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