C++ GridCurve类代码示例

/**
 * Exceptions
 *
 */
bool testExceptions(const string &filename)
{

  GridCurve<KhalimskySpaceND<2> > c; //grid curve

  trace.info() << endl;
  trace.info() << "Trying to read bad file: " << filename << endl;
  
  ifstream instream; // input stream
  instream.open (filename.c_str(), ifstream::in);

  try {
    c.initFromVectorStream(instream);
    trace.info() << "no exception catched!?" << endl;
    return false;
  }  catch (DGtal::ConnectivityException& e) {
    trace.info() << e.what() << endl;
    return true;
  } catch (DGtal::InputException& e) {
    trace.info() << e.what() << endl;
    return true;
  } catch (exception& e) {
    trace.info() << e.what() << endl;
    return true;
  } 
}

bool testIOGridCurve(const string& filename)
{

  unsigned int d = KSpace::Point::dimension;
  GridCurve<KSpace> c; //grid curve

//////////////////////////////////////////
  trace.info() << endl;
  trace.info() << "Reading GridCurve d=" << d << " "; 
  
  ifstream instream; // input stream
  instream.open (filename.c_str(), ifstream::in);

  c.initFromVectorStream(instream);

  trace.info() << "(" << c.size() << ") elts" << std::endl; 
  trace.info() << c << endl;

///////////////////////////////////////////
  std::stringstream s; 
  s << "gridcurve" << d << ".dat"; 

  trace.info() << "Writing GridCurve d=" << d << " in " << s.str() << endl;

  ofstream outstream(s.str().c_str()); //output stream
  if (!outstream.is_open()) return false;
  else {
    c.writeVectorToStream(outstream);
  }
  outstream.close();

  return true;
}

/**
 * Display
 *
 */
bool testDrawGridCurve(const string &filename)
{

  GridCurve<KhalimskySpaceND<2> > c; //grid curve

  trace.info() << endl;
  trace.info() << "Displaying GridCurve " << endl;
  
  //reading grid curve
  fstream inputStream;
  inputStream.open (filename.c_str(), ios::in);
  c.initFromVectorStream(inputStream); 
  inputStream.close();

  //displaying it
  Board2D aBoard;
  aBoard.setUnit(Board2D::UCentimeter);
  aBoard << c; 
  aBoard.saveEPS( "GridCurve.eps", Board2D::BoundingBox, 5000 );
#ifdef WITH_CAIRO
  aBoard.saveCairo("GridCurve-cairo.pdf", Board2D::CairoPDF, Board2D::BoundingBox, 5000);
#endif

  return true;
}

/**
 * Recogition of randomly generated digital circles
 */
bool testRecognition()
{

  typedef KhalimskySpaceND<2,int> KSpace; 
  GridCurve<KSpace> c; 
  
  unsigned int nbok = 0;
  unsigned int nb = 0;

  trace.beginBlock ( "Recognition" );
  
  for (unsigned int i = 0; i < 50; ++i)
  {
    //generate digital circle
    double cx = (rand()%100 ) / 100.0;
    double cy = (rand()%100 ) / 100.0;
    double radius = (rand()%100 )+100;
    c = ballGenerator<KSpace>( cx, cy, radius, ((i%2)==1) ); 
    trace.info() << " #ball #" << i << " c(" << cx << "," << cy << ") r=" << radius << endl; 
    
    //range
    typedef GridCurve<KSpace>::IncidentPointsRange Range; 
    Range r = c.getIncidentPointsRange();
    
    //recognition
    typedef Range::ConstIterator ConstIterator; //iterator
    StabbingCircleComputer<ConstIterator> s;
    longestSegment(s,r.begin(),r.end()); 

    //checking if the circle is separating
    bool flag = true;
    typedef CircleFrom3Points<KSpace::Point> Circle; 
    typedef functors::Point2ShapePredicate<Circle,false,true> 
      FirstInCirclePred; 
    typedef functors::Point2ShapePredicate<Circle,true,true> 
      SecondInCirclePred; 
    for (ConstIterator it = s.begin(); ((it != s.end()) && flag) ; ++it)
    {
      FirstInCirclePred p1( s.getSeparatingCircle() ); 
      SecondInCirclePred p2( s.getSeparatingCircle() ); 
      flag = ( p1(it->first)&&p2(it->second) ); 
    }
    
    //conclusion
    nbok += flag ? 1 : 0; 
    nb++;
  }

  trace.endBlock();
  
  trace.info() << "(" << nbok << "/" << nb << ") " << endl;
  return nbok == nb;
}

void
ballGenerator(const int& size, double aCx, double aCy, double aR, GridCurve<TKSpace>& gc)
{

  ASSERT( aR < (double) size ); 
 
  // Types
  typedef TKSpace KSpace;  
  typedef typename KSpace::SCell SCell;
  typedef typename KSpace::Space Space;  
  typedef typename Space::Point Point;

  KSpace K;
  bool ok = K.init( Point(-size,-size), Point(size,size), true );
  if ( ! ok )
    {
      std::cerr << " error in creating KSpace." << std::endl;
    }
  try 
    {
      BallPredicate<Point> dig(aCx, aCy, aR); 
      // 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 );
      gc.initFromVector(points); 
    }
  catch ( InputException e )
    {
      std::cerr << " error in finding a bel." << std::endl;
    }
}

/**
 * Open/Closed
 *
 */
bool testIsOpen(const string &filename, const bool& aFlag)
{

  trace.info() << endl;
  trace.info() << "Open/Closed test" << endl;

  GridCurve<KhalimskySpaceND<2> > c; //grid curve

  ifstream instream; // input stream
  instream.open (filename.c_str(), ifstream::in);
  c.initFromVectorStream(instream);

  trace.info() << c.isOpen() << " == " << aFlag << endl;

  return (c.isOpen() == aFlag);
}

int main( int argc, char** argv )
{
  trace.beginBlock ( "Testing class StabbingCircleComputer" );
  trace.info() << "Args:";
  for ( int i = 0; i < argc; ++i )
    trace.info() << " " << argv[ i ];
  trace.info() << endl;

  bool res; 
  
  {//concept checking
    testStabbingCircleComputerConceptChecking();
  }
  
  {//basic operations
    typedef KhalimskySpaceND<2,int> KSpace; 

    GridCurve<KSpace> c, rc; 
    c = ballGenerator<KSpace>(0,0,6,false); 
    std::vector<PointVector<2,int> > rv; 
    rc = ballGenerator<KSpace>(0,0,6,true); 

    res = testStabbingCircleComputer(c)
  && drawingTestStabbingCircleComputer(c, "CCW")
  && drawingTestStabbingCircleComputer(rc, "CW"); 
  }
  
  {//recognition 
    res = res && testRecognition();
  }
  
  {//segmentations
    std::string filename = testPath + "samples/sinus2D4.dat";
    ifstream instream; // input stream
    instream.open (filename.c_str(), ifstream::in);
    
    typedef KhalimskySpaceND<2,int> KSpace; 
    GridCurve<KSpace> c; //grid curve
    c.initFromVectorStream(instream);
    
    res = res && testSegmentation(c); 
  }
  
  trace.emphase() << ( res ? "Passed." : "Error." ) << endl;
  trace.endBlock();
  return res ? 0 : 1;
}

/**
 * Example of a test. To be completed.
 *
 */
bool testTrueLocalEstimator(const std::string &filename)
{
  trace.info() << "Reading GridCurve " << endl;
  ifstream instream; // input stream
  instream.open (filename.c_str(), ifstream::in);
  typedef KhalimskySpaceND<2> Kspace; //space
  GridCurve<Kspace> c; 
  c.initFromVectorStream(instream); //building grid curve
  typedef GridCurve<Kspace >::PointsRange Range;//range
  Range r = c.getPointsRange();//building range

  
  typedef Ball2D<Z2i::Space> Shape;
  typedef GridCurve<KhalimskySpaceND<2> >::PointsRange Range;
  typedef Range::ConstIterator ConstIteratorOnPoints;
  typedef ParametricShapeCurvatureFunctor< Shape > Curvature;
  typedef ParametricShapeTangentFunctor< Shape > Tangent;
  typedef ParametricShapeArcLengthFunctor< Shape > Length;

  Shape ball(Z2i::Point(0,0), 30);

   
  TrueLocalEstimatorOnPoints< ConstIteratorOnPoints, Shape, Curvature  >  curvatureEstimator;
  TrueLocalEstimatorOnPoints< ConstIteratorOnPoints, Shape, Tangent  >  tangentEstimator;
  TrueGlobalEstimatorOnPoints< ConstIteratorOnPoints, Shape, Length  >  lengthEstimator;

  curvatureEstimator.init( 1, r.begin(), r.end(), &ball, true);
  tangentEstimator.init( 1, r.begin(), r.end(), &ball, true);
 

  ConstIteratorOnPoints it = r.begin();
  //  ConstIteratorOnPoints it2 = r.begin()+15;
  ConstIteratorOnPoints it2 = it;
  for (  int compteur = 0; compteur < 15; ++compteur ) ++it2;
  lengthEstimator.init( 1, it, it2, &ball, true);
  
  
  trace.info() << "Current point = "<<*it<<std::endl;
  trace.info() << "Current point+15 = "<<*it2<<std::endl;
  trace.info() << "Eval curvature (begin, h=1) = "<< curvatureEstimator.eval(it2)<<std::endl;
  trace.info() << "Eval tangent (begin, h=1) = "<< tangentEstimator.eval(it2)<<std::endl;
  trace.info() << "Eval length ( h=1) = "<< lengthEstimator.eval(it,it2)<<std::endl;
  
  return true;

}

int main( int argc, char** argv )
{
  trace.beginBlock ( "Testing class GeometricalDSS" );
  trace.info() << "Args:";
  for ( int i = 0; i < argc; ++i )
    trace.info() << " " << argv[ i ];
  trace.info() << endl;

  bool res; 
  
  {//concept checking
    testGeometricalDSSConceptChecking();
  }
  
  {//basic operations
    std::string filename = testPath + "samples/DSS.dat";
    ifstream instream; // input stream
    instream.open (filename.c_str(), ifstream::in);
    
    typedef KhalimskySpaceND<2,int> KSpace; 
    GridCurve<KSpace> c; //grid curve
    c.initFromVectorStream(instream);

    res = testGeometricalDSS(c)
  && drawingTestGeometricalDSS(c); 
  }
  
  {//segmentations
    std::string filename = testPath + "samples/sinus2D4.dat";
    ifstream instream; // input stream
    instream.open (filename.c_str(), ifstream::in);
    
    typedef KhalimskySpaceND<2,int> KSpace; 
    GridCurve<KSpace> c; //grid curve
    c.initFromVectorStream(instream);
    
    res = res && testSegmentation(c); 
  }
  
  trace.emphase() << ( res ? "Passed." : "Error." ) << endl;
  trace.endBlock();
  return res ? 0 : 1;
}

/**
 * Test 
 *
 */
bool testFP(string filename)
{


  trace.info() << endl;
  trace.info() << "Reading GridCurve from " << filename << endl;
  
  ifstream instream; // input stream
  instream.open (filename.c_str(), ifstream::in);

  //range of points
  typedef int Coordinate; 
  typedef KhalimskySpaceND<2,Coordinate> Kspace; //space
  GridCurve<Kspace> c; //building grid curve
  c.initFromVectorStream(instream);
  typedef GridCurve<Kspace >::PointsRange Range;//range
  Range r = c.getPointsRange();//building range

  typedef Range::ConstIterator ConstIterator;//constIterator

  //faithful polyon
  trace.info() << "Building FP (process digital curve as"; 
  trace.info() << ( (c.isClosed())?"closed":"open" ) << ")" << endl;

  bool res = true; 
  if (c.isClosed())
    {
      typedef FP<Range::ConstCirculator,Coordinate,4> FP; 
      FP theFP( r.c(), r.c() );
      res = theFP.isValid();       
    }
  else 
    {
      typedef FP<Range::ConstIterator,Coordinate,4> FP; 
      FP theFP( r.begin(), r.end() );
      res = theFP.isValid(); 
    }
  return res;
 
}

/**
 * Applying test on a given data file 
 *
 */
bool testEval(string filename)
{

  trace.info() << endl;
  trace.info() << "Reading GridCurve from " << filename << endl;
  
  ifstream instream; // input stream
  instream.open (filename.c_str(), ifstream::in);
  typedef KhalimskySpaceND<2> Kspace; //space
  GridCurve<Kspace> c; //building grid curve
  c.initFromVectorStream(instream);
  typedef GridCurve<Kspace >::PointsRange Range;//range
  Range r = c.getPointsRange();//building range

  trace.info() << "Building Estimator (process range as"; 
  trace.info() << ( (c.isClosed())?"closed":"open" ) << ")" << endl;

  if (c.isClosed())
    return test(r.c(), r.c()); 
  else 
    return test(r.begin(), r.end()); 

}

int main()
{

  std::string filename = testPath + "samples/DSS.dat";
  ifstream instream; // input stream
  instream.open (filename.c_str(), ifstream::in);
  
  typedef KhalimskySpaceND<2,int> KSpace; 
  GridCurve<KSpace> c; //grid curve
  c.initFromVectorStream(instream);


  trace.beginBlock("Simple preimage test");

  typedef StraightLineFrom2Points<GridCurve<KSpace>::Point> StraightLine;
  StraightLine aStraightLine; //instance of straight line
  typedef Preimage2D<StraightLine> Preimage2D;

  GridCurve<KSpace>::IncidentPointsRange r = c.getIncidentPointsRange(); //range
  GridCurve<KSpace>::IncidentPointsRange::ConstReverseIterator it (r.rbegin()); //iterators
  GridCurve<KSpace>::IncidentPointsRange::ConstReverseIterator itEnd (r.rend()); 

  //preimage computation
  Preimage2D thePreimage(it->first, it->second, aStraightLine);
  ++it; 
  while ( (it != itEnd) &&
              (thePreimage.addBack(it->first, it->second)) )
  {
  trace.info() << (it - r.rbegin()) << endl << thePreimage << endl;
    ++it;
  }

  trace.endBlock();

  return 0;
}

int main( int argc, char** argv )
{
  trace.beginBlock ( "Testing class GridCurve" );
  trace.info() << "Args:";
  for ( int i = 0; i < argc; ++i )
    trace.info() << " " << argv[ i ];
  trace.info() << endl;


  std::string sinus2D4 = testPath + "samples/sinus2D4.dat";
  std::string polyg2D = testPath + "samples/polyg2D.dat";
  std::string sinus3D = testPath + "samples/sinus3D.dat";
  std::string emptyFile = testPath + "samples/emptyFile.dat";
  std::string square = testPath + "samples/smallSquare.dat";

  typedef KhalimskySpaceND<2> K2;
  typedef KhalimskySpaceND<3> K3;

///////// general tests
  bool res = testIOGridCurve<K2>(sinus2D4)
    && testIOGridCurve<K3>(sinus3D)
    && testExceptions(sinus3D)
    && testExceptions(polyg2D)
    && testExceptions(emptyFile)
    && testDrawGridCurve(sinus2D4)
    && testIsOpen(sinus2D4,true)
    && testIsOpen(square,false); 

/////////// ranges test
  typedef GridCurve<K2> GridCurve;

testRangeConceptChecking<GridCurve::SCellsRange>();
testRangeConceptChecking<GridCurve::SCellsRange>();
testRangeConceptChecking<GridCurve::PointsRange>();
testRangeConceptChecking<GridCurve::MidPointsRange>();
testRangeConceptChecking<GridCurve::ArrowsRange>();
testRangeConceptChecking<GridCurve::InnerPointsRange>();
testRangeConceptChecking<GridCurve::OuterPointsRange>();
testRangeConceptChecking<GridCurve::IncidentPointsRange>();

  //reading grid curve
  GridCurve c; 
  fstream inputStream;
  inputStream.open (square.c_str(), ios::in);
  c.initFromVectorStream(inputStream);
  inputStream.close();

  res = res 
    && testRange<GridCurve::SCellsRange>(c.getSCellsRange())
    && testRange<GridCurve::PointsRange>(c.getPointsRange())
    && testRange<GridCurve::MidPointsRange>(c.getMidPointsRange())
    && testPairsRange<GridCurve::ArrowsRange>(c.getArrowsRange())
    && testRange<GridCurve::InnerPointsRange>(c.getInnerPointsRange())
    && testRange<GridCurve::OuterPointsRange>(c.getOuterPointsRange())
    && testPairsRange<GridCurve::IncidentPointsRange>(c.getIncidentPointsRange())
    && testRange<GridCurve::CodesRange>(c.getCodesRange())
;

  res = res 
    && testDisplayRange<GridCurve::SCellsRange>(c.getSCellsRange())
    && testDisplayRange<GridCurve::PointsRange>(c.getPointsRange())
    && testDisplayRange<GridCurve::MidPointsRange>(c.getMidPointsRange())
    && testDisplayRange<GridCurve::ArrowsRange>(c.getArrowsRange())
    && testDisplayRange<GridCurve::InnerPointsRange>(c.getInnerPointsRange())
    && testDisplayRange<GridCurve::OuterPointsRange>(c.getOuterPointsRange())
    && testDisplayRange<GridCurve::IncidentPointsRange>(c.getIncidentPointsRange())
    && testDisplayRange<GridCurve::CodesRange>(c.getCodesRange())
;

  res = res 
    && testDrawRange<GridCurve::SCellsRange>(c.getSCellsRange(),"1cells","Grid")
    && testDrawRange<GridCurve::PointsRange>(c.getPointsRange(),"Points","Paving")
    && testDrawRange<GridCurve::MidPointsRange>(c.getMidPointsRange(),"MidPoints","Paving")
    && testDrawRange<GridCurve::ArrowsRange>(c.getArrowsRange(),"Arrows","Paving")
    && testDrawRange<GridCurve::InnerPointsRange>(c.getInnerPointsRange(),"InnerPoints","Grid")
    && testDrawRange<GridCurve::OuterPointsRange>(c.getOuterPointsRange(),"OuterPoints","Grid")
    && testDrawRange<GridCurve::IncidentPointsRange>(c.getIncidentPointsRange(),"IncidentPoints","Grid")
;

//////////////////////

  trace.emphase() << ( res ? "Passed." : "Error." ) << endl;
  trace.endBlock();
  
  return res ? 0 : 1;
}

bool
lengthEstimators( const std::string & /*name*/,
      Shape & aShape, 
      double h )
{
  // Types
  typedef typename Space::Point Point;
  typedef typename Space::Vector Vector;
  typedef typename Space::RealPoint RealPoint;
  typedef typename Space::Integer Integer;
  typedef HyperRectDomain<Space> Domain;
  typedef KhalimskySpaceND<Space::dimension,Integer> KSpace;
  typedef typename KSpace::SCell SCell;
  typedef typename GridCurve<KSpace>::PointsRange PointsRange;
  typedef typename GridCurve<KSpace>::ArrowsRange ArrowsRange;

  // Digitizer
  GaussDigitizer<Space,Shape> dig;  
  dig.attach( aShape ); // attaches the shape.
  Vector vlow(-1,-1); Vector vup(1,1);
  dig.init( aShape.getLowerBound()+vlow, aShape.getUpperBound()+vup, h ); 
  Domain domain = dig.getDomain();

  // Create cellular space
  KSpace K;
  bool ok = K.init( dig.getLowerBound(), dig.getUpperBound(), true );
  if ( ! ok )
    {
      std::cerr << "[lengthEstimators]"
    << " 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 );
    // Create GridCurve
    GridCurve<KSpace> gridcurve;
    gridcurve.initFromVector( points );
    // Ranges
    ArrowsRange ra = gridcurve.getArrowsRange(); 
    PointsRange rp = gridcurve.getPointsRange(); 


    // Estimations
    typedef typename PointsRange::ConstIterator ConstIteratorOnPoints; 
    typedef ParametricShapeArcLengthFunctor< Shape > Length;
    TrueGlobalEstimatorOnPoints< ConstIteratorOnPoints, Shape, Length  >  trueLengthEstimator;
    trueLengthEstimator.init( h, rp.begin(), rp.end(), &aShape, gridcurve.isClosed());

    L1LengthEstimator< typename ArrowsRange::ConstCirculator > l1length;
    DSSLengthEstimator< typename PointsRange::ConstCirculator > DSSlength;
    MLPLengthEstimator< typename PointsRange::ConstIterator > MLPlength;
    FPLengthEstimator< typename PointsRange::ConstIterator > FPlength;
    BLUELocalLengthEstimator< typename ArrowsRange::ConstIterator > BLUElength;
    RosenProffittLocalLengthEstimator< typename ArrowsRange::ConstIterator > RosenProffittlength;
  
    // Output
    double trueValue = trueLengthEstimator.eval();
    double l1, blue, rosen,dss,mlp,fp;
    double Tl1, Tblue, Trosen,Tdss,Tmlp,Tfp;
    
    Clock c;

    //Length evaluation & timing
    c.startClock();
    l1length.init(h, ra.c(), ra.c());
    l1 = l1length.eval();
    Tl1 = c.stopClock();
    
    c.startClock();
    BLUElength.init(h, ra.begin(), ra.end(), gridcurve.isClosed());
    blue = BLUElength.eval();
    Tblue = c.stopClock();
    
    c.startClock();
    RosenProffittlength.init(h, ra.begin(), ra.end(), gridcurve.isClosed());
    rosen = RosenProffittlength.eval();
    Trosen = c.stopClock();
    
    c.startClock();
    DSSlength.init(h, rp.c(), rp.c());
    dss = DSSlength.eval();
    Tdss = c.stopClock();
    
    c.startClock();
    MLPlength.init(h, rp.begin(), rp.end(), gridcurve.isClosed());
    mlp = MLPlength.eval();
    Tmlp = c.stopClock();

    c.startClock();
    FPlength.init(h, rp.begin(), rp.end(), gridcurve.isClosed());
    fp = FPlength.eval();
    Tfp = c.stopClock();

    std::cout << std::setprecision( 15 ) << h << " " << rp.size() << " " << trueValue 
   << " " << l1
//.........这里部分代码省略.........

bool testCompareEstimator(const std::string &name, Shape & aShape, double h)
{
  using namespace Z2i;

  trace.beginBlock ( ( "Testing CompareEstimator on digitization of "
           + name ). c_str() );
  
  // Creates a digitizer on the window (xLow, xUp).
  typedef Space::RealPoint RealPoint;
  RealPoint xLow( -10.0, -10.0 );
  RealPoint xUp( 10.0, 10.0 );
  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 << "[testCompareEstimators]"
    << " error in creating KSpace." << std::endl;
    }
  else
    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 );
      // Create GridCurve
      GridCurve<KSpace> gridcurve;
      gridcurve.initFromVector( points );
      typedef GridCurve<KhalimskySpaceND<2> >::PointsRange Range;
      typedef Range::ConstIterator ConstIteratorOnPoints;
      Range r = gridcurve.getPointsRange();//building range

      unsigned int nb = 0; 
      unsigned int nbok = 0; 
      //curvature
      typedef ParametricShapeCurvatureFunctor< Shape > Curvature;
      typedef TrueLocalEstimatorOnPoints< ConstIteratorOnPoints, Shape, Curvature  >  TrueCurvature;
      TrueCurvature curvatureEstimator;
      TrueCurvature curvatureEstimatorBis;
      curvatureEstimator.init( h, r.begin(), r.end() );
      curvatureEstimator.attach( &aShape ); 
      curvatureEstimatorBis.init( h, r.begin(), r.end() );
      curvatureEstimatorBis.attach( &aShape ); 

      typedef CompareLocalEstimators< TrueCurvature, TrueCurvature> Comparator;

      trace.info()<< "True curvature comparison at "<< *r.begin() << " = "
		  << Comparator::compare(curvatureEstimator,curvatureEstimatorBis, r.begin())
		  << std::endl;
      
      typename Comparator::OutputStatistic error
	=Comparator::compare(curvatureEstimator, curvatureEstimatorBis, 
			     r.begin(),
			     r.end());
      
      trace.info() << "Nb samples= "<< error.samples()<<std::endl;
      trace.info() << "Error mean= "<< error.mean()<<std::endl;
      trace.info() << "Error max= "<< error.max()<<std::endl;
      nbok += ( (error.samples() == r.size())&&(error.max() == 0) )?1:0; 
      nb++;
      trace.info() << nbok << "/" << nb << std::endl; 

      //tangents
      typedef ParametricShapeTangentFunctor< Shape > Tangent;
      typedef TrueLocalEstimatorOnPoints< ConstIteratorOnPoints, Shape, Tangent  >  TrueTangent;

      typedef ArithmeticalDSS<ConstIteratorOnPoints,KSpace::Integer,4> 
  SegmentComputer;
      typedef TangentFromDSSEstimator<SegmentComputer> Functor;
      typedef MostCenteredMaximalSegmentEstimator<SegmentComputer,Functor> 
  MSTangentEstimator;

      SegmentComputer sc;
      Functor f; 
      
      TrueTangent tang1;
      MSTangentEstimator tang2(sc, f); 
    
      tang1.init( h, r.begin(), r.end() );
      tang1.attach( &aShape ); 
      tang2.init( h, r.begin(), r.end() );
      
      typedef CompareLocalEstimators< TrueTangent, MSTangentEstimator> ComparatorTan;

      trace.info()<< "Tangent comparison at "<< *r.begin() << " = " 
		  << ComparatorTan::compareVectors( tang1, tang2, r.begin())
		  << std::endl; 
      
      typename ComparatorTan::OutputVectorStatistic error2
	=ComparatorTan::compareVectors(tang1, tang2, 
//.........这里部分代码省略.........