Population
Image processing library in C++
Watershed 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:

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

Mat3UI8 img;
img.loadFromDirectory("/home/vincent/Desktop/ImageProject/SedimentaryRock/slice/","Data","png");
Rock.gif

Here, you load all png files in the folder /home/vincent/Desktop/ImageProject/SedimentaryRock/slice/ with the basename Data and the extension png 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 :)

rock1.png

This matrix of rock 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();
rock2.png

The grey level inside the pore-space (the blackest phase), the oil (the middle gray phase) and the granular (the lighest phase) is quite homogeneous but you cannot apply the threshold segmenetation due to the interface artefact when you have three phase. You have some voxels at the interface between the blackest phase and the lighest phase with a middle gray. So, if you apply the threshold segmentation for the middle gray phase, these voxels will be assigned badly to this phase. You can see this artefact in the below figure where some blue lines corresponding to the oil phase are located at the interface between the pore phase and the grain phase.

rock3.png
Oil segmentation with interface artefact

For this matrix, we will apply the watershed segmentation as presented in the chapter 3 of Phd thesis http://tel.archives-ouvertes.fr/tel-00516939/en/. The main of this segmentation is the localization of a seed inside each phase with these constraints for each seed:

To create the seeds, you apply the threshold operation in order to respect these previous constraints. To find the threshold values, you can use this manual tool:

int lowvalue,highvalue;
Application::thresholdSelection(imgfilter,lowvalue,highvalue);

When these values are found, you can apply the threshold operation

Mat3UI8 grain= Processing::threshold(imgfilter,155);
Mat3UI8 oil = Processing::threshold(imgfilter,70,110);
oil = Processing::openingRegionGrowing(oil,2);//To remove the interface artefact
Mat3UI8 air = Processing::threshold(imgfilter,0,40);

For the oil seed, you apply the opening operator to remove the interface artefact. Then to merge the seeds in a single label matrix, you do that:

Mat3UI8 seed = Processing::labelMerge(grain,oil);
seed = Processing::labelMerge(seed,air);
Visualization::labelForeground(seed,imgfilter).display();//check the good seed localization
rock4.png
Seed Localization

Then, you create the topographic surface as the magnitude of the gradient matrix. The Deriche's gradient removes the fluctuation at small scales for a smooth gradient.

Mat3UI8 gradient = Processing::gradientMagnitudeDeriche(img,1.5);

Finally, you apply the watershed transformation on the topographic surface with the seed matrix.

Mat3UI8 water = Processing::watershed(seed,gradient);
Visualization::labelAverageColor(water,img).display();
rock5.png

Now, you can comment the crop at the beginning to apply the segmentation procedure on the full image to get this matrix at the end:

Segmentation.gif
Segmented matrix

C++ code

#include"Population.h"//Single header
using namespace pop;//Population namespace
int main(){
Mat3UI8 img;
img.load(POP_PROJECT_SOURCE_DIR+std::string("/image/rock3d.pgm"));
//img.loadFromDirectory("/home/vincent/Desktop/tomo/","tomo2048","png");//if the 3d image is a stack of 2d images
img = img(Vec3I32(0,0,0),Vec3I32(128,128,128));
// img.display();
Mat3UI8 grain= Processing::threshold(imgfilter,155);
Mat3UI8 oil = Processing::threshold(imgfilter,70,110);
oil = Processing::openingRegionGrowing(oil,2);//To remove the interface artefact
Mat3UI8 air = Processing::threshold(imgfilter,0,40);
Mat3UI8 seed = Processing::labelMerge(grain,oil);
seed = Processing::labelMerge(seed,air);
// Visualization::labelForeground(seed,imgfilter).display();//check the good seed localization
Mat3UI8 water = Processing::watershed(seed,gradient);
grain = Processing::labelFromSingleSeed(water,grain);
grain=grain/2;
oil = oil/4;
Scene3d scene;
scene.display();
}

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 img = Mat3UI8()
11 img.load(PathPop+"/image/rock3d.pgm")
12 img = img(Vec3I32(0,0,0),Vec3I32(64,64,64))
13 pde = PDE()
14 imgfilter= pde.nonLinearAnisotropicDiffusionDericheFast(img)
15 app=Application()
16 #app.thresholdSelection(imgfilter)
17 proc = Processing()
18 grain= proc.threshold(imgfilter,155)
19 oil = proc.threshold(imgfilter,70,110)
20 oil = proc.openingRegionGrowing(oil,2)
21 air = proc.threshold(imgfilter,0,40)
22 seed = proc.labelMerge(grain,oil)
23 seed = proc.labelMerge(seed,air)
24 visu =Visualization()
25 #visu.labelForeground(seed,imgfilter).display()
26 gradient = proc.gradientMagnitudeDeriche(img,1.5)
27 water = proc.watershed(seed,gradient)
28 grain = proc.labelFromSingleSeed(water,grain)
29 grain=grain/2
30 oil = proc.labelFromSingleSeed(water,oil)
31 oil = oil/4
32 scene = Scene3d()
33 visu.marchingCube(scene,grain)
34 visu.marchingCube(scene,oil)
35 visu.lineCube(scene,img)
36 scene.display()
37