DGtalTools 2.0.0
Loading...
Searching...
No Matches
2dLocalEstimators.cpp
1
39#include <iostream>
40#include <iomanip>
41#include <vector>
42#include <string>
43#include <fstream>
44
45#include "CLI11.hpp"
46
47#include "DGtal/base/Common.h"
48#include "DGtal/base/Clock.h"
49
50//shapes
51#include "DGtal/shapes/ShapeFactory.h"
52#include "DGtal/shapes/Shapes.h"
53#include "DGtal/helpers/StdDefs.h"
54#include "DGtal/topology/helpers/Surfaces.h"
55#include "DGtal/topology/DigitalSurface.h"
56
57//Digitizer
58#include "DGtal/shapes/GaussDigitizer.h"
59#include "DGtal/geometry/curves/GridCurve.h"
60#include "DGtal/topology/LightImplicitDigitalSurface.h"
61#include "DGtal/graph/DepthFirstVisitor.h"
62#include "DGtal/graph/GraphVisitorRange.h"
63#include "DGtal/geometry/volumes/KanungoNoise.h"
64
65
66//Estimators
67#include "DGtal/geometry/curves/estimation/TrueLocalEstimatorOnPoints.h"
68#include "DGtal/geometry/curves/estimation/TrueGlobalEstimatorOnPoints.h"
69#include "DGtal/geometry/curves/estimation/ParametricShapeCurvatureFunctor.h"
70#include "DGtal/geometry/curves/estimation/ParametricShapeTangentFunctor.h"
71#include "DGtal/geometry/curves/estimation/ParametricShapeArcLengthFunctor.h"
72
73#include "DGtal/geometry/curves/BinomialConvolver.h"
74#include "DGtal/geometry/curves/estimation/MostCenteredMaximalSegmentEstimator.h"
75#include "DGtal/geometry/curves/estimation/SegmentComputerEstimators.h"
76#include "DGtal/geometry/curves/ArithmeticalDSSComputer.h"
77#include "DGtal/geometry/curves/StabbingCircleComputer.h"
78
79#include "DGtal/images/ImageHelper.h"
80#include "DGtal/geometry/surfaces/estimation/IIGeometricFunctors.h"
81#include "DGtal/geometry/surfaces/estimation/IntegralInvariantVolumeEstimator.h"
82
83#include "DGtal/kernel/BasicPointFunctors.h"
84
85using namespace DGtal;
86
87
88
169std::vector<std::string> shapes2D;
170std::vector<std::string> shapesDesc;
171std::vector<std::string> shapesParam1;
172std::vector<std::string> shapesParam2;
173std::vector<std::string> shapesParam3;
174std::vector<std::string> shapesParam4;
175
176template< typename RealPoint >
177struct OptionsIntegralInvariant
178{
179 double alpha; // <! Alpha parameter for the convolution kernel. 1/3 by default
180 double radius; // <! Radius of the convolution kernel.
181 RealPoint center; // <! Center of the shape.
182 bool lambda_optimized;
183};
184
185
190void createList()
191{
192 shapes2D.push_back("ball");
193 shapesDesc.push_back("Ball for the Euclidean metric.");
194 shapesParam1.push_back("--radius [-R]");
195 shapesParam2.push_back("");
196 shapesParam3.push_back("");
197 shapesParam4.push_back("");
198
199 shapes2D.push_back("square");
200 shapesDesc.push_back("square (no signature).");
201 shapesParam1.push_back("--width [-w]");
202 shapesParam2.push_back("");
203 shapesParam3.push_back("");
204 shapesParam4.push_back("");
205
206 shapes2D.push_back("lpball");
207 shapesDesc.push_back("Ball for the l_power metric (no signature).");
208 shapesParam1.push_back("--radius [-R],");
209 shapesParam2.push_back("--power [-p]");
210 shapesParam3.push_back("");
211 shapesParam4.push_back("");
212
213 shapes2D.push_back("flower");
214 shapesDesc.push_back("Flower with k petals with radius ranging from R+/-v.");
215 shapesParam1.push_back("--radius [-R],");
216 shapesParam2.push_back("--varsmallradius [-v],");
217 shapesParam3.push_back("--k [-k],");
218 shapesParam4.push_back("--phi");
219
220 shapes2D.push_back("ngon");
221 shapesDesc.push_back("Regular k-gon.");
222 shapesParam1.push_back("--radius [-R],");
223 shapesParam2.push_back("--k [-k],");
224 shapesParam3.push_back("--phi");
225 shapesParam4.push_back("");
226
227 shapes2D.push_back("accflower");
228 shapesDesc.push_back("Accelerated Flower with k petals.");
229 shapesParam1.push_back("--radius [-R],");
230 shapesParam2.push_back("--varsmallradius [-v],");
231 shapesParam3.push_back("--k [-k],");
232 shapesParam4.push_back("--phi");
233
234 shapes2D.push_back("ellipse");
235 shapesDesc.push_back("Ellipse.");
236 shapesParam1.push_back("--axis1 [-A],");
237 shapesParam2.push_back("--axis2 [-a],");
238 shapesParam3.push_back("--phi");
239 shapesParam4.push_back("");
240
241
242}
243
248void displayList()
249{
250 trace.emphase()<<"2D Shapes:"<<std::endl;
251 for(unsigned int i=0; i<shapes2D.size(); ++i)
252 trace.info()<<"\t"<<shapes2D[i]<<"\t"
253 <<shapesDesc[i]<<std::endl
254 <<"\t\tRequired parameter(s): "
255 << shapesParam1[i]<<" "
256 << shapesParam2[i]<<" "
257 << shapesParam3[i]<<" "
258 << shapesParam4[i]<<std::endl;
259
260}
261
262
271unsigned int checkAndReturnIndex(const std::string &shapeName)
272{
273 unsigned int pos=0;
274
275 while ((pos < shapes2D.size()) && (shapes2D[pos] != shapeName))
276 pos++;
277
278 if (pos == shapes2D.size())
279 {
280 trace.error() << "The specified shape has not found.";
281 trace.info() << std::endl;
282 exit(1);
283 }
284
285 return pos;
286}
287
293void missingParam(std::string param)
294{
295 trace.error() <<" Parameter: "<<param<<" is required.";
296 trace.info()<<std::endl;
297 exit(1);
298}
299
306void estimationError(int currentSize, int expectedSize)
307{
308 if (currentSize != expectedSize)
309 {
310 trace.error() << " error in the estimation"
311 << " (got " << currentSize << " values"
312 << " instead of " << expectedSize << ")";
313 trace.info() << std::endl;
314 exit(1);
315 }
316
317}
318
329template <typename Estimator, typename ConstIterator, typename OutputIterator>
330void
331estimation( Estimator & estimator, double h,
332 const ConstIterator& itb, const ConstIterator& ite, const OutputIterator& ito )
333{
334 Clock c;
335 c.startClock();
336 estimator.eval( itb, ite, ito, h );
337 double time = c.stopClock();
338 std::cout << "# Time: " << time << std::endl;
339}
340
341
346template< typename ConstIteratorOnPoints, typename RPoint >
347unsigned int suggestedSizeIntegralInvariant( const double h,
348 const RPoint& center,
349 const ConstIteratorOnPoints& itb,
350 const ConstIteratorOnPoints& ite )
351{
352 ConstIteratorOnPoints it = itb;
353 RPoint p( *it );
354 RPoint distance = p - center;
355 auto minRadius = distance.norm();
356 ++it;
357
358 for ( ; it != ite; ++it )
359 {
360 p = *it;
361 distance = p - center;
362 if ( distance.norm() < minRadius )
363 {
364 minRadius = distance.norm();
365 }
366 }
367
368 return minRadius * h;
369}
370
371
384template <typename Space, typename Shape>
385bool
386computeLocalEstimations( const std::string & filename,
387 const Shape& aShape,
388 const double & h,
389 struct OptionsIntegralInvariant< Z2i::RealPoint > optionsII,
390 const std::string & options,
391 const std::string & properties,
392 const std::string & outShape,
393 double noiseLevel = 0.0 )
394{
395 // Types
396 typedef typename Space::Vector Vector;
397 typedef typename Space::RealPoint RealPoint;
398 typedef typename Space::Integer Integer;
399 typedef HyperRectDomain<Space> Domain;
400 typedef KhalimskySpaceND<Space::dimension,Integer> KSpace;
401 typedef typename KSpace::SCell SCell;
402 typedef GaussDigitizer<Space,Shape> Digitizer;
403 typedef KanungoNoise< Digitizer, Z2i::Domain > KanungoPredicate;
404
405 bool withNoise = ( noiseLevel <= 0.0 ) ? false : true;
406 /*if( withNoise )
407 noiseLevel = std::pow(noiseLevel, h);*/
408
409 ASSERT (( noiseLevel < 1.0 ));
410
411 bool tangent = ( properties.at( 0 ) != '0' ) ? true : false;
412 bool curvature = ( properties.at( 1 ) != '0' ) ? true : false;
413
414 // Digitizer
415 Digitizer* dig = new Digitizer();
416 dig->attach( aShape ); // attaches the shape.
417 Vector vlow(-1,-1); Vector vup(1,1);
418 dig->init( aShape.getLowerBound()+vlow, aShape.getUpperBound()+vup, h );
419 Domain domain = dig->getDomain();
420
421 //Noise
422
423 Clock c;
424
425 // Create cellular space
426 KSpace K;
427 bool ok = K.init( dig->getLowerBound(), dig->getUpperBound(), true );
428 if ( ! ok )
429 {
430 std::cerr << "[2dLocalEstimators]"
431 << " error in creating KSpace." << std::endl;
432 return false;
433 }
434 try {
435
436 // Extracts shape boundary
437 SurfelAdjacency< KSpace::dimension > SAdj( true );
438 SCell bel;
439 std::vector< SCell > points;
440
441 KanungoPredicate *noisifiedObject;
442 if ( withNoise )
443 {
444 noisifiedObject = new KanungoPredicate( *dig, domain, noiseLevel );
445 bel = Surfaces< KSpace >::findABel( K, *noisifiedObject, 10000 );
446 Surfaces< KSpace >::track2DBoundary( points, K, SAdj, *noisifiedObject, bel );
447
448 double minsize = dig->getUpperBound()[0] - dig->getLowerBound()[0];
449 while( points.size() < 2 * minsize )
450 {
451 points.clear();
452 bel = Surfaces< KSpace >::findABel( K, *noisifiedObject, 10000 );
453 Surfaces< KSpace >::track2DBoundary( points, K, SAdj, *noisifiedObject, bel );
454 }
455 }
456 else
457 {
458 bel = Surfaces< KSpace >::findABel( K, *dig, 10000 );
459 Surfaces< KSpace >::track2DBoundary( points, K, SAdj, *dig, bel );
460 }
461
462 // Create GridCurve
463 GridCurve< KSpace > gridcurve;
464 gridcurve.initFromSCellsVector( points );
465 if(outShape != "")
466 {
467 std::ofstream outS;
468 outS.open(outShape.c_str());
469 for(const auto &p : points)
470 {
471
472 Dimension track = *( K.sDirs( p ) );
473 SCell pointel = K.sIndirectIncident( p, track );
474 outS << K.sCoords( pointel )[0] << " " << K.sCoords( pointel )[1] << std::endl;
475 }
476 outS.close();
477 }
478 // Ranges
479 typedef typename GridCurve< KSpace >::MidPointsRange PointsRange;
480 PointsRange pointsRange = gridcurve.getMidPointsRange();
481
482 // Estimations
483 if (gridcurve.isClosed())
484 {
485 if (options.at(0) != '0')
486 {
487 if( tangent )
488 {
489 char full_filename[360];
490 sprintf( full_filename, "%s%s", filename.c_str(), "_True_tangeant.dat" );
491 std::ofstream file( full_filename );
492
493 file << "# h = " << h << std::endl;
494 file << "# True tangents computation" << std::endl;
495 file << "# range size = " << pointsRange.size() << std::endl;
496 if ( withNoise )
497 {
498 file << "# noise level (init) = " << noiseLevel/h << std::endl;
499 file << "# noise level (current) = " << noiseLevel << std::endl;
500 }
501
502 std::ostream_iterator< RealPoint > out_it( file, "\n" );
503
504 typedef ParametricShapeTangentFunctor< Shape > TangentFunctor;
505 typedef typename PointsRange::ConstCirculator C;
506 TrueLocalEstimatorOnPoints< C, Shape, TangentFunctor >
507 trueTangentEstimator;
508 trueTangentEstimator.attach( aShape );
509 estimation( trueTangentEstimator, h,
510 pointsRange.c(), pointsRange.c(),
511 out_it );
512
513 file.close();
514
515 }
516
517 if( curvature )
518 {
519 char full_filename[360];
520 sprintf( full_filename, "%s%s", filename.c_str(), "_True_curvature.dat" );
521 std::ofstream file( full_filename );
522
523 file << "# h = " << h << std::endl;
524 file << "# True curvature computation" << std::endl;
525 file << "# range size = " << pointsRange.size() << std::endl;
526 if ( withNoise )
527 {
528 file << "# noise level (init) = " << noiseLevel/h << std::endl;
529 file << "# noise level (current) = " << noiseLevel << std::endl;
530 }
531
532 std::ostream_iterator< double > out_it( file, "\n" );
533
534 typedef ParametricShapeCurvatureFunctor< Shape > CurvatureFunctor;
535 typedef typename PointsRange::ConstCirculator C;
536 TrueLocalEstimatorOnPoints< C, Shape, CurvatureFunctor >
537 trueCurvatureEstimator;
538 trueCurvatureEstimator.attach( aShape );
539 estimation( trueCurvatureEstimator, h,
540 pointsRange.c(), pointsRange.c(),
541 out_it );
542
543 file.close();
544 }
545 }
546
547 // Maximal Segments
548 if (options.at(1) != '0')
549 {
550 if( tangent )
551 {
552 char full_filename[360];
553 sprintf( full_filename, "%s%s", filename.c_str(), "_MDSS_tangeant.dat" );
554 std::ofstream file( full_filename );
555
556 file << "# h = " << h << std::endl;
557 file << "# Most centered maximal DSS tangent estimation" << std::endl;
558 file << "# range size = " << pointsRange.size() << std::endl;
559 if ( withNoise )
560 {
561 file << "# noise level (init) = " << noiseLevel/h << std::endl;
562 file << "# noise level (current) = " << noiseLevel << std::endl;
563 }
564
565 std::ostream_iterator< RealPoint > out_it( file, "\n" );
566
567 typedef typename GridCurve< KSpace >::PointsRange PointsRange2;
568 PointsRange2 pointsRange2 = gridcurve.getPointsRange();
569
570 typedef typename PointsRange2::ConstCirculator C;
571 typedef ArithmeticalDSSComputer< C, int, 4 > SegmentComputer;
572 typedef TangentFromDSSEstimator<SegmentComputer> SCFunctor;
573 SegmentComputer sc;
574 SCFunctor f;
575
576
577 MostCenteredMaximalSegmentEstimator<SegmentComputer,SCFunctor> MDSSTangentEstimator(sc, f);
578 estimation( MDSSTangentEstimator, h,
579 pointsRange2.c(), pointsRange2.c(),
580 out_it );
581
582 file.close();
583 }
584 if( curvature )
585 {
586 c.startClock();
587
588 char full_filename[360];
589 sprintf( full_filename, "%s%s", filename.c_str(), "_MDSSl_curvature.dat" );
590 std::ofstream file( full_filename );
591
592 file << "# h = " << h << std::endl;
593 file << "# Most centered maximal DSS (length) curvature estimation" << std::endl;
594 file << "# range size = " << pointsRange.size() << std::endl;
595 if ( withNoise )
596 {
597 file << "# noise level (init) = " << noiseLevel/h << std::endl;
598 file << "# noise level (current) = " << noiseLevel << std::endl;
599 }
600
601 std::ostream_iterator< double > out_it( file, "\n" );
602
603 typedef typename GridCurve< KSpace >::PointsRange PointsRange2;
604 PointsRange2 pointsRange2 = gridcurve.getPointsRange();
605
606 typedef typename PointsRange2::ConstCirculator C;
607 typedef ArithmeticalDSSComputer< C, int, 4 > SegmentComputer;
608 typedef CurvatureFromDSSLengthEstimator<SegmentComputer> SCFunctor;
609 SegmentComputer sc;
610 SCFunctor f;
611 MostCenteredMaximalSegmentEstimator<SegmentComputer,SCFunctor> MDSSCurvatureEstimator(sc, f);
612
613 estimation( MDSSCurvatureEstimator, h,
614 pointsRange2.c(), pointsRange2.c(),
615 out_it );
616
617 file.close();
618
619
620 memset(&full_filename[0], 0, sizeof(full_filename));
621 sprintf( full_filename, "%s%s", filename.c_str(), "_MDSSlw_curvature.dat" );
622 file.open( full_filename );
623
624 file << "# h = " << h << std::endl;
625 file << "# Most centered maximal DSS (length & width) curvature estimation" << std::endl;
626 file << "# range size = " << pointsRange.size() << std::endl;
627 if ( withNoise )
628 {
629 file << "# noise level (init) = " << noiseLevel/h << std::endl;
630 file << "# noise level (current) = " << noiseLevel << std::endl;
631 }
632
633 std::ostream_iterator< double > out_it2( file, "\n" );
634
635 typedef CurvatureFromDSSEstimator<SegmentComputer> SCFunctor2;
636 SegmentComputer sc2;
637 SCFunctor2 f2;
638 MostCenteredMaximalSegmentEstimator<SegmentComputer,SCFunctor2> MDSSCurvatureEstimator2(sc2, f2);
639 estimation( MDSSCurvatureEstimator2, h,
640 pointsRange2.c(), pointsRange2.c(),
641 out_it2 );
642
643 double time = c.stopClock();
644 file << "# Time: " << time << std::endl;
645
646 file.close();
647
648 }
649 }
650
651 //Maximal circular arcs
652 if (options.at(2) != '0')
653 {
654 if( tangent )
655 {
656 char full_filename[360];
657 sprintf( full_filename, "%s%s", filename.c_str(), "_MDCA_tangent.dat" );
658 std::ofstream file( full_filename );
659
660 file << "# h = " << h << std::endl;
661 file << "# Most centered maximal DCA tangents estimation" << std::endl;
662 file << "# range size = " << pointsRange.size() << std::endl;
663 if ( withNoise )
664 {
665 file << "# noise level (init) = " << noiseLevel/h << std::endl;
666 file << "# noise level (current) = " << noiseLevel << std::endl;
667 }
668
669 std::ostream_iterator< RealPoint > out_it( file, "\n" );
670
671 typedef typename GridCurve<KSpace>::IncidentPointsRange Range;
672 typedef typename Range::ConstCirculator C;
673 Range r = gridcurve.getIncidentPointsRange();
674 typedef StabbingCircleComputer<C> SegmentComputer;
675 typedef TangentFromDCAEstimator<SegmentComputer> SCFunctor;
676 SegmentComputer sc;
677 SCFunctor f;
678 MostCenteredMaximalSegmentEstimator<SegmentComputer,SCFunctor> MDCATangentEstimator(sc, f);
679 estimation( MDCATangentEstimator, h,
680 r.c(), r.c(),
681 out_it );
682
683 file.close();
684 }
685
686 if( curvature )
687 {
688 c.startClock();
689
690 char full_filename[360];
691 sprintf( full_filename, "%s%s", filename.c_str(), "_MDCA_curvature.dat" );
692 std::ofstream file( full_filename );
693
694 file << "# h = " << h << std::endl;
695 file << "# Most centered maximal DCA curvature estimation" << std::endl;
696 file << "# range size = " << pointsRange.size() << std::endl;
697 if ( withNoise )
698 {
699 file << "# noise level (init) = " << noiseLevel/h << std::endl;
700 file << "# noise level (current) = " << noiseLevel << std::endl;
701 }
702
703 std::ostream_iterator< double > out_it( file, "\n" );
704
705 typedef typename GridCurve<KSpace>::IncidentPointsRange Range;
706 typedef typename Range::ConstCirculator C;
707 Range r = gridcurve.getIncidentPointsRange();
708 typedef StabbingCircleComputer<C> SegmentComputer;
709 typedef CurvatureFromDCAEstimator<SegmentComputer, false> SCFunctor;
710 SegmentComputer sc;
711 SCFunctor f;
712 MostCenteredMaximalSegmentEstimator<SegmentComputer,SCFunctor> MDCACurvatureEstimator(sc, f);
713 estimation( MDCACurvatureEstimator, h,
714 r.c(), r.c(),
715 out_it );
716
717 double time = c.stopClock();
718 file << "# Time: " << time << std::endl;
719
720 file.close();
721 }
722 }
723
724 //Binomial convolver
725 if (options.at(3) != '0')
726 {
727 if( tangent )
728 {
729 c.startClock();
730
731 char full_filename[360];
732 sprintf( full_filename, "%s%s", filename.c_str(), "_BC_tangeant.dat" );
733 std::ofstream file( full_filename );
734
735 file << "# h = " << h << std::endl;
736 file << "# Tangents estimation from binomial convolution" << std::endl;
737 file << "# range size = " << pointsRange.size() << std::endl;
738 if ( withNoise )
739 {
740 file << "# noise level (init) = " << noiseLevel/h << std::endl;
741 file << "# noise level (current) = " << noiseLevel << std::endl;
742 }
743
744 typedef typename PointsRange::ConstIterator I;
745 typedef BinomialConvolver<I, double> MyBinomialConvolver;
746 file << "# mask size = " <<
747 MyBinomialConvolver::suggestedSize( h, pointsRange.begin(), pointsRange.end() ) << std::endl;
748
749 typedef TangentFromBinomialConvolverFunctor< MyBinomialConvolver, RealPoint >
750 TangentBCFct;
751 BinomialConvolverEstimator< MyBinomialConvolver, TangentBCFct> BCTangentEstimator;
752
753 std::ostream_iterator< RealPoint > out_it( file, "\n" );
754
755 BCTangentEstimator.init( h, pointsRange.begin(), pointsRange.end(), true );
756 BCTangentEstimator.eval( pointsRange.begin(), pointsRange.end(), out_it );
757
758 double time = c.stopClock();
759 file << "# Time: " << time << std::endl;
760
761 file.close();
762 }
763
764 if( curvature )
765 {
766 c.startClock();
767
768 char full_filename[360];
769 sprintf( full_filename, "%s%s", filename.c_str(), "_BC_curvature.dat" );
770 std::ofstream file( full_filename );
771
772 file << "# h = " << h << std::endl;
773 file << "# Curvature estimation from binomial convolution" << std::endl;
774 file << "# range size = " << pointsRange.size() << std::endl;
775 if ( withNoise )
776 {
777 file << "# noise level (init) = " << noiseLevel/h << std::endl;
778 file << "# noise level (current) = " << noiseLevel << std::endl;
779 }
780
781 typedef typename PointsRange::ConstIterator I;
782 typedef BinomialConvolver<I, double> MyBinomialConvolver;
783 file << "# mask size = " <<
784 MyBinomialConvolver::suggestedSize( h, pointsRange.begin(), pointsRange.end() ) << std::endl;
785
786 std::ostream_iterator< double > out_it( file, "\n" );
787
788 typedef CurvatureFromBinomialConvolverFunctor< MyBinomialConvolver, double >
789 CurvatureBCFct;
790 BinomialConvolverEstimator< MyBinomialConvolver, CurvatureBCFct> BCCurvatureEstimator;
791
792 BCCurvatureEstimator.init( h, pointsRange.begin(), pointsRange.end(), true );
793 BCCurvatureEstimator.eval( pointsRange.begin(), pointsRange.end(), out_it );
794
795 double time = c.stopClock();
796 file << "# Time: " << time << std::endl;
797
798 file.close();
799 }
800 }
801
803 //Integral Invariants
804 if (options.at(4) != '0')
805 {
806 if( curvature )
807 {
808 c.startClock();
809
810 char full_filename[360];
811 sprintf( full_filename, "%s%s", filename.c_str(), "_II_curvature.dat" );
812 std::ofstream file( full_filename );
813
814 file << "# h = " << h << std::endl;
815 file << "# Integral Invariant curvature estimation" << std::endl;
816 file << "# range size = " << pointsRange.size() << std::endl;
817 if ( withNoise )
818 {
819 file << "# noise level (init) = " << noiseLevel/h << std::endl;
820 file << "# noise level (current) = " << noiseLevel << std::endl;
821 }
822
823 if( optionsII.radius <= 0.0 )
824 {
825 optionsII.radius = suggestedSizeIntegralInvariant( h, optionsII.center, pointsRange.begin(), pointsRange.end() );
826 file << "# Estimated radius: " << optionsII.radius << std::endl;
827 }
828 double re = optionsII.radius * std::pow( h, optionsII.alpha );
829 file << "# full kernel (digital) size (with alpha = " << optionsII.alpha << ") = " <<
830 re / h << std::endl;
831
832 std::ostream_iterator< double > out_it( file, "\n" );
833
834 if ( withNoise )
835 {
836 typedef LightImplicitDigitalSurface< KSpace, KanungoPredicate > LightImplicitDigSurface;
837 typedef DigitalSurface< LightImplicitDigSurface > DigSurface;
838
839 LightImplicitDigSurface LightImplDigSurf( K, *noisifiedObject, SAdj, bel );
840 DigSurface surf( LightImplDigSurf );
841
842 typedef DepthFirstVisitor< DigSurface > Visitor;
843 typedef GraphVisitorRange< Visitor > VisitorRange;
844 typedef typename VisitorRange::ConstIterator VisitorConstIterator;
845
846 VisitorRange range( new Visitor( surf, *surf.begin() ));
847 VisitorConstIterator ibegin = range.begin();
848 VisitorConstIterator iend = range.end();
849
850 typedef functors::IICurvatureFunctor<Z2i::Space> MyIICurvatureFunctor;
851 typedef IntegralInvariantVolumeEstimator< KSpace, KanungoPredicate, MyIICurvatureFunctor > MyIICurvatureEstimator;
852
853 MyIICurvatureFunctor curvatureFunctor;
854 curvatureFunctor.init( h, re );
855
856 MyIICurvatureEstimator curvatureEstimator( curvatureFunctor );
857 curvatureEstimator.attach( K, *noisifiedObject );
858 curvatureEstimator.setParams( re/h );
859 curvatureEstimator.init( h, ibegin, iend );
860
861 curvatureEstimator.eval( ibegin, iend, out_it );
862 }
863 else
864 {
865 typedef LightImplicitDigitalSurface< KSpace, Digitizer > LightImplicitDigSurface;
866 typedef DigitalSurface< LightImplicitDigSurface > DigSurface;
867
868 LightImplicitDigSurface LightImplDigSurf( K, *dig, SAdj, bel );
869 DigSurface surf( LightImplDigSurf );
870
871 typedef DepthFirstVisitor< DigSurface > Visitor;
872 typedef GraphVisitorRange< Visitor > VisitorRange;
873 typedef typename VisitorRange::ConstIterator VisitorConstIterator;
874
875 VisitorRange range( new Visitor( surf, *surf.begin() ));
876 VisitorConstIterator ibegin = range.begin();
877 VisitorConstIterator iend = range.end();
878
879 typedef functors::IICurvatureFunctor<Z2i::Space> MyIICurvatureFunctor;
880 typedef IntegralInvariantVolumeEstimator< KSpace, Digitizer, MyIICurvatureFunctor > MyIICurvatureEstimator;
881
882 MyIICurvatureFunctor curvatureFunctor;
883 curvatureFunctor.init( h, re );
884
885 MyIICurvatureEstimator curvatureEstimator( curvatureFunctor );
886 curvatureEstimator.attach( K, *dig );
887 curvatureEstimator.setParams( re/h );
888 curvatureEstimator.init( h, ibegin, iend );
889
890 curvatureEstimator.eval( ibegin, iend, out_it );
891 }
892
893 double time = c.stopClock();
894 file << "# Time: " << time << std::endl;
895
896 file.close();
897 }
898 }
899
900 //delete noisifiedObject;
901 delete dig;
902 }
903 else
904 {
905 //delete noisifiedObject;
906 delete dig;
907 std::cerr << "[computeLocalEstimations]"
908 << " error: open digital curve found." << std::endl;
909 return false;
910 }
911 }
912 catch ( InputException e )
913 {
914 std::cerr << "[computeLocalEstimations]"
915 << " error in finding a bel." << std::endl;
916 return false;
917 }
918 return true;
919}
920
921
923
924
925int main( int argc, char** argv )
926{
927 // parse command line CLI ----------------------------------------------
928 CLI::App app;
929 std::string shapeName;
930 std::string filename;
931 double radius, kernelradius;
932 double power {2.0};
933 double smallradius {5};
934 double varsmallradius {5};
935 double cx {0.0}, cy {0.0};
936 double h {1.0};
937 unsigned int k {3};
938 double phi {0.0};
939 double width {10.0};
940 double axis1, axis2;
941 double alpha {1.0/3.0};
942 double noiseLevel {0.0};
943 std::string properties {"11"};
944 std::string outShape {""};
945 bool lambda {false};
946 std::string options {"1000"};
947
948 app.description("Compares local estimators on implicit shapes using DGtal library.\n Typical use example:\n \t 2dlocalEstimators --output <output> --shape <shapeName> [required parameters] --estimators <binaryWord> --properties <binaryWord>\n");
949 auto listOpt = app.add_flag("--list,-l","List all available shapes");
950 auto outputOpt = app.add_option("--output,-o", filename, "Output")->required();
951 auto shapeNameOpt = app.add_option("--shape,-s", shapeName, "Shape name")->required();
952 auto radiusOpt = app.add_option("--radius,-R", radius, "Radius of the shape" );
953 auto kernelradiusOpt = app.add_option("--kernelradius,-K", radius, "Radius of the convolution kernel (Integral invariants estimators)");
954 auto alphaOpt = app.add_option("--alpha", alpha, "Alpha parameter for Integral Invariant computation");
955 auto axis1Opt = app.add_option("--axis1,-A", axis1, "Half big axis of the shape (ellipse)" );
956 auto axis2Opt = app.add_option("--axis2,-a", axis2, "Half small axis of the shape (ellipse)" );
957 auto smallradiusOpt = app.add_option("--smallradius,-r", smallradius, "Small radius of the shape");
958 auto varsmallradiusOpt = app.add_option("--varsmallradius,-v", varsmallradius, "Variable small radius of the shape" );
959 auto kOpt = app.add_option("-k", k, "Number of branches or corners the shape (default 3)" );
960 auto phiOpt = app.add_option("--phi", phi, "Phase of the shape (in radian)" );
961 auto widthOpt = app.add_option("--width,-w", width, "Width of the shape" );
962 auto powerOpt = app.add_option("--power,-p", power, "Power of the metric" );
963 app.add_option("--center_x,-x", cx, "x-coordinate of the shape center" );
964 app.add_option("--center_y,-y", cy, "y-coordinate of the shape center" );
965 app.add_option("--gridstep,-g", h, "Gridstep for the digitization" );
966 app.add_option("--noise,-n", noiseLevel, "Level of noise to perturb the shape");
967 app.add_option("--properties", properties, "the i-th property is disabled iff there is a 0 at position i");
968 app.add_option("--estimators,-e", options, "the i-th estimator is disabled iff there is a 0 at position i");
969 app.add_option("--exportShape,-E", outShape, "Exports the contour of the source shape as a sequence of discrete points (.sdp)");
970 app.add_option("--lambda", lambda, "Use the shape to get a better approximation of the surface (optional)");
971
972 app.get_formatter()->column_width(40);
973 CLI11_PARSE(app, argc, argv);
974 // END parse command line using CLI
975
976 //List creation
977 createList();
978
979 if ( listOpt->count() > 0 )
980 {
981 displayList();
982 return 0;
983 }
984
985 unsigned int nb = 4; //number of available methods
986 if (options.size() < nb)
987 {
988 trace.error() << " At least " << nb
989 << " characters are required "
990 << " with option --estimators.";
991 trace.info() << std::endl;
992 exit(1);
993 }
994
995 nb = 2; //number of available properties
996 if (properties.size() < nb)
997 {
998 trace.error() << " At least " << nb
999 << " characters are required "
1000 << " with option --properties.";
1001 trace.info() << std::endl;
1002 exit(1);
1003 }
1004
1005 //We check that the shape is known
1006 unsigned int id = checkAndReturnIndex(shapeName);
1007
1008 // standard types
1009 typedef Z2i::Space Space;
1010 typedef Space::RealPoint RealPoint;
1011
1012 RealPoint center( cx, cy );
1013
1014 struct OptionsIntegralInvariant< RealPoint > optII;
1015 optII.radius = kernelradius;
1016 optII.alpha = alpha;
1017 optII.lambda_optimized = lambda;
1018 optII.center = center;
1019
1020 if (id ==0)
1021 {
1022 if (radiusOpt->count()==0) missingParam("--radius");
1023 //if (!(vm.count("kernelradius"))) missingParam("--kernelradius");
1024
1025 Ball2D<Space> ball( center, radius );
1026 computeLocalEstimations<Space>( filename, ball, h, optII, options, properties, outShape, noiseLevel );
1027 }
1028 else if (id ==1)
1029 {
1030 //if (widthOpt->count()==0) missingParam("--width");
1031
1032 ImplicitHyperCube<Space> object(Z2i::Point(0,0), width/2);
1033 trace.error()<< "Not available.";
1034 trace.info()<<std::endl;
1035 }
1036 else if (id ==2)
1037 {
1038 //if (powerOpt->count()==0) missingParam("--power");
1039 if (radiusOpt->count()==0) missingParam("--radius");
1040
1041 ImplicitRoundedHyperCube<Space> ball( Z2i::Point(0,0), radius, power );
1042 trace.error()<< "Not available.";
1043 trace.info()<<std::endl;
1044 }
1045 else if (id ==3)
1046 {
1047 //if (varsmallradiusOpt->count()==0) missingParam("--varsmallradius");
1048 if (radiusOpt->count()==0) missingParam("--radius");
1049 //if (kOpt->count()==0) missingParam("--k");
1050 //if (phiOpt->count()==0) missingParam("--phi");
1051 //if (!(vm.count("kernelradius"))) missingParam("--kernelradius");
1052
1053 Flower2D<Space> flower( center, radius, varsmallradius, k, phi );
1054 computeLocalEstimations<Space>( filename, flower, h, optII, options, properties, outShape, noiseLevel );
1055 }
1056 else if (id ==4)
1057 {
1058 if (radiusOpt->count()==0) missingParam("--radius");
1059 //if (kOpt->count()==0) missingParam("--k");
1060 //if (phiOpt->count()==0) missingParam("--phi");
1061 //if (!(vm.count("kernelradius"))) missingParam("--kernelradius");
1062
1063 NGon2D<Space> object( center, radius, k, phi );
1064 computeLocalEstimations<Space>( filename, object, h, optII, options, properties, outShape, noiseLevel );
1065 }
1066 else if (id ==5)
1067 {
1068 //if (varsmallradiusOpt->count()==0) missingParam("--varsmallradius");
1069 if (radiusOpt->count()==0) missingParam("--radius");
1070 //if (kOpt->count()==0) missingParam("--k");
1071 //if (phiOpt->count()==0) missingParam("--phi");
1072 //if (!(vm.count("kernelradius"))) missingParam("--kernelradius");
1073
1074 AccFlower2D<Space> accflower( center, radius, varsmallradius, k, phi );
1075 computeLocalEstimations<Space>( filename, accflower, h, optII, options, properties, outShape, noiseLevel );
1076 }
1077 else if (id ==6)
1078 {
1079 if (axis1Opt->count()==0) missingParam("--axis1");
1080 if (axis2Opt->count()==0) missingParam("--axis2");
1081 //if (phiOpt->count()==0) missingParam("--phi");
1082 //if (!(vm.count("kernelradius"))) missingParam("--kernelradius");
1083
1084
1085 Ellipse2D<Space> ellipse( center, axis1, axis2, phi );
1086 computeLocalEstimations<Space>( filename, ellipse, h, optII, options, properties, outShape, noiseLevel );
1087 }
1088}
DGtal
Definition ATu0v1.h:57