Population
Image processing library in C++
Classes

Matrix In -> Measure (2-point correlation function, REV, histogram,...) More...

Classes

class  pop::Analysis
 Analyse a 2D/3D matrix. More...
 

Detailed Description

Matrix In -> Measure (2-point correlation function, REV, histogram,...)

After segmentation, space is partitioned into phases. One challenging problem deals with the ability to describe the geometrical organisation of these phases. Quantitative knowledge of this organisation is important in many fields, for instance, in porous media, for the understanding of the role of geometric confinement in adsorption, condensation, transport, reaction processes. As seen in my PhD, an analysis can be divided in two categories~:

This class provides a collection of algorithms to analyse a 2d/3d matrix. A simple analysis can be as follows:

#include"data/mat/MatN.h"
#include"algorithm/Analysis.h"
using namespace pop;
int main(){
Mat2UI8 lena;
lena.load((std::string(POP_PROJECT_SOURCE_DIR)+"/image/Lena.bmp").c_str());
std::cout<<Analysis::maxValue(lena)<<std::endl;
return 0;
}

Most of the algorithms return a matrix containing the information. To plot this information, you can convert it in DistributionRegularStep and display it as follows:

Mat2F32 mverpore = Analysis::REVPorosity(grain,VecN<3,F32>(grain.getDomain())*0.5,200);
DistributionRegularStep dverpor(mverpore);
dverpor.display();
mverpore.saveAscii("ver.txt");

Because this displayer is not nice, an another solution is to save the ascii data and to open it with gnuplot or matlab.
This code presents a extended analysis of a 3d matrix (to have the 3d visualization, uncomment the lines of code):

Mat3UI8 porespace;
porespace.load("../image/spinodal.pgm");
porespace = pop::Processing::greylevelRemoveEmptyValue(porespace);//the porespace label is now 1 (before 255)
//porespace = porespace(Vec3I32(0,0,0),Vec3I32(50,50,50));//for the first execution, you can test the process in small sample
spinodal.png
Spinodal structure (see pop::Visualization::marchingCubeLevelSet)
//############ Representative Elementary Volume CHECKING ###########################
Mat2F32 mverpore = Analysis::REVPorosity(porespace,VecN<3,F32>(porespace.getDomain())*0.5,200);
mverpore.saveAscii("spinodal_VER.m");
spinodal_VER.png
The sample size is twice time larger than the stablilzation of the porosity measurement -> so REV oke
//############ METRIC ###################
Mat2F32 mhisto = Analysis::histogram(porespace);
std::cout<<"porespace fraction:"<<mhisto(1,1)<<std::endl;
Mat2F32 mchord = Analysis::chord(porespace);
mchord.saveAscii("spinodal_chord.m");
mchord = mchord.deleteCol(1);
DistributionRegularStep dchord_solid(mchord);
std::cout<<"Charateristic length of solid space "<<Statistics::moment(dchord_solid,1,0,300,1)<<std::endl;// \sum_{i=0}^{300} i*d(i)=27.2
spinodal_chord.png
Peak following by exponantial decrease (for the meaning of this signature see Chap 2, p 37. In Handbook of Porous Media. P. Levitz)
Mat2F32 mcorre= Analysis::correlation(porespace);
mcorre.saveAscii("spinodal_corre.m");
spinodal_corr.png
No long-range oscillion=No periodic structure and REV oke as in REVPorosity
Mat3UI8 fdistance;
Mat2F32 mldistance= Analysis::ldistance(porespace,2,fdistance);//euclidean distance only in Population library ;-)
mldistance.saveAscii("spinodal_ldistance.m");
Mat3RGBUI8 dcorregrad=Visualization::labelToRandomRGB(fdistance);//random color
Scene3d scene;
Visualization::plane(scene,dcorregrad,dcorregrad.getDomain()(0)/2,0);
Visualization::plane(scene,dcorregrad,dcorregrad.getDomain()(1)/2,1);
Visualization::plane(scene,dcorregrad,dcorregrad.getDomain()(2)/2,2);
Visualization::lineCube(scene,dcorregrad);
scene.display();
spinodal_ldistance_pore.png
see Pore-Size Probability Density Function in Torquato's book
Mat3UI8 fgranulo;
//granulo of the solid space
Mat2F32 mlgranulo= Analysis::granulometryMatheron(porespace,2,fgranulo);//quite long algorithm in euclidean distance
mlgranulo.saveAscii("spinodal_granulo.m");
Mat3RGBUI8 dcorregrad=Visualization::labelToRandomRGB(fgranulo);//random color
Scene3d scene;
Visualization::plane(scene,dcorregrad,dcorregrad.getDomain()(0)/2,0);
Visualization::plane(scene,dcorregrad,dcorregrad.getDomain()(1)/2,1);
Visualization::plane(scene,dcorregrad,dcorregrad.getDomain()(2)/2,2);
Visualization::lineCube(scene,dcorregrad);
scene.display();
spinodal_granulo.png
see X-ray microtomography characterisation of the changes in statistical homogeneity of an unsaturated sand during imbibitions
Mat2F32 mgeometrical = Analysis::geometricalTortuosity(porespace);
mgeometrical.saveAscii("spinodal_geometrical.m");

We get 1.19465 in 3 directions // see Jeulin's paper estimation of tortuosity and reconstruction of geodesic paths in 3d

//############ TOPOLOGY #######################
F32 euler = Analysis::eulerPoincare(porespace);
std::ofstream out("spinodal_euler.m");
out<<euler;//euler
out.close();
Mat2F32 mpercolationopening = Analysis::percolationOpening(porespace,2);//->charactertic length related to permeability
mpercolationopening.saveAscii("spinodal_percolationopening.m");//output is 6 in three direction-> the structure sill percolates after an opening of ball of radius of 6 but not with a radius of size 7
Mat2F32 mpercolationerosion = Analysis::percolationErosion(porespace,2);
mpercolationerosion.saveAscii("spinodal_percolationerosion.m");//output is 5 in three direction-> the structure sill percolates after an erosion of ball of radius of 5 but not with a radius of size 6
Mat3UI8 porespace_hole= pop::Processing::holeFilling(porespace);
Mat3UI8 skeleton= Analysis::thinningAtConstantTopology(porespace_hole,"../file/topo24.dat");
Scene3d scene;
scene.display();
spinodal_skeleton.png
Topological skeleton
std::pair<Mat3UI8,Mat3UI8> vertex_edge = Analysis::fromSkeletonToVertexAndEdge (skeleton);
Mat3UI32 verteces = pop::Processing::clusterToLabel(vertex_edge.first,0);
Mat3UI32 edges = pop::Processing::clusterToLabel(vertex_edge.second,0);
scene.display();
spinodal_edge.png
Labelisation of the edges of the topological skeleton
int tore;
GraphAdjencyList<Vec3I32> g = Analysis::linkEdgeVertex(verteces,edges,tore);
Scene3d scene;
scene.display();
spinodal_graph.png
Topological graph
std::cout<<euler/g.sizeVertex()<<std::endl;//N_3/alpha_0 normalised topogical characteristic (connectivity number in my phd)
//############ PHYSICAL ###################
Mat2F32 mdiffusion = PDE::randomWalk(porespace);
mdiffusion.saveAscii("spinodal_self_diffusion.m");
spinodal_self_diffusion.png
Coefficient of self diffusion
MatN<3,Vec3F32> velocityfield;
Mat2F32 K(3,3);
VecF32 kx = PDE::permeability(porespace,velocityfield,0,0.05);//permeability in x-direction
VecF32 ky = PDE::permeability(porespace,velocityfield,1,0.05);//permeability in y-direction
VecF32 kz = PDE::permeability(porespace,velocityfield,2,0.05);//permeability in z-direction
//merge the results in the permeability matrix
K.setCol(0,kx);
K.setCol(1,ky);
K.setCol(2,kz);
//display the norm of the last valocity field
Mat3F32 velocityfield_norm(velocityfield);
ForEachDomain3D(x,velocityfield)
{
velocityfield_norm(x)=normValue(velocityfield(x));
}
Mat3RGBUI8 velocityfield_norm_grad= pop::Visualization::labelToRGBGradation(velocityfield_norm);
Scene3d scene;
Visualization::plane(scene,velocityfield_norm_grad,velocityfield_norm_grad.getDomain()(0)/2,0);
Visualization::plane(scene,velocityfield_norm_grad,velocityfield_norm_grad.getDomain()(1)/2,1);
Visualization::plane(scene,velocityfield_norm_grad,velocityfield_norm_grad.getDomain()(2)/2,2);
Visualization::lineCube(scene,velocityfield_norm_grad);
scene.display();
spinodal_permeability.png
Amplitude of the velocity field