35 #include "DGtal/base/Common.h"
36 #include "DGtal/base/Clock.h"
37 #include "DGtal/helpers/StdDefs.h"
39 #include <boost/program_options/options_description.hpp>
40 #include <boost/program_options/parsers.hpp>
41 #include <boost/program_options/variables_map.hpp>
44 #include "DGtal/shapes/implicit/ImplicitBall.h"
45 #include "DGtal/math/MPolynomial.h"
46 #include "DGtal/io/readers/MPolynomialReader.h"
47 #include "DGtal/shapes/implicit/ImplicitPolynomial3Shape.h"
50 #include "DGtal/shapes/GaussDigitizer.h"
51 #include "DGtal/topology/LightImplicitDigitalSurface.h"
52 #include "DGtal/topology/DigitalSurface.h"
53 #include "DGtal/graph/DepthFirstVisitor.h"
54 #include "DGtal/geometry/volumes/KanungoNoise.h"
55 #include "DGtal/topology/CanonicSCellEmbedder.h"
60 #include "DGtal/kernel/BasicPointFunctors.h"
61 #include "DGtal/geometry/surfaces/FunctorOnCells.h"
63 #include "DGtal/geometry/curves/estimation/TrueLocalEstimatorOnPoints.h"
64 #include "DGtal/geometry/surfaces/estimation/IIGeometricFunctors.h"
65 #include "DGtal/geometry/surfaces/estimation/IntegralInvariantVolumeEstimator.h"
66 #include "DGtal/geometry/surfaces/estimation/IntegralInvariantCovarianceEstimator.h"
68 #include "DGtal/geometry/surfaces/estimation/LocalEstimatorFromSurfelFunctorAdapter.h"
69 #include "DGtal/geometry/surfaces/estimation/estimationFunctors/MongeJetFittingMeanCurvatureEstimator.h"
70 #include "DGtal/geometry/surfaces/estimation/estimationFunctors/MongeJetFittingGaussianCurvatureEstimator.h"
71 #include "DGtal/geometry/surfaces/estimation/estimationFunctors/MongeJetFittingPrincipalCurvaturesEstimator.h"
73 using namespace DGtal;
76 typedef std::pair<double,double> PrincipalCurvatures;
78 template <
typename Shape,
typename KSpace,
typename ConstIterator,
typename OutputIterator >
80 estimateTrueMeanCurvatureQuantity(
const ConstIterator & it_begin,
81 const ConstIterator & it_end,
82 OutputIterator & output,
90 Embedder embedder( K );
91 RealPoint currentRealPoint;
93 for ( ConstIterator it = it_begin; it != it_end; ++it )
95 currentRealPoint = embedder( *it_begin ) * h;
96 *output = aShape->meanCurvature( currentRealPoint );
101 template <
typename Shape,
typename KSpace,
typename ConstIterator,
typename OutputIterator >
103 estimateTrueGaussianCurvatureQuantity(
const ConstIterator & it_begin,
104 const ConstIterator & it_end,
105 OutputIterator & output,
113 Embedder embedder( K );
114 RealPoint currentRealPoint;
116 for ( ConstIterator it = it_begin; it != it_end; ++it )
118 currentRealPoint = embedder( *it_begin ) * h;
119 *output = aShape->gaussianCurvature( currentRealPoint );
124 template <
typename Shape,
typename KSpace,
typename ConstIterator,
typename OutputIterator >
126 estimateTruePrincipalCurvaturesQuantity(
const ConstIterator & it_begin,
127 const ConstIterator & it_end,
128 OutputIterator & output,
136 Embedder embedder( K );
137 RealPoint currentRealPoint;
139 for ( ConstIterator it = it_begin; it != it_end; ++it )
141 currentRealPoint = embedder( *it_begin ) * h;
143 aShape->principalCurvatures( currentRealPoint, k1, k2 );
144 PrincipalCurvatures result;
152 template <
typename Space,
typename Shape>
154 compareShapeEstimators(
const std::string & filename,
155 const Shape * aShape,
156 const typename Space::RealPoint & border_min,
157 const typename Space::RealPoint & border_max,
159 const double & radius_kernel,
160 const double & alpha,
161 const std::string & options,
162 const std::string & properties,
163 const bool & lambda_optimized,
164 double noiseLevel = 0.0 )
166 typedef typename Space::RealPoint RealPoint;
172 bool withNoise = ( noiseLevel <= 0.0 ) ?
false :
true;
174 ASSERT (( noiseLevel < 1.0 ));
176 DigitalShape* dshape =
new DigitalShape();
177 dshape->attach( *aShape );
178 dshape->init( border_min, border_max, h );
181 if ( ! K.
init( dshape->getLowerBound(), dshape->getUpperBound(), true ) )
183 std::cerr <<
"[3dLocalEstimators] error in creating KSpace." << std::endl;
194 typedef typename MyDigitalSurface::ConstIterator ConstIterator;
198 typedef typename VisitorRange::ConstIterator VisitorConstIterator;
204 KanungoPredicate * noisifiedObject =
new KanungoPredicate( *dshape, dshape->getDomain(), noiseLevel );
207 MyDigitalSurface surf ( *boundary );
209 double minsize = dshape->getUpperBound()[0] - dshape->getLowerBound()[0];
210 unsigned int tries = 0;
211 while( surf.size() < 2 * minsize || tries > 150 )
216 surf = MyDigitalSurface( *boundary );
222 std::cerr <<
"Can't found a proper bel. So .... I ... just ... kill myself." << std::endl;
226 VisitorRange * range;
227 VisitorConstIterator ibegin;
228 VisitorConstIterator iend;
234 if( options.at( 0 ) !=
'0' )
237 if( properties.at( 0 ) !=
'0' )
241 char full_filename[360];
242 sprintf( full_filename,
"%s%s", filename.c_str(),
"_True_mean.dat" );
243 std::ofstream file( full_filename );
244 file <<
"# h = " << h << std::endl;
245 file <<
"# True Mean Curvature estimation" << std::endl;
247 std::ostream_iterator< double > out_it_true_mean( file,
"\n" );
249 range =
new VisitorRange(
new Visitor( surf, *surf.begin() ));
250 ibegin = range->begin();
255 estimateTrueMeanCurvatureQuantity( ibegin,
263 file <<
"# time = " << TTrueMeanCurv << std::endl;
272 if( properties.at( 1 ) !=
'0' )
276 char full_filename[360];
277 sprintf( full_filename,
"%s%s", filename.c_str(),
"_True_gaussian.dat" );
278 std::ofstream file( full_filename );
279 file <<
"# h = " << h << std::endl;
280 file <<
"# True Gaussian Curvature estimation" << std::endl;
282 std::ostream_iterator< double > out_it_true_gaussian( file,
"\n" );
284 range =
new VisitorRange(
new Visitor( surf, *surf.begin() ));
285 ibegin = range->begin();
290 estimateTrueGaussianCurvatureQuantity( ibegin,
292 out_it_true_gaussian,
297 double TTrueGaussianCurv = c.
stopClock();
298 file <<
"# time = " << TTrueGaussianCurv << std::endl;
307 if( properties.at( 2 ) !=
'0' )
311 char full_filename[360];
312 sprintf( full_filename,
"%s%s", filename.c_str(),
"_True_principal.dat" );
313 std::ofstream file( full_filename );
314 file <<
"# h = " << h << std::endl;
315 file <<
"# True Gaussian Curvature estimation" << std::endl;
317 std::ostream_iterator< std::string > out_it_true_pc( file,
"\n" );
319 std::vector<PrincipalCurvatures> v_results;
320 std::back_insert_iterator< std::vector<PrincipalCurvatures> > bkIt(v_results);
321 range =
new VisitorRange(
new Visitor( surf, *surf.begin() ));
322 ibegin = range->begin();
327 estimateTruePrincipalCurvaturesQuantity( ibegin,
334 for(
unsigned int ii = 0; ii < v_results.size(); ++ii )
336 std::stringstream ss;
337 ss << v_results[ii].first <<
" " << v_results[ii].second;
338 *out_it_true_pc = ss.str();
342 file <<
"# time = " << TTruePrincCurv << std::endl;
352 double re = radius_kernel * std::pow( h, alpha );
355 if( options.at( 1 ) !=
'0' )
358 if( properties.at( 0 ) !=
'0' )
362 char full_filename[360];
363 sprintf( full_filename,
"%s%s", filename.c_str(),
"_II_mean.dat" );
364 std::ofstream file( full_filename );
365 file <<
"# h = " << h << std::endl;
366 file <<
"# Mean Curvature estimation from the Integral Invariant" << std::endl;
367 file <<
"# computed kernel radius = " << re << std::endl;
369 range =
new VisitorRange(
new Visitor( surf, *surf.begin() ));
370 ibegin = range->begin();
378 MyIICurvatureFunctor curvatureFunctor;
379 curvatureFunctor.
init( h, re );
381 MyIICurvatureEstimator curvatureEstimator( curvatureFunctor );
382 curvatureEstimator.attach( K, *noisifiedObject );
383 curvatureEstimator.setParams( re/h );
384 curvatureEstimator.init( h, ibegin, iend );
386 std::ostream_iterator< double > out_it_ii_mean( file,
"\n" );
387 curvatureEstimator.eval( ibegin, iend, out_it_ii_mean );
390 file <<
"# time = " << TIIMeanCurv << std::endl;
399 if( properties.at( 1 ) !=
'0' )
403 char full_filename[360];
404 sprintf( full_filename,
"%s%s", filename.c_str(),
"_II_gaussian.dat" );
405 std::ofstream file( full_filename );
406 file <<
"# h = " << h << std::endl;
407 file <<
"# Gaussian Curvature estimation from the Integral Invariant" << std::endl;
408 file <<
"# computed kernel radius = " << re << std::endl;
410 range =
new VisitorRange(
new Visitor( surf, *surf.begin() ));
411 ibegin = range->begin();
419 MyIICurvatureFunctor curvatureFunctor;
420 curvatureFunctor.
init( h, re );
422 MyIICurvatureEstimator curvatureEstimator( curvatureFunctor );
423 curvatureEstimator.attach( K, *noisifiedObject );
424 curvatureEstimator.setParams( re/h );
425 curvatureEstimator.init( h, ibegin, iend );
427 std::ostream_iterator< double > out_it_ii_gaussian( file,
"\n" );
428 curvatureEstimator.eval( ibegin, iend, out_it_ii_gaussian );
431 file <<
"# time = " << TIIGaussCurv << std::endl;
440 if( properties.at( 2 ) !=
'0' )
444 char full_filename[360];
445 sprintf( full_filename,
"%s%s", filename.c_str(),
"_II_principal.dat" );
446 std::ofstream file( full_filename );
447 file <<
"# h = " << h << std::endl;
448 file <<
"# Gaussian Curvature estimation from the Integral Invariant" << std::endl;
449 file <<
"# computed kernel radius = " << re << std::endl;
451 range =
new VisitorRange(
new Visitor( surf, *surf.begin() ));
452 ibegin = range->begin();
460 MyIICurvatureFunctor curvatureFunctor;
461 curvatureFunctor.
init( h, re );
463 MyIICurvatureEstimator curvatureEstimator( curvatureFunctor );
464 curvatureEstimator.attach( K, *noisifiedObject );
465 curvatureEstimator.setParams( re/h );
466 curvatureEstimator.init( h, ibegin, iend );
468 std::vector<PrincipalCurvatures> v_results;
469 std::back_insert_iterator< std::vector<PrincipalCurvatures> > bkIt(v_results);
471 curvatureEstimator.eval( ibegin, iend, bkIt );
473 std::ostream_iterator< std::string > out_it_ii_principal( file,
"\n" );
474 for(
unsigned int ii = 0; ii < v_results.size(); ++ii )
476 std::stringstream ss;
477 ss << v_results[ii].first <<
" " << v_results[ii].second;
478 *out_it_ii_principal = ss.str();
479 ++out_it_ii_principal;
483 file <<
"# time = " << TIIGaussCurv << std::endl;
493 if( options.at( 2 ) !=
'0' )
496 if( properties.at( 0 ) !=
'0' )
504 FunctorMean estimatorH( embedder, h );
505 ConvFunctor convFunc(1.0);
507 reporterH.attach( surf );
508 reporterH.setParams(
Z3i::l2Metric, estimatorH, convFunc, re/h );
511 range =
new VisitorRange(
new Visitor( surf, *surf.begin() ));
512 ibegin = range->begin();
515 reporterH.init( h , ibegin, iend );
517 char full_filename[360];
518 sprintf( full_filename,
"%s%s", filename.c_str(),
"_MongeJetFitting_mean.dat" );
519 std::ofstream file( full_filename );
520 file <<
"# h = " << h << std::endl;
521 file <<
"# Mean Curvature estimation from CGAL Monge from and Jet Fitting" << std::endl;
522 file <<
"# computed kernel radius = " << re << std::endl;
523 std::ostream_iterator< double > out_it_monge_mean( file,
"\n" );
526 range =
new VisitorRange(
new Visitor( surf, *surf.begin() ));
527 ibegin = range->begin();
531 reporterH.eval(ibegin, iend, out_it_monge_mean);
533 file <<
"# time = " << TMongeMeanCurv << std::endl;
542 if( properties.at( 1 ) !=
'0' )
550 FunctorGaussian estimatorK( embedder, h );
551 ConvFunctor convFunc(1.0);
553 reporterK.attach( surf );
554 reporterK.setParams(
Z3i::l2Metric, estimatorK, convFunc, re/h );
557 range =
new VisitorRange(
new Visitor( surf, *surf.begin() ));
558 ibegin = range->begin();
561 reporterK.init( h , ibegin, iend );
567 range =
new VisitorRange(
new Visitor( surf, *surf.begin() ));
568 ibegin = range->begin();
571 char full_filename[360];
572 sprintf( full_filename,
"%s%s", filename.c_str(),
"_MongeJetFitting_gaussian.dat" );
573 std::ofstream file( full_filename );
574 file <<
"# h = " << h << std::endl;
575 file <<
"# Gaussian Curvature estimation from CGAL Monge from and Jet Fitting" << std::endl;
576 file <<
"# computed kernel radius = " << re << std::endl;
577 std::ostream_iterator< double > out_it_monge_gaussian( file,
"\n" );
578 reporterK.eval(ibegin, iend , out_it_monge_gaussian);
580 file <<
"# time = " << TMongeGaussCurv << std::endl;
589 if( properties.at( 2 ) !=
'0' )
597 FunctorPrincCurv estimatorK( embedder, h );
598 ConvFunctor convFunc(1.0);
600 reporterK.attach( surf );
601 reporterK.setParams(
Z3i::l2Metric, estimatorK, convFunc, re/h );
605 range =
new VisitorRange(
new Visitor( surf, *surf.begin() ));
606 ibegin = range->begin();
608 reporterK.init( h , ibegin, iend );
610 char full_filename[360];
611 sprintf( full_filename,
"%s%s", filename.c_str(),
"_MongeJetFitting_principal.dat" );
612 std::ofstream file( full_filename );
613 file <<
"# h = " << h << std::endl;
614 file <<
"# Gaussian Curvature estimation from CGAL Monge from and Jet Fitting" << std::endl;
615 file <<
"# computed kernel radius = " << re << std::endl;
616 std::ostream_iterator< std::string > out_it_monge_principal( file,
"\n" );
618 std::vector<PrincipalCurvatures> v_results;
619 std::back_insert_iterator< std::vector<PrincipalCurvatures> > bkIt(v_results);
622 range =
new VisitorRange(
new Visitor( surf, *surf.begin() ));
623 ibegin = range->begin();
626 reporterK.eval(ibegin, iend , bkIt);
628 for(
unsigned int ii = 0; ii < v_results.size(); ++ii )
630 std::stringstream ss;
631 ss << v_results[ii].first <<
" " << v_results[ii].second;
632 *out_it_monge_principal = ss.str();
633 ++out_it_monge_principal;
637 file <<
"# time = " << TMongeGaussCurv << std::endl;
650 typedef typename MyDigitalSurface::ConstIterator ConstIterator;
654 typedef typename VisitorRange::ConstIterator VisitorConstIterator;
662 MyDigitalSurface surf ( boundary );
664 VisitorRange * range;
665 VisitorConstIterator ibegin;
666 VisitorConstIterator iend;
668 unsigned int cntIn = 0;
671 if( dshape->operator ()(*it))
677 std::cout <<
"boundary:" << surf.size() << std::endl;
678 std::cout <<
"full:" << cntIn << std::endl;
684 if( options.at( 0 ) !=
'0' )
687 if( properties.at( 0 ) !=
'0' )
691 char full_filename[360];
692 sprintf( full_filename,
"%s%s", filename.c_str(),
"_True_mean.dat" );
693 std::ofstream file( full_filename );
694 file <<
"# h = " << h << std::endl;
695 file <<
"# True Mean Curvature estimation" << std::endl;
697 std::ostream_iterator< double > out_it_true_mean( file,
"\n" );
699 range =
new VisitorRange(
new Visitor( surf, *surf.begin() ));
700 ibegin = range->begin();
705 estimateTrueMeanCurvatureQuantity( ibegin,
713 file <<
"# time = " << TTrueMeanCurv << std::endl;
722 if( properties.at( 1 ) !=
'0' )
726 char full_filename[360];
727 sprintf( full_filename,
"%s%s", filename.c_str(),
"_True_gaussian.dat" );
728 std::ofstream file( full_filename );
729 file <<
"# h = " << h << std::endl;
730 file <<
"# True Gaussian Curvature estimation" << std::endl;
732 std::ostream_iterator< double > out_it_true_gaussian( file,
"\n" );
734 range =
new VisitorRange(
new Visitor( surf, *surf.begin() ));
735 ibegin = range->begin();
740 estimateTrueGaussianCurvatureQuantity( ibegin,
742 out_it_true_gaussian,
747 double TTrueGaussianCurv = c.
stopClock();
748 file <<
"# time = " << TTrueGaussianCurv << std::endl;
758 if( properties.at( 2 ) !=
'0' )
762 char full_filename[360];
763 sprintf( full_filename,
"%s%s", filename.c_str(),
"_True_principal.dat" );
764 std::ofstream file( full_filename );
765 file <<
"# h = " << h << std::endl;
766 file <<
"# True Gaussian Curvature estimation" << std::endl;
768 std::ostream_iterator< std::string > out_it_true_pc( file,
"\n" );
770 std::vector<PrincipalCurvatures> v_results;
771 std::back_insert_iterator< std::vector<PrincipalCurvatures> > bkIt(v_results);
773 range =
new VisitorRange(
new Visitor( surf, *surf.begin() ));
774 ibegin = range->begin();
779 estimateTruePrincipalCurvaturesQuantity( ibegin,
787 for(
unsigned int ii = 0; ii < v_results.size(); ++ii )
789 std::stringstream ss;
790 ss << v_results[ii].first <<
" " << v_results[ii].second;
791 *out_it_true_pc = ss.str();
796 file <<
"# time = " << TTruePrincCurv << std::endl;
806 double re = radius_kernel * std::pow( h, alpha );
809 if( options.at( 1 ) !=
'0' )
812 if( properties.at( 0 ) !=
'0' )
816 char full_filename[360];
817 sprintf( full_filename,
"%s%s", filename.c_str(),
"_II_mean.dat" );
818 std::ofstream file( full_filename );
819 file <<
"# h = " << h << std::endl;
820 file <<
"# Mean Curvature estimation from the Integral Invariant" << std::endl;
821 file <<
"# computed kernel radius = " << re << std::endl;
823 range =
new VisitorRange(
new Visitor( surf, *surf.begin() ));
824 ibegin = range->begin();
832 MyIICurvatureFunctor curvatureFunctor;
833 curvatureFunctor.
init( h, re );
835 MyIICurvatureEstimator curvatureEstimator( curvatureFunctor );
836 curvatureEstimator.attach( K, *dshape );
837 curvatureEstimator.setParams( re/h );
838 curvatureEstimator.init( h, ibegin, iend );
840 std::ostream_iterator< double > out_it_ii_mean( file,
"\n" );
841 curvatureEstimator.eval( ibegin, iend, out_it_ii_mean );
844 file <<
"# time = " << TIIMeanCurv << std::endl;
853 if( properties.at( 1 ) !=
'0' )
857 char full_filename[360];
858 sprintf( full_filename,
"%s%s", filename.c_str(),
"_II_gaussian.dat" );
859 std::ofstream file( full_filename );
860 file <<
"# h = " << h << std::endl;
861 file <<
"# Gaussian Curvature estimation from the Integral Invariant" << std::endl;
862 file <<
"# computed kernel radius = " << re << std::endl;
864 range =
new VisitorRange(
new Visitor( surf, *surf.begin() ));
865 ibegin = range->begin();
873 MyIICurvatureFunctor curvatureFunctor;
874 curvatureFunctor.
init( h, re );
876 MyIICurvatureEstimator curvatureEstimator( curvatureFunctor );
877 curvatureEstimator.attach( K, *dshape );
878 curvatureEstimator.setParams( re/h );
879 curvatureEstimator.init( h, ibegin, iend );
881 std::ostream_iterator< double > out_it_ii_gaussian( file,
"\n" );
882 curvatureEstimator.eval( ibegin, iend, out_it_ii_gaussian );
885 file <<
"# time = " << TIIGaussCurv << std::endl;
894 if( properties.at( 2 ) !=
'0' )
898 char full_filename[360];
899 sprintf( full_filename,
"%s%s", filename.c_str(),
"_II_principal.dat" );
900 std::ofstream file( full_filename );
901 file <<
"# h = " << h << std::endl;
902 file <<
"# Gaussian Curvature estimation from the Integral Invariant" << std::endl;
903 file <<
"# computed kernel radius = " << re << std::endl;
905 range =
new VisitorRange(
new Visitor( surf, *surf.begin() ));
906 ibegin = range->begin();
914 MyIICurvatureFunctor curvatureFunctor;
915 curvatureFunctor.
init( h, re );
917 MyIICurvatureEstimator curvatureEstimator( curvatureFunctor );
918 curvatureEstimator.attach( K, *dshape );
919 curvatureEstimator.setParams( re/h );
920 curvatureEstimator.init( h, ibegin, iend );
922 std::vector<PrincipalCurvatures> v_results;
923 std::back_insert_iterator< std::vector<PrincipalCurvatures> > bkIt(v_results);
925 curvatureEstimator.eval( ibegin, iend, bkIt );
927 std::ostream_iterator< std::string > out_it_ii_principal( file,
"\n" );
928 for(
unsigned int ii = 0; ii < v_results.size(); ++ii )
930 std::stringstream ss;
931 ss << v_results[ii].first <<
" " << v_results[ii].second;
932 *out_it_ii_principal = ss.str();
933 ++out_it_ii_principal;
937 file <<
"# time = " << TIIGaussCurv << std::endl;
947 if( options.at( 2 ) !=
'0' )
950 if( properties.at( 0 ) !=
'0' )
958 FunctorMean estimatorH( embedder, h );
959 ConvFunctor convFunc(1.0);
961 reporterH.attach( surf );
962 reporterH.setParams(
Z3i::l2Metric, estimatorH, convFunc, re/h );
965 range =
new VisitorRange(
new Visitor( surf, *surf.begin() ));
966 ibegin = range->begin();
968 reporterH.init( h , ibegin, iend );
970 char full_filename[360];
971 sprintf( full_filename,
"%s%s", filename.c_str(),
"_MongeJetFitting_mean.dat" );
972 std::ofstream file( full_filename );
973 file <<
"# h = " << h << std::endl;
974 file <<
"# Mean Curvature estimation from CGAL Monge from and Jet Fitting" << std::endl;
975 file <<
"# computed kernel radius = " << re << std::endl;
976 std::ostream_iterator< double > out_it_monge_mean( file,
"\n" );
979 range =
new VisitorRange(
new Visitor( surf, *surf.begin() ));
980 ibegin = range->begin();
983 reporterH.eval(ibegin, iend , out_it_monge_mean);
985 file <<
"# time = " << TMongeMeanCurv << std::endl;
994 if( properties.at( 1 ) !=
'0' )
1002 FunctorGaussian estimatorK( embedder, h );
1003 ConvFunctor convFunc(1.0);
1004 ReporterK reporterK;
1005 reporterK.attach( surf );
1006 reporterK.setParams(
Z3i::l2Metric, estimatorK, convFunc, re/h );
1010 range =
new VisitorRange(
new Visitor( surf, *surf.begin() ));
1011 ibegin = range->begin();
1012 iend = range->end();
1013 reporterK.init( h , ibegin, iend );
1016 range =
new VisitorRange(
new Visitor( surf, *surf.begin() ));
1017 ibegin = range->begin();
1018 iend = range->end();
1020 char full_filename[360];
1021 sprintf( full_filename,
"%s%s", filename.c_str(),
"_MongeJetFitting_gaussian.dat" );
1022 std::ofstream file( full_filename );
1023 file <<
"# h = " << h << std::endl;
1024 file <<
"# Gaussian Curvature estimation from CGAL Monge from and Jet Fitting" << std::endl;
1025 file <<
"# computed kernel radius = " << re << std::endl;
1026 std::ostream_iterator< double > out_it_monge_gaussian( file,
"\n" );
1027 reporterK.eval(ibegin, iend , out_it_monge_gaussian);
1029 file <<
"# time = " << TMongeGaussCurv << std::endl;
1038 if( properties.at( 2 ) !=
'0' )
1046 FunctorPrincCurv estimatorK( embedder, h );
1047 ConvFunctor convFunc(1.0);
1048 ReporterK reporterK;
1049 reporterK.attach(surf);
1050 reporterK.setParams(
Z3i::l2Metric, estimatorK, convFunc, re/h);
1054 range =
new VisitorRange(
new Visitor( surf, *surf.begin() ));
1055 ibegin = range->begin();
1056 iend = range->end();
1057 reporterK.init( h , ibegin, iend );
1060 range =
new VisitorRange(
new Visitor( surf, *surf.begin() ));
1061 ibegin = range->begin();
1062 iend = range->end();
1064 char full_filename[360];
1065 sprintf( full_filename,
"%s%s", filename.c_str(),
"_MongeJetFitting_principal.dat" );
1066 std::ofstream file( full_filename );
1067 file <<
"# h = " << h << std::endl;
1068 file <<
"# Gaussian Curvature estimation from CGAL Monge from and Jet Fitting" << std::endl;
1069 file <<
"# computed kernel radius = " << re << std::endl;
1070 std::ostream_iterator< std::string > out_it_monge_principal( file,
"\n" );
1072 std::vector<PrincipalCurvatures> v_results;
1073 std::back_insert_iterator< std::vector<PrincipalCurvatures> > bkIt(v_results);
1075 reporterK.eval(ibegin, iend , bkIt);
1077 for(
unsigned int ii = 0; ii < v_results.size(); ++ii )
1079 std::stringstream ss;
1080 ss << v_results[ii].first <<
" " << v_results[ii].second;
1081 *out_it_monge_principal = ss.str();
1082 ++out_it_monge_principal;
1086 file <<
"# time = " << TMongeGaussCurv << std::endl;
1099 std::cerr <<
"[estimatorCurvatureComparator3D]"
1104 << border_min[0] <<
" "
1105 << border_max[0] <<
" "
1107 << radius_kernel <<
" "
1108 << lambda_optimized <<
" "
1125 void missingParam( std::string param )
1127 trace.
error() <<
" Parameter: " << param <<
" is required.";
1132 namespace po = boost::program_options;
1134 int main(
int argc,
char** argv )
1139 #error You need to have activated CGAL (WITH_CGAL) to include this file.
1142 #error You need to have activated EIGEN (WITH_EIGEN) to include this file.
1146 po::options_description general_opt(
"Allowed options are");
1147 general_opt.add_options()
1148 (
"help,h",
"display this message")
1149 (
"shape,s", po::value< std::string >(),
"Shape")
1150 (
"output,o", po::value< std::string >(),
"Output file")
1151 (
"radius,r", po::value< double >(),
"Kernel radius for IntegralInvariant" )
1152 (
"alpha", po::value<double>()->default_value(1.0/3.0),
"Alpha parameter for Integral Invariant computation" )
1153 (
"h", po::value< double >(),
"Grid step" )
1154 (
"minAABB,a", po::value< double >()->default_value( -10.0 ),
"Min value of the AABB bounding box (domain)" )
1155 (
"maxAABB,A", po::value< double >()->default_value( 10.0 ),
"Max value of the AABB bounding box (domain)" )
1156 (
"noise,n", po::value<double>()->default_value(0.0),
"Level of noise to perturb the shape" )
1157 (
"lambda,l", po::value< bool >()->default_value(
false ),
"Use the shape to get a better approximation of the surface (optional)" )
1158 (
"properties", po::value<std::string>()->default_value(
"110"),
"the i-th property is disabled iff there is a 0 at position i" )
1159 (
"estimators,e", po::value< std::string >()->default_value(
"110"),
"the i-th estimator is disabled iff there is a 0 at position i" );
1162 bool parseOK =
true;
1163 po::variables_map vm;
1166 po::store( po::parse_command_line( argc, argv, general_opt ), vm );
1168 catch(
const std::exception & ex )
1171 trace.
info() <<
"Error checking program options: " << ex.what() << std::endl;
1174 if( !parseOK || vm.count(
"help") || argc <= 1 )
1176 trace.
info()<<
"Compare local estimators on implicit shapes using DGtal library" <<std::endl
1177 <<
"Basic usage: "<<std::endl
1178 <<
"\t3dlocalEstimators --shape <shape> --h <h> --radius <radius> --estimators <binaryWord> --output <output>"<<std::endl
1180 <<
"Below are the different available families of estimators: " << std::endl
1181 <<
"\t - Integral Invariant Mean" << std::endl
1182 <<
"\t - Integral Invariant Gaussian" << std::endl
1183 <<
"\t - Monge Jet Fitting Mean" << std::endl
1184 <<
"\t - Monge Jet Fitting Gaussian" << std::endl
1186 <<
"The i-th family of estimators is enabled if the i-th character of the binary word is not 0. "
1187 <<
"The default binary word is '1100'. This means that the first family of estimators, "
1188 <<
"ie. Integral Invariant, is enabled, whereas the next ones are disabled. "
1189 <<
"Below are the different available properties: " << std::endl
1190 <<
"\t - Mean Curvature" << std::endl
1191 <<
"\t - Gaussian Curvature" << std::endl
1192 <<
"\t - k1/k2" << std::endl
1197 if (!(vm.count(
"output"))) missingParam(
"--output");
1198 if (!(vm.count(
"shape"))) missingParam(
"--shape");
1199 if (!(vm.count(
"h"))) missingParam(
"--h");
1200 if (!(vm.count(
"radius"))) missingParam(
"--radius");
1203 std::string file_export = vm[
"output"].as< std::string >();
1205 std::string options = vm[
"estimators"].as< std::string >();
1206 if (options.size() < nb)
1209 <<
" characters are required "
1210 <<
" with option --estimators.";
1214 double h = vm[
"h"].as<
double >();
1215 double radius = vm[
"radius"].as<
double >();
1216 double alpha = vm[
"alpha"].as<
double >();
1217 std::string poly_str = vm[
"shape"].as< std::string >();
1218 bool lambda_optimized = vm[
"lambda"].as<
bool >();
1219 double noiseLevel = vm[
"noise"].as<
double>();
1222 std::string properties = vm[
"properties"].as<std::string>();
1223 if (properties.size() < nb)
1226 <<
" characters are required "
1227 <<
" with option --properties.";
1235 RealPoint border_min( vm[
"minAABB"].as< double >(), vm[
"minAABB"].as< double >(), vm[
"minAABB"].as< double >() );
1236 RealPoint border_max( vm[
"maxAABB"].as< double >(), vm[
"maxAABB"].as< double >(), vm[
"maxAABB"].as< double >() );
1245 Polynomial3Reader reader;
1246 std::string::const_iterator iter = reader.read( poly, poly_str.begin(), poly_str.end() );
1247 if ( iter != poly_str.end() )
1249 std::cerr <<
"ERROR: I read only <"
1250 << poly_str.substr( 0, iter - poly_str.begin() )
1251 <<
">, and I built P=" << poly << std::endl;
1255 ImplicitShape* shape =
new ImplicitShape( poly );
1257 compareShapeEstimators< Z3i::Space, ImplicitShape > (
1260 border_min, border_max,
void beginBlock(const std::string &keyword="")
PointVector< dim, double > RealPoint
static SCell findABel(const KSpace &K, const PointPredicate &pp, unsigned int nbtries=1000)
void init(const double _h, SurfelConstIterator itb, SurfelConstIterator ite)
bool init(const Point &lower, const Point &upper, bool isClosed)
Trace trace(traceWriterTerm)
void init(const double _h, SurfelConstIterator itb, SurfelConstIterator ite)
static const L2Metric l2Metric