Granular Geomechanics Group

The script is to get the local void ratio distribution by using delaunay triangulation

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

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 = threeDarray(41:(A-40), 41:(B-40),41:(C-40));
threeDarray = logical(threeDarray);

% GLOBAL VIOD RATIO CALCULATION

n=length(find(threeDarray==0));

GVR=n/((A-80)*(B-80)*(C-80)-n);

% 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 whole 2D images

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

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

NumOfPts=numel(X);

%REMOVE BOUNDARY POINTS TO GET REGULAR TETRAHEDRONS FOR CALCULATING LVR

N = [X(1:(NumOfPts)),Y(1:(NumOfPts)),Z(1:(NumOfPts))];

Tri = delaunayn(N);

figure
tetramesh(Tri,N);
axis equal;
axis off;
camorbit(20,0);

% MAPPING EACH VOXEL TO CARTESIAN COORDINATES

for i=1:(numberOftifs-80)
    [x{i},y{i}]=find(threeDarray(:,:,i));
    z{i}=i*ones(numel(x{1,i}),1);
end

XX=cat(1,x{1,1},x{1,2});
YY=cat(1,y{1,1},y{1,2});
ZZ=cat(1,z{1,1},z{1,2});

for k=3:(numberOftifs-80)
XX=cat(1,XX,x{1,k});
YY=cat(1,YY,y{1,k});
ZZ=cat(1,ZZ,z{1,k});
end

% TO DETERMINE THE POINTS NO. INSIDE EACH TETRAHEDRON

DT = delaunayTriangulation(N);

PtsCoordinates=[XX,YY,ZZ];
ParticleID = pointLocation(DT,PtsCoordinates);
ParticleID(isnan(ParticleID))=[];

ParticleID_Unique = unique(ParticleID);

NumberOfparticle =histc(ParticleID,ParticleID_Unique);
ParticleSpace = NumberOfparticle;

% TO OBTAIN THE VOLUME OF EACH TETRAHEDRON

V = zeros(size(Tri,1),1);

for i=1:size(Tri,1)
    V(i) = 1/6*abs(dot(N(Tri(i,2),:)-N(Tri(i,1),:), ...
    cross(N(Tri(i,3),:)-N(Tri(i,1),:),N(Tri(i,4),:)-N(Tri(i,1),:))));
end

% LOCAL VOID RATIO DISTRIBUTION
ParticleID_Unique = 1:numel(ParticleSpace);

figure

LV_RATIO = (V(ParticleID_Unique)-ParticleSpace)./ParticleSpace;

V_tetraheron = V(ParticleID_Unique);

Glob_Void = sum(ParticleSpace.*LV_RATIO)./sum(ParticleSpace);

LV_RATIO = LV_RATIO(LV_RATIO>0.1);
E = sort(LV_RATIO);
E = E(E<10);
E = E(1:(numel(E)));

emin = min(E);
emax = max(E);
ee=linspace(emin,emax,100);
nn=histc(E,ee);
nn=nn/length(LV_RATIO)*100;
bar(ee,nn);

toc

* Sample images

             1. Delaunay triangulation of Simple cubic

           2. Local void ratio distribution of HCP sample