0
$\begingroup$

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. enter image description here

$\endgroup$
  • $\begingroup$ I think your code was written for the case of an ellipsoid whose axes are the same as cartesian axes. $\endgroup$ – Intelligenti pauca Jun 28 '17 at 20:28
  • $\begingroup$ That may be the case. Every reference I look at gives the ellipsoid equation assuming the ellipsoid is centered on the origin and the axis are the coordinate system axis. How do I incorporate the axis vectors into my equation then? $\endgroup$ – Josh Jun 28 '17 at 20:35
  • $\begingroup$ Substitute in your code pos_norm.X with the dot product of pos.norm with axis C vector, and so on for pos_norm.Y and pos_norm.Z. After that, 1/Tparam should give the relative radial position. $\endgroup$ – Intelligenti pauca Jun 28 '17 at 20:43
1
$\begingroup$

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. $$

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ I edited my post to give more information and a concrete example. I am not sure that I understand what I am really computing with Tparam, and thus how 1 / Tparam gives the correct answer? $\endgroup$ – Josh Jun 28 '17 at 20:20
  • $\begingroup$ I edited my answer: hope it is clear now. $\endgroup$ – Intelligenti pauca Jun 29 '17 at 7:20
  • $\begingroup$ I was able to get my code working from what you told me yesterday, but it is great to see it more formalized, and I think that I understand now. Thank you very much with all of your help! $\endgroup$ – Josh Jun 29 '17 at 16:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.