DGtalTools  0.9.3
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::Integer Integer;
286  typedef HyperRectDomain<Space> Domain;
287  typedef KhalimskySpaceND<Space::dimension,Integer> KSpace;
288  typedef typename KSpace::SCell SCell;
289  typedef typename GridCurve<KSpace>::PointsRange PointsRange;
290  typedef typename GridCurve<KSpace>::ArrowsRange ArrowsRange;
291 
292  // Digitizer
293  GaussDigitizer<Space,Shape> dig;
294  dig.attach( aShape ); // attaches the shape.
295  Vector vlow(-1,-1); Vector vup(1,1);
296  dig.init( aShape.getLowerBound()+vlow, aShape.getUpperBound()+vup, h );
297  Domain domain = dig.getDomain();
298 
299  // Create cellular space
300  KSpace K;
301  bool ok = K.init( dig.getLowerBound(), dig.getUpperBound(), true );
302  if ( ! ok )
303  {
304  std::cerr << "[lengthEstimators]"
305  << " error in creating KSpace." << std::endl;
306  return false;
307  }
308  try {
309  // Extracts shape boundary
310  SurfelAdjacency<KSpace::dimension> SAdj( true );
311  SCell bel = Surfaces<KSpace>::findABel( K, dig, 10000 );
312  // Getting the consecutive surfels of the 2D boundary
313  std::vector<Point> points;
314  Surfaces<KSpace>::track2DBoundaryPoints( points, K, SAdj, dig, bel );
315  // Create GridCurve
316  GridCurve<KSpace> gridcurve;
317  gridcurve.initFromVector( points );
318  // Ranges
319  ArrowsRange ra = gridcurve.getArrowsRange();
320  PointsRange rp = gridcurve.getPointsRange();
321 
322 
323  // Estimations
324  typedef typename PointsRange::ConstIterator ConstIteratorOnPoints;
325  typedef ParametricShapeArcLengthFunctor< Shape > Length;
326  TrueGlobalEstimatorOnPoints< ConstIteratorOnPoints, Shape, Length > trueLengthEstimator;
327  trueLengthEstimator.init( h, rp.begin(), rp.end(), &aShape, gridcurve.isClosed());
328 
329  L1LengthEstimator< typename ArrowsRange::ConstCirculator > l1length;
330  DSSLengthEstimator< typename PointsRange::ConstCirculator > DSSlength;
331  MLPLengthEstimator< typename PointsRange::ConstIterator > MLPlength;
332  FPLengthEstimator< typename PointsRange::ConstIterator > FPlength;
333  BLUELocalLengthEstimator< typename ArrowsRange::ConstIterator > BLUElength;
334  RosenProffittLocalLengthEstimator< typename ArrowsRange::ConstIterator > RosenProffittlength;
335 
336  // Output
337  double trueValue = trueLengthEstimator.eval();
338  double l1, blue, rosen,dss,mlp,fp;
339  double Tl1, Tblue, Trosen,Tdss,Tmlp,Tfp;
340 
341  Clock c;
342 
343  //Length evaluation & timing
344  c.startClock();
345  l1length.init(h, ra.c(), ra.c());
346  l1 = l1length.eval();
347  Tl1 = c.stopClock();
348 
349  c.startClock();
350  BLUElength.init(h, ra.begin(), ra.end(), gridcurve.isClosed());
351  blue = BLUElength.eval();
352  Tblue = c.stopClock();
353 
354  c.startClock();
355  RosenProffittlength.init(h, ra.begin(), ra.end(), gridcurve.isClosed());
356  rosen = RosenProffittlength.eval();
357  Trosen = c.stopClock();
358 
359  c.startClock();
360  DSSlength.init(h, rp.c(), rp.c());
361  dss = DSSlength.eval();
362  Tdss = c.stopClock();
363 
364  c.startClock();
365  MLPlength.init(h, rp.begin(), rp.end(), gridcurve.isClosed());
366  mlp = MLPlength.eval();
367  Tmlp = c.stopClock();
368 
369  c.startClock();
370  FPlength.init(h, rp.begin(), rp.end(), gridcurve.isClosed());
371  fp = FPlength.eval();
372  Tfp = c.stopClock();
373 
374  std::cout << std::setprecision( 15 ) << h << " " << rp.size() << " " << trueValue
375  << " " << l1
376  << " " << blue
377  << " " << rosen
378  << " " << dss
379  << " " << mlp
380  << " " << fp
381  << " " << Tl1
382  << " " << Tblue
383  << " " << Trosen
384  << " " << Tdss
385  << " " << Tmlp
386  << " " << Tfp
387  << std::endl;
388  return true;
389  }
390  catch ( InputException e )
391  {
392  std::cerr << "[lengthEstimators]"
393  << " error in finding a bel." << std::endl;
394  return false;
395  }
396 }
398 namespace po = boost::program_options;
399 
400 int main( int argc, char** argv )
401 {
402  // parse command line ----------------------------------------------
403  po::options_description general_opt("Allowed options are: ");
404  general_opt.add_options()
405  ("help,h", "display this message")
406  ("shape,s", po::value<std::string>(), "Shape name")
407  ("list,l", "List all available shapes")
408  ("radius,R", po::value<double>(), "Radius of the shape" )
409  ("axis1,A", po::value<double>(), "Half big axis of the shape (ellipse)" )
410  ("axis2,a", po::value<double>(), "Half small axis of the shape (ellipse)" )
411  ("smallradius,r", po::value<double>()->default_value(5), "Small radius of the shape" )
412  ("varsmallradius,v", po::value<double>()->default_value(5), "Variable small radius of the shape" )
413  ("k,k", po::value<unsigned int>()->default_value(3), "Number of branches or corners the shape" )
414  ("phi", po::value<double>()->default_value(0.0), "Phase of the shape (in radian)" )
415  ("width,w", po::value<double>()->default_value(10.0), "Width of the shape" )
416  ("power,p", po::value<double>()->default_value(2.0), "Power of the metric (double)" )
417  ("hMin", po::value<double>()->default_value(0.0001), "Minimum value for the grid step h (double)" )
418  ("steps", po::value<int>()->default_value(32), "Number of multigrid steps between 1 and hMin (integer)" );
419 
420  bool parseOK=true;
421  po::variables_map vm;
422  try{
423  po::store(po::parse_command_line(argc, argv, general_opt), vm);
424  }catch(const std::exception& ex){
425  parseOK=false;
426  trace.info()<< "Error checking program options: "<< ex.what()<< std::endl;
427  }
428  po::notify(vm);
429  if(!parseOK || vm.count("help")||argc<=1)
430  {
431  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
432  << "\tLengthEstimators [options] --shape <shapeName>"<<std::endl
433  << general_opt << "\n";
434  return 0;
435  }
436 
437  //List creation
438  createList();
439 
440  if (vm.count("list"))
441  {
442  displayList();
443  return 0;
444  }
445 
446  //Parse options
447  if (!(vm.count("shape"))) missingParam("--shape");
448  std::string shapeName = vm["shape"].as<std::string>();
449  double hMin = vm["hMin"].as<double>();
450  int nbSteps = vm["steps"].as<int>();
451 
452  //We check that the shape is known
453  unsigned int id = checkAndRetrunIndex(shapeName);
454 
455 
457  std::cout << "#h nbPoints true-length L1 BLUE RosenProffit "
458  << "DSS MLP FP Time-L1 Time-BLUE Time-RosenProffitt "
459  << "Time-DSS Time-MLP Time-FP" <<std::endl;
460  std::cout << "# timings are given in msec." <<std::endl;
461 
462  double h = 1;
463  double step = exp( log(hMin) / (double)nbSteps);
464  while (h > hMin) {
465 
466  if (id ==0)
467  {
468  if (!(vm.count("radius"))) missingParam("--radius");
469  double radius = vm["radius"].as<double>();
470 
471  Ball2D<Z2i::Space> ball(Z2i::Point(0,0), radius);
472 
473  lengthEstimators<Ball2D<Z2i::Space>,Z2i::Space>("ball",ball,h);
474  }
475  else
476  if (id ==1)
477  {
478  if (!(vm.count("width"))) missingParam("--width");
479  double width = vm["width"].as<double>();
480 
481  ImplicitHyperCube<Z2i::Space> object(Z2i::Point(0,0), width/2);
482 
483  trace.error()<< "Not available.";
484  trace.info()<<std::endl;
485  }
486  else
487  if (id ==2)
488  {
489  if (!(vm.count("power"))) missingParam("--power");
490  if (!(vm.count("radius"))) missingParam("--radius");
491  double radius = vm["radius"].as<double>();
492  double power = vm["power"].as<double>();
493 
494  ImplicitRoundedHyperCube<Z2i::Space> ball(Z2i::Point(0,0), radius, power);
495 
496  trace.error()<< "Not available.";
497  trace.info()<<std::endl;
498  }
499  else
500  if (id ==3)
501  {
502  if (!(vm.count("varsmallradius"))) missingParam("--varsmallradius");
503  if (!(vm.count("radius"))) missingParam("--radius");
504  if (!(vm.count("k"))) missingParam("--k");
505  if (!(vm.count("phi"))) missingParam("--phi");
506  double radius = vm["radius"].as<double>();
507  double varsmallradius = vm["varsmallradius"].as<double>();
508  unsigned int k = vm["k"].as<unsigned int>();
509  double phi = vm["phi"].as<double>();
510 
511  Flower2D<Z2i::Space> flower(Z2i::Point(0,0), radius, varsmallradius,k,phi);
512 
513  lengthEstimators<Flower2D<Z2i::Space>,Z2i::Space>("flower",flower,h);
514  }
515  else
516  if (id ==4)
517  {
518  if (!(vm.count("radius"))) missingParam("--radius");
519  if (!(vm.count("k"))) missingParam("--k");
520  if (!(vm.count("phi"))) missingParam("--phi");
521  double radius = vm["radius"].as<double>();
522  unsigned int k = vm["k"].as<unsigned int>();
523  double phi = vm["phi"].as<double>();
524 
525  NGon2D<Z2i::Space> object(Z2i::Point(0,0), radius,k,phi);
526 
527  lengthEstimators<NGon2D<Z2i::Space>,Z2i::Space>("NGon",object,h);
528 
529  }
530  else
531  if (id ==5)
532  {
533  if (!(vm.count("varsmallradius"))) missingParam("--varsmallradius");
534  if (!(vm.count("radius"))) missingParam("--radius");
535  if (!(vm.count("k"))) missingParam("--k");
536  if (!(vm.count("phi"))) missingParam("--phi");
537  double radius = vm["radius"].as<double>();
538  double varsmallradius = vm["varsmallradius"].as<double>();
539  unsigned int k = vm["k"].as<unsigned int>();
540  double phi = vm["phi"].as<double>();
541 
542  AccFlower2D<Z2i::Space> flower(Z2i::Point(0,0), radius, varsmallradius,k,phi);
543  lengthEstimators<AccFlower2D<Z2i::Space>,Z2i::Space>("accFlower",flower,h);
544 
545  }
546  else
547  //if (id ==6)
548  {
549  if (!(vm.count("axis1"))) missingParam("--axis1");
550  if (!(vm.count("axis2"))) missingParam("--axis2");
551  if (!(vm.count("phi"))) missingParam("--phi");
552  double a1 = vm["axis1"].as<double>();
553  double a2 = vm["axis2"].as<double>();
554  double phi = vm["phi"].as<double>();
555 
556  Ellipse2D<Z2i::Space> ell(Z2i::Point(0,0), a1, a2,phi);
557 
558  lengthEstimators<Ellipse2D<Z2i::Space>,Z2i::Space>("Ellipse",ell,h);
559 
560  }
561 
562  h = h * step;
563  }
564  return 0;
565 }
DGtal
Definition: ATu0v1.h:56
Companion project associated to DGtal ; Generated on Mon Feb 13 2017 07:58:04 for DGtalTools by   doxygen 1.8.10