How to find the maximum and minimum diameter from a set of 3D data points ?
24 views (last 30 days)
Show older comments
Hi,
I have a set of data points (x,y,z) that I have imported on to Matlab. It is of a slightly rotated ellipsoidal shape. I did a 3D scatter plot of these and its shown in the attached screenshot image. The coordinates of these are attached in the text file "cell.txt" where columns 2,3 and 4 representing the X,Y,Z axis respectively.
I applied a mesh and a surface over these points. These are shown in the Matlab code file I included.
I would like to measure the maximum and minimum diameter of this 3D ellipse in different planes (ie XY, YZ, ZX planes) using Matlab.
An example is shown in image "XY plane" where the ellipsoid is in XY plane. I want to measure the maximum and minimum diameter of this ellipsoid.
How could this be carried out?
Thank you.
3 Comments
KSSV
on 5 Feb 2016
As you said you want to find the a, b for the ellipses in XY,YZ and XZ plane. I guess the equation shall work. Say ellipse is in XYplane. pick all points for which Z = 0 ; use equation and solve for a, b.
Answers (3)
John BG
on 5 Feb 2016
There is literature that shows how to find minimum ellypsoid volume, for instance http://compgeom.com/~piyush/papers/emve.pdf
find the smallest ellypsoid WITHIN the samples.
Since you did not mention 'minimum', you want to find the best fit, understanding your samples (are not, I assume your data correct, but) show bit noisy compared to the ellypsoid model, and that some samples are above the ideal surface attempting to approximate your data, and some other samples are below the ellypsoid surface
have a look to MATLAB CENTRAL function for fitting standard volumes to 3D shapes http://uk.mathworks.com/matlabcentral/mlc-downloads/downloads/submissions/48913/versions/4/screenshot.jpg
you may want to spend a bit of time understanding Guiseppe's IEEE paper: Approximation of n-Dimensional Datq using spherical and ellipsoidal primitives http://staff.polito.it/giuseppe.calafiore/Documenti/Papers/Ellipsoidal%20Approximation_TSMC-02.pdf
There is MATLAB CENTRAL function http://uk.mathworks.com/matlabcentral/fileexchange/9542-minimum-volume-enclosing-ellipsoid/content/MinVolEllipse.m
that calculates A, the matrix that defines ellypsoids, copied the function at the end of this text bear in mind that an ellypsoid may be pointing on any direction and centered on any space point, c is the centre: (P_i - c)' * A * (P_i - c) <= 1 and so can your shape, either your -1.35 centering requires a bit more code, or you may want to consider working with this general equation of the ellypsid (P_i - c)' * A * (P_i - c)
from: https://www.wolframalpha.com/input/?i=ellysoid
x(u,v)=a*cos(u)*sin(v) y(u,v)=b*sin(u)*sin(v) z(u,v)=c*cos(v) u within [0:u_step:2*pi] v within [0:v_step:pi]
Volume=4/3*pi*a*b*c
Surface = 2 pi (c^2+b sqrt(a^2-c^2) E(am(sn^(-1)(sqrt(a^2-c^2)/a|(a^2 (b^2-c^2))/(b^2 (a^2-c^2))), (a^2 (b^2-c^2))/(b^2 (a^2-c^2)))|(a^2 (b^2-c^2))/(b^2 (a^2-c^2)))+(b c^2 sn^(-1)(sqrt(a^2-c^2)/a|(a^2 (b^2-c^2))/(b^2 (a^2-c^2))))/sqrt(a^2-c^2))
E:EllipticE am: JacobiSN
The equivalent 2D elipse has parameters that need no further refinement x(t)=a*cos(t) y(t)=b*sin(t)
t within [0 2*pi]
A=pi*a*b
arc length L = 4 a E(1-b^2/a^2)
but Wolfram's equation for the surface of the ellypsoid has operators E() and am() that seem to refer to 'EllipticE' perhaphs for Elliptic enclosure and am() for JacobiSN. Since it is possible to approximate your 3D data with a volume, with the mesh you have used in your supplied example, or directly with the class alphaShape that generates a volume object directly from the 3D points:
Vol1=alphaShape(x1,y1,z1)
plot(Vol1)
So you best shot to find the 3 diameters, assuming your data is centered, is to find an ellypsoid that has as close as possible volume and surface to Vol1
vol_target=volume(Vol1)
x1_min=min(x1);x1_max=max(x1)
y1_min=min(y1);y1_max=max(y1)
z1_min=min(z1);z1_max=max(z1)
1st approach, ellypsoid centered on the median of your samples
assuming you have correctly centered your data, that you probably
haven't bear in mind
a0=abs(x1_max-x1_min)/2
b0=abs(y1_max-y1_min)/2
c0=abs(z1_max-z1_min)/2
volume_approx=4/3*pi*a0*b0*c0
a brute but sometimes effective solution would be:
a0_step=a0/100
a0_range=[a0-a0*.1 a0+a0*.1 ]
b0_step=b0/100
b0_range=[b0-b0*.1 b0+b0*.1 ]
c0_step=c0/100
c0_range=[c0-c0*.1 c0+c0*.1 ]
err=0
for a1=a0_range(1):a0_step:a0_range(2)
for b1=b0_range(1):b0_step:b0_range(2)
for c1=c0_range(1):c0_step:c0_range(2)
v1=4/3*pi*a1*b1*c1
err1=abs(v1-vol_target)/vol_target*100
err=[err err1]
end
end
end
this triple loop is more time consuming than expected, had to cut it
the error err looks like this:
despite the error showing modulated saw within a saw, err shows a mild up-hill slope
find(err(err<0.05))
ans = 2 3 4 5 6 7 8 9 10 11 12 13
find(err(err<0.02))
ans = 2 3 4 5
find(err(err<0.015))
ans = 2 3 4 5
find(err(err<0.013))
ans = []
find(err(err<0.013)) and lower hits the bottom of my sweep window, ans = []
So, an ellypsoid approximating your data, assuming your data correctly centered, that probably is not, would be
a0=abs(x1_max-x1_min)/2
b0=abs(y1_max-y1_min)/2
c0=abs(z1_max-z1_min)/2
a0 = 0.004999500000000
b0 = 0.005000000000000
c0 = 0.004998500000000
Please bear in mind that these lines only aim at helping you on the right direction please read Guiseppe's paper or any other literature you find useful to first center your data, or work with MinVolEllipse.m and then slightly increase it so not all your samples are above the surface of the approximating shape.
Note I have not included another triple loop you may want to test for the surface, then crossing [a,b,c] values with low volume error with [a,b,c] showing low area error values, you may get a bit more precise result.
If you find the ellypsoid surface formula with clearly explained E() and am() operators, to plug it directly to MATLAB script, would you mind to let me know?
Hope it helps, if you find this answer useful please give thumbs up, thanks in advance
John BG
on 5 Feb 2016
Edited: Walter Roberson
on 5 Feb 2016
function [A , c] = MinVolEllipse(P, tolerance)
% [A , c] = MinVolEllipse(P, tolerance)
% Finds the minimum volume enclsing ellipsoid (MVEE) of a set of data
% points stored in matrix P. The following optimization problem is solved:
%
% minimize log(det(A))
% subject to (P_i - c)' * A * (P_i - c) <= 1
%
% in variables A and c, where P_i is the i-th column of the matrix P.
% The solver is based on Khachiyan Algorithm, and the final solution
% is different from the optimal value by the pre-spesified amount of 'tolerance'.
%
% inputs:
%---------
% P : (d x N) dimnesional matrix containing N points in R^d.
% tolerance : error in the solution with respect to the optimal value.
%
% outputs:
%---------
% A : (d x d) matrix of the ellipse equation in the 'center form':
% (x-c)' * A * (x-c) = 1
% c : 'd' dimensional vector as the center of the ellipse.
%
% example:
% --------
% P = rand(5,100);
% [A, c] = MinVolEllipse(P, .01)
%
% To reduce the computation time, work with the boundary points only:
%
% K = convhulln(P');
% K = unique(K(:));
% Q = P(:,K);
% [A, c] = MinVolEllipse(Q, .01)
%
%
% Nima Moshtagh (nima@seas.upenn.edu)
% University of Pennsylvania
%
% December 2005
% UPDATE: Jan 2009
%%%%%%%%%%%%%%%%%%%%%Solving the Dual problem%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% ---------------------------------
% data points
% -----------------------------------
[d N] = size(P);
Q = zeros(d+1,N);
Q(1:d,:) = P(1:d,1:N);
Q(d+1,:) = ones(1,N);
% initializations
% -----------------------------------
count = 1;
err = 1;
u = (1/N) * ones(N,1); % 1st iteration
% Khachiyan Algorithm
% -----------------------------------
while err > tolerance,
X = Q * diag(u) * Q'; % X = \sum_i ( u_i * q_i * q_i') is a (d+1)x(d+1) matrix
M = diag(Q' * inv(X) * Q); % M the diagonal vector of an NxN matrix
[maximum j] = max(M);
step_size = (maximum - d -1)/((d+1)*(maximum-1));
new_u = (1 - step_size)*u ;
new_u(j) = new_u(j) + step_size;
count = count + 1;
err = norm(new_u - u);
u = new_u;
end
%%%%%%%%%%%%%%%%%%%Computing the Ellipse parameters%%%%%%%%%%%%%%%%%%%%%%
% Finds the ellipse equation in the 'center form':
% (x-c)' * A * (x-c) = 1
% It computes a dxd matrix 'A' and a d dimensional vector 'c' as the center
% of the ellipse.
U = diag(u);
% the A matrix for the ellipse
% --------------------------------------------
A = (1/d) * inv(P * U * P' - (P * u)*(P*u)' );
% center of the ellipse
% --------------------------------------------
c = P * u;
2 Comments
Walter Roberson
on 7 Feb 2016
You used %.1f formats. That limits the results to one digit after the decimal point. Use %g instead of %.1f
kamel troudi
on 23 Mar 2016
hey, the function (MinVolEllipse.m) works well with numbers almost points lower at 5000, is there a solution to work with a large number of points. Thank you
0 Comments
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!