testLengthEstimators.cpp

#include <cmath>
#include "DGtal/base/Common.h"
 
//space / domain
#include "DGtal/kernel/SpaceND.h"
#include "DGtal/kernel/domains/HyperRectDomain.h"
#include "DGtal/topology/KhalimskySpaceND.h"
 
//shape and digitizer
#include "DGtal/shapes/ShapeFactory.h"
#include "DGtal/shapes/Shapes.h"
#include "DGtal/helpers/StdDefs.h"
#include "DGtal/topology/helpers/Surfaces.h"
 
#include "DGtal/shapes/GaussDigitizer.h"
#include "DGtal/geometry/curves/GridCurve.h"
 
//estimators
#include "DGtal/geometry/curves/estimation/TrueLocalEstimatorOnPoints.h"
#include "DGtal/geometry/curves/estimation/TrueGlobalEstimatorOnPoints.h"
#include "DGtal/geometry/curves/estimation/ParametricShapeCurvatureFunctor.h"
#include "DGtal/geometry/curves/estimation/ParametricShapeTangentFunctor.h"
#include "DGtal/geometry/curves/estimation/ParametricShapeArcLengthFunctor.h"
 
#include "DGtal/geometry/curves/estimation/L1LengthEstimator.h"
#include "DGtal/geometry/curves/estimation/TwoStepLocalLengthEstimator.h"
#include "DGtal/geometry/curves/estimation/BLUELocalLengthEstimator.h"
#include "DGtal/geometry/curves/estimation/RosenProffittLocalLengthEstimator.h"
 
#include "DGtal/geometry/curves/estimation/MLPLengthEstimator.h"
#include "DGtal/geometry/curves/estimation/FPLengthEstimator.h"
#include "DGtal/geometry/curves/estimation/DSSLengthEstimator.h"
 
#include "ConfigTest.h"
 
#include "DGtal/io/boards/Board2D.h"

 
using namespace std;
using namespace DGtal;
 

// Functions for testing Length Estimator classes.
 
bool testLengthEstimatorsOnBall(double radius, double h)
{
 
  // Types
  typedef SpaceND<2,int> Space;  
  typedef Ball2D<Space> Shape;
  typedef Space::Point Point;
  typedef Space::RealPoint RealPoint;
  typedef Space::Integer Integer;
  typedef HyperRectDomain<Space> Domain;
  typedef KhalimskySpaceND<Space::dimension,Integer> KSpace;
  typedef KSpace::SCell SCell;
  typedef GridCurve<KSpace>::PointsRange PointsRange;
  typedef GridCurve<KSpace>::ArrowsRange ArrowsRange;
 
 
  //Forme
  Shape aShape(Point(0,0), radius);
 
  trace.info() << "#ball created, r=" << radius << endl;
 
  // Window for the estimation
  RealPoint xLow ( -radius-1, -radius-1 );
  RealPoint xUp( radius+1, radius+1 );
  GaussDigitizer<Space,Shape> dig;  
  dig.attach( aShape ); // attaches the shape.
  dig.init( xLow, xUp, h ); 
  // The domain size is given by the digitizer according to the window
  // and the step.
  Domain domain = dig.getDomain();
  // Create cellular space
  KSpace K;
  bool ok = K.init( dig.getLowerBound(), dig.getUpperBound(), true );
  if ( ! ok )
    {
      std::cerr << " "
                << " error in creating KSpace." << std::endl;
      return false;
    }
 
  SurfelAdjacency<KSpace::dimension> SAdj( true );
  SCell bel; 
  try 
    {
      bel = Surfaces<KSpace>::findABel( K, dig, 10000 );
    }    
  catch ( InputException& e )
    {
      std::cerr << " "
                << " error in finding a bel." << std::endl;
      return false;
    }
 
    // Getting the consecutive surfels of the 2D boundary
    std::vector<Point> points;
    Surfaces<KSpace>::track2DBoundaryPoints( points, K, SAdj, dig, bel );
    trace.info() << "#tracking..." << endl;
    // Create GridCurve
    GridCurve<KSpace> gridcurve;
    gridcurve.initFromPointsVector( points );
    trace.info() << "#grid curve created, h=" << h << endl;
 
    //ranges
    ArrowsRange ra = gridcurve.getArrowsRange(); 
    PointsRange rp = gridcurve.getPointsRange(); 

    double trueValue = M_PI*2*radius;
    L1LengthEstimator< GridCurve<KSpace>::PointsRange::ConstCirculator > l1length;
    TwoStepLocalLengthEstimator< GridCurve<KSpace>::ArrowsRange::ConstIterator > locallength(1.0,sqrt(2.0));
    BLUELocalLengthEstimator< GridCurve<KSpace>::ArrowsRange::ConstIterator > BLUElength;
    RosenProffittLocalLengthEstimator< GridCurve<KSpace>::ArrowsRange::ConstIterator > RosenProffittlength;
    DSSLengthEstimator< GridCurve<KSpace>::PointsRange::ConstCirculator > DSSlength;
    MLPLengthEstimator< GridCurve<KSpace>::PointsRange::ConstCirculator > MLPlength;
    FPLengthEstimator< GridCurve<KSpace>::PointsRange::ConstCirculator > FPlength;
 
    trace.info() << "#Estimations" <<std::endl;
    trace.info() << "#h true naive 1-sqrt(2) BLUE RosenProffitt DSS MLP FP " <<std::endl;
    trace.info() << h << " " << trueValue  
                 << " " << l1length.eval(rp.c(), rp.c(), h) 
                 << " " << locallength.eval(ra.begin(), ra.end(), h) 
                 << " " << BLUElength.eval(ra.begin(), ra.end(), h) 
                 << " " << RosenProffittlength.eval(ra.begin(), ra.end(), h)
                 << " " << DSSlength.eval(rp.c(), rp.c(), h) 
                 << " " << MLPlength.eval(rp.c(), rp.c(), h)
                 << " " << FPlength.eval(rp.c(), rp.c(), h)
     << std::endl;
 
  return true;
}
 
 
 
bool testDisplay(double radius, double h)
{
 
  // Types
  typedef Ball2D<Z2i::Space> Shape;
  typedef Z2i::Space::Point Point;
  typedef Z2i::Space::RealPoint RealPoint;
  typedef Z2i::Space::Integer Integer;
  typedef HyperRectDomain<Z2i::Space> Domain;
  typedef KhalimskySpaceND<Z2i::Space::dimension,Integer> KSpace;
  typedef KSpace::SCell SCell;
  typedef GridCurve<KSpace>::PointsRange PointsRange;
  typedef GridCurve<KSpace>::ArrowsRange ArrowsRange;
  typedef GridCurve<KSpace>::SCellsRange SCellsRange;
 
 
  //Forme
  Shape aShape(Point(0,0), radius);
 
  trace.info() << "#ball created, r=" << radius << endl;
 
  // Window for the estimation
  RealPoint xLow ( -radius-1, -radius-1 );
  RealPoint xUp( radius+1, radius+1 );
  GaussDigitizer<Z2i::Space,Shape> dig;  
  dig.attach( aShape ); // attaches the shape.
  dig.init( xLow, xUp, h ); 
  // The domain size is given by the digitizer according to the window
  // and the step.
  Domain domain = dig.getDomain();
  // Create cellular space
  KSpace K;
  bool ok = K.init( dig.getLowerBound(), dig.getUpperBound(), true );
  if ( ! ok )
    {
      std::cerr << " "
                << " error in creating KSpace." << std::endl;
      return false;
    }
  try {
      
    // Extracts shape boundary
    SurfelAdjacency<KSpace::dimension> SAdj( true );
    SCell bel = Surfaces<KSpace>::findABel( K, dig, 10000 );
    // Getting the consecutive surfels of the 2D boundary
    std::vector<Point> points;
    Surfaces<KSpace>::track2DBoundaryPoints( points, K, SAdj, dig, bel );
    trace.info() << "# tracking..." << endl;
    // Create GridCurve
    GridCurve<KSpace> gridcurve;
    gridcurve.initFromPointsVector( points );
    trace.info() << "#grid curve created, h=" << h << endl;
    
    //ranges
    ArrowsRange ra = gridcurve.getArrowsRange(); 
    PointsRange rp = gridcurve.getPointsRange(); 
    SCellsRange rc = gridcurve.getSCellsRange(); 
    
    //Explicit reshaping for drawing purposes
    Z2i::DigitalSet set(domain);
    Shapes<Z2i::Domain>::euclideanShaper( set, Shape(Point(0,0), radius/h));
         
    Board2D board;
    
    board << domain << set;
    board.saveSVG( "Ranges-Set.svg" );
    
    board.clear();
    board << domain;
    board << rp; 
    // for( PointsRange::ConstIterator it =rp.begin(), ite=rp.end();
    //      it != ite; ++it)
    //   board << (*it);
    board.saveSVG( "Ranges-Points.svg" );
    
 
    board.clear();
    board << domain;
    board << rc; 
    // for( SCellsRange::ConstIterator it =rc.begin(), ite=rc.end();
    //      it != ite; ++it)
    //   board << (*it);
    board.saveSVG( "Ranges-SCells.svg" );
    
 
    board.clear();
    board << domain;
    board << ra; 
    // Z2i::Space::Vector shift;
    // board.setPenColor( Color::Black );
    // for(  ArrowsRange::ConstIterator it = ra.begin(), itend = ra.end();
    //       it != itend;   
    //       ++it)
    //   {
    //     shift =   (*it).second ;
    //     draw(board, shift, (*it).first );
    //   }
    board.saveSVG( "Ranges-Arrows.svg" );
 
  }    
  catch ( InputException& e )
    {
      std::cerr << " "
                << " error in finding a bel." << std::endl;
      return false;
    }
 
 
 
  return true;
}

// Standard services - public :
 
int main( int argc, char** argv )
{
  trace.beginBlock ( "Testing class LengthEstimators" );
  trace.info() << "Args:";
  for ( int i = 0; i < argc; ++i )
    trace.info() << " " << argv[ i ];
  trace.info() << endl;
 
  double r = 5; 
  bool res = testLengthEstimatorsOnBall(r,1)
    && testLengthEstimatorsOnBall(r,0.1)
    && testLengthEstimatorsOnBall(r,0.01)
    && testLengthEstimatorsOnBall(r,0.001)
    && testDisplay(r,0.9);
  ;
 
  trace.emphase() << ( res ? "Passed." : "Error." ) << endl;
  trace.endBlock();
  return res ? 0 : 1;
}
//                                                                           //