website articles
distance to an ellipse
Often, ellipses are usually not much more difficult than disks/circles, and one can compute lots of their properties analytically and easily, like their bounding boxes and intersections. However, some properties turn out to be a bit more complicated to compute than one would expect. One such property is the distance from a point to an ellipse, in 2D!


The maths




So, for a given point in the plane and a parametric ellipse e centered in the origin with radiuses a and b in the x and y axes , the distance from z to the ellipse is



which expands to



The closest point will be given by



which expands in this case to



This calls for the change of variable , so that



which expands into a quartic equation with the following coefficients:



or by dividing the equation by the leading coefficient (the only condition to be able to do so is that the ellipse is not, in fact, a circle), and by then renaming ,



From this point on we proceed to resolve the quartic equation analytically with the standard method. Due to some of the symmetries in the coefficients, we get some pretty nice simplifications though. For example, the resolvent cubic equation loses its linear term and then depresses to



with






The code




Due to all the symmetries mentioned, and the fact that we know we are looking only for real roots (possibly only one, if we are measuring the distance from a point outside the elliposid), the code gets quite compact. Also, we use the symmetry of the geometrical problem and we force our point to be on the first quadrant of the coordinate system (x>0 and y>0).

float sdEllipse( in vec2 z, in vec2 ab )
{
    vec3 p = abs( z ); if( p.x > p.y ){ p=p.yx; ab=ab.yx; }
	
    float l = ab.y*ab.y - ab.x*ab.x;
    float m = ab.x*p.x/l; float m2 = m*m;
    float n = ab.y*p.y/l; float n2 = n*n;
    float c = (m2 + n2 - 1.0)/3.0; float c3 = c*c*c;
    float q = c3 + m2*n2*2.0;
    float d = c3 + m2*n2;
    float g = m + m*n2;

    float co;

    if( d<0.0 )
    {
        float p = acos(q/c3)/3.0;
        float s = cos(p);
        float t = sin(p)*sqrt(3.0);
        float rx = sqrt( -c*(s + t + 2.0) + m2 );
        float ry = sqrt( -c*(s - t + 2.0) + m2 );
        co = ( ry + sign(l)*rx + abs(g)/(rx*ry) - m)/2.0;
    }
    else
    {
        float h = 2.0*m*n*sqrt( d );
        float s = sign(q+h)*pow( abs(q+h), 1.0/3.0 );
        float u = sign(q-h)*pow( abs(q-h), 1.0/3.0 );
        float rx = -s - u - c*4.0 + 2.0*m2;
        float ry = (s - u)*sqrt(3.0);
        float rm = sqrt( rx*rx + ry*ry );
        float p = ry/sqrt(rm-rx);
        co = (p + 2.0*g/rm - m)/2.0;
    }

    float si = sqrt( 1.0 - co*co );
 
    vec2 closestPoint = vec2( ab.x*co, ab.y*si );
	
    return length(closestPoint - p ) * sign(p.y-closestPoint.y);
}


An interactive version (mouse mouse over ellipse to deform it) and reference code (click in the title of the viewport) can be found here: