The script is to get the cartisian coordinates of the mass center of each particle by inporting the binary images.

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

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',8.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
ParticleLabel_1 = ParticleLabel;

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);

     % Obtaining Coordinates of each particle centroid
     Sum_X{j} = sum(X{j});
     Sum_Y{j} = sum(Y{j});
     Sum_Z{j} = sum(Z{j});

     % Mass center
     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

% Write center of mass coordinates to xls file
COM = xlswrite('CenterOfMass.xlsx',[XXX',YYY',ZZZ']);

toc