Pore Analysis

using System;
using ImageDev;
using IOLink;
using IOFormat;

namespace T06_02_PoreAnalysis
{
    class Program
    {
        static void Main( string[] args )
        {
            // Initialize the ImageDev library if not done
            if ( Initialization.IsInitialized() == false )
                Initialization.Init();

            // Open a standard tif file and display the image properties
            ImageView imageInput = IOFormat.ViewIO.ReadImage( @"Data/images/shale.am" ) as ImageView;

            // Preprocess the input with a bilateral filter for denoising
            ImageView imageDenoised =
                Processing.BilateralFilter3d( imageInput, 3, 3, 3, 20, BilateralFilter3d.FilterMode.BILATERAL )
                    as ImageView;

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

            // Slightly erode these markers in order to avoid contact with the high level markers
            ImageView imageLowEro =
                Processing.Erosion3d( imageLowBin, 1, Erosion3d.Neighborhood.CONNECTIVITY_26 ) as ImageView;

            // Initialize the marker image (label 1)
            ImageView imageLowLab =
                Processing.ConvertImage( imageLowEro, ConvertImage.OutputType.LABEL_16_BIT ) as ImageView;

            // Threshold the high graylevels (mineral markers)
            ImageView imageHighBin = Processing.ThresholdingByCriterion(
                imageDenoised, ThresholdingByCriterion.ComparisonCriterion.GREATER_THAN, 90 ) as ImageView;

            // Merge high intensity markers with label 2
            ImageView imageMarkers = Processing.AddObjectToLabel( imageHighBin, imageLowLab, 2 ) as ImageView;

            // Extract the gradient to build the relief image for watershed
            var resultGrad = Processing.GradientOperator3d( imageDenoised,
                                                            GradientOperator3d.GradientOperator.CANNY,
                                                            GradientOperator3d.GradientMode.AMPLITUDE_MAXIMUM,
                                                            60 );

            // Perform the Watershed to detect pore boundaries
            ImageView imageWatershed =
                Processing.MarkerBasedWatershed3d( resultGrad.outputAmplitudeImage,
                                                   imageMarkers,
                                                   MarkerBasedWatershed3d.AlgorithmMode.REPEATABLE,
                                                   MarkerBasedWatershed3d.OutputType.LINES,
                                                   MarkerBasedWatershed3d.Neighborhood.CONNECTIVITY_26 ) as ImageView;

            // Fill all pores in each 2D slice of the volume
            Data.SetDimensionalInterpretation( imageWatershed, ImageTypeId.IMAGE_SEQUENCE );
            ImageView imageFilled = Processing.FillHoles2d( imageWatershed, FillHoles2d.Neighborhood.CONNECTIVITY_4 );
            Data.SetDimensionalInterpretation( imageFilled, ImageTypeId.VOLUME );

            // Remove pores touching the image border
            ImageView imageInside =
                Processing.BorderKill( imageFilled, BorderKill.Neighborhood.CONNECTIVITY_26 ) as ImageView;

            // Assign labels to the pores
            ImageView imagePores = Processing.Labeling3d( imageInside,
                                                          Labeling3d.LabelType.LABEL_16_BIT,
                                                          Labeling3d.Neighborhood.CONNECTIVITY_26 ) as ImageView;

            // STEP 2: Pore quantification
            // Compute the pore volume percentage
            var poreVolume = Processing.IntensityIntegral3d( imageInside );
            Console.WriteLine(
                "  - Pore volume fraction = " + ( 100.0 * poreVolume.volumeFraction() ).ToString( "0.00" ) + " %" );

            // Compute the pore number
            var poreCount = Processing.ObjectCount( imagePores );
            Console.WriteLine( "  - Pore number = " + poreCount.outputMeasurement.count( 0 ) );

            // Define the analysis features to be computed
            AnalysisMsr analysis = new AnalysisMsr();
            var diameter = analysis.Select( NativeMeasurements.EquivalentDiameter );
            var volume = analysis.Select( NativeMeasurements.Volume3d );
            var shape = analysis.Select( NativeMeasurements.InverseSphericity3d );

            // Launch the feature extraction on the segmented image
            Processing.LabelAnalysis( imagePores, null, analysis );
            Console.WriteLine("Pore\t" + diameter.Name() + "\t" + volume.Name() + "\t" + shape.Name() );
            // Print the analysis results for 5% of the pores
            for ( int i = 0; i < ( int )( analysis.LabelCount() / 20 ); i++ )
            {
                Console.WriteLine( ( i + 1 ) + "\t\t" + diameter.Value( i ).ToString( "0.00" ) + "\t\t\t" +
                                   volume.Value( i ).ToString( "0.00" ) + "\t\t" +
                                   shape.Value( i ).ToString( "0.00" ) );
            }

            // Save the segmented image with IOFormat
            IOFormat.ViewIO.WriteView( imagePores, @"T06_02_output.tif" );

            // Close the ImageDev library
            Initialization.Finish();
        }
    }
}