### The maths

So, for a given point in the plane z={x,y} and a parametric ellipse

*e*centered in the origin with radii

*a*and

*b*in the

*x*and

*y*axes e(ω) = { a⋅cos ω, b⋅sin ω }, the distance from

*z*to the ellipse is

s

^{2}(ω) = |e(ω)-z|

^{2}

which expands to

s

^{2}(ω) = |e(ω)z|

^{2}+ |z|

^{2}- 2<e(ω),z> = a

^{2}⋅cos

^{2}ω + b

^{2}⋅sin

^{2}ω + x

^{2}+ y

^{2}- 2ax⋅cos ω - 2by⋅sin ω

The closest point will be given by

which expands in this case to

(b

^{2}-a

^{2})⋅sin ω⋅cos ω + ax⋅sin ω - by⋅cos ω = 0

This calls for the change of variable λ=cos(ω), so that

which expands into a quartic equation with the following coefficients:

or by dividing the equation by the leading coefficient k

_{4}(the only condition to be able to do so is that the ellipse is not, in fact, a circle), and by then renaming ,

k

_{4}= 1

k

_{3}= 2m

k

_{2}= m

^{2}+ n

^{2}- 1

k

_{1}= -2m

k

_{0}= -m

^{2}

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

z

^{3}- 3Qz + 2R = 0

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: