#include <ImageDev/ImageDev.h>
#include <ioformat/IOFormat.h>
#include <string.h>

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

int
main( int argc, char* argv[] )
{
    int status = 0;

    try
    {
        // 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 = gradientMagnitude3d( imageDenoised,
                                               GradientMagnitude3d::GAUSSIAN,
                                               { 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
        auto imageWatershed = markerBasedWatershed3d( resultGrad,
                                                      imageMarkers,
                                                      MarkerBasedWatershed3d::REPEATABLE,
                                                      MarkerBasedWatershed3d::LINES,
                                                      MarkerBasedWatershed3d::CONNECTIVITY_26 );

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

        // Remove pores touching the image border
        auto imageInside = killBorder3d( imageFilled, KillBorder3d::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, "T06_02_output.tif" );

        std::cout << "This example ran successfully." << std::endl;
    }
    catch ( const imagedev::Exception& error )
    {
        // Print potential exception in the standard output
        std::cerr << "ImageDev exception: " << error.what() << std::endl;
        status = -1;
    }

    // ImageDev library finalization
    imagedev::finish();

    // Check if we must ask for an enter key to close the program
    if ( !( ( argc == 2 ) && strcmp( argv[1], "--no-stop-at-end" ) == 0 ) )
        std::cout << "Press Enter key to close this window." << std::endl, getchar();

    return status;
}
