REllipse.getVectorTo(p) exceptions (FS#2564)

[CVH]

Andrew,

REllipse.getVectorTo(p) or defining the nearest point on an ellipse from any arbitrary point may fail or is not very accurate.

This has major implications, as getVectorTo() is an important resource on which many other resources are based.

• The distance to
  A perpendicular point (One, while there are 2-4)
  isTangent
  isOnShape
  ... and so on.

Several related resources solve this with an increased (but fixed) tolerance ... Up to 1e-4 or 100000 times RS.PointTolerance

A) Incorrect results: FS#2564
The method typically fails for points on the major Axis and then strictly inside the evolute of the ellipse.

This can be reproduced in the attached DXF example.
- Ensure that layer 0 is the current active layer (All other are locked)
- Start 'Line from 2 Points' (LI) in segment mode and select any point 3 or the center point as first point.
- Activate the perpendicular snapping option (SU) and indicate anywhere near the ellipse shape.

FalsePerpendicularOnEllipse.dxf

(117.67 KiB) Downloaded 270 times

The preview will be a line towards the nearest major point and that is most incorrect.
The correct but also dual solution in green for 3 and C are visualized in the related layers under 'X_Solutions'.
Toggle the layer visibility to see them.

For the center it should be more than obvious that both minor points are the nearest perpendicular points.
For all points on the major axis, strictly inside the evolute, there are 2 mirrored solutions, both are equally 'far' or equally 'near'.
The problem here is that on calling getVectorTo() for any shape type, none or one singular result is expected.
A misconception when handling Circles, Ellipses or Arcs.

This can be solved with an extra flag for these cases: singular = true/false.
By default, true, returning a singular result or none on duality (Similar as with a circle center).
For distance related resources, true, only returning 1 of both solutions at equal distance.
False for special cases, for example for perpendicular points although still incomplete for the ellipse case.


B) Less accurate:
getVectorTo() for an ellipse is based on the Newton's method to converge to a solution within 32 (costly) iterations to less than dEpsilon=1.0e-8.
A common used method but there are 'Practical considerations' pointed out in the WIKI.

Verifying that the returned RVector is on the ellipse shape is usually NOT very satisfactory.
For any point on an generalized ellipse solving x²/a²+y²/b² must be equal to 1.0 (Equation of an ellipse with a/b = major/minor radius)
Values between 1.0 and 10.0 or even larger are not uncommon.
This error will trickle down to all resources that are based on this key resource.

Any difference from 1.0 cannot be used as a measurable tolerance but very near 1.0 should be the outcome for any point on the ellipse.
How close to 1.0? Essentially as close a possible whiting the number system.


Solution:
Extensively testing a 'Simple Method for Distance to Ellipse' by Carl Chatfield for more than 6 months now.
Essentially based on the fact that any normal to an ellipse is also tangent to its evolute.

GitHub repository (MIT License): https://github.com/0xfaded/ellipse_demo ... atfield.io
Exploiting the 'trig-free modification' by Adrian Stephens is a very cost effective approach that does not compromise on very high accuracy.
What in turn is later adopted in the GitHub repository.

I implemented this as JS with max 12 iterations halting when X and Y have converged to within 1e-15 for a unified ellipse.
Typically less than 4 iterations are required and 'onShape' is by default as true as can be.
Scaling this back (up) maintains the relative accuracy and larger coordinate values ​​cannot be more accurate than that (Give or take).

There are some counter indications, the method uses 1/4 ellipse what means that both semi-axes are the limit cases.
More than 12 iterations are required for points nearing the cusps of the evolute.
It is also less stable converging very near the axes.

This is pre-solved as one special case using the equation of any normal to an ellipse for a point at angle t on the ellipse.
ax/cos(t)-by/sin(t) = a²-b² reduces to: cos(t) = ax/(a²-b²) for y = 0 (Major axis) and |x| < (a²-b²)/a.
Further away from the center than |x| the single solution is the nearest major point.
Otherwise cos(t) is well defined, thus sin(t) (+/-), thus the point on the ellipse and its mirror are known exactly.
=> (X = a * cos(t), Y = b * +/-sin(t)) for any point (|x|<(a²-b²)/a, y≈0).

This approach spares us the trouble of starting to guess ... And converge.
I transferred this from a method to define all normal points, exploiting it for both axes.
Solving a Quartic equation was also less stable near the axes ... See some results in the related layers under 'X_Solutions (All)'.
This tool also returns the 'far-side' normals or perpendicular points for snapping, tangent circles/arcs and so on.
Related topics: Exceptions with AT (ArcTPR.js), Auxiliary shapes, ...


In the perspective that a circle (arc) is a special case of an ellipse (arc), then the 'far-side' perpendicular point is also omitted.
See feature request FS#2621. Related topic by 'A user': Perpendicular / Tangential Snap

PerpendicularToArc.png (604.04 KiB) Viewed 29261 times

Adding a copy of ArcTPR-issues.png for a 'far-side' perpendicular point P' on a circle and the 'far-side' perpendicular points ABC.

ArcTPR-Issue.png (24.07 KiB) Viewed 29261 times

Code will be added in due time ... Only as JS, I think that the conversion to C++ should not be an issue.

Regards,
CVH

[CVH]

As promised and as proof of concept.

Start a new document and run the attached code using: Misc .. 'Run Script' XC:

POC_NearestPointOnEllipse.js

(20.43 KiB) Downloaded 311 times

An ellipse is drawn with 20 trial points (magenta).
From each test point the normal to the nearest point on the ellipse is generated.
Each normal as line segment is validated and that info is included as custom properties.
This data is listed in the Property Editor when selecting a 'Normal'.

The segment is always created from a solution of the alternative method.
In red when the standard resource REllipse.getVectorTo(point) would fail.
-> Typical for points on the major axis and inside the evolute.
In cyan when the standard resource returns something valid.
-> Most likely less accurate.

'Passing' or not for the standard resource is verified as within a larger tolerance of 1 unit for X and Y.
'On ellipse' is verified in two distinct ways:
-> x²/a²+y²/b² must be equal to 1.0, anything within 1e-14 is near perfect.
-> The sum of the distances to the foci must be 2a ≈ 87.8, anything within 1e-13 is near perfect.
Orientation (0-2Pi) of the result is verified in regards with the bisector of two lines to the foci, anything within 1e-14 is near perfect.

Except where the standard resource fails, its largest errors occur at a larger change in curvature.
E.g.: See details for the segment related to the most left-lower point.


getNearestPointsOnFullREllipse2D() is the part that finds the nearest point on an ellipse.
-> Nothing with invalid data
-> The ellipse center when Major radius <= 1.0e-6
-> RCircle method when almost circular (ratio >= 0.999999)
-> Both minor points when the point is near the center
-> 3 cases for a point very near the ellipse major axis
->-> Degenerated ellipse (Minor axis <= 1.0e-6)
->-> Outside the evolute
->-> None of the above
-> 1 case for a point very near the ellipse minor axis
-> Solution by converging for none of the above (Line 163 - 251) Alternative code adapted to QCAD.

getVectorToREllipse2D is the JS analog of REllipse::getVectorTo() with the extra flag in case of a duality.

validateNormalPoint is a helper function for the validation in this POC, comparing results with the original resource.

main() is the driver code that creates the document shapes including their validation by the above.
Calling getVectorToREllipse2D based on getNearestPointsOnFullREllipse2D() for each of the 20 points.

The alternative method is found to be less complex and faster, more exact and does not fail.
The downside is that dual solutions exists while getVectorTo() is geared to a single solution, valid or not.

Regards,
CVH