Distance to an ellipse
A wonderful dive into the maths of ellipses.
Finding the distance between a point and an ellipse involves solving a quartic equation. The ellipse can intersect an ellipse at most four times.
The problem is to find the circle that intersects the ellipse once.
An iterative method is relatively easy to implement, lets see.
An oversized circle will intersect the ellipse at least two times. The iteration to minimise the distance between point a on the ellipse and p is as follows:
Choose a point a on the ellipse. Compute the distance r between a and p . On the circle centred at p with radius r, find the other intersection b. Set a’ to the midpoint of a and b on the ellipse.
Visually a’ is obviously a better estimate than a. This process is still tricky to compute but we can approximate to improve the situation.
Consider the previous algorithm when the ellipse is a circle, a’ then trivially lies on the line between c and p. It can also be seen that a’ is optimal. The idea is to approximate the ellipse as a circle, but first we must consider curvature.
Every point on the ellipse has a perpendicular normal. Also every point has a radius of curvature. A tighter curve will result in a smaller radius of curvature, while a slight curve will yield a larger radius of curvature. The centre of curvature of some point e on the ellipse is the centre of the circle with radius equal to the radius of curvature and perpendicular to the ellipse at e ( 3 of these poonts are shown below). In other words, it is the centre of the circle that locally approximates the ellipse at e. The centre of curvatures of all points can be plotted to form the evolute of the ellipse, shown as the blue line below.
Mathematically the parametric definition is best here, so for the ellipse we have,
x = a * cos(t)
y = b * sin(t)
The ellipse evolute also has a parametric form.
ev_x = (aa-bb)/a * cos^3(t)
ev_y = (bb-aa)/b * sin^3(t)
This is useful for approximating the ellipse as a circle, since we know the centre of curvature for any t.
Now when the ellipse is a circle, a’ lies on the line between the circle centre c and p. So here a’ lies on the line between the centre of curvature c and p.
Intersecting a line with an ellipse is a much simpler problem. Below we can see how to define two vectors r and q. Very loosely speaking, r is the radius of curvature and q intersects the ellipse where we would like the radius of curvature to have been. Since we using our circle approximation, we can compute the arc length between a and a’ as the angle between r and q times the radius of curvature.
\[ \large sin(\frac{Δc}{|r|}) = \frac{(r \times q)}{|r||q|} \]
Since Δc is small, we can drop the sine.
Now we need to determine how much to vary t by to achieve the same arc length Δc on the ellipse. With a bit of calculus:
dx/dt = -a * sin(t)
dy/dt = b * cos(t)
dc/dt = sqrt( (dx/dt)^2 + (dy/dt)^2 )
Δc/Δt ≈ sqrt(a^2 * sin^2(t) + b^2 * cos^2(t))
Δt ≈ Δc / sqrt(a^2 * sin^2(t) + b^2 * cos^2(t))
Δt ≈ Δc / sqrt(a^2 - a^2 * cos^2(t) + b^2 - b^2 * sin^2(t))
Δt ≈ Δc / sqrt(a^2 + b^2 - x^2 - y^2)
a’x = a * cos(t + Δt)
a’y = b * sin(t + Δt)
def solve(semi_major, semi_minor, p):
px = abs(p[0])
py = abs(p[1])
t = math.pi / 4 # a guess
a = semi_major
b = semi_minor
for x in range(0, 3):
x = a * math.cos(t)
y = b * math.sin(t)
ex = (a*a - b*b) * math.cos(t)**3 / a
ey = (b*b - a*a) * math.sin(t)**3 / b
rx = x - ex
ry = y - ey
qx = px - ex
qy = py - ey
r = math.hypot(ry, rx)
q = math.hypot(qy, qx)
delta_c = r * math.asin((rx*qy - ry*qx)/(r*q))
delta_t = delta_c / math.sqrt(a*a + b*b - x*x - y*y)
t += delta_t
t = min(math.pi/2, max(0, t))
return (math.copysign(x, p[0]), math.copysign(y, p[1]))Three iterations of this produces excellent results.
Why is this useful ?
Well, when I generate dummy circle dose profiles , initially I had my gamma ΔD and Δx being the same. But what if we need to determine gamma for 3 cGy and 1 mm. Rather than changing the gamma criteria we rescale (by their relative ratios) the distribution from our evaluation point, so that the profile changes from a circle to an ellipse and then we can use this function to determine the individual gamma indices.