From 792f62f41f674bfd2b4e7b6aa37af13666313e06 Mon Sep 17 00:00:00 2001 From: Kerautret Date: Sun, 14 Jun 2020 19:09:17 +0200 Subject: [PATCH 1/6] const alias use --- estimators/3dLocalEstimators.cpp | 42 ++++++++++++++++---------------- estimators/CMakeLists.txt | 2 +- 2 files changed, 22 insertions(+), 22 deletions(-) diff --git a/estimators/3dLocalEstimators.cpp b/estimators/3dLocalEstimators.cpp index 8afc1282..9fe48400 100644 --- a/estimators/3dLocalEstimators.cpp +++ b/estimators/3dLocalEstimators.cpp @@ -173,12 +173,12 @@ compareShapeEstimators( const std::string & filename, ASSERT (( noiseLevel < 1.0 )); // Digitizer - DigitalShape* dshape = new DigitalShape(); - dshape->attach( *aShape ); - dshape->init( border_min, border_max, h ); + DigitalShape dshape; + dshape.attach( *aShape ); + dshape.init( border_min, border_max, h ); KSpace K; - if ( ! K.init( dshape->getLowerBound(), dshape->getUpperBound(), true ) ) + if ( ! K.init( dshape.getLowerBound(), dshape.getUpperBound(), true ) ) { std::cerr << "[3dLocalEstimators] error in creating KSpace." << std::endl; return false; @@ -201,18 +201,19 @@ compareShapeEstimators( const std::string & filename, // typedef FunctorOnCells< MyPointFunctor, KSpace > MySpelFunctor; // Extracts shape boundary - KanungoPredicate * noisifiedObject = new KanungoPredicate( *dshape, dshape->getDomain(), noiseLevel ); - SCell bel = Surfaces< KSpace >::findABel( K, *noisifiedObject, 10000 ); - Boundary * boundary = new Boundary( K, *noisifiedObject, SurfelAdjacency< KSpace::dimension >( true ), bel ); + auto d = dshape.getDomain(); + KanungoPredicate noisifiedObject ( dshape, d, noiseLevel ); + SCell bel = Surfaces< KSpace >::findABel( K, noisifiedObject, 10000 ); + Boundary * boundary = new Boundary( K, noisifiedObject, SurfelAdjacency< KSpace::dimension >( true ), bel ); MyDigitalSurface surf ( *boundary ); - double minsize = dshape->getUpperBound()[0] - dshape->getLowerBound()[0]; + double minsize = dshape.getUpperBound()[0] - dshape.getLowerBound()[0]; unsigned int tries = 0; while( surf.size() < 2 * minsize || tries > 150 ) { delete boundary; - bel = Surfaces< KSpace >::findABel( K, *noisifiedObject, 10000 ); - boundary = new Boundary( K, *noisifiedObject, SurfelAdjacency< KSpace::dimension >( true ), bel ); + bel = Surfaces< KSpace >::findABel( K, noisifiedObject, 10000 ); + boundary = new Boundary( K, noisifiedObject, SurfelAdjacency< KSpace::dimension >( true ), bel ); surf = MyDigitalSurface( *boundary ); ++tries; } @@ -379,7 +380,7 @@ compareShapeEstimators( const std::string & filename, curvatureFunctor.init( h, re ); MyIICurvatureEstimator curvatureEstimator( curvatureFunctor ); - curvatureEstimator.attach( K, *noisifiedObject ); + curvatureEstimator.attach( K, noisifiedObject ); curvatureEstimator.setParams( re/h ); curvatureEstimator.init( h, ibegin, iend ); @@ -420,7 +421,7 @@ compareShapeEstimators( const std::string & filename, curvatureFunctor.init( h, re ); MyIICurvatureEstimator curvatureEstimator( curvatureFunctor ); - curvatureEstimator.attach( K, *noisifiedObject ); + curvatureEstimator.attach( K, noisifiedObject ); curvatureEstimator.setParams( re/h ); curvatureEstimator.init( h, ibegin, iend ); @@ -461,7 +462,7 @@ compareShapeEstimators( const std::string & filename, curvatureFunctor.init( h, re ); MyIICurvatureEstimator curvatureEstimator( curvatureFunctor ); - curvatureEstimator.attach( K, *noisifiedObject ); + curvatureEstimator.attach( K, noisifiedObject ); curvatureEstimator.setParams( re/h ); curvatureEstimator.init( h, ibegin, iend ); @@ -657,8 +658,8 @@ compareShapeEstimators( const std::string & filename, // typedef FunctorOnCells< MyPointFunctor, KSpace > MySpelFunctor; // Extracts shape boundary - SCell bel = Surfaces::findABel ( K, *dshape, 10000 ); - Boundary boundary( K, *dshape, SurfelAdjacency< KSpace::dimension >( true ), bel ); + SCell bel = Surfaces::findABel ( K, dshape, 10000 ); + Boundary boundary( K, dshape, SurfelAdjacency< KSpace::dimension >( true ), bel ); MyDigitalSurface surf ( boundary ); VisitorRange * range; @@ -666,9 +667,9 @@ compareShapeEstimators( const std::string & filename, VisitorConstIterator iend; unsigned int cntIn = 0; - for( typename Z3i::Domain::ConstIterator it = dshape->getDomain().begin(), ite = dshape->getDomain().end(); it != ite; ++it ) + for( typename Z3i::Domain::ConstIterator it = dshape.getDomain().begin(), ite = dshape.getDomain().end(); it != ite; ++it ) { - if( dshape->operator ()(*it)) + if( dshape.operator ()(*it)) { ++cntIn; } @@ -833,7 +834,7 @@ compareShapeEstimators( const std::string & filename, curvatureFunctor.init( h, re ); MyIICurvatureEstimator curvatureEstimator( curvatureFunctor ); - curvatureEstimator.attach( K, *dshape ); + curvatureEstimator.attach( K, dshape ); curvatureEstimator.setParams( re/h ); curvatureEstimator.init( h, ibegin, iend ); @@ -874,7 +875,7 @@ compareShapeEstimators( const std::string & filename, curvatureFunctor.init( h, re ); MyIICurvatureEstimator curvatureEstimator( curvatureFunctor ); - curvatureEstimator.attach( K, *dshape ); + curvatureEstimator.attach( K, dshape ); curvatureEstimator.setParams( re/h ); curvatureEstimator.init( h, ibegin, iend ); @@ -915,7 +916,7 @@ compareShapeEstimators( const std::string & filename, curvatureFunctor.init( h, re ); MyIICurvatureEstimator curvatureEstimator( curvatureFunctor ); - curvatureEstimator.attach( K, *dshape ); + curvatureEstimator.attach( K, dshape ); curvatureEstimator.setParams( re/h ); curvatureEstimator.init( h, ibegin, iend ); @@ -1112,7 +1113,6 @@ compareShapeEstimators( const std::string & filename, return false; } - delete dshape; return true; } diff --git a/estimators/CMakeLists.txt b/estimators/CMakeLists.txt index 5e96a9b8..33874927 100644 --- a/estimators/CMakeLists.txt +++ b/estimators/CMakeLists.txt @@ -36,7 +36,7 @@ endif ( WITH_VISU3D_QGLVIEWER ) if ( WITH_CGAL ) SET(CGAL_TESTS_SRC -# 3dLocalEstimators + 3dLocalEstimators ) FOREACH(FILE ${CGAL_TESTS_SRC}) add_executable(${FILE} ${FILE}) From 30764e1e79adc90293804c0b0c17b6aa6d94f843 Mon Sep 17 00:00:00 2001 From: Kerautret Date: Sat, 18 Jul 2020 14:51:24 +0200 Subject: [PATCH 2/6] Fix compilation following testMonge of DGtal --- estimators/3dLocalEstimators.cpp | 60 +++++++++++++++----------------- 1 file changed, 28 insertions(+), 32 deletions(-) diff --git a/estimators/3dLocalEstimators.cpp b/estimators/3dLocalEstimators.cpp index 9fe48400..259d9485 100644 --- a/estimators/3dLocalEstimators.cpp +++ b/estimators/3dLocalEstimators.cpp @@ -55,10 +55,10 @@ #include "DGtal/topology/CanonicSCellEmbedder.h" - //Estimators #include "DGtal/kernel/BasicPointFunctors.h" #include "DGtal/geometry/surfaces/FunctorOnCells.h" +#include "DGtal/geometry/volumes/distance/LpMetric.h" #include "DGtal/geometry/curves/estimation/TrueLocalEstimatorOnPoints.h" #include "DGtal/geometry/surfaces/estimation/IIGeometricFunctors.h" @@ -133,7 +133,7 @@ estimateTruePrincipalCurvaturesQuantity( const ConstIterator & it_begin, typedef typename KSpace::Space::RealPoint RealPoint; typedef CanonicSCellEmbedder< KSpace > Embedder; - Embedder embedder( K ); + Embedder embedder( K ); RealPoint currentRealPoint; for ( ConstIterator it = it_begin; it != it_end; ++it ) @@ -497,16 +497,17 @@ compareShapeEstimators( const std::string & filename, if( properties.at( 0 ) != '0' ) { trace.beginBlock( "Monge mean curvature" ); - + typedef LpMetric L2Metric; typedef functors::MongeJetFittingMeanCurvatureEstimator > FunctorMean; typedef functors::ConstValue< double > ConvFunctor; - typedef LocalEstimatorFromSurfelFunctorAdapter ReporterH; + typedef LocalEstimatorFromSurfelFunctorAdapter ReporterH; CanonicSCellEmbedder embedder( K ); FunctorMean estimatorH( embedder, h ); ConvFunctor convFunc(1.0); + L2Metric l2(2.0); ReporterH reporterH; reporterH.attach( surf ); - reporterH.setParams( Z3i::l2Metric, estimatorH, convFunc, re/h ); + reporterH.setParams( l2, estimatorH, convFunc, re/h ); c.startClock(); range = new VisitorRange( new Visitor( surf, *surf.begin() )); @@ -527,8 +528,6 @@ compareShapeEstimators( const std::string & filename, range = new VisitorRange( new Visitor( surf, *surf.begin() )); ibegin = range->begin(); iend = range->end(); - //typename ReporterH::SurfelConstIterator aabegin = surf.begin(); - //typename ReporterH::SurfelConstIterator aaend = surf.end(); reporterH.eval(ibegin, iend, out_it_monge_mean); double TMongeMeanCurv = c.stopClock(); file << "# time = " << TMongeMeanCurv << std::endl; @@ -546,13 +545,15 @@ compareShapeEstimators( const std::string & filename, typedef functors::MongeJetFittingGaussianCurvatureEstimator > FunctorGaussian; typedef functors::ConstValue< double > ConvFunctor; - typedef LocalEstimatorFromSurfelFunctorAdapter ReporterK; + typedef LpMetric L2Metric; + typedef LocalEstimatorFromSurfelFunctorAdapter ReporterK; CanonicSCellEmbedder embedder( K ); FunctorGaussian estimatorK( embedder, h ); ConvFunctor convFunc(1.0); ReporterK reporterK; reporterK.attach( surf ); - reporterK.setParams( Z3i::l2Metric, estimatorK, convFunc, re/h ); + L2Metric l2(2.0); + reporterK.setParams( l2, estimatorK, convFunc, re/h ); c.startClock(); range = new VisitorRange( new Visitor( surf, *surf.begin() )); @@ -561,9 +562,6 @@ compareShapeEstimators( const std::string & filename, reporterK.init( h , ibegin, iend ); - //typename ReporterK::SurfelConstIterator aaabegin = surf.begin(); - //typename ReporterK::SurfelConstIterator aaaend = surf.end(); - delete range; range = new VisitorRange( new Visitor( surf, *surf.begin() )); ibegin = range->begin(); @@ -590,16 +588,17 @@ compareShapeEstimators( const std::string & filename, if( properties.at( 2 ) != '0' ) { trace.beginBlock( "Monge Principal Curvature" ); - + typedef LpMetric L2Metric; typedef functors::MongeJetFittingPrincipalCurvaturesEstimator > FunctorPrincCurv; typedef functors::ConstValue< double > ConvFunctor; - typedef LocalEstimatorFromSurfelFunctorAdapter ReporterK; + typedef LocalEstimatorFromSurfelFunctorAdapter ReporterK; CanonicSCellEmbedder embedder( K ); FunctorPrincCurv estimatorK( embedder, h ); ConvFunctor convFunc(1.0); + L2Metric l2(2.0); ReporterK reporterK; reporterK.attach( surf ); - reporterK.setParams( Z3i::l2Metric, estimatorK, convFunc, re/h ); + reporterK.setParams( l2, estimatorK, convFunc, re/h ); c.startClock(); @@ -617,13 +616,12 @@ compareShapeEstimators( const std::string & filename, std::ostream_iterator< std::string > out_it_monge_principal( file, "\n" ); std::vector v_results; - std::back_insert_iterator< std::vector > bkIt(v_results); + std::back_insert_iterator > bkIt(v_results); delete range; range = new VisitorRange( new Visitor( surf, *surf.begin() )); ibegin = range->begin(); iend = range->end(); - reporterK.eval(ibegin, iend , bkIt);//out_it_monge_principal); for(unsigned int ii = 0; ii < v_results.size(); ++ii ) @@ -654,9 +652,6 @@ compareShapeEstimators( const std::string & filename, typedef GraphVisitorRange< Visitor > VisitorRange; typedef typename VisitorRange::ConstIterator VisitorConstIterator; - // typedef PointFunctorFromPointPredicateAndDomain< DigitalShape, Z3i::Domain, unsigned int > MyPointFunctor; - // typedef FunctorOnCells< MyPointFunctor, KSpace > MySpelFunctor; - // Extracts shape boundary SCell bel = Surfaces::findABel ( K, dshape, 10000 ); Boundary boundary( K, dshape, SurfelAdjacency< KSpace::dimension >( true ), bel ); @@ -951,16 +946,17 @@ compareShapeEstimators( const std::string & filename, if( properties.at( 0 ) != '0' ) { trace.beginBlock( "Monge mean curvature" ); - + typedef LpMetric L2Metric; typedef functors::MongeJetFittingMeanCurvatureEstimator > FunctorMean; typedef functors::ConstValue< double > ConvFunctor; - typedef LocalEstimatorFromSurfelFunctorAdapter ReporterH; + typedef LocalEstimatorFromSurfelFunctorAdapter ReporterH; CanonicSCellEmbedder embedder( K ); FunctorMean estimatorH( embedder, h ); ConvFunctor convFunc(1.0); + L2Metric l2(2.0); ReporterH reporterH; reporterH.attach( surf ); - reporterH.setParams( Z3i::l2Metric, estimatorH, convFunc, re/h ); + reporterH.setParams( l2, estimatorH, convFunc, re/h ); c.startClock(); range = new VisitorRange( new Visitor( surf, *surf.begin() )); @@ -995,16 +991,17 @@ compareShapeEstimators( const std::string & filename, if( properties.at( 1 ) != '0' ) { trace.beginBlock( "Monge Gaussian curvature" ); - + typedef LpMetric L2Metric; typedef functors::MongeJetFittingGaussianCurvatureEstimator > FunctorGaussian; typedef functors::ConstValue< double > ConvFunctor; - typedef LocalEstimatorFromSurfelFunctorAdapter ReporterK; + typedef LocalEstimatorFromSurfelFunctorAdapter ReporterK; CanonicSCellEmbedder embedder( K ); FunctorGaussian estimatorK( embedder, h ); ConvFunctor convFunc(1.0); ReporterK reporterK; reporterK.attach( surf ); - reporterK.setParams( Z3i::l2Metric, estimatorK, convFunc, re/h ); + L2Metric l2(2.0); + reporterK.setParams( l2, estimatorK, convFunc, re/h ); c.startClock(); @@ -1039,19 +1036,19 @@ compareShapeEstimators( const std::string & filename, if( properties.at( 2 ) != '0' ) { trace.beginBlock( "Monge Principal Curvature" ); + typedef LpMetric L2Metric; typedef functors::MongeJetFittingPrincipalCurvaturesEstimator > FunctorPrincCurv; typedef functors::ConstValue< double > ConvFunctor; - typedef LocalEstimatorFromSurfelFunctorAdapter ReporterK; + typedef LocalEstimatorFromSurfelFunctorAdapter ReporterK; CanonicSCellEmbedder embedder( K ); FunctorPrincCurv estimatorK( embedder, h ); ConvFunctor convFunc(1.0); ReporterK reporterK; reporterK.attach(surf); - reporterK.setParams(Z3i::l2Metric, estimatorK, convFunc, re/h); - + L2Metric l2(2.0); + reporterK.setParams(l2, estimatorK, convFunc, re/h); c.startClock(); - range = new VisitorRange( new Visitor( surf, *surf.begin() )); ibegin = range->begin(); iend = range->end(); @@ -1071,8 +1068,7 @@ compareShapeEstimators( const std::string & filename, std::ostream_iterator< std::string > out_it_monge_principal( file, "\n" ); std::vector v_results; - std::back_insert_iterator< std::vector > bkIt(v_results); - + std::back_insert_iterator > bkIt(v_results); reporterK.eval(ibegin, iend , bkIt);//out_it_monge_principal); for(unsigned int ii = 0; ii < v_results.size(); ++ii ) From 3ab21c86094b35b71b8a322b47af0665aa886227 Mon Sep 17 00:00:00 2001 From: Kerautret Date: Sat, 18 Jul 2020 15:26:14 +0200 Subject: [PATCH 3/6] changelog --- ChangeLog.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/ChangeLog.md b/ChangeLog.md index 71e63ac1..e4d6b4f8 100644 --- a/ChangeLog.md +++ b/ChangeLog.md @@ -39,6 +39,8 @@ - *estimators* - volSurfaceRegularization now in the "make install" command. (David Coeurjolly, [#376](https://github.com/DGtal-team/DGtalTools/pull/376)) + - 3dLocalEstimators now include in the main build and fix compilation issues #382. + (Bertrand Kerautret [#383](https://github.com/DGtal-team/DGtalTools/pull/383)) - *converters* - volBoundary2obj improved using new Shortcuts helpers. From 09ac9109224d2ce767726c4e66187a3516678208 Mon Sep 17 00:00:00 2001 From: Kerautret Date: Sat, 18 Jul 2020 15:29:39 +0200 Subject: [PATCH 4/6] changelog --- ChangeLog.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/ChangeLog.md b/ChangeLog.md index e4d6b4f8..e62a905d 100644 --- a/ChangeLog.md +++ b/ChangeLog.md @@ -21,7 +21,6 @@ - volAddBorder: Add an option that sets zero value to domain boundary voxels without changing the domain extent (Bertrand Kertautret [#371](https://github.com/DGtal-team/DGtalTools/pull/371)) - - Fix a wrong error message that appears when using the tool (wrong IO error) (Bertrand Kerautret [#368](https://github.com/DGtal-team/DGtalTools/pull/368)) @@ -39,9 +38,10 @@ - *estimators* - volSurfaceRegularization now in the "make install" command. (David Coeurjolly, [#376](https://github.com/DGtal-team/DGtalTools/pull/376)) - - 3dLocalEstimators now include in the main build and fix compilation issues #382. + - 3dLocalEstimators now include in the main build and fix compilation issues + [#382](https://github.com/DGtal-team/DGtalTools/issues/382). (Bertrand Kerautret [#383](https://github.com/DGtal-team/DGtalTools/pull/383)) - + - *converters* - volBoundary2obj improved using new Shortcuts helpers. (Bertrand Kerautret From 40c12eb8017a244af9379252f1545f2f66debc28 Mon Sep 17 00:00:00 2001 From: Kerautret Date: Sat, 18 Jul 2020 16:46:25 +0200 Subject: [PATCH 5/6] indent --- estimators/3dLocalEstimators.cpp | 2261 +++++++++++++++--------------- 1 file changed, 1129 insertions(+), 1132 deletions(-) diff --git a/estimators/3dLocalEstimators.cpp b/estimators/3dLocalEstimators.cpp index 259d9485..75b191d5 100644 --- a/estimators/3dLocalEstimators.cpp +++ b/estimators/3dLocalEstimators.cpp @@ -78,1039 +78,1036 @@ typedef std::pair PrincipalCurvatures; template < typename Shape, typename KSpace, typename ConstIterator, typename OutputIterator > void estimateTrueMeanCurvatureQuantity( const ConstIterator & it_begin, - const ConstIterator & it_end, - OutputIterator & output, - const KSpace & K, - const double & h, - Shape * aShape ) + const ConstIterator & it_end, + OutputIterator & output, + const KSpace & K, + const double & h, + Shape * aShape ) { - typedef typename KSpace::Space::RealPoint RealPoint; - typedef CanonicSCellEmbedder< KSpace > Embedder; - - Embedder embedder( K ); - RealPoint currentRealPoint; - - for ( ConstIterator it = it_begin; it != it_end; ++it ) - { - currentRealPoint = embedder( *it_begin ) * h; - *output = aShape->meanCurvature( currentRealPoint ); - ++output; - } + typedef typename KSpace::Space::RealPoint RealPoint; + typedef CanonicSCellEmbedder< KSpace > Embedder; + + Embedder embedder( K ); + RealPoint currentRealPoint; + + for ( ConstIterator it = it_begin; it != it_end; ++it ) + { + currentRealPoint = embedder( *it_begin ) * h; + *output = aShape->meanCurvature( currentRealPoint ); + ++output; + } } template < typename Shape, typename KSpace, typename ConstIterator, typename OutputIterator > void estimateTrueGaussianCurvatureQuantity( const ConstIterator & it_begin, - const ConstIterator & it_end, - OutputIterator & output, - const KSpace & K, - const double & h, - Shape * aShape ) + const ConstIterator & it_end, + OutputIterator & output, + const KSpace & K, + const double & h, + Shape * aShape ) { - typedef typename KSpace::Space::RealPoint RealPoint; - typedef CanonicSCellEmbedder< KSpace > Embedder; - - Embedder embedder( K ); - RealPoint currentRealPoint; - - for ( ConstIterator it = it_begin; it != it_end; ++it ) - { - currentRealPoint = embedder( *it_begin ) * h; - *output = aShape->gaussianCurvature( currentRealPoint ); - ++output; - } + typedef typename KSpace::Space::RealPoint RealPoint; + typedef CanonicSCellEmbedder< KSpace > Embedder; + + Embedder embedder( K ); + RealPoint currentRealPoint; + + for ( ConstIterator it = it_begin; it != it_end; ++it ) + { + currentRealPoint = embedder( *it_begin ) * h; + *output = aShape->gaussianCurvature( currentRealPoint ); + ++output; + } } template < typename Shape, typename KSpace, typename ConstIterator, typename OutputIterator > void estimateTruePrincipalCurvaturesQuantity( const ConstIterator & it_begin, - const ConstIterator & it_end, - OutputIterator & output, - const KSpace & K, - const double & h, - Shape * aShape ) + const ConstIterator & it_end, + OutputIterator & output, + const KSpace & K, + const double & h, + Shape * aShape ) { - typedef typename KSpace::Space::RealPoint RealPoint; - typedef CanonicSCellEmbedder< KSpace > Embedder; - + typedef typename KSpace::Space::RealPoint RealPoint; + typedef CanonicSCellEmbedder< KSpace > Embedder; + Embedder embedder( K ); - RealPoint currentRealPoint; - - for ( ConstIterator it = it_begin; it != it_end; ++it ) - { - currentRealPoint = embedder( *it_begin ) * h; - double k1, k2; - aShape->principalCurvatures( currentRealPoint, k1, k2 ); - PrincipalCurvatures result; - result.first = k1; - result.second = k2; - *output = result; - ++output; - } + RealPoint currentRealPoint; + + for ( ConstIterator it = it_begin; it != it_end; ++it ) + { + currentRealPoint = embedder( *it_begin ) * h; + double k1, k2; + aShape->principalCurvatures( currentRealPoint, k1, k2 ); + PrincipalCurvatures result; + result.first = k1; + result.second = k2; + *output = result; + ++output; + } } template bool compareShapeEstimators( const std::string & filename, - const Shape * aShape, - const typename Space::RealPoint & border_min, - const typename Space::RealPoint & border_max, - const double & h, - const double & radius_kernel, - const double & alpha, - const std::string & options, - const std::string & properties, - const bool & lambda_optimized, - double noiseLevel = 0.0 ) + const Shape * aShape, + const typename Space::RealPoint & border_min, + const typename Space::RealPoint & border_max, + const double & h, + const double & radius_kernel, + const double & alpha, + const std::string & options, + const std::string & properties, + const bool & lambda_optimized, + double noiseLevel = 0.0 ) { - typedef typename Space::RealPoint RealPoint; - typedef GaussDigitizer< Z3i::Space, Shape > DigitalShape; - typedef Z3i::KSpace KSpace; - typedef typename KSpace::SCell SCell; - typedef typename KSpace::Surfel Surfel; - - bool withNoise = ( noiseLevel <= 0.0 ) ? false : true; - - ASSERT (( noiseLevel < 1.0 )); - // Digitizer - DigitalShape dshape; + typedef typename Space::RealPoint RealPoint; + typedef GaussDigitizer< Z3i::Space, Shape > DigitalShape; + typedef Z3i::KSpace KSpace; + typedef typename KSpace::SCell SCell; + typedef typename KSpace::Surfel Surfel; + + bool withNoise = ( noiseLevel <= 0.0 ) ? false : true; + + ASSERT (( noiseLevel < 1.0 )); + // Digitizer + DigitalShape dshape; dshape.attach( *aShape ); dshape.init( border_min, border_max, h ); - - KSpace K; + + KSpace K; if ( ! K.init( dshape.getLowerBound(), dshape.getUpperBound(), true ) ) + { + std::cerr << "[3dLocalEstimators] error in creating KSpace." << std::endl; + return false; + } + + try + { + if ( withNoise ) { - std::cerr << "[3dLocalEstimators] error in creating KSpace." << std::endl; + typedef KanungoNoise< DigitalShape, Z3i::Domain > KanungoPredicate; + typedef LightImplicitDigitalSurface< KSpace, KanungoPredicate > Boundary; + typedef DigitalSurface< Boundary > MyDigitalSurface; + typedef typename MyDigitalSurface::ConstIterator ConstIterator; + + typedef DepthFirstVisitor< MyDigitalSurface > Visitor; + typedef GraphVisitorRange< Visitor > VisitorRange; + typedef typename VisitorRange::ConstIterator VisitorConstIterator; + + // typedef PointFunctorFromPointPredicateAndDomain< KanungoPredicate, Z3i::Domain, unsigned int > MyPointFunctor; + // typedef FunctorOnCells< MyPointFunctor, KSpace > MySpelFunctor; + + // Extracts shape boundary + auto d = dshape.getDomain(); + KanungoPredicate noisifiedObject ( dshape, d, noiseLevel ); + SCell bel = Surfaces< KSpace >::findABel( K, noisifiedObject, 10000 ); + Boundary * boundary = new Boundary( K, noisifiedObject, SurfelAdjacency< KSpace::dimension >( true ), bel ); + MyDigitalSurface surf ( *boundary ); + + double minsize = dshape.getUpperBound()[0] - dshape.getLowerBound()[0]; + unsigned int tries = 0; + while( surf.size() < 2 * minsize || tries > 150 ) + { + delete boundary; + bel = Surfaces< KSpace >::findABel( K, noisifiedObject, 10000 ); + boundary = new Boundary( K, noisifiedObject, SurfelAdjacency< KSpace::dimension >( true ), bel ); + surf = MyDigitalSurface( *boundary ); + ++tries; + } + + if( tries > 150 ) + { + std::cerr << "Can't found a proper bel. So .... I ... just ... kill myself." << std::endl; return false; - } - - try - { - if ( withNoise ) + } + + VisitorRange * range; + VisitorConstIterator ibegin; + VisitorConstIterator iend; + + // Estimations + Clock c; + + // True + if( options.at( 0 ) != '0' ) + { + // True Mean Curvature + if( properties.at( 0 ) != '0' ) { - typedef KanungoNoise< DigitalShape, Z3i::Domain > KanungoPredicate; - typedef LightImplicitDigitalSurface< KSpace, KanungoPredicate > Boundary; - typedef DigitalSurface< Boundary > MyDigitalSurface; - typedef typename MyDigitalSurface::ConstIterator ConstIterator; - - typedef DepthFirstVisitor< MyDigitalSurface > Visitor; - typedef GraphVisitorRange< Visitor > VisitorRange; - typedef typename VisitorRange::ConstIterator VisitorConstIterator; - - // typedef PointFunctorFromPointPredicateAndDomain< KanungoPredicate, Z3i::Domain, unsigned int > MyPointFunctor; - // typedef FunctorOnCells< MyPointFunctor, KSpace > MySpelFunctor; - - // Extracts shape boundary - auto d = dshape.getDomain(); - KanungoPredicate noisifiedObject ( dshape, d, noiseLevel ); - SCell bel = Surfaces< KSpace >::findABel( K, noisifiedObject, 10000 ); - Boundary * boundary = new Boundary( K, noisifiedObject, SurfelAdjacency< KSpace::dimension >( true ), bel ); - MyDigitalSurface surf ( *boundary ); - - double minsize = dshape.getUpperBound()[0] - dshape.getLowerBound()[0]; - unsigned int tries = 0; - while( surf.size() < 2 * minsize || tries > 150 ) - { - delete boundary; - bel = Surfaces< KSpace >::findABel( K, noisifiedObject, 10000 ); - boundary = new Boundary( K, noisifiedObject, SurfelAdjacency< KSpace::dimension >( true ), bel ); - surf = MyDigitalSurface( *boundary ); - ++tries; - } - - if( tries > 150 ) - { - std::cerr << "Can't found a proper bel. So .... I ... just ... kill myself." << std::endl; - return false; - } - - VisitorRange * range; - VisitorConstIterator ibegin; - VisitorConstIterator iend; - - // Estimations - Clock c; - - // True - if( options.at( 0 ) != '0' ) - { - // True Mean Curvature - if( properties.at( 0 ) != '0' ) - { - trace.beginBlock( "True mean curvature" ); - - char full_filename[360]; - sprintf( full_filename, "%s%s", filename.c_str(), "_True_mean.dat" ); - std::ofstream file( full_filename ); - file << "# h = " << h << std::endl; - file << "# True Mean Curvature estimation" << std::endl; - - std::ostream_iterator< double > out_it_true_mean( file, "\n" ); - - range = new VisitorRange( new Visitor( surf, *surf.begin() )); - ibegin = range->begin(); - iend = range->end(); - - c.startClock(); - - estimateTrueMeanCurvatureQuantity( ibegin, - iend, - out_it_true_mean, - K, - h, - aShape ); - - double TTrueMeanCurv = c.stopClock(); - file << "# time = " << TTrueMeanCurv << std::endl; - - file.close(); - delete range; - - trace.endBlock(); - } - - // True Gaussian Curvature - if( properties.at( 1 ) != '0' ) - { - trace.beginBlock( "True Gausian curvature" ); - - char full_filename[360]; - sprintf( full_filename, "%s%s", filename.c_str(), "_True_gaussian.dat" ); - std::ofstream file( full_filename ); - file << "# h = " << h << std::endl; - file << "# True Gaussian Curvature estimation" << std::endl; - - std::ostream_iterator< double > out_it_true_gaussian( file, "\n" ); - - range = new VisitorRange( new Visitor( surf, *surf.begin() )); - ibegin = range->begin(); - iend = range->end(); - - c.startClock(); - - estimateTrueGaussianCurvatureQuantity( ibegin, - iend, - out_it_true_gaussian, - K, - h, - aShape ); - - double TTrueGaussianCurv = c.stopClock(); - file << "# time = " << TTrueGaussianCurv << std::endl; - - file.close(); - delete range; - - trace.endBlock(); - } - - // True Principal Curvatures - if( properties.at( 2 ) != '0' ) - { - trace.beginBlock( "True principal curvatures" ); - - char full_filename[360]; - sprintf( full_filename, "%s%s", filename.c_str(), "_True_principal.dat" ); - std::ofstream file( full_filename ); - file << "# h = " << h << std::endl; - file << "# True Gaussian Curvature estimation" << std::endl; - - std::ostream_iterator< std::string > out_it_true_pc( file, "\n" ); - - std::vector v_results; - std::back_insert_iterator< std::vector > bkIt(v_results); - range = new VisitorRange( new Visitor( surf, *surf.begin() )); - ibegin = range->begin(); - iend = range->end(); - - c.startClock(); - - estimateTruePrincipalCurvaturesQuantity( ibegin, - iend, - bkIt,//out_it_true_pc, - K, - h, - aShape ); - - for(unsigned int ii = 0; ii < v_results.size(); ++ii ) - { - std::stringstream ss; - ss << v_results[ii].first << " " << v_results[ii].second; - *out_it_true_pc = ss.str(); - ++out_it_true_pc; - } - double TTruePrincCurv = c.stopClock(); - file << "# time = " << TTruePrincCurv << std::endl; - - file.close(); - - delete range; - - trace.endBlock(); - } - } - - double re = radius_kernel * std::pow( h, alpha ); // to obtains convergence results, re must follow the rule re=kh^(1/3) - - // II - if( options.at( 1 ) != '0' ) - { - // Integral Invariant Mean Curvature - if( properties.at( 0 ) != '0' ) - { - trace.beginBlock( "II mean curvature" ); - - char full_filename[360]; - sprintf( full_filename, "%s%s", filename.c_str(), "_II_mean.dat" ); - std::ofstream file( full_filename ); - file << "# h = " << h << std::endl; - file << "# Mean Curvature estimation from the Integral Invariant" << std::endl; - file << "# computed kernel radius = " << re << std::endl; - - range = new VisitorRange( new Visitor( surf, *surf.begin() )); - ibegin = range->begin(); - iend = range->end(); - - c.startClock(); - - typedef functors::IIMeanCurvature3DFunctor MyIICurvatureFunctor; - typedef IntegralInvariantVolumeEstimator< KSpace, KanungoPredicate, MyIICurvatureFunctor > MyIICurvatureEstimator; - - MyIICurvatureFunctor curvatureFunctor; - curvatureFunctor.init( h, re ); - - MyIICurvatureEstimator curvatureEstimator( curvatureFunctor ); - curvatureEstimator.attach( K, noisifiedObject ); - curvatureEstimator.setParams( re/h ); - curvatureEstimator.init( h, ibegin, iend ); - - std::ostream_iterator< double > out_it_ii_mean( file, "\n" ); - curvatureEstimator.eval( ibegin, iend, out_it_ii_mean ); - - double TIIMeanCurv = c.stopClock(); - file << "# time = " << TIIMeanCurv << std::endl; - file.close(); - - delete range; - - trace.endBlock(); - } - - // Integral Invariant Gaussian Curvature - if( properties.at( 1 ) != '0' ) - { - trace.beginBlock( "II Gaussian curvature" ); - - char full_filename[360]; - sprintf( full_filename, "%s%s", filename.c_str(), "_II_gaussian.dat" ); - std::ofstream file( full_filename ); - file << "# h = " << h << std::endl; - file << "# Gaussian Curvature estimation from the Integral Invariant" << std::endl; - file << "# computed kernel radius = " << re << std::endl; - - range = new VisitorRange( new Visitor( surf, *surf.begin() )); - ibegin = range->begin(); - iend = range->end(); - - c.startClock(); - - typedef functors::IIGaussianCurvature3DFunctor MyIICurvatureFunctor; - typedef IntegralInvariantCovarianceEstimator< KSpace, KanungoPredicate, MyIICurvatureFunctor > MyIICurvatureEstimator; - - MyIICurvatureFunctor curvatureFunctor; - curvatureFunctor.init( h, re ); - - MyIICurvatureEstimator curvatureEstimator( curvatureFunctor ); - curvatureEstimator.attach( K, noisifiedObject ); - curvatureEstimator.setParams( re/h ); - curvatureEstimator.init( h, ibegin, iend ); - - std::ostream_iterator< double > out_it_ii_gaussian( file, "\n" ); - curvatureEstimator.eval( ibegin, iend, out_it_ii_gaussian ); - - double TIIGaussCurv = c.stopClock(); - file << "# time = " << TIIGaussCurv << std::endl; - file.close(); - - delete range; - - trace.endBlock(); - } - - // Integral Invariant Principal Curvatures - if( properties.at( 2 ) != '0' ) - { - trace.beginBlock( "II Principal curvatures" ); - - char full_filename[360]; - sprintf( full_filename, "%s%s", filename.c_str(), "_II_principal.dat" ); - std::ofstream file( full_filename ); - file << "# h = " << h << std::endl; - file << "# Gaussian Curvature estimation from the Integral Invariant" << std::endl; - file << "# computed kernel radius = " << re << std::endl; - - range = new VisitorRange( new Visitor( surf, *surf.begin() )); - ibegin = range->begin(); - iend = range->end(); - - c.startClock(); - - typedef functors::IIPrincipalCurvatures3DFunctor MyIICurvatureFunctor; - typedef IntegralInvariantCovarianceEstimator< KSpace, KanungoPredicate, MyIICurvatureFunctor > MyIICurvatureEstimator; - - MyIICurvatureFunctor curvatureFunctor; - curvatureFunctor.init( h, re ); - - MyIICurvatureEstimator curvatureEstimator( curvatureFunctor ); - curvatureEstimator.attach( K, noisifiedObject ); - curvatureEstimator.setParams( re/h ); - curvatureEstimator.init( h, ibegin, iend ); - - std::vector v_results; - std::back_insert_iterator< std::vector > bkIt(v_results); - - curvatureEstimator.eval( ibegin, iend, bkIt ); - - std::ostream_iterator< std::string > out_it_ii_principal( file, "\n" ); - for( unsigned int ii = 0; ii < v_results.size(); ++ii ) - { - std::stringstream ss; - ss << v_results[ii].first << " " << v_results[ii].second; - *out_it_ii_principal = ss.str(); - ++out_it_ii_principal; - } - - double TIIGaussCurv = c.stopClock(); - file << "# time = " << TIIGaussCurv << std::endl; - file.close(); - - delete range; - - trace.endBlock(); - } - } - - // Monge - if( options.at( 2 ) != '0' ) - { - // Monge Mean Curvature - if( properties.at( 0 ) != '0' ) - { - trace.beginBlock( "Monge mean curvature" ); - typedef LpMetric L2Metric; - typedef functors::MongeJetFittingMeanCurvatureEstimator > FunctorMean; - typedef functors::ConstValue< double > ConvFunctor; - typedef LocalEstimatorFromSurfelFunctorAdapter ReporterH; - CanonicSCellEmbedder embedder( K ); - FunctorMean estimatorH( embedder, h ); - ConvFunctor convFunc(1.0); - L2Metric l2(2.0); - ReporterH reporterH; - reporterH.attach( surf ); - reporterH.setParams( l2, estimatorH, convFunc, re/h ); - - c.startClock(); - range = new VisitorRange( new Visitor( surf, *surf.begin() )); - ibegin = range->begin(); - iend = range->end(); - - reporterH.init( h , ibegin, iend ); - - char full_filename[360]; - sprintf( full_filename, "%s%s", filename.c_str(), "_MongeJetFitting_mean.dat" ); - std::ofstream file( full_filename ); - file << "# h = " << h << std::endl; - file << "# Mean Curvature estimation from CGAL Monge from and Jet Fitting" << std::endl; - file << "# computed kernel radius = " << re << std::endl; - std::ostream_iterator< double > out_it_monge_mean( file, "\n" ); - - delete range; - range = new VisitorRange( new Visitor( surf, *surf.begin() )); - ibegin = range->begin(); - iend = range->end(); - reporterH.eval(ibegin, iend, out_it_monge_mean); - double TMongeMeanCurv = c.stopClock(); - file << "# time = " << TMongeMeanCurv << std::endl; - file.close(); - delete range; - - - trace.endBlock(); - } - - // Monge Gaussian Curvature - if( properties.at( 1 ) != '0' ) - { - trace.beginBlock( "Monge Gaussian curvature" ); - - typedef functors::MongeJetFittingGaussianCurvatureEstimator > FunctorGaussian; - typedef functors::ConstValue< double > ConvFunctor; - typedef LpMetric L2Metric; - typedef LocalEstimatorFromSurfelFunctorAdapter ReporterK; - CanonicSCellEmbedder embedder( K ); - FunctorGaussian estimatorK( embedder, h ); - ConvFunctor convFunc(1.0); - ReporterK reporterK; - reporterK.attach( surf ); - L2Metric l2(2.0); - reporterK.setParams( l2, estimatorK, convFunc, re/h ); - c.startClock(); - - range = new VisitorRange( new Visitor( surf, *surf.begin() )); - ibegin = range->begin(); - iend = range->end(); - - reporterK.init( h , ibegin, iend ); - - delete range; - range = new VisitorRange( new Visitor( surf, *surf.begin() )); - ibegin = range->begin(); - iend = range->end(); - - char full_filename[360]; - sprintf( full_filename, "%s%s", filename.c_str(), "_MongeJetFitting_gaussian.dat" ); - std::ofstream file( full_filename ); - file << "# h = " << h << std::endl; - file << "# Gaussian Curvature estimation from CGAL Monge from and Jet Fitting" << std::endl; - file << "# computed kernel radius = " << re << std::endl; - std::ostream_iterator< double > out_it_monge_gaussian( file, "\n" ); - reporterK.eval(ibegin, iend , out_it_monge_gaussian); - double TMongeGaussCurv = c.stopClock(); - file << "# time = " << TMongeGaussCurv << std::endl; - file.close(); - delete range; - - - trace.endBlock(); - } - - // Monge Principal Curvatures - if( properties.at( 2 ) != '0' ) - { - trace.beginBlock( "Monge Principal Curvature" ); - typedef LpMetric L2Metric; - typedef functors::MongeJetFittingPrincipalCurvaturesEstimator > FunctorPrincCurv; - typedef functors::ConstValue< double > ConvFunctor; - typedef LocalEstimatorFromSurfelFunctorAdapter ReporterK; - CanonicSCellEmbedder embedder( K ); - FunctorPrincCurv estimatorK( embedder, h ); - ConvFunctor convFunc(1.0); - L2Metric l2(2.0); - ReporterK reporterK; - reporterK.attach( surf ); - reporterK.setParams( l2, estimatorK, convFunc, re/h ); - - c.startClock(); - - range = new VisitorRange( new Visitor( surf, *surf.begin() )); - ibegin = range->begin(); - iend = range->end(); - reporterK.init( h , ibegin, iend ); - - char full_filename[360]; - sprintf( full_filename, "%s%s", filename.c_str(), "_MongeJetFitting_principal.dat" ); - std::ofstream file( full_filename ); - file << "# h = " << h << std::endl; - file << "# Gaussian Curvature estimation from CGAL Monge from and Jet Fitting" << std::endl; - file << "# computed kernel radius = " << re << std::endl; - std::ostream_iterator< std::string > out_it_monge_principal( file, "\n" ); - - std::vector v_results; - std::back_insert_iterator > bkIt(v_results); - - delete range; - range = new VisitorRange( new Visitor( surf, *surf.begin() )); - ibegin = range->begin(); - iend = range->end(); - reporterK.eval(ibegin, iend , bkIt);//out_it_monge_principal); - - for(unsigned int ii = 0; ii < v_results.size(); ++ii ) - { - std::stringstream ss; - ss << v_results[ii].first << " " << v_results[ii].second; - *out_it_monge_principal = ss.str(); - ++out_it_monge_principal; - } - - double TMongeGaussCurv = c.stopClock(); - file << "# time = " << TMongeGaussCurv << std::endl; - file.close(); - delete range; - - - trace.endBlock(); - } - } + trace.beginBlock( "True mean curvature" ); + + char full_filename[360]; + sprintf( full_filename, "%s%s", filename.c_str(), "_True_mean.dat" ); + std::ofstream file( full_filename ); + file << "# h = " << h << std::endl; + file << "# True Mean Curvature estimation" << std::endl; + + std::ostream_iterator< double > out_it_true_mean( file, "\n" ); + + range = new VisitorRange( new Visitor( surf, *surf.begin() )); + ibegin = range->begin(); + iend = range->end(); + + c.startClock(); + + estimateTrueMeanCurvatureQuantity( ibegin, + iend, + out_it_true_mean, + K, + h, + aShape ); + + double TTrueMeanCurv = c.stopClock(); + file << "# time = " << TTrueMeanCurv << std::endl; + + file.close(); + delete range; + + trace.endBlock(); } - else // no noise + + // True Gaussian Curvature + if( properties.at( 1 ) != '0' ) { - typedef LightImplicitDigitalSurface< KSpace, DigitalShape > Boundary; - typedef DigitalSurface< Boundary > MyDigitalSurface; - typedef typename MyDigitalSurface::ConstIterator ConstIterator; - - typedef DepthFirstVisitor< MyDigitalSurface > Visitor; - typedef GraphVisitorRange< Visitor > VisitorRange; - typedef typename VisitorRange::ConstIterator VisitorConstIterator; - - // Extracts shape boundary - SCell bel = Surfaces::findABel ( K, dshape, 10000 ); - Boundary boundary( K, dshape, SurfelAdjacency< KSpace::dimension >( true ), bel ); - MyDigitalSurface surf ( boundary ); - - VisitorRange * range; - VisitorConstIterator ibegin; - VisitorConstIterator iend; - - unsigned int cntIn = 0; - for( typename Z3i::Domain::ConstIterator it = dshape.getDomain().begin(), ite = dshape.getDomain().end(); it != ite; ++it ) - { - if( dshape.operator ()(*it)) - { - ++cntIn; - } - } - - std::cout << "boundary:" << surf.size() << std::endl; - std::cout << "full:" << cntIn << std::endl; - - // Estimations - Clock c; - - // True - if( options.at( 0 ) != '0' ) - { - // True Mean Curvature - if( properties.at( 0 ) != '0' ) - { - trace.beginBlock( "True mean curvature" ); - - char full_filename[360]; - sprintf( full_filename, "%s%s", filename.c_str(), "_True_mean.dat" ); - std::ofstream file( full_filename ); - file << "# h = " << h << std::endl; - file << "# True Mean Curvature estimation" << std::endl; - - std::ostream_iterator< double > out_it_true_mean( file, "\n" ); - - range = new VisitorRange( new Visitor( surf, *surf.begin() )); - ibegin = range->begin(); - iend = range->end(); - - c.startClock(); - - estimateTrueMeanCurvatureQuantity( ibegin, - iend, - out_it_true_mean, - K, - h, - aShape ); - - double TTrueMeanCurv = c.stopClock(); - file << "# time = " << TTrueMeanCurv << std::endl; - - file.close(); - delete range; - - trace.endBlock(); - } - - // True Gaussian Curvature - if( properties.at( 1 ) != '0' ) - { - trace.beginBlock( "True Gaussian curvature" ); - - char full_filename[360]; - sprintf( full_filename, "%s%s", filename.c_str(), "_True_gaussian.dat" ); - std::ofstream file( full_filename ); - file << "# h = " << h << std::endl; - file << "# True Gaussian Curvature estimation" << std::endl; - - std::ostream_iterator< double > out_it_true_gaussian( file, "\n" ); - - range = new VisitorRange( new Visitor( surf, *surf.begin() )); - ibegin = range->begin(); - iend = range->end(); - - c.startClock(); - - estimateTrueGaussianCurvatureQuantity( ibegin, - iend, - out_it_true_gaussian, - K, - h, - aShape ); - - double TTrueGaussianCurv = c.stopClock(); - file << "# time = " << TTrueGaussianCurv << std::endl; - - file.close(); - - delete range; - - trace.endBlock(); - } - - // True Principal Curvatures - if( properties.at( 2 ) != '0' ) - { - trace.beginBlock( "True principal curvatures" ); - - char full_filename[360]; - sprintf( full_filename, "%s%s", filename.c_str(), "_True_principal.dat" ); - std::ofstream file( full_filename ); - file << "# h = " << h << std::endl; - file << "# True Gaussian Curvature estimation" << std::endl; - - std::ostream_iterator< std::string > out_it_true_pc( file, "\n" ); - - std::vector v_results; - std::back_insert_iterator< std::vector > bkIt(v_results); - - range = new VisitorRange( new Visitor( surf, *surf.begin() )); - ibegin = range->begin(); - iend = range->end(); - - c.startClock(); - - estimateTruePrincipalCurvaturesQuantity( ibegin, - iend, - bkIt,// out_it_true_pc, - K, - h, - aShape ); - - - for(unsigned int ii = 0; ii < v_results.size(); ++ii ) - { - std::stringstream ss; - ss << v_results[ii].first << " " << v_results[ii].second; - *out_it_true_pc = ss.str(); - ++out_it_true_pc; - } - - double TTruePrincCurv = c.stopClock(); - file << "# time = " << TTruePrincCurv << std::endl; - - file.close(); - - delete range; - - trace.endBlock(); - } - } - - double re = radius_kernel * std::pow( h, alpha ); // to obtains convergence results, re must follow the rule re=kh^(1/3) - - // II - if( options.at( 1 ) != '0' ) - { - // Integral Invariant Mean Curvature - if( properties.at( 0 ) != '0' ) - { - trace.beginBlock( "II mean curvature" ); - - char full_filename[360]; - sprintf( full_filename, "%s%s", filename.c_str(), "_II_mean.dat" ); - std::ofstream file( full_filename ); - file << "# h = " << h << std::endl; - file << "# Mean Curvature estimation from the Integral Invariant" << std::endl; - file << "# computed kernel radius = " << re << std::endl; - - range = new VisitorRange( new Visitor( surf, *surf.begin() )); - ibegin = range->begin(); - iend = range->end(); - - c.startClock(); - - typedef functors::IIMeanCurvature3DFunctor MyIICurvatureFunctor; - typedef IntegralInvariantVolumeEstimator< KSpace, DigitalShape, MyIICurvatureFunctor > MyIICurvatureEstimator; - - MyIICurvatureFunctor curvatureFunctor; - curvatureFunctor.init( h, re ); - - MyIICurvatureEstimator curvatureEstimator( curvatureFunctor ); - curvatureEstimator.attach( K, dshape ); - curvatureEstimator.setParams( re/h ); - curvatureEstimator.init( h, ibegin, iend ); - - std::ostream_iterator< double > out_it_ii_mean( file, "\n" ); - curvatureEstimator.eval( ibegin, iend, out_it_ii_mean ); - - double TIIMeanCurv = c.stopClock(); - file << "# time = " << TIIMeanCurv << std::endl; - file.close(); - - delete range; - - trace.endBlock(); - } - - // Integral Invariant Gaussian Curvature - if( properties.at( 1 ) != '0' ) - { - trace.beginBlock( "II Gaussian curvature" ); - - char full_filename[360]; - sprintf( full_filename, "%s%s", filename.c_str(), "_II_gaussian.dat" ); - std::ofstream file( full_filename ); - file << "# h = " << h << std::endl; - file << "# Gaussian Curvature estimation from the Integral Invariant" << std::endl; - file << "# computed kernel radius = " << re << std::endl; - - range = new VisitorRange( new Visitor( surf, *surf.begin() )); - ibegin = range->begin(); - iend = range->end(); - - c.startClock(); - - typedef functors::IIGaussianCurvature3DFunctor MyIICurvatureFunctor; - typedef IntegralInvariantCovarianceEstimator< KSpace, DigitalShape, MyIICurvatureFunctor > MyIICurvatureEstimator; - - MyIICurvatureFunctor curvatureFunctor; - curvatureFunctor.init( h, re ); - - MyIICurvatureEstimator curvatureEstimator( curvatureFunctor ); - curvatureEstimator.attach( K, dshape ); - curvatureEstimator.setParams( re/h ); - curvatureEstimator.init( h, ibegin, iend ); - - std::ostream_iterator< double > out_it_ii_gaussian( file, "\n" ); - curvatureEstimator.eval( ibegin, iend, out_it_ii_gaussian ); - - double TIIGaussCurv = c.stopClock(); - file << "# time = " << TIIGaussCurv << std::endl; - file.close(); - - delete range; - - trace.endBlock(); - } - - // Integral Invariant Principal Curvatures - if( properties.at( 2 ) != '0' ) - { - trace.beginBlock( "II Principal curvatures" ); - - char full_filename[360]; - sprintf( full_filename, "%s%s", filename.c_str(), "_II_principal.dat" ); - std::ofstream file( full_filename ); - file << "# h = " << h << std::endl; - file << "# Gaussian Curvature estimation from the Integral Invariant" << std::endl; - file << "# computed kernel radius = " << re << std::endl; - - range = new VisitorRange( new Visitor( surf, *surf.begin() )); - ibegin = range->begin(); - iend = range->end(); - - c.startClock(); - - typedef functors::IIPrincipalCurvatures3DFunctor MyIICurvatureFunctor; - typedef IntegralInvariantCovarianceEstimator< KSpace, DigitalShape, MyIICurvatureFunctor > MyIICurvatureEstimator; - - MyIICurvatureFunctor curvatureFunctor; - curvatureFunctor.init( h, re ); - - MyIICurvatureEstimator curvatureEstimator( curvatureFunctor ); - curvatureEstimator.attach( K, dshape ); - curvatureEstimator.setParams( re/h ); - curvatureEstimator.init( h, ibegin, iend ); - - std::vector v_results; - std::back_insert_iterator< std::vector > bkIt(v_results); - - curvatureEstimator.eval( ibegin, iend, bkIt ); - - std::ostream_iterator< std::string > out_it_ii_principal( file, "\n" ); - for( unsigned int ii = 0; ii < v_results.size(); ++ii ) - { - std::stringstream ss; - ss << v_results[ii].first << " " << v_results[ii].second; - *out_it_ii_principal = ss.str(); - ++out_it_ii_principal; - } - - double TIIGaussCurv = c.stopClock(); - file << "# time = " << TIIGaussCurv << std::endl; - file.close(); - - delete range; - - trace.endBlock(); - } - } - - // Monge - if( options.at( 2 ) != '0' ) - { - // Monge Mean Curvature - if( properties.at( 0 ) != '0' ) - { - trace.beginBlock( "Monge mean curvature" ); - typedef LpMetric L2Metric; - typedef functors::MongeJetFittingMeanCurvatureEstimator > FunctorMean; - typedef functors::ConstValue< double > ConvFunctor; - typedef LocalEstimatorFromSurfelFunctorAdapter ReporterH; - CanonicSCellEmbedder embedder( K ); - FunctorMean estimatorH( embedder, h ); - ConvFunctor convFunc(1.0); - L2Metric l2(2.0); - ReporterH reporterH; - reporterH.attach( surf ); - reporterH.setParams( l2, estimatorH, convFunc, re/h ); - c.startClock(); - - range = new VisitorRange( new Visitor( surf, *surf.begin() )); - ibegin = range->begin(); - iend = range->end(); - reporterH.init( h , ibegin, iend ); - - char full_filename[360]; - sprintf( full_filename, "%s%s", filename.c_str(), "_MongeJetFitting_mean.dat" ); - std::ofstream file( full_filename ); - file << "# h = " << h << std::endl; - file << "# Mean Curvature estimation from CGAL Monge from and Jet Fitting" << std::endl; - file << "# computed kernel radius = " << re << std::endl; - std::ostream_iterator< double > out_it_monge_mean( file, "\n" ); - - delete range; - range = new VisitorRange( new Visitor( surf, *surf.begin() )); - ibegin = range->begin(); - iend = range->end(); - - reporterH.eval(ibegin, iend , out_it_monge_mean); - double TMongeMeanCurv = c.stopClock(); - file << "# time = " << TMongeMeanCurv << std::endl; - file.close(); - delete range; - - - trace.endBlock(); - } - - // Monge Gaussian Curvature - if( properties.at( 1 ) != '0' ) - { - trace.beginBlock( "Monge Gaussian curvature" ); - typedef LpMetric L2Metric; - typedef functors::MongeJetFittingGaussianCurvatureEstimator > FunctorGaussian; - typedef functors::ConstValue< double > ConvFunctor; - typedef LocalEstimatorFromSurfelFunctorAdapter ReporterK; - CanonicSCellEmbedder embedder( K ); - FunctorGaussian estimatorK( embedder, h ); - ConvFunctor convFunc(1.0); - ReporterK reporterK; - reporterK.attach( surf ); - L2Metric l2(2.0); - reporterK.setParams( l2, estimatorK, convFunc, re/h ); - - c.startClock(); - - range = new VisitorRange( new Visitor( surf, *surf.begin() )); - ibegin = range->begin(); - iend = range->end(); - reporterK.init( h , ibegin, iend ); - - delete range; - range = new VisitorRange( new Visitor( surf, *surf.begin() )); - ibegin = range->begin(); - iend = range->end(); - - char full_filename[360]; - sprintf( full_filename, "%s%s", filename.c_str(), "_MongeJetFitting_gaussian.dat" ); - std::ofstream file( full_filename ); - file << "# h = " << h << std::endl; - file << "# Gaussian Curvature estimation from CGAL Monge from and Jet Fitting" << std::endl; - file << "# computed kernel radius = " << re << std::endl; - std::ostream_iterator< double > out_it_monge_gaussian( file, "\n" ); - reporterK.eval(ibegin, iend , out_it_monge_gaussian); - double TMongeGaussCurv = c.stopClock(); - file << "# time = " << TMongeGaussCurv << std::endl; - file.close(); - delete range; - - - trace.endBlock(); - } - - // Monge Principal Curvatures - if( properties.at( 2 ) != '0' ) - { - trace.beginBlock( "Monge Principal Curvature" ); - typedef LpMetric L2Metric; - - typedef functors::MongeJetFittingPrincipalCurvaturesEstimator > FunctorPrincCurv; - typedef functors::ConstValue< double > ConvFunctor; - typedef LocalEstimatorFromSurfelFunctorAdapter ReporterK; - CanonicSCellEmbedder embedder( K ); - FunctorPrincCurv estimatorK( embedder, h ); - ConvFunctor convFunc(1.0); - ReporterK reporterK; - reporterK.attach(surf); - L2Metric l2(2.0); - reporterK.setParams(l2, estimatorK, convFunc, re/h); - c.startClock(); - range = new VisitorRange( new Visitor( surf, *surf.begin() )); - ibegin = range->begin(); - iend = range->end(); - reporterK.init( h , ibegin, iend ); - - delete range; - range = new VisitorRange( new Visitor( surf, *surf.begin() )); - ibegin = range->begin(); - iend = range->end(); - - char full_filename[360]; - sprintf( full_filename, "%s%s", filename.c_str(), "_MongeJetFitting_principal.dat" ); - std::ofstream file( full_filename ); - file << "# h = " << h << std::endl; - file << "# Gaussian Curvature estimation from CGAL Monge from and Jet Fitting" << std::endl; - file << "# computed kernel radius = " << re << std::endl; - std::ostream_iterator< std::string > out_it_monge_principal( file, "\n" ); - - std::vector v_results; - std::back_insert_iterator > bkIt(v_results); - reporterK.eval(ibegin, iend , bkIt);//out_it_monge_principal); - - for(unsigned int ii = 0; ii < v_results.size(); ++ii ) - { - std::stringstream ss; - ss << v_results[ii].first << " " << v_results[ii].second; - *out_it_monge_principal = ss.str(); - ++out_it_monge_principal; - } - - double TMongeGaussCurv = c.stopClock(); - file << "# time = " << TMongeGaussCurv << std::endl; - file.close(); - delete range; - - - trace.endBlock(); - } - - } + trace.beginBlock( "True Gausian curvature" ); + + char full_filename[360]; + sprintf( full_filename, "%s%s", filename.c_str(), "_True_gaussian.dat" ); + std::ofstream file( full_filename ); + file << "# h = " << h << std::endl; + file << "# True Gaussian Curvature estimation" << std::endl; + + std::ostream_iterator< double > out_it_true_gaussian( file, "\n" ); + + range = new VisitorRange( new Visitor( surf, *surf.begin() )); + ibegin = range->begin(); + iend = range->end(); + + c.startClock(); + + estimateTrueGaussianCurvatureQuantity( ibegin, + iend, + out_it_true_gaussian, + K, + h, + aShape ); + + double TTrueGaussianCurv = c.stopClock(); + file << "# time = " << TTrueGaussianCurv << std::endl; + + file.close(); + delete range; + + trace.endBlock(); + } + + // True Principal Curvatures + if( properties.at( 2 ) != '0' ) + { + trace.beginBlock( "True principal curvatures" ); + + char full_filename[360]; + sprintf( full_filename, "%s%s", filename.c_str(), "_True_principal.dat" ); + std::ofstream file( full_filename ); + file << "# h = " << h << std::endl; + file << "# True Gaussian Curvature estimation" << std::endl; + + std::ostream_iterator< std::string > out_it_true_pc( file, "\n" ); + + std::vector v_results; + std::back_insert_iterator< std::vector > bkIt(v_results); + range = new VisitorRange( new Visitor( surf, *surf.begin() )); + ibegin = range->begin(); + iend = range->end(); + + c.startClock(); + + estimateTruePrincipalCurvaturesQuantity( ibegin, + iend, + bkIt,//out_it_true_pc, + K, + h, + aShape ); + + for(unsigned int ii = 0; ii < v_results.size(); ++ii ) + { + std::stringstream ss; + ss << v_results[ii].first << " " << v_results[ii].second; + *out_it_true_pc = ss.str(); + ++out_it_true_pc; + } + double TTruePrincCurv = c.stopClock(); + file << "# time = " << TTruePrincCurv << std::endl; + + file.close(); + + delete range; + + trace.endBlock(); + } + } + + double re = radius_kernel * std::pow( h, alpha ); // to obtains convergence results, re must follow the rule re=kh^(1/3) + + // II + if( options.at( 1 ) != '0' ) + { + // Integral Invariant Mean Curvature + if( properties.at( 0 ) != '0' ) + { + trace.beginBlock( "II mean curvature" ); + + char full_filename[360]; + sprintf( full_filename, "%s%s", filename.c_str(), "_II_mean.dat" ); + std::ofstream file( full_filename ); + file << "# h = " << h << std::endl; + file << "# Mean Curvature estimation from the Integral Invariant" << std::endl; + file << "# computed kernel radius = " << re << std::endl; + + range = new VisitorRange( new Visitor( surf, *surf.begin() )); + ibegin = range->begin(); + iend = range->end(); + + c.startClock(); + + typedef functors::IIMeanCurvature3DFunctor MyIICurvatureFunctor; + typedef IntegralInvariantVolumeEstimator< KSpace, KanungoPredicate, MyIICurvatureFunctor > MyIICurvatureEstimator; + + MyIICurvatureFunctor curvatureFunctor; + curvatureFunctor.init( h, re ); + + MyIICurvatureEstimator curvatureEstimator( curvatureFunctor ); + curvatureEstimator.attach( K, noisifiedObject ); + curvatureEstimator.setParams( re/h ); + curvatureEstimator.init( h, ibegin, iend ); + + std::ostream_iterator< double > out_it_ii_mean( file, "\n" ); + curvatureEstimator.eval( ibegin, iend, out_it_ii_mean ); + + double TIIMeanCurv = c.stopClock(); + file << "# time = " << TIIMeanCurv << std::endl; + file.close(); + + delete range; + + trace.endBlock(); + } + + // Integral Invariant Gaussian Curvature + if( properties.at( 1 ) != '0' ) + { + trace.beginBlock( "II Gaussian curvature" ); + + char full_filename[360]; + sprintf( full_filename, "%s%s", filename.c_str(), "_II_gaussian.dat" ); + std::ofstream file( full_filename ); + file << "# h = " << h << std::endl; + file << "# Gaussian Curvature estimation from the Integral Invariant" << std::endl; + file << "# computed kernel radius = " << re << std::endl; + + range = new VisitorRange( new Visitor( surf, *surf.begin() )); + ibegin = range->begin(); + iend = range->end(); + + c.startClock(); + + typedef functors::IIGaussianCurvature3DFunctor MyIICurvatureFunctor; + typedef IntegralInvariantCovarianceEstimator< KSpace, KanungoPredicate, MyIICurvatureFunctor > MyIICurvatureEstimator; + + MyIICurvatureFunctor curvatureFunctor; + curvatureFunctor.init( h, re ); + + MyIICurvatureEstimator curvatureEstimator( curvatureFunctor ); + curvatureEstimator.attach( K, noisifiedObject ); + curvatureEstimator.setParams( re/h ); + curvatureEstimator.init( h, ibegin, iend ); + + std::ostream_iterator< double > out_it_ii_gaussian( file, "\n" ); + curvatureEstimator.eval( ibegin, iend, out_it_ii_gaussian ); + + double TIIGaussCurv = c.stopClock(); + file << "# time = " << TIIGaussCurv << std::endl; + file.close(); + + delete range; + + trace.endBlock(); + } + + // Integral Invariant Principal Curvatures + if( properties.at( 2 ) != '0' ) + { + trace.beginBlock( "II Principal curvatures" ); + + char full_filename[360]; + sprintf( full_filename, "%s%s", filename.c_str(), "_II_principal.dat" ); + std::ofstream file( full_filename ); + file << "# h = " << h << std::endl; + file << "# Gaussian Curvature estimation from the Integral Invariant" << std::endl; + file << "# computed kernel radius = " << re << std::endl; + + range = new VisitorRange( new Visitor( surf, *surf.begin() )); + ibegin = range->begin(); + iend = range->end(); + + c.startClock(); + + typedef functors::IIPrincipalCurvatures3DFunctor MyIICurvatureFunctor; + typedef IntegralInvariantCovarianceEstimator< KSpace, KanungoPredicate, MyIICurvatureFunctor > MyIICurvatureEstimator; + + MyIICurvatureFunctor curvatureFunctor; + curvatureFunctor.init( h, re ); + + MyIICurvatureEstimator curvatureEstimator( curvatureFunctor ); + curvatureEstimator.attach( K, noisifiedObject ); + curvatureEstimator.setParams( re/h ); + curvatureEstimator.init( h, ibegin, iend ); + + std::vector v_results; + std::back_insert_iterator< std::vector > bkIt(v_results); + + curvatureEstimator.eval( ibegin, iend, bkIt ); + + std::ostream_iterator< std::string > out_it_ii_principal( file, "\n" ); + for( unsigned int ii = 0; ii < v_results.size(); ++ii ) + { + std::stringstream ss; + ss << v_results[ii].first << " " << v_results[ii].second; + *out_it_ii_principal = ss.str(); + ++out_it_ii_principal; + } + + double TIIGaussCurv = c.stopClock(); + file << "# time = " << TIIGaussCurv << std::endl; + file.close(); + + delete range; + + trace.endBlock(); } + } + + // Monge + if( options.at( 2 ) != '0' ) + { + // Monge Mean Curvature + if( properties.at( 0 ) != '0' ) + { + trace.beginBlock( "Monge mean curvature" ); + typedef LpMetric L2Metric; + typedef functors::MongeJetFittingMeanCurvatureEstimator > FunctorMean; + typedef functors::ConstValue< double > ConvFunctor; + typedef LocalEstimatorFromSurfelFunctorAdapter ReporterH; + CanonicSCellEmbedder embedder( K ); + FunctorMean estimatorH( embedder, h ); + ConvFunctor convFunc(1.0); + L2Metric l2(2.0); + ReporterH reporterH; + reporterH.attach( surf ); + reporterH.setParams( l2, estimatorH, convFunc, re/h ); + + c.startClock(); + range = new VisitorRange( new Visitor( surf, *surf.begin() )); + ibegin = range->begin(); + iend = range->end(); + + reporterH.init( h , ibegin, iend ); + + char full_filename[360]; + sprintf( full_filename, "%s%s", filename.c_str(), "_MongeJetFitting_mean.dat" ); + std::ofstream file( full_filename ); + file << "# h = " << h << std::endl; + file << "# Mean Curvature estimation from CGAL Monge from and Jet Fitting" << std::endl; + file << "# computed kernel radius = " << re << std::endl; + std::ostream_iterator< double > out_it_monge_mean( file, "\n" ); + + delete range; + range = new VisitorRange( new Visitor( surf, *surf.begin() )); + ibegin = range->begin(); + iend = range->end(); + reporterH.eval(ibegin, iend, out_it_monge_mean); + double TMongeMeanCurv = c.stopClock(); + file << "# time = " << TMongeMeanCurv << std::endl; + file.close(); + delete range; + + + trace.endBlock(); + } + + // Monge Gaussian Curvature + if( properties.at( 1 ) != '0' ) + { + trace.beginBlock( "Monge Gaussian curvature" ); + + typedef functors::MongeJetFittingGaussianCurvatureEstimator > FunctorGaussian; + typedef functors::ConstValue< double > ConvFunctor; + typedef LpMetric L2Metric; + typedef LocalEstimatorFromSurfelFunctorAdapter ReporterK; + CanonicSCellEmbedder embedder( K ); + FunctorGaussian estimatorK( embedder, h ); + ConvFunctor convFunc(1.0); + ReporterK reporterK; + reporterK.attach( surf ); + L2Metric l2(2.0); + reporterK.setParams( l2, estimatorK, convFunc, re/h ); + c.startClock(); + + range = new VisitorRange( new Visitor( surf, *surf.begin() )); + ibegin = range->begin(); + iend = range->end(); + + reporterK.init( h , ibegin, iend ); + + delete range; + range = new VisitorRange( new Visitor( surf, *surf.begin() )); + ibegin = range->begin(); + iend = range->end(); + + char full_filename[360]; + sprintf( full_filename, "%s%s", filename.c_str(), "_MongeJetFitting_gaussian.dat" ); + std::ofstream file( full_filename ); + file << "# h = " << h << std::endl; + file << "# Gaussian Curvature estimation from CGAL Monge from and Jet Fitting" << std::endl; + file << "# computed kernel radius = " << re << std::endl; + std::ostream_iterator< double > out_it_monge_gaussian( file, "\n" ); + reporterK.eval(ibegin, iend , out_it_monge_gaussian); + double TMongeGaussCurv = c.stopClock(); + file << "# time = " << TMongeGaussCurv << std::endl; + file.close(); + delete range; + + + trace.endBlock(); + } + + // Monge Principal Curvatures + if( properties.at( 2 ) != '0' ) + { + trace.beginBlock( "Monge Principal Curvature" ); + typedef LpMetric L2Metric; + typedef functors::MongeJetFittingPrincipalCurvaturesEstimator > FunctorPrincCurv; + typedef functors::ConstValue< double > ConvFunctor; + typedef LocalEstimatorFromSurfelFunctorAdapter ReporterK; + CanonicSCellEmbedder embedder( K ); + FunctorPrincCurv estimatorK( embedder, h ); + ConvFunctor convFunc(1.0); + L2Metric l2(2.0); + ReporterK reporterK; + reporterK.attach( surf ); + reporterK.setParams( l2, estimatorK, convFunc, re/h ); + + c.startClock(); + + range = new VisitorRange( new Visitor( surf, *surf.begin() )); + ibegin = range->begin(); + iend = range->end(); + reporterK.init( h , ibegin, iend ); + + char full_filename[360]; + sprintf( full_filename, "%s%s", filename.c_str(), "_MongeJetFitting_principal.dat" ); + std::ofstream file( full_filename ); + file << "# h = " << h << std::endl; + file << "# Gaussian Curvature estimation from CGAL Monge from and Jet Fitting" << std::endl; + file << "# computed kernel radius = " << re << std::endl; + std::ostream_iterator< std::string > out_it_monge_principal( file, "\n" ); + + std::vector v_results; + std::back_insert_iterator > bkIt(v_results); + + delete range; + range = new VisitorRange( new Visitor( surf, *surf.begin() )); + ibegin = range->begin(); + iend = range->end(); + reporterK.eval(ibegin, iend , bkIt);//out_it_monge_principal); + + for(unsigned int ii = 0; ii < v_results.size(); ++ii ) + { + std::stringstream ss; + ss << v_results[ii].first << " " << v_results[ii].second; + *out_it_monge_principal = ss.str(); + ++out_it_monge_principal; + } + + double TMongeGaussCurv = c.stopClock(); + file << "# time = " << TMongeGaussCurv << std::endl; + file.close(); + delete range; + + + trace.endBlock(); + } + } } - catch ( InputException e ) + else // no noise { - std::cerr << "[estimatorCurvatureComparator3D]" - << " error." - << e.what() - << " " - << filename << " " - << border_min[0] << " " - << border_max[0] << " " - << h << " " - << radius_kernel << " " - << lambda_optimized << " " - << options << " " - << alpha << " " - << std::endl; - return false; + typedef LightImplicitDigitalSurface< KSpace, DigitalShape > Boundary; + typedef DigitalSurface< Boundary > MyDigitalSurface; + typedef typename MyDigitalSurface::ConstIterator ConstIterator; + + typedef DepthFirstVisitor< MyDigitalSurface > Visitor; + typedef GraphVisitorRange< Visitor > VisitorRange; + typedef typename VisitorRange::ConstIterator VisitorConstIterator; + + // Extracts shape boundary + SCell bel = Surfaces::findABel ( K, dshape, 10000 ); + Boundary boundary( K, dshape, SurfelAdjacency< KSpace::dimension >( true ), bel ); + MyDigitalSurface surf ( boundary ); + + VisitorRange * range; + VisitorConstIterator ibegin; + VisitorConstIterator iend; + + unsigned int cntIn = 0; + for( typename Z3i::Domain::ConstIterator it = dshape.getDomain().begin(), ite = dshape.getDomain().end(); it != ite; ++it ) + { + if( dshape.operator ()(*it)) + { + ++cntIn; + } + } + + std::cout << "boundary:" << surf.size() << std::endl; + std::cout << "full:" << cntIn << std::endl; + + // Estimations + Clock c; + + // True + if( options.at( 0 ) != '0' ) + { + // True Mean Curvature + if( properties.at( 0 ) != '0' ) + { + trace.beginBlock( "True mean curvature" ); + + char full_filename[360]; + sprintf( full_filename, "%s%s", filename.c_str(), "_True_mean.dat" ); + std::ofstream file( full_filename ); + file << "# h = " << h << std::endl; + file << "# True Mean Curvature estimation" << std::endl; + + std::ostream_iterator< double > out_it_true_mean( file, "\n" ); + + range = new VisitorRange( new Visitor( surf, *surf.begin() )); + ibegin = range->begin(); + iend = range->end(); + + c.startClock(); + + estimateTrueMeanCurvatureQuantity( ibegin, + iend, + out_it_true_mean, + K, + h, + aShape ); + + double TTrueMeanCurv = c.stopClock(); + file << "# time = " << TTrueMeanCurv << std::endl; + + file.close(); + delete range; + + trace.endBlock(); + } + + // True Gaussian Curvature + if( properties.at( 1 ) != '0' ) + { + trace.beginBlock( "True Gaussian curvature" ); + + char full_filename[360]; + sprintf( full_filename, "%s%s", filename.c_str(), "_True_gaussian.dat" ); + std::ofstream file( full_filename ); + file << "# h = " << h << std::endl; + file << "# True Gaussian Curvature estimation" << std::endl; + + std::ostream_iterator< double > out_it_true_gaussian( file, "\n" ); + + range = new VisitorRange( new Visitor( surf, *surf.begin() )); + ibegin = range->begin(); + iend = range->end(); + + c.startClock(); + + estimateTrueGaussianCurvatureQuantity( ibegin, + iend, + out_it_true_gaussian, + K, + h, + aShape ); + + double TTrueGaussianCurv = c.stopClock(); + file << "# time = " << TTrueGaussianCurv << std::endl; + + file.close(); + + delete range; + + trace.endBlock(); + } + + // True Principal Curvatures + if( properties.at( 2 ) != '0' ) + { + trace.beginBlock( "True principal curvatures" ); + + char full_filename[360]; + sprintf( full_filename, "%s%s", filename.c_str(), "_True_principal.dat" ); + std::ofstream file( full_filename ); + file << "# h = " << h << std::endl; + file << "# True Gaussian Curvature estimation" << std::endl; + + std::ostream_iterator< std::string > out_it_true_pc( file, "\n" ); + + std::vector v_results; + std::back_insert_iterator< std::vector > bkIt(v_results); + + range = new VisitorRange( new Visitor( surf, *surf.begin() )); + ibegin = range->begin(); + iend = range->end(); + + c.startClock(); + + estimateTruePrincipalCurvaturesQuantity( ibegin, + iend, + bkIt,// out_it_true_pc, + K, + h, + aShape ); + + + for(unsigned int ii = 0; ii < v_results.size(); ++ii ) + { + std::stringstream ss; + ss << v_results[ii].first << " " << v_results[ii].second; + *out_it_true_pc = ss.str(); + ++out_it_true_pc; + } + + double TTruePrincCurv = c.stopClock(); + file << "# time = " << TTruePrincCurv << std::endl; + + file.close(); + + delete range; + + trace.endBlock(); + } + } + + double re = radius_kernel * std::pow( h, alpha ); // to obtains convergence results, re must follow the rule re=kh^(1/3) + + // II + if( options.at( 1 ) != '0' ) + { + // Integral Invariant Mean Curvature + if( properties.at( 0 ) != '0' ) + { + trace.beginBlock( "II mean curvature" ); + char full_filename[360]; + sprintf( full_filename, "%s%s", filename.c_str(), "_II_mean.dat" ); + std::ofstream file( full_filename ); + file << "# h = " << h << std::endl; + file << "# Mean Curvature estimation from the Integral Invariant" << std::endl; + file << "# computed kernel radius = " << re << std::endl; + + range = new VisitorRange( new Visitor( surf, *surf.begin() )); + ibegin = range->begin(); + iend = range->end(); + + c.startClock(); + + typedef functors::IIMeanCurvature3DFunctor MyIICurvatureFunctor; + typedef IntegralInvariantVolumeEstimator< KSpace, DigitalShape, MyIICurvatureFunctor > MyIICurvatureEstimator; + + MyIICurvatureFunctor curvatureFunctor; + curvatureFunctor.init( h, re ); + + MyIICurvatureEstimator curvatureEstimator( curvatureFunctor ); + curvatureEstimator.attach( K, dshape ); + curvatureEstimator.setParams( re/h ); + curvatureEstimator.init( h, ibegin, iend ); + + std::ostream_iterator< double > out_it_ii_mean( file, "\n" ); + curvatureEstimator.eval( ibegin, iend, out_it_ii_mean ); + + double TIIMeanCurv = c.stopClock(); + file << "# time = " << TIIMeanCurv << std::endl; + file.close(); + + delete range; + + trace.endBlock(); + } + + // Integral Invariant Gaussian Curvature + if( properties.at( 1 ) != '0' ) + { + trace.beginBlock( "II Gaussian curvature" ); + + char full_filename[360]; + sprintf( full_filename, "%s%s", filename.c_str(), "_II_gaussian.dat" ); + std::ofstream file( full_filename ); + file << "# h = " << h << std::endl; + file << "# Gaussian Curvature estimation from the Integral Invariant" << std::endl; + file << "# computed kernel radius = " << re << std::endl; + + range = new VisitorRange( new Visitor( surf, *surf.begin() )); + ibegin = range->begin(); + iend = range->end(); + + c.startClock(); + + typedef functors::IIGaussianCurvature3DFunctor MyIICurvatureFunctor; + typedef IntegralInvariantCovarianceEstimator< KSpace, DigitalShape, MyIICurvatureFunctor > MyIICurvatureEstimator; + + MyIICurvatureFunctor curvatureFunctor; + curvatureFunctor.init( h, re ); + + MyIICurvatureEstimator curvatureEstimator( curvatureFunctor ); + curvatureEstimator.attach( K, dshape ); + curvatureEstimator.setParams( re/h ); + curvatureEstimator.init( h, ibegin, iend ); + + std::ostream_iterator< double > out_it_ii_gaussian( file, "\n" ); + curvatureEstimator.eval( ibegin, iend, out_it_ii_gaussian ); + + double TIIGaussCurv = c.stopClock(); + file << "# time = " << TIIGaussCurv << std::endl; + file.close(); + + delete range; + + trace.endBlock(); + } + + // Integral Invariant Principal Curvatures + if( properties.at( 2 ) != '0' ) + { + trace.beginBlock( "II Principal curvatures" ); + + char full_filename[360]; + sprintf( full_filename, "%s%s", filename.c_str(), "_II_principal.dat" ); + std::ofstream file( full_filename ); + file << "# h = " << h << std::endl; + file << "# Gaussian Curvature estimation from the Integral Invariant" << std::endl; + file << "# computed kernel radius = " << re << std::endl; + + range = new VisitorRange( new Visitor( surf, *surf.begin() )); + ibegin = range->begin(); + iend = range->end(); + + c.startClock(); + + typedef functors::IIPrincipalCurvatures3DFunctor MyIICurvatureFunctor; + typedef IntegralInvariantCovarianceEstimator< KSpace, DigitalShape, MyIICurvatureFunctor > MyIICurvatureEstimator; + + MyIICurvatureFunctor curvatureFunctor; + curvatureFunctor.init( h, re ); + + MyIICurvatureEstimator curvatureEstimator( curvatureFunctor ); + curvatureEstimator.attach( K, dshape ); + curvatureEstimator.setParams( re/h ); + curvatureEstimator.init( h, ibegin, iend ); + + std::vector v_results; + std::back_insert_iterator< std::vector > bkIt(v_results); + + curvatureEstimator.eval( ibegin, iend, bkIt ); + + std::ostream_iterator< std::string > out_it_ii_principal( file, "\n" ); + for( unsigned int ii = 0; ii < v_results.size(); ++ii ) + { + std::stringstream ss; + ss << v_results[ii].first << " " << v_results[ii].second; + *out_it_ii_principal = ss.str(); + ++out_it_ii_principal; + } + + double TIIGaussCurv = c.stopClock(); + file << "# time = " << TIIGaussCurv << std::endl; + file.close(); + + delete range; + + trace.endBlock(); + } + } + + // Monge + if( options.at( 2 ) != '0' ) + { + // Monge Mean Curvature + if( properties.at( 0 ) != '0' ) + { + trace.beginBlock( "Monge mean curvature" ); + typedef LpMetric L2Metric; + typedef functors::MongeJetFittingMeanCurvatureEstimator > FunctorMean; + typedef functors::ConstValue< double > ConvFunctor; + typedef LocalEstimatorFromSurfelFunctorAdapter ReporterH; + CanonicSCellEmbedder embedder( K ); + FunctorMean estimatorH( embedder, h ); + ConvFunctor convFunc(1.0); + L2Metric l2(2.0); + ReporterH reporterH; + reporterH.attach( surf ); + reporterH.setParams( l2, estimatorH, convFunc, re/h ); + c.startClock(); + + range = new VisitorRange( new Visitor( surf, *surf.begin() )); + ibegin = range->begin(); + iend = range->end(); + reporterH.init( h , ibegin, iend ); + + char full_filename[360]; + sprintf( full_filename, "%s%s", filename.c_str(), "_MongeJetFitting_mean.dat" ); + std::ofstream file( full_filename ); + file << "# h = " << h << std::endl; + file << "# Mean Curvature estimation from CGAL Monge from and Jet Fitting" << std::endl; + file << "# computed kernel radius = " << re << std::endl; + std::ostream_iterator< double > out_it_monge_mean( file, "\n" ); + + delete range; + range = new VisitorRange( new Visitor( surf, *surf.begin() )); + ibegin = range->begin(); + iend = range->end(); + + reporterH.eval(ibegin, iend , out_it_monge_mean); + double TMongeMeanCurv = c.stopClock(); + file << "# time = " << TMongeMeanCurv << std::endl; + file.close(); + delete range; + + + trace.endBlock(); + } + + // Monge Gaussian Curvature + if( properties.at( 1 ) != '0' ) + { + trace.beginBlock( "Monge Gaussian curvature" ); + typedef LpMetric L2Metric; + typedef functors::MongeJetFittingGaussianCurvatureEstimator > FunctorGaussian; + typedef functors::ConstValue< double > ConvFunctor; + typedef LocalEstimatorFromSurfelFunctorAdapter ReporterK; + CanonicSCellEmbedder embedder( K ); + FunctorGaussian estimatorK( embedder, h ); + ConvFunctor convFunc(1.0); + ReporterK reporterK; + reporterK.attach( surf ); + L2Metric l2(2.0); + reporterK.setParams( l2, estimatorK, convFunc, re/h ); + + c.startClock(); + + range = new VisitorRange( new Visitor( surf, *surf.begin() )); + ibegin = range->begin(); + iend = range->end(); + reporterK.init( h , ibegin, iend ); + + delete range; + range = new VisitorRange( new Visitor( surf, *surf.begin() )); + ibegin = range->begin(); + iend = range->end(); + + char full_filename[360]; + sprintf( full_filename, "%s%s", filename.c_str(), "_MongeJetFitting_gaussian.dat" ); + std::ofstream file( full_filename ); + file << "# h = " << h << std::endl; + file << "# Gaussian Curvature estimation from CGAL Monge from and Jet Fitting" << std::endl; + file << "# computed kernel radius = " << re << std::endl; + std::ostream_iterator< double > out_it_monge_gaussian( file, "\n" ); + reporterK.eval(ibegin, iend , out_it_monge_gaussian); + double TMongeGaussCurv = c.stopClock(); + file << "# time = " << TMongeGaussCurv << std::endl; + file.close(); + delete range; + + + trace.endBlock(); + } + + // Monge Principal Curvatures + if( properties.at( 2 ) != '0' ) + { + trace.beginBlock( "Monge Principal Curvature" ); + typedef LpMetric L2Metric; + + typedef functors::MongeJetFittingPrincipalCurvaturesEstimator > FunctorPrincCurv; + typedef functors::ConstValue< double > ConvFunctor; + typedef LocalEstimatorFromSurfelFunctorAdapter ReporterK; + CanonicSCellEmbedder embedder( K ); + FunctorPrincCurv estimatorK( embedder, h ); + ConvFunctor convFunc(1.0); + ReporterK reporterK; + reporterK.attach(surf); + L2Metric l2(2.0); + reporterK.setParams(l2, estimatorK, convFunc, re/h); + c.startClock(); + range = new VisitorRange( new Visitor( surf, *surf.begin() )); + ibegin = range->begin(); + iend = range->end(); + reporterK.init( h , ibegin, iend ); + + delete range; + range = new VisitorRange( new Visitor( surf, *surf.begin() )); + ibegin = range->begin(); + iend = range->end(); + + char full_filename[360]; + sprintf( full_filename, "%s%s", filename.c_str(), "_MongeJetFitting_principal.dat" ); + std::ofstream file( full_filename ); + file << "# h = " << h << std::endl; + file << "# Gaussian Curvature estimation from CGAL Monge from and Jet Fitting" << std::endl; + file << "# computed kernel radius = " << re << std::endl; + std::ostream_iterator< std::string > out_it_monge_principal( file, "\n" ); + + std::vector v_results; + std::back_insert_iterator > bkIt(v_results); + reporterK.eval(ibegin, iend , bkIt);//out_it_monge_principal); + + for(unsigned int ii = 0; ii < v_results.size(); ++ii ) + { + std::stringstream ss; + ss << v_results[ii].first << " " << v_results[ii].second; + *out_it_monge_principal = ss.str(); + ++out_it_monge_principal; + } + + double TMongeGaussCurv = c.stopClock(); + file << "# time = " << TMongeGaussCurv << std::endl; + file.close(); + delete range; + + + trace.endBlock(); + } + + } } - - - return true; + } + catch ( InputException e ) + { + std::cerr << "[estimatorCurvatureComparator3D]" + << " error." + << e.what() + << " " + << filename << " " + << border_min[0] << " " + << border_max[0] << " " + << h << " " + << radius_kernel << " " + << lambda_optimized << " " + << options << " " + << alpha << " " + << std::endl; + return false; + } + return true; } /** @@ -1120,147 +1117,147 @@ compareShapeEstimators( const std::string & filename, */ void missingParam( std::string param ) { - trace.error() << " Parameter: " << param << " is required."; - trace.info() << std::endl; - exit( 1 ); + trace.error() << " Parameter: " << param << " is required."; + trace.info() << std::endl; + exit( 1 ); } namespace po = boost::program_options; int main( int argc, char** argv ) { - - + + #ifndef WITH_CGAL #error You need to have activated CGAL (WITH_CGAL) to include this file. #endif #ifndef WITH_EIGEN #error You need to have activated EIGEN (WITH_EIGEN) to include this file. #endif - - // parse command line ---------------------------------------------- - po::options_description general_opt("Allowed options are"); - general_opt.add_options() - ("help,h", "display this message") - ("shape,s", po::value< std::string >(), "Shape") - ("output,o", po::value< std::string >(), "Output file") - ("radius,r", po::value< double >(), "Kernel radius for IntegralInvariant" ) - ("alpha", po::value()->default_value(1.0/3.0), "Alpha parameter for Integral Invariant computation" ) - ("h", po::value< double >(), "Grid step" ) - ("minAABB,a", po::value< double >()->default_value( -10.0 ), "Min value of the AABB bounding box (domain)" ) - ("maxAABB,A", po::value< double >()->default_value( 10.0 ), "Max value of the AABB bounding box (domain)" ) - ("noise,n", po::value()->default_value(0.0), "Level of noise to perturb the shape" ) - ("lambda,l", po::value< bool >()->default_value( false ), "Use the shape to get a better approximation of the surface (optional)" ) - ("properties", po::value()->default_value("110"), "the i-th property is disabled iff there is a 0 at position i" ) - ("estimators,e", po::value< std::string >()->default_value("110"), "the i-th estimator is disabled iff there is a 0 at position i" ); - - - bool parseOK = true; - po::variables_map vm; - try - { - po::store( po::parse_command_line( argc, argv, general_opt ), vm ); - } - catch( const std::exception & ex ) - { - parseOK = false; - trace.info() << "Error checking program options: " << ex.what() << std::endl; - } - po::notify( vm ); - if( !parseOK || vm.count("help") || argc <= 1 ) - { - trace.info()<< "Compare local estimators on implicit shapes using DGtal library" < --h --radius --estimators --output "<(); - int nb = 3; - std::string options = vm["estimators"].as< std::string >(); - if (options.size() < nb) - { - trace.error() << " At least " << nb - << " characters are required " - << " with option --estimators."; - trace.info() << std::endl; - exit(1); - } - double h = vm["h"].as< double >(); - double radius = vm["radius"].as< double >(); - double alpha = vm["alpha"].as< double >(); - std::string poly_str = vm["shape"].as< std::string >(); - bool lambda_optimized = vm["lambda"].as< bool >(); - double noiseLevel = vm["noise"].as(); - - nb = 3; //number of available properties - std::string properties = vm["properties"].as(); - if (properties.size() < nb) - { - trace.error() << " At least " << nb - << " characters are required " - << " with option --properties."; - trace.info() << std::endl; - exit(1); - } - - typedef Z3i::Space::RealPoint RealPoint; - typedef Z3i::Space::RealPoint::Coordinate Ring; - - RealPoint border_min( vm["minAABB"].as< double >(), vm["minAABB"].as< double >(), vm["minAABB"].as< double >() ); - RealPoint border_max( vm["maxAABB"].as< double >(), vm["maxAABB"].as< double >(), vm["maxAABB"].as< double >() ); - - /// Construction of the polynomial shape - - typedef MPolynomial< 3, Ring > Polynomial3; - typedef MPolynomialReader<3, Ring> Polynomial3Reader; - typedef ImplicitPolynomial3Shape ImplicitShape; - - Polynomial3 poly; - Polynomial3Reader reader; - std::string::const_iterator iter = reader.read( poly, poly_str.begin(), poly_str.end() ); - if ( iter != poly_str.end() ) - { - std::cerr << "ERROR: I read only <" - << poly_str.substr( 0, iter - poly_str.begin() ) - << ">, and I built P=" << poly << std::endl; - return 1; - } - - ImplicitShape* shape = new ImplicitShape( poly ); - - compareShapeEstimators< Z3i::Space, ImplicitShape > ( - file_export, - shape, - border_min, border_max, - h, - radius, - alpha, - options, - properties, - lambda_optimized, - noiseLevel ); - - delete shape; + + // parse command line ---------------------------------------------- + po::options_description general_opt("Allowed options are"); + general_opt.add_options() + ("help,h", "display this message") + ("shape,s", po::value< std::string >(), "Shape") + ("output,o", po::value< std::string >(), "Output file") + ("radius,r", po::value< double >(), "Kernel radius for IntegralInvariant" ) + ("alpha", po::value()->default_value(1.0/3.0), "Alpha parameter for Integral Invariant computation" ) + ("h", po::value< double >(), "Grid step" ) + ("minAABB,a", po::value< double >()->default_value( -10.0 ), "Min value of the AABB bounding box (domain)" ) + ("maxAABB,A", po::value< double >()->default_value( 10.0 ), "Max value of the AABB bounding box (domain)" ) + ("noise,n", po::value()->default_value(0.0), "Level of noise to perturb the shape" ) + ("lambda,l", po::value< bool >()->default_value( false ), "Use the shape to get a better approximation of the surface (optional)" ) + ("properties", po::value()->default_value("110"), "the i-th property is disabled iff there is a 0 at position i" ) + ("estimators,e", po::value< std::string >()->default_value("110"), "the i-th estimator is disabled iff there is a 0 at position i" ); + + + bool parseOK = true; + po::variables_map vm; + try + { + po::store( po::parse_command_line( argc, argv, general_opt ), vm ); + } + catch( const std::exception & ex ) + { + parseOK = false; + trace.info() << "Error checking program options: " << ex.what() << std::endl; + } + po::notify( vm ); + if( !parseOK || vm.count("help") || argc <= 1 ) + { + trace.info()<< "Compare local estimators on implicit shapes using DGtal library" < --h --radius --estimators --output "<(); + int nb = 3; + std::string options = vm["estimators"].as< std::string >(); + if (options.size() < nb) + { + trace.error() << " At least " << nb + << " characters are required " + << " with option --estimators."; + trace.info() << std::endl; + exit(1); + } + double h = vm["h"].as< double >(); + double radius = vm["radius"].as< double >(); + double alpha = vm["alpha"].as< double >(); + std::string poly_str = vm["shape"].as< std::string >(); + bool lambda_optimized = vm["lambda"].as< bool >(); + double noiseLevel = vm["noise"].as(); + + nb = 3; //number of available properties + std::string properties = vm["properties"].as(); + if (properties.size() < nb) + { + trace.error() << " At least " << nb + << " characters are required " + << " with option --properties."; + trace.info() << std::endl; + exit(1); + } + + typedef Z3i::Space::RealPoint RealPoint; + typedef Z3i::Space::RealPoint::Coordinate Ring; + + RealPoint border_min( vm["minAABB"].as< double >(), vm["minAABB"].as< double >(), vm["minAABB"].as< double >() ); + RealPoint border_max( vm["maxAABB"].as< double >(), vm["maxAABB"].as< double >(), vm["maxAABB"].as< double >() ); + + /// Construction of the polynomial shape + + typedef MPolynomial< 3, Ring > Polynomial3; + typedef MPolynomialReader<3, Ring> Polynomial3Reader; + typedef ImplicitPolynomial3Shape ImplicitShape; + + Polynomial3 poly; + Polynomial3Reader reader; + std::string::const_iterator iter = reader.read( poly, poly_str.begin(), poly_str.end() ); + if ( iter != poly_str.end() ) + { + std::cerr << "ERROR: I read only <" + << poly_str.substr( 0, iter - poly_str.begin() ) + << ">, and I built P=" << poly << std::endl; + return 1; + } + + ImplicitShape* shape = new ImplicitShape( poly ); + + compareShapeEstimators< Z3i::Space, ImplicitShape > ( + file_export, + shape, + border_min, border_max, + h, + radius, + alpha, + options, + properties, + lambda_optimized, + noiseLevel ); + + delete shape; } From 6e8858319a51b00adbd389eb3cb0733e874b6117 Mon Sep 17 00:00:00 2001 From: Kerautret Date: Sat, 18 Jul 2020 17:34:42 +0200 Subject: [PATCH 6/6] add doc in 3dlocalEstimators --- ChangeLog.md | 3 +- doc/estimators.dox | 2 +- estimators/3dLocalEstimators.cpp | 57 +++++++++++++++++++++++++++++++- 3 files changed, 59 insertions(+), 3 deletions(-) diff --git a/ChangeLog.md b/ChangeLog.md index e62a905d..6d434c93 100644 --- a/ChangeLog.md +++ b/ChangeLog.md @@ -38,7 +38,8 @@ - *estimators* - volSurfaceRegularization now in the "make install" command. (David Coeurjolly, [#376](https://github.com/DGtal-team/DGtalTools/pull/376)) - - 3dLocalEstimators now include in the main build and fix compilation issues + - 3dLocalEstimators now included in the main build and fix compilation issues + and documentation added. [#382](https://github.com/DGtal-team/DGtalTools/issues/382). (Bertrand Kerautret [#383](https://github.com/DGtal-team/DGtalTools/pull/383)) diff --git a/doc/estimators.dox b/doc/estimators.dox index 5ec489db..06bec656 100644 --- a/doc/estimators.dox +++ b/doc/estimators.dox @@ -18,9 +18,9 @@ @section estimators_Doc Estimators - - @ref Doc2dLocalEstimators : compares local estimators on implicit shapes using DGtal Library. - @ref Doc3dCurveTangentEstimator : estimates the tangent vector to a set of 3D integer points. + - @ref Doc3dLocalEstimators : compares local estimators on implicit shapes using DGtal library. - @ref curvatureBC : estimates curvature using a binomial convolver. - @ref curvatureMCMS : estimates curvature using length of most centered segment computers. - @ref curvatureScaleSpaceBCC : generates the Curvature Scale Space image using a binomial convolver based estimator. diff --git a/estimators/3dLocalEstimators.cpp b/estimators/3dLocalEstimators.cpp index 75b191d5..1d0af947 100644 --- a/estimators/3dLocalEstimators.cpp +++ b/estimators/3dLocalEstimators.cpp @@ -73,8 +73,60 @@ using namespace DGtal; using namespace functors; -typedef std::pair PrincipalCurvatures; + +/** + @page Doc3dLocalEstimators 3dLocalEstimators + + @brief Compares local estimators on implicit shapes using DGtal library. + + @b Usage: 3dLocalEstimators [options] --shape --h --radius --estimators --output + + Below are the different available families of estimators: + - Integral Invariant Mean + - Integral Invariant Gaussian + - Monge Jet Fitting Mean + - Monge Jet Fitting Gaussian + + @b Allowed @b options @b are : + @code + -h [ --help ] display this message + -s [ --shape ] arg Shape + -o [ --output ] arg Output file + -r [ --radius ] arg Kernel radius for IntegralInvariant + --alpha arg (=0.33333333333333331) Alpha parameter for Integral Invariant + computation + --h arg Grid step + -a [ --minAABB ] arg (=-10) Min value of the AABB bounding box + (domain) + -A [ --maxAABB ] arg (=10) Max value of the AABB bounding box + (domain) + -n [ --noise ] arg (=0) Level of noise to perturb the shape + -l [ --lambda ] arg (=0) Use the shape to get a better + approximation of the surface (optional) + --properties arg (=110) the i-th property is disabled iff there is + a 0 at position i + -e [ --estimators ] arg (=110) the i-th estimator is disabled iff there + is a 0 at position i + @endcode + + @b Example: + The following example will estimate the curvature of an implicit cone shape and will produce six resulting files (toto_II_gaussian.dat,toto_II_mean.dat, toto_MongeJetFitting_gaussian.dat, toto_MongeJetFitting_mean.dat, toto_True_gaussian.dat, toto_True_mean.dat) + @code + ./estimators/3dlocalEstimators --shape "z^2-x^2-y^2" --output result --h 0.4 --radius 1.0 + + @endcode + You can check other example of implicit shapes like the followinf: + - whitney : x^2-y*z^2 + - 4lines : x*y*(y-x)*(y-z*x) + - cone : z^2-x^2-y^2 + - simonU : x^2-z*y^2+x^4+y^4 + - cayley3 : 4*(x^2 + y^2 + z^2) + 16*x*y*z - 1 + - crixxi : -0.9*(y^2+z^2-1)^2-(x^2+y^2-1)^3 +**/ + + +typedef std::pair PrincipalCurvatures; template < typename Shape, typename KSpace, typename ConstIterator, typename OutputIterator > void estimateTrueMeanCurvatureQuantity( const ConstIterator & it_begin, @@ -1183,6 +1235,9 @@ int main( int argc, char** argv ) << "\t - Mean Curvature" << std::endl << "\t - Gaussian Curvature" << std::endl << "\t - k1/k2" << std::endl + << std::endl + << general_opt + << "Example: \n ./estimators/3dlocalEstimators --shape \"z^2-x^2-y^2\" --output result --h 0.4 --radius 1.0" << std::endl; return 0; }