curvatureMCMS.cpp

#include <iostream>

#include "DGtal/base/Common.h"

#include "DGtal/shapes/ShapeFactory.h"
#include "DGtal/shapes/Shapes.h"
#include "DGtal/helpers/StdDefs.h"
#include "DGtal/topology/helpers/Surfaces.h"

//Grid curve
#include "DGtal/geometry/curves/FreemanChain.h"
#include "DGtal/geometry/curves/GridCurve.h"

//Estimators
#include "DGtal/geometry/curves/ArithmeticalDSSComputer.h"
#include "DGtal/geometry/curves/estimation/MostCenteredMaximalSegmentEstimator.h"
#include "DGtal/geometry/curves/estimation/SegmentComputerEstimators.h"

#include "CLI11.hpp"

#include <vector>
#include <string>
#include <iomanip>


using namespace DGtal;
using namespace std;


/*
  Curvature estimation from the length of the most centered maximal segments
  @param h grid step
  @param itb begin iterator on points
  @param ite end iterator on points
  @param ito output iterator on curvature
 */
template<typename I, typename O>
void estimationFromLength( double h, const I& itb, const I& ite, const O& ito )
{
  typedef ArithmeticalDSSComputer<I,int,4> SegmentComputer;
  SegmentComputer sc;
  typedef CurvatureFromDSSLengthEstimator<SegmentComputer> SCEstimator;
  SCEstimator f; 
  typedef MostCenteredMaximalSegmentEstimator<SegmentComputer,SCEstimator> Estimator;
  Estimator CurvatureEstimator(sc, f);

  CurvatureEstimator.eval( itb, ite, ito, h ); 
}

/*
  Curvature estimation from the length and width
  of the most centered maximal segments
  @param h grid step
  @param itb begin iterator on points
  @param ite end iterator on points
  @param ito output iterator on curvature
 */
template<typename I, typename O>
void estimationFromLengthAndWidth( double h, const I& itb, const I& ite, const O& ito )
{
  typedef ArithmeticalDSSComputer<I,int,4> SegmentComputer;
  SegmentComputer sc;
  typedef CurvatureFromDSSEstimator<SegmentComputer> SCEstimator;
  SCEstimator f; 
  typedef MostCenteredMaximalSegmentEstimator<SegmentComputer,SCEstimator> Estimator;
  Estimator CurvatureEstimator(sc, f);

  CurvatureEstimator.eval( itb, ite, ito, h ); 
}



int main( int argc, char** argv )
{
  // parse command line CLI ----------------------------------------------
  CLI::App app;
  string fileName;
  double h {1.0};

  app.description("Estimates curvature using length of most centered segment computers.\n Typical use example:\n \t curvatureMCMS [options] --input  <fileName>\n");
  app.add_option("-i,--input, 1",fileName,"Input FreemanChain file name")->required()->check(CLI::ExistingFile);
  app.add_option("--GridStep", h, "Grid step",true);
  
  app.get_formatter()->column_width(40);
  CLI11_PARSE(app, argc, argv);
  // END parse command line using CLI ----------------------------------------------  

  typedef Z2i::Space Space; 
  typedef Space::Point Point; 
  typedef Space::Integer Integer;  
  typedef Z2i::KSpace KSpace; 
  typedef FreemanChain<Integer> FreemanChain; 

  vector< FreemanChain > vectFcs =  PointListReader< Point >::getFreemanChainsFromFile<Integer> (fileName);  

  for(unsigned int i=0; i<vectFcs.size(); i++){

    // Freeman chain
    FreemanChain fc = vectFcs.at(i); 
    // Create GridCurve
    GridCurve<> gridcurve;
    gridcurve.initFromPointsRange( fc.begin(), fc.end() );

    cout << "# grid curve " << i+1 << "/" 
   << gridcurve.size() << " "
   << ( (gridcurve.isClosed())?"closed":"open" ) << endl;

    // Create range of incident points
    typedef GridCurve<KSpace>::PointsRange Range;
    Range r = gridcurve.getPointsRange();//building range

    // Estimation
    cout << "# Curvature estimation from maximal segments" << endl; 
    std::vector<double> estimations1; 
    std::vector<double> estimations2; 
    if (gridcurve.isOpen())
      { 
  cout << "# open grid curve" << endl;
  estimationFromLength( h, r.begin(), r.end(), back_inserter(estimations1) ); 
  estimationFromLengthAndWidth( h, r.begin(), r.end(), back_inserter(estimations2) ); 
      }
    else
      { 
  cout << "# closed grid curve" << endl;
  estimationFromLength( h, r.c(), r.c(), back_inserter(estimations1) ); 
  estimationFromLengthAndWidth( h, r.c(), r.c(), back_inserter(estimations2) ); 
      }

    // Output
    cout << "# id curvatureFromLength curvatureFromLengthAndWidth" << endl;  
    unsigned int j = 0;
    for ( Range::ConstIterator it = r.begin(), itEnd = r.end();
    it != itEnd; ++it, ++j ) {
cout << j << setprecision( 15 )
     << " " << estimations1[ j ]
     << " " << estimations2[ j ] << endl;
    }

 }

 
  return 0;
}