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