Message Boards Message Boards

0
|
10562 Views
|
9 Replies
|
5 Total Likes
View groups...
Share
Share this post:

Measure the diameters of an ellipsoidal shape?

Posted 8 years ago

Hi,

I have a set of data points (x,y,z). It has an ellipsoidal shape and I would like to measure the maximum and minimum diameters of this shape using Mathematica, but I am not sure how to carry this out. Can you please help me?

Thanks.

Attachments:
POSTED BY: SAMM Hill
9 Replies
Posted 8 years ago

Here is another approach. The equation for an ellipse centered at (xc, yc, zc} with semimajor axes {a, b, c} and oriented with respect to the data coordinate system by Euler angles {[Alpha], [Beta], [Gamma]} is described by matrix methods. Given the equation LHS ==1, the sum of (LHS-1)^2 over the data is used as an error function. FindMinimum is used to determine the values of the parameters which minimize the error. This determines not only the semimajor axes, but also the center and orientation of the best-fit ellipsoid.

(* Equation for an ellipse centered at (xc,yc,zc} with semimajor axes \
{a,b,c} and oriented with respect to the data coordinate system by \
Euler angles {\[Alpha],\[Beta],\[Gamma]} is described by matrix \
methods *)

(* rotation of semimajor axes *)
rot = EulerMatrix[{\[Alpha], \[Beta], \[Gamma]}];

(* the transpose *)
rotT = Transpose@rot;

(* diagonalized form *)
amat = DiagonalMatrix[{1/a^2, 1/b^2, 1/c^2}];

(* coordinate vectors *)
cmat = {{x - xc}, {y - yc}, {z - zc}};

(* the transpose *)
cmatT = Transpose@cmat;

(* this forms the left hand side of LHS \[Equal]1 *)
lhs = Simplify[cmatT.rotT.amat.rot.cmat][[1, 1]];
In[10]:= (* the square error from LHS \[Equal]1 for a single {x,y,z} \
point and given parameter set *)
squareError[{x_, y_, z_}, xc_, yc_, zc_, a_, b_, 
   c_, \[Alpha]_, \[Beta]_, \[Gamma]_] = (lhs - 1)^2;

(* find the parameter set which minimizes the error sum over the data \
*)
{residual, fit} = FindMinimum[
  {
   Plus @@ (squareError[#, xc, yc, zc, a, b, 
        c, \[Alpha], \[Beta], \[Gamma]] &) /@ data1,
   0 <= \[Alpha] < Pi, 0 <= \[Beta] < Pi, 0 <= \[Gamma] < Pi
   },
  {{xc, 2}, {yc, 0}, {zc, 0}, {a, 2}, {b, 2}, {c, 2}, {\[Alpha], 
    0}, {\[Beta], 0}, {\[Gamma], 0}}]

(* Out[11]= {0.044751062, {xc -> 2.0074857, yc -> 0.15231692, zc -> -0.10929044, a -> 4.311946, b -> 3.5987792, c -> 3.5830354, [Alpha] -> 0.00068328857, [Beta] -> 0.31168585, [Gamma] -> 0.85806301}} *)

In[12]:= plot = 
 Show[ListPointPlot3D[data1, AxesLabel -> {"X", "Y", "Z"}, 
   PlotStyle -> Red, BoxRatios -> Automatic], 
  ContourPlot3D[
   Evaluate[1 == lhs /. fit], {x, -5, 6}, {y, -5, 5}, {z, -5, 5}, 
   ContourStyle -> Opacity[.3], Mesh -> None]]

enter image description here

Attachments:
POSTED BY: David Keith
Posted 8 years ago

And an improved visualization:

lines = data1 /. {x_, y_, z_} -> Line[{{xc, yc, zc}, {x, y, z}}] /. 
   fit;

p2 = Show[plot, Graphics3D[lines]]

enter image description here

POSTED BY: David Keith

I'm not sure what BoundedRegion does internally but it has a "FastEllipse" form specification. This is not a 'fit' as Daniel does (right?). i.e. with a fit you have data that is above/below the fit (or inside/outside) but this does not do that, it is always 'outside' of the data.

Depending what you are after this might be a alternative/better solution.

POSTED BY: Sander Huisman
Posted 8 years ago

Geometric tool states least squares fitting method for ellipsoid with documentation and C code. https://www.geometrictools.com/Documentation/LeastSquaresFitting.pdf

POSTED BY: Frederick Wu

Yes, but those approaches are sometimes harder than need be (using nonlinear fitting where there is a formulation for a linear LS fit).

POSTED BY: Daniel Lichtblau

One way to approximate it is to get a best-fit ellipse by singular values decomposition. Since only the stretch factors are of interest we just compute the singular values. If I recall correctly it is their square roots that give the stretch. The method below requires that we center at the origin, hence the subtracting of the mean.

dataCentered = Map[# - Mean[data1] &, data1];
Sqrt[SingularValueList[dataCentered]]

(* Out[34]= {4.5038878147, 4.19597091008, 4.18475712383} *)

To check this we normalize and see if we get something that looks like the unit sphere.

dataNormed = Map[#/Sqrt[svals] &, dataCentered];
ListPointPlot3D[dataNormed, BoxRatios -> {1, 1, 1}]

enter image description here

POSTED BY: Daniel Lichtblau
Posted 8 years ago

Hi I know this is not an efficient way to do it but I did this way. Note that data1=[{}]; is wrong notation, I changed to data1={{}};

 In[72]:= data3 = 
      Table[EuclideanDistance[data1[[j]], data1[[i]]], {i, 
        Length@data1}, {j, Length@data1}];
In[73]:= Max@data3

Out[73]= 8.50821

In[74]:= Position[data3, Max@data3]

Out[74]= {{5, 68}, {68, 5}}

In[75]:= data1[[5]]

Out[75]= {-1.47388, 2.66873, -0.915369}

In[76]:= data1[[68]]

Out[76]= {5.29808, -1.95733, 1.34955}



Show[HighlightMesh[ConvexHullMesh[data1], Style[2, Opacity[0.5]]], 
 Graphics3D[{Red, Line[{data1[[5]], data1[[68]]}]}], 
 Graphics3D[Point[data1]]]
Show[ListPointPlot3D[data1], 
 Graphics3D[{Red, Line[{data1[[5]], data1[[68]]}]}]]
POSTED BY: Okkes Dulgerci
Posted 8 years ago

I find your way is easiest to follow as I am only looking for an approximate length. I am trying to figure out now how to find the diameter of the other dimension of the ellipsoid.

POSTED BY: SAMM Hill
Posted 8 years ago

Hi, It is hard for me to find it. What I can suggest is you need to split your data in two piece along the line (assume you have a plane that slice the ellipsoid in two half). Then find max the distance between two data. and similarly for other dimension..

BTW I modified @Marco Thiel ' s the code for this problem and here it is.

farpoints = 
  Take[data1[[#]] &@Flatten@Position[#, Max[#]] &@ DistanceMatrix[data1], 2];
distance = EuclideanDistance[First@farpoints , Last@farpoints ];
Show[ListPointPlot3D[{data1}], 
 Graphics3D[{{PointSize[Large], Red, 
    Point[#] & /@farpoints }, {Thick, Blue, Line[  farpoints ]}}], 
 ImageSize -> 500, 
 PlotLabel -> Framed@Style["distance = " <> ToString[distance], 20]]
POSTED BY: Okkes Dulgerci
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract