# Grain(Particle) Size Distribution (PSD) Analysis with Mathematica

I would like to do an analysis of the grain size distribution of a Monte Carlo Grain Growth simulation I implemented based on the nice example by Rituraj Nandan (see Link). The result of this calculation is a 2D integer matrix, where every point represents a certain grain orientation. Thus a cluster of similiar values represents a grain. Converting this information into a graphical representation by assigning a similiar color to every similiar integer gives an image which shows the grain distribution of the simulation clearly. Here is an example of such an simulation result (I converted the initial implementation in Java into Mathematica syntax): The code used for generating this image is as follows (see above link to the java implementation for the details of the physics behind it):

meshgen[xmax_:0.01,ymax:0.01,nx_:60,ny_:60]:=Module[{hx,hy},
hx=xmax/nx;
hy=ymax/ny;
Return[{{xmax,ymax},{nx,ny},{hx,hy},{Table[i hx,{i,0,nx-1}],Table[j hy,{j,0,ny-1}]}}];
];

temperature[mesh_,torch_,vel_]:=Module[{x,y,xc,yc,nx,ny,rad,al,t,
q0=500(*net power*),
d=0.002(*thickness of plate*),
\[Rho]=7800(*density*),
cp=600(*specifc heat capacity*),
k=45(*thermal conductivity*),
t0=298.0(*Initial temperature*)},(*quasi-steady thermal profile*)
al=k/(\[Rho] cp);(*thermal diffusivity*)
{nx,ny}=mesh[];
{x,y}=mesh[[-1]];
{xc,yc}=torch;
rad=Outer[Norm[{#1,#2}]&,x-xc,y-yc];
t=Table[t0+q0/(d 2 \[Pi] k) E^(-(1/2)vel (x[[i]]-xc)/al) BesselK[0,1/2 vel rad[[i,j]]/al],{j,1,ny},{i,1,nx}];
Return[t];
];

tmc[mesh_,t_,vel_]:=Module[
{K1=1.01,(*Model constant*)
n1=0.42,(*Model constant*)
A=1,(*Accommodation probability*)
\[CapitalNu]=1.53 10^19,(*Average number of atoms per unit area at the grain boundary in #/m\.b2*)
Na=6.022 10^23,(*Avogadro's number*)
Vm=7.11 10^-6, (* ??? *)
h=6.62 10^-34(*Planck's constant m\.b2kg/s*),
R=8.314, (*Gas Constant J/molK*)
\[CapitalDelta]Sf=9.48, (*Activation entropy of grain boundary migration J/molK*)
Q=93400, (*Activation enthalpy for grain boundary migration in J/mol*)
\[Gamma]=1.77, (*Grain boundary energy J/m\.b2*)
\[Lambda], (*Grid Size*)
L0,  (*Initial Grain Size*)
dt,time,expt,C1,x,nx},(*monte-carlo time at each location (i,j)*)
x=mesh[[-1,1]];
nx=mesh[[2,1]];
\[Lambda]=mesh[[3,1]];
L0=mesh[[3,1]];
C1=(4 \[Gamma] A \[CapitalNu] Vm^2)/(Na^2 h) E^(\[CapitalDelta]Sf/R);
time=x/vel;
dt=Append[Differences[time],time[[-1]]];
expt=Table[Exp[-Q/(R*t[[;;,i]])]*dt[[i]],{i,1,nx}];
expt=Reverse[Accumulate[Reverse[expt]]];
Return[(Sqrt[C1 expt+L0^2]/(K1 \[Lambda])-1/K1)^(1/n1)];
];

energypart[indx_,{nx_,ny_},reflect_:True]:=
Module[{ii,jj,it,ib,jl,jr},(*get part specification of nearest neighbours of location (ii,jj). If reflect is set True then at the array boundary reflecting boundary conditions are assumed. Otherwise periodic boundary conditions are used*)
ii=indx[];jj=indx[];
(*rows of top and bottom neighbours*)
it=If[reflect,If[ii==1,1,ii-1],If[ii==1,ny,ii-1]];
ib=If[reflect,If[ii==ny,ny,ii+1],If[ii==ny,1,ii+1]];
(*rows of left and right neighbours*)
(*rows of top and bottom neighbours*)
jl=If[reflect,If[jj==1,1,jj-1],If[jj==1,nx,jj-1]];
jr=If[reflect,If[jj==nx,nx,jj+1],If[jj==nx,1,jj+1]];
{it;;ib,jl;;jr}
];

pottsenergy[osub_]:=Total[Map[1-KroneckerDelta[#,osub[[2,2]]]&,osub,{2}],2];(*Potts-Function of 3x3 Array with relevant grain in the center at {2,2}. The total of this array gives the Energy of the grain at {2,2}*)

energy[o_,indx_,nmax_,reflect_:True]:=pottsenergy[o[[Sequence@@energypart[indx,nmax,reflect]]]];(*Select relevant block with center element and all neighbours and apply Potts-Function to it. The result of this function gives the Energy of the grain at position indx*)

RemoveNoise[frame_,limit_:2]:=Module[{emap,nx,ny,result,i},
{nx,ny}=Dimensions[frame];(*Set all grain size orientations for grains which just consist of one or two pixels (set by limit) (E=8) equal to -1 (with the color function defined below they will be set to black*)
emap=MapIndexed[energy[frame,#2,{nx,ny}]&,frame,{2}];
i=5;
emap[[1,;;]]=Nest[ReplacePart[#,Position[#,i--]->8]&,emap[[1,;;]],limit];(*Correct for lower number of neighbours at boundaries*)
i=5;
emap[[nx,;;]]=Nest[ReplacePart[#,Position[#,i--]->8]&,emap[[nx,;;]],limit];
i=5;
emap[[;;,ny]]=Nest[ReplacePart[#,Position[#,i--]->8]&,emap[[;;,ny]],limit];
i=5;
emap[[;;,1]]=Nest[ReplacePart[#,Position[#,i--]->8]&,emap[[;;,1]],limit];
i=8;
result=Nest[ReplacePart[#,Position[emap,i--]->-1]&,frame,limit];
Return[result];
];

grainsize[mesh_,o_]:=Module[{dist=Table[0,{10}],avg,size,sum=0,num=0,nx,ny,hx,hy,i,j,k,xmax,ymax},(*measurement of grain size distribution and average grain size*)
{xmax,ymax}=mesh[];
{nx,ny}=mesh[];
{hx,hy}=mesh[];
For[i=1,i<=nx,i+=5,(*measurement of average grain size using mean lineal intercept method*)
For[j=1,j<=ny,j++,
If[o[[i,j]]!=o[[i,Mod[j+ny-1,ny,1]]],(*periodic boundary*)
num+=1
];
]
];
avg=nx ymax/num;(*average grain size*)
For[i=1,i<=nx,i++,(*grain size distribution*)
For[j=1,j<=ny,j++,
If[o[[i,j]]==o[[i,Mod[j+ny-1,ny,1]]],
sum+=1
,
size=sum hy/avg;
sum=0;
For[k=0,k<10,k++,
If[(size<If[k<9,k+1,1000]/2)&&size>=k/2,dist[[k+1]]+=1];
];
];
];
];
Print["Grain size distribution : "];
For[i=0,i<10,i++,Print[i/2//N," - ",If[i<9,i+1,1000]/2//N,"    ",dist[[i+1]]]];
Print["Average grain size : ",avg];
Return[avg];
];


The above routine grainsize is the one still to be optimized since there I did not do much conversion so far.

Finally with the following statements the Graph above can be generated:

mesh = meshgen[0.01, 0.01, 500, 500];
temp = temperature[mesh, {0.0070001, yc = 0.0050001}, 0.005];
tmcs = tmc[mesh, temp, 0.005];
o = Table[RandomInteger[{0, 31}], {mesh[[2, 1]]}, {mesh[[2, 2]]}];
frames = {};
If [! FileExistsQ [NotebookDirectory [] <> "Graingrowth.dmp"],
For[i = 1, i < 250, i++,
tempc = PrintTemporary["Frame " <> ToString[i] <> "/250"];
If[i > 1,
tempp = PrintTemporary[
ArrayPlot[RemoveNoise[frames[[-1]]], PlotRange -> {-1, 31},
ColorFunction -> Function[v, If[v == 0, Black, Hue[v]]]]]];
on = potts[mesh, o, tmcs, 5000000];
AppendTo[frames, on\[Transpose]];
o = on;
If[i > 1, NotebookDelete[tempp]];
NotebookDelete[tempc];
];
AppendTo[frames, o\[Transpose]];
DumpSave[NotebookDirectory[] <> "Graingrowth.dmp", frames];
,
Get [NotebookDirectory[] <> "Graingrowth.dmp"];
];
ArrayPlot[RemoveNoise[frames[[-1]]], PlotRange -> {-1, 31}, ColorFunction -> Function[v, If[v == 0, Black, Hue[v]]]]
Export[NotebookDirectory[] <> "FinalGrainStructure.png", %];ColorFunction -> Function[v, If[v == 0, Black, Hue[v]]]]


Of course there are multiple sources on the web of how to do PSD analysis, however I did not find any details about the potential of using Mathematica for this application. There is a very interesting Reference in the Wolfram Library, but it could not manage to get a hard or soft copy of this work.

I guess the statistical routines within Mathematica would make it perfectly suitable for such an analysis.

Any ideas/help/advice?

• Colorize@dataMatrix might do something close to what you want for the visualisation and perhaps MorphlogicalComponents may also be of use. – image_doctor Nov 10 '12 at 9:44
• "I converted the initial implementation in Java into Mathematica syntax..." - could you maybe include this conversion you speak of in your post? – J. M.'s ennui Nov 10 '12 at 10:05
• @image_doctor: Thanks for the inputs! MorphologicalComponents will come in handy for further statistical evaluation on the extracted grain boundaries I suppose. With respect to Colorize I'll check if this routine is fastern than my ColorFunction approach. Do you have any inputs with respect to the statistical side of the PSD analysis? – Rainer Nov 10 '12 at 13:01
• @J.M.: I'll post the code soon. The current version is very different from the original JAVA implementation, because I tried to optimize the code for speed. Initially the mathematica code was around 50x slower than the java code. Now I guess it is just a few times slower. – Rainer Nov 10 '12 at 13:04
• You might also be interested in DeleteSmallComponents or ComponentMeasurements for removing small grains. – image_doctor Nov 11 '12 at 10:02

## 1 Answer

UPDATE: I found the answer with the help of the inputs by image_doctor and the following post:

As an example I use the final, noise-corrected image shown above for further processing.

With the help of ComponentMeasurements I first separate all grains with same orientation (color) from all grains with different colors.

directedgrains = ComponentMeasurements[img, "Mask"][[;; , 2]];


The part specification just selects the sparse arrays giving the shapes of the extracted grains.

Now with MorphologicalComponents I separate all grains with same orientation from each other and thread over the list directedgrains generated by ComponentMeasurement

subgrains = MorphologicalComponents[Normal[#]] & /@ directedgrains;


With the structure subgrains I can now do all sort of statistical analysis like

Histogram[Flatten[(ComponentMeasurements[#, "Circularity"] & /@subgrains)[[;; , ;; , 2]]]]


which yields a histogram of the circularity of the grains or the "Most Wanted" Information of the PSD with

Histogram[2 Flatten[(ComponentMeasurements[#, "EquivalentDiskRadius"] & /@subgrains)[[;; , ;; , 2]]]]


which finally yields the Histogram According to theory a lognormal distribution is expected. Ignoring the first bar in the histogram, which is most probably due to remaining statistical noise, the distribution looks rather similiar to a lognormal shape.

Thanks again to image_doctor to for the very useful hints!