**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 resole 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 looses 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 measuting the distance from a point outside the elliposide), 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 interatcive version (mouse mouse over ellipse to deform it) and referecne code (click in the title of the viewport) can be found here: