3 ms·

>I had a problem which required finding the minimum distance between an ellipse and a point. If that had been my day job ... OT: I had the problem as part of m

>I had a problem which required finding the minimum distance between an ellipse and a point. If that had been my day job ... OT: I had the problem as part of my 1st day job in 1990 on Parasolid. IIRC, there was more than one way of converting the problem into solving a quartic polynomial, which is solvable by radicals and so doesn't need Newton's method. If the point is inside the evolute of the ellipse then the quartic has four real solutions, which gives two local minima and two local maxima of distance. If the point is outside or on the evolute of the ellipse then the quartic has two real solutions, which gives one local minima and one local maxima of distance. Having said that, all methods using quartic polynomials could be unstable if the ellipse was pathological, for example very long and thin. Maybe that would be a time to revert to Newton's method.

0xfaded 9y ago
There are certainly multiple ways of doing it. I'm curious, do you remember if the specific quartic had a nice "trick" to the formulation? I played with it a little and couldn't see anything that would cancel in a nice way. I try to avoid using quadratic solutions in code because of numeric instability issues. Granted I never implemented the direct solution, but the number of opportunities for catastrophic cancellation gave me shivers.
- delhanty 9y ago
  I don't think there was any trick. One approach is to use the rational quadratic parametrization of the ellipse and substitute that into the squared distance function to get a quartic. >I try to avoid using quadratic solutions in code because of numeric instability issues. Sensible! Back then the quartic solver had been written already and extensively tested by someone else - probably for line X torus intersection. Even so there were numerical problems introduced because the quartic was in power basis and so the coefficients didn't have geometric meaning. I guess phenomena like Wilkinson's polynomial could occur. [0] If I was doing it today, I would probably proceed as follows: 1. special case for the point exactly on an axis. 2. otherwise, choose the quadrant containing the point and represent that as a rational quadratic Bezier curve. This parametrizes the quadrant from 0 to 1 and is guaranteed to contain the global minimum. 3. substitute the quadratic Bezier curve into the square distance function to get a Bernstein basis quartic polynomial. 4. differentiate to get a Bernstein basic cubic. 5. solve the roots of the cubic numerically in the range 0 to 1: either Newton's or clipping or some other hybrid numerical method. A good reference is [1]. [0] https://en.wikipedia.org/wiki/Wilkinson%27s_polynomial https://en.wikipedia.org/wiki/Wilkinson%27s_polynomial [1] Shape Interrogation for Computer Aided Design and Manufacture by Nichola M. Patrikalakis and Takashi Maekawa
  - mkl 9y ago
    > substitute that into the squared distance function How the heck did "" get in there? > 5. solve the roots of the cubic numerically in the range 0 to 1: either Newton's or clipping or some other hybrid numerical method. Is there a reason not to use the cubic formula in this situation?
    - delhanty 9y ago
      >How the heck did "" get in there? Something like (E(t)-P).(E(t)-P) where E is the ellipse and P is the point. >Is there a reason not to use the cubic formula in this situation? Numerical stability in edge cases, which will occur for code exercised as extensively as Parasolid's. I guess you could a Bernstein version of the cubic formula, though I've not seen one derived in the wild. It would be nice to have as well as Bernstein quartic formula. The Bernstein quadratic is easy.

0xfaded 9y ago

There are certainly multiple ways of doing it. I'm curious, do you remember if the specific quartic had a nice "trick" to the formulation? I played with it a little and couldn't see anything that would cancel in a nice way. I try to avoid using quadratic solutions in code because of numeric instability issues. Granted I never implemented the direct solution, but the number of opportunities for catastrophic cancellation gave me shivers.

delhanty 9y ago
I don't think there was any trick. One approach is to use the rational quadratic parametrization of the ellipse and substitute that into the squared distance function to get a quartic. >I try to avoid using quadratic solutions in code because of numeric instability issues. Sensible! Back then the quartic solver had been written already and extensively tested by someone else - probably for line X torus intersection. Even so there were numerical problems introduced because the quartic was in power basis and so the coefficients didn't have geometric meaning. I guess phenomena like Wilkinson's polynomial could occur. [0] If I was doing it today, I would probably proceed as follows: 1. special case for the point exactly on an axis. 2. otherwise, choose the quadrant containing the point and represent that as a rational quadratic Bezier curve. This parametrizes the quadrant from 0 to 1 and is guaranteed to contain the global minimum. 3. substitute the quadratic Bezier curve into the square distance function to get a Bernstein basis quartic polynomial. 4. differentiate to get a Bernstein basic cubic. 5. solve the roots of the cubic numerically in the range 0 to 1: either Newton's or clipping or some other hybrid numerical method. A good reference is [1]. [0] https://en.wikipedia.org/wiki/Wilkinson%27s_polynomial https://en.wikipedia.org/wiki/Wilkinson%27s_polynomial [1] Shape Interrogation for Computer Aided Design and Manufacture by Nichola M. Patrikalakis and Takashi Maekawa
- mkl 9y ago
  > substitute that into the squared distance function How the heck did "" get in there? > 5. solve the roots of the cubic numerically in the range 0 to 1: either Newton's or clipping or some other hybrid numerical method. Is there a reason not to use the cubic formula in this situation?
  - delhanty 9y ago
    >How the heck did "" get in there? Something like (E(t)-P).(E(t)-P) where E is the ellipse and P is the point. >Is there a reason not to use the cubic formula in this situation? Numerical stability in edge cases, which will occur for code exercised as extensively as Parasolid's. I guess you could a Bernstein version of the cubic formula, though I've not seen one derived in the wild. It would be nice to have as well as Bernstein quartic formula. The Bernstein quadratic is easy.