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