The script is get the sphericity of each particle.

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

clc
clear all

tic


listOfimages=dir('TenSpheres\*.tif');
numberOftifs=numel(listOfimages);

for i=1:numberOftifs
    imageName=strcat('TenSpheres\','data',num2str(i),'r.tif');
    I{i}=imread(imageName);
end

% 2D stacks to 3D array
ThreeDImage=cat(3,I{1,1},I{1,2});

for i=3:numberOftifs
    ThreeDImage=cat(3,ThreeDImage,I{1,i});
end

[A,B,C]=size(ThreeDImage);

ThreeDImage = logical(ThreeDImage);

% Image Erosion to separate different particles

SE = strel('disk',16.0);
ThreeDImage = imerode(ThreeDImage,SE);

% Get Number of particles

ParticleLabel = bwlabeln(ThreeDImage,6);
ParticleLabelNO = unique(ParticleLabel);
ParticleLabelNO(ParticleLabelNO==0)=[];

% Get coordinates of each particle surface and each particle center

for j=1:numel(ParticleLabelNO)
    ParticleLabel_1 = ParticleLabel;
    ParticleLabel_1(ParticleLabel~=ParticleLabelNO(j)) = 0;% Remove other particle voxels
    DiffPart{j} = ParticleLabel_1; % Number of voxels for each particle
    % Edges{i} = edge(DiffPart{i});% Number of voxels for each particel surface

    for iii=1:C
       Edges{j}(:,:,iii) = edge(DiffPart{j}(:,:,iii),'sobel');
    end
 
    % Obtaining Coordinates of each particle surface
    [row{j},col{j}] = find(Edges{j});

    X{j} = row{j};
    Y{j} = mod(col{j},B);
    Z{j} = ceil(col{j}./B);

    %%% Mass center Via Troditional Method
    XXX(j) = sum(X{j})/nnz(Edges{j}); %non zero numbers of the surface.
    YYY(j) = sum(Y{j})/nnz(Edges{j});
    ZZZ(j) = sum(Z{j})/nnz(Edges{j});
end

% Sphericity of each particle
for j=1:numel(XXX)
     R(j,1) = (3*nnz(DiffPart{j})/4/pi)^(1/3);
     Sphere_area{j} = 4*pi*R(j,1)^2;
     S(j,1) = Sphere_area{j}/nnz(Edges{j})/sqrt(2);
end

%Save sphericities of each particle to a xls file
Sphericity = xlswrite('sphericities.xlsx',S);

toc