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