The script is to get the particle orientation of closed particles

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

clc
clear all

tic

NoOfParticles=dir('ClosedSurface\*.csv');
numberOfcsvs=numel(NoOfParticles);

for i=1:numberOfcsvs
csvName=strcat('ClosedSurface\','surfacevoxel-',num2str(i,'%02i'),'.csv');
I{i}=xlsread(csvName);
x{i} = I{i}(:,1);
y{i} = I{i}(:,2);
z{i} = I{i}(:,3);
% Mass Center
XXX(i) = sum(x{i})/numel(x{i}); %non zero numbers of the surface.
YYY(i) = sum(y{i})/numel(y{i});
ZZZ(i) = sum(z{i})/numel(z{i});

% Transfer to around origin

x{i} = x{i} - XXX(i);
y{i} = y{i} - YYY(i);
z{i} = z{i} - ZZZ(i);

%%% Ellipsoid fitting
% fit ellipsoid in the form Ax^2 + By^2 + Cz^2 + 2Dxy + 2Exz +
% 2Fyz + 2Gx + 2Hy + 2Iz = 1

D{i} = [ x{i} .* x{i}, y{i} .* y{i}, z{i} .* z{i}, 2 * x{i} .* y{i}, ...
2 * x{i} .* z{i}, 2 * y{i} .* z{i}, 2 * x{i}, 2 * y{i}, 2 * z{i} ];

% solve the normal system of equations
v{i} = ( D{i}' * D{i} ) \ ( D{i}' * ones( size( x{i}, 1 ), 1 ) );

% Algebraic form of the ellipsoid
E{i} = [ v{i}(1) v{i}(4) v{i}(5) v{i}(7); ...
v{i}(4) v{i}(2) v{i}(6) v{i}(8); ...
v{i}(5) v{i}(6) v{i}(3) v{i}(9); ...
v{i}(7) v{i}(8) v{i}(9) -1 ];
% find the center of the ellipsoid
center{i} = -E{i}( 1:3, 1:3 ) \ [ v{i}(7); v{i}(8); v{i}(9) ];
% Corresponding translation matrix

T = eye( 4 );
T( 4, 1:3 ) = center{i}';
% translate to the center
R{i} = T * E{i} * T';
% solve the eigenproblem
% eigenvectors corresponding to the semi axis directions
[ evecs{i},evals{i} ] = eig( R{i}( 1:3, 1:3 ) / -R{i}( 4, 4 ) );
% Three semi-axis
Rad{i} = sqrt( 1 ./ diag( evals{i} ) );

% Particle orientations of all the particles
ParticleOrientation(i,:)= evecs{i}(1,:);
% Particle Semi-Axes of all the particles
ParticleSemiAxes(i,:)= Rad{i}';
end

% Plot arrows of particle orientations

for i=1:numberOfcsvs

     x = ParticleOrientation(i,1);
     y = ParticleOrientation(i,2);
     z = ParticleOrientation(i,3);

     plot3([0;x],[0;y],[0;z]);
     hold on
     axis off
     axis equal

     % Plot arrows
     alpha = 0.1;
     beta = 0.3;

     hu = [x-alpha*(x+beta*(y+eps)); x; x-alpha*(x-beta*(y+eps))];
     hv = [y-alpha*(y-beta*(x+eps)); y; y-alpha*(y+beta*(x+eps))];
     hw = [z-alpha*z;z;z-alpha*z];

     plot3(hu(:),hv(:),hw(:));

end

toc

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

     Arrows showing particle orientations of all the read particles