Population
Image processing library in C++
Threshold segmentation

A prerequisite of a quantitative matrix analysis is the phase partition (segmentation) that is the transformation of the level matrix to the labelled matrix where each label corresponds to a phase of the material. The error of the quantitative analysis depends on the match between the segmented structure and the real structure. Basically, segmentation is decomposed in two parts:

filtering: to increase the signal to noise ratio without removing some structural features like the presence of cracks. The library includes the classical morphological filters and also the efficient non-linear anisotropic diffusion.
segmentation: to convert the grey-level matrix to a labelled matrix. When the contrast-to-noise ratio between the phases of the material are enough good, the classical threshold segmentation can be used. But, in many cases, the matrix is hampered by an low signal-to-noise ratio, this library includes a methodology based on a growing process starting from seeds on the topographic surface.

When the contrast-to-noise ratio between the phases of the material are enough good, the classical threshold segmentation can be used.

To illustrate the threshold segmentation, we will show you an example with a 3D matrix.
First, you load the 3d matrix:

#include"Population.h"//Single header
using namespace pop;//Population namespace
int main(){
Mat3UI8 img;
img.loadFromDirectory("/home/vincent/Desktop/WorkSegmentation/sand/","500-755","pgm");

Here, you load all pgm files in the folder /home/vincent/Desktop/WorkSegmentation/sand/ with the basename 500-755 and the extension pgm to create the 3D matrix.
In order to create a prototype and calibrate it rapidly, you extract a sub-matrix to reduce the computational time of process:

img = img(Vec3I32(0,0,0),Vec3I32(64,64,64));
img.display();

When the prototype is finished, you comment this line to process the full matrix and you go to the coffee room :)

sand1.png

This matrix of sand material is hampered by noise (see the above matrix). You reduce it by applying the non linear anisotropic diffusion:

Mat3UI8 imgfilter= PDE::nonLinearAnisotropicDiffusion(img);
imgfilter.display();
sand2.png

Now the grey level inside the pore-space and the sand grains is quite homogeneous, you can apply the threshold segmentation. To find the threshold value, you plot the histogram:

Mat2F32 m = Analysis::histogram(img);
DistributionRegularStep d(m);
DistributionDisplay::display(d,d.getXmin(),d.getXmax());
sand3.png

In this histogram, the minimum value between the peaks populated by the pore phase and the granular phase correspond to the argument 155. Therefore you will threshold the matrix with this value:

Mat3UI8 grain= Processing::threshold(imgfilter,155);//or you can apply Processing::thresholdOtsuMethod(imgfilter);
grain.display();
Mat3RGBUI8 color= Visualization::labelForegroundBoundary(grain,img);
color.display();
grain.save("/home/vincent/Desktop/WorkSegmentation/sand/grain.pgm");

To check the agreement between the segmentation and the original matrix, you draw the boundary of the segmented phase on the original matrix.

sand4.png

C++ code

#include"Population.h"//Single header
#include"data/notstable/CharacteristicCluster.h"
using namespace pop;//Population namespace
int main(int argc, char *argv[]){
Mat3UI8 img;
img.load(POP_PROJECT_SOURCE_DIR+std::string("/image/rock3d.pgm"));
img = img(Vec3I32(0,0,0),Vec3I32(50,50,50));
// img.display();
Mat3UI8 imgfilter=Processing::median(img,2,2);
// imgfilter.display();
DistributionRegularStep d(m);
//DistributionDisplay::display(d,d.getXmin(),d.getXmax());
int threshold;
Mat3UI8 grain= Processing::thresholdOtsuMethod(imgfilter,threshold);
// grain.display();
color.display();
grain.save("grain.pgm");
}

Pyton code

1 import sys
2 import os
3 PathPop= ""
4 if os.path.isfile(PathPop+"population.py")==0:
5  print "set the variable PathPop to the path where you compile population, for instance D:\Users/vtariel/Desktop/ANV/Population-build/. This folder must contain population.py"
6  sys.exit(-1)
7 sys.path.append(PathPop)
8 from population import *
9 
10 
11 img = Mat3UI8()
12 img.load(PathPop+"/image/rock3d.pgm")
13 # img.loadFromDirectory("/home/vincent/Desktop/WorkSegmentation/sand/","500-755","pgm") #to load a stack of images
14 img = img(Vec3I32(0,0,0),Vec3I32(64,64,64))
15 img.display()
16 pde = PDE()
17 imgfilter= pde.nonLinearAnisotropicDiffusionDericheFast(img)
18 proc = Processing()
19 grain= proc.thresholdOtsuMethod(imgfilter)
20 grain.display()
21 
22 visu = Visualization()
23 color= visu.labelForegroundBoundary(grain,img)
24 color.display()
25 grain.save(PathPop+"/image/grain.pgm")
26 
27 scene = Scene3d()
28 visu.marchingCube(scene,grain)
29 visu.lineCube(scene,img)
30 scene.display()
31