35 #include <boost/program_options/options_description.hpp> 36 #include <boost/program_options/parsers.hpp> 37 #include <boost/program_options/variables_map.hpp> 39 #include "DGtal/base/Common.h" 40 #include "DGtal/base/CountedPtr.h" 41 #include "DGtal/base/CountedConstPtrOrConstPtr.h" 42 #include "DGtal/helpers/StdDefs.h" 43 #include "DGtal/math/Statistic.h" 44 #include "DGtal/math/MPolynomial.h" 45 #include "DGtal/topology/LightImplicitDigitalSurface.h" 46 #include "DGtal/graph/DepthFirstVisitor.h" 47 #include "DGtal/graph/GraphVisitorRange.h" 48 #include "DGtal/geometry/surfaces/estimation/CNormalVectorEstimator.h" 49 #include "DGtal/geometry/surfaces/estimation/VoronoiCovarianceMeasureOnDigitalSurface.h" 50 #include "DGtal/geometry/surfaces/estimation/VCMDigitalSurfaceLocalEstimator.h" 51 #include "DGtal/geometry/surfaces/estimation/TrueDigitalSurfaceLocalEstimator.h" 52 #include "DGtal/geometry/surfaces/estimation/IIGeometricFunctors.h" 53 #include "DGtal/geometry/surfaces/estimation/IntegralInvariantCovarianceEstimator.h" 54 #include "DGtal/geometry/volumes/KanungoNoise.h" 55 #include "DGtal/shapes/GaussDigitizer.h" 56 #include "DGtal/shapes/ShapeGeometricFunctors.h" 57 #include "DGtal/shapes/implicit/ImplicitPolynomial3Shape.h" 58 #include "DGtal/images/ImageContainerBySTLVector.h" 59 #include "DGtal/images/SimpleThresholdForegroundPredicate.h" 60 #include "DGtal/io/readers/MPolynomialReader.h" 61 #include "DGtal/io/colormaps/GradientColorMap.h" 62 #ifdef WITH_VISU3D_QGLVIEWER 63 #include "DGtal/io/viewers/Viewer3D.h" 64 #include "DGtal/io/Display3DFactory.h" 68 using namespace DGtal;
69 namespace po = boost::program_options;
240 template <
typename SCell,
typename RealVector>
241 struct GradientMapAdapter {
242 typedef std::map<SCell,RealVector> SCell2RealVectorMap;
243 typedef SCell Argument;
244 typedef RealVector Value;
245 GradientMapAdapter( ConstAlias<SCell2RealVectorMap> map )
247 RealVector operator()(
const Argument& arg )
const 249 typename SCell2RealVectorMap::const_iterator it = myMap->find( arg );
250 if ( it != myMap->end() )
return it->second;
251 else return RealVector();
253 CountedConstPtrOrConstPtr<SCell2RealVectorMap> myMap;
256 template <
typename SCellEmbedder>
257 struct SCellEmbedderWithNormal :
public SCellEmbedder
259 using SCellEmbedder::space;
260 using SCellEmbedder::operator();
261 typedef typename SCellEmbedder::KSpace KSpace;
262 typedef typename SCellEmbedder::SCell SCell;
263 typedef typename SCellEmbedder::RealPoint RealPoint;
264 typedef SCell Argument;
265 typedef RealPoint Value;
266 typedef typename KSpace::Space::RealVector RealVector;
267 typedef std::map<SCell,RealVector> SCell2RealVectorMap;
268 typedef GradientMapAdapter<SCell,RealVector> GradientMap;
270 SCellEmbedderWithNormal( ConstAlias<SCellEmbedder> embedder,
271 ConstAlias<SCell2RealVectorMap> map )
272 : SCellEmbedder( embedder ), myMap( map )
275 GradientMap gradientMap()
const 277 return GradientMap( myMap );
280 CountedConstPtrOrConstPtr<SCell2RealVectorMap> myMap;
283 template <
typename DigitalSurface,
285 void exportNOFFSurface(
const DigitalSurface& surface,
286 const Estimator& estimator,
287 std::ostream& output )
289 typedef typename DigitalSurface::KSpace KSpace;
290 typedef typename DigitalSurface::ConstIterator ConstIterator;
291 typedef typename DigitalSurface::Surfel Surfel;
292 typedef typename KSpace::SCell SCell;
293 typedef typename Estimator::Quantity Quantity;
294 const KSpace& ks = surface.container().space();
295 std::map<Surfel,Quantity> normals;
296 for ( ConstIterator it = surface.begin(), itE = surface.end(); it != itE; ++it )
298 Quantity n_est = estimator.eval( it );
299 normals[ *it ] = n_est;
301 CanonicSCellEmbedder<KSpace> surfelEmbedder( ks );
302 typedef SCellEmbedderWithNormal< CanonicSCellEmbedder<KSpace> > Embedder;
303 Embedder embedder( surfelEmbedder, normals );
304 surface.exportAs3DNOFF( output, embedder );
311 template <
typename KSpace,
312 typename ImplicitShape,
314 typename TrueEstimator,
316 void computeEstimation
317 (
const po::variables_map& vm,
319 const ImplicitShape& shape,
320 const Surface& surface,
321 TrueEstimator& true_estimator,
322 Estimator& estimator )
324 typedef typename Surface::ConstIterator ConstIterator;
325 typedef typename Surface::Surfel Surfel;
326 typedef typename Estimator::Quantity Quantity;
327 typedef double Scalar;
328 typedef DepthFirstVisitor< Surface > Visitor;
329 typedef GraphVisitorRange< Visitor > VisitorRange;
330 typedef typename VisitorRange::ConstIterator VisitorConstIterator;
332 std::string fname = vm[
"output" ].as<std::string>();
333 string nameEstimator = vm[
"estimator" ].as<
string>();
334 trace.beginBlock(
"Computing " + nameEstimator +
"estimations." );
335 CountedPtr<VisitorRange> range(
new VisitorRange(
new Visitor( surface, *(surface.begin()) )) );
336 std::vector<Quantity> n_estimations;
337 estimator.eval( range->begin(), range->end(), std::back_inserter( n_estimations ) );
338 trace.info() <<
"- nb estimations = " << n_estimations.size() << std::endl;
341 trace.beginBlock(
"Computing ground truth." );
342 range = CountedPtr<VisitorRange>(
new VisitorRange(
new Visitor( surface, *(surface.begin()) )) );
343 std::vector<Quantity> n_true_estimations;
344 true_estimator.eval( range->begin(), range->end(), std::back_inserter( n_true_estimations ) );
345 trace.info() <<
"- nb estimations = " << n_true_estimations.size() << std::endl;
348 trace.beginBlock(
"Correcting orientations." );
349 ASSERT( n_estimations.size() == n_true_estimations.size() );
350 for (
unsigned int i = 0; i < n_estimations.size(); ++i )
351 if ( n_estimations[ i ].dot( n_true_estimations[ i ] ) < 0 )
352 n_estimations[ i ] = -n_estimations[ i ];
355 DGtal::GradientColorMap<double> grad( 0.0, 40.0 );
358 grad.addColor( DGtal::Color( 128, 128, 255 ) );
359 grad.addColor( DGtal::Color( 128, 255, 255 ) );
360 grad.addColor( DGtal::Color( 128, 255, 128 ) );
361 grad.addColor( DGtal::Color( 255, 255, 128 ) );
362 grad.addColor( DGtal::Color( 255, 255, 0 ) );
363 grad.addColor( DGtal::Color( 255, 128, 0 ) );
364 grad.addColor( DGtal::Color( 255, 0, 0 ) );
365 grad.addColor( DGtal::Color( 128, 0, 0 ) );
366 grad.addColor( DGtal::Color( 128, 128, 128 ) );
368 if ( vm.count(
"angle-deviation-stats" ) )
370 trace.beginBlock(
"Computing angle deviation error stats." );
371 std::ostringstream adev_sstr;
372 adev_sstr << fname <<
"-" << nameEstimator <<
"-angle-deviation-" 373 << estimator.h() <<
".txt";
374 DGtal::Statistic<Scalar> adev_stat;
376 range = CountedPtr<VisitorRange>(
new VisitorRange(
new Visitor( surface, *(surface.begin()) )) );
377 for ( VisitorConstIterator it = range->begin(), itE = range->end(); it != itE; ++it, ++i )
379 Quantity n_est = n_estimations[ i ];
380 Quantity n_true_est = n_true_estimations[ i ];
381 Scalar angle_error = acos( n_est.dot( n_true_est ) );
382 adev_stat.addValue( angle_error );
384 adev_stat.terminate();
385 std::ofstream adev_output( adev_sstr.str().c_str() );
386 adev_output <<
"# Average error X of the absolute angle between two vector estimations." << std::endl;
387 adev_output <<
"# h L1 L2 Loo E[X] Var[X] Min[X] Max[X] Nb[X]" << std::endl;
388 adev_output << estimator.h()
389 <<
" " << adev_stat.mean()
390 <<
" " << sqrt( adev_stat.unbiasedVariance()
391 + adev_stat.mean()*adev_stat.mean() )
392 <<
" " << adev_stat.max()
393 <<
" " << adev_stat.mean()
394 <<
" " << adev_stat.unbiasedVariance()
395 <<
" " << adev_stat.min()
396 <<
" " << adev_stat.max()
397 <<
" " << adev_stat.samples()
402 if ( vm[
"export" ].as<string>() !=
"None" )
404 trace.beginBlock(
"Exporting cell geometry." );
405 std::ostringstream export_sstr;
406 export_sstr << fname <<
"-" << nameEstimator <<
"-cells-" 407 << estimator.h() <<
".txt";
408 std::ofstream export_output( export_sstr.str().c_str() );
409 export_output <<
"# ImaGene viewer (viewSetOfSurfels) file format for displaying cells." << std::endl;
410 bool adev = vm[
"export" ].as<
string>() ==
"AngleDeviation";
412 range = CountedPtr<VisitorRange>(
new VisitorRange(
new Visitor( surface, *(surface.begin()) )) );
413 for ( VisitorConstIterator it = range->begin(), itE = range->end(); it != itE; ++it, ++i )
415 Quantity n_est = n_estimations[ i ];
416 Quantity n_true_est = n_true_estimations[ i ];
417 Scalar angle_error = acos( n_est.dot( n_true_est ) )*180.0 / 3.14159625;
421 <<
" " << min( 1023, max( 512+K.sKCoord( s, 0 ), 0 ) )
422 <<
" " << min( 1023, max( 512+K.sKCoord( s, 1 ), 0 ) )
423 <<
" " << min( 1023, max( 512+K.sKCoord( s, 2 ), 0 ) )
424 <<
" " << K.sSign( s );
426 if ( adev ) c = grad( max( 0.0, min( angle_error, 40.0 ) ) );
427 export_output <<
" " << ((double) c.red() / 255.0 )
428 <<
" " << ((
double) c.green() / 255.0 )
429 << " " << ((
double) c.blue() / 255.0 );
430 export_output << " " << n_est[ 0 ] << " " << n_est[ 1 ]
431 << " " << n_est[ 2 ] <<
std::endl;
433 export_output.close();
436 if ( vm.count( "normals" ) )
438 trace.beginBlock(
"Exporting cells normals." );
439 std::ostringstream export_sstr;
440 export_sstr << fname <<
"-" << nameEstimator <<
"-normals-" 441 << estimator.h() <<
".txt";
442 std::ofstream export_output( export_sstr.str().c_str() );
443 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;
445 range = CountedPtr<VisitorRange>(
new VisitorRange(
new Visitor( surface, *(surface.begin()) )) );
446 for ( VisitorConstIterator it = range->begin(), itE = range->end(); it != itE; ++it, ++i )
448 Quantity n_est = n_estimations[ i ];
449 Quantity n_true_est = n_true_estimations[ i ];
452 << K.sKCoord( s, 0 ) <<
" " << K.sKCoord( s, 1 ) <<
" " << K.sKCoord( s, 2 )
453 <<
" " << K.sSign( s )
454 <<
" " << n_est[ 0 ] <<
" " << n_est[ 1 ] <<
" " << n_est[ 2 ]
455 <<
" " << n_true_est[ 0 ] <<
" " << n_true_est[ 1 ] <<
" " << n_true_est[ 2 ]
458 export_output.close();
461 if ( vm.count(
"noff" ) )
463 trace.beginBlock(
"Exporting NOFF file." );
464 std::ostringstream export_sstr;
465 export_sstr << fname <<
"-" << nameEstimator <<
"-noff-" 466 << estimator.h() <<
".off";
467 std::ofstream export_output( export_sstr.str().c_str() );
468 std::map<Surfel,Quantity> normals;
470 range = CountedPtr<VisitorRange>(
new VisitorRange(
new Visitor( surface, *(surface.begin()) )) );
471 for ( VisitorConstIterator it = range->begin(), itE = range->end(); it != itE; ++it, ++i )
473 Quantity n_est = n_estimations[ i ];
474 normals[ *it ] = n_est;
476 CanonicSCellEmbedder<KSpace> surfelEmbedder( K );
477 typedef SCellEmbedderWithNormal< CanonicSCellEmbedder<KSpace> > Embedder;
478 Embedder embedder( surfelEmbedder, normals );
479 surface.exportAs3DNOFF( export_output, embedder );
480 export_output.close();
483 #ifdef WITH_VISU3D_QGLVIEWER 484 if ( vm[
"view" ].as<string>() !=
"None" )
486 typedef typename KSpace::Space Space;
487 typedef Viewer3D<Space,KSpace> MyViewever3D;
488 typedef Display3DFactory<Space,KSpace> MyDisplay3DFactory;
490 char name[] =
"Viewer";
494 QApplication application( argc, argv );
495 MyViewever3D viewer( K );
497 viewer << SetMode3D( s.className(),
"Basic" );
498 trace.beginBlock(
"Viewing surface." );
499 bool adev = vm[
"view" ].as<
string>() ==
"AngleDeviation";
502 range = CountedPtr<VisitorRange>(
new VisitorRange(
new Visitor( surface, *(surface.begin()) )) );
503 for ( VisitorConstIterator it = range->begin(), itE = range->end(); it != itE; ++it, ++i )
505 Quantity n_est = n_estimations[ i ];
506 Quantity n_true_est = n_true_estimations[ i ];
507 Scalar angle_error = acos( n_est.dot( n_true_est ) )*180.0 / 3.14159625;
510 if ( adev ) c = grad( max( 0.0, min( angle_error, 40.0 ) ) );
511 viewer.setFillColor( c );
512 MyDisplay3DFactory::drawOrientedSurfelWithNormal( viewer, s, n_est,
false );
515 viewer << MyViewever3D::updateDisplay;
522 template <
typename KSpace,
523 typename ImplicitShape,
525 typename KernelFunction,
526 typename PointPredicate>
528 (
const po::variables_map& vm,
530 const ImplicitShape& shape,
531 const Surface& surface,
532 const KernelFunction& chi,
533 const PointPredicate& ptPred )
535 string nameEstimator = vm[
"estimator" ].as<
string>();
536 double h = vm[
"gridstep"].as<
double>();
537 typedef functors::ShapeGeometricFunctors::ShapeNormalVectorFunctor<ImplicitShape> NormalFunctor;
538 typedef TrueDigitalSurfaceLocalEstimator<KSpace, ImplicitShape, NormalFunctor> TrueEstimator;
539 TrueEstimator true_estimator;
540 true_estimator.attach( shape );
541 true_estimator.setParams( K, NormalFunctor(), 20, 0.1, 0.01 );
542 true_estimator.init( h, surface.begin(), surface.end() );
543 if ( nameEstimator ==
"True" )
545 trace.beginBlock(
"Chosen estimator is: True." );
546 typedef TrueDigitalSurfaceLocalEstimator<KSpace, ImplicitShape, NormalFunctor> Estimator;
548 estimator.attach( shape );
549 estimator.setParams( K, NormalFunctor(), 20, 0.1, 0.01 );
550 estimator.init( h, surface.begin(), surface.end() );
552 computeEstimation( vm, K, shape, surface, true_estimator, estimator );
554 else if ( nameEstimator ==
"VCM" )
556 trace.beginBlock(
"Chosen estimator is: VCM." );
557 typedef typename KSpace::Space Space;
558 typedef typename Surface::DigitalSurfaceContainer SurfaceContainer;
559 typedef ExactPredicateLpSeparableMetric<Space,2> Metric;
560 typedef VoronoiCovarianceMeasureOnDigitalSurface<SurfaceContainer,Metric,
561 KernelFunction> VCMOnSurface;
562 typedef functors::VCMNormalVectorFunctor<VCMOnSurface> NormalFunctor;
563 typedef VCMDigitalSurfaceLocalEstimator<SurfaceContainer,Metric,
564 KernelFunction, NormalFunctor> VCMNormalEstimator;
565 int embedding = vm[
"embedding"].as<
int>();
566 Surfel2PointEmbedding embType = embedding == 0 ? Pointels :
567 embedding == 1 ? InnerSpel : OuterSpel;
568 double R = vm[
"R-radius"].as<
double>();
569 double r = vm[
"r-radius"].as<
double>();
570 double t = vm[
"trivial-radius"].as<
double>();
571 double alpha = vm[
"alpha"].as<
double>();
572 if ( alpha != 0.0 ) R *= pow( h, alpha-1.0 );
573 if ( alpha != 0.0 ) r *= pow( h, alpha-1.0 );
574 trace.info() <<
"- R=" << R <<
" r=" << r <<
" t=" << t << std::endl;
575 VCMNormalEstimator estimator;
576 estimator.attach( surface );
577 estimator.setParams( embType, R, r, chi, t, Metric(),
true );
578 estimator.init( h, surface.begin(), surface.end() );
580 computeEstimation( vm, K, shape, surface, true_estimator, estimator );
582 else if ( nameEstimator ==
"II" )
584 trace.beginBlock(
"Chosen estimator is: II." );
585 typedef typename KSpace::Space Space;
586 typedef HyperRectDomain<Space> Domain;
587 typedef ImageContainerBySTLVector< Domain, bool> Image;
588 typedef typename Domain::ConstIterator DomainConstIterator;
589 typedef functors::SimpleThresholdForegroundPredicate<Image> ThresholdedImage;
590 typedef functors::IINormalDirectionFunctor<Space> IINormalFunctor;
591 typedef IntegralInvariantCovarianceEstimator<KSpace, ThresholdedImage, IINormalFunctor> IINormalEstimator;
592 double r = vm[
"r-radius"].as<
double>();
593 double alpha = vm[
"alpha"].as<
double>();
594 if ( alpha != 0.0 ) r *= pow( h, alpha-1.0 );
595 trace.info() <<
" r=" << r << std::endl;
596 trace.beginBlock(
"Preparing characteristic set." );
597 Domain domain( K.lowerBound(), K.upperBound() );
598 Image image( domain );
599 for ( DomainConstIterator it = domain.begin(), itE = domain.end(); it != itE; ++it )
601 image.setValue( *it, ptPred( *it ) );
604 trace.beginBlock(
"Initialize II estimator." );
605 ThresholdedImage thresholdedImage( image,
false );
606 IINormalEstimator ii_estimator( K, thresholdedImage );
607 ii_estimator.setParams( r );
608 ii_estimator.init( h, surface.begin(), surface.end() );
611 computeEstimation( vm, K, shape, surface, true_estimator, ii_estimator );
616 template <
typename KSpace,
617 typename ImplicitShape,
619 typename PointPredicate>
621 (
const po::variables_map& vm,
623 const ImplicitShape& shape,
624 const Surface& surface,
625 const PointPredicate& ptPred )
627 string kernel = vm[
"kernel" ].as<
string>();
628 double h = vm[
"gridstep"].as<
double>();
629 double r = vm[
"r-radius"].as<
double>();
630 double alpha = vm[
"alpha"].as<
double>();
631 if ( alpha != 0.0 ) r *= pow( h, alpha-1.0 );
632 if ( kernel ==
"hat" ) {
633 typedef typename KSpace::Point Point;
634 typedef functors::HatPointFunction<Point,double> KernelFunction;
635 KernelFunction chi_r( 1.0, r );
636 chooseEstimator( vm, K, shape, surface, chi_r, ptPred );
637 }
else if ( kernel ==
"ball" ) {
638 typedef typename KSpace::Point Point;
639 typedef functors::BallConstantPointFunction<Point,double> KernelFunction;
640 KernelFunction chi_r( 1.0, r );
641 chooseEstimator( vm, K, shape, surface, chi_r, ptPred );
645 template <
typename KSpace,
646 typename ImplicitShape,
647 typename ImplicitDigitalShape >
649 (
const po::variables_map& vm,
651 const ImplicitShape& shape,
652 const ImplicitDigitalShape& dshape )
655 typedef double Scalar;
656 Scalar noiseLevel = vm[
"noise" ].as<
double>();
657 if ( noiseLevel == 0.0 )
659 trace.beginBlock(
"Make digital surface..." );
660 typedef LightImplicitDigitalSurface<KSpace,ImplicitDigitalShape> SurfaceContainer;
661 typedef DigitalSurface< SurfaceContainer > Surface;
662 typedef typename Surface::Surfel Surfel;
663 SurfelAdjacency< KSpace::dimension > surfAdj(
true );
666 bel = Surfaces<KSpace>::findABel( K, dshape, 10000 );
667 }
catch (DGtal::InputException e) {
668 trace.error() <<
"ERROR Unable to find bel." << std::endl;
671 SurfaceContainer* surfaceContainer =
new SurfaceContainer( K, dshape, surfAdj, bel );
672 CountedPtr<Surface> ptrSurface(
new Surface( surfaceContainer ) );
673 trace.info() <<
"- surface component has " << ptrSurface->size() <<
" surfels." << std::endl;
675 chooseKernel( vm, K, shape, *ptrSurface, dshape );
679 trace.beginBlock(
"Make digital surface..." );
680 typedef typename ImplicitDigitalShape::Domain Domain;
681 typedef KanungoNoise< ImplicitDigitalShape, Domain > KanungoPredicate;
682 typedef LightImplicitDigitalSurface< KSpace, KanungoPredicate > SurfaceContainer;
683 typedef DigitalSurface< SurfaceContainer > Surface;
684 typedef typename Surface::Surfel Surfel;
685 SurfelAdjacency< KSpace::dimension > surfAdj(
true );
687 const Domain shapeDomain = dshape.getDomain();
688 KanungoPredicate* noisified_dshape =
new KanungoPredicate( dshape, shapeDomain, noiseLevel );
690 CountedPtr<Surface> ptrSurface;
691 double minsize = dshape.getUpperBound()[0] - dshape.getLowerBound()[0];
692 unsigned int nb_surfels = 0;
693 unsigned int tries = 0;
696 bel = Surfaces<KSpace>::findABel( K, *noisified_dshape, 10000 );
697 }
catch (DGtal::InputException e) {
698 trace.error() <<
"ERROR Unable to find bel." << std::endl;
701 SurfaceContainer* surfaceContainer =
new SurfaceContainer( K, *noisified_dshape, surfAdj, bel );
702 ptrSurface = CountedPtr<Surface>(
new Surface( surfaceContainer ) );
703 nb_surfels = ptrSurface->size();
704 }
while ( ( nb_surfels < 2 * minsize ) && ( tries++ < 150 ) );
707 trace.error() <<
"ERROR cannot find a proper bel in a big enough component." << std::endl;
710 trace.info() <<
"- surface component has " << nb_surfels <<
" surfels." << std::endl;
712 chooseKernel( vm, K, shape, *ptrSurface, *noisified_dshape );
719 int main(
int argc,
char** argv )
722 namespace po = boost::program_options;
723 po::options_description general_opt(
"Allowed options are: ");
724 general_opt.add_options()
725 (
"help,h",
"display this message")
726 (
"polynomial,p", po::value<string>(),
"the implicit polynomial whose zero-level defines the shape of interest." )
727 (
"noise,N", po::value<double>()->default_value( 0.0 ),
"the Kanungo noise level l=arg, with l^d the probability that a point at distance d is flipped inside/outside." )
728 (
"minAABB,a", po::value<double>()->default_value( -10.0 ),
"the min value of the AABB bounding box (domain)" )
729 (
"maxAABB,A", po::value<double>()->default_value( 10.0 ),
"the max value of the AABB bounding box (domain)" )
730 (
"gridstep,g", po::value< double >()->default_value( 1.0 ),
"the gridstep that defines the digitization (often called h). " )
731 (
"estimator,e", po::value<string>()->default_value(
"True" ),
"the chosen normal estimator: True | VCM | II | Trivial" )
732 (
"R-radius,R", po::value<double>()->default_value( 5 ),
"the constant for parameter R in R(h)=R h^alpha (VCM)." )
733 (
"r-radius,r", po::value<double>()->default_value( 3 ),
"the constant for parameter r in r(h)=r h^alpha (VCM,II,Trivial)." )
734 (
"kernel,k", po::value<string>()->default_value(
"hat" ),
"the function chi_r, either hat or ball." )
735 (
"alpha", po::value<double>()->default_value( 0.0 ),
"the parameter alpha in r(h)=r h^alpha (VCM)." )
736 (
"trivial-radius,t", po::value<double>()->default_value( 3 ),
"the parameter t defining the radius for the Trivial estimator. Also used for reorienting the VCM." )
737 (
"embedding,E", po::value<int>()->default_value( 0 ),
"the surfel -> point embedding for VCM estimator: 0: Pointels, 1: InnerSpel, 2: OuterSpel." )
738 (
"output,o", po::value<string>()->default_value(
"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." )
739 (
"angle-deviation-stats,S",
"computes angle deviation error and outputs them in file <basename>-angle-deviation-<gridstep>.txt, as specified by -o <basename>." )
740 (
"export,x", po::value<string>()->default_value(
"None" ),
"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" )
741 (
"normals,n",
"outputs every surfel, its estimated normal, and the ground truth normal in file <basename>-normals-<gridstep>.txt, as specified by -o <basename>." )
742 (
"noff,O",
"exports the digital surface with normals as NOFF file <basename>-noff-<gridstep>.off, as specified by -o <basename>.." )
743 #ifdef WITH_VISU3D_QGLVIEWER
744 (
"view,V", po::value<string>()->default_value(
"None" ),
"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." )
748 po::variables_map vm;
750 po::store(po::parse_command_line(argc, argv, general_opt), vm);
751 }
catch(
const exception& ex){
753 cerr <<
"Error checking program options: "<< ex.what()<< endl;
756 if( ! parseOK || vm.count(
"help") || ! vm.count(
"polynomial" ) )
758 if ( ! vm.count(
"polynomial" ) )
759 cerr <<
"Need parameter --polynomial / -p" << endl;
761 cerr <<
"Usage: " << argv[0] <<
" -p <polynomial> [options]\n" 762 <<
"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." 764 << general_opt <<
"\n";
765 cerr <<
"Example of use:\n" 766 <<
"./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
767 <<
"Example of implicit surfaces (specified by -p):" << endl
768 <<
" - ellipse : 90-3*x^2-2*y^2-z^2 " << endl
769 <<
" - torus : -1*(x^2+y^2+z^2+6*6-2*2)^2+4*6*6*(x^2+y^2) " << endl
770 <<
" - rcube : 6561-x^4-y^4-z^4" << endl
771 <<
" - goursat : 8-0.03*x^4-0.03*y^4-0.03*z^4+2*x^2+2*y^2+2*z^2" << endl
772 <<
" - distel : 10000-(x^2+y^2+z^2+1000*(x^2+y^2)*(x^2+z^2)*(y^2+z^2))" << endl
773 <<
" - leopold : 100-(x^2*y^2*z^2+4*x^2+4*y^2+3*z^2)" << endl
774 <<
" - diabolo : x^2-(y^2+z^2)^2" << endl
775 <<
" - heart : -1*(x^2+2.25*y^2+z^2-1)^3+x^2*z^3+0.1125*y^2*z^3" << endl
776 <<
" - crixxi : -0.9*(y^2+z^2-1)^2-(x^2+y^2-1)^3" << endl << endl;
777 cerr <<
"Implemented estimators (specified by -e):" << endl
778 <<
" - True : supposed to be the ground truth for computations. Of course, it is only approximations." << endl
779 <<
" - VCM : normal estimator by digital Voronoi Covariance Matrix. Radii parameters are given by -R, -r." << endl
780 <<
" - II : normal estimator by Integral Invariants. Radius parameter is given by -r." << endl
781 <<
" - Trivial : the normal obtained by average trivial surfel normals in a ball neighborhood. Radius parameter is given by -r." << endl
783 cerr <<
"Note:" << endl
784 <<
" - 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
785 <<
" - 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;
789 trace.beginBlock(
"Make implicit shape..." );
790 typedef Z3i::Space Space;
791 typedef double Scalar;
792 typedef MPolynomial< 3, Scalar > Polynomial3;
793 typedef MPolynomialReader<3, Scalar> Polynomial3Reader;
794 typedef ImplicitPolynomial3Shape<Space> ImplicitShape;
795 string poly_str = vm[
"polynomial" ].as<
string>();
797 Polynomial3Reader reader;
798 string::const_iterator iter = reader.read( poly, poly_str.begin(), poly_str.end() );
799 if ( iter != poly_str.end() )
801 trace.error() <<
"ERROR reading polynomial: I read only <" << poly_str.substr( 0, iter - poly_str.begin() )
802 <<
">, and I built P=" << poly << std::endl;
805 CountedPtr<ImplicitShape> shape(
new ImplicitShape( poly ) );
808 trace.beginBlock(
"Make implicit digital shape..." );
809 typedef Z3i::KSpace KSpace;
810 typedef KSpace::Point Point;
811 typedef Space::RealPoint RealPoint;
812 typedef GaussDigitizer< Space, ImplicitShape > ImplicitDigitalShape;
813 typedef ImplicitDigitalShape::Domain Domain;
814 Scalar min_x = vm[
"minAABB" ].as<
double>();
815 Scalar max_x = vm[
"maxAABB" ].as<
double>();
816 Scalar h = vm[
"gridstep" ].as<
double>();
817 RealPoint p1( min_x, min_x, min_x );
818 RealPoint p2( max_x, max_x, max_x );
819 CountedPtr<ImplicitDigitalShape> dshape(
new ImplicitDigitalShape() );
820 dshape->attach( *shape );
821 dshape->init( p1, p2, h );
822 Domain domain = dshape->getDomain();
824 K.init( domain.lowerBound(), domain.upperBound(), true );
825 trace.info() <<
"- domain is " << domain << std::endl;
828 chooseSurface( vm, K, *shape, *dshape );