ImageDev

Pore Analysis

This example chains several classical image processing algorithms on a three-dimensional digital rocks image. It extracts several features characterizing the pores contained in shale rock volume.
The input image is a 3D sample of shale rock made up of minerals, organic matter, and pores.

<b>Figure 1.</b> Original image with (red) pores, (green) organic-matter grains, and (blue) minerals
Figure 1. Original image with (red) pores, (green) organic-matter grains, and (blue) minerals

First, this example applies the BilateralFilter3d to denoise the input image while preserving the edges of the different materials.

<b>Figure 2.</b> Denoised image
Figure 2. Denoised image

A first thresholding operation is done to isolate the darker objects corresponding to the pores.

<b>Figure 3.</b> Rough binarization of the pores
Figure 3. Rough binarization of the pores

An erosion is performed to isolate the pores from the organic matters and build pore markers.

<b>Figure 4.</b> Eroded pore markers
Figure 4. Eroded pore markers

A second thresholding is applied to detect the minerals.

<b>Figure 5.</b> Rough binarization of the minerals
Figure 5. Rough binarization of the minerals

Then, the AddObjectToLabel algorithm is applied to create a label image that merges pore markers as label 1 with background markers as label 2.

<b>Figure 6.</b> Merged markers: (light blue) pores and (dark blue) organic matters and minerals
Figure 6. Merged markers: (light blue) pores and (dark blue) organic matters and minerals

The denoised image is processed by the GradientOperator3d algorithm to reveal the edges of the different phases.

<b>Figure 7.</b> Gradient of the denoised image
Figure 7. Gradient of the denoised image

This gradient image can be used as a landscape image for the watershed algorithm. MarkerBasedWatershed3d is applied on this image with the merged markers image to detect the separation surfaces between pores and background.
The FillHoles2d algorithm is applied on each 2D slice of the resulting volume to fill the pores. A BorderKill operation is also applied to remove pores that are not entirely visible in the image.
<b>Figure 8.</b> Result after filling pores and rejecting pores intersecting the image borders
Figure 8. Result after filling pores and rejecting pores intersecting the image borders

The image is labeled to assign the same intensity to all voxels belonging to the same connected pore.

<b>Figure 9.</b> Original image and labeled pores
Figure 9. Original image and labeled pores

Some global measurements are computed to extract the number of detected pores and the phase fraction represented by the pores.

- Pore volume fraction = 1.84%
- Pore number = 114


Finally, three measurements are computed for each pore: its volume, the diameter of a sphere of the same volume, and a shape factor indicating if the pore is spherical or not. The corresponding results are displayed in the console for the first pores of the volume.

Pore EquivalentDiameter Volume3d InverseSphericity3d
1 11.43 781.90 2.12
2 13.42 1266.16 1.72
3 6.91 172.88 1.14
4 9.27 558.92 1.89 

#include <ImageDev/Data/ImageViewHelper.h>
#include <ImageDev/Data/NativeMeasurements.h>
#include <ImageDev/ImageDev.h>
#include <ioformat/IOFormat.h>

using namespace imagedev;
using namespace ioformat;
using namespace iolink;

int
main()
{
    // Initialize the ImageDev library if not done
    if ( imagedev::isInitialized() == false )
        imagedev::init();

    // Open a standard tif file and display the image properties
    auto imageInput = readImage( std::string( IMAGEDEVDATA_IMAGES_FOLDER ) + "shale.am" );

    // Preprocess the input with a bilateral filter for denoising
    auto imageDenoised = bilateralFilter3d( imageInput, 3, 3, 3, 20, BilateralFilter3d::BILATERAL );

    // STEP 1: Pore segmentation with the watershed algorithm
    // Threshold the low graylevels (pore markers)
    auto imageLowBin =
        thresholdingByCriterion( imageDenoised, ThresholdingByCriterion::ComparisonCriterion::LESS_THAN, 40 );

    // Slightly erode these markers in order to avoid contact with the high level markers
    auto imageLowEro = erosion3d( imageLowBin, 1, Erosion3d::CONNECTIVITY_26 );

    // Initialize the marker image (label 1)
    auto imageLowLab = convertImage( imageLowEro, ConvertImage::LABEL_16_BIT );

    // Threshold the high graylevels (mineral markers)
    auto imageHighBin =
        thresholdingByCriterion( imageDenoised, ThresholdingByCriterion::ComparisonCriterion::GREATER_THAN, 90 );

    // Merge high intensity markers with label 2
    auto imageMarkers = addObjectToLabel( imageHighBin, imageLowLab, 2 );

    // Extract the gradient to build the relief image for watershed
    auto resultGrad =
        gradientOperator3d( imageDenoised, GradientOperator3d::CANNY, GradientOperator3d::AMPLITUDE_MAXIMUM, 60 );

    // Perform the Watershed to detect pore boundaries
    auto imageWatershed = markerBasedWatershed3d( resultGrad.outputAmplitudeImage,
                                                  imageMarkers,
                                                  MarkerBasedWatershed3d::REPEATABLE,
                                                  MarkerBasedWatershed3d::LINES,
                                                  MarkerBasedWatershed3d::CONNECTIVITY_26 );

    // Fill all pores in each 2D slice of the volume
    setDimensionalInterpretation( imageWatershed, ImageTypeId::IMAGE_SEQUENCE );
    auto imageFilled = fillHoles2d( imageWatershed, FillHoles2d::CONNECTIVITY_4 );
    setDimensionalInterpretation( imageFilled, ImageTypeId::VOLUME );

    // Remove pores touching the image border
    auto imageInside = borderKill( imageFilled, BorderKill::Neighborhood::CONNECTIVITY_26 );

    // Assign labels to the pores
    auto imagePores = labeling3d( imageInside, Labeling3d::LABEL_16_BIT, Labeling3d::CONNECTIVITY_26 );

    // STEP 2: Pore quantification
    // Compute the pore volume percentage
    auto poreVolume = intensityIntegral3d( imageInside );
    std::cout << "  - Pore volume fraction = " << ( 100.0 * poreVolume->volumeFraction() ) << "%" << std::endl;

    // Compute the pore number
    auto poreCount = objectCount( imagePores );
    std::cout << "  - Pore number = " << poreCount.outputMeasurement->count( 0 ) << std::endl;

    // Define the analysis features to be computed
    AnalysisMsr::Ptr analysis = AnalysisMsr::Ptr( new AnalysisMsr() );
    auto diameter = analysis->select( NativeMeasurements::equivalentDiameter );
    auto volume = analysis->select( NativeMeasurements::volume3d );
    auto shape = analysis->select( NativeMeasurements::inverseSphericity3d );

    // Launch the feature extraction on the segmented image
    labelAnalysis( imagePores, NULL, analysis );
    std::cout << "Pore\t" << diameter->name() << "\t" << volume->name() << "\t" << shape->name() << std::endl;
    // Print the analysis results for 5% of the pores
    for ( int i = 0; i < ( int )( analysis->labelCount() / 20 ); i++ )
    {
        std::cout << ( i + 1 ) << "\t\t" << diameter->value( i ) << "\t\t\t" << volume->value( i ) << "\t\t"
                  << shape->value( i ) << std::endl;
    }

    // Save the segmented image with IOFormat
    writeView( imagePores, R"(T06_02_output.tif)" );

    // Close the ImageDev library
    imagedev::finish();

    return 0;
}


See also