DGtalTools  0.9.2
lengthEstimators.cpp
1 
32 #include <cmath>
34 #include <iostream>
35 #include <iomanip>
36 #include <vector>
37 #include <string>
38 
39 #include <boost/program_options/options_description.hpp>
40 #include <boost/program_options/parsers.hpp>
41 #include <boost/program_options/variables_map.hpp>
42 
43 #include "DGtal/base/Common.h"
44 #include "DGtal/base/Clock.h"
45 
46 //space / domain
47 #include "DGtal/kernel/SpaceND.h"
48 #include "DGtal/kernel/domains/HyperRectDomain.h"
49 #include "DGtal/topology/KhalimskySpaceND.h"
50 
51 //shape and digitizer
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 
57 #include "DGtal/shapes/GaussDigitizer.h"
58 #include "DGtal/geometry/curves/GridCurve.h"
59 
60 //estimators
61 #include "DGtal/geometry/curves/estimation/TrueLocalEstimatorOnPoints.h"
62 #include "DGtal/geometry/curves/estimation/TrueGlobalEstimatorOnPoints.h"
63 #include "DGtal/geometry/curves/estimation/ParametricShapeArcLengthFunctor.h"
64 
65 #include "DGtal/geometry/curves/estimation/L1LengthEstimator.h"
66 #include "DGtal/geometry/curves/estimation/BLUELocalLengthEstimator.h"
67 #include "DGtal/geometry/curves/estimation/RosenProffittLocalLengthEstimator.h"
68 #include "DGtal/geometry/curves/estimation/MLPLengthEstimator.h"
69 #include "DGtal/geometry/curves/estimation/FPLengthEstimator.h"
70 #include "DGtal/geometry/curves/estimation/DSSLengthEstimator.h"
71 
72 
73 using namespace DGtal;
74 
75 
76 
153 std::vector<std::string> shapes2D;
154 std::vector<std::string> shapesDesc;
155 std::vector<std::string> shapesParam1;
156 std::vector<std::string> shapesParam2;
157 std::vector<std::string> shapesParam3;
158 std::vector<std::string> shapesParam4;
159 
160 
165 void createList()
166 {
167  shapes2D.push_back("ball");
168  shapesDesc.push_back("Ball for the Euclidean metric.");
169  shapesParam1.push_back("--radius [-R]");
170  shapesParam2.push_back("");
171  shapesParam3.push_back("");
172  shapesParam4.push_back("");
173 
174  shapes2D.push_back("square");
175  shapesDesc.push_back("square (no signature).");
176  shapesParam1.push_back("--width [-w]");
177  shapesParam2.push_back("");
178  shapesParam3.push_back("");
179  shapesParam4.push_back("");
180 
181  shapes2D.push_back("lpball");
182  shapesDesc.push_back("Ball for the l_power metric (no signature).");
183  shapesParam1.push_back("--radius [-R],");
184  shapesParam2.push_back("--power [-p]");
185  shapesParam3.push_back("");
186  shapesParam4.push_back("");
187 
188  shapes2D.push_back("flower");
189  shapesDesc.push_back("Flower with k petals.");
190  shapesParam1.push_back("--radius [-R],");
191  shapesParam2.push_back("--varsmallradius [-v],");
192  shapesParam3.push_back("--k [-k],");
193  shapesParam4.push_back("--phi");
194 
195  shapes2D.push_back("ngon");
196  shapesDesc.push_back("Regular k-gon.");
197  shapesParam1.push_back("--radius [-R],");
198  shapesParam2.push_back("--k [-k],");
199  shapesParam3.push_back("--phi");
200  shapesParam4.push_back("");
201 
202  shapes2D.push_back("accflower");
203  shapesDesc.push_back("Accelerated Flower with k petals.");
204  shapesParam1.push_back("--radius [-R],");
205  shapesParam2.push_back("--varsmallradius [-v],");
206  shapesParam3.push_back("--k [-k],");
207  shapesParam4.push_back("--phi");
208 
209  shapes2D.push_back("ellipse");
210  shapesDesc.push_back("Ellipse.");
211  shapesParam1.push_back("--axis1 [-A],");
212  shapesParam2.push_back("--axis2 [-a],");
213  shapesParam3.push_back("--phi");
214  shapesParam4.push_back("");
215 
216 
217 }
218 
223 void displayList()
224 {
225  trace.emphase()<<"2D Shapes:"<<std::endl;
226  for(unsigned int i=0; i<shapes2D.size(); ++i)
227  trace.info()<<"\t"<<shapes2D[i]<<"\t"
228  <<shapesDesc[i]<<std::endl
229  <<"\t\tRequired parameter(s): "
230  << shapesParam1[i]<<" "
231  << shapesParam2[i]<<" "
232  << shapesParam3[i]<<" "
233  << shapesParam4[i]<<std::endl;
234 
235 }
236 
237 
246 unsigned int checkAndRetrunIndex(const std::string &shapeName)
247 {
248  unsigned int pos=0;
249 
250  while ((pos < shapes2D.size()) && (shapes2D[pos] != shapeName))
251  pos++;
252 
253  if (pos == shapes2D.size())
254  {
255  trace.error() << "The specified shape has not found.";
256  trace.info()<<std::endl;
257  exit(1);
258  }
259 
260  return pos;
261 }
262 
268 void missingParam(std::string param)
269 {
270  trace.error() <<" Parameter: "<<param<<" is required..";
271  trace.info()<<std::endl;
272  exit(1);
273 }
275 
276 template <typename Shape, typename Space>
277 bool
278 lengthEstimators( const std::string & /*name*/,
279  Shape & aShape,
280  double h )
281 {
282  // Types
283  typedef typename Space::Point Point;
284  typedef typename Space::Vector Vector;
285  typedef typename Space::RealPoint RealPoint;
286  typedef typename Space::Integer Integer;
287  typedef HyperRectDomain<Space> Domain;
288  typedef KhalimskySpaceND<Space::dimension,Integer> KSpace;
289  typedef typename KSpace::SCell SCell;
290  typedef typename GridCurve<KSpace>::PointsRange PointsRange;
291  typedef typename GridCurve<KSpace>::ArrowsRange ArrowsRange;
292 
293  // Digitizer
294  GaussDigitizer<Space,Shape> dig;
295  dig.attach( aShape ); // attaches the shape.
296  Vector vlow(-1,-1); Vector vup(1,1);
297  dig.init( aShape.getLowerBound()+vlow, aShape.getUpperBound()+vup, h );
298  Domain domain = dig.getDomain();
299 
300  // Create cellular space
301  KSpace K;
302  bool ok = K.init( dig.getLowerBound(), dig.getUpperBound(), true );
303  if ( ! ok )
304  {
305  std::cerr << "[lengthEstimators]"
306  << " error in creating KSpace." << std::endl;
307  return false;
308  }
309  try {
310  // Extracts shape boundary
311  SurfelAdjacency<KSpace::dimension> SAdj( true );
312  SCell bel = Surfaces<KSpace>::findABel( K, dig, 10000 );
313  // Getting the consecutive surfels of the 2D boundary
314  std::vector<Point> points;
315  Surfaces<KSpace>::track2DBoundaryPoints( points, K, SAdj, dig, bel );
316  // Create GridCurve
317  GridCurve<KSpace> gridcurve;
318  gridcurve.initFromVector( points );
319  // Ranges
320  ArrowsRange ra = gridcurve.getArrowsRange();
321  PointsRange rp = gridcurve.getPointsRange();
322 
323 
324  // Estimations
325  typedef typename PointsRange::ConstIterator ConstIteratorOnPoints;
326  typedef ParametricShapeArcLengthFunctor< Shape > Length;
327  TrueGlobalEstimatorOnPoints< ConstIteratorOnPoints, Shape, Length > trueLengthEstimator;
328  trueLengthEstimator.init( h, rp.begin(), rp.end(), &aShape, gridcurve.isClosed());
329 
330  L1LengthEstimator< typename ArrowsRange::ConstCirculator > l1length;
331  DSSLengthEstimator< typename PointsRange::ConstCirculator > DSSlength;
332  MLPLengthEstimator< typename PointsRange::ConstIterator > MLPlength;
333  FPLengthEstimator< typename PointsRange::ConstIterator > FPlength;
334  BLUELocalLengthEstimator< typename ArrowsRange::ConstIterator > BLUElength;
335  RosenProffittLocalLengthEstimator< typename ArrowsRange::ConstIterator > RosenProffittlength;
336 
337  // Output
338  double trueValue = trueLengthEstimator.eval();
339  double l1, blue, rosen,dss,mlp,fp;
340  double Tl1, Tblue, Trosen,Tdss,Tmlp,Tfp;
341 
342  Clock c;
343 
344  //Length evaluation & timing
345  c.startClock();
346  l1length.init(h, ra.c(), ra.c());
347  l1 = l1length.eval();
348  Tl1 = c.stopClock();
349 
350  c.startClock();
351  BLUElength.init(h, ra.begin(), ra.end(), gridcurve.isClosed());
352  blue = BLUElength.eval();
353  Tblue = c.stopClock();
354 
355  c.startClock();
356  RosenProffittlength.init(h, ra.begin(), ra.end(), gridcurve.isClosed());
357  rosen = RosenProffittlength.eval();
358  Trosen = c.stopClock();
359 
360  c.startClock();
361  DSSlength.init(h, rp.c(), rp.c());
362  dss = DSSlength.eval();
363  Tdss = c.stopClock();
364 
365  c.startClock();
366  MLPlength.init(h, rp.begin(), rp.end(), gridcurve.isClosed());
367  mlp = MLPlength.eval();
368  Tmlp = c.stopClock();
369 
370  c.startClock();
371  FPlength.init(h, rp.begin(), rp.end(), gridcurve.isClosed());
372  fp = FPlength.eval();
373  Tfp = c.stopClock();
374 
375  std::cout << std::setprecision( 15 ) << h << " " << rp.size() << " " << trueValue
376  << " " << l1
377  << " " << blue
378  << " " << rosen
379  << " " << dss
380  << " " << mlp
381  << " " << fp
382  << " " << Tl1
383  << " " << Tblue
384  << " " << Trosen
385  << " " << Tdss
386  << " " << Tmlp
387  << " " << Tfp
388  << std::endl;
389  return true;
390  }
391  catch ( InputException e )
392  {
393  std::cerr << "[lengthEstimators]"
394  << " error in finding a bel." << std::endl;
395  return false;
396  }
397 }
399 namespace po = boost::program_options;
400 
401 int main( int argc, char** argv )
402 {
403  // parse command line ----------------------------------------------
404  po::options_description general_opt("Allowed options are: ");
405  general_opt.add_options()
406  ("help,h", "display this message")
407  ("shape,s", po::value<std::string>(), "Shape name")
408  ("list,l", "List all available shapes")
409  ("radius,R", po::value<double>(), "Radius of the shape" )
410  ("axis1,A", po::value<double>(), "Half big axis of the shape (ellipse)" )
411  ("axis2,a", po::value<double>(), "Half small axis of the shape (ellipse)" )
412  ("smallradius,r", po::value<double>()->default_value(5), "Small radius of the shape" )
413  ("varsmallradius,v", po::value<double>()->default_value(5), "Variable small radius of the shape" )
414  ("k,k", po::value<unsigned int>()->default_value(3), "Number of branches or corners the shape" )
415  ("phi", po::value<double>()->default_value(0.0), "Phase of the shape (in radian)" )
416  ("width,w", po::value<double>()->default_value(10.0), "Width of the shape" )
417  ("power,p", po::value<double>()->default_value(2.0), "Power of the metric (double)" )
418  ("hMin", po::value<double>()->default_value(0.0001), "Minimum value for the grid step h (double)" )
419  ("steps", po::value<int>()->default_value(32), "Number of multigrid steps between 1 and hMin (integer)" );
420 
421  bool parseOK=true;
422  po::variables_map vm;
423  try{
424  po::store(po::parse_command_line(argc, argv, general_opt), vm);
425  }catch(const std::exception& ex){
426  parseOK=false;
427  trace.info()<< "Error checking program options: "<< ex.what()<< std::endl;
428  }
429  po::notify(vm);
430  if(!parseOK || vm.count("help")||argc<=1)
431  {
432  trace.info()<< "Generate multigrid length estimations of paramteric shapes using DGtal library. It will output length estimations (and timings) using several algorithms for decreasing grid steps." <<std::endl << "Basic usage: "<<std::endl
433  << "\tLengthEstimators [options] --shape <shapeName>"<<std::endl
434  << general_opt << "\n";
435  return 0;
436  }
437 
438  //List creation
439  createList();
440 
441  if (vm.count("list"))
442  {
443  displayList();
444  return 0;
445  }
446 
447  //Parse options
448  if (!(vm.count("shape"))) missingParam("--shape");
449  std::string shapeName = vm["shape"].as<std::string>();
450  double hMin = vm["hMin"].as<double>();
451  int nbSteps = vm["steps"].as<int>();
452 
453  //We check that the shape is known
454  unsigned int id = checkAndRetrunIndex(shapeName);
455 
456 
458  std::cout << "#h nbPoints true-length L1 BLUE RosenProffit "
459  << "DSS MLP FP Time-L1 Time-BLUE Time-RosenProffitt "
460  << "Time-DSS Time-MLP Time-FP" <<std::endl;
461  std::cout << "# timings are given in msec." <<std::endl;
462 
463  double h = 1;
464  double step = exp( log(hMin) / (double)nbSteps);
465  while (h > hMin) {
466 
467  if (id ==0)
468  {
469  if (!(vm.count("radius"))) missingParam("--radius");
470  double radius = vm["radius"].as<double>();
471 
472  Ball2D<Z2i::Space> ball(Z2i::Point(0,0), radius);
473 
474  lengthEstimators<Ball2D<Z2i::Space>,Z2i::Space>("ball",ball,h);
475  }
476  else
477  if (id ==1)
478  {
479  if (!(vm.count("width"))) missingParam("--width");
480  double width = vm["width"].as<double>();
481 
482  ImplicitHyperCube<Z2i::Space> object(Z2i::Point(0,0), width/2);
483 
484  trace.error()<< "Not available.";
485  trace.info()<<std::endl;
486  }
487  else
488  if (id ==2)
489  {
490  if (!(vm.count("power"))) missingParam("--power");
491  if (!(vm.count("radius"))) missingParam("--radius");
492  double radius = vm["radius"].as<double>();
493  double power = vm["power"].as<double>();
494 
495  ImplicitRoundedHyperCube<Z2i::Space> ball(Z2i::Point(0,0), radius, power);
496 
497  trace.error()<< "Not available.";
498  trace.info()<<std::endl;
499  }
500  else
501  if (id ==3)
502  {
503  if (!(vm.count("varsmallradius"))) missingParam("--varsmallradius");
504  if (!(vm.count("radius"))) missingParam("--radius");
505  if (!(vm.count("k"))) missingParam("--k");
506  if (!(vm.count("phi"))) missingParam("--phi");
507  double radius = vm["radius"].as<double>();
508  double varsmallradius = vm["varsmallradius"].as<double>();
509  unsigned int k = vm["k"].as<unsigned int>();
510  double phi = vm["phi"].as<double>();
511 
512  Flower2D<Z2i::Space> flower(Z2i::Point(0,0), radius, varsmallradius,k,phi);
513 
514  lengthEstimators<Flower2D<Z2i::Space>,Z2i::Space>("flower",flower,h);
515  }
516  else
517  if (id ==4)
518  {
519  if (!(vm.count("radius"))) missingParam("--radius");
520  if (!(vm.count("k"))) missingParam("--k");
521  if (!(vm.count("phi"))) missingParam("--phi");
522  double radius = vm["radius"].as<double>();
523  unsigned int k = vm["k"].as<unsigned int>();
524  double phi = vm["phi"].as<double>();
525 
526  NGon2D<Z2i::Space> object(Z2i::Point(0,0), radius,k,phi);
527 
528  lengthEstimators<NGon2D<Z2i::Space>,Z2i::Space>("NGon",object,h);
529 
530  }
531  else
532  if (id ==5)
533  {
534  if (!(vm.count("varsmallradius"))) missingParam("--varsmallradius");
535  if (!(vm.count("radius"))) missingParam("--radius");
536  if (!(vm.count("k"))) missingParam("--k");
537  if (!(vm.count("phi"))) missingParam("--phi");
538  double radius = vm["radius"].as<double>();
539  double varsmallradius = vm["varsmallradius"].as<double>();
540  unsigned int k = vm["k"].as<unsigned int>();
541  double phi = vm["phi"].as<double>();
542 
543  AccFlower2D<Z2i::Space> flower(Z2i::Point(0,0), radius, varsmallradius,k,phi);
544  lengthEstimators<AccFlower2D<Z2i::Space>,Z2i::Space>("accFlower",flower,h);
545 
546  }
547  else
548  //if (id ==6)
549  {
550  if (!(vm.count("axis1"))) missingParam("--axis1");
551  if (!(vm.count("axis2"))) missingParam("--axis2");
552  if (!(vm.count("phi"))) missingParam("--phi");
553  double a1 = vm["axis1"].as<double>();
554  double a2 = vm["axis2"].as<double>();
555  double phi = vm["phi"].as<double>();
556 
557  Ellipse2D<Z2i::Space> ell(Z2i::Point(0,0), a1, a2,phi);
558 
559  lengthEstimators<Ellipse2D<Z2i::Space>,Z2i::Space>("Ellipse",ell,h);
560 
561  }
562 
563  h = h * step;
564  }
565  return 0;
566 }
Companion project associated to DGtal ; Generated on Mon Jun 27 2016 12:21:49 for DGtalTools by   doxygen 1.8.11