Pore Analysis

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

namespace T06_02_PoreAnalysis
{
    class Program
    {
        static void Main(string[] args)
        {
            int status = 0;

            try
            {
                // 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.GradientMagnitude3d(imageDenoised,
                                               GradientMagnitude3d.GradientOperator.GAUSSIAN,
                                               new double[] { 0.9, 0.9, 0.9 },
                                               GradientMagnitude3d.FilterMode.RECURSIVE,
                                               GradientMagnitude3d.OutputType.SAME_AS_INPUT,
                                               GradientMagnitude3d.GradientMode.AMPLITUDE_EUCLIDEAN);

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

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

                // Remove pores touching the image border
                ImageView imageInside =
                    Processing.KillBorder3d(imageFilled, KillBorder3d.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");

            }
            catch (Exception error)
            {
                // Print potential exception in the standard output
                System.Console.WriteLine("HelloImageDev exception: " + error.ToString());
                status = -1;
            }

            // ImageDev library finalization
            Initialization.Finish();

            // Check if we must ask for an enter key to close the program
            if (!((args.Length >= 1) && (args[0] == "--no-stop-at-end")))
            {
                System.Console.WriteLine("Press Enter key to close this window.");
                System.Console.ReadKey();
            }

            System.Environment.Exit(status);
        }
    }
}