tangentBC.cpp

#include <cmath>
#include <iostream>
#include <vector>
#include <string>

#include "CLI11.hpp"

#include "DGtal/base/Common.h"
#include "DGtal/helpers/StdDefs.h"
#include "DGtal/io/readers/PointListReader.h"

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

//Estimators
#include "DGtal/geometry/curves/BinomialConvolver.h"

using namespace DGtal;


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

  app.description("Estimates tangent using a binomial convolver.\n Typical use example:\n \t tangentBC [options] --input  <fileName>\n");
  auto filenameOpt = app.add_option("--input,-i,1",fileName,"input file name: FreemanChain (.fc) or a sequence of discrete points (.sdp).")->required()->check(CLI::ExistingFile);
  app.add_option("--GridStep", h, "Grid step (default 1.0)", true);
   
  app.get_formatter()->column_width(40);
  CLI11_PARSE(app, argc, argv);
  // END parse command line using CLI ----------------------------------------------  

  if(filenameOpt->count()>0){
    std::string extension =  fileName.substr( fileName.find_last_of(".") + 1 );
    bool isSDP = extension == "sdp";
    typedef Z2i::Space Space; 
    typedef Space::Point Point; 
    typedef PointVector<2, double> RealPoint; 
    typedef Space::Integer Integer;  
    typedef FreemanChain<Integer> FreemanChain; 
    typedef std::vector< Point > Storage;
    typedef Storage::const_iterator ConstIteratorOnPoints; 

    std::vector< FreemanChain > vectFcs;
    if(!isSDP)
      {
        vectFcs =   
          PointListReader< Point >:: getFreemanChainsFromFile<Integer> (fileName); 
      }
    for(unsigned int i=0; i<vectFcs.size() || (i==0 && isSDP); i++){
      Storage vectPts; 
      bool isClosed;
      if(!isSDP)
        {
          isClosed = vectFcs.at(i).isClosed(); 
          std::cout << "# grid curve " << i << "/" << vectFcs.size() << " "
                    << ( (isClosed)?"closed":"open" ) << std::endl;
          FreemanChain::getContourPoints( vectFcs.at(i), vectPts ); 
        }
      else
        {
          vectPts = PointListReader<Z2i::Point>::getPointsFromFile(fileName);
          Z2i::Point pf =vectPts[0];
          Z2i::Point pl =vectPts[vectPts.size()-1];
          isClosed = (pf[0]-pl[0])+(pf[1]-pl[1]) <= 1;
        }
      
      // Binomial
      std::cout << "# Curvature estimation from binomial convolution" << std::endl;
      typedef BinomialConvolver<ConstIteratorOnPoints, double> MyBinomialConvolver;
      std::cout << "# mask size = " << 
      MyBinomialConvolver::suggestedSize( h, vectPts.begin(), vectPts.end() ) << std::endl;
      typedef 
        TangentFromBinomialConvolverFunctor< MyBinomialConvolver, RealPoint >
        TangentBCFct;
      BinomialConvolverEstimator< MyBinomialConvolver, TangentBCFct> 
        BCTangentEstimator;
      
      BCTangentEstimator.init( h, vectPts.begin(), vectPts.end(), isClosed );

      std::vector<RealPoint> tangents( vectPts.size() ); 
      BCTangentEstimator.eval( vectPts.begin(), vectPts.end(), 
             tangents.begin() ); 

      // Output
      std::cout << "# id tangent.x tangent.y angle(atan2(y,x)) x y" << std::endl;  
      unsigned int j = 0;
      for ( ConstIteratorOnPoints 
        it = vectPts.begin(), it_end = vectPts.end();
      it != it_end; ++it, ++j ) 
  {
    double x = tangents[ j ][ 0 ];
    double y = tangents[ j ][ 1 ];
    std::cout << j << std::setprecision( 15 )
              << " " << x << " " << y 
              << " " << atan2( y, x ) << " " << (*it)[0] << " " << (*it)[1]
              << std::endl;
  }

    }

  }
  return 0;
}