The script is to do ellipsoil fitting based on the particle surface voxel cartisian coordinates

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc
clear all

tic

A = xlsread('particle.xlsx');
x = A(:,1);
y = A(:,2);
z = A(:,3);

% Mass Center
XXX = sum(x)/numel(x); %non zero numbers of the surface.
YYY = sum(y)/numel(y);
ZZZ = sum(z)/numel(z);

% Transfer to around origin
x = x - XXX;
y = y - YYY;
z = z - ZZZ;

% Ellipsoid fitting
% fit ellipsoid in the form Ax^2 + By^2 + Cz^2 + 2Dxy + 2Exz + 2Fyz + 2Gx + 2Hy + 2Iz = 1
D = [ x .* x, y .* y, z .* z, 2 * x .* y,2 * x .* z, 2 * y .* z, ...
2 * x, ...
2 * y, ...
2 * z ];
% solve the normal system of equations
v = ( D' * D ) \ ( D' * ones( size( x, 1 ), 1 ) );
% Algebraic form of the Ellipsoid
E = [ v(1) v(4) v(5) v(7);v(4) v(2) v(6) v(8); v(5) v(6) v(3) v(9); ...
v(7) v(8) v(9) -1 ];
% find the center of the ellipsoid
center = -E( 1:3, 1:3 ) \ [ v(7); v(8); v(9) ];
% form the corresponding translation matrix
T = eye( 4 );
T( 4, 1:3 ) = center';
% translate to the center
R = T * E * T';
% solve the eigenproblem
% eigenvectors corresponding to the semi axis directions
[evecs,evals] = eig( R( 1:3, 1:3 ) / -R( 4, 4 ) );
% Three semi-axis
Rad = sqrt( 1 ./ diag( evals ) );

% [ center, radii, evecs, v ] = ellipsoid_fit( [x y z ] );
fprintf( 'Original Mass center: %.3g %.3g %.3g\n', XXX,YYY,ZZZ );
fprintf( 'Final Mass center: 0,0,0\n' );
fprintf( 'Ellipsoid center: %.3g %.3g %.3g\n', center );
fprintf( 'Ellipsoid semi axis : %.3g %.3g %.3g\n', Rad );
fprintf( 'Ellipsoid Eigenvectors :\n' );
fprintf( '%.3g %.3g %.3g\n%.3g %.3g %.3g\n%.3g %.3g %.3g\n', ...
evecs(1), evecs(2), evecs(3), evecs(4), evecs(5), evecs(6), ...
evecs(7), evecs(8), evecs(9) );
fprintf( 'Algebraic form :\n' );
fprintf( '%.3g ', v );
fprintf( '\n' );

% draw data
figure
plot3( x, y, z, '.b' );
axis equal;
axis off;
hold on;

figure
plot3( x, y, z, '.g' );
axis equal;
axis off;
hold on;

%Fitted ellipsoid
maxx = max(x);
maxy = max(y);
maxz = max(z);
minx = min(x);
miny = min(y);
minz = min(z);

[x,y,z] = meshgrid(minx:maxx,miny:maxy,minz:maxz);

Ellipsoid = v(1) *x.*x + v(2) * y.*y + v(3) * z.*z + ...
2*v(4) *x.*y + 2*v(5)*x.*z + 2*v(6) * y.*z + ...
2*v(7) *x + 2*v(8)*y + 2*v(9) * z;
p = patch( isosurface( x, y, z, Ellipsoid, 1 ) );
set( p, 'FaceColor', 'Magenta', 'EdgeColor', 'none' );
view( -70, 40 );
axis vis3d;
camlight;
lighting phong;

toc

 

%%%%

       1. Particle surface voxel cartisian coordinates

     2. Fitted ellipsoid