Find the relative radial position of a point within an ellipsoid

Question

I have parameters for a general arbitrary ellipsoid in a 3D cartesian coordinate system:

• ellipsoid centered on (x0, y0, z0)
• axis A with vector (xa, ya, za) and length la
• axis B with vector (xb, yb, zb) and length lb
• axis C with vector (xc, yc, zc) and length lc

I would like to find the relative radial position of a 3D point with coordinates position within the ellipsoid (0 == center; 1 == periphery).

Currently, I am doing the following:

var pos_norm = position - subsurface.CenterOfMass;
var Tparam = Math.Sqrt(1 /
                          ( (Math.Pow(pos_norm.X, 2) / Math.Pow(lc, 2))
                          + (Math.Pow(pos_norm.Y, 2) / Math.Pow(lb, 2))
                          + (Math.Pow(pos_norm.Z, 2) / Math.Pow(la, 2))
                          ));
var surfacePoint = pos_norm * Tparam;
var sdist = pos_norm.Distance(surfacePoint);
var eradius = surfacePoint.Distance(Vector3.Zero);
return 1 - (sdist / eradius);

This seems to generally work, but I am observing some cases where this seems to produce incorrect radial positions, almost inverted. I wonder if I am getting into trouble by ignoring the ellipsoid axis vectors. How do I correctly perform this calculation in the general case? Thanks.

EDIT: The code above was inherited, and is purported to compute a point on the ellipsoid surface that lies on the line intersecting the ellipsoid center and the query point, and finally takes the ratio of the segments on either side of the query point as the radial position.

Here is a concrete example with numbers:

• ellipsoid center at (14.329, 10.476, 3.241)
• axis A with vector (-0.024, -0.028, -0.999) and length 2.896
• axis B with vector (0.735, 0.677, -0.036) and length 4.096
• axis C with vector (0.678, -0.735, 0.004) and length 10.224
• query point at (14.000, 10.000, 3.241)

using the code aboveTparam = 6.831 which puts the point on the ellipsoid at (-2.247, -3.252, -1.646). So then sdist = 21.475 and eradius = 4.282 and thus the return value is -4.015. The expected return value should be close to zero.

However the surface point computed is clearly way off, and nowhere near the ellipsoid surface. In the image below, the transparent blue is an isosurface from which the ellipsoid parameters have been determined, with the center as the dark blue square. The sphere within the surface is the query point in the concrete example above, and the sphere outside the surface is what I computed as the ellipsoid surface point surfacePoint.

Answer 1

You have an ellipsoid, centered at $O$, with mutually orthogonal axes of length $l_a$, $l_b$, $l_c$, whose directions are given by three unit vectors $\vec{n_a}$, $\vec{n_b}$, $\vec{n_c}$. Then a point $Q$ lies on the ellipsoid if the following equation holds: $$ {(\vec{r}\cdot\vec{n_a})^2\over l_a^2}+ {(\vec{r}\cdot\vec{n_b})^2\over l_b^2}+ {(\vec{r}\cdot\vec{n_c})^2\over l_c^2}=1, $$ where I set $\vec r=(Q-O)$ and $\cdot$ denotes the dot product between two vectors: $\vec{a}\cdot\vec b=a_xb_x+a_yb_y+a_zb_z$.

If we now have a point $P$ inside the ellipsoid, all points $Q$ of ray $OP$ can be pameterized as $Q-O=t(P-O)$, with $t\ge0$. For $t=1$, in particular, you get point $P$ itself, and for some $t>1$ the corresponding $Q$ lies on the ellipsoid.

To find this value of $t$, set $\vec p=P-O$, substitute $\vec r=t\vec p$ into the above equation and solve for $t$. The result is:

$$ {1\over t}=\sqrt{ {(\vec{p}\cdot\vec{n_a})^2\over l_a^2}+ {(\vec{p}\cdot\vec{n_b})^2\over l_b^2}+ {(\vec{p}\cdot\vec{n_c})^2\over l_c^2}. } $$ As you are interested in the ratio $PO/QO$, from the above definitions it follows that this is equal to $1/t$ given above.

With your concrete example, for instance: $$ p=(14.000, 10.000, 3.241) - (14.329, 10.476, 3.241)=(-0.329, -0.476, 0),\\ \vec{p}\cdot\vec{n_a}=0.021224,\quad \vec{p}\cdot\vec{n_b}=-0.564067,\quad \vec{p}\cdot\vec{n_c}=0.126798,\\ 1/t=0.138463. $$