# Distance from a Point to an Ellipse, an Ellipsoid, or a Hyperellipsoid

David Eberly

<figure><img src="image_1.png" /></figure>

## To cite this version:

David Eberly. Distance from a Point to an Ellipse, an Ellipsoid, or a Hyperellipsoid. Geometric Tools. 2020. [hal-04946422]

HAL Id: hal-04946422

https://hal.science/hal-04946422v1

Submitted on 13 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Distance from a Point to an Ellipse, an Ellipsoid, or a Hyperellipsoid

David Eberly, Geometric Tools, Redmond WA 98052
https://www.geometrictools.com/

This work is licensed under the Creative Commons Attribution 4.0 International License. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/ or send a letter to Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.

Created: June 28, 2013

Last Modified: September 11, 2020

## Contents

1 Introduction 3

2 Distance from a Point to an Ellipse 3

2.1 The Closest Point's Normal is Directed Toward the Query Point 3

2.2 The Case of a Circle 4

2.3 The Query Point is the Origin 4

2.4 The Query Point is on the Vertical Axis 4

2.5 The Query Point is on the Horizontal Axis 5

2.6 The Query Point is Strictly in the First Quadrant 5

2.7 A Summary of the Mathematical Algorithm 6

2.8 Robust Root Finders 7

2.8.1 Bisection Method 7

2.8.2 Newton's Method 8

2.8.3 Conversion to a Polynomial Equation 10

2.9 A Robust Implementation 11

3 Distance from a Point to an Ellipsoid 13

3.1 The Closest Point's Normal is Directed Toward the Query Point 13

3.2 The Case of a Sphere 14

3.3 The Case of an Oblate Spheroid 14

3.4 The Case of a Prolate Spheroid 14
<table><tr><td>3.5</td><td>The Query Point is the Origin</td><td>14</td></tr><tr><td>3.6</td><td>The Query Point has y<sub>2</sub> &gt; 0</td><td>14</td></tr><tr><td>3.7</td><td>The Query Point has y<sub>2</sub> = 0</td><td>16</td></tr><tr><td>3.8</td><td>A Summary of the Mathematical Algorithm</td><td>16</td></tr><tr><td>3.9</td><td>Robust Root Finding and Conversion to a Polynomial</td><td>18</td></tr><tr><td>3.10</td><td>A Robust Implementation</td><td>18</td></tr><tr><td>4</td><td>Distance from a Point to a Hyperellipsoid</td><td>20</td></tr><tr><td>4.1</td><td>The Closest Point's Normal is Directed Toward the Query Point</td><td>20</td></tr><tr><td>4.2</td><td>The Key Observations</td><td>21</td></tr><tr><td>4.3</td><td>A Robust Implementation</td><td>22</td></tr></table>
# 1 Introduction

This document describes an algorithm for computing the distance from a point to an ellipse (2D), from a point to an ellipsoid (3D), and from a point to a hyperellipsoid (any dimension).

# 2 Distance from a Point to an Ellipse

A general ellipse in 2D is represented by a center point C, an orthonormal set of axis-direction vectors {U₀, U₁}, and associated extents eᵢ with e₀ ≥ e₁ > 0. The ellipse points are

$$ \mathbf{P} = \mathbf{C} + x_0 \mathbf{U}_0 + x_1 \mathbf{U}_1 \quad (1) $$

where

$$ \left(\frac{x_0}{e_0}\right)^2 + \left(\frac{x_1}{e_1}\right)^2 = 1 \quad (2) $$

If e₀ = e₁, then the ellipse is a circle with center C and radius e₀. The orthonormality of the axis directions and Equation (1) imply $x_i = U_i \cdot (\mathbf{P} - \mathbf{C})$. Substituting this into Equation (2) we obtain

$$ (\mathbf{P} - \mathbf{C})^{\mathrm{T}} M (\mathbf{P} - \mathbf{C}) = 1 \quad (3) $$

where $M = RDR^{\mathrm{T}}$, R is an orthogonal matrix whose columns are U₀ and U₁, and D is a diagonal matrix whose diagonal entries are 1/e₀² and 1/e₁².

The problem is to compute the distance from a point Q to the ellipse. It is sufficient to solve this problem in the coordinate system of the ellipse; that is, represent $Q = C + y_0U_0 + y_1U_1$. The distance from Q to the closest point P on the ellipse as defined by Equation (3) is the same as the distance from $Y = (y_0, y_1)$ to the closest point $X = (x_0, x_1)$ on the standard ellipse of Equation (2).

We may additionally use symmetry to simplify the construction. It is sufficient to consider the case when $(y_0, y_1)$ is in the first quadrant: $y_0 \ge 0$ and $y_1 \ge 0$. For example, if $(y_0, y_1)$ is in the second quadrant where $y_0 < 0$, and if $(x_0, x_1)$ is the closest point which must be in the second quadrant, then $(-y_0, y_1)$ is in the first quadrant and $(-x_0, x_1)$ is the closest ellipse point which must be in the first quadrant. If we transform the query point into the first quadrant by sign changes on the components, we can construct the closest point in the first quadrant, and then undo the sign changes to obtain the result in the original quadrant.

## 2.1 The Closest Point's Normal is Directed Toward the Query Point

A parameterization of the standard ellipse is $X(\theta) = (e_0 \cos \theta, e_1 \sin \theta)$ for $\theta \in [0, 2\pi)$. The squared distance from Y to any point on the ellipse is

$$ F(\theta) = |X(\theta) - Y|^2 \quad (4) $$

This is a nonnegative, periodic, and differentiable function; it must have a global minimum occurring at an angle for which the first-order derivative is zero,

$$ F'(\theta) = 2(X(\theta) - Y) \cdot X'(\theta) = 0 \quad (5) $$

For the derivative to be zero, the vectors $(X(\theta) - Y)$ and $X'(\theta)$ must be perpendicular. The vector $X'(\theta)$ is tangent to the ellipse at $X(\theta)$. This implies that the vector from Y to the closest ellipse point X is normal to
the curve at X. Using the implicit form of the ellipse, namely, $G(x_0, x_1) = (x_0/e_0)^2 + (x_1/e_1)^2 - 1 = 0$, half the gradient of $G(x_0, x_1)$ is a normal vector to the ellipse at $(x_0, x_1)$, so $(y_0, y_1) - (x_0, x_1) = t \nabla G(x_0, x_1)/2 = t(x_0/e_0^2, x_1/e_1^2)$ for some scalar $t$, or

$$y_0 = x_0 \left( 1 + \frac{t}{e_0^2} \right), \quad y_1 = x_1 \left( 1 + \frac{t}{e_1^2} \right) \qquad (6)$$

If $(y_0, y_1)$ is outside the ellipse, it is necessary that $t > 0$. If $(y_0, y_1)$ is inside the ellipse, it is necessary that $t < 0$. If $(y_0, y_1)$ is already on the ellipse, then $t = 0$ and the distance is zero.

## 2.2 The Case of a Circle

If $e_0 = e_1$, then the ellipse is a circle. The origin $(0,0)$ has infinitely many closest circle points (all of them), but clearly the distance from the origin to the circle is $e_0$. The closest circle point to a point $(y_0, y_1) \neq (0,0)$ is $(x_0, x_1) = e_0(y_0, y_1)/|(y_0, y_1)|$. Equation (6) is consistent with this, because it implies $(x_0, x_1) = (e_0^2/(t + e_0^2))(y_0, y_1)$ for some $t$; that is, $(x_0, x_1)$ is parallel to $(y_0, y_1)$ and must have length $e_0$. It is easily shown that $t = -e_0^2 + e_0\sqrt{y_0^2 + y_1^2}$.

For the remainder of Section 2, we assume $e_0 > e_1$ and that all analysis is restricted to the first quadrant.

## 2.3 The Query Point is the Origin

Let $y_0 = 0$ and $y_1 = 0$. Equation (6) becomes $0 = x_0(1 + t/e_0^2)$ and $0 = x_1(1 + t/e_1^2)$. We have four cases to consider.

* If $x_0 = 0$ and $x_1 = 0$, the point is not on the ellipse; this case may be discarded.

* If $x_0 = 0$ and $x_1 \neq 0$, then $t = -e_1^2$ and the only constraint on $x_1$ is that $(x_0, x_1)$ be a point on the ellipse, which means $x_1 = e_1$. The candidate closest point is $(0, e_1)$.

* If $x_0 \neq 0$ and $x_1 = 0$, then $t = -e_0^2$ and the only constraint on $x_0$ is that $(x_0, x_1)$ be a point on the ellipse, which means $x_0 = e_0$. The candidate closest point is $(e_0, 0)$.

* If $x_0 \neq 0$ and $x_1 \neq 0$, then $t = -e_0^2$ and $t = -e_1^2$, which is a contradiction because $e_0 \neq e_1$; this case may be discarded.

The only candidate ellipse points in the first quadrant closest to $(0,0)$ are $(e_0, 0)$ and $(0, e_1)$. Of these two, $(0, e_1)$ is closer. In summary: The closest ellipse point in the first quadrant to $(0,0)$ is the point $(0, e_1)$ with distance $d = e_1$.

## 2.4 The Query Point is on the Vertical Axis

Let $y_0 = 0$ and $y_1 > 0$. Equation (6) becomes $0 = x_0(1 + t/e_0^2)$ and $y_1 = x_1(1 + t/e_1^2)$. We have two cases to consider.

* If $x_0 = 0$, then for $(x_0, x_1)$ to be on the ellipse, we need $x_1 = e_1$. The candidate closest point is $(0, e_1)$.
* If $x_0 \neq 0$, then $t = -e_0^2$ and $y_1 = x_1(1 - e_0^2/e_1^2) \le 0$, where the last inequality is a consequence of $x_1 \ge 0$ and $e_0 > e_1$. This contradicts the assumption that $y_1 > 0$.

The only candidate ellipse point is $(0, e_1)$. In summary: The closest in the first quadrant to $(0, y_1)$ for $y_1 > 0$ is $(0, e_1)$ with distance $d = |y_1 - e_1|$.

## 2.5 The Query Point is on the Horizontal Axis

Let $y_0 > 0$ and $y_1 = 0$. Equation (6) becomes $y_0 = x_0(1 + t/e_0^2)$ and $0 = x_1(1 + t/e_1^2)$. We have two subcases to consider.

* If $x_1 = 0$, then for $(x_0, x_1)$ to be on the ellipse, we need $x_0 = e_0$. The candidate closest point is $(e_0, 0)$. The squared distance to this candidate is

$$d_0^2 = (y_0 - e_0)^2 \qquad (7)$$

* If $x_1 \neq 0$, then $t = -e_1^2$ and $y_0 = x_0(1 - e_1^2/e_0^2)$. It follows that $x_0 = e_0^2 y_0 / (e_0^2 - e_1^2) \le e_0$. The last inequality is a consequence of $(x_0, x_1)$ being on the ellipse. The implication is that $y_0 \le (e_0^2 - e_1^2)/e_0 < e_0$. Notice that this limits the construction to points $(y_0, 0)$ inside the ellipse. For points $(y_0, 0)$ with $y_0 \ge (e_0^2 - e_1^2)/e_0$, which includes points inside and outside the ellipse, the closest point is necessarily $(e_0, 0)$. When $y_0 \le (e_0^2 - e_1^2)/e_0$, for $(x_0, x_1)$ to be on the ellipse, we may solve for $x_1 = e_1 \sqrt{1 - (x_0/e_0)^2}$. The squared distance from the candidate $(x_0, x_1)$ to $(y_0, 0)$ is

$$d_1^2 = (x_0 - y_0)^2 + x_1^2 = e_1^2 \left( 1 - \frac{y_0^2}{e_0^2 - e_1^2} \right) \qquad (8)$$

The two candidates are $(e_0, 0)$ and $(x_0, x_1)$. We need to determine which of these is closer to $(y_0, 0)$. It may be shown that

$$d_0^2 - d_1^2 = \frac{(e_0 y_0 + e_1^2 - e_0^2)^2}{e_0^2 - e_1^2} \ge 0 \qquad (9)$$

This implies that $(x_0, x_1)$ is the closer point. In summary: The closest point in the first quadrant to $(y_0, 0)$ for $0 < y_0 < (e_0^2 - e_1^2)/e_0$ is $(x_0, x_1)$ with $x_0 = e_0^2 y_0 / (e_0^2 - e_1^2)$, $x_1 = e_1 \sqrt{1 - (x_0/e_0)^2}$, and distance $d = \sqrt{(x_0 - y_0)^2 + x_1^2} = e_1 \sqrt{1 - y_0^2 / (e_0^2 - e_1^2)}$. The closest point to $(y_0, 0)$ for $y_0 \ge e_0 - e_1^2/e_0$ is $(e_0, 0)$ with distance $d = |y_0 - e_0|$.

## 2.6 The Query Point is Strictly in the First Quadrant

Let $y_0 > 0$ and $y_1 > 0$. Equation (6) has solution

$$x_0 = \frac{e_0^2 y_0}{t + e_0^2}, \quad x_1 = \frac{e_1^2 y_1}{t + e_1^2} \qquad (10)$$

for some scalar $t$. We know that the closest point in the first quadrant requires $x_0 \ge 0$ and $x_1 \ge 0$, which implies $t > -e_0^2$ and $t > -e_1^2$. Because $e_0 > e_1$, it is enough to analyze only those $t$-values for which $t > -e_1^2$.
Substitute Equation (10) in Equation (2) to obtain

$$F(t) = \left( \frac{e_0 y_0}{t + e_0^2} \right)^2 + \left( \frac{e_1 y_1}{t + e_1^2} \right)^2 - 1 = 0 \qquad (11)$$

The first equality defines the function $F(t)$ with domain $(-e_1^2, \infty)$. The candidates for the ellipse point $(x_0, x_1)$ closest to $(y_0, y_1)$ are generated by the roots $t$ to $F(t) = 0$.

The first-order and second-order derivatives of $F$ are

$$F'(t) = -\frac{2e_0^2 y_0^2}{(t + e_0^2)^3} - \frac{2e_1^2 y_1^2}{(t + e_1^2)^3}, \quad F''(t) = \frac{6e_0^2 y_0^2}{(t + e_0^2)^4} + \frac{6e_1^2 y_1^2}{(t + e_1^2)^4} \qquad (12)$$

For $t > -e_1^2$ we have the conditions $F'(t) < 0$ and $F''(t) > 0$. Also observe that

$$\lim_{t \to -e_1^2^+} F(t) = +\infty, \quad \lim_{t \to \infty} F(t) = -1 \qquad (13)$$

The first expression is a one-sided limit where $t$ approaches $e_1^2$ through values larger than $e_1^2$. We have shown that $F(t)$ is a strictly decreasing function for $t \in (-e_1^2, +\infty)$ that is initially positive, then becomes negative. Consequently it has a unique root on the specified domain. Figure 1 shows a typical graph of $F(t)$.

<figure><img src="image_3.png" /><figcaption>Figure 1. A typical graph of $F(t)$ for $t > -e_1^2$. The unique root $\bar{t}$ is shown.</figcaption></figure>

The domain $(-e_1^2, \infty)$ contains 0 and $F(0) = (y_0/e_0)^2 + (y_1/e_1)^2 - 1$. When $(y_0, y_1)$ is inside the ellipse we see that $F(0) < 0$ and $\bar{t} < 0$. When $(y_0, y_1)$ is outside the ellipse we see that $F(0) > 0$ and $\bar{t} > 0$. When $(y_0, y_1)$ is on the ellipse we see that $F(0) = 0$ and $\bar{t} = 0$.

## 2.7 A Summary of the Mathematical Algorithm

Listing 1 summarizes the algorithm for computing the ellipse point $(x_0, x_1)$ closest to the query point $(y_0, y_1)$ in the first quadrant and for computing the distance to the closest point. The preconditions are $e_0 \ge e_1 > 0$, $y_0 \ge 0$, and $y_1 \ge 0$.

**Listing 1.** Pseudocode for computing the closest ellipse point and distance to a query point. It is required that $e_0 \ge e_1 > 0$, $y_0 \ge 0$, and $y_1 \ge 0$.
```
Real DistancePointEllipse(Real e0, Real e1, Real y0, Real y1, Real& x0, Real& x1)
{
    Real distance;
    if (y1 > 0)
    {
        if (y0 > 0)
        {
            Compute the unique root tbar of F(t) on (-e1*e1, +infinity);
            x0 = e0*e0*y0/(tbar + e0*e0);
            x1 = e1*e1*y1/(tbar + e1*e1);
            distance = sqrt((x0 - y0)*(x0 - y0) + (x1 - y1)*(x1 - y1));
        }
        else // y0 == 0
        {
            x0 = 0;
            x1 = e1;
            distance = abs(y1 - e1);
        }
    }
    else // y1 == 0
    {
        if (y0 < (e0*e0 - e1*e1) / e0)
        {
            x0 = e0*e0*y0/(e0*e0 - e1*e1);
            x1 = e1*sqrt(1 - (x0/e0)*(x0/e0));
            distance = sqrt((x0 - y0)*(x0 - y0) + x1*x1);
        }
        else
        {
            x0 = e0;
            x1 = 0;
            distance = abs(y0 - e0);
        }
    }
    return distance;
}
```

## 2.8 Robust Root Finders

Various methods are discussed for computing the unique root $\bar{t} \in (-e_1^2, +\infty)$ of $F(t)$. The first uses bisection and is the most robust. The second uses Newton's method, but is has numerical problems when $y_1$ is nearly zero. The third uses a conversion to a polynomial equation, but this approach has its own problems. The conclusion is that bisection is the safest and most robust method to use.

### 2.8.1 Bisection Method

You may locate the unique root $\bar{t} \in (-e_1^2, +\infty)$ of $F(t)$ using bisection. A finite-length bounding interval $[t_0, t_1]$ that contains $\bar{t}$ is required. The minimum of the bounding interval may be chosen as

$$t_0 = -e_1^2 + e_1 y_1 > -e_1^2 \qquad (14)$$

The function value at that parameter is

$$F(t_0) = \left( \frac{e_0 y_0}{t_0 + e_0^2} \right)^2 > 0 \qquad (15)$$

Because $t + e_0^2 > t + e_1^2$, we can see that

$$F(t) = \left(\frac{e_0 y_0}{t + e_0^2}\right)^2 + \left(\frac{e_1 y_1}{t + e_1^2}\right)^2 - 1 < \left(\frac{e_0 y_0}{t + e_0^2}\right)^2 + \left(\frac{e_1 y_1}{t + e_1^2}\right)^2 - 1 = \frac{e_0^2 y_0^2 + e_1^2 y_1^2}{(t + e_1^2)^2} - 1 \qquad (16)$$
The right-hand side of this equation is negative for $t \ge t_1$, where

$$t_1 = -e_1^2 + \sqrt{e_0^2 y_0^2 + e_1^2 y_1^2} \qquad (17)$$

that is, $F(t) < 0$ for $t \ge t_1$. In summary, a bounding interval for the root $\bar{t}$ is $[t_0, t_1]$, where $t_0$ is defined by Equation (14) and $t_1$ is defined by Equation (17).

You may bisect until you reach your own stopping criterion, say, $|t_1 - t_0| < \varepsilon$ for a specified small number $\varepsilon > 0$. Our implementation bisects until $F = 0$ exactly as a floating-point number or until the midpoint $(t_0 + t_1)/2$ is equal to an endpoint due to rounding. The maximum number of iterations is 149 for float (32-bit numbers) or 1074 for double (64-bit numbers). Generally, the maximum number is

$$\max\text{iterations} = \text{std::numeric\_limits\langle Real \rangle::digits} - \text{std::numeric\_limits\langle Real \rangle::min\_exponent}; \qquad (18)$$

where Real is a floating-point type.

### 2.8.2 Newton's Method

The condition $F''(t) > 0$ says that $F$ is a convex function. Such functions are ideally suited for the application of Newton's Method for root finding. Given an initial guess $t_0$, the iterates are

$$t_{n+1} = t_n - \frac{F(t_n)}{F'(t_n)}, \quad n \ge 0 \qquad (19)$$

It is important to choose an initial guess for which the method converges. Newton's method has an intuitive geometric appeal to it. The next value $t_{n+1}$ is computed by determining where the tangent line to the graph at $(t_n, F(t_n))$ intersects the t-axis. The intersection point is $(t_{n+1}, 0)$. If we choose an initial guess $t_0 < \bar{t}$, the tangent line to $(t_0, F(t_0))$ intersects the t-axis at $(t_1, 0)$ where $t_0 < t_1 < \bar{t}$. Figure 2 illustrates this.

<figure><img src="image_4.png" /><figcaption>Figure 2. An initial guess $t_0$ to the left of $\bar{t}$ guarantees $t_0 < t_1 < \bar{t}$.</figcaption></figure>

If we were instead to choose an initial guess $t_0 > \bar{t}$, the new iterate satisfies $t_1 < \bar{t}$, but potentially $t_1 < -e_1^2$ which is outside the domain of interest, namely $t \in (-e_1^2, \infty)$. Figure 3 illustrates this.

<figure><img src="image_5.png" /><figcaption>Figure 3. An initial guess $t_0$ to the right of $\bar{t}$ guarantees $t_1 < \bar{t}$ but does not guarantee $t_1 > -e_1^2$.</figcaption></figure>
To avoid this potential problem, it is better to choose an initial guess to the left of the root. Any $t_0 > -e_1^2$ for which $F(t_0) > 0$ will work. In particular, you can choose $t_0 = -e_1^2 + e_1y_1$, in which case $F(t_0) = [e_0y_0/(e_1y_1 + e_0^2 - e_1^2)]^2 > 0$.

The stopping criterion for Newton's Method typically involves testing the values $F(t_n)$ for closeness to zero. For robustness, testing is recommended to determine whether progress is made in the domain; that is, if the difference $|t_{n+1} - t_n|$ of consecutive iterates is sufficiently small, then the iteration terminates.

As mentioned previously, you need to be careful when using Newton's Method. At first it is mathematically comforting to have a convex function $F(t)$ for $t \in (-e_1^2, \infty)$, knowing that there is a unique root. However, there are potential problems when $y_1$ is nearly zero. When $y_0 > 0$ and $y_1 > 0$, Figure 4 shows typical graphs of $F(t)$.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Typical graphs of $F(t)$ when $y_0 > 0$ and $y_1 > 0$.</figcaption></figure>

The three possibilities are based on where $(y_0, y_1)$ is located with respect to the *evolute* of the ellipse. Define $G(y_0, y_1) = (e_0y_0)^{2/3} + (e_1y_1)^{2/3} - (e_0^2 - e_1^2)^{2/3}$. The evolute is the level set defined by $G(y_0, y_1) = 0$. The left image of Figure 4 corresponds to $(y_0, y_1)$ outside the evolute, so $G(y_0, y_1) > 0$. The middle image corresponds to $(y_0, y_1)$ on the evolute itself, so $G(y_0, y_1) = 0$. The right image corresponds to $(y_0, y_1)$ inside the evolute, so $G(y_0, y_1) < 0$. You can generate the evolute equation by eliminating $t$ from the simultaneous equations $F(t) = 0$ and $F'(t) = 0$.

Imagine choosing positive values of $y_1$ closer and closer to zero. The vertical asymptotes at $t = -e_0^2$ and $t = -e_1^2$ remain so, but the graph of $F(t)$ near the asymptote $t = -e_1^2$ changes: it begins to hug the vertical asymptote. In the limiting sense as $y_1$ approaches zero, the graph is illustrated by Figure 5.
**Figure 5.** The limiting graph of $F(t)$ when $y_0 > 0$ and $y_1 \to 0^+$. Notice that there is effectively a discontinuity at $t = -e_1^2$, and the largest root jumps from $\bar{t}$ to $\hat{t}$.

<figure><img src="image_7.png" /></figure>

Just setting $y_1 = 0$, we have $F(t) = (e_0y_0/(t + e_0^2))^2 - 1$ and the largest root of $F(t)$ is $\hat{t} \in (-e_0^2, -e_1^2)$. As a function of $y_1$, the largest root of $F(t)$ is discontinuous in the sense that

$$ -e_1^2 = \lim_{y_1 \to 0^+} \bar{t}(y_1) \neq \hat{t} = -e_0^2 + e_0y_0 \quad (20) $$

Thus, the graph of $F(t)$ has a topological change at $y_1 = 0$, which can cause numerical problems with Newton's Method. For the initial guess $t_0 = -e_1^2 + e_1y_1$ when $y_1 > 0$,

$$ F'(t_0) = - \frac{2e_0^2 y_0^2}{(e_1 y_1 + e_0^2 - e_1^2)^3} - \frac{2}{e_1 y_1} \quad (21) $$

which is a negative number of large magnitude when $y_1$ is nearly zero. The number can be so large that the floating-point representation is infinite. The next iterate $t_1 = t_0 - F(t_0)/F'(t_0)$ is computed in floating-point arithmetic to be $t_0$ because $F(t_0) > 0$ as a floating-point number and $F(t_0)/F'(t_0)$ is zero as a floating-point number. The algorithm makes no progress towards computing the root.

### 2.8.3 Conversion to a Polynomial Equation

An approach for ellipsoids—but applies to ellipses as well—is described in [1]. The roots of the following quartic polynomial may be computed, leading to candidates for the closest ellipse point,

$$ P(t) = (t + e_0^2)^2(t + e_1^2)^4 F(t) = e_0^2 y_0^2 (t + e_1^2)^2 + e_1^2 y_1^2 (t + e_0^2)^2 - (t + e_0^2)^2 (t + e_1^2) \quad (22) $$

When $y_1 > 0$, the roots of $P(t)$ are the same as those for $F(t)$, and the largest root of $P(t)$ is the same as the largest root $\bar{t}$ for $F(t)$. However, when $y_1 = 0$, $P(t)$ has a double root at $t = -e_1^2$ but $F(t)$ does not have a root at $t = -e_1^2$.

The Graphics Gems IV article mentions that the largest root of $P(t)$ is the one that corresponds to the closest point. This is not true when $y_1 = 0$, as mentioned in the previous paragraph. In this case, $P(t) = (t + e_1^2)^2[e_0^2 y_0^2 - (t + e_0^2)^2]$. The largest root is $\bar{t} = -e_1^2$ and another root is $\hat{t} = -e_0^2 + e_0y_0$, as shown in Figure 5. It is $\hat{t}$ that determines the closest ellipse point, not $\bar{t}$.
The article also mentions that Newton's Method may be used, starting with an initial guess greater than the largest root. When $(y_0, y_1)$ is inside the ellipse, then $\bar{t} < 0$ and so the initial guess is chosen to be $t_0 = 0$. When $(y_0, y_1)$ is outside the ellipse, a geometric and algebraic argument is used to construct an initial guess $t_0 = e_0\sqrt{y_0^2 + y_1^2}$. No analysis is provided for the graph of $P(t)$ at the largest root. In particular, there are no guarantees that $P(t)$ is convex at that root, so the convergence of the Newton iterates is not guaranteed. Even so, when $y_1 = 0$, the iterates might converge to the root $\bar{t} = -e_1^2$, which is not the correct root that determines the closest ellipse point.

Empirical evidence suggests that Newton's Method applied to $P(t)$ to find its largest root suffers from numerical problems when $y_0$ or $y_1$ is nearly zero. For example, if $y_0$ is positive and if $y_1$ is a very small positive number, then $\bar{t}$ is the largest root of $P(t)$. It is nearly a double root, which can cause problems with Newton's Method. Specifically, if a function $G(t)$ has a double root at $\bar{t}$, Newton's Method must be slightly modified to $t_{i+1} = t_i - 2G(t_i)/G'(t_i)$. Generally, if $G(t)$ has a root of multiplicity $m$, then Newton's Method should be $t_{i+1} = t_i - mG(t_i)/G'(t_i)$.

And finally, the initial guess suggested in the Graphics Gems IV article is also quite large. From the discussion of the bisection method, $t_1 = -e_1 + \sqrt{e_0^2 y_0^2 + e_1^2 y_1^2}$ is an estimate greater than the maximum root but smaller than $t_0 = e_0\sqrt{y_0^2 + y_1^2}$. To see this, this,

$$ (t_1 + e_1)^2 = e_0^2 y_0^2 + e_1^2 y_1^2 = t_0^2 - (e_0^2 - e_1^2)y_1^2 < t_0^2 \quad (23) $$

which implies $t_1 < t_1 + e_1 < t_0$.

## 2.9 A Robust Implementation

An implementation that uses the robust bisection method is described next. A change of variables helps keep the magnitude of the numbers to a moderate size, $s = t/e_1^2$. The query point components are also scaled by the ellipse extents to keep the inputs on the order of 1 when the query is executed for points near the ellipse, $z_i = y_i/e_i$. Define the ratio of squared extents $r_0 = e_0^2/e_1^2 > 1$. The function $F(t)$ is transformed to

$$ G(s) = \left( \frac{r_0 z_0}{s + r_0} \right)^2 + \left( \frac{z_1}{s + 1} \right)^2 - 1 \quad (24) $$

for $s \in (-1, +\infty)$. The root of $G(s) = 0$ is $\bar{s} = \bar{t}/e_1^2$ and it is in the interval $[s_0, s_1]$, where $s_0 = -1 + z_1$ and $s_1 = -1 + \sqrt{r_0^2 z_0^2 + z_1^2}$.

Listing 2 is the robust implementation of the algorithm. The maximum number of iterations is provided by Equation (18).
<caption>**Listing 2.** Pseudocode for robustly computing the closest ellipse point and distance to a query point. It is required that $e_0 \ge e_1 > 0$, $y_0 \ge 0$, and $y_1 \ge 0$.</caption>

```c
Real GetRoot(Real r0, Real z0, Real z1, Real g)
{
    Real n0 = r0*z0;
    Real s0 = z1 - 1, s1 = (g < 0 ? 0 : RobustLength(n0, z1) - 1);
    Real s = 0;
    for (int i = 0; i < maxIterations; ++i)
    {
        s = (s0 + s1) / 2;
        if (s == s0 || s == s1) { break; }
        Real ratio0 = n0/(s + r0), ratio1 = z1/(s + 1);
        g = Sqr(ratio0) + Sqr(ratio1) - 1;
        if (g > 0) { s0 = s; } else if (g < 0) { s1 = s; } else { break; }
    }
    return s;
}

Real DistancePointEllipse(Real e0, Real e1, Real y0, Real y1, Real& x0, Real& x1)
{
    Real distance;
    if (y1 > 0)
    {
        if (y0 > 0)
        {
            Real z0 = y0/e0, z1 = y1/e1, g = Sqr(z0) + Sqr(z1) - 1;
            if (g != 0)
            {
                Real r0 = Sqr(e0/e1), sbar = GetRoot(r0, z0, z1, g);
                x0 = r0*y0/(sbar + r0); x1 = y1/(sbar + 1);
                distance = sqrt(Sqr(x0 - y0) + Sqr(x1 - y1));
            }
            else
            {
                x0 = y0; x1 = y1; distance = 0;
            }
        }
        else // y0 == 0
        {
            x0 = 0; x1 = e1; distance = abs(y1 - e1);
        }
    }
    else // y1 == 0
    {
        Real numer0 = e0*y0, denom0 = Sqr(e0) - Sqr(e1);
        if (numer0 < denom0)
        {
            Real xde0 = numer0/denom0;
            x0 = e0*xde0; x1 = e1*sqrt(1 - xde0*xde0);
            distance = sqrt(Sqr(x0 - y0) + Sqr(x1));
        }
        else
        {
            x0 = e0; x1 = 0; distance = abs(y0 - e0);
        }
    }
    return distance;
}
```

The function Sqr(t) simply returns $t^2$. The function RobustLength(v0,v1) computes the length of the input vector (v0, v1) by avoiding floating-point overflow that could occur normally when computing $v_0^2 + v_1^2$. If $|v_0| = \max\{|v_0|, |v_1|\}$, then $\sqrt{v_0^2 + v_1^2} = |v_0|\sqrt{1 + (v_1/v_0)^2}$. If $|v_1| = \max\{|v_0|, |v_1|\}$, then $\sqrt{v_0^2 + v_1^2} = |v_1|\sqrt{1 + (v_0/v_1)^2}$.
## 3 Distance from a Point to an Ellipsoid

A general ellipsoid in 3D is represented by a center point C, an orthonormal set of axis-direction vectors {U₀, U₁, U₂}, and associated extents eᵢ with e₀ ≥ e₁ ≥ e₂ > 0. The ellipsoid points are

$$ \mathbf{P} = \mathbf{C} + x_0\mathbf{U}_0 + x_1\mathbf{U}_1 + x_2\mathbf{U}_2 \quad (25) $$

where

$$ \left(\frac{x_0}{e_0}\right)^2 + \left(\frac{x_1}{e_1}\right)^2 + \left(\frac{x_2}{e_2}\right)^2 = 1 \quad (26) $$

If e₀ = e₁ = e₂, then the ellipsoid is a sphere with center C and radius e₀. If e₀ = e₁ > e₂, the ellipsoid is said to be an *oblate spheroid*. If e₀ > e₁ = e₂, the ellipsoid is said to be a *prolate spheroid*. The orthonormality of the axis directions and Equation (25) imply $x_i = U_i \cdot (\mathbf{P} - \mathbf{C})$. Substituting this into Equation (26) we obtain

$$ (\mathbf{P} - \mathbf{C})^{\mathrm{T}} M (\mathbf{P} - \mathbf{C}) = 1 \quad (27) $$

where M = RDRT, R is an orthogonal matrix whose columns are U₀, U₁, and U₂ and D is a diagonal matrix whose diagonal entries are 1/e₂₀, 1/e₂₁, and 1/e₂₂.

The problem is to compute the distance from a point Q to the ellipsoid. It is sufficient to solve this problem in the coordinate system of the ellipsoid; that is, represent Q = C + y₀U₀ + y₁U₁ + y₂U₂. The distance from Q to the closest point P on the ellipsoid as defined by Equation (25) is the same as the distance from Y = (y₀, y₁, y₂) to the closest point X = (x₀, x₁, x₂) on the standard ellipsoid of Equation (26).

As in the 2D ellipse problem, we may additionally use symmetry to simplify the construction. It is sufficient to consider the case when (y₀, y₁, y₂) is in the first octant: y₀ ≥ 0, y₁ ≥ 0, and y₂ ≥ 0.

### 3.1 The Closest Point's Normal is Directed Toward the Query Point

A parameterization of the standard ellipsoid is $X(\theta, \phi) = (e_0 \cos \theta \sin \phi, e_1 \sin \theta \sin \phi, e_2 \cos \phi)$ for $\theta \in [0, 2\pi)$ and $\phi \in [0, \pi]$. The squared distance from Y to any point on the ellipsoid is

$$ F(\theta, \phi) = |X(\theta, \phi) - Y|^2 \quad (28) $$

This is a nonnegative, doubly periodic, and differentiable function; it must have a global minimum occurring at angles for which the first-order partial derivatives are zero,

$$ \frac{\partial F}{\partial \theta} = 2(\mathbf{X}(\theta, \phi) - \mathbf{Y}) \cdot \frac{\partial \mathbf{X}}{\partial \theta} = 0, \quad \frac{\partial F}{\partial \phi} = 2(\mathbf{X}(\theta, \phi) - \mathbf{Y}) \cdot \frac{\partial \mathbf{X}}{\partial \phi} = 0 \quad (29) $$

For the derivatives to be zero, the vector $(\mathbf{X}(\theta, \phi) - \mathbf{Y})$ must be perpendicular to the tangent vectors $\partial\mathbf{X}/\partial\theta$ and $\partial\mathbf{X}/\partial\phi$. This implies that the vector from Y to the closest ellipsoid point X must be normal to the surface at X. Using the implicit form of the ellipsoid, namely, $G(x_0, x_1, x_2) = (x_0/e_0)^2 + (x_1/e_1)^2 + (x_2/e_2)^2 - 1$, half the gradient of $G(x_0, x_1, x_2)$ is a normal vector to the ellipsoid at $(x_0, x_1, x_2)$, so we have $(y_0, y_1, y_2) - (x_0, x_1, x_2) = t\nabla G(x_0, x_1, x_2)/2 = t(x_0/e_0^2, x_1/e_1^2, x_2/e_2^2)$ for some scalar t, or

$$ y_0 = x_0 \left( 1 + \frac{t}{e_0^2} \right), \quad y_1 = x_1 \left( 1 + \frac{t}{e_1^2} \right), \quad y_2 = x_2 \left( 1 + \frac{t}{e_2^2} \right) \quad (30) $$

If $(y_0, y_1, y_2)$ is outside the ellipsoid, it is necessary that $t > 0$. If $(y_0, y_1, y_2)$ is inside the ellipsoid, it is necessary that $t < 0$. If $(y_0, y_1, y_2)$ is already on the ellipsoid, then $t = 0$ and the distance is zero.
## 3.2 The Case of a Sphere

If $e_0 = e_1 = e_2$, then the ellipsoid is a sphere. The origin $(0,0,0)$ has infinitely many closest sphere points (all of them), but clearly the distance from the origin to the sphere is $e_0$. The closest sphere point to $(y_0, y_1, y_2) \neq (0,0,0)$ is $(x_0, x_1, x_2) = e_0(y_0, y_1, y_2)/|(y_0, y_1, y_2)|$. Equation (30) is consistent with this, because it implies $(x_0, x_1, x_2) = (e_0^2/(t+e_0^2))(y_0, y_1, y_2)$ for some $t$; that is, $(x_0, x_1, x_2)$ is parallel to $(y_0, y_1, y_2)$ and must have length $e_0$. It is easily shown that $t = -e_0^2 + e_0\sqrt{y_0^2 + y_1^2 + y_2^2}$.

## 3.3 The Case of an Oblate Spheroid

If $e_0 = e_1 > e_2$, then the ellipsoid is an oblate spheroid. The standard ellipse equation reduces to $(x_0^2 + x_1^2)/e_0^2 + x_2^2/e_2^2 = 1$, which is the equation of an ellipse in the $(r, x_2)$-plane where $r = \sqrt{x_0^2 + x_1^2}$. An implementation can make the reduction by mapping the query point $(y_0, y_1, y_2)$ to $(\sqrt{y_0^2 + y_1^2}, y_2)$ and using the point-ellipse algorithm discussed in Section 2.

## 3.4 The Case of a Prolate Spheroid

If $e_0 > e_1 = e_2$, then the ellipsoid is a prolate spheroid. The standard ellipse equation reduces to $x_0^2/e_0^2+(x_1^2+x_2^2/e_2^2) = 1$, which is the equation of an ellipse in the $(x_0, r)$-plane where $r = \sqrt{x_1^2 + x_2^2}$. An implementation can make the reduction by mapping the query point $(y_0, y_1, y_2)$ to $(y_0, \sqrt{y_1^2 + y_2^2}$ and using the point-ellipse algorithm discussed in Section 2.

For the remainder of Section 3, we assume that $e_0 > e_1 > e_2$ and that all analysis is restricted to the first octant.

## 3.5 The Query Point is the Origin

Equation (30) becomes $0 = x_i(1 + t/e_i^2)$ for $i = 0, 1, 2$. The construction of candidate points is similar to that shown in Section 2.3. There are eight subcases to consider. The subcase $x_0 = x_1 = x_2 = 0$ is discarded because the point is not on the ellipsoid. Any two subcases that require $t = -e_i^2$ and $t = -e_j^2$ with $i \neq j$ are contradictory because $e_i \neq e_j$; there are four such subcases. The remaining three subcases lead to candidates $(e_0, 0, 0)$, $(0, e_1, 0)$, and $(0, 0, e_2)$. Of these, $(0, 0, e_2)$ is closest to the origin. In summary: The closest ellipsoid point in the first octant to $(0, 0, 0)$ is $(0, 0, e_2)$ with distance $d = e_2$.

## 3.6 The Query Point has $y_2 > 0$

Equation (30) leads to $y_2 = x_2(1 + t/e_2^2)$. Because $y_2 > 0$ and the search for closest point is in the first octant, it must be that $x_2 > 0$ and $1 + t/e_2^2 > 0$. We have four cases to consider depending on whether $y_0$ or $y_1$ are zero or positive.

* Let $y_0 = 0$ and $y_1 = 0$; then $0 = x_0(1 + t/e_0^2)$ and $0 = x_1(1 + t/e_1^2)$. If $x_0 > 0$, then $t = -e_0^2$ and $y_2 = x_2(1 - e_0^2/e_1^2) < 0$. The last inequality is due to $e_0 > e_1$, but is a contradiction because $y_2 > 0$. Similarly, if $x_1 > 0$, then $t = -e_1^2$ and $y_2 = x_2(1 - e_1^2/e_2^2) < 0$, also a contradiction. It must be that $x_0 = 0$ and $x_1 = 0$, so the closest ellipsoid point is $(x_0, x_1, x_2) = (0, 0, e_2)$.
* Let $y_0 > 0$ and $y_1 = 0$; then $y_0 = x_0(1 + t/e_0^2)$, in which case $x_0 > 0$. We saw previously that $y_1 = 0$ forces $x_1 = 0$. This reduces the 3D point-ellipsoid problem to the 2D point-ellipse problem where the ellipse is $(x_0/e_0)^2 + (x_2/e_2)^2 = 1$ with query point $(y_0, y_2)$ whose components are both positive; see Section 2.6.

* Let $y_1 > 0$ and $y_0 = 0$; then $y_1 = x_1(1 + t/e_1^2)$, in which case $x_1 > 0$. We saw previously that $y_0 = 0$ forces $x_0 = 0$. This reduces the 3D point-ellipsoid problem to the 2D point-ellipse problem where the ellipse is $(x_1/e_1)^2 + (x_2/e_2)^2 = 1$ with query point $(y_1, y_2)$ whose components are both positive.; see Section 2.6.

* Let $y_0 > 0$ and $y_1 > 0$; then Equation (30) is solved for $x_i = e_i^2 y_i / (t + e_i^2)$ for $0 \le i \le 2$. The algorithm for computing the closest ellipsoid point is similar to that of Section 2.6. Substituting the $x_i$ into Equation (26), we obtain

$$F(t) = \left( \frac{e_0 y_0}{t + e_0^2} \right)^2 + \left( \frac{e_1 y_1}{t + e_1^2} \right)^2 + \left( \frac{e_2 y_2}{t + e_2^2} \right)^2 - 1 = 0 \quad (31)$$

The first equality defines the function $F(t)$ with domain $(-e_2^2, \infty)$. The candidates for the ellipsoid point $(x_0, x_1, x_2)$ closest to $(y_0, y_1, y_2)$ are generated by the roots $t$ to $F(t) = 0$.

The first-order and second-order derivatives of $F$ are

$$F'(t) = -\frac{2e_0^2 y_0^2}{(t + e_0^2)^3} - \frac{2e_1^2 y_1^2}{(t + e_1^2)^3} - \frac{2e_2^2 y_2^2}{(t + e_2^2)^3}, \quad F''(t) = \frac{6e_0^2 y_0^2}{(t + e_0^2)^4} + \frac{6e_1^2 y_1^2}{(t + e_1^2)^4} + \frac{6e_2^2 y_2^2}{(t + e_2^2)^4} \quad (32)$$

We know that $y_0 > 0$, $y_1 > 0$, $y_2 > 0$, and $t > -e_2^2$. These conditions imply $F'(t) < 0$ and $F''(t) > 0$ for $t > -e_2^2$. Observe that

$$\lim_{t \to -e_2^2^+} F(t) = +\infty, \quad \lim_{t \to \infty} F(t) = -1 \quad (33)$$

The first expression is a one-sided limit where $t$ approaches $e_2^2$ through values larger than $e_2^2$. We have shown that $F(t)$ is a strictly decreasing function for $t \in (-e_2^2, +\infty)$ that is initially positive, then becomes negative. Consequently it has a unique root on the specified domain. Figure 6 shows a typical graph of $F(t)$.

<figure><img src="image_8.png" /><figcaption>**Figure 6.** A typical graph of $F(t)$ for $t > -e_2^2$. The unique root $\bar{t}$ is shown.</figcaption></figure>

Observe that the domain $(-e_2^2, \infty)$ contains 0 and that $F(0) = (y_0/e_0)^2 + (y_1/e_1)^2 + (y_2/e_2)^2 - 1$. When $(y_0, y_1, y_2)$ is inside the ellipsoid we see that $F(0) < 0$ and $\bar{t} < 0$. When $(y_0, y_1, y_2)$ is outside
the ellipsoid we see that $F(0) > 0$ and $\bar{t} > 0$. When $(y_0, y_1, y_2)$ is on the ellipse we see that $F(0) = 0$ and $\bar{t} = 0$.

For all $t$, a typical graph of $F(t)$ is shown in Figure 7.

<figure><img src="image_9.png" /><figcaption>Figure 7. A typical graph of $F(t)$ when $y_0 > 0$, $y_1 > 0$, and $y_2 > 0$.</figcaption></figure>

The illustration shows $F(t)$ with 4 roots. However, it is possible for $F(t)$ to have up to 6 roots.

## 3.7 The Query Point has $y_2 = 0$

The construction is analogous to that of Section 2.5. Equation (30) states that $y_0 = x_0(1 + t/e_0^2)$, $y_1 = x_1(1 + t/e_1^2)$, and $0 = x_2(1 + t/e_2^2)$. The last equation has two cases.

* Let $x_2 = 0$. For $(x_0, x_1, 0)$ to be on the ellipsoid, we need $(x_0/e_0)^2 + (x_1/e_1)^2 = 1$. The 3D point-ellipsoid problem is therefore reduced to the 2D point-ellipse problem with query point $(y_0, y_1)$.

* Let $x_2 > 0$; then $t = -e_2^2$, $y_0 = x_0(1 - e_2^2/e_0^2)$, and $y_1 = x_1(1 - e_2^2/e_1^2)$. It follows that $x_0 = e_0^2 y_0 / (e_0^2 - e_2^2)$ and $x_1 = e_1^2 y_1 / (e_1^2 - e_2^2)$. Because $(x_0, x_1, x_2)$ is on the ellipsoid, we know that

$$1 \ge 1 - \left(\frac{x_2}{e_2}\right)^2 = \left(\frac{x_0}{e_0}\right)^2 + \left(\frac{x_1}{e_1}\right)^2 = \left(\frac{e_0 y_0}{e_0^2 - e_2^2}\right)^2 + \left(\frac{e_1 y_1}{e_1^2 - e_2^2}\right)^2 \quad (34)$$

The implication is that $(y_0, y_1)$ is contained by a subellipse of the $x_0x_1$-domain for the hemiellipsoid $x_2 \ge 0$. For $(y_0, y_1)$ outside this subellipse, the closest point is necessarily on the ellipse $(x_0/e_0)^2 + (x_1/e_1)^2 = 1$. When $(y_0, y_1)$ is inside the subellipse, for $(x_0, x_1, x_2)$ to be on the ellipsoid, we may solve for $x_2 = e_2 \sqrt{1 - (x_0/e_0)^2 - (x_1/e_1)^2}$.

## 3.8 A Summary of the Mathematical Algorithm

Listing 3 summarizes the algorithm for computing the ellipsoid point $(x_0, x_1, x_2)$ closest to the query point $(y_0, y_1, y_2)$ in the first octant and for computing the distance to the closest point. The preconditions are $e_0 \ge e_1 \ge e_2 > 0$, $y_0 \ge 0$, $y_1 \ge 0$, and $y_2 \ge 0$.
<caption>**Listing 3.** Pseudocode for computing the closest ellipsoid point and distance to a query point. It is required that $e_0 \ge e_1 \ge e_2 > 0$, $y_0 \ge 0$, $y_1 \ge 0$, and $y_2 \ge 0$.</caption>

```c
Real DistancePointEllipsoid(Real e0, Real e1, Real e2, Real y0, Real y1, Real y2,
                             Real& x0, Real& x1, Real& x2)
{
    Real distance;
    if (y2 > 0)
    {
        if (y1 > 0)
        {
            if (y0 > 0)
            {
                Compute the unique root tbar of F(t) on (-e2*e2, +infinity);
                x0 = e0*e0*y0/(tbar + e0*e0);
                x1 = e1*e1*y1/(tbar + e1*e1);
                x2 = e2*e2*y2/(tbar + e2*e2);
                distance = sqrt((x0 - y0)*(x0 - y0) + (x1 - y1)*(x1 - y1) + (x2 - y2)*(x2 - y2));
            }
            else // y0 == 0
            {
                x0 = 0;
                distance = DistancePointEllipse(e1, e2, y1, y2, x1, x2);
            }
        }
        else // y1 == 0
        {
            x1 = 0;
            if (y0 > 0)
            {
                distance = DistancePointEllipse(e0, e2, y0, y2, x0, x2);
            }
            else // y0 == 0
            {
                x0 = 0;
                x2 = e2;
                distance = abs(y2 - e2);
            }
        }
    }
    else // y2 == 0
    {
        Real denom0 = e0*e0 - e2*e2, denom1 = e1*e1 - e2*e2;
        Real numer0 = e0*y0, numer1 = e1*y1;
        bool computed = false;
        if (numer0 < denom0 && numer1 < denom1)
        {
            Real xde0 = numer0/denom0, xdel = numer1/denom1;
            Real xde0sqr = xde0*xde0, xdelsqr = xdel*xdel;
            Real discr = 1 - xde0sqr - xdelsqr;
            if (discr > 0)
            {
                x0 = e0*xde0;
                x1 = e1*xdel;
                x2 = e2*sqrt(discr);
                distance = sqrt((x0 - y0)*(x0 - y0) + (x1 - y1)*(x1 - y1) + x2*x2);
                computed = true;
            }
        }
        if (!computed)
        {
            x2 = 0;
            distance = DistancePointEllipse(e0, e1, y0, y1, x0, x1);
        }
    }
    return distance;
}
```
## 3.9 Robust Root Finding and Conversion to a Polynomial

The ideas of Section 2.8 apply equally well in the point-ellipsoid problem. The bisection method is robust. The finite bounding interval for $\bar{t}$ is $[t_0, t_1]$ with $t_0 = -e_2^2 + e_2y_2$ and $t_1 = -e_2^2 + \sqrt{(e_0y_0)^2 + (e_1y_1)^2 + (e_2y_2)^2}$.

Newton's Method may be used with initial guess $t_0$, but the same numerical issues must be addressed when $y_2$ is nearly zero. In the limit as $y_2$ approaches to zero, the vertical asymptote disappears and there is a discontinuity in the root $\bar{t}$ when viewed as a function of $y_2$.

The conversion to a polynomial of degree 6 is suggested in [1], but has the same problems mentioned for the point-ellipse problem. The roots of the polynomial may be computed, leading to candidates for the closest ellipsoid point,

$$ 
\begin{aligned}
P(t) &= (t + e_0^2)^2(t + e_1^2)^2(t + e_2^2)^2 F(t) \\
&= e_0^2 y_0^2 (t + e_1^2)^2 (t + e_2^2)^2 + e_1^2 y_1^2 (t + e_0^2)^2 (t + e_2^2)^2 + e_2^2 y_2^2 (t + e_0^2)^2 (t + e_1^2)^2
\end{aligned} 
\quad (35) $$

When $y_2 > 0$, the roots of $P(t)$ are the same as those for $F(t)$, and the largest root of $P(t)$ is the same as the largest root of $F(t)$. However, when $y_2 = 0$, $P(t)$ has a double root at $t = -e_2^2$ but $F(t)$ does not have a root at $t = -e_2^2$. Just as in the point-ellipse distance algorithm, this causes numerical problems when $y_2$ is nearly zero. In fact, numerical problems occur when any of the $y_i$ are nearly zero.

The Graphics Gems IV article suggests using Newton's Method starting with an initial guess greater than the largest root of $P(t)$, namely, $t_0 = e_0 \sqrt{y_0^2 + y_1^2 + y_2^2}$. As in Section 2.8.3, the convergence is not guaranteed, especially when $y_2$ is nearly zero. The discussion about double roots in that section applies here as well, and if you choose to use Newton's Method for $P(t)$, a better initial guess is $t_1 = -e_1 + \sqrt{e_0^2 y_0^2 + e_1^2 y_1^2 + e_2^2 y_2^2}$. It may be shown that $t_1 < t_1 + e_2 < t_0$.

## 3.10 A Robust Implementation

An implementation that uses the robust bisection method is described next. A change of variables helps keep the magnitude of the numbers to a moderate size, $s = t/e_2^2$. Te query point components are also scaled by the ellipsoid extents to keep the inputs on the order of 1 when the query is executed for points near the ellipsoid, $z_i = y_i/e_i$. Define the ratios of squared extents $r_0 = e_0^2/e_2^2$ and $r_1 = e_1^2/e_2^2$. The function $F(t)$ is transformed to

$$ G(s) = \left( \frac{r_0 z_0}{s + r_0} \right)^2 + \left( \frac{r_1 z_1}{s + r_1} \right)^2 + \left( \frac{z_2}{s + 1} \right)^2 - 1 \quad (36) $$

for $s \in (-1, +\infty)$. The root of $G(s) = 0$ is $\bar{s} = \bar{t}/e_2^2$ and it is in the interval $[s_0, s_1]$, where $s_0 = -1 + z_2$ and $s_1 = -1 + \sqrt{r_0^2 z_0^2 + r_1^2 z_1^2 + z_2^2}$.

Listing 4 is the robust implementation of the algorithm. The maximum number of iterations is provided by Equation (18).

**Listing 4.** Pseudocode for robustly computing the closest ellipsoid point and distance to a query point. It is required that $e_0 \ge e_1 \ge e_2 > 0$, $y_0 \ge 0$, $y_1 \ge 0$, and $y_2 \ge 0$.

```
Real GetRoot(Real r0, Real r1, Real z0, Real z1, Real z2, Real g)
{
    Real n0 = r0*z0, n1 = r1*z1;
```
```c
Real s0 = z2 - 1, s1 = (g < 0 ? 0 : RobustLength(n0, n1, z2) - 1);
Real s = 0;
for (int i = 0; i < maxIterations; ++i)
{
    s = (s0 + s1) / 2;
    if (s == s0 || s == s1) { break; }
    Real ratio0 = n0/(s + r0), ratio1 = n1/(s + r1), ratio2 = z2/(s + 1);
    g = Sqr(ratio0) + Sqr(ratio1) + Sqr(ratio2) - 1;
    if (g > 0) { s0 = s; } else if (g < 0) { s1 = s; } break; }
return s;
}

Real DistancePointEllipsoid(Real e0, Real e1, Real e2, Real y0, Real y1, Real y2,
Real& x0, Real& x1, Real& x2)
{
    Real distance;
    if (y2 > 0)
    {
        if (y1 > 0)
        {
            if (y0 > 0)
            {
                Real z0 = y0/e0, z1 = y1/e1, z2 = y2/e2;
                Real g = Sqr(z0) + Sqr(z1) + Sqr(z2) - 1;
                if (g != 0)
                {
                    Real r0 = Sqr(e0/e2), r1 = Sqr(e1/e2);
                    Real sbar = GetRoot(r0, r1, z0, z1, z2, g);
                    x0 = r0*y0/(sbar + r0); x1 = r1*y1/(sbar + r1); x2 = y2/(sbar + 1);
                    distance = sqrt(Sqr(x0 - y0) + Sqr(x1 - y1) + Sqr(x2 - y2));
                }
                else
                {
                    x0 = y0; x1 = y1; x2 = y2; distance = 0;
                }
            }
            else // y0 == 0
            {
                x0 = 0; distance = DistancePointEllipse(e1, e2, y1, y2, x1, x2);
            }
        }
        else // y1 == 0
        {
            if (y0 > 0)
            {
                x1 = 0; distance = DistancePointEllipse(e0, e2, y0, y2, x0, x2);
            }
            else // y0 == 0
            {
                x0 = 0; x1 = 0; x2 = e2; distance = abs(y2 - e2);
            }
        }
    }
    else // y2 == 0
    {
        Real denom0 = e0*0 - e2*e2, denom1 = e1*e1 - e2*e2, numer0 = e0*y0, numer1 = e1*y1;
        bool computed = false;
        if (numer0 < denom0 && numer1 < denom1)
        {
            Real xde0 = numer0/denom0, xde1 = numer1/denom1;
            Real xde0sqr = xde0*xde0, xde1sqr = xde1*xde1;
            Real discr = 1 - xde0sqr - xde1sqr;
            if (discr > 0)
            {
                x0 = e0*xde0; x1 = e1*xde1; x2 = e2*sqrt(discr);
                distance = sqrt((x0 - y0)*(x0 - y0) + (x1 - y1)*(x1 - y1) + x2*x2);
                computed = true;
            }
        }
        if (!computed)
        {
```
```c
        x2 = 0; distance = DistancePointEllipse(e0, e1, y0, y1, x0, x1);
    }
}
return distance;
}
```

The function Sqr(t) simply returns t². The function RobustLength(v0,v1,v2) computes the length of the input vector V = (v₀, v₁, v₂) by avoiding floating-point overflow that could occur normally when computing v₀² + v₁² + v₂². If |vₖ| = max{|v₀|, |v₁|, |v₂|}, then |V| = |vₖ||V/vₖ|.

The code calls the point-ellipse distance function when the point-ellipsoid query reduces to a 2D problem. As we will see in the last section of this document, the recursion on dimension is not necessary.

# 4 Distance from a Point to a Hyperellipsoid

A general hyperellipsoid in n dimensions is represented by a center point C, an orthonormal set of axis-direction vectors {Uᵢ}n-1ᵢ=0, and associated extents eᵢ with eⱼ ≥ eⱼ₊₁ > 0 for all j. The hyperellipsoid points are

$$ \mathbf{P} = \mathbf{C} + \sum_{i=0}^{n-1} x_i \mathbf{U}_i \qquad (37) $$

where

$$ \sum_{i=0}^{n-1} \left( \frac{x_i}{e_i} \right)^2 = 1 \qquad (38) $$

As with an ellipse or an ellipsoid, special cases occur when some of the extents are equal. The implementation handles these correctly, so these are not required to have special code to deal with them. The orthonormality of the axis directions and Equation (37) imply xᵢ = Uᵢ · (P - C). Substituting this into Equation (38), we obtain

$$ (\mathbf{P} - \mathbf{C})^{\mathrm{T}} M (\mathbf{P} - \mathbf{C}) = 1 \qquad (39) $$

where M = RDRT, R is an orthogonal matrix whose columns are the Uᵢ vectors, and D is a diagonal matrix whose diagonal entries are the 1/eᵢ² values.

The problem is to compute the distance from a point Q to the hyperellipsoid. It is sufficient to solve this problem in the coordinate system of the hyperellipsoid; that is, represent Q = C + Σn-1ᵢ=0 yᵢUᵢ. The distance from Q to the closest point P on the hyperellipsoid is the same as the distance from Y = (y₀, ..., yₙ₋₁) to the closest point X = (x₀, ..., xₙ₋₁) on the standard hyperellipsoid of Equation (38).

We may additionally use symmetry to simplify the construction. It is sufficient to consider the case when yᵢ ≥ 0 for all i.

## 4.1 The Closest Point's Normal is Directed Toward the Query Point

A parameterization of the standard hyperellipsoid involves n-1 parameters Θ = (θ₀, ..., θₙ₋₂), say, X(Θ). The actual form of the components can use a generalization of spherical coordinates. The squared distance from Y to any point on the hyperellipsoid is

$$ F(\Theta) = |\mathbf{X}(\Theta) - \mathbf{Y}|^2 \qquad (40) $$
This is a nonnegative differentiable function that is periodic in each of its independent variables; it must have a global minimum occurring at parameters for which the first-order partial derivatives are zero,

$$ \frac{\partial F}{\partial \theta_i} = 2(\mathbf{X}(\mathbf{\Theta}) - \mathbf{Y}) \cdot \frac{\partial \mathbf{X}}{\partial \theta_i} = 0 \qquad (41) $$

for all i. For the derivatives to be zero, the vector $(\mathbf{X}(\mathbf{\Theta}) - \mathbf{Y})$ must be perpendicular to the tangent vectors $\partial\mathbf{X}/\partial\theta_i$. This implies the vector from $\mathbf{Y}$ to the closest hyperellipsoid point $\mathbf{X}$ must be normal to the surface at $\mathbf{X}$. Using the implicit form of the ellipsoid, namely, $G(\mathbf{X}) = \mathbf{X}^T D \mathbf{X} - 1$, half the gradient of $G(\mathbf{X})$ is a normal vector to the hyperellipsoid at $\mathbf{X}$, so we have $\mathbf{Y} - \mathbf{X} = t \nabla G(\mathbf{X}) / 2 = t D \mathbf{X}$ for some scalar $t$, or

$$ y_i = x_i \left( 1 + \frac{t}{e_i^2} \right) \qquad (42) $$

for all i. If $\mathbf{Y}$ is outside the hyperellipsoid, it is necessary that $t > 0$. If $\mathbf{Y}$ is inside the hyperellipsoid, it is necessary that $t < 0$. If $\mathbf{Y}$ is already on the hyperellipsoid, then $t = 0$ and the distance is zero.

## 4.2 The Key Observations

In the point-ellipse and point-ellipsoid distance algorithms, an important observation is that the main branch- ing depends on whether the last component of $\mathbf{Y}$ is positive or zero.

* When $y_{n-1} > 0$, notice that whenever $y_i = 0$ for any $i < n - 1$, the corresponding closest point has component $x_i = 0$.

* When $y_{n-1} = 0$, notice that two cases occur. The projection of $\mathbf{Y}$ onto the coordinate hyperplane $y_{n-1} = 0$ is either inside or outside a special hyperellipsoid living in $n - 1$ dimensions (the original hyperellipsoid lives in $n$ dimensions). If inside the special hyperellipsoid, we can compute the corre- sponding $x_i$ directly from Equation (42), and then compute $x_{n-1} = \sqrt{1 - \sum_{i=0}^{n-2} (x_i / e_i)^2}$. If outside the special hyperellipsoid, we have the recursion in dimension—the query reduces to point-hyperellipsoid within the hyperplane $y_{n-1} = 0$.

This suggests that we should compute a vector $\mathbf{Y}'$ from $\mathbf{Y}$ by storing all the positive components, construct a vector $\mathbf{e}'$ that stores the extents corresponding to those positive components, and then reduce the problem either to a bisection in a smaller-dimensional space or to a direct computation of the components in a smaller-dimensional space.

Let there be $p$ positive components of $\mathbf{Y}$ located at the indices $i_0$ through $i_{p-1}$. If $y_{n-1} > 0$, then $i_{p-1} = n-1$. Construct $\mathbf{Y}' = (y_{i_0}, \dots, y_{i_{p-1}})$. Let $\mathbf{e}' = (e_{i_0}, \dots, e_{i_{p-1}})$ and $\mathbf{X}' = (x_{i_0}, \dots, x_{i_{p-1}})$. The function to bisect is

$$ F(t) = \sum_{j=0}^{p-1} \left( \frac{e_{ij} y_{ij}}{t + e_{ij}^2} \right)^2 - 1 \qquad (43) $$

for $t \in (-e_{i_{p-1}}^2, +\infty)$. The pseudocode is shown in Listing 5.

**Listing 5.** Pseudocode for computing the distance from a point to a hyperellipsoid. The extent vector must have nonincreasing components and the query point must have nonnegative components.