How to find the maximum and minimum diameter from a set of 3D data points ?

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.

[KSSV]

As it is an ellipse, it will obey the equation: (x/a)^2+(y/b)^2 = 1. Where a, b are the major and minor radii. By substituting the points (x,y) (points are known to you) cant you frame equations and solve for a, b?

[Meera V]

Thank you Dr Kolukula for your reply.

I understand the equation, however, I need to identify the maximum and minimum points before I could measure the distance using that equation or 3D vectors etc.

The "ellipsoid" is a sphere which has undergone shear and compressive load deformation. I attempted to find the maximum point's coordinates using the scatter plot and calculated the distance. However the minimum diameter is proving to be tricky as the shape has been rotated/deformed due to shear. Therefore I can not identify where they lie in order to get their coordinates.

Would Matlab be able to identify the minimum and maximum points on the deformed ellipsoid?

Thanks.

[KSSV]

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.

Answer 1

[John BG]

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 [e-mail]

[Meera V]

Thank you John,

I will try what you suggested and read the paper.

I understand what you are saying, that it does not form a smooth shape. These data points are from a CAE program which the simulation were run on and they are the surface points of the deformed sphere.

Thanks.

Answer 2

[John BG]

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 ([e-mail])
% 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;

[Meera V]

I came across this previous post where Image Analyst suggested something that could be helpful to me. I changed some parts of the codes to use the X,Y,Z coordinates of my own sphere and to calculate the distance between the points using the 3D vector. The code I used is attached with this post.

However, the code calculates the distance of the furthest point from each of the points. Would it be possible to calculate the furthest distance of the points from the origin.

I.e:

At the moment, the code outputs the results as:

 Point #761 at (-0.0, 0.0) is farthest away from point #1574 at (0.0, -0.0) and is at a distance of 0.010000

• Can we define the "away from point" as the origin? So the distances could be calculated from that central origin to the other points. This should allow me to identify which point is the furthest and which point is the closest. From this I can calculate the minimum and maximum diameter?
• The command window outputs the coordinates of the points in 2D

     (eg _Point #761 at (-0.0, 0.0)_)

Can it be changed so it displays the 3D coordinates? I changed the code to include the Z axis when calculating the distance but can not change the output display in the command window.

• In addition to this, can the number of decimal places be increased in the output answer? Since the distances are very small, it require higher number of decimal places. I tried

format long

but it did not increase the decimal places in the results.

Could the problems I am facing, be resolved?

Ps. I am not sure how I can tag Image Analyst so they can see this post. Can someone tag Image Analyst?

Thank you.

[Walter Roberson]

You used %.1f formats. That limits the results to one digit after the decimal point. Use %g instead of %.1f

Answer 3

[kamel troudi]

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