DGtalTools  1.2.0
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 
85 using namespace DGtal;
86 
87 
88 
168 std::vector<std::string> shapes2D;
169 std::vector<std::string> shapesDesc;
170 std::vector<std::string> shapesParam1;
171 std::vector<std::string> shapesParam2;
172 std::vector<std::string> shapesParam3;
173 std::vector<std::string> shapesParam4;
174 
175 template< typename RealPoint >
176 struct 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 
189 void 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 
247 void 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 
270 unsigned 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 
292 void missingParam(std::string param)
293 {
294  trace.error() <<" Parameter: "<<param<<" is required.";
295  trace.info()<<std::endl;
296  exit(1);
297 }
298 
305 void 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 
328 template <typename Estimator, typename ConstIterator, typename OutputIterator>
329 void
330 estimation( 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 
345 template< typename ConstIteratorOnPoints, typename RPoint >
346 unsigned 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 
383 template <typename Space, typename Shape>
384 bool
385 computeLocalEstimations( 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 
924 int 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)", true);
953  auto alphaOpt = app.add_option("--alpha", alpha, "Alpha parameter for Integral Invariant computation", true);
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", true);
957  auto varsmallradiusOpt = app.add_option("--varsmallradius,-v", varsmallradius, "Variable small radius of the shape", true );
958  auto kOpt = app.add_option("-k", k, "Number of branches or corners the shape (default 3)", true );
959  auto phiOpt = app.add_option("--phi", phi, "Phase of the shape (in radian)", true );
960  auto widthOpt = app.add_option("--width,-w", width, "Width of the shape", true );
961  auto powerOpt = app.add_option("--power,-p", power, "Power of the metric", true );
962  app.add_option("--center_x,-x", cx, "x-coordinate of the shape center", true );
963  app.add_option("--center_y,-y", cy, "y-coordinate of the shape center", true );
964  app.add_option("--gridstep,-g", h, "Gridstep for the digitization", true );
965  app.add_option("--noise,-n", noiseLevel, "Level of noise to perturb the shape", true);
966  app.add_option("--properties", properties, "the i-th property is disabled iff there is a 0 at position i", true);
967  app.add_option("--estimators,-e", options, "the i-th estimator is disabled iff there is a 0 at position i", true);
968  app.add_option("--exportShape,-E", outShape, "Exports the contour of the source shape as a sequence of discrete points (.sdp)", true);
969  app.add_option("--lambda", lambda, "Use the shape to get a better approximation of the surface (optional)", true);
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
Companion project associated to DGtal ; Generated on Wed Jun 2 2021 15:34:51 for DGtalTools by   doxygen 1.9.1