Thread Subject: Shortest distance from point to ellipsoid surface

I have a ellipsoid "defined" at a point E. I have its shape information stored:

elpsd_props{i,1}(:,:) = [centroid xyz coords;
                                  unit semi-axis vector a;
                                  unit semi-axis vector b;
                                  unit semi-axis vector c;
                                  semi-axis lengths a,b,c]

I have a matrix called astrd_int which holds a massive grid 3-D points. All of these points lie outside of the ellipsoid. I call these interior points (we're inside another shape).

I wish to find the minimum distance from each interior point to the surface of the ellipsoid without the use of generating points on the surface of the ellipsoid.

I already find all of the minimum distances by generating the points on the ellipsoid, calculating distances, and taking the minimum, for each point. I don't wish to do this because the points generated by the ellipsoid function are finite, so if my grid of 3-D points becomes very dense (many many points), it's likely that their minimum distance will "link" to an ellipsoid surface point which is farther away than the point on the real surface. I could control the number of points generated but even so I would rather have a "true" solution, not one that relied on generated surface points.

Is there a better way to do this without using ellipsoid surface points?

Define "massive" - are we talking about tens of millions of voxels
here? Personally I'd use bwdist (in the Image Processing Toolbox) and
just get a quantized value. That should probably be good enough - if
it's not accurate enough explain why not and what you're doing -
what's the big picture, what are you really trying to solve?

If you really need a perfect analytical solution using trigonometry
and geometry then you'll probably have to hope Roger gives it to you,
because it's way more work than I want to put into it.

Thank you for replying,

I apologize for not specifying; my grid of 3-D points usually contains between 25,000 and 50,000 "interior points". Usually I'll generate (200+1)^2 points using the mesh command, and then they're "filtered" by checking if they lie inside of an asteroid.

I'm improving a self-made finite element analysis program for simulating gravitational fields of asteroids. The primary goal is to use the fewest sphere and/or ellipsoid elements as possible; hence, the goal becomes defining the largest possible elements inside of the asteroid.

I've aced the sphere element technique without needing to generate sphere surface points, but I'm trying to optimize my ellipsoid element technique.

I haven't used the Image Processing Toolbox before, so I'll have to dig around and familiar myself with bwdist.

I'm trying to find the largest possible ellipsoid inside of an asteroid (of arbitrary shape: convex, concave, etc.). I need to assign to each interior point its minimum distance from itself to a previously-defined ellipsoid surface. It would take a while to explain why I need these distances lol.

The vectors of and distances from the ellipsoid center to all interior points are calculated for other uses. If those vectors are the same vectors along which the minimum distance would be found, but if they were, I'm sure that they could be parametrized and solved in conjunction with the ellipsoid surface equation.

And yes Roger has helped me a TON before, and I truly appreciate all of your help!

Robert:
It sure sounds like bwdist will do what you want, as long as you can
live with the quantization error. It finds the Euclidean distance
between every voxel in the interior and the voxel closest to it that
lies on the enclosing surface. But it sounds like you already have
these distances because you say they were "calculated for other uses"
so maybe bwdist won't add anything useful beyond what you have.

 It sounds like you're quantized already anyway because you're talking
about finding the distance between interior points and the surface as
if there are a finite number of those interior points, when if it
weren't quantized, there would be an infinite number of points
inside. But again, this assumes you'd use bwdist which may not be
necessary.

And your volume isn't yet massive so you could probably even increase
the resolution more (tinier voxels and more of them) before you run
out of memory.

John D'Errico used to have some nice routines about finding the
largest or smallest shapes (circles, rectangles) encapsulating or
within a cluster of points, but, very unfortunately, he removed them
from the File Exchange.

Sorry I can't help more but my math skills are pale in comparison to
Roger's. Maybe Roger will replicate all John's old submissions ;-)

If you don't want to fit the *largest* ellipsoid *inside* your
irregular surface but are willing to live with a *fit* of an ellipsoid
to the surface instead (where the fitted ellipsoid is inside the
suface about half the time and outside the surface about half the
time), then regionprops can do that.

One final suggestion, it sounds like you might be interested in
quadtree decomposition, which was pretty much made for this kind of
situation - look it up in the Image Processing Toolbox help. I'd be
willing to gamble that is what you really might need.
-ImageAnalyst

"Robert Phillips" <phillir1@my.erau.edu> wrote in message <ivd95k$6h2$1@newscl01ah.mathworks.com>...

> Is there a better way to do this without using ellipsoid surface points?

See my FEX:

http://www.mathworks.com/matlabcentral/fileexchange/27711-euclidian-projection-on-ellipsoid-and-conic

It might be slow, but the method gives exact solution(s).

Bruno

"Robert Phillips" <phillir1@my.erau.edu> wrote in message <ivd95k$6h2$1@newscl01ah.mathworks.com>...
> I have a ellipsoid "defined" at a point E. I have its shape information stored:
>
> elpsd_props{i,1}(:,:) = [centroid xyz coords;
> unit semi-axis vector a;
> unit semi-axis vector b;
> unit semi-axis vector c;
> semi-axis lengths a,b,c]
>
> I have a matrix called astrd_int which holds a massive grid 3-D points. All of these points lie outside of the ellipsoid. I call these interior points (we're inside another shape).
>
> I wish to find the minimum distance from each interior point to the surface of the ellipsoid without the use of generating points on the surface of the ellipsoid.
- - - - - - - - - -
  Finding the minimum distance from an arbitrary point to the surface of a specified ellipsoid is not a problem with a simple solution. I have done part of the work and left the rest to you, Robert.

  The part I have not done is to translate and rotate the ellipsoid, along with the set of "grid" points, so that the ellipsoid's new coordinate satisfy a simpler equation of the form:

 x^2/a + y^2/b + z^2/c - 1 = 0

where a, b, and c are the squares of the three semi-axis lengths. (I have used their squares here to simplify the form of the analysis that follows.)

  The part I have worked out for you is, having done the above transformation for the ellipsoid, then for each (transformed) point (x0,y0,z0), we must find the minimum distance to a point (x,y,z) which lies on the ellipsoid. That is, we are to minimize (x-x0)^2+(y-y0)^2+(z-z0)^2 subject to the constraint of the above equation. This is a problem that calls for finding a Lagrange multiplier.

  In analysis it is known that extremums must have the property that the set of partial derivatives of the function to be minimized must be proportional to the corresponding partial derivatives of the constraint function. Hence we must have 2*(x-x0),2*(y-y0),2*(z-z0) proportional to 2*x/a,2*y/b,2*z/c. If we call the proportionality factor (the Lagrange multiplier) t, this gives the equations

 t = a*(x-x0)/x = b*(y-y0)/y = c*(z-z0)/z

which leads to

 x = a/(a-t)*x0, y = b/(b-t)*y0, z = c/(c-t)*z0

Substituting these into the ellipsoid equation gives

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

and from this we obtain

 (a-t)^2*(b-t)^2*(c-t)^2 =
 a*(b-t)^2*(c-t)^2*x0^2+b*(c-t)^2*(a-t)^2*y0^2+c*(a-t)^2*(b-t)^2*z0^2

  This is a sixth order polynomial equation in the multiplier t, and we can solve it using matlab's roots function. I have worked out its seven coefficients, courtesy of matlab's symbolic feature, as follows:

 c1 = 1;
 c2 = -2*(a+b+c);
 c3 = a^2+b^2+c^2+4*(b*c+c*a+a*b)-(a*x0^2+b*y0^2+c*z0^2);
 c4 = 2*(-a^2*(b+c)-b^2*(c+a)-c^2*(a+b) ...
      -4*a*b*c+a*(b+c)*x0^2+b*(c+a)*y0^2 +c*(a+b)*z0^2);
 c5 = b^2*c^2+c^2*a^2+a^2*b^2+4*a*b*c*(a+b+c) ...
      -a*(b^2+c^2)*x0^2-b*(c^2+a^2)*y0^2-c*(a^2+b^2)*z0^2 ...
      -4*a*b*c*(x0^2+y0^2+z0^2);
 c6 = 2*a*b*c*((b+c)*x0^2+(c+a)*y0^2+(a+b)*z0^2-b*c-c*a-a*b);
 c7 = a*b*c*(a*b*c-b*c*x0^2-c*a*y0^2-a*b*z0^2);

(I am sure with a little work you could further simplify these expressions. I got tired messing around with them.)

  The rest of the code could look something like this:

 t = roots([c1,c2,c3,c4,c5,c6,c7]); % Find all six roots
 p = imag(t)==0; % Determine which are the real roots
 t = t(p); % Narrow t vector down to only real roots
 x = a./(a-t)*x0; % Find corresponding x,y,z coordinates
 y = b./(b-t)*y0;
 z = c./(c-t)*z0;
 d = sqrt(min((x-x0).^2+(y-y0).^2+(z-z0).^2)); % Choose the min. dist.

If you want the coordinates of the closest point, you can use the second argument of 'min' combined with p to select those coordinates.

  In theory it would have been possible to solve the problem directly without performing the simplifying translation/rotation, but as you can imagine, the complexity of the result would become unimaginably worse and the required processing to solve the resulting sixth order equation might very well outweigh the processing involved in making the transformation. For that reason (and because it would have taken me a lot more time to solve,) I have chosen this simpler method. I hope you will not have difficulty in working out the translation and rotation stuff.

Roger Stafford

"Robert Phillips" <phillir1@my.erau.edu> wrote in message <ivd95k$6h2$1@newscl01ah.mathworks.com>...
> I have a ellipsoid "defined" at a point E. I have its shape information stored:
>
> elpsd_props{i,1}(:,:) = [centroid xyz coords;
> unit semi-axis vector a;
> unit semi-axis vector b;
> unit semi-axis vector c;
> semi-axis lengths a,b,c]
>
> I have a matrix called astrd_int which holds a massive grid 3-D points. All of these points lie outside of the ellipsoid. I call these interior points (we're inside another shape).
>
> I wish to find the minimum distance from each interior point to the surface of the ellipsoid without the use of generating points on the surface of the ellipsoid.
>
> I already find all of the minimum distances by generating the points on the ellipsoid, calculating distances, and taking the minimum, for each point. I don't wish to do this because the points generated by the ellipsoid function are finite, so if my grid of 3-D points becomes very dense (many many points), it's likely that their minimum distance will "link" to an ellipsoid surface point which is farther away than the point on the real surface. I could control the number of points generated but even so I would rather have a "true" solution, not one that relied on generated surface points.
>
> Is there a better way to do this without using ellipsoid surface points?

I had a funny feeling this was doable using Lagrange
multipliers (but not trivial.) So before I had a chance
to work it out, I did a quick search. Look what popped
up on top:

http://videoprocessing.ucsd.edu/~stanleychan/publication/unpublished/Ellipse.pdf

John

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <ivep55$3b$1@newscl01ah.mathworks.com>...
> I had a funny feeling this was doable using Lagrange
> multipliers (but not trivial.) So before I had a chance
> to work it out, I did a quick search. Look what popped
> up on top:
>
> http://videoprocessing.ucsd.edu/~stanleychan/publication/unpublished/Ellipse.pdf
>
> John
- - - - - - - - - - -
  To John: Hurrah! You are back again! We have missed you, John. It stands to reason it would take something like this interesting kind of problem to tempt you.

  To Robert: You might be interested in reading the cssm thread

 http://www.mathworks.com/matlabcentral/newsreader/view_thread/282755

started May 21, 2010 about the same problem. It inspired Bruno to write the FEX contribution he made about that time. I had completely forgotten about that thread.

Roger Stafford

"Roger Stafford" wrote in message <ivfb8h$rmv$1@newscl01ah.mathworks.com>...
>
>
> started May 21, 2010 about the same problem. It inspired Bruno to write the FEX contribution he made about that time. I had completely forgotten about that thread.
>

Right and on the idea from you Roger, the same you explained above that transforms the problem into polynomial equation.

In the FEX I have work out to so that the function that can handle any dimension and any type of conic hypersurface.

I have an impression that numerical roundoff would have more effect in higher dimension and possibly Newton method could be more effective.

Bruno

Thank you all for the overwhelming amount of help! I'll be reading through all of this and your recommended threads and resources. I will report back with my findings and solutions soon.

I've finally had the time to work more on this. I've gone through step-by-step and obtained the polynomial equation and its seven coefficients (although it seems my coefficients are the negatives of the ones provided by Roger Stafford). I wasn't familiar with how Lagrange multipliers work so it took a bit to absorb it.

Before defining the polynomial's seven coefficients and taking roots etc., I'm simply "interpreting" each ellipsoid as being translated and rotated. In my head it wouldn't make sense to literally transform it because there's technically "nothing" to transform. Am I correct in this thought process?


Also, I'm having difficulty with the roots command (I haven't worked much with the symbolic toolbox). I receive an error:
"??? Undefined function or method 'isfinite' for input
arguments of type 'sym'.

Error in ==> roots at 27
if ~all(isfinite(c))"
Naturally, I have c1,c2,c3,c4,c5,c6,c7 defined as 'sym'. Am I missing a step here? Am I supposed to define the coefficients in some other way?


Finally, I am having trouble understanding the process after obtaining:
d = sqrt(min((x-x0).^2+(y-y0).^2+(z-z0).^2)).
Because the ellipsoid was translated and rotated to the origin (at least interpreted that way, right?), naturally, I have to translate and rotate the solution appropriately (using the centroid E and its rotation matrix [unit vectors E_a; E_b; E_c], right?). An ideal solution would be to have the vector from the point of interest to the ellipsoid's surface point as well as the distance itself

I'm hoping you can validate my steps...

- Obtain 'd' and the corresponding minimum-distance point (xe,ye,ze)
- Calculate the vector from the point of interest x0,y0,z0 to xe,ye,ze
- Rotate the vector by the ellipsoid's rotation matrix [E_a;E_b;E_c] (not sure if I multiply or divide the vector by the matrix?)
- Translate the vector by adding the components of 'E' to the vector
- Then adjust the shortest distance (take the magnitude of the transformed vector?)

Again I really appreciate all the help you've offered!!

"Robert Phillips" <roll24dive7@gmail.com> wrote in message
news:j5e0fb$7es$1@newscl01ah.mathworks.com...
> I've finally had the time to work more on this. I've gone through
> step-by-step and obtained the polynomial equation and its seven
> coefficients (although it seems my coefficients are the negatives of the
> ones provided by Roger Stafford). I wasn't familiar with how Lagrange
> multipliers work so it took a bit to absorb it.
>
> Before defining the polynomial's seven coefficients and taking roots etc.,
> I'm simply "interpreting" each ellipsoid as being translated and rotated.
> In my head it wouldn't make sense to literally transform it because
> there's technically "nothing" to transform. Am I correct in this thought
> process?
>
>
> Also, I'm having difficulty with the roots command (I haven't worked much
> with the symbolic toolbox). I receive an error:
> "??? Undefined function or method 'isfinite' for input
> arguments of type 'sym'.
>
> Error in ==> roots at 27
> if ~all(isfinite(c))"
> Naturally, I have c1,c2,c3,c4,c5,c6,c7 defined as 'sym'. Am I missing a
> step here? Am I supposed to define the coefficients in some other way?

ROOTS is intended to be used with _numeric_ vectors of polynomial
coefficients; it will not work with _symbolic_ vectors of polynomial
coefficients. Since you have a symbolic vector of polynomial coefficients,
use POLY2SYM and pass the resulting polynomial into SOLVE.

*snip*

--
Steve Lord
slord@mathworks.com
To contact Technical Support use the Contact Us link on
http://www.mathworks.com

Thank you for the help; I ended up defining the coefficients as calculated numeric values and then using a for-loop to take ROOTS for each point (I'm calculating minimum distances from points in a Nx3 matrix to the ellipsoid surface).

I've checked my coefficients and roots, and they seem to be correct. However, I am in fact having trouble with the transformation procedure...

N is the number of points that are used for calculating minimum distances to the ellipsoid surface. N = size(astrd_int_2R,1).

a_sqr, b_sqr, & c_sqr are a^2, b^2, & c^2, respectively, and a>b>c.

rm is the "ellipsoid rotation matrix", where its 3 rows are normalized and orthogonal semiaxis vectors, in order a; b; c. So essentially it's [ax,ay,az; bx,by,bz; cx,cy,cz].

ctrd is the "ellipsoid centroid", [Ex,Ey,Ez].

I have the coefficients C1 through C7 (which match exactly to Roger's I believe), all of which are Nx1 size, which are used in my loop. They are ordered such that C1 corresponds to t^6, C2 to t^5, and so on until C7=1 for t^0.

    % For all interior points inside 2*max(R(E-1)):
    for j=1:N
        
        % Calculate t_roots.
        t_roots = roots([C1(j,1),C2(j,1),C3(j,1),C4(j,1),...
                                C5(j,1),C6(j,1),C7(j,1)]);
        
        % Take real t_roots.
        t_roots = t_roots(imag(t_roots)==0);
        
        % Calculate ellipsoid surface xyz coordinates corresponding to
        % t_roots.
        x_e_rel = (a_sqr*astrd_int_2R(j,1))./(a_sqr-t_roots);
        y_e_rel = (b_sqr*astrd_int_2R(j,2))./(b_sqr-t_roots);
        z_e_rel = (c_sqr*astrd_int_2R(j,3))./(c_sqr-t_roots);
        
        % Assign e_rel to xyz coordinates corresponding to point with
        % minimum distance.
        [~,idx_mn] = min(sum([x_e_rel - astrd_int_2R(j,1),...
                                          y_e_rel - astrd_int_2R(j,2),...
                                          z_e_rel - astrd_int_2R(j,3)].^2,2));
        e_rel(j,:) = [x_e_rel(idx_mn),...
                         y_e_rel(idx_mn),...
                         z_e_rel(idx_mn)];
                
    end % End all interior points inside 2*max(R(E-1)).

And I perform a transformation, *which I'm not sure is appropriate / correct*:
e_new = [sum(bsxfun(@times,e_rel,(rm(:,1)')),2) + ctrd(1,1)...
               sum(bsxfun(@times,e_rel,(rm(:,2)')),2) + ctrd(1,2)...
               sum(bsxfun(@times,e_rel,(rm(:,3)')),2) + ctrd(1,3)];

which is (or at least should be) equivalent to:

x_e_new = e_rel(:,1)*rm(1,1) +...
                 e_rel(:,2)*rm(2,1) +...
                 e_rel(:,3)*rm(3,1) + ctrd(1,1);
y_e_new = e_rel(:,1)*rm(1,2) +...
                 e_rel(:,2)*rm(2,2) +...
                 e_rel(:,3)*rm(3,2) + ctrd(1,2);
z_e_new = e_rel(:,1)*rm(1,3) +...
                 e_rel(:,2)*rm(2,3) +...
                 e_rel(:,3)*rm(3,3) + ctrd(1,3);
e_new = [x_e_new, y_e_new, z_e_new];

When I use scatter3 to visualize the e_rel points and e_new points (along with other objects), they look (very) wrong. The e_new points do not seem to lie on the surface of the ellipsoid (whose surface points are properly visualized in the same figure). Also the e_rel and e_new points are very small in magnitude, too.

I figured that, since the ellipsoid surface points are (at least interpreted to be) oriented along the coordinate system's axes and centered at the coordinate system's origin, then it is appropriate to transform the coordinates e_new according to the ellipsoid rotation matrix and centroid.

Is there something that I'm doing wrongly here?

For reference, the real t_roots for the first for-loop iteration are
-444.64
50.1011
29.2212
19.3991
and the t_roots selected for minimum distance was -444.64


"Steven_Lord" <slord@mathworks.com> wrote in message <j5fed5$f6h$1@newscl01ah.mathworks.com>...
>
>
> "Robert Phillips" <roll24dive7@gmail.com> wrote in message
> news:j5e0fb$7es$1@newscl01ah.mathworks.com...
> > I've finally had the time to work more on this. I've gone through
> > step-by-step and obtained the polynomial equation and its seven
> > coefficients (although it seems my coefficients are the negatives of the
> > ones provided by Roger Stafford). I wasn't familiar with how Lagrange
> > multipliers work so it took a bit to absorb it.
> >
> > Before defining the polynomial's seven coefficients and taking roots etc.,
> > I'm simply "interpreting" each ellipsoid as being translated and rotated.
> > In my head it wouldn't make sense to literally transform it because
> > there's technically "nothing" to transform. Am I correct in this thought
> > process?
> >
> >
> > Also, I'm having difficulty with the roots command (I haven't worked much
> > with the symbolic toolbox). I receive an error:
> > "??? Undefined function or method 'isfinite' for input
> > arguments of type 'sym'.
> >
> > Error in ==> roots at 27
> > if ~all(isfinite(c))"
> > Naturally, I have c1,c2,c3,c4,c5,c6,c7 defined as 'sym'. Am I missing a
> > step here? Am I supposed to define the coefficients in some other way?
>
> ROOTS is intended to be used with _numeric_ vectors of polynomial
> coefficients; it will not work with _symbolic_ vectors of polynomial
> coefficients. Since you have a symbolic vector of polynomial coefficients,
> use POLY2SYM and pass the resulting polynomial into SOLVE.
>
> *snip*
>
> --
> Steve Lord
> slord@mathworks.com
> To contact Technical Support use the Contact Us link on
> http://www.mathworks.com

"Robert Phillips" <roll24dive7@gmail.com> wrote in message <j5ld06$blu$1@newscl01ah.mathworks.com>...

>
> Is there something that I'm doing wrongly here?
>

Why don't you use the FEX, at least to check the results against each other?

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <j5lein$gmf$1@newscl01ah.mathworks.com>...
> "Robert Phillips" <roll24dive7@gmail.com> wrote in message <j5ld06$blu$1@newscl01ah.mathworks.com>...
>
> >
> > Is there something that I'm doing wrongly here?
> >
>
> Why don't you use the FEX, at least to check the results against each other?
>
> Bruno

Bruno, I've been reading through the related posts and through your FEX; it's just a bit much for me to absorb at the moment.

But as an update, I reversed the order of my C1 through C7 input arguments for ROOTS, and the problem seems to be solved! I don't understand why COEFFS would return my results in reverse order for use with ROOTS... Maybe I messed up somewhere.

More updates to come.

Thank you,
Robby