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 (118 KB)
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): GitHub - 0xfaded/ellipse_demo at blog.chatfield.io · GitHub
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

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

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.4 KB)
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

[andrew]

I let Claude look into this today and here’s its reply:

Thanks CVH — I went through the proof of concept and it is now implemented in REllipse::getVectorTo().

Your diagnosis is confirmed, and the problem is worse than “not very accurate”. Measured against a 60 digit reference (bisection on (a·u/(t+a²))² + (b·v/(t+b²))² = 1, which is strictly monotonic for t > -b²), the old Newton iteration does not merely lose precision — from its starting point t0 = b(v-b) it converges to a wrong root for points near the major axis and near the evolute. Worst error relative to the major radius:

┌─────────────────────┬─────────┬─────────┐
│       regime        │   old   │   new   │
├─────────────────────┼─────────┼─────────┤
│ general             │ 5.5e-05 │ 5.9e-12 │
├─────────────────────┼─────────┼─────────┤
│ near the curve      │ 1.8e-06 │ 7.2e-13 │
├─────────────────────┼─────────┼─────────┤
│ near the evolute    │ 6.5e-04 │ 1.6e-11 │
├─────────────────────┼─────────┼─────────┤
│ near the major axis │ 3.0     │ 2.2e-09 │
├─────────────────────┼─────────┼─────────┤
│ near the minor axis │ 4.8e-09 │ 1.0e-09 │
├─────────────────────┼─────────┼─────────┤
│ eccentric           │ 2.4e-01 │ 3.8e-13 │
├─────────────────────┼─────────┼─────────┤
│ near circular       │ 1.0e-06 │ 1.0e-12 │
├─────────────────────┼─────────┼─────────┤
│ far away            │ 4.8e-09 │ 1.1e-10 │
└─────────────────────┴─────────┴─────────┘

A worst case of 3.0 means the returned point was three major radii away from the ellipse. 4600 random cases over 8 regimes, random centre and rotation, both columns from the actual compiled REllipse.

Two things were added on top of your iteration:

A Newton refinement of the parametric angle on h(t) = (b²-a²)/2·sin(2t) + a|u|sin(t) - b|v|cos(t), the derivative of the squared distance. Your iteration is the one regime where the old code was competitive: for eccentric ellipses it is off by up to 5.6e-4. Most of that turned out to be non-convergence within the 12 iterations rather than a precision limit — 24 iterations already gets to 4.3e-9 — but two refinement steps take it to 3.8e-13 and cost less. So the iteration supplies the globally correct starting point and the refinement supplies the precision.

Cases neither version handled. Degenerate ellipses (ratio 0, or major radius 0) returned an invalid vector; now the collapsed segment resp. the centre is used. Ellipses with a ratio greater than one took the axis shortcut that is only correct for a ratio below one; the two axes are now swapped for the whole computation (worst deviation from the ellipse equation 4.2e-05 → 3.3e-13).

Your handling of points in line with the major axis, including the two mirrored solutions, is kept as it was — it is needed for ellipse arcs, where the mirrored point may still be on the arc. I made its tolerance relative to the major radius so it behaves the same for small and large ellipses. For arcs, the result is valid in exactly the same cases as before (4000 random arcs, zero divergences) and never further from the point: against a brute force scan the old code overshot by up to 1.2e-8, the new one by 0.

I did not need the near circular shortcut — the iteration converges on a circle in a single step. Which is just as well, because that branch references an undefined majorR (the variable is called a), so it would throw if it were ever reached. Worth fixing in your copy.

Cost is 192ns → 268ns per call, which seems a fair trade.