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