vol2normalField.cpp

#include <iostream>
#include <iterator>
#include "DGtal/base/Common.h"
#include "DGtal/topology/CanonicDigitalSurfaceEmbedder.h"
#include "DGtal/topology/DigitalSurface.h"
#include "DGtal/topology/DigitalSetBoundary.h"
#include "DGtal/topology/ImplicitDigitalSurface.h"
#include "DGtal/topology/LightImplicitDigitalSurface.h"
#include "DGtal/topology/ExplicitDigitalSurface.h"
#include "DGtal/topology/LightExplicitDigitalSurface.h"
#include "DGtal/graph/BreadthFirstVisitor.h"
#include "DGtal/topology/helpers/FrontierPredicate.h"
#include "DGtal/topology/helpers/BoundaryPredicate.h"
#include "DGtal/graph/CUndirectedSimpleLocalGraph.h"
#include "DGtal/graph/CUndirectedSimpleGraph.h"

#include "DGtal/io/readers/VolReader.h"
#include "DGtal/images/imagesSetsUtils/SetFromImage.h"
#include "DGtal/images/SimpleThresholdForegroundPredicate.h"
#include "DGtal/images/ImageSelector.h"
#include "DGtal/shapes/Shapes.h"
#include "DGtal/helpers/StdDefs.h"
#include "DGtal/kernel/CanonicEmbedder.h"

#include "DGtal/geometry/surfaces/estimation/CNormalVectorEstimator.h"
#include "DGtal/geometry/surfaces/estimation/BasicConvolutionWeights.h"
#include "DGtal/geometry/surfaces/estimation/LocalConvolutionNormalVectorEstimator.h"
#include "DGtal/geometry/surfaces/estimation/DigitalSurfaceEmbedderWithNormalVectorEstimator.h"

#include "CLI11.hpp"


using namespace std;
using namespace DGtal;
using namespace Z3i;

void missingParam ( std::string param )
{
    trace.error() <<" Parameter: "<<param<<" is required..";
    trace.info() <<std::endl;
    exit ( 1 );
}


int main ( int argc, char**argv )
{

    // parse command line CLI ----------------------------------------------
    CLI::App app;
    std::string filename;
    std::string outputFileName;
    unsigned int level {0};
    double sigma {5.0};
    unsigned int neighborhood {10};
    double normExport {1.0};

    app.description("Generates normal vector field from a vol file using DGtal library.\n Typical use example:\n \t vol2normalField[options] --input <volFileName> --o <outputFileName>\n");
    app.add_option("-i,--input,1",filename,"Input vol file.")->required()->check(CLI::ExistingFile);
    app.add_option("-o,--output,2",outputFileName,"Output file.")->required();
    app.add_option("--level,-l",level,"Iso-level for the surface construction (default 0).",true);
    app.add_option("--sigma,-s", sigma,"Sigma parameter of the Gaussian kernel (default 5.0).",true);
    auto expOpt = app.add_flag("--exportOriginAndExtremity", "exports the origin and extremity of the vector fields when exporting the vector field in TXT format (useful to be displayed in other viewer like meshViewer).");
    app.add_option("--vectorsNorm,-N", normExport, "set the norm of the exported vectors in TXT format (when the extremity points are exported with --exportOriginAndExtremity). By using a negative value you will invert the direction of the vectors (default 1.0).",true); 
    app.add_option("--neighborhood,-n", neighborhood,"Size of the neighborhood for the convolution (distance on surfel graph, default 10).",true);  

    app.get_formatter()->column_width(40);
    CLI11_PARSE(app, argc, argv);
    // END parse command line using CLI ----------------------------------------------  

    typedef ImageSelector < Z3i::Domain, unsigned char>::Type Image;
    Image image = VolReader<Image>::importVol ( filename );

    trace.info() <<image<<std::endl;

    functors::SimpleThresholdForegroundPredicate<Image> simplePredicate ( image, level );

    KSpace ks;
    bool space_ok = ks.init ( image.domain().lowerBound(),
                              image.domain().upperBound(), true );
    if ( !space_ok )
    {
        trace.error() << "Error in the Khamisky space construction."<<std::endl;
        return 2;
    }

    typedef SurfelAdjacency<KSpace::dimension> MySurfelAdjacency;
    MySurfelAdjacency surfAdj ( true ); // interior in all directions.

    //Set up digital surface.
    typedef LightImplicitDigitalSurface<KSpace, functors::SimpleThresholdForegroundPredicate<Image>  > MyDigitalSurfaceContainer;
    typedef DigitalSurface<MyDigitalSurfaceContainer> MyDigitalSurface;
    SCell bel = Surfaces<KSpace>::findABel ( ks, simplePredicate );

    MyDigitalSurfaceContainer* ptrSurfContainer =
        new MyDigitalSurfaceContainer ( ks, simplePredicate, surfAdj, bel );
    MyDigitalSurface digSurf ( ptrSurfContainer ); // acquired
    MyDigitalSurface::ConstIterator it = digSurf.begin();

    // Embedder definition
    typedef CanonicDigitalSurfaceEmbedder<MyDigitalSurface> SurfaceEmbedder;
    SurfaceEmbedder surfaceEmbedder ( digSurf );

    //Convolution kernel
    deprecated::GaussianConvolutionWeights < MyDigitalSurface::Size > Gkernel ( sigma );

    //Estimator definition
    typedef deprecated::LocalConvolutionNormalVectorEstimator  < MyDigitalSurface,
                                                     deprecated::GaussianConvolutionWeights< MyDigitalSurface::Size>  > MyGaussianEstimator;
    BOOST_CONCEPT_ASSERT ( ( concepts::CNormalVectorEstimator< MyGaussianEstimator > ) );
    MyGaussianEstimator myNormalEstimatorG ( digSurf, Gkernel );

    // Embedder definition
    typedef DigitalSurfaceEmbedderWithNormalVectorEstimator<SurfaceEmbedder,MyGaussianEstimator> SurfaceEmbedderWithGaussianNormal;
    SurfaceEmbedderWithGaussianNormal mySurfelEmbedderG ( surfaceEmbedder, myNormalEstimatorG );

    // Compute normal vector field and displays it.
    myNormalEstimatorG.init ( 1.0, neighborhood );

    trace.info() << "Generating the NOFF surface "<< std::endl;
    ofstream out2 ( ( outputFileName + ".off" ).c_str() );
    if ( out2.good() )
      digSurf.exportAs3DNOFF ( out2 ,mySurfelEmbedderG );
    out2.close();

    trace.info() << "Generating the polar coordinates file"<< std::endl;
    ofstream out3 ( ( outputFileName + ".txt" ).c_str() );
    if ( out3.good() )
    {
        MyGaussianEstimator::Quantity res;
        for ( MyDigitalSurface::ConstIterator it =digSurf.begin(),
                itend = digSurf.end(); it != itend; ++it )
        {
            res = myNormalEstimatorG.eval ( it );
            //We output Theta - Phi
            out3<< acos ( res [2] ) *180.0/M_PI <<"  " << ( atan2 ( res [1], res [0] ) + M_PI ) *180.0/M_PI;

            if (expOpt->count()>0)
              {
                res *= normExport;
                out3 << " " << mySurfelEmbedderG(*it)[0]
                     << " " << mySurfelEmbedderG(*it)[1]
                     << " " << mySurfelEmbedderG(*it)[2] << " " 
                     <<  mySurfelEmbedderG(*it)[0]+res[0] << " "
                     << mySurfelEmbedderG(*it)[1]+res[1] <<  " "
                     << mySurfelEmbedderG(*it)[2]+res[2];
              
              }

              out3 <<std::endl;
        }
    }
    out3.close();

    return 0;
}


//                                                                           //