DGtalTools  0.9.2
3dLocalEstimators.cpp
1 
30 #include <iostream>
32 #include <iomanip>
33 #include <string>
34 
35 #include "DGtal/base/Common.h"
36 #include "DGtal/base/Clock.h"
37 #include "DGtal/helpers/StdDefs.h"
38 
39 #include <boost/program_options/options_description.hpp>
40 #include <boost/program_options/parsers.hpp>
41 #include <boost/program_options/variables_map.hpp>
42 
43 //shapes
44 #include "DGtal/shapes/implicit/ImplicitBall.h"
45 #include "DGtal/math/MPolynomial.h"
46 #include "DGtal/io/readers/MPolynomialReader.h"
47 #include "DGtal/shapes/implicit/ImplicitPolynomial3Shape.h"
48 
49 //Digitizer
50 #include "DGtal/shapes/GaussDigitizer.h"
51 #include "DGtal/topology/LightImplicitDigitalSurface.h"
52 #include "DGtal/topology/DigitalSurface.h"
53 #include "DGtal/graph/DepthFirstVisitor.h"
54 #include "DGtal/geometry/volumes/KanungoNoise.h"
55 #include "DGtal/topology/CanonicSCellEmbedder.h"
56 
57 
58 
59 //Estimators
60 #include "DGtal/kernel/BasicPointFunctors.h"
61 #include "DGtal/geometry/surfaces/FunctorOnCells.h"
62 
63 #include "DGtal/geometry/curves/estimation/TrueLocalEstimatorOnPoints.h"
64 #include "DGtal/geometry/surfaces/estimation/IIGeometricFunctors.h"
65 #include "DGtal/geometry/surfaces/estimation/IntegralInvariantVolumeEstimator.h"
66 #include "DGtal/geometry/surfaces/estimation/IntegralInvariantCovarianceEstimator.h"
67 
68 #include "DGtal/geometry/surfaces/estimation/LocalEstimatorFromSurfelFunctorAdapter.h"
69 #include "DGtal/geometry/surfaces/estimation/estimationFunctors/MongeJetFittingMeanCurvatureEstimator.h"
70 #include "DGtal/geometry/surfaces/estimation/estimationFunctors/MongeJetFittingGaussianCurvatureEstimator.h"
71 #include "DGtal/geometry/surfaces/estimation/estimationFunctors/MongeJetFittingPrincipalCurvaturesEstimator.h"
72 
73 using namespace DGtal;
74 using namespace functors;
75 
76 typedef std::pair<double,double> PrincipalCurvatures;
77 
78 template < typename Shape, typename KSpace, typename ConstIterator, typename OutputIterator >
79 void
80 estimateTrueMeanCurvatureQuantity( const ConstIterator & it_begin,
81  const ConstIterator & it_end,
82  OutputIterator & output,
83  const KSpace & K,
84  const double & h,
85  Shape * aShape )
86 {
87  typedef typename KSpace::Space::RealPoint RealPoint;
88  typedef CanonicSCellEmbedder< KSpace > Embedder;
89 
90  Embedder embedder( K );
91  RealPoint currentRealPoint;
92 
93  for ( ConstIterator it = it_begin; it != it_end; ++it )
94  {
95  currentRealPoint = embedder( *it_begin ) * h;
96  *output = aShape->meanCurvature( currentRealPoint );
97  ++output;
98  }
99 }
100 
101 template < typename Shape, typename KSpace, typename ConstIterator, typename OutputIterator >
102 void
103 estimateTrueGaussianCurvatureQuantity( const ConstIterator & it_begin,
104  const ConstIterator & it_end,
105  OutputIterator & output,
106  const KSpace & K,
107  const double & h,
108  Shape * aShape )
109 {
110  typedef typename KSpace::Space::RealPoint RealPoint;
111  typedef CanonicSCellEmbedder< KSpace > Embedder;
112 
113  Embedder embedder( K );
114  RealPoint currentRealPoint;
115 
116  for ( ConstIterator it = it_begin; it != it_end; ++it )
117  {
118  currentRealPoint = embedder( *it_begin ) * h;
119  *output = aShape->gaussianCurvature( currentRealPoint );
120  ++output;
121  }
122 }
123 
124 template < typename Shape, typename KSpace, typename ConstIterator, typename OutputIterator >
125 void
126 estimateTruePrincipalCurvaturesQuantity( const ConstIterator & it_begin,
127  const ConstIterator & it_end,
128  OutputIterator & output,
129  const KSpace & K,
130  const double & h,
131  Shape * aShape )
132 {
133  typedef typename KSpace::Space::RealPoint RealPoint;
134  typedef CanonicSCellEmbedder< KSpace > Embedder;
135 
136  Embedder embedder( K );
137  RealPoint currentRealPoint;
138 
139  for ( ConstIterator it = it_begin; it != it_end; ++it )
140  {
141  currentRealPoint = embedder( *it_begin ) * h;
142  double k1, k2;
143  aShape->principalCurvatures( currentRealPoint, k1, k2 );
144  PrincipalCurvatures result;
145  result.first = k1;
146  result.second = k2;
147  *output = result;
148  ++output;
149  }
150 }
151 
152 template <typename Space, typename Shape>
153 bool
154 compareShapeEstimators( const std::string & filename,
155  const Shape * aShape,
156  const typename Space::RealPoint & border_min,
157  const typename Space::RealPoint & border_max,
158  const double & h,
159  const double & radius_kernel,
160  const double & alpha,
161  const std::string & options,
162  const std::string & properties,
163  const bool & lambda_optimized,
164  double noiseLevel = 0.0 )
165 {
166  typedef typename Space::RealPoint RealPoint;
167  typedef GaussDigitizer< Z3i::Space, Shape > DigitalShape;
168  typedef Z3i::KSpace KSpace;
169  typedef typename KSpace::SCell SCell;
170  typedef typename KSpace::Surfel Surfel;
171 
172  bool withNoise = ( noiseLevel <= 0.0 ) ? false : true;
173 
174  ASSERT (( noiseLevel < 1.0 ));
175  // Digitizer
176  DigitalShape* dshape = new DigitalShape();
177  dshape->attach( *aShape );
178  dshape->init( border_min, border_max, h );
179 
180  KSpace K;
181  if ( ! K.init( dshape->getLowerBound(), dshape->getUpperBound(), true ) )
182  {
183  std::cerr << "[3dLocalEstimators] error in creating KSpace." << std::endl;
184  return false;
185  }
186 
187  try
188  {
189  if ( withNoise )
190  {
191  typedef KanungoNoise< DigitalShape, Z3i::Domain > KanungoPredicate;
192  typedef LightImplicitDigitalSurface< KSpace, KanungoPredicate > Boundary;
193  typedef DigitalSurface< Boundary > MyDigitalSurface;
194  typedef typename MyDigitalSurface::ConstIterator ConstIterator;
195 
196  typedef DepthFirstVisitor< MyDigitalSurface > Visitor;
197  typedef GraphVisitorRange< Visitor > VisitorRange;
198  typedef typename VisitorRange::ConstIterator VisitorConstIterator;
199 
200  // typedef PointFunctorFromPointPredicateAndDomain< KanungoPredicate, Z3i::Domain, unsigned int > MyPointFunctor;
201  // typedef FunctorOnCells< MyPointFunctor, KSpace > MySpelFunctor;
202 
203  // Extracts shape boundary
204  KanungoPredicate * noisifiedObject = new KanungoPredicate( *dshape, dshape->getDomain(), noiseLevel );
205  SCell bel = Surfaces< KSpace >::findABel( K, *noisifiedObject, 10000 );
206  Boundary * boundary = new Boundary( K, *noisifiedObject, SurfelAdjacency< KSpace::dimension >( true ), bel );
207  MyDigitalSurface surf ( *boundary );
208 
209  double minsize = dshape->getUpperBound()[0] - dshape->getLowerBound()[0];
210  unsigned int tries = 0;
211  while( surf.size() < 2 * minsize || tries > 150 )
212  {
213  delete boundary;
214  bel = Surfaces< KSpace >::findABel( K, *noisifiedObject, 10000 );
215  boundary = new Boundary( K, *noisifiedObject, SurfelAdjacency< KSpace::dimension >( true ), bel );
216  surf = MyDigitalSurface( *boundary );
217  ++tries;
218  }
219 
220  if( tries > 150 )
221  {
222  std::cerr << "Can't found a proper bel. So .... I ... just ... kill myself." << std::endl;
223  return false;
224  }
225 
226  VisitorRange * range;
227  VisitorConstIterator ibegin;
228  VisitorConstIterator iend;
229 
230  // Estimations
231  Clock c;
232 
233  // True
234  if( options.at( 0 ) != '0' )
235  {
236  // True Mean Curvature
237  if( properties.at( 0 ) != '0' )
238  {
239  trace.beginBlock( "True mean curvature" );
240 
241  char full_filename[360];
242  sprintf( full_filename, "%s%s", filename.c_str(), "_True_mean.dat" );
243  std::ofstream file( full_filename );
244  file << "# h = " << h << std::endl;
245  file << "# True Mean Curvature estimation" << std::endl;
246 
247  std::ostream_iterator< double > out_it_true_mean( file, "\n" );
248 
249  range = new VisitorRange( new Visitor( surf, *surf.begin() ));
250  ibegin = range->begin();
251  iend = range->end();
252 
253  c.startClock();
254 
255  estimateTrueMeanCurvatureQuantity( ibegin,
256  iend,
257  out_it_true_mean,
258  K,
259  h,
260  aShape );
261 
262  double TTrueMeanCurv = c.stopClock();
263  file << "# time = " << TTrueMeanCurv << std::endl;
264 
265  file.close();
266  delete range;
267 
268  trace.endBlock();
269  }
270 
271  // True Gaussian Curvature
272  if( properties.at( 1 ) != '0' )
273  {
274  trace.beginBlock( "True Gausian curvature" );
275 
276  char full_filename[360];
277  sprintf( full_filename, "%s%s", filename.c_str(), "_True_gaussian.dat" );
278  std::ofstream file( full_filename );
279  file << "# h = " << h << std::endl;
280  file << "# True Gaussian Curvature estimation" << std::endl;
281 
282  std::ostream_iterator< double > out_it_true_gaussian( file, "\n" );
283 
284  range = new VisitorRange( new Visitor( surf, *surf.begin() ));
285  ibegin = range->begin();
286  iend = range->end();
287 
288  c.startClock();
289 
290  estimateTrueGaussianCurvatureQuantity( ibegin,
291  iend,
292  out_it_true_gaussian,
293  K,
294  h,
295  aShape );
296 
297  double TTrueGaussianCurv = c.stopClock();
298  file << "# time = " << TTrueGaussianCurv << std::endl;
299 
300  file.close();
301  delete range;
302 
303  trace.endBlock();
304  }
305 
306  // True Principal Curvatures
307  if( properties.at( 2 ) != '0' )
308  {
309  trace.beginBlock( "True principal curvatures" );
310 
311  char full_filename[360];
312  sprintf( full_filename, "%s%s", filename.c_str(), "_True_principal.dat" );
313  std::ofstream file( full_filename );
314  file << "# h = " << h << std::endl;
315  file << "# True Gaussian Curvature estimation" << std::endl;
316 
317  std::ostream_iterator< std::string > out_it_true_pc( file, "\n" );
318 
319  std::vector<PrincipalCurvatures> v_results;
320  std::back_insert_iterator< std::vector<PrincipalCurvatures> > bkIt(v_results);
321  range = new VisitorRange( new Visitor( surf, *surf.begin() ));
322  ibegin = range->begin();
323  iend = range->end();
324 
325  c.startClock();
326 
327  estimateTruePrincipalCurvaturesQuantity( ibegin,
328  iend,
329  bkIt,//out_it_true_pc,
330  K,
331  h,
332  aShape );
333 
334  for(unsigned int ii = 0; ii < v_results.size(); ++ii )
335  {
336  std::stringstream ss;
337  ss << v_results[ii].first << " " << v_results[ii].second;
338  *out_it_true_pc = ss.str();
339  ++out_it_true_pc;
340  }
341  double TTruePrincCurv = c.stopClock();
342  file << "# time = " << TTruePrincCurv << std::endl;
343 
344  file.close();
345 
346  delete range;
347 
348  trace.endBlock();
349  }
350  }
351 
352  double re = radius_kernel * std::pow( h, alpha ); // to obtains convergence results, re must follow the rule re=kh^(1/3)
353 
354  // II
355  if( options.at( 1 ) != '0' )
356  {
357  // Integral Invariant Mean Curvature
358  if( properties.at( 0 ) != '0' )
359  {
360  trace.beginBlock( "II mean curvature" );
361 
362  char full_filename[360];
363  sprintf( full_filename, "%s%s", filename.c_str(), "_II_mean.dat" );
364  std::ofstream file( full_filename );
365  file << "# h = " << h << std::endl;
366  file << "# Mean Curvature estimation from the Integral Invariant" << std::endl;
367  file << "# computed kernel radius = " << re << std::endl;
368 
369  range = new VisitorRange( new Visitor( surf, *surf.begin() ));
370  ibegin = range->begin();
371  iend = range->end();
372 
373  c.startClock();
374 
375  typedef functors::IIMeanCurvature3DFunctor<Z3i::Space> MyIICurvatureFunctor;
376  typedef IntegralInvariantVolumeEstimator< KSpace, KanungoPredicate, MyIICurvatureFunctor > MyIICurvatureEstimator;
377 
378  MyIICurvatureFunctor curvatureFunctor;
379  curvatureFunctor.init( h, re );
380 
381  MyIICurvatureEstimator curvatureEstimator( curvatureFunctor );
382  curvatureEstimator.attach( K, *noisifiedObject );
383  curvatureEstimator.setParams( re/h );
384  curvatureEstimator.init( h, ibegin, iend );
385 
386  std::ostream_iterator< double > out_it_ii_mean( file, "\n" );
387  curvatureEstimator.eval( ibegin, iend, out_it_ii_mean );
388 
389  double TIIMeanCurv = c.stopClock();
390  file << "# time = " << TIIMeanCurv << std::endl;
391  file.close();
392 
393  delete range;
394 
395  trace.endBlock();
396  }
397 
398  // Integral Invariant Gaussian Curvature
399  if( properties.at( 1 ) != '0' )
400  {
401  trace.beginBlock( "II Gaussian curvature" );
402 
403  char full_filename[360];
404  sprintf( full_filename, "%s%s", filename.c_str(), "_II_gaussian.dat" );
405  std::ofstream file( full_filename );
406  file << "# h = " << h << std::endl;
407  file << "# Gaussian Curvature estimation from the Integral Invariant" << std::endl;
408  file << "# computed kernel radius = " << re << std::endl;
409 
410  range = new VisitorRange( new Visitor( surf, *surf.begin() ));
411  ibegin = range->begin();
412  iend = range->end();
413 
414  c.startClock();
415 
416  typedef functors::IIGaussianCurvature3DFunctor<Z3i::Space> MyIICurvatureFunctor;
417  typedef IntegralInvariantCovarianceEstimator< KSpace, KanungoPredicate, MyIICurvatureFunctor > MyIICurvatureEstimator;
418 
419  MyIICurvatureFunctor curvatureFunctor;
420  curvatureFunctor.init( h, re );
421 
422  MyIICurvatureEstimator curvatureEstimator( curvatureFunctor );
423  curvatureEstimator.attach( K, *noisifiedObject );
424  curvatureEstimator.setParams( re/h );
425  curvatureEstimator.init( h, ibegin, iend );
426 
427  std::ostream_iterator< double > out_it_ii_gaussian( file, "\n" );
428  curvatureEstimator.eval( ibegin, iend, out_it_ii_gaussian );
429 
430  double TIIGaussCurv = c.stopClock();
431  file << "# time = " << TIIGaussCurv << std::endl;
432  file.close();
433 
434  delete range;
435 
436  trace.endBlock();
437  }
438 
439  // Integral Invariant Principal Curvatures
440  if( properties.at( 2 ) != '0' )
441  {
442  trace.beginBlock( "II Principal curvatures" );
443 
444  char full_filename[360];
445  sprintf( full_filename, "%s%s", filename.c_str(), "_II_principal.dat" );
446  std::ofstream file( full_filename );
447  file << "# h = " << h << std::endl;
448  file << "# Gaussian Curvature estimation from the Integral Invariant" << std::endl;
449  file << "# computed kernel radius = " << re << std::endl;
450 
451  range = new VisitorRange( new Visitor( surf, *surf.begin() ));
452  ibegin = range->begin();
453  iend = range->end();
454 
455  c.startClock();
456 
457  typedef functors::IIPrincipalCurvatures3DFunctor<Z3i::Space> MyIICurvatureFunctor;
458  typedef IntegralInvariantCovarianceEstimator< KSpace, KanungoPredicate, MyIICurvatureFunctor > MyIICurvatureEstimator;
459 
460  MyIICurvatureFunctor curvatureFunctor;
461  curvatureFunctor.init( h, re );
462 
463  MyIICurvatureEstimator curvatureEstimator( curvatureFunctor );
464  curvatureEstimator.attach( K, *noisifiedObject );
465  curvatureEstimator.setParams( re/h );
466  curvatureEstimator.init( h, ibegin, iend );
467 
468  std::vector<PrincipalCurvatures> v_results;
469  std::back_insert_iterator< std::vector<PrincipalCurvatures> > bkIt(v_results);
470 
471  curvatureEstimator.eval( ibegin, iend, bkIt );
472 
473  std::ostream_iterator< std::string > out_it_ii_principal( file, "\n" );
474  for( unsigned int ii = 0; ii < v_results.size(); ++ii )
475  {
476  std::stringstream ss;
477  ss << v_results[ii].first << " " << v_results[ii].second;
478  *out_it_ii_principal = ss.str();
479  ++out_it_ii_principal;
480  }
481 
482  double TIIGaussCurv = c.stopClock();
483  file << "# time = " << TIIGaussCurv << std::endl;
484  file.close();
485 
486  delete range;
487 
488  trace.endBlock();
489  }
490  }
491 
492  // Monge
493  if( options.at( 2 ) != '0' )
494  {
495  // Monge Mean Curvature
496  if( properties.at( 0 ) != '0' )
497  {
498  trace.beginBlock( "Monge mean curvature" );
499 
500  typedef functors::MongeJetFittingMeanCurvatureEstimator<Surfel, CanonicSCellEmbedder<KSpace> > FunctorMean;
501  typedef functors::ConstValue< double > ConvFunctor;
502  typedef LocalEstimatorFromSurfelFunctorAdapter<typename MyDigitalSurface::DigitalSurfaceContainer, Z3i::L2Metric, FunctorMean, ConvFunctor> ReporterH;
503  CanonicSCellEmbedder<KSpace> embedder( K );
504  FunctorMean estimatorH( embedder, h );
505  ConvFunctor convFunc(1.0);
506  ReporterH reporterH;
507  reporterH.attach( surf );
508  reporterH.setParams( Z3i::l2Metric, estimatorH, convFunc, re/h );
509 
510  c.startClock();
511  range = new VisitorRange( new Visitor( surf, *surf.begin() ));
512  ibegin = range->begin();
513  iend = range->end();
514 
515  reporterH.init( h , ibegin, iend );
516 
517  char full_filename[360];
518  sprintf( full_filename, "%s%s", filename.c_str(), "_MongeJetFitting_mean.dat" );
519  std::ofstream file( full_filename );
520  file << "# h = " << h << std::endl;
521  file << "# Mean Curvature estimation from CGAL Monge from and Jet Fitting" << std::endl;
522  file << "# computed kernel radius = " << re << std::endl;
523  std::ostream_iterator< double > out_it_monge_mean( file, "\n" );
524 
525  delete range;
526  range = new VisitorRange( new Visitor( surf, *surf.begin() ));
527  ibegin = range->begin();
528  iend = range->end();
529  //typename ReporterH::SurfelConstIterator aabegin = surf.begin();
530  //typename ReporterH::SurfelConstIterator aaend = surf.end();
531  reporterH.eval(ibegin, iend, out_it_monge_mean);
532  double TMongeMeanCurv = c.stopClock();
533  file << "# time = " << TMongeMeanCurv << std::endl;
534  file.close();
535  delete range;
536 
537 
538  trace.endBlock();
539  }
540 
541  // Monge Gaussian Curvature
542  if( properties.at( 1 ) != '0' )
543  {
544  trace.beginBlock( "Monge Gaussian curvature" );
545 
546  typedef functors::MongeJetFittingGaussianCurvatureEstimator<Surfel, CanonicSCellEmbedder<KSpace> > FunctorGaussian;
547  typedef functors::ConstValue< double > ConvFunctor;
548  typedef LocalEstimatorFromSurfelFunctorAdapter<typename MyDigitalSurface::DigitalSurfaceContainer, Z3i::L2Metric, FunctorGaussian, ConvFunctor> ReporterK;
549  CanonicSCellEmbedder<KSpace> embedder( K );
550  FunctorGaussian estimatorK( embedder, h );
551  ConvFunctor convFunc(1.0);
552  ReporterK reporterK;
553  reporterK.attach( surf );
554  reporterK.setParams( Z3i::l2Metric, estimatorK, convFunc, re/h );
555  c.startClock();
556 
557  range = new VisitorRange( new Visitor( surf, *surf.begin() ));
558  ibegin = range->begin();
559  iend = range->end();
560 
561  reporterK.init( h , ibegin, iend );
562 
563  //typename ReporterK::SurfelConstIterator aaabegin = surf.begin();
564  //typename ReporterK::SurfelConstIterator aaaend = surf.end();
565 
566  delete range;
567  range = new VisitorRange( new Visitor( surf, *surf.begin() ));
568  ibegin = range->begin();
569  iend = range->end();
570 
571  char full_filename[360];
572  sprintf( full_filename, "%s%s", filename.c_str(), "_MongeJetFitting_gaussian.dat" );
573  std::ofstream file( full_filename );
574  file << "# h = " << h << std::endl;
575  file << "# Gaussian Curvature estimation from CGAL Monge from and Jet Fitting" << std::endl;
576  file << "# computed kernel radius = " << re << std::endl;
577  std::ostream_iterator< double > out_it_monge_gaussian( file, "\n" );
578  reporterK.eval(ibegin, iend , out_it_monge_gaussian);
579  double TMongeGaussCurv = c.stopClock();
580  file << "# time = " << TMongeGaussCurv << std::endl;
581  file.close();
582  delete range;
583 
584 
585  trace.endBlock();
586  }
587 
588  // Monge Principal Curvatures
589  if( properties.at( 2 ) != '0' )
590  {
591  trace.beginBlock( "Monge Principal Curvature" );
592 
593  typedef functors::MongeJetFittingPrincipalCurvaturesEstimator<Surfel, CanonicSCellEmbedder<KSpace> > FunctorPrincCurv;
594  typedef functors::ConstValue< double > ConvFunctor;
595  typedef LocalEstimatorFromSurfelFunctorAdapter<typename MyDigitalSurface::DigitalSurfaceContainer, Z3i::L2Metric, FunctorPrincCurv, ConvFunctor> ReporterK;
596  CanonicSCellEmbedder<KSpace> embedder( K );
597  FunctorPrincCurv estimatorK( embedder, h );
598  ConvFunctor convFunc(1.0);
599  ReporterK reporterK;
600  reporterK.attach( surf );
601  reporterK.setParams( Z3i::l2Metric, estimatorK, convFunc, re/h );
602 
603  c.startClock();
604 
605  range = new VisitorRange( new Visitor( surf, *surf.begin() ));
606  ibegin = range->begin();
607  iend = range->end();
608  reporterK.init( h , ibegin, iend );
609 
610  char full_filename[360];
611  sprintf( full_filename, "%s%s", filename.c_str(), "_MongeJetFitting_principal.dat" );
612  std::ofstream file( full_filename );
613  file << "# h = " << h << std::endl;
614  file << "# Gaussian Curvature estimation from CGAL Monge from and Jet Fitting" << std::endl;
615  file << "# computed kernel radius = " << re << std::endl;
616  std::ostream_iterator< std::string > out_it_monge_principal( file, "\n" );
617 
618  std::vector<PrincipalCurvatures> v_results;
619  std::back_insert_iterator< std::vector<PrincipalCurvatures> > bkIt(v_results);
620 
621  delete range;
622  range = new VisitorRange( new Visitor( surf, *surf.begin() ));
623  ibegin = range->begin();
624  iend = range->end();
625 
626  reporterK.eval(ibegin, iend , bkIt);//out_it_monge_principal);
627 
628  for(unsigned int ii = 0; ii < v_results.size(); ++ii )
629  {
630  std::stringstream ss;
631  ss << v_results[ii].first << " " << v_results[ii].second;
632  *out_it_monge_principal = ss.str();
633  ++out_it_monge_principal;
634  }
635 
636  double TMongeGaussCurv = c.stopClock();
637  file << "# time = " << TMongeGaussCurv << std::endl;
638  file.close();
639  delete range;
640 
641 
642  trace.endBlock();
643  }
644  }
645  }
646  else // no noise
647  {
648  typedef LightImplicitDigitalSurface< KSpace, DigitalShape > Boundary;
649  typedef DigitalSurface< Boundary > MyDigitalSurface;
650  typedef typename MyDigitalSurface::ConstIterator ConstIterator;
651 
652  typedef DepthFirstVisitor< MyDigitalSurface > Visitor;
653  typedef GraphVisitorRange< Visitor > VisitorRange;
654  typedef typename VisitorRange::ConstIterator VisitorConstIterator;
655 
656  // typedef PointFunctorFromPointPredicateAndDomain< DigitalShape, Z3i::Domain, unsigned int > MyPointFunctor;
657  // typedef FunctorOnCells< MyPointFunctor, KSpace > MySpelFunctor;
658 
659  // Extracts shape boundary
660  SCell bel = Surfaces<KSpace>::findABel ( K, *dshape, 10000 );
661  Boundary boundary( K, *dshape, SurfelAdjacency< KSpace::dimension >( true ), bel );
662  MyDigitalSurface surf ( boundary );
663 
664  VisitorRange * range;
665  VisitorConstIterator ibegin;
666  VisitorConstIterator iend;
667 
668  unsigned int cntIn = 0;
669  for( typename Z3i::Domain::ConstIterator it = dshape->getDomain().begin(), ite = dshape->getDomain().end(); it != ite; ++it )
670  {
671  if( dshape->operator ()(*it))
672  {
673  ++cntIn;
674  }
675  }
676 
677  std::cout << "boundary:" << surf.size() << std::endl;
678  std::cout << "full:" << cntIn << std::endl;
679 
680  // Estimations
681  Clock c;
682 
683  // True
684  if( options.at( 0 ) != '0' )
685  {
686  // True Mean Curvature
687  if( properties.at( 0 ) != '0' )
688  {
689  trace.beginBlock( "True mean curvature" );
690 
691  char full_filename[360];
692  sprintf( full_filename, "%s%s", filename.c_str(), "_True_mean.dat" );
693  std::ofstream file( full_filename );
694  file << "# h = " << h << std::endl;
695  file << "# True Mean Curvature estimation" << std::endl;
696 
697  std::ostream_iterator< double > out_it_true_mean( file, "\n" );
698 
699  range = new VisitorRange( new Visitor( surf, *surf.begin() ));
700  ibegin = range->begin();
701  iend = range->end();
702 
703  c.startClock();
704 
705  estimateTrueMeanCurvatureQuantity( ibegin,
706  iend,
707  out_it_true_mean,
708  K,
709  h,
710  aShape );
711 
712  double TTrueMeanCurv = c.stopClock();
713  file << "# time = " << TTrueMeanCurv << std::endl;
714 
715  file.close();
716  delete range;
717 
718  trace.endBlock();
719  }
720 
721  // True Gaussian Curvature
722  if( properties.at( 1 ) != '0' )
723  {
724  trace.beginBlock( "True Gaussian curvature" );
725 
726  char full_filename[360];
727  sprintf( full_filename, "%s%s", filename.c_str(), "_True_gaussian.dat" );
728  std::ofstream file( full_filename );
729  file << "# h = " << h << std::endl;
730  file << "# True Gaussian Curvature estimation" << std::endl;
731 
732  std::ostream_iterator< double > out_it_true_gaussian( file, "\n" );
733 
734  range = new VisitorRange( new Visitor( surf, *surf.begin() ));
735  ibegin = range->begin();
736  iend = range->end();
737 
738  c.startClock();
739 
740  estimateTrueGaussianCurvatureQuantity( ibegin,
741  iend,
742  out_it_true_gaussian,
743  K,
744  h,
745  aShape );
746 
747  double TTrueGaussianCurv = c.stopClock();
748  file << "# time = " << TTrueGaussianCurv << std::endl;
749 
750  file.close();
751 
752  delete range;
753 
754  trace.endBlock();
755  }
756 
757  // True Principal Curvatures
758  if( properties.at( 2 ) != '0' )
759  {
760  trace.beginBlock( "True principal curvatures" );
761 
762  char full_filename[360];
763  sprintf( full_filename, "%s%s", filename.c_str(), "_True_principal.dat" );
764  std::ofstream file( full_filename );
765  file << "# h = " << h << std::endl;
766  file << "# True Gaussian Curvature estimation" << std::endl;
767 
768  std::ostream_iterator< std::string > out_it_true_pc( file, "\n" );
769 
770  std::vector<PrincipalCurvatures> v_results;
771  std::back_insert_iterator< std::vector<PrincipalCurvatures> > bkIt(v_results);
772 
773  range = new VisitorRange( new Visitor( surf, *surf.begin() ));
774  ibegin = range->begin();
775  iend = range->end();
776 
777  c.startClock();
778 
779  estimateTruePrincipalCurvaturesQuantity( ibegin,
780  iend,
781  bkIt,// out_it_true_pc,
782  K,
783  h,
784  aShape );
785 
786 
787  for(unsigned int ii = 0; ii < v_results.size(); ++ii )
788  {
789  std::stringstream ss;
790  ss << v_results[ii].first << " " << v_results[ii].second;
791  *out_it_true_pc = ss.str();
792  ++out_it_true_pc;
793  }
794 
795  double TTruePrincCurv = c.stopClock();
796  file << "# time = " << TTruePrincCurv << std::endl;
797 
798  file.close();
799 
800  delete range;
801 
802  trace.endBlock();
803  }
804  }
805 
806  double re = radius_kernel * std::pow( h, alpha ); // to obtains convergence results, re must follow the rule re=kh^(1/3)
807 
808  // II
809  if( options.at( 1 ) != '0' )
810  {
811  // Integral Invariant Mean Curvature
812  if( properties.at( 0 ) != '0' )
813  {
814  trace.beginBlock( "II mean curvature" );
815 
816  char full_filename[360];
817  sprintf( full_filename, "%s%s", filename.c_str(), "_II_mean.dat" );
818  std::ofstream file( full_filename );
819  file << "# h = " << h << std::endl;
820  file << "# Mean Curvature estimation from the Integral Invariant" << std::endl;
821  file << "# computed kernel radius = " << re << std::endl;
822 
823  range = new VisitorRange( new Visitor( surf, *surf.begin() ));
824  ibegin = range->begin();
825  iend = range->end();
826 
827  c.startClock();
828 
829  typedef functors::IIMeanCurvature3DFunctor<Z3i::Space> MyIICurvatureFunctor;
830  typedef IntegralInvariantVolumeEstimator< KSpace, DigitalShape, MyIICurvatureFunctor > MyIICurvatureEstimator;
831 
832  MyIICurvatureFunctor curvatureFunctor;
833  curvatureFunctor.init( h, re );
834 
835  MyIICurvatureEstimator curvatureEstimator( curvatureFunctor );
836  curvatureEstimator.attach( K, *dshape );
837  curvatureEstimator.setParams( re/h );
838  curvatureEstimator.init( h, ibegin, iend );
839 
840  std::ostream_iterator< double > out_it_ii_mean( file, "\n" );
841  curvatureEstimator.eval( ibegin, iend, out_it_ii_mean );
842 
843  double TIIMeanCurv = c.stopClock();
844  file << "# time = " << TIIMeanCurv << std::endl;
845  file.close();
846 
847  delete range;
848 
849  trace.endBlock();
850  }
851 
852  // Integral Invariant Gaussian Curvature
853  if( properties.at( 1 ) != '0' )
854  {
855  trace.beginBlock( "II Gaussian curvature" );
856 
857  char full_filename[360];
858  sprintf( full_filename, "%s%s", filename.c_str(), "_II_gaussian.dat" );
859  std::ofstream file( full_filename );
860  file << "# h = " << h << std::endl;
861  file << "# Gaussian Curvature estimation from the Integral Invariant" << std::endl;
862  file << "# computed kernel radius = " << re << std::endl;
863 
864  range = new VisitorRange( new Visitor( surf, *surf.begin() ));
865  ibegin = range->begin();
866  iend = range->end();
867 
868  c.startClock();
869 
870  typedef functors::IIGaussianCurvature3DFunctor<Z3i::Space> MyIICurvatureFunctor;
871  typedef IntegralInvariantCovarianceEstimator< KSpace, DigitalShape, MyIICurvatureFunctor > MyIICurvatureEstimator;
872 
873  MyIICurvatureFunctor curvatureFunctor;
874  curvatureFunctor.init( h, re );
875 
876  MyIICurvatureEstimator curvatureEstimator( curvatureFunctor );
877  curvatureEstimator.attach( K, *dshape );
878  curvatureEstimator.setParams( re/h );
879  curvatureEstimator.init( h, ibegin, iend );
880 
881  std::ostream_iterator< double > out_it_ii_gaussian( file, "\n" );
882  curvatureEstimator.eval( ibegin, iend, out_it_ii_gaussian );
883 
884  double TIIGaussCurv = c.stopClock();
885  file << "# time = " << TIIGaussCurv << std::endl;
886  file.close();
887 
888  delete range;
889 
890  trace.endBlock();
891  }
892 
893  // Integral Invariant Principal Curvatures
894  if( properties.at( 2 ) != '0' )
895  {
896  trace.beginBlock( "II Principal curvatures" );
897 
898  char full_filename[360];
899  sprintf( full_filename, "%s%s", filename.c_str(), "_II_principal.dat" );
900  std::ofstream file( full_filename );
901  file << "# h = " << h << std::endl;
902  file << "# Gaussian Curvature estimation from the Integral Invariant" << std::endl;
903  file << "# computed kernel radius = " << re << std::endl;
904 
905  range = new VisitorRange( new Visitor( surf, *surf.begin() ));
906  ibegin = range->begin();
907  iend = range->end();
908 
909  c.startClock();
910 
911  typedef functors::IIPrincipalCurvatures3DFunctor<Z3i::Space> MyIICurvatureFunctor;
912  typedef IntegralInvariantCovarianceEstimator< KSpace, DigitalShape, MyIICurvatureFunctor > MyIICurvatureEstimator;
913 
914  MyIICurvatureFunctor curvatureFunctor;
915  curvatureFunctor.init( h, re );
916 
917  MyIICurvatureEstimator curvatureEstimator( curvatureFunctor );
918  curvatureEstimator.attach( K, *dshape );
919  curvatureEstimator.setParams( re/h );
920  curvatureEstimator.init( h, ibegin, iend );
921 
922  std::vector<PrincipalCurvatures> v_results;
923  std::back_insert_iterator< std::vector<PrincipalCurvatures> > bkIt(v_results);
924 
925  curvatureEstimator.eval( ibegin, iend, bkIt );
926 
927  std::ostream_iterator< std::string > out_it_ii_principal( file, "\n" );
928  for( unsigned int ii = 0; ii < v_results.size(); ++ii )
929  {
930  std::stringstream ss;
931  ss << v_results[ii].first << " " << v_results[ii].second;
932  *out_it_ii_principal = ss.str();
933  ++out_it_ii_principal;
934  }
935 
936  double TIIGaussCurv = c.stopClock();
937  file << "# time = " << TIIGaussCurv << std::endl;
938  file.close();
939 
940  delete range;
941 
942  trace.endBlock();
943  }
944  }
945 
946  // Monge
947  if( options.at( 2 ) != '0' )
948  {
949  // Monge Mean Curvature
950  if( properties.at( 0 ) != '0' )
951  {
952  trace.beginBlock( "Monge mean curvature" );
953 
954  typedef functors::MongeJetFittingMeanCurvatureEstimator<Surfel, CanonicSCellEmbedder<KSpace> > FunctorMean;
955  typedef functors::ConstValue< double > ConvFunctor;
956  typedef LocalEstimatorFromSurfelFunctorAdapter<typename MyDigitalSurface::DigitalSurfaceContainer, Z3i::L2Metric, FunctorMean, ConvFunctor> ReporterH;
957  CanonicSCellEmbedder<KSpace> embedder( K );
958  FunctorMean estimatorH( embedder, h );
959  ConvFunctor convFunc(1.0);
960  ReporterH reporterH;
961  reporterH.attach( surf );
962  reporterH.setParams( Z3i::l2Metric, estimatorH, convFunc, re/h );
963  c.startClock();
964 
965  range = new VisitorRange( new Visitor( surf, *surf.begin() ));
966  ibegin = range->begin();
967  iend = range->end();
968  reporterH.init( h , ibegin, iend );
969 
970  char full_filename[360];
971  sprintf( full_filename, "%s%s", filename.c_str(), "_MongeJetFitting_mean.dat" );
972  std::ofstream file( full_filename );
973  file << "# h = " << h << std::endl;
974  file << "# Mean Curvature estimation from CGAL Monge from and Jet Fitting" << std::endl;
975  file << "# computed kernel radius = " << re << std::endl;
976  std::ostream_iterator< double > out_it_monge_mean( file, "\n" );
977 
978  delete range;
979  range = new VisitorRange( new Visitor( surf, *surf.begin() ));
980  ibegin = range->begin();
981  iend = range->end();
982 
983  reporterH.eval(ibegin, iend , out_it_monge_mean);
984  double TMongeMeanCurv = c.stopClock();
985  file << "# time = " << TMongeMeanCurv << std::endl;
986  file.close();
987  delete range;
988 
989 
990  trace.endBlock();
991  }
992 
993  // Monge Gaussian Curvature
994  if( properties.at( 1 ) != '0' )
995  {
996  trace.beginBlock( "Monge Gaussian curvature" );
997 
998  typedef functors::MongeJetFittingGaussianCurvatureEstimator<Surfel, CanonicSCellEmbedder<KSpace> > FunctorGaussian;
999  typedef functors::ConstValue< double > ConvFunctor;
1000  typedef LocalEstimatorFromSurfelFunctorAdapter<typename MyDigitalSurface::DigitalSurfaceContainer, Z3i::L2Metric, FunctorGaussian, ConvFunctor> ReporterK;
1001  CanonicSCellEmbedder<KSpace> embedder( K );
1002  FunctorGaussian estimatorK( embedder, h );
1003  ConvFunctor convFunc(1.0);
1004  ReporterK reporterK;
1005  reporterK.attach( surf );
1006  reporterK.setParams( Z3i::l2Metric, estimatorK, convFunc, re/h );
1007 
1008  c.startClock();
1009 
1010  range = new VisitorRange( new Visitor( surf, *surf.begin() ));
1011  ibegin = range->begin();
1012  iend = range->end();
1013  reporterK.init( h , ibegin, iend );
1014 
1015  delete range;
1016  range = new VisitorRange( new Visitor( surf, *surf.begin() ));
1017  ibegin = range->begin();
1018  iend = range->end();
1019 
1020  char full_filename[360];
1021  sprintf( full_filename, "%s%s", filename.c_str(), "_MongeJetFitting_gaussian.dat" );
1022  std::ofstream file( full_filename );
1023  file << "# h = " << h << std::endl;
1024  file << "# Gaussian Curvature estimation from CGAL Monge from and Jet Fitting" << std::endl;
1025  file << "# computed kernel radius = " << re << std::endl;
1026  std::ostream_iterator< double > out_it_monge_gaussian( file, "\n" );
1027  reporterK.eval(ibegin, iend , out_it_monge_gaussian);
1028  double TMongeGaussCurv = c.stopClock();
1029  file << "# time = " << TMongeGaussCurv << std::endl;
1030  file.close();
1031  delete range;
1032 
1033 
1034  trace.endBlock();
1035  }
1036 
1037  // Monge Principal Curvatures
1038  if( properties.at( 2 ) != '0' )
1039  {
1040  trace.beginBlock( "Monge Principal Curvature" );
1041 
1042  typedef functors::MongeJetFittingPrincipalCurvaturesEstimator<Surfel, CanonicSCellEmbedder<KSpace> > FunctorPrincCurv;
1043  typedef functors::ConstValue< double > ConvFunctor;
1044  typedef LocalEstimatorFromSurfelFunctorAdapter<typename MyDigitalSurface::DigitalSurfaceContainer, Z3i::L2Metric, FunctorPrincCurv, ConvFunctor> ReporterK;
1045  CanonicSCellEmbedder<KSpace> embedder( K );
1046  FunctorPrincCurv estimatorK( embedder, h );
1047  ConvFunctor convFunc(1.0);
1048  ReporterK reporterK;
1049  reporterK.attach(surf);
1050  reporterK.setParams(Z3i::l2Metric, estimatorK, convFunc, re/h);
1051 
1052  c.startClock();
1053 
1054  range = new VisitorRange( new Visitor( surf, *surf.begin() ));
1055  ibegin = range->begin();
1056  iend = range->end();
1057  reporterK.init( h , ibegin, iend );
1058 
1059  delete range;
1060  range = new VisitorRange( new Visitor( surf, *surf.begin() ));
1061  ibegin = range->begin();
1062  iend = range->end();
1063 
1064  char full_filename[360];
1065  sprintf( full_filename, "%s%s", filename.c_str(), "_MongeJetFitting_principal.dat" );
1066  std::ofstream file( full_filename );
1067  file << "# h = " << h << std::endl;
1068  file << "# Gaussian Curvature estimation from CGAL Monge from and Jet Fitting" << std::endl;
1069  file << "# computed kernel radius = " << re << std::endl;
1070  std::ostream_iterator< std::string > out_it_monge_principal( file, "\n" );
1071 
1072  std::vector<PrincipalCurvatures> v_results;
1073  std::back_insert_iterator< std::vector<PrincipalCurvatures> > bkIt(v_results);
1074 
1075  reporterK.eval(ibegin, iend , bkIt);//out_it_monge_principal);
1076 
1077  for(unsigned int ii = 0; ii < v_results.size(); ++ii )
1078  {
1079  std::stringstream ss;
1080  ss << v_results[ii].first << " " << v_results[ii].second;
1081  *out_it_monge_principal = ss.str();
1082  ++out_it_monge_principal;
1083  }
1084 
1085  double TMongeGaussCurv = c.stopClock();
1086  file << "# time = " << TMongeGaussCurv << std::endl;
1087  file.close();
1088  delete range;
1089 
1090 
1091  trace.endBlock();
1092  }
1093 
1094  }
1095  }
1096  }
1097  catch ( InputException e )
1098  {
1099  std::cerr << "[estimatorCurvatureComparator3D]"
1100  << " error."
1101  << e.what()
1102  << " "
1103  << filename << " "
1104  << border_min[0] << " "
1105  << border_max[0] << " "
1106  << h << " "
1107  << radius_kernel << " "
1108  << lambda_optimized << " "
1109  << options << " "
1110  << alpha << " "
1111  << std::endl;
1112  return false;
1113  }
1114 
1115  delete dshape;
1116 
1117  return true;
1118 }
1119 
1125 void missingParam( std::string param )
1126 {
1127  trace.error() << " Parameter: " << param << " is required.";
1128  trace.info() << std::endl;
1129  exit( 1 );
1130 }
1131 
1132 namespace po = boost::program_options;
1133 
1134 int main( int argc, char** argv )
1135 {
1136 
1137 
1138 #ifndef WITH_CGAL
1139 #error You need to have activated CGAL (WITH_CGAL) to include this file.
1140 #endif
1141 #ifndef WITH_EIGEN
1142 #error You need to have activated EIGEN (WITH_EIGEN) to include this file.
1143 #endif
1144 
1145  // parse command line ----------------------------------------------
1146  po::options_description general_opt("Allowed options are");
1147  general_opt.add_options()
1148  ("help,h", "display this message")
1149  ("shape,s", po::value< std::string >(), "Shape")
1150  ("output,o", po::value< std::string >(), "Output file")
1151  ("radius,r", po::value< double >(), "Kernel radius for IntegralInvariant" )
1152  ("alpha", po::value<double>()->default_value(1.0/3.0), "Alpha parameter for Integral Invariant computation" )
1153  ("h", po::value< double >(), "Grid step" )
1154  ("minAABB,a", po::value< double >()->default_value( -10.0 ), "Min value of the AABB bounding box (domain)" )
1155  ("maxAABB,A", po::value< double >()->default_value( 10.0 ), "Max value of the AABB bounding box (domain)" )
1156  ("noise,n", po::value<double>()->default_value(0.0), "Level of noise to perturb the shape" )
1157  ("lambda,l", po::value< bool >()->default_value( false ), "Use the shape to get a better approximation of the surface (optional)" )
1158  ("properties", po::value<std::string>()->default_value("110"), "the i-th property is disabled iff there is a 0 at position i" )
1159  ("estimators,e", po::value< std::string >()->default_value("110"), "the i-th estimator is disabled iff there is a 0 at position i" );
1160 
1161 
1162  bool parseOK = true;
1163  po::variables_map vm;
1164  try
1165  {
1166  po::store( po::parse_command_line( argc, argv, general_opt ), vm );
1167  }
1168  catch( const std::exception & ex )
1169  {
1170  parseOK = false;
1171  trace.info() << "Error checking program options: " << ex.what() << std::endl;
1172  }
1173  po::notify( vm );
1174  if( !parseOK || vm.count("help") || argc <= 1 )
1175  {
1176  trace.info()<< "Compare local estimators on implicit shapes using DGtal library" <<std::endl
1177  << "Basic usage: "<<std::endl
1178  << "\t3dlocalEstimators --shape <shape> --h <h> --radius <radius> --estimators <binaryWord> --output <output>"<<std::endl
1179  << std::endl
1180  << "Below are the different available families of estimators: " << std::endl
1181  << "\t - Integral Invariant Mean" << std::endl
1182  << "\t - Integral Invariant Gaussian" << std::endl
1183  << "\t - Monge Jet Fitting Mean" << std::endl
1184  << "\t - Monge Jet Fitting Gaussian" << std::endl
1185  << std::endl
1186  << "The i-th family of estimators is enabled if the i-th character of the binary word is not 0. "
1187  << "The default binary word is '1100'. This means that the first family of estimators, "
1188  << "ie. Integral Invariant, is enabled, whereas the next ones are disabled. "
1189  << "Below are the different available properties: " << std::endl
1190  << "\t - Mean Curvature" << std::endl
1191  << "\t - Gaussian Curvature" << std::endl
1192  << "\t - k1/k2" << std::endl
1193  << std::endl;
1194  return 0;
1195  }
1196 
1197  if (!(vm.count("output"))) missingParam("--output");
1198  if (!(vm.count("shape"))) missingParam("--shape");
1199  if (!(vm.count("h"))) missingParam("--h");
1200  if (!(vm.count("radius"))) missingParam("--radius");
1201 
1202 
1203  std::string file_export = vm["output"].as< std::string >();
1204  int nb = 3;
1205  std::string options = vm["estimators"].as< std::string >();
1206  if (options.size() < nb)
1207  {
1208  trace.error() << " At least " << nb
1209  << " characters are required "
1210  << " with option --estimators.";
1211  trace.info() << std::endl;
1212  exit(1);
1213  }
1214  double h = vm["h"].as< double >();
1215  double radius = vm["radius"].as< double >();
1216  double alpha = vm["alpha"].as< double >();
1217  std::string poly_str = vm["shape"].as< std::string >();
1218  bool lambda_optimized = vm["lambda"].as< bool >();
1219  double noiseLevel = vm["noise"].as<double>();
1220 
1221  nb = 3; //number of available properties
1222  std::string properties = vm["properties"].as<std::string>();
1223  if (properties.size() < nb)
1224  {
1225  trace.error() << " At least " << nb
1226  << " characters are required "
1227  << " with option --properties.";
1228  trace.info() << std::endl;
1229  exit(1);
1230  }
1231 
1232  typedef Z3i::Space::RealPoint RealPoint;
1233  typedef Z3i::Space::RealPoint::Coordinate Ring;
1234 
1235  RealPoint border_min( vm["minAABB"].as< double >(), vm["minAABB"].as< double >(), vm["minAABB"].as< double >() );
1236  RealPoint border_max( vm["maxAABB"].as< double >(), vm["maxAABB"].as< double >(), vm["maxAABB"].as< double >() );
1237 
1239 
1240  typedef MPolynomial< 3, Ring > Polynomial3;
1241  typedef MPolynomialReader<3, Ring> Polynomial3Reader;
1242  typedef ImplicitPolynomial3Shape<Z3i::Space> ImplicitShape;
1243 
1244  Polynomial3 poly;
1245  Polynomial3Reader reader;
1246  std::string::const_iterator iter = reader.read( poly, poly_str.begin(), poly_str.end() );
1247  if ( iter != poly_str.end() )
1248  {
1249  std::cerr << "ERROR: I read only <"
1250  << poly_str.substr( 0, iter - poly_str.begin() )
1251  << ">, and I built P=" << poly << std::endl;
1252  return 1;
1253  }
1254 
1255  ImplicitShape* shape = new ImplicitShape( poly );
1256 
1257  compareShapeEstimators< Z3i::Space, ImplicitShape > (
1258  file_export,
1259  shape,
1260  border_min, border_max,
1261  h,
1262  radius,
1263  alpha,
1264  options,
1265  properties,
1266  lambda_optimized,
1267  noiseLevel );
1268 
1269  delete shape;
1270 }
Companion project associated to DGtal ; Generated on Mon Jun 27 2016 12:21:49 for DGtalTools by   doxygen 1.8.11