website articles
distance to an ellipse - 2009
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 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

s2(ω) = |e(ω)-z|2

which expands to

s2(ω) = |e(ω)z|2 + |z|2 - 2<e(ω),z> = a2⋅cos2ω + b2⋅sin2ω + x2 + y2 - 2ax⋅cos ω - 2by⋅sin ω

The closest point will be given by



which expands in this case to

(b2-a2)⋅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 k4 (the only condition to be able to do so is that the ellipse is not, in fact, a circle), and by then renaming ,

k4 = 1
k3 = 2m
k2 = m2 + n2 - 1
k1 = -2m
k0 = -m2



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

z3 - 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: