This script is to get the centroid points of each particle.

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

clc
clear all

tic

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

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

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

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

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

threeDarray = logical(threeDarray);

% Remove sample surface
threeDarray = logical(threeDarray(2:(A-1),2:(B-1),2:(C-1)));

%%% GET THE CENTRIOD VOXEL COORDINATES OF EACH PARTICLE

c = numel(threeDarray(1,:,1)); %Columes of one whole 2D images

% Image Ultimate erosion
M = bwulterode(threeDarray,'euclidean');
[row,col] = find(M);

X = row;
Y = mod(col,c);
Z = ceil(col./c);

%%% The commented lines below are considering three different cases that
%%% centroid points can not be obtained by ultimate erosion

% % Unique value for 'col' to erode uneroded surface to lines
% [col,ia,ic] = unique(col);
% X = row(ia);
% Y = mod(col,c);
% Z = ceil(col./c);
%
% %Case 1, if the uneroded line in x-direction
% [YZ,a1,b1] = unique([Y,Z],'rows','stable');
% X = X(a1);
% Y = YZ(:,1);
% Z = YZ(:,2);
%
% %Case 2, if the uneroded line in y-direction
% [XZ,a2,b2] = unique([X,Z],'rows','stable');
% X = XZ(:,1);
% Y = Y(a2);
% Z = XZ(:,2);
%
% %Case 3, if the uneroded line in z-direction
% [XY,a3,b3] = unique([X,Y],'rows','stable');
% X = XY(:,1);
% Y = XY(:,2);
% Z = Z(a3);

toc