curvatureBC.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/BinomialConvolver.h"


#include <boost/program_options/options_description.hpp>
#include <boost/program_options/parsers.hpp>
#include <boost/program_options/variables_map.hpp>

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

using namespace DGtal;



namespace po = boost::program_options;

int main( int argc, char** argv )
{
  // parse command line ----------------------------------------------
  po::options_description general_opt("Allowed options are: ");
  general_opt.add_options()
    ("help,h", "display this message")
    ("input,i", po::value<std::string>(), "input FreemanChain file name")
    ("GridStep,step", po::value<double>()->default_value(1.0), "Grid step");
  
  
  bool parseOK=true;
  po::variables_map vm;
  try{
    po::store(po::parse_command_line(argc, argv, general_opt), vm);  
  }catch(const std::exception& ex){
    parseOK=false;
    trace.info()<< "Error checking program options: "<< ex.what()<< std::endl;
  }
  po::notify(vm);    
  if(!parseOK || vm.count("help")||argc<=1 || (!(vm.count("input"))) )
    {
      trace.info()<< "Curvature using a binomial convolver " <<std::endl << "Basic usage: "<<std::endl
      << "\t curvatureBC [options] --input  <fileName> "<<std::endl
      << general_opt << "\n";
      return 0;
    }
  
  
  double h = vm["GridStep"].as<double>();  


 
  if(vm.count("input")){
    std::string fileName = vm["input"].as<std::string>();

    typedef Z2i::Space Space; 
    typedef Space::Point Point; 
    typedef Space::Integer Integer;  
    typedef FreemanChain<Integer> FreemanChain; 
    typedef std::vector< Point > Storage;
    typedef Storage::const_iterator ConstIteratorOnPoints; 

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

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

      bool isClosed = vectFcs.at(i).isClosed(); 
      std::cout << "# grid curve " << i+1 << "/" << vectFcs.size() << " "
      << ( (isClosed)?"closed":"open" ) << std::endl;

      Storage vectPts; 
      FreemanChain::getContourPoints( vectFcs.at(i), vectPts ); 

      // 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 CurvatureFromBinomialConvolverFunctor< MyBinomialConvolver, double >
      CurvatureBCFct;
      BinomialConvolverEstimator< MyBinomialConvolver, CurvatureBCFct> BCCurvatureEstimator;

      BCCurvatureEstimator.init( h, vectPts.begin(), vectPts.end(), isClosed );

      std::vector <double> curvatures( vectPts.size() ); 
      BCCurvatureEstimator.eval( vectPts.begin(), vectPts.end(), curvatures.begin() ); 

      // Output
      std::cout << "# id curvature" << std::endl;  
      unsigned int j = 0;
      for ( ConstIteratorOnPoints it = vectPts.begin(), it_end = vectPts.end();
      it != it_end; ++it, ++j ) {
  std::cout << j << std::setprecision( 15 )
       << " " << curvatures[ j ] << std::endl;
      }

   }

 }
 
  return 0;
}