This post provides an approach for calculating the minimum distance between two ellipses using gradient descent to optimize the square distance equation between points on the ellipses.

Ellipse Diagram

Figure 1: Ellipse Geometry

Equation Development

The ellipse geometry is shown in Figure 1. The parametric equations for the local coordinates of an ellipse, i.e. the coordinates of the ellipse centered at (0, 0) with no rotation, are

\[x^* = (a \cos t, b \cos t)\]

From this, the global coordinates for the ellipse are

\[x = r x^* + x_0\]

where the rotation matrix is

\[r = \begin{bmatrix} \cos \phi & -\sin \phi \\ \sin \phi & \cos \phi \end{bmatrix}\]

To find the minimum distance, we can optimize the square distance between the ellipse coordinates, $ x_1 $ and $ x_2 $:

\[D^2 = (x_1 - x_2) \cdot (x_1 - x_2)\]

Gradient descent may be used to find the minima of this function. For this, the derivatives with respect to the parameters $ t_1 $ and $ t_2 $ are required. The derivatives are:

\[\frac{d D^2}{d t_1} = 2 \frac{d x_1}{d t_1} \cdot (x_1 - x_2)\] \[\frac{d D^2}{d t_2} = -2 \frac{d x_2}{d t_2} \cdot (x_1 - x_2)\] \[\frac{d x_1}{d t_1} = r_1 \begin{bmatrix} -a_1 \sin t_1 \\ b_1 \cos t_1 \end{bmatrix}\] \[\frac{d x_2}{d t_2} = r_2 \begin{bmatrix} -a_2 \sin t_2 \\ b_2 \cos t_2 \end{bmatrix}\]

Finally, for the gradient descent method, iterations on $ t_1 $ and $ t_2 $ may be taken as

\[\begin{bmatrix} t_1 \\ t_2 \end{bmatrix}_{n+1} = \begin{bmatrix} t_1 \\ t_2 \end{bmatrix}_n - \alpha \begin{bmatrix} \frac{d D^2}{d t_1} \\ \frac{d D^2}{d t_2} \end{bmatrix}\]

For ellipses whose centers are not contained within each other, reasonable initial values for $ t_1 $ and $ t_2 $ are

\[t_1 = \text{atan2}\bigg(\frac{y_2 - y_1}{x_2 - x_1}\bigg) - \phi_1\] \[t_2 = \text{atan2}\bigg(\frac{y_1 - y_2}{x_1 - x_2}\bigg) - \phi_2\]

For ellipses whose centers are contained within each other, however, this value is not adequate to prevent the solution from converging to another local solution. Generating polygon vertices around the ellipses and using the angles corresponding to the minimum distance between the ellipse vertices resolves this issue, though this solution may be somewhat inelegant.

Import Statements

The following imports the Python packages that will be used for implementing the algorithms. This post uses the numpy package for array operations and matplotlib package for plotting results.

%matplotlib notebook
import math
import numpy as np
import matplotlib.pyplot as plt

Ellipse Class Implementation

An ellipse class implementing the developed equations into the find_distance method is shown below.

class Ellipse:
    def __init__(self, x, y, width, height, angle = 0):
        self.x = x
        self.y = y
        self.width = width
        self.height = height
        self.angle = angle
        
    def rotation_matrix(self):
        """
        Returns the rotation matrix for the ellipse's rotation.
        """
        a = math.cos(self.angle)
        b = math.sin(self.angle)
        return np.array([[a, -b], [b, a]])
    
    def get_point(self, angle):
        """
        Returns the point on the ellipse at the specified local angle.
        """
        r = self.rotation_matrix()
        xe = 0.5 * self.width * math.cos(angle)
        ye = 0.5 * self.height * math.sin(angle)
        return np.dot(r, [xe, ye]) + [self.x, self.y]
    
    def get_points(self, count):
        """
        Returns an array of points around the ellipse in the specified count.
        """
        t = np.linspace(0, 2 * math.pi, count)
        xe = 0.5 * self.width * np.cos(t)
        ye = 0.5 * self.height * np.sin(t)
        r = self.rotation_matrix()
        return np.dot(np.column_stack([xe, ye]), r.T) + [self.x, self.y]
    
    def contains(self, x):
        """
        Returns true if the point is contained in the ellipse.
        """
        x = np.asarray(x)
        a = 0.5 * self.width
        b = 0.5 * self.height
        r = self.rotation_matrix()
        xe = np.dot(r.T, x)
        return (xe[0] / a)**2 + (xe[1] / b)**2 <= 1
    
    def find_distance(self, other, tolerance = 1e-8, max_iterations = 10000, learning_rate = 0.04, decay_rate = 0.002):
        """
        Finds the minimum distance to another ellipse using gradient descent.
        """
        r1 = self.rotation_matrix()
        r2 = other.rotation_matrix()
        dx0 = np.array([self.x - other.x, self.y - other.y])
        t1 = math.atan2(-dx0[1], -dx0[0]) - self.angle
        t2 = math.atan2(dx0[1], dx0[0]) - other.angle
        t = np.array([t1, t2])
        a1 = 0.5 * self.width
        b1 = 0.5 * self.height
        a2 = 0.5 * other.width
        b2 = 0.5 * other.height
        iterations = 0
        error = tolerance
        errors = []
        
        # If ellipse center is contained in another ellipse, use vertices to determine initial angles
        if self.contains([other.x, other.y]) or other.contains([self.x, self.y]):
            ts = np.linspace(0, 2 * math.pi, 24, endpoint = False)
            ts = np.array(np.meshgrid(ts, ts)).T.reshape(-1, 2)
            xe = a1 * np.cos(ts[:,0])
            ye = b1 * np.sin(ts[:,0])
            x1 = np.dot(np.column_stack([xe, ye]), r1.T) + [self.x, self.y]
            xe = a2 * np.cos(ts[:,1])
            ye = b2 * np.sin(ts[:,1])
            x2 = np.dot(np.column_stack([xe, ye]), r2.T) + [other.x, other.y]
            t = ts[np.argmin(np.linalg.norm(x1 - x2, axis = 1))]
        
        while error >= tolerance and iterations < max_iterations:
            cost = np.cos(t)
            sint = np.sin(t)
            x1 = np.dot(r1, [a1 * cost[0], b1 * sint[0]])
            x2 = np.dot(r2, [a2 * cost[1], b2 * sint[1]])
            dx1 = np.dot(r1, [-a1 * sint[0], b1 * cost[0]])
            dx2 = np.dot(r2, [-a2 * sint[1], b2 * cost[1]])
            delta = x1 - x2 + dx0
            dp = np.array([2 * np.dot(dx1, delta), -2 * np.dot(dx2, delta)])
            alpha = learning_rate * math.exp(-decay_rate * iterations)
            t -= alpha * dp
            error = np.max(np.abs(dp))
            errors.append(error)
            iterations += 1
            
        errors = np.array(errors)
        y = np.linalg.norm(x1 - x2 + dx0)
        success = error < tolerance and iterations < max_iterations
        return dict(x = t, y = y, error = error, iterations = iterations, success = success, errors = errors)

Method Evaluation

In order to evaluate the method, the following creates an ellipse and several sample ellipses to which distances may be evaluated. As shown, the method successfully converges for all samples.

ellipse = Ellipse(0, 0, 6, 10, np.deg2rad(15))

ellipses = [
    Ellipse(0, 0, 6, 10, np.deg2rad(15)),
    Ellipse(0, 0, 6, 10, np.deg2rad(195)),
    Ellipse(0, 0, 6, 9, np.deg2rad(15)),
    Ellipse(0, 0, 6, 9, np.deg2rad(195)),
    Ellipse(0, 0, 5, 10, np.deg2rad(15)),
    Ellipse(0, 0, 5, 10, np.deg2rad(195)),
    Ellipse(0, 0, 6, 10, np.deg2rad(-15)),
    Ellipse(0, 0, 6, 10, np.deg2rad(165)),
    Ellipse(7, 7, 5, 9, np.deg2rad(15)),
    Ellipse(7, 7, 5, 9, np.deg2rad(195)),
    Ellipse(-7, -7, 7, 4, np.deg2rad(0)),
    Ellipse(-7, -7, 7, 4, np.deg2rad(180)),
    Ellipse(-7, 7, 10, 4, np.deg2rad(0)),
    Ellipse(-7, 7, 10, 4, np.deg2rad(180)),
    Ellipse(7, -7, 4, 6, np.deg2rad(45)),
    Ellipse(7, -7, 4, 6, np.deg2rad(225)),
    Ellipse(-8, 0, 4, 4, np.deg2rad(0)),
    Ellipse(-8, 0, 4, 4, np.deg2rad(15)),
    Ellipse(0.5, -0.5, 2, 4, np.deg2rad(0)),
]

solutions = [ellipse.find_distance(x) for x in ellipses]
print(f"Method Successful: {np.all([x['success'] for x in solutions])}")

Method Successful: True

Method Convergence

A plot of the method convergence is shown below. For the selected learning schedule and parameters, the method converges relatively quickly for most of the samples, with some noise demonstrated in the early iterations of a few samples, likely due to overshooting. An adaptive learning rate might fair better but is beyond the scope of this post.

fig = plt.figure()
ax = fig.add_subplot(111, title = "Error vs. Iteration")

for solution in solutions:
    x = solution["errors"]
    ax.plot(np.arange(len(x)), x, "-")

Visual Verification

The following code plots the results for visual verification. The shortest distance vectors, shown in black, appear correct for all of the ellipses in the sample.

fig = plt.figure()
ax = fig.add_subplot(111, title = "Ellipse Distances", aspect = "equal")
ellipse_points = ellipse.get_points(100)
ax.plot(ellipse_points[:,0], ellipse_points[:,1], "k-", linewidth = 5)

for solution, other in zip(solutions, ellipses):
    x1 = other.get_points(100)
    x2 = np.array([ellipse.get_point(solution["x"][0]), other.get_point(solution["x"][1])])
    ax.plot(x1[:,0], x1[:,1], "-")
    ax.plot(x2[:,0], x2[:,1], "k-")
    ax.plot(x2[:,0], x2[:,1], "k.")

Conclusion

From the above results, the following conclusions can be drawn:

  • The developed equations and gradient descent method were successful is finding the minimum distance between the sample ellipses.
  • An adaptive learning rate might improve the convergence rate of the method.

Edits

  • 3/1/24 - Added missing dx0 term to returned distance calculation per issue discovered by wergil81.

Read Next