36

I have an ellipse, defined by Center Point, radiusX and radiusY, and I have a Point. I want to find the point on the ellipse that is closest to the given point. In the illustration below, that would be S1.

graph1

Now I already have code, but there is a logical error somewhere in it, and I seem to be unable to find it. I broke the problem down to the following code example:

#include <vector>
#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <math.h>

using namespace std;

void dostuff();

int main()
{
    dostuff();
    return 0;
}

typedef std::vector<cv::Point> vectorOfCvPoints;

void dostuff()
{

    const double ellipseCenterX = 250;
    const double ellipseCenterY = 250;
    const double ellipseRadiusX = 150;
    const double ellipseRadiusY = 100;

    vectorOfCvPoints datapoints;

    for (int i = 0; i < 360; i+=5)
    {
        double angle = i / 180.0 * CV_PI;
        double x = ellipseRadiusX * cos(angle);
        double y = ellipseRadiusY * sin(angle);
        x *= 1.4;
        y *= 1.4;
        x += ellipseCenterX;
        y += ellipseCenterY;
        datapoints.push_back(cv::Point(x,y));
    }

    cv::Mat drawing = cv::Mat::zeros( 500, 500, CV_8UC1 );

    for (int i = 0; i < datapoints.size(); i++)
    {
        const cv::Point & curPoint = datapoints[i];
        const double curPointX = curPoint.x;
        const double curPointY = curPoint.y * -1; //transform from image coordinates to geometric coordinates

        double angleToEllipseCenter = atan2(curPointY - ellipseCenterY * -1, curPointX - ellipseCenterX); //ellipseCenterY * -1 for transformation to geometric coords (from image coords)

        double nearestEllipseX = ellipseCenterX + ellipseRadiusX * cos(angleToEllipseCenter);
        double nearestEllipseY = ellipseCenterY * -1 + ellipseRadiusY * sin(angleToEllipseCenter); //ellipseCenterY * -1 for transformation to geometric coords (from image coords)


        cv::Point center(ellipseCenterX, ellipseCenterY);
        cv::Size axes(ellipseRadiusX, ellipseRadiusY);
        cv::ellipse(drawing, center, axes, 0, 0, 360, cv::Scalar(255));
        cv::line(drawing, curPoint, cv::Point(nearestEllipseX,nearestEllipseY*-1), cv::Scalar(180));

    }
    cv::namedWindow( "ellipse", CV_WINDOW_AUTOSIZE );
    cv::imshow( "ellipse", drawing );
    cv::waitKey(0);
}

It produces the following image:

snapshot1

You can see that it actually finds "near" points on the ellipse, but it are not the "nearest" points. What I intentionally want is this: (excuse my poor drawing)

snapshot2

would you extent the lines in the last image, they would cross the center of the ellipse, but this is not the case for the lines in the previous image.
I hope you get the picture. Can anyone tell me what I am doing wrong?

user2950911
  • 873
  • 3
  • 12
  • 19
  • it will be easier if You just describe your method for finding the point than the actual code – rank1 Apr 09 '14 at 10:30

9 Answers9

43

Consider a bounding circle around the given point (c, d), which passes through the nearest point on the ellipse. From the diagram it is clear that the closest point is such that a line drawn from it to the given point must be perpendicular to the shared tangent of the ellipse and circle. Any other points would be outside the circle and so must be further away from the given point.

enter image description here

So the point you are looking for is not the intersection between the line and the ellipse, but the point (x, y) in the diagram.

Gradient of tangent:

enter image description here

Gradient of line:

enter image description here

Condition for perpedicular lines - product of gradients = -1:

enter image description here

enter image description here

enter image description here

When rearranged and substituted into the equation of your ellipse...

enter image description here

...this will give two nasty quartic (4th-degree polynomial) equations in terms of either x or y. AFAIK there are no general analytical (exact algebraic) methods to solve them. You could try an iterative method - look up the Newton-Raphson iterative root-finding algorithm.

Take a look at this very good paper on the subject: http://www.spaceroots.org/documents/distance/distance-to-ellipse.pdf

Sorry for the incomplete answer - I totally blame the laws of mathematics and nature...

EDIT: oops, i seem to have a and b the wrong way round in the diagram xD

  • 1
    hell, you are right. That solves the problem of my faulty implementation, not by finding the error, but by exposing my approach as being faulty in the first place. Thanks – user2950911 Apr 09 '14 at 12:04
  • @user2950911 Glad to be of help –  Apr 09 '14 at 12:05
  • @willywonka_dailyblah would the following also work? Create a transformation matrix which scales the ellipse to a circle, while keeping the center the same. Applying the transformation matrix to the point p1. Calculating the intersection s1 of the point to the center of the circle. Transforming p1 and s1 back and calculating their distance. – mfuchs Jul 31 '15 at 12:39
  • 1
    @mat69 I like your thinking, but no. The transformation matrix you're thinking of is a scaling matrix, which is a linear transformation. For a circle the closest point would indeed be the intersection he assumed - but when the inverse transformation is applied, all points on that line would be scaled by the same amount - leading to the wrong answer he obtained _again_. –  Jul 31 '15 at 15:07
  • 4
    Quartic polynomials do have exact solutions (they're the highest degree to have them, see Galois Theory). They're non-trivial, though, so it might be worth investigating if the polynomial is reducible. – MSalters Feb 16 '16 at 14:19
  • 1
    @MSalters you're right, but doing so isn't worth the computational effort compared to iterative convergence –  Feb 18 '16 at 20:20
  • This ans is not readable in dark mode... :( – echefede Jun 21 '21 at 20:36
  • On the opposite, there *is* an analytical formula to solve the quartics. Well, it is a little complicated (involving the resolution of a cubic and three quadratics), but workable anyway. – Yves Daoust Apr 26 '23 at 08:43
22

There is a relatively simple numerical method with better convergence than Newtons Method. I have a blog post about why it works http://wet-robots.ghost.io/simple-method-for-distance-to-ellipse/

This implementation works without any trig functions:

def solve(semi_major, semi_minor, p):  
    px = abs(p[0])
    py = abs(p[1])

    tx = 0.707
    ty = 0.707

    a = semi_major
    b = semi_minor

    for x in range(0, 3):
        x = a * tx
        y = b * ty

        ex = (a*a - b*b) * tx**3 / a
        ey = (b*b - a*a) * ty**3 / b

        rx = x - ex
        ry = y - ey

        qx = px - ex
        qy = py - ey

        r = math.hypot(ry, rx)
        q = math.hypot(qy, qx)

        tx = min(1, max(0, (qx * r / q + ex) / a))
        ty = min(1, max(0, (qy * r / q + ey) / b))
        t = math.hypot(ty, tx)
        tx /= t 
        ty /= t 

    return (math.copysign(a * tx, p[0]), math.copysign(b * ty, p[1]))

Convergence

Credit to Adrian Stephens for the Trig-Free Optimization.

Johannes
  • 6,232
  • 9
  • 43
  • 59
user364952
  • 570
  • 4
  • 10
  • 2
    This approach is excellent. An optimization can be made that completely avoids any trig functions (all credit to the issue posted on the github repo for this blog post: https://github.com/0xfaded/ellipse_demo/issues/1 ). I made an implementation that uses this optimization in C# for Unity3D: https://gist.github.com/JohannesMP/777bdc8e84df6ddfeaa4f0ddb1c7adb3 – Johannes Jul 13 '18 at 06:16
  • Code doesn't seem to be complete after your edit. E.g. `tx` is not defined before its first use. – JHBonarius Jul 13 '18 at 07:15
  • Whoops, didn't copy the optimization correctly. Thanks for checking @JHBonarius – Johannes Jul 14 '18 at 01:26
  • 1
    Should this line be changed from "for x in range(0, 3):" to "for i in range(0, 3):"? – quarkz Jan 07 '22 at 09:00
4

Here is the code translated to C# implemented from this paper to solve for the ellipse: http://www.geometrictools.com/Documentation/DistancePointEllipseEllipsoid.pdf

Note that this code is untested - if you find any errors let me know.

    //Pseudocode for robustly computing the closest ellipse point and distance to a query point. It
    //is required that e0 >= e1 > 0, y0 >= 0, and y1 >= 0.
    //e0,e1 = ellipse dimension 0 and 1, where 0 is greater and both are positive.
    //y0,y1 = initial point on ellipse axis (center of ellipse is 0,0)
    //x0,x1 = intersection point

    double GetRoot ( double r0 , double z0 , double z1 , double g )
    {
        double n0 = r0*z0;
        double s0 = z1 - 1; 
        double s1 = ( g < 0 ? 0 : Math.Sqrt(n0*n0+z1*z1) - 1 ) ;
        double s = 0;
        for ( int i = 0; i < maxIter; ++i ){
            s = ( s0 + s1 ) / 2 ;
            if ( s == s0 || s == s1 ) {break; }
            double ratio0 = n0 /( s + r0 );
            double ratio1 = z1 /( s + 1 );
            g = ratio0*ratio0 + ratio1*ratio1 - 1 ;
            if (g > 0) {s0 = s;} else if (g < 0) {s1 = s ;} else {break ;}
        }
        return s;
    }
    double DistancePointEllipse( double e0 , double e1 , double y0 , double y1 , out double x0 , out double x1)
    {
        double distance;
        if ( y1 > 0){
            if ( y0 > 0){
                double z0 = y0 / e0; 
                double z1 = y1 / e1; 
                double g = z0*z0+z1*z1 - 1;
                if ( g != 0){
                    double r0 = (e0/e1)*(e0/e1);
                    double sbar = GetRoot(r0 , z0 , z1 , g);
                    x0 = r0 * y0 /( sbar + r0 );
                    x1 = y1 /( sbar + 1 );
                    distance = Math.Sqrt( (x0-y0)*(x0-y0) + (x1-y1)*(x1-y1) );
                    }else{
                        x0 = y0; 
                        x1 = y1;
                        distance = 0;
                    }
                }
                else // y0 == 0
                    x0 = 0 ; x1 = e1 ; distance = Math.Abs( y1 - e1 );
        }else{ // y1 == 0
            double numer0 = e0*y0 , denom0 = e0*e0 - e1*e1;
            if ( numer0 < denom0 ){
                    double xde0 = numer0/denom0;
                    x0 = e0*xde0 ; x1 = e1*Math.Sqrt(1 - xde0*xde0 );
                    distance = Math.Sqrt( (x0-y0)*(x0-y0) + x1*x1 );
                }else{
                    x0 = e0; 
                    x1 = 0; 
                    distance = Math.Abs( y0 - e0 );
            }
        }
        return distance;
    }
johnml1135
  • 4,519
  • 3
  • 23
  • 18
3

The following python code implements the equations described at "Distance from a Point to an Ellipse" and uses newton's method to find the roots and from that the closest point on the ellipse to the point.

Unfortunately, as can be seen from the example, it seems to only be accurate outside the ellipse. Within the ellipse weird things happen.

from math import sin, cos, atan2, pi, fabs


def ellipe_tan_dot(rx, ry, px, py, theta):
    '''Dot product of the equation of the line formed by the point
    with another point on the ellipse's boundary and the tangent of the ellipse
    at that point on the boundary.
    '''
    return ((rx ** 2 - ry ** 2) * cos(theta) * sin(theta) -
            px * rx * sin(theta) + py * ry * cos(theta))


def ellipe_tan_dot_derivative(rx, ry, px, py, theta):
    '''The derivative of ellipe_tan_dot.
    '''
    return ((rx ** 2 - ry ** 2) * (cos(theta) ** 2 - sin(theta) ** 2) -
            px * rx * cos(theta) - py * ry * sin(theta))


def estimate_distance(x, y, rx, ry, x0=0, y0=0, angle=0, error=1e-5):
    '''Given a point (x, y), and an ellipse with major - minor axis (rx, ry),
    its center at (x0, y0), and with a counter clockwise rotation of
    `angle` degrees, will return the distance between the ellipse and the
    closest point on the ellipses boundary.
    '''
    x -= x0
    y -= y0
    if angle:
        # rotate the points onto an ellipse whose rx, and ry lay on the x, y
        # axis
        angle = -pi / 180. * angle
        x, y = x * cos(angle) - y * sin(angle), x * sin(angle) + y * cos(angle)

    theta = atan2(rx * y, ry * x)
    while fabs(ellipe_tan_dot(rx, ry, x, y, theta)) > error:
        theta -= ellipe_tan_dot(
            rx, ry, x, y, theta) / \
            ellipe_tan_dot_derivative(rx, ry, x, y, theta)

    px, py = rx * cos(theta), ry * sin(theta)
    return ((x - px) ** 2 + (y - py) ** 2) ** .5

Here's an example:

rx, ry = 12, 35  # major, minor ellipse axis
x0 = y0 = 50  # center point of the ellipse
angle = 45  # ellipse's rotation counter clockwise
sx, sy = s = 100, 100  # size of the canvas background

dist = np.zeros(s)
for x in range(sx):
    for y in range(sy):
        dist[x, y] = estimate_distance(x, y, rx, ry, x0, y0, angle)

plt.imshow(dist.T, extent=(0, sx, 0, sy), origin="lower")
plt.colorbar()
ax = plt.gca()
ellipse = Ellipse(xy=(x0, y0), width=2 * rx, height=2 * ry, angle=angle,
                  edgecolor='r', fc='None', linestyle='dashed')
ax.add_patch(ellipse)
plt.show()

Which generates rotated ellipse an ellipse and the distance from the boundary of the ellipse as a heat map. As can be seen, at the boundary the distance is zero (deep blue).

Matt
  • 1,327
  • 1
  • 12
  • 21
  • Be warned that the paper has some mistakes. Specifically, in the 3D section, the formula given for (P-X) dot dX/dtheta is incorrect. I verified by hand and again with Mathematica. There also seem to be factors of "r" missing in various places -- maybe less of a problem since you can set r to 1 without loss of generality. – Paul Du Bois Dec 16 '16 at 22:13
  • The code above is for 2D which is correct, though. Even for non-1 r. – Matt Dec 17 '16 at 02:25
  • Regarding the 2D case, the paper conflates "solution to (delta dot tangent)==0" with "closest point". There are 4 solutions inside the ellipse. What your image shows is that the paper's suggestion for initial theta sometimes leads to the incorrect root. It would probably be interesting to plot the direction of the solution found rather than the distance. – Paul Du Bois Dec 17 '16 at 23:36
  • But outside the ellipse it's valid, right? I guess that's really what I needed it for. – Matt Dec 20 '16 at 05:08
2

Given an ellipse E in parametric form and a point P

eqn

the square of the distance between P and E(t) is

eqn

The minimum must satisfy

eqn

Using the trigonometric identities

eqn

and substituting

enter image description here

yields the following quartic equation:

enter image description here

Here's an example C function that solves the quartic directly and computes sin(t) and cos(t) for the nearest point on the ellipse:

void nearest(double a, double b, double x, double y, double *ecos_ret, double *esin_ret) {
    double ax = fabs(a*x);
    double by = fabs(b*y);
    double r  = b*b - a*a;
    double c, d;
    int switched = 0;

    if (ax <= by) {
        if (by == 0) {
            if (r >= 0) { *ecos_ret = 1; *esin_ret = 0; }
            else        { *ecos_ret = 0; *esin_ret = 1; }
            return;
        }
        c = (ax - r) / by;
        d = (ax + r) / by;
    } else {
        c = (by + r) / ax;
        d = (by - r) / ax;
        switched = 1;
    }

    double cc = c*c;
    double D0 = 12*(c*d + 1);      // *-4
    double D1 = 54*(d*d - cc);     // *4
    double D  = D1*D1 + D0*D0*D0;  // *16

    double St;
    if (D < 0) {
        double t = sqrt(-D0);             // *2
        double phi = acos(D1 / (t*t*t));
        St = 2*t*cos((1.0/3)*phi);        // *2
    } else {
        double Q = cbrt(D1 + sqrt(D));    // *2
        St = Q - D0 / Q;                  // *2
    }

    double p    = 3*cc;                          // *-2
    double SS   = (1.0/3)*(p + St);              // *4
    double S    = sqrt(SS);                      // *2
    double q    = 2*cc*c + 4*d;                  // *2
    double l    = sqrt(p - SS + q / S) - S - c;  // *2
    double ll   = l*l;                           // *4
    double ll4  = ll + 4;                        // *4
    double esin = (4*l)    / ll4;
    double ecos = (4 - ll) / ll4;

    if (switched) {
        double t = esin;
        esin = ecos;
        ecos = t;
    }

    *ecos_ret = copysign(ecos, a*x);
    *esin_ret = copysign(esin, b*y);
}

Try it online!

nwellnhof
  • 32,319
  • 7
  • 89
  • 113
0

You just need to calculate the intersection of the line [P1,P0] to your elipse which is S1.

If the line equeation is:

enter image description here

and the elipse equesion is:

elipse equesion

than the values of S1 will be:

enter image description here

Now you just need to calculate the distance between S1 to P1 , the formula (for A,B points) is: distance

Idan Yehuda
  • 524
  • 7
  • 21
  • How can 2D point coordinates have an undefined leading sign? That would leave 4 possible points. I cant quite follow how you have derived the S2 coordinate equations, also it seems you have inserted the wrong image for the line equation, which i assume is something like y=mx+t or y=((y1-y0)/(x1-x0)) * (x-x0) + y0 – user2950911 Apr 09 '14 at 11:13
  • 2
    As user3235832 showed, this is wrong as it's not going to lie on a straight line from the point. – Matt Dec 17 '16 at 02:26
0

I've solved the distance issue via focal points.

For every point on the ellipse

r1 + r2 = 2*a0

where

r1 - Euclidean distance from the given point to focal point 1

r2 - Euclidean distance from the given point to focal point 2

a0 - semimajor axis length

I can also calculate the r1 and r2 for any given point which gives me another ellipse that this point lies on that is concentric to the given ellipse. So the distance is

d = Abs((r1 + r2) / 2 - a0)

  • 1
    This post suggests that your answer doesn't give the correct result: https://math.stackexchange.com/questions/1897587/is-minimum-distance-between-points-on-concentric-ellipses-constant – user664303 Jul 01 '19 at 23:53
0

As propposed by user3235832 you shall solve quartic equation to find the normal to the ellipse (https://www.mathpages.com/home/kmath505/kmath505.htm). With good initial value only few iterations are needed (I use it myself). As an initial value I use S1 from your picture.

Hynek
  • 1
0

The fastest method I guess is http://wwwf.imperial.ac.uk/~rn/distance2ellipse.pdf Which has been mentioned also by Matt but as he found out the method doesn't work very well inside of ellipse.

The problem is the theta initialization. I proposed an stable initialization:

  1. Find the intersection of ellipse and horizontal line passing the point.
  2. Find the other intersection using vertical line.
  3. Choose the one that is closer the point.
  4. Calculate the initial angle based on that point.

I got good results with no issue inside and outside:

As you can see in the following image it just iterated about 3 times to reach 1e-8. Close to axis it is 1 iteration.

The C++ code is here:

double initialAngle(double a, double b, double x, double y) {
    auto abs_x = fabs(x);
    auto abs_y = fabs(y);
    bool isOutside = false;
    if (abs_x > a || abs_y > b) isOutside = true;

    double xd, yd;
    if (!isOutside) {
        xd = sqrt((1.0 - y * y / (b * b)) * (a * a));
        if (abs_x > xd)
            isOutside = true;
        else {
            yd = sqrt((1.0 - x * x / (a * a)) * (b * b));
            if (abs_y > yd)
                isOutside = true;
        }
    }
    double t;
    if (isOutside)
        t = atan2(a * y, b * x); //The point is outside of ellipse
    else {
        //The point is inside
        if (xd < yd) {
            if (x < 0) xd = -xd;
            t = atan2(y, xd);
        }
        else {
            if (y < 0) yd = -yd;
            t = atan2(yd, x);
        }
    }
    return t;
}


double distanceToElipse(double a, double b, double x, double y, int maxIter = 10, double maxError = 1e-5) {
    //std::cout <<"p="<< x << "," << y << std::endl;

    
    auto a2mb2 = a * a - b * b;
    double t = initialAngle(a, b, x, y);
    auto ct = cos(t);
    auto st = sin(t);

    int i;
    double err;
    for (i = 0; i < maxIter; i++) {
        auto f = a2mb2 * ct * st - x * a * st + y * b * ct;
        auto fp = a2mb2 * (ct * ct - st * st) - x * a * ct - y * b * st;
        auto t2 = t - f / fp;
        err = fabs(t2 - t);
        //std::cout << i + 1 << " " << err << std::endl;
        t = t2;
        ct = cos(t);
        st = sin(t);
        if (err < maxError) break;
    }
    auto dx = a * ct - x;
    auto dy = b * st - y;
    //std::cout << a * ct << "," << b * st << std::endl;
    return   sqrt(dx * dx + dy * dy);
}