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