DGtalTools  1.2.0
lengthEstimators.cpp
1 
33 #include <cmath>
34 #include <iostream>
35 #include <iomanip>
36 #include <vector>
37 #include <string>
38 
39 #include "CLI11.hpp"
40 
41 #include "DGtal/base/Common.h"
42 #include "DGtal/base/Clock.h"
43 
44 //space / domain
45 #include "DGtal/kernel/SpaceND.h"
46 #include "DGtal/kernel/domains/HyperRectDomain.h"
47 #include "DGtal/topology/KhalimskySpaceND.h"
48 
49 //shape and digitizer
50 #include "DGtal/shapes/ShapeFactory.h"
51 #include "DGtal/shapes/Shapes.h"
52 #include "DGtal/helpers/StdDefs.h"
53 #include "DGtal/topology/helpers/Surfaces.h"
54 
55 #include "DGtal/shapes/GaussDigitizer.h"
56 #include "DGtal/geometry/curves/GridCurve.h"
57 
58 //estimators
59 #include "DGtal/geometry/curves/estimation/TrueLocalEstimatorOnPoints.h"
60 #include "DGtal/geometry/curves/estimation/TrueGlobalEstimatorOnPoints.h"
61 #include "DGtal/geometry/curves/estimation/ParametricShapeArcLengthFunctor.h"
62 
63 #include "DGtal/geometry/curves/estimation/L1LengthEstimator.h"
64 #include "DGtal/geometry/curves/estimation/BLUELocalLengthEstimator.h"
65 #include "DGtal/geometry/curves/estimation/RosenProffittLocalLengthEstimator.h"
66 #include "DGtal/geometry/curves/estimation/MLPLengthEstimator.h"
67 #include "DGtal/geometry/curves/estimation/FPLengthEstimator.h"
68 #include "DGtal/geometry/curves/estimation/DSSLengthEstimator.h"
69 
70 using namespace DGtal;
71 
72 
73 
149 std::vector<std::string> shapes2D;
150 std::vector<std::string> shapesDesc;
151 std::vector<std::string> shapesParam1;
152 std::vector<std::string> shapesParam2;
153 std::vector<std::string> shapesParam3;
154 std::vector<std::string> shapesParam4;
155 
156 
161 void createList()
162 {
163  shapes2D.push_back("ball");
164  shapesDesc.push_back("Ball for the Euclidean metric.");
165  shapesParam1.push_back("--radius [-R]");
166  shapesParam2.push_back("");
167  shapesParam3.push_back("");
168  shapesParam4.push_back("");
169 
170  shapes2D.push_back("square");
171  shapesDesc.push_back("square (no signature).");
172  shapesParam1.push_back("--width [-w]");
173  shapesParam2.push_back("");
174  shapesParam3.push_back("");
175  shapesParam4.push_back("");
176 
177  shapes2D.push_back("lpball");
178  shapesDesc.push_back("Ball for the l_power metric (no signature).");
179  shapesParam1.push_back("--radius [-R],");
180  shapesParam2.push_back("--power [-p]");
181  shapesParam3.push_back("");
182  shapesParam4.push_back("");
183 
184  shapes2D.push_back("flower");
185  shapesDesc.push_back("Flower with k petals.");
186  shapesParam1.push_back("--radius [-R],");
187  shapesParam2.push_back("--varsmallradius [-v],");
188  shapesParam3.push_back("--k [-k],");
189  shapesParam4.push_back("--phi");
190 
191  shapes2D.push_back("ngon");
192  shapesDesc.push_back("Regular k-gon.");
193  shapesParam1.push_back("--radius [-R],");
194  shapesParam2.push_back("--k [-k],");
195  shapesParam3.push_back("--phi");
196  shapesParam4.push_back("");
197 
198  shapes2D.push_back("accflower");
199  shapesDesc.push_back("Accelerated Flower with k petals.");
200  shapesParam1.push_back("--radius [-R],");
201  shapesParam2.push_back("--varsmallradius [-v],");
202  shapesParam3.push_back("--k [-k],");
203  shapesParam4.push_back("--phi");
204 
205  shapes2D.push_back("ellipse");
206  shapesDesc.push_back("Ellipse.");
207  shapesParam1.push_back("--axis1 [-A],");
208  shapesParam2.push_back("--axis2 [-a],");
209  shapesParam3.push_back("--phi");
210  shapesParam4.push_back("");
211 
212 
213 }
214 
219 void displayList()
220 {
221  trace.emphase()<<"2D Shapes:"<<std::endl;
222  for(unsigned int i=0; i<shapes2D.size(); ++i)
223  trace.info()<<"\t"<<shapes2D[i]<<"\t"
224  << shapesDesc[i]<<std::endl
225  << "\t\tRequired parameter(s): "
226  << shapesParam1[i]<<" "
227  << shapesParam2[i]<<" "
228  << shapesParam3[i]<<" "
229  << shapesParam4[i]<<std::endl;
230 
231 }
232 
233 
242 unsigned int checkAndRetrunIndex(const std::string &shapeName)
243 {
244  unsigned int pos=0;
245 
246  while ((pos < shapes2D.size()) && (shapes2D[pos] != shapeName))
247  pos++;
248 
249  if (pos == shapes2D.size())
250  {
251  trace.error() << "The specified shape has not found.";
252  trace.info()<<std::endl;
253  exit(1);
254  }
255 
256  return pos;
257 }
258 
264 void missingParam(std::string param)
265 {
266  trace.error() <<" Parameter: "<<param<<" is required..";
267  trace.info()<<std::endl;
268  exit(1);
269 }
271 
272 template <typename Shape, typename Space>
273 bool
274 lengthEstimators( const std::string & /*name*/,
275  Shape & aShape,
276  double h )
277 {
278  // Types
279  typedef typename Space::Point Point;
280  typedef typename Space::Vector Vector;
281  typedef typename Space::Integer Integer;
282  typedef HyperRectDomain<Space> Domain;
283  typedef KhalimskySpaceND<Space::dimension,Integer> KSpace;
284  typedef typename KSpace::SCell SCell;
285  typedef typename GridCurve<KSpace>::PointsRange PointsRange;
286  typedef typename GridCurve<KSpace>::ArrowsRange ArrowsRange;
287 
288  // Digitizer
289  GaussDigitizer<Space,Shape> dig;
290  dig.attach( aShape ); // attaches the shape.
291  Vector vlow(-1,-1); Vector vup(1,1);
292  dig.init( aShape.getLowerBound()+vlow, aShape.getUpperBound()+vup, h );
293  Domain domain = dig.getDomain();
294 
295  // Create cellular space
296  KSpace K;
297  bool ok = K.init( dig.getLowerBound(), dig.getUpperBound(), true );
298  if ( ! ok )
299  {
300  std::cerr << "[lengthEstimators]"
301  << " error in creating KSpace." << std::endl;
302  return false;
303  }
304  try
305  {
306  // Extracts shape boundary
307  SurfelAdjacency<KSpace::dimension> SAdj( true );
308  SCell bel = Surfaces<KSpace>::findABel( K, dig, 10000 );
309  // Getting the consecutive surfels of the 2D boundary
310  std::vector<Point> points;
311  Surfaces<KSpace>::track2DBoundaryPoints( points, K, SAdj, dig, bel );
312  // Create GridCurve
313  GridCurve<KSpace> gridcurve;
314  gridcurve.initFromVector( points );
315  // Ranges
316  ArrowsRange ra = gridcurve.getArrowsRange();
317  PointsRange rp = gridcurve.getPointsRange();
318 
319  // Estimations
320  typedef typename PointsRange::ConstIterator ConstIteratorOnPoints;
321  typedef ParametricShapeArcLengthFunctor< Shape > Length;
322  TrueGlobalEstimatorOnPoints< ConstIteratorOnPoints, Shape, Length > trueLengthEstimator;
323  trueLengthEstimator.attach( aShape );
324 
325  L1LengthEstimator< typename PointsRange::ConstCirculator > l1length;
326  DSSLengthEstimator< typename PointsRange::ConstCirculator > DSSlength;
327  MLPLengthEstimator< typename PointsRange::ConstCirculator > MLPlength;
328  FPLengthEstimator< typename PointsRange::ConstCirculator > FPlength;
329  BLUELocalLengthEstimator< typename ArrowsRange::ConstCirculator > BLUElength;
330  RosenProffittLocalLengthEstimator< typename ArrowsRange::ConstCirculator > RosenProffittlength;
331 
332  // Output
333  double trueValue = trueLengthEstimator.eval();
334  double l1, blue, rosen,dss,mlp,fp;
335  double Tl1, Tblue, Trosen,Tdss,Tmlp,Tfp;
336 
337  Clock c;
338 
339  //Length evaluation & timing
340  c.startClock();
341  l1 = l1length.eval( rp.c(), rp.c(), h );
342  Tl1 = c.stopClock();
343 
344  c.startClock();
345  blue = BLUElength.eval( ra.c(), ra.c(), h );
346  Tblue = c.stopClock();
347 
348  c.startClock();
349  rosen = RosenProffittlength.eval( ra.c(), ra.c(), h );
350  Trosen = c.stopClock();
351 
352  c.startClock();
353  dss = DSSlength.eval( rp.c(), rp.c(), h );
354  Tdss = c.stopClock();
355 
356  c.startClock();
357  mlp = MLPlength.eval( rp.c(), rp.c(), h );
358  Tmlp = c.stopClock();
359 
360  c.startClock();
361  fp = FPlength.eval( rp.c(), rp.c(), h );
362  Tfp = c.stopClock();
363 
364  std::cout << std::setprecision( 15 ) << h << " " << rp.size() << " " << trueValue
365  << " " << l1
366  << " " << blue
367  << " " << rosen
368  << " " << dss
369  << " " << mlp
370  << " " << fp
371  << " " << Tl1
372  << " " << Tblue
373  << " " << Trosen
374  << " " << Tdss
375  << " " << Tmlp
376  << " " << Tfp
377  << std::endl;
378  return true;
379  }
380  catch ( InputException e )
381  {
382  std::cerr << "[lengthEstimators]"
383  << " error in finding a bel." << std::endl;
384  return false;
385  }
386 }
388 
389 int main( int argc, char** argv )
390 {
391  // parse command line CLI ----------------------------------------------
392  CLI::App app;
393  std::string shapeName;
394  double hMin {0.0001};
395  int nbSteps {32};
396  double radius;
397  double power {2.0};
398  double smallradius {5};
399  double varsmallradius {5};
400  unsigned int k {3};
401  double phi {0.0};
402  double width {10.0};
403  double axis1, axis2;
404 
405  app.description("Generates multigrid contours of 2d digital shapes using DGtal library.\n Typical use example:\n \t contourGenerator --shape <shapeName> [requiredParam] [otherOptions]\n");
406  auto listOpt = app.add_flag("--list,-l","List all available shapes");
407  auto shapeNameOpt = app.add_option("--shape,-s", shapeName, "Shape name");
408  auto radiusOpt = app.add_option("--radius,-R", radius, "Radius of the shape" );
409  auto axis1Opt = app.add_option("--axis1,-A", axis1, "Half big axis of the shape (ellipse)" );
410  auto axis2Opt = app.add_option("--axis2,-a", axis2, "Half small axis of the shape (ellipse)" );
411  auto smallradiusOpt = app.add_option("--smallradius,-r", smallradius, "Small radius of the shape (default 5)", true);
412  auto varsmallradiusOpt = app.add_option("--varsmallradius,-v", varsmallradius, "Variable small radius of the shape (default 5)", true );
413  auto kOpt = app.add_option("-k", k, "Number of branches or corners the shape (default 3)", true );
414  auto phiOpt = app.add_option("--phi", phi, "Phase of the shape (in radian, default 0.0)", true );
415  auto widthOpt = app.add_option("--width,-w", width, "Width of the shape (default 10.0)", true );
416  auto powerOpt = app.add_option("--power,-p", power, "Power of the metric (default 2.0)", true );
417  app.add_option("--hMin", hMin, "Minimum value for the grid step h (double, default 0.0001)", true );
418  app.add_option("--steps", nbSteps, "Number of multigrid steps between 1 and hMin (integer, default 32)", true );
419 
420  app.get_formatter()->column_width(40);
421  CLI11_PARSE(app, argc, argv);
422  // END parse command line using CLI ----------------------------------------------
423 
424  //List creation
425  createList();
426 
427  if ( listOpt->count() > 0 )
428  {
429  displayList();
430  return 0;
431  }
432 
433  //Parse options
434  if(shapeNameOpt->count()==0) missingParam("--shape");
435 
436  //We check that the shape is known
437  unsigned int id = checkAndRetrunIndex(shapeName);
438 
439 
441  std::cout << "#h nbPoints true-length L1 BLUE RosenProffit "
442  << "DSS MLP FP Time-L1 Time-BLUE Time-RosenProffitt "
443  << "Time-DSS Time-MLP Time-FP" << std::endl;
444  std::cout << "# timings are given in msec." << std::endl;
445 
446  double h = 1;
447  double step = exp( log(hMin) / (double)nbSteps);
448  while (h > hMin) {
449 
450  if (id ==0)
451  {
452  if (radiusOpt->count()==0) missingParam("--radius");
453 
454  Ball2D<Z2i::Space> ball(Z2i::Point(0,0), radius);
455 
456  lengthEstimators<Ball2D<Z2i::Space>,Z2i::Space>("ball",ball,h);
457  }
458  else
459  if (id ==1)
460  {
461  //if (widthOpt->count()==0) missingParam("--width");
462 
463  ImplicitHyperCube<Z2i::Space> object(Z2i::Point(0,0), width/2);
464 
465  trace.error()<< "Not available.";
466  trace.info()<<std::endl;
467  }
468  else
469  if (id ==2)
470  {
471  //if (powerOpt->count()==0) missingParam("--power");
472  if (radiusOpt->count()==0) missingParam("--radius");
473 
474  ImplicitRoundedHyperCube<Z2i::Space> ball(Z2i::Point(0,0), radius, power);
475 
476  trace.error()<< "Not available.";
477  trace.info()<<std::endl;
478  }
479  else
480  if (id ==3)
481  {
482  //if (varsmallradiusOpt->count()==0) missingParam("--varsmallradius");
483  if (radiusOpt->count()==0) missingParam("--radius");
484  //if (kOpt->count()==0) missingParam("--k");
485  //if (phiOpt->count()==0) missingParam("--phi");
486 
487  Flower2D<Z2i::Space> flower(Z2i::Point(0,0), radius, varsmallradius,k,phi);
488 
489  lengthEstimators<Flower2D<Z2i::Space>,Z2i::Space>("flower",flower,h);
490  }
491  else
492  if (id ==4)
493  {
494  if (radiusOpt->count()==0) missingParam("--radius");
495  //if (kOpt->count()==0) missingParam("--k");
496  //if (phiOpt->count()==0) missingParam("--phi");
497 
498  NGon2D<Z2i::Space> object(Z2i::Point(0,0), radius,k,phi);
499 
500  lengthEstimators<NGon2D<Z2i::Space>,Z2i::Space>("NGon",object,h);
501 
502  }
503  else
504  if (id ==5)
505  {
506  //if (varsmallradiusOpt->count()==0) missingParam("--varsmallradius");
507  if (radiusOpt->count()==0) missingParam("--radius");
508  //if (kOpt->count()==0) missingParam("--k");
509  //if (phiOpt->count()==0) missingParam("--phi");
510 
511  AccFlower2D<Z2i::Space> flower(Z2i::Point(0,0), radius, varsmallradius,k,phi);
512  lengthEstimators<AccFlower2D<Z2i::Space>,Z2i::Space>("accFlower",flower,h);
513 
514  }
515  else
516  //if (id ==6)
517  {
518  if (axis1Opt->count()==0) missingParam("--axis1");
519  if (axis2Opt->count()==0) missingParam("--axis2");
520  //if (phiOpt->count()==0) missingParam("--phi");
521 
522  Ellipse2D<Z2i::Space> ell(Z2i::Point(0,0), axis1, axis2,phi);
523 
524  lengthEstimators<Ellipse2D<Z2i::Space>,Z2i::Space>("Ellipse",ell,h);
525  }
526 
527  h = h * step;
528  }
529  return 0;
530 }
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