This script is to get moment of inertia 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
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

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

     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

% Moment of inertia

clear i
clear j

for i=1:numel(XXX)
   %suppose each voxel mass is 1
   M{i}= ones(1,numel(X{i}));
   mt{i} = sum(M{i});

 for j=1:numel(X{i})
    Ig{i}{j} = [0.0,0.0,0.0;
                     0.0,0.0,0.0;
                     0.0,0.0,0.0];

     x{i}(j)=(X{i}(j,1)-XXX(i));
      y{i}(j)=(Y{i}(j,1)-YYY(i));
      z{i}(j)=(Z{i}(j,1)-ZZZ(i));

          Ig{i}{j} = Ig{i}{j}+[(y{i}(j))^2+(z{i}(j))^2,-x{i}(j)*y{i}(j),-x{i}(j)*z{i}(j);
                           -y{i}(j)*x{i}(j), (x{i}(j))^2+(z{i}(j))^2, -y{i}(j)*z{i}(j);
                           -x{i}(j)*z{i}(j), -y{i}(j)*z{i}(j), (y{i}(j))^2+(x{i}(j))^2];


  end
end

clear i
clear j

% Moment of Inertia of each particle
  for i=1:numel(Ig)
       MomentOfInertia{i} = Ig{1,1}{1,1};
       for j= 2:numel(Ig{i})
          MomentOfInertia{i} = plus(MomentOfInertia{i},Ig{1,i}{1,j});
       end
  end


toc