I don’t think those equations are right… I just tried them now in the 2D case, and they didn’t give the right answer at all. (I used an ellipse with a = 2 and b = 1, and (x0, y0) = (2, 1.5), and plotted the parametric curve according to the equation you posted. It goes through a point on the ellipse, but definitely not point on the ellipse nearest to (2, 1.5).)

But, assuming that some typo was made, I believe I can give a partial explanation:

Imagine taking the point (x0, y0, z0) and drawing a small sphere around it, then gradually expanding that sphere till it just touches the ellipsoid. The point where it touches the ellipsoid is the place you’re looking for, and since the sphere and ellipsoid are just touching at that point, they’re tangent there. This means the ellipsoid’s normal at that point is in the opposite direction to the sphere’s normal - and since the sphere’s normal points directly away from its center, it follows the ellipsoid’s normal points directly towards the sphere center, aka (x0, y0, z0).

So we know the location we want is the point on the ellipsoid where its normal points toward (x0, y0, z0). Now, if you define an ellipsoid as a scalar field f(x,y,z) where f = 0 for points on the surface, f < 0 inside and f > 0 outside, then the normal of an ellipsoid is given by the gradient of f.

Now, here is where I’m not sure what your author is doing anymore - he seems to be trying to construct a 1D curve that (somehow) will go through the right point on the ellipsoid. At this point it makes more sense to me to parameterize the ellipsoid as a function P(u,v) mapping u,v to some point on the ellipsoid surface. Then, what we want to find are the u,v for which

(x0, y0, z0) - P(u,v) = k * gradientF(P(u,v))

for some k, i.e. the vector from the point on the surface to (x0, y0, z0) is in the same direction as the gradient at that point. This can be re-expressed as a function that normalizes both vectors and takes their dot product:

G(u,v) = normalize((x0, y0, z0) - P(u,v)) dot normalize(gradientF(P(u,v)))

Now, this function will have a maximum exactly at the u,v coordinates where you’ll find the point you’re looking for. You have to apply a two-dimensional maximization now, but that’s not so bad as the function is very well behaved, there are no local maxima. You could even probably solve it analytically, though it would be messy (and possibly the numerical solution would perform better).

To determine the closest point on an ellipsoid’s surface to any given other point (x0, y0, z0), I have found a method of parameterizing the ellipsoid equation by a parameter t, so that:

f(t) = a² * x0² / (a² + 2t)² + b² * y0² / (b² + 2t)² + c² * z0² / (c² + 2t)² - 1

Which you will find is simply the standard equation, where x, y and z have been substituted by:

x(t) = x0 * a² / (a² + 2t)

y(t) = y0 * b² / (b² + 2t)

z(t) = z0 * c² / (c² + 2t)

Since the 1 has been pulled to the other side, the root of f(t) will be the value for t at which the components x(t), y(t) and z(t) of a new point will be such that the new point is on the ellipsoid surface. (Since if the equation is 1, this is the case.) This point is the one we are searching for.

To find this root, a simple Newton-Raphson approximation is applied to the function. I’ve tested this out, and it works just fine. However, I’m not sure why x(t), y(t) and z(t) are defined the way they are.

It is said in the source of this approach that it is done this way because “the derivative of the ellipsoid equation at the nearest point on the ellipsoid to the skin point represents a vector between the skin point and its nearest point”. This is one thing I don’t grasp. The other is that the three functions are called components of a “parametric line function”, though when I input different values for t and visualized this, the resulting path in 3d space was clearly a curve, not a straight line.

The author then just describes the approximation as moving step by step on the line/curve created by the three functions towards the t value that gives a point almost on the ellipsoid surface. Another important property of the functions he mentions is that for t = 0, the resulting point on the line/curve is exactly the point we want to find the nearest point for (as one can easily see).

So, all in all, I’d like someone to tell me how the author got to these three equations.