DGtalTools  1.2.0
generic3dNormalEstimators.cpp
1 
32 #include <iostream>
33 #include <fstream>
34 #include <sstream>
35 
36 #include "CLI11.hpp"
37 
38 #include "DGtal/base/Common.h"
39 #include "DGtal/base/CountedPtr.h"
40 #include "DGtal/base/CountedConstPtrOrConstPtr.h"
41 #include "DGtal/helpers/StdDefs.h"
42 #include "DGtal/math/Statistic.h"
43 #include "DGtal/math/MPolynomial.h"
44 #include "DGtal/topology/LightImplicitDigitalSurface.h"
45 #include "DGtal/graph/DepthFirstVisitor.h"
46 #include "DGtal/graph/GraphVisitorRange.h"
47 #include "DGtal/geometry/surfaces/estimation/CNormalVectorEstimator.h"
48 #include "DGtal/geometry/surfaces/estimation/VoronoiCovarianceMeasureOnDigitalSurface.h"
49 #include "DGtal/geometry/surfaces/estimation/VCMDigitalSurfaceLocalEstimator.h"
50 #include "DGtal/geometry/surfaces/estimation/TrueDigitalSurfaceLocalEstimator.h"
51 #include "DGtal/geometry/surfaces/estimation/IIGeometricFunctors.h"
52 #include "DGtal/geometry/surfaces/estimation/IntegralInvariantCovarianceEstimator.h"
53 #include "DGtal/geometry/volumes/KanungoNoise.h"
54 #include "DGtal/shapes/GaussDigitizer.h"
55 #include "DGtal/shapes/ShapeGeometricFunctors.h"
56 #include "DGtal/shapes/implicit/ImplicitPolynomial3Shape.h"
57 #include "DGtal/images/ImageContainerBySTLVector.h"
58 #include "DGtal/images/SimpleThresholdForegroundPredicate.h"
59 #include "DGtal/io/readers/MPolynomialReader.h"
60 #include "DGtal/io/colormaps/GradientColorMap.h"
61 #ifdef WITH_VISU3D_QGLVIEWER
62 #include "DGtal/io/viewers/Viewer3D.h"
63 #include "DGtal/io/Display3DFactory.h"
64 #endif
65 
66 using namespace std;
67 using namespace DGtal;
68 
69 
70 
199 struct AllParams{
200  double noise {0.0};
201  double minAABB {-10.0};
202  double maxAABB {10.0};
203  double gridstep {1.0};
204  double R_radius {5.0};
205  double r_radius {3.0};
206  double alpha {0.0};
207  double trivial_radius {3.0};
208  int embedding {0};
209  std::string polynomials;
210  std::string estimator {"True"};
211  std::string kernel {"hat"};
212  std::string output {"output"};
213  std::string exportX {"None"};
214  std::string view {"None"};
215  bool angle_deviation_stats;
216  bool normals;
217  bool noff;
218 };
219 
220 
221 
222 template <typename SCell, typename RealVector>
223 struct GradientMapAdapter {
224  typedef std::map<SCell,RealVector> SCell2RealVectorMap;
225  typedef SCell Argument;
226  typedef RealVector Value;
227  GradientMapAdapter( ConstAlias<SCell2RealVectorMap> map )
228  : myMap( map ) {}
229  RealVector operator()( const Argument& arg ) const
230  {
231  typename SCell2RealVectorMap::const_iterator it = myMap->find( arg );
232  if ( it != myMap->end() ) return it->second;
233  else return RealVector();
234  }
235  CountedConstPtrOrConstPtr<SCell2RealVectorMap> myMap;
236 };
237 
238 template <typename SCellEmbedder>
239 struct SCellEmbedderWithNormal : public SCellEmbedder
240 {
241  using SCellEmbedder::space;
242  using SCellEmbedder::operator();
243  typedef typename SCellEmbedder::KSpace KSpace;
244  typedef typename SCellEmbedder::SCell SCell;
245  typedef typename SCellEmbedder::RealPoint RealPoint;
246  typedef SCell Argument;
247  typedef RealPoint Value;
248  typedef typename KSpace::Space::RealVector RealVector;
249  typedef std::map<SCell,RealVector> SCell2RealVectorMap;
250  typedef GradientMapAdapter<SCell,RealVector> GradientMap;
251 
252  SCellEmbedderWithNormal( ConstAlias<SCellEmbedder> embedder,
253  ConstAlias<SCell2RealVectorMap> map )
254  : SCellEmbedder( embedder ), myMap( map )
255  {}
256 
257  GradientMap gradientMap() const
258  {
259  return GradientMap( myMap );
260  }
261 
262  CountedConstPtrOrConstPtr<SCell2RealVectorMap> myMap;
263 };
264 
265 template <typename DigitalSurface,
266  typename Estimator>
267 void exportNOFFSurface( const DigitalSurface& surface,
268  const Estimator& estimator,
269  std::ostream& output )
270 {
271  typedef typename DigitalSurface::KSpace KSpace;
272  typedef typename DigitalSurface::ConstIterator ConstIterator;
273  typedef typename DigitalSurface::Surfel Surfel;
274  typedef typename Estimator::Quantity Quantity;
275  const KSpace& ks = surface.container().space();
276  std::map<Surfel,Quantity> normals;
277  for ( ConstIterator it = surface.begin(), itE = surface.end(); it != itE; ++it )
278  {
279  Quantity n_est = estimator.eval( it );
280  normals[ *it ] = n_est;
281  }
282  CanonicSCellEmbedder<KSpace> surfelEmbedder( ks );
283  typedef SCellEmbedderWithNormal< CanonicSCellEmbedder<KSpace> > Embedder;
284  Embedder embedder( surfelEmbedder, normals );
285  surface.exportAs3DNOFF( output, embedder );
286 }
287 
288 
292 template <typename KSpace,
293  typename ImplicitShape,
294  typename Surface,
295  typename TrueEstimator,
296  typename Estimator>
297 void computeEstimation
298 ( const AllParams &params, //< command-line parameters
299  const KSpace& K, //< cellular grid space
300  const ImplicitShape& shape, //< implicit shape "ground truth"
301  const Surface& surface, //< digital surface approximating shape
302  TrueEstimator& true_estimator, //< "ground truth" estimator
303  Estimator& estimator ) //< an initialized estimator
304 {
305  typedef typename Surface::Surfel Surfel;
306  typedef typename Estimator::Quantity Quantity;
307  typedef double Scalar;
308  typedef DepthFirstVisitor< Surface > Visitor;
309  typedef GraphVisitorRange< Visitor > VisitorRange;
310  typedef typename VisitorRange::ConstIterator VisitorConstIterator;
311 
312  std::string fname = params.output;
313  string nameEstimator = params.estimator;
314  trace.beginBlock( "Computing " + nameEstimator + "estimations." );
315  CountedPtr<VisitorRange> range( new VisitorRange( new Visitor( surface, *(surface.begin()) )) );
316  std::vector<Quantity> n_estimations;
317  estimator.eval( range->begin(), range->end(), std::back_inserter( n_estimations ) );
318  trace.info() << "- nb estimations = " << n_estimations.size() << std::endl;
319  trace.endBlock();
320 
321  trace.beginBlock( "Computing ground truth." );
322  range = CountedPtr<VisitorRange>( new VisitorRange( new Visitor( surface, *(surface.begin()) )) );
323  std::vector<Quantity> n_true_estimations;
324  true_estimator.eval( range->begin(), range->end(), std::back_inserter( n_true_estimations ) );
325  trace.info() << "- nb estimations = " << n_true_estimations.size() << std::endl;
326  trace.endBlock();
327 
328  trace.beginBlock( "Correcting orientations." );
329  ASSERT( n_estimations.size() == n_true_estimations.size() );
330  for ( unsigned int i = 0; i < n_estimations.size(); ++i )
331  if ( n_estimations[ i ].dot( n_true_estimations[ i ] ) < 0 )
332  n_estimations[ i ] = -n_estimations[ i ];
333  trace.endBlock();
334 
335  DGtal::GradientColorMap<double> grad( 0.0, 40.0 );
336  // 0 metallic blue, 5 light cyan, 10 light green, 15 light
337  // yellow, 20 yellow, 25 orange, 30 red, 35, dark red, 40- grey
338  grad.addColor( DGtal::Color( 128, 128, 255 ) ); // 0
339  grad.addColor( DGtal::Color( 128, 255, 255 ) ); // 5
340  grad.addColor( DGtal::Color( 128, 255, 128 ) ); // 10
341  grad.addColor( DGtal::Color( 255, 255, 128 ) ); // 15
342  grad.addColor( DGtal::Color( 255, 255, 0 ) ); // 20
343  grad.addColor( DGtal::Color( 255, 128, 0 ) ); // 25
344  grad.addColor( DGtal::Color( 255, 0, 0 ) ); // 30
345  grad.addColor( DGtal::Color( 128, 0, 0 ) ); // 35
346  grad.addColor( DGtal::Color( 128, 128, 128 ) ); // 40
347 
348  if ( params.angle_deviation_stats )
349  {
350  trace.beginBlock( "Computing angle deviation error stats." );
351  std::ostringstream adev_sstr;
352  adev_sstr << fname << "-" << nameEstimator << "-angle-deviation-"
353  << estimator.h() << ".txt";
354  DGtal::Statistic<Scalar> adev_stat;
355  unsigned int i = 0;
356  range = CountedPtr<VisitorRange>( new VisitorRange( new Visitor( surface, *(surface.begin()) )) );
357  for ( VisitorConstIterator it = range->begin(), itE = range->end(); it != itE; ++it, ++i )
358  {
359  Quantity n_est = n_estimations[ i ];
360  Quantity n_true_est = n_true_estimations[ i ];
361  Scalar angle_error = acos( n_est.dot( n_true_est ) );
362  adev_stat.addValue( angle_error );
363  }
364  adev_stat.terminate();
365  std::ofstream adev_output( adev_sstr.str().c_str() );
366  adev_output << "# Average error X of the absolute angle between two vector estimations." << std::endl;
367  adev_output << "# h L1 L2 Loo E[X] Var[X] Min[X] Max[X] Nb[X]" << std::endl;
368  adev_output << estimator.h()
369  << " " << adev_stat.mean() // L1
370  << " " << sqrt( adev_stat.unbiasedVariance()
371  + adev_stat.mean()*adev_stat.mean() ) // L2
372  << " " << adev_stat.max() // Loo
373  << " " << adev_stat.mean() // E[X] (=L1)
374  << " " << adev_stat.unbiasedVariance() // Var[X]
375  << " " << adev_stat.min() // Min[X]
376  << " " << adev_stat.max() // Max[X]
377  << " " << adev_stat.samples() // Nb[X]
378  << std::endl;
379  adev_output.close();
380  trace.endBlock();
381  }
382  if ( params.exportX != "None" )
383  {
384  trace.beginBlock( "Exporting cell geometry." );
385  std::ostringstream export_sstr;
386  export_sstr << fname << "-" << nameEstimator << "-cells-"
387  << estimator.h() << ".txt";
388  std::ofstream export_output( export_sstr.str().c_str() );
389  export_output << "# ImaGene viewer (viewSetOfSurfels) file format for displaying cells." << std::endl;
390  bool adev = params.exportX == "AngleDeviation";
391  unsigned int i = 0;
392  range = CountedPtr<VisitorRange>( new VisitorRange( new Visitor( surface, *(surface.begin()) )) );
393  for ( VisitorConstIterator it = range->begin(), itE = range->end(); it != itE; ++it, ++i )
394  {
395  Quantity n_est = n_estimations[ i ];
396  Quantity n_true_est = n_true_estimations[ i ];
397  Scalar angle_error = acos( n_est.dot( n_true_est ) )*180.0 / 3.14159625;
398  Surfel s = *it;
399  export_output
400  << "CellN"
401  << " " << min( 1023, max( 512+K.sKCoord( s, 0 ), 0 ) )
402  << " " << min( 1023, max( 512+K.sKCoord( s, 1 ), 0 ) )
403  << " " << min( 1023, max( 512+K.sKCoord( s, 2 ), 0 ) )
404  << " " << K.sSign( s );
405  Color c = grad( 0 );
406  if ( adev ) c = grad( max( 0.0, min( angle_error, 40.0 ) ) );
407  export_output << " " << ((double) c.red() / 255.0 )
408  << " " << ((double) c.green() / 255.0 )
409  << " " << ((double) c.blue() / 255.0 );
410  export_output << " " << n_est[ 0 ] << " " << n_est[ 1 ]
411  << " " << n_est[ 2 ] << std::endl;
412  }
413  export_output.close();
414  trace.endBlock();
415  }
416  if ( params.normals )
417  {
418  trace.beginBlock( "Exporting cells normals." );
419  std::ostringstream export_sstr;
420  export_sstr << fname << "-" << nameEstimator << "-normals-"
421  << estimator.h() << ".txt";
422  std::ofstream export_output( export_sstr.str().c_str() );
423  export_output << "# kx ky kz sign n_est[0] n_est[1] n_est[2] n_true[0] n_true[1] n_true[2]" << std::endl;
424  unsigned int i = 0;
425  range = CountedPtr<VisitorRange>( new VisitorRange( new Visitor( surface, *(surface.begin()) )) );
426  for ( VisitorConstIterator it = range->begin(), itE = range->end(); it != itE; ++it, ++i )
427  {
428  Quantity n_est = n_estimations[ i ];
429  Quantity n_true_est = n_true_estimations[ i ];
430  Surfel s = *it;
431  export_output
432  << K.sKCoord( s, 0 ) << " " << K.sKCoord( s, 1 ) << " " << K.sKCoord( s, 2 )
433  << " " << K.sSign( s )
434  << " " << n_est[ 0 ] << " " << n_est[ 1 ] << " " << n_est[ 2 ]
435  << " " << n_true_est[ 0 ] << " " << n_true_est[ 1 ] << " " << n_true_est[ 2 ]
436  << std::endl;
437  }
438  export_output.close();
439  trace.endBlock();
440  }
441  if ( params.noff )
442  {
443  trace.beginBlock( "Exporting NOFF file." );
444  std::ostringstream export_sstr;
445  export_sstr << fname << "-" << nameEstimator << "-noff-"
446  << estimator.h() << ".off";
447  std::ofstream export_output( export_sstr.str().c_str() );
448  std::map<Surfel,Quantity> normals;
449  unsigned int i = 0;
450  range = CountedPtr<VisitorRange>( new VisitorRange( new Visitor( surface, *(surface.begin()) )) );
451  for ( VisitorConstIterator it = range->begin(), itE = range->end(); it != itE; ++it, ++i )
452  {
453  Quantity n_est = n_estimations[ i ];
454  normals[ *it ] = n_est;
455  }
456  CanonicSCellEmbedder<KSpace> surfelEmbedder( K );
457  typedef SCellEmbedderWithNormal< CanonicSCellEmbedder<KSpace> > Embedder;
458  Embedder embedder( surfelEmbedder, normals );
459  surface.exportAs3DNOFF( export_output, embedder );
460  export_output.close();
461  trace.endBlock();
462  }
463 #ifdef WITH_VISU3D_QGLVIEWER
464  if ( params.exportX != "None" )
465  {
466  typedef typename KSpace::Space Space;
467  typedef Viewer3D<Space,KSpace> MyViewever3D;
468  typedef Display3DFactory<Space,KSpace> MyDisplay3DFactory;
469  int argc = 1;
470  char name[] = "Viewer";
471  char* argv[ 1 ];
472  argv[ 0 ] = name;
473  Surfel s;
474  QApplication application( argc, argv );
475  MyViewever3D viewer( K );
476  viewer.show();
477  viewer << SetMode3D( s.className(), "Basic" );
478  trace.beginBlock( "Viewing surface." );
479  bool adev = params.view == "AngleDeviation";
480 
481  unsigned int i = 0;
482  range = CountedPtr<VisitorRange>( new VisitorRange( new Visitor( surface, *(surface.begin()) )) );
483  for ( VisitorConstIterator it = range->begin(), itE = range->end(); it != itE; ++it, ++i )
484  {
485  Quantity n_est = n_estimations[ i ];
486  Quantity n_true_est = n_true_estimations[ i ];
487  Scalar angle_error = acos( n_est.dot( n_true_est ) )*180.0 / 3.14159625;
488  s = *it;
489  Color c = grad( 0 );
490  if ( adev ) c = grad( max( 0.0, min( angle_error, 40.0 ) ) );
491  viewer.setFillColor( c );
492  MyDisplay3DFactory::drawOrientedSurfelWithNormal( viewer, s, n_est, false );
493  }
494  trace.endBlock();
495  viewer << MyViewever3D::updateDisplay;
496  application.exec();
497  }
498 #endif
499 
500 }
501 
502 template <typename KSpace,
503  typename ImplicitShape,
504  typename Surface,
505  typename KernelFunction,
506  typename PointPredicate>
507 void chooseEstimator
508 ( const AllParams &params, //< command-line parameters
509  const KSpace& K, //< cellular grid space
510  const ImplicitShape& shape, //< implicit shape "ground truth"
511  const Surface& surface, //< digital surface approximating shape
512  const KernelFunction& chi, //< the kernel function
513  const PointPredicate& ptPred ) //< analysed implicit digital shape as a PointPredicate
514 {
515  string nameEstimator = params.estimator;
516  double h = params.gridstep;
517  typedef functors::ShapeGeometricFunctors::ShapeNormalVectorFunctor<ImplicitShape> NormalFunctor;
518  typedef TrueDigitalSurfaceLocalEstimator<KSpace, ImplicitShape, NormalFunctor> TrueEstimator;
519  TrueEstimator true_estimator;
520  true_estimator.attach( shape );
521  true_estimator.setParams( K, NormalFunctor(), 20, 0.1, 0.01 );
522  true_estimator.init( h, surface.begin(), surface.end() );
523  if ( nameEstimator == "True" )
524  {
525  trace.beginBlock( "Chosen estimator is: True." );
526  typedef TrueDigitalSurfaceLocalEstimator<KSpace, ImplicitShape, NormalFunctor> Estimator;
527  Estimator estimator;
528  estimator.attach( shape );
529  estimator.setParams( K, NormalFunctor(), 20, 0.1, 0.01 );
530  estimator.init( h, surface.begin(), surface.end() );
531  trace.endBlock();
532  computeEstimation( params, K, shape, surface, true_estimator, estimator );
533  }
534  else if ( nameEstimator == "VCM" )
535  {
536  trace.beginBlock( "Chosen estimator is: VCM." );
537  typedef typename KSpace::Space Space;
538  typedef typename Surface::DigitalSurfaceContainer SurfaceContainer;
539  typedef ExactPredicateLpSeparableMetric<Space,2> Metric;
540  typedef VoronoiCovarianceMeasureOnDigitalSurface<SurfaceContainer,Metric,
541  KernelFunction> VCMOnSurface;
542  typedef functors::VCMNormalVectorFunctor<VCMOnSurface> NormalFunctor;
543  typedef VCMDigitalSurfaceLocalEstimator<SurfaceContainer,Metric,
544  KernelFunction, NormalFunctor> VCMNormalEstimator;
545  int embedding = params.embedding;
546  Surfel2PointEmbedding embType = embedding == 0 ? Pointels :
547  embedding == 1 ? InnerSpel : OuterSpel;
548  double R = params.R_radius;
549  double r = params.r_radius;
550  double t = params.trivial_radius;
551  double alpha = params.alpha;
552  if ( alpha != 0.0 ) R *= pow( h, alpha-1.0 );
553  if ( alpha != 0.0 ) r *= pow( h, alpha-1.0 );
554  trace.info() << "- R=" << R << " r=" << r << " t=" << t << std::endl;
555  VCMNormalEstimator estimator;
556  estimator.attach( surface );
557  estimator.setParams( embType, R, r, chi, t, Metric(), true );
558  estimator.init( h, surface.begin(), surface.end() );
559  trace.endBlock();
560  computeEstimation( params, K, shape, surface, true_estimator, estimator );
561  }
562  else if ( nameEstimator == "II" )
563  {
564  trace.beginBlock( "Chosen estimator is: II." );
565  typedef typename KSpace::Space Space;
566  typedef HyperRectDomain<Space> Domain;
567  typedef ImageContainerBySTLVector< Domain, bool> Image;
568  typedef typename Domain::ConstIterator DomainConstIterator;
569  typedef functors::SimpleThresholdForegroundPredicate<Image> ThresholdedImage;
570  typedef functors::IINormalDirectionFunctor<Space> IINormalFunctor;
571  typedef IntegralInvariantCovarianceEstimator<KSpace, ThresholdedImage, IINormalFunctor> IINormalEstimator;
572  double r = params.r_radius;
573  double alpha = params.alpha;
574  if ( alpha != 0.0 ) r *= pow( h, alpha-1.0 );
575  trace.info() << " r=" << r << std::endl;
576  trace.beginBlock( "Preparing characteristic set." );
577  Domain domain( K.lowerBound(), K.upperBound() );
578  Image image( domain );
579  for ( DomainConstIterator it = domain.begin(), itE = domain.end(); it != itE; ++it )
580  {
581  image.setValue( *it, ptPred( *it ) );
582  }
583  trace.endBlock();
584  trace.beginBlock( "Initialize II estimator." );
585  ThresholdedImage thresholdedImage( image, false );
586  IINormalEstimator ii_estimator( K, thresholdedImage );
587  ii_estimator.setParams( r );
588  ii_estimator.init( h, surface.begin(), surface.end() );
589  trace.endBlock();
590  trace.endBlock();
591  computeEstimation( params, K, shape, surface, true_estimator, ii_estimator );
592  }
593 
594 }
595 
596 template <typename KSpace,
597  typename ImplicitShape,
598  typename Surface,
599  typename PointPredicate>
600 void chooseKernel
601 ( const AllParams& params, //< command-line parameters
602  const KSpace& K, //< cellular grid space
603  const ImplicitShape& shape, //< implicit shape "ground truth"
604  const Surface& surface, //< digital surface approximating shape
605  const PointPredicate& ptPred ) //< analysed implicit digital shape as a PointPredicate
606 {
607  string kernel = params.kernel;
608  double h = params.gridstep;
609  double r = params.r_radius;
610  double alpha = params.alpha;
611  if ( alpha != 0.0 ) r *= pow( h, alpha-1.0 );
612  if ( kernel == "hat" ) {
613  typedef typename KSpace::Point Point;
614  typedef functors::HatPointFunction<Point,double> KernelFunction;
615  KernelFunction chi_r( 1.0, r );
616  chooseEstimator( params, K, shape, surface, chi_r, ptPred );
617  } else if ( kernel == "ball" ) {
618  typedef typename KSpace::Point Point;
619  typedef functors::BallConstantPointFunction<Point,double> KernelFunction;
620  KernelFunction chi_r( 1.0, r );
621  chooseEstimator( params, K, shape, surface, chi_r, ptPred );
622  }
623 }
624 
625 template <typename KSpace,
626  typename ImplicitShape,
627  typename ImplicitDigitalShape >
628 int chooseSurface
629 ( const AllParams &params, //< command-line parameters
630  const KSpace& K, //< cellular grid space
631  const ImplicitShape& shape, //< implicit shape "ground truth"
632  const ImplicitDigitalShape& dshape ) //< analysed implicit digital shape
633 {
634  // Selecting a model of surface depending on noise / not noise.
635  typedef double Scalar;
636  if ( params.noise == 0.0 )
637  { // no noise
638  trace.beginBlock( "Make digital surface..." );
639  typedef LightImplicitDigitalSurface<KSpace,ImplicitDigitalShape> SurfaceContainer;
640  typedef DigitalSurface< SurfaceContainer > Surface;
641  typedef typename Surface::Surfel Surfel;
642  SurfelAdjacency< KSpace::dimension > surfAdj( true );
643  Surfel bel;
644  try {
645  bel = Surfaces<KSpace>::findABel( K, dshape, 10000 );
646  } catch (DGtal::InputException e) {
647  trace.error() << "ERROR Unable to find bel." << std::endl;
648  return 3;
649  }
650  SurfaceContainer* surfaceContainer = new SurfaceContainer( K, dshape, surfAdj, bel );
651  CountedPtr<Surface> ptrSurface( new Surface( surfaceContainer ) ); // acquired
652  trace.info() << "- surface component has " << ptrSurface->size() << " surfels." << std::endl;
653  trace.endBlock();
654  chooseKernel( params, K, shape, *ptrSurface, dshape );
655  }
656  else
657  { // noise
658  trace.beginBlock( "Make digital surface..." );
659  typedef typename ImplicitDigitalShape::Domain Domain;
660  typedef KanungoNoise< ImplicitDigitalShape, Domain > KanungoPredicate;
661  typedef LightImplicitDigitalSurface< KSpace, KanungoPredicate > SurfaceContainer;
662  typedef DigitalSurface< SurfaceContainer > Surface;
663  typedef typename Surface::Surfel Surfel;
664  SurfelAdjacency< KSpace::dimension > surfAdj( true );
665  Surfel bel;
666  const Domain shapeDomain = dshape.getDomain();
667  KanungoPredicate* noisified_dshape = new KanungoPredicate( dshape, shapeDomain, params.noise );
668  // We have to search for a big connected component.
669  CountedPtr<Surface> ptrSurface;
670  double minsize = dshape.getUpperBound()[0] - dshape.getLowerBound()[0];
671  unsigned int nb_surfels = 0;
672  unsigned int tries = 0;
673  do {
674  try { // Search initial bel
675  bel = Surfaces<KSpace>::findABel( K, *noisified_dshape, 10000 );
676  } catch (DGtal::InputException e) {
677  trace.error() << "ERROR Unable to find bel." << std::endl;
678  return 3;
679  }
680  SurfaceContainer* surfaceContainer = new SurfaceContainer( K, *noisified_dshape, surfAdj, bel );
681  ptrSurface = CountedPtr<Surface>( new Surface( surfaceContainer ) ); // acquired
682  nb_surfels = ptrSurface->size();
683  } while ( ( nb_surfels < 2 * minsize ) && ( tries++ < 150 ) );
684  if( tries >= 150 )
685  {
686  trace.error() << "ERROR cannot find a proper bel in a big enough component." << std::endl;
687  return 4;
688  }
689  trace.info() << "- surface component has " << nb_surfels << " surfels." << std::endl;
690  trace.endBlock();
691  chooseKernel( params, K, shape, *ptrSurface, *noisified_dshape );
692  }
693  return 0;
694 }
695 
696 
698 int main( int argc, char** argv )
699 {
700 
701  AllParams allParams;
702 
703  // parse command line ----------------------------------------------
704  CLI::App app;
705  std::stringstream ss_descr;
706  ss_descr << "Computes a normal vector field over a digitized 3D implicit surface for several estimators (II|VCM|Trivial|True), specified with -e. You may add Kanungo noise with option -N. These estimators are compared with ground truth. You may then: 1) visualize the normals or the angle deviations with -V (if WITH_QGL_VIEWER is enabled), 2) outputs them as a list of cells/estimations with -n, 3) outputs them as a ImaGene file with -O, 4) outputs them as a NOFF file with -O, 5) computes estimation statistics with option -S.";
707 
708  ss_descr << "Example of use:\n"
709  << "./generic3dNormalEstimator -p \"90-3*x^2-2*y^2-z^2\" -o VCM-ellipse -a -10 -A 10 -e VCM -R 3 -r 3 -t 2 -E 0 -x Normals" << endl << endl
710  << "Example of implicit surfaces (specified by -p):" << endl
711  << " - ellipse : 90-3*x^2-2*y^2-z^2 " << endl
712  << " - torus : -1*(x^2+y^2+z^2+6*6-2*2)^2+4*6*6*(x^2+y^2) " << endl
713  << " - rcube : 6561-x^4-y^4-z^4" << endl
714  << " - goursat : 8-0.03*x^4-0.03*y^4-0.03*z^4+2*x^2+2*y^2+2*z^2" << endl
715  << " - distel : 10000-(x^2+y^2+z^2+1000*(x^2+y^2)*(x^2+z^2)*(y^2+z^2))" << endl
716  << " - leopold : 100-(x^2*y^2*z^2+4*x^2+4*y^2+3*z^2)" << endl
717  << " - diabolo : x^2-(y^2+z^2)^2" << endl
718  << " - heart : -1*(x^2+2.25*y^2+z^2-1)^3+x^2*z^3+0.1125*y^2*z^3" << endl
719  << " - crixxi : -0.9*(y^2+z^2-1)^2-(x^2+y^2-1)^3" << endl << endl;
720  ss_descr << "Implemented estimators (specified by -e):" << endl
721  << " - True : supposed to be the ground truth for computations. Of course, it is only approximations." << endl
722  << " - VCM : normal estimator by digital Voronoi Covariance Matrix. Radii parameters are given by -R, -r." << endl
723  << " - II : normal estimator by Integral Invariants. Radius parameter is given by -r." << endl
724  << " - Trivial : the normal obtained by average trivial surfel normals in a ball neighborhood. Radius parameter is given by -r." << endl
725  << endl;
726  ss_descr << "Note:" << endl
727  << " - This is a normal *direction* evaluator more than a normal vector evaluator. Orientations of normals are deduced from ground truth. This is due to the fact that II and VCM only estimates normal directions." << endl
728  << " - This tool only analyses one surface component, and one that contains at least as many surfels as the width of the digital bounding box. This is required when analysing noisy data, where a lot of the small components are spurious. The drawback is that you cannot analyse the normals on a surface with several components." << endl;
729 
730  app.description(ss_descr.str());
731 
732  app.add_option("--polynomial,-p",allParams.polynomials, "the implicit polynomial whose zero-level defines the shape of interest.") ->required();
733  app.add_option("--noise,-N", allParams.noise,"the Kanungo noise level l=arg, with l^d the probability that a point at distance d is flipped inside/outside.", true );
734  app.add_option("--minAABB,-a", allParams.minAABB, "the min value of the AABB bounding box (domain)", true );
735  app.add_option("--maxAABB,-A", allParams.maxAABB, "the max value of the AABB bounding box (domain)", true );
736  app.add_option("--gridstep,-g", allParams.gridstep, "the gridstep that defines the digitization (often called h).", true);
737  app.add_option("--estimator,-e", allParams.estimator, "the chosen normal estimator: True | VCM | II | Trivial", true)
738  -> check(CLI::IsMember({"True", "VCM", "II", "Trivial"}));
739  app.add_option("--R-radius,-R", allParams.R_radius, "the constant for parameter R in R(h)=R h^alpha (VCM).",true);
740  app.add_option("--r-radius,-r", allParams.r_radius, "the constant for parameter r in r(h)=r h^alpha (VCM,II,Trivial).", true);
741  app.add_option("--kernel,-k", allParams.kernel, "the function chi_r, either hat or ball.", true);
742  app.add_option("--alpha", allParams.alpha, "the parameter alpha in r(h)=r h^alpha (VCM).", true);
743  app.add_option("--trivial-radius,-t", allParams.trivial_radius, "the parameter t defining the radius for the Trivial estimator. Also used for reorienting the VCM.", true);
744  app.add_option("--embedding,-E",allParams.embedding, "the surfel -> point embedding for VCM estimator: 0: Pointels, 1: InnerSpel, 2: OuterSpel.", true);
745  app.add_option("--output,-o", allParams.output, "the output basename. All generated files will have the form <arg>-*, for instance <arg>-angle-deviation-<gridstep>.txt, <arg>-normals-<gridstep>.txt, <arg>-cells-<gridstep>.txt, <arg>-noff-<gridstep>.off.", true);
746  app.add_flag("--angle-deviation-stats,-S", allParams.angle_deviation_stats, "computes angle deviation error and outputs them in file <basename>-angle-deviation-<gridstep>.txt, as specified by -o <basename>.");
747 
748  app.add_option("--export,-x",allParams.exportX, "exports surfel normals which can be viewed with ImaGene tool 'viewSetOfSurfels' in file <basename>-cells-<gridstep>.txt, as specified by -o <basename>. Parameter <arg> is None|Normals|AngleDeviation. The color depends on the angle deviation in degree: 0 metallic blue, 5 light cyan, 10 light green, 15 light yellow, 20 yellow, 25 orange, 30 red, 35, dark red, 40- grey", true );
749  app.add_flag("--normals,-n", allParams.normals, "outputs every surfel, its estimated normal, and the ground truth normal in file <basename>-normals-<gridstep>.txt, as specified by -o <basename>.");
750  app.add_flag("--noff,-O", allParams.noff, "exports the digital surface with normals as NOFF file <basename>-noff-<gridstep>.off, as specified by -o <basename>..");
751 #ifdef WITH_VISU3D_QGLVIEWER
752  app.add_option("--view,-V", allParams.view, "view the digital surface with normals. Parameter <arg> is None|Normals|AngleDeviation. The color depends on the angle deviation in degree: 0 metallic blue, 5 light cyan, 10 light green, 15 light yellow, 20 yellow, 25 orange, 30 red, 35, dark red, 40- grey.");
753 #endif
754 
755 
756  app.get_formatter()->column_width(40);
757  CLI11_PARSE(app, argc, argv);
758  // END parse command line using CLI ----------------------------------------------
759 
760 
761 
762  trace.beginBlock( "Make implicit shape..." );
763  typedef Z3i::Space Space;
764  typedef double Scalar;
765  typedef MPolynomial< 3, Scalar > Polynomial3;
766  typedef MPolynomialReader<3, Scalar> Polynomial3Reader;
767  typedef ImplicitPolynomial3Shape<Space> ImplicitShape;
768  Polynomial3 poly;
769  Polynomial3Reader reader;
770  string::const_iterator iter = reader.read( poly, allParams.polynomials.begin(), allParams.polynomials.end() );
771  if ( iter != allParams.polynomials.end() )
772  {
773  trace.error() << "ERROR reading polynomial: I read only <" << allParams.polynomials.substr( 0, iter - allParams.polynomials.begin() )
774  << ">, and I built P=" << poly << std::endl;
775  return 2;
776  }
777  CountedPtr<ImplicitShape> shape( new ImplicitShape( poly ) ); // smart pointer
778  trace.endBlock();
779 
780  trace.beginBlock( "Make implicit digital shape..." );
781  typedef Z3i::KSpace KSpace;
782  typedef Space::RealPoint RealPoint;
783  typedef GaussDigitizer< Space, ImplicitShape > ImplicitDigitalShape;
784  typedef ImplicitDigitalShape::Domain Domain;
785  Scalar min_x = allParams.minAABB;
786  Scalar max_x = allParams.maxAABB;
787  Scalar h = allParams.gridstep;
788  RealPoint p1( min_x, min_x, min_x );
789  RealPoint p2( max_x, max_x, max_x );
790  CountedPtr<ImplicitDigitalShape> dshape( new ImplicitDigitalShape() );
791  dshape->attach( *shape );
792  dshape->init( p1, p2, h );
793  Domain domain = dshape->getDomain();
794  KSpace K;
795  K.init( domain.lowerBound(), domain.upperBound(), true );
796  trace.info() << "- domain is " << domain << std::endl;
797  trace.endBlock();
798 
799  chooseSurface( allParams, K, *shape, *dshape );
800 
801  return 0;
802 }
803 
Definition: ATu0v1.h:57