DGtal  0.9.2
testFMM.cpp
1 
29 #include <iostream>
31 #include <iomanip>
32 #include <functional>
33 
34 #include "DGtal/base/Common.h"
35 
36 #include "DGtal/kernel/SpaceND.h"
37 #include "DGtal/kernel/domains/HyperRectDomain.h"
38 #include "DGtal/kernel/BasicPointPredicates.h"
39 #include "DGtal/kernel/domains/DomainPredicate.h"
40 #include "DGtal/kernel/sets/DigitalSetFromMap.h"
41 #include "DGtal/images/ImageContainerBySTLMap.h"
42 #include "DGtal/images/SimpleThresholdForegroundPredicate.h"
43 
44 //DT
45 #include "DGtal/images/ImageSelector.h"
46 #include "DGtal/geometry/volumes/distance/DistanceTransformation.h"
47 #include "DGtal/geometry/volumes/distance/ExactPredicateLpSeparableMetric.h"
48 
49 //FMM
50 #include "DGtal/geometry/volumes/distance/FMM.h"
51 
52 //Display
53 #include "DGtal/io/colormaps/HueShadeColorMap.h"
54 #include "DGtal/io/boards/Board2D.h"
55 
56 //shape and digitizer
57 #include "DGtal/shapes/ShapeFactory.h"
58 #include "DGtal/shapes/Shapes.h"
59 #include "DGtal/topology/helpers/Surfaces.h"
60 #include "DGtal/shapes/GaussDigitizer.h"
61 #include "DGtal/geometry/curves/GridCurve.h"
62 
64 
65 using namespace std;
66 using namespace DGtal;
67 using namespace DGtal::functors;
68 
70 //
71 template <typename TImage, typename TSet, int norm>
72 struct DistanceTraits
73 {
74  typedef LInfLocalDistance<TImage, TSet> Distance;
75 };
76 //partial specialization
77 template <typename TImage, typename TSet>
78 struct DistanceTraits<TImage, TSet, 1>
79 {
80  typedef L1LocalDistance<TImage, TSet> Distance;
81 };
82 
84 // digital circle generator
85 template <typename TPoint>
86 class BallPredicate
87 {
88 public:
89  typedef TPoint Point;
90 
91 public:
92 
93  BallPredicate(double aCx, double aCy, double aR ):
94  myCx(aCx), myCy(aCy), myR(aR)
95  { ASSERT(myR > 0); };
96 
97  bool operator()(const TPoint& aPoint) const
98  {
99  double d = std::sqrt( std::pow( (myCx-aPoint[0] ), 2)
100  + std::pow( (myCy-aPoint[1] ), 2) );
101  if (d <= myR) return true;
102  else return false;
103  };
104 private:
105  double myCx, myCy, myR;
106 };
107 
108 template <typename TPoint>
109 class BallFunctor
110 {
111 public:
112  typedef TPoint Point;
113  typedef double Value;
114 public:
115 
116  BallFunctor(double aCx, double aCy, double aR ):
117  myCx(aCx), myCy(aCy), myR(aR)
118  { ASSERT(myR > 0); };
119 
120  Value operator()(const TPoint& aPoint) const
121  {
122  double d = std::sqrt( std::pow( (myCx-aPoint[0] ), 2)
123  + std::pow( (myCy-aPoint[1] ), 2) );
124  return (d - myR);
125  };
126 private:
127  double myCx, myCy, myR;
128 };
129 
130 
131 template<typename TKSpace>
132 void
133 ballGenerator(const int& size, double aCx, double aCy, double aR, GridCurve<TKSpace>& gc)
134 {
135 
136  ASSERT( aR < (double) size );
137 
138  // Types
139  typedef TKSpace KSpace;
140  typedef typename KSpace::SCell SCell;
141  typedef typename KSpace::Space Space;
142  typedef typename Space::Point Point;
143 
144  KSpace K;
145  bool ok = K.init( Point(-size,-size), Point(size,size), true );
146  if ( ! ok )
147  {
148  std::cerr << " error in creating KSpace." << std::endl;
149  }
150  try
151  {
152  BallPredicate<Point> dig(aCx, aCy, aR);
153  // Extracts shape boundary
155  SCell bel = Surfaces<KSpace>::findABel( K, dig, 10000 );
156  // Getting the consecutive surfels of the 2D boundary
157  std::vector<Point> points;
158  Surfaces<KSpace>::track2DBoundaryPoints( points, K, SAdj, dig, bel );
159  gc.initFromVector(points);
160  }
161  catch ( InputException e )
162  {
163  std::cerr << " error in finding a bel." << std::endl;
164  }
165 }
166 
167 template< typename TIterator >
168 void draw( const TIterator& itb, const TIterator& ite, const int& size, std::string basename)
169 {
170  typedef typename std::iterator_traits<TIterator>::value_type Pair;
171  typedef typename Pair::first_type Point;
172  HueShadeColorMap<unsigned char, 2> colorMap(0,3*size);
173 
174  Board2D b;
176 
177  TIterator it = itb;
178  for ( ; it != ite; ++it)
179  {
180  Point p = it->first;
181  b << CustomStyle( p.className(), new CustomFillColor( colorMap( it->second) ) );
182  b << p;
183  }
184 
185  std::stringstream s;
186  s << basename << ".eps";
187  b.saveEPS(s.str().c_str());
188 }
189 
190 
195 bool testDisplayDT2d(int size, int area, double distance)
196 {
197 
198  static const DGtal::Dimension dimension = 2;
199 
200  //Domain
202  typedef Domain::Point Point;
203  Domain d(Point::diagonal(-size), Point::diagonal(size));
205 
206  //Image and set
208  Image map( d, (size*size) );
209  map.setValue( Point::diagonal(0), 0.0 );
210  typedef DigitalSetFromMap<Image> Set;
211  Set set(map);
212 
213  //computation
214  trace.beginBlock ( "Display 2d FMM results " );
215 
217  FMM fmm(map, set, dp, area, distance);
218  fmm.compute();
219  trace.info() << fmm << std::endl;
220 
221  trace.endBlock();
222 
223  //display
224  std::stringstream s;
225  s << "DTFrom2dPt-" << size << "-" << area << "-" << distance;
226  draw(map.begin(), map.end(), size, s.str());
227 
228  return fmm.isValid();
229 }
230 
231 
233 bool accuracyTest(int size)
234 {
235 
236  static const DGtal::Dimension dimension = 2;
237 
238  //Domain
240  typedef Domain::Point Point;
241  Domain d(Point::diagonal(-size), Point::diagonal(size));
242  double h = 1.0/(double)size;
243 
244  //predicate
245  int radius = (size/2);
246  typedef BallPredicate<Point> Predicate;
247  Predicate predicate( 0, 0, radius );
248 
249  trace.beginBlock ( "test of accuracy" );
250  trace.info() << " # circle of radius 0.5 "
251  << "digitized in a square of size 1 "
252  << "at step h=" << h << endl;
253 
254  //Digital circle generation
255  typedef KhalimskySpaceND< dimension, int > KSpace;
256  KSpace K; K.init( Point::diagonal(-size), Point::diagonal(size), true);
258  KSpace::SCell bel = Surfaces<KSpace>::findABel( K, predicate, 10000 );
259  std::vector<KSpace::SCell> vSCells;
260  Surfaces<KSpace>::track2DBoundary( vSCells, K, SAdj, predicate, bel );
261 
262  double diff1, diff2, diff3 = 0.0;
263  {
264  //Image and set
266  typedef DigitalSetFromMap<Image> Set;
267  Image map( d );
268  Set set(map);
269 
270  //initialisation
272  typedef FMM<Image, Set, Predicate > FMM;
274 
276  FMM::initFromBelsRange( K,
277  vSCells.begin(), vSCells.end(),
278  map, set, 0.5 );
280 
281  //computation
283  FMM fmm(map, set, predicate);
284  fmm.compute();
285  trace.info() << fmm << std::endl;
287 
288  //max
289  double truth = radius*h;
290  double found = ( std::max(std::abs(fmm.max()),std::abs(fmm.min())) )*h;
291  double diff = std::abs(found-truth);
292  trace.info() << " # radius (low accuracy)" << std::endl;
293  trace.info() << " # truth: " << truth << std::endl;
294  trace.info() << " # found: " << found << std::endl;
295  trace.info() << " # diff.: " << diff << std::endl;
296 
297  diff1 = diff;
298  }
299 
300  {
301  //Image and set
303  typedef DigitalSetFromMap<Image> Set;
304  Image map( d );
305  Set set(map);
306 
307  //initialisation
308  typedef FMM<Image, Set, Predicate > FMM;
309 
310  typedef BallFunctor<Point> Functor;
311  Functor functor( 0, 0, radius );
313  FMM::initFromBelsRange( K,
314  vSCells.begin(), vSCells.end(),
315  functor, map, set );
317 
318  //computation
319  FMM fmm(map, set, predicate);
320  fmm.compute();
321  trace.info() << fmm << std::endl;
322 
323  //max
324  double truth = radius*h;
325  double found = ( std::max(std::abs(fmm.max()),std::abs(fmm.min())) )*h;
326  double diff = std::abs(found-truth);
327  trace.info() << " # radius (medium accuracy)" << std::endl;
328  trace.info() << " # truth: " << truth << std::endl;
329  trace.info() << " # found: " << found << std::endl;
330  trace.info() << " # diff.: " << diff << std::endl;
331 
332  diff2 = diff;
333  }
334 
335  {
336  //Image and set
338  typedef DigitalSetFromMap<Image> Set;
339  Image map( d );
340  Set set(map);
341 
342  //initialisation
344  typedef L2SecondOrderLocalDistance<Image, Set> Distance;
347 
348  typedef BallFunctor<Point> Functor;
349  Functor functor( 0, 0, radius );
350 
351  FMM::initFromBelsRange( K,
352  vSCells.begin(), vSCells.end(),
353  functor, map, set );
354 
355  //computation
356  Distance distance(map, set);
357  FMM fmm(map, set, predicate, distance);
358  fmm.compute();
359  trace.info() << fmm << std::endl;
360 
361  //max
362  double truth = radius*h;
363  double found = ( std::max(std::abs(fmm.max()),std::abs(fmm.min())) )*h;
364  double diff = std::abs(found-truth);
365  trace.info() << " # radius (high accuracy)" << std::endl;
366  trace.info() << " # truth: " << truth << std::endl;
367  trace.info() << " # found: " << found << std::endl;
368  trace.info() << " # diff.: " << diff << std::endl;
369 
370  diff3 = diff;
371  }
372 
373  trace.endBlock();
374 
375  return ( (diff1 >= diff2)&&(diff2 >= diff3) );
376 }
377 
382 bool testDisplayDT3d(int size, int area, double distance)
383 {
384 
385  static const DGtal::Dimension dimension = 3;
386 
387  //Domain
389  typedef Domain::Point Point;
390  Domain d(Point::diagonal(-size), Point::diagonal(size));
392 
393  //Image and set
395  Image map( d, 0.0 );
396  map.setValue( Point::diagonal(0), 0.0 );
397  typedef DigitalSetFromMap<Image> Set;
398  Set set(map);
399 
400  //computation
401  trace.beginBlock ( "Display 3d FMM results " );
402 
404  FMM fmm(map, set, dp, area, distance);
405  fmm.compute();
406  trace.info() << fmm << std::endl;
407 
408  trace.endBlock();
409 
410  { //display
411  HueShadeColorMap<unsigned char, 2> colorMap(0,2*size);
412 
413  Board2D b;
415 
416  Domain::ConstIterator it = d.begin();
417  for ( ; it != d.end(); ++it)
418  {
419  Point p3 = *it;
420  if (p3[2] == 0)
421  {
422  PointVector<2,Point::Coordinate> p2(p3[0], p3[1]);
423  b << CustomStyle( p2.className(),
424  new CustomFillColor( colorMap(map(p3)) ) )
425  << p2;
426  }
427  }
428 
429  std::stringstream s;
430  s << "DTFrom3dPt-" << size << "-" << area << "-" << distance
431  << ".eps";
432  b.saveEPS(s.str().c_str());
433  }
434 
435  return fmm.isValid();
436 }
437 
438 bool testDisplayDTFromCircle(int size)
439 {
440 
441  static const DGtal::Dimension dimension = 2;
442 
443  //Domain
445  typedef Domain::Point Point;
446  Domain d(Point::diagonal(-size), Point::diagonal(size));
448 
449  //Image and set
451  typedef DigitalSetFromMap<Image> Set;
452 
453  //Digital circle generation
454  typedef KhalimskySpaceND< dimension, int > KSpace;
455  GridCurve<KSpace> gc;
456  double radius = (rand()%size);
457  trace.info() << " #ball c(" << 0 << "," << 0 << ") r=" << radius << endl;
458  ballGenerator<KSpace>( size, 0, 0, radius, gc );
459 
460 
461  unsigned int nbok = 0;
462  unsigned int nb = 0;
463 
464  double dmaxInt = 0;
465  trace.beginBlock ( "Interior " );
466  {
467  typedef BallPredicate<Point> Predicate;
468  typedef FMM<Image, Set, Predicate > FMM;
469 
470  //init
471  Image map( d );
472  Set set(map);
474  FMM::initFromPointsRange(r.begin(), r.end(), map, set, 0.5);
475 
476  //computation
477  Predicate bp(0,0,radius);
478  FMM fmm(map, set, bp);
479  fmm.compute();
480  trace.info() << fmm << std::endl;
481  nbok += (fmm.isValid()?1:0);
482  trace.info() << nbok << "/" << ++nb << std::endl;
483 
484  //max
485  dmaxInt = fmm.getMax();
486 
487  //display
488  std::stringstream s;
489  s << "DTInCircle-" << size;
490  draw(map.begin(), map.end(), size, s.str());
491 
492  }
493  trace.endBlock();
494 
495  double dmaxExt = 0;
496  trace.beginBlock ( "Exterior " );
497  {
498  typedef NotPointPredicate<BallPredicate<Point> > PointPredicate;
499  typedef BinaryPointPredicate<PointPredicate,
500  DomainPredicate<Domain> > Predicate;
501  typedef FMM<Image, Set, Predicate > FMM;
502 
503  //init
504  Image map( d );
505  Set set(map);
507  FMM::initFromPointsRange(r.begin(), r.end(), map, set, 0.5);
508 
509  //computation
510  BallPredicate<Point> bp(0,0,radius);
511  PointPredicate nbp( bp );
512  Predicate pred( nbp, dp, andBF2 );
513  FMM fmm(map, set, pred);
514  fmm.compute();
515  trace.info() << fmm << std::endl;
516  nbok += (fmm.isValid()?1:0);
517  trace.info() << nbok << "/" << ++nb << std::endl;
518 
519  //max
520  dmaxExt = fmm.getMax();
521 
522  //display
523  std::stringstream s;
524  s << "DTOutCircle-" << size;
525  draw(map.begin(), map.end(), size, s.str());
526  }
527  trace.endBlock();
528 
529  double dmin = 0; //2*size*size;
530  double dmax = 0;
531  trace.beginBlock ( "Both " );
532  {
533  typedef DomainPredicate<Domain> Predicate;
534  typedef FMM<Image, Set, Predicate > FMM;
535 
536  //init
537  Image map( d );
538  Set set(map);
540  FMM::initFromIncidentPointsRange(r.begin(), r.end(), map, set, 0.5);
541 
542  //computation
543  FMM fmm(map, set, dp);
544  fmm.compute();
545  trace.info() << fmm << std::endl;
546  nbok += (fmm.isValid()?1:0);
547  trace.info() << nbok << "/" << ++nb << std::endl;
548 
549  //min, max
550  dmin = fmm.getMin();
551  dmax = fmm.getMax();
552 
553  //display
554  std::stringstream s;
555  s << "DTfromCircle-" << size;
556  draw(map.begin(), map.end(), size, s.str());
557  }
558  trace.endBlock();
559 
560  trace.beginBlock ( "Comparison " );
561  {
562  double epsilon = 0.0001;
563  nbok += ( ( (std::abs(-dmaxInt - dmin) < epsilon)
564  && (std::abs(dmaxExt - dmax) < epsilon) )?1:0);
565  trace.info() << nbok << "/" << ++nb << std::endl;
566  }
567  trace.endBlock();
568 
569  return (nb == nbok);
570 
571 }
572 
573 
578 template<Dimension dim, int norm>
579 bool testComparison(int size, int area, double dist)
580 {
581 
582  static const DGtal::Dimension dimension = dim;
583 
584  //Domain
586  typedef typename Domain::Point Point;
587  Domain d(Point::diagonal(-size), Point::diagonal(size));
589 
590  //Image and set for FMM
592  Image map( d );
593  map.setValue( Point::diagonal(0), 0);
594  typedef DigitalSetBySTLSet<Domain> Set;
595  Set set( d );
596  set.insert( Point::diagonal(0) );
597 
598  //Image for separable DT
599  typedef Image Image2;
600  Image2 image ( d );
601  typename Domain::Iterator dit = d.begin();
602  typename Domain::Iterator ditEnd = d.end();
603  for ( ; dit != ditEnd; ++dit)
604  {
605  image.setValue(*dit, 128);
606  }
607  image.setValue(Point::diagonal(0), 0);
608 
609  //computation
610  trace.beginBlock ( " FMM computation " );
611 
612  typedef typename DistanceTraits<Image,Set,norm>::Distance Distance;
613  typedef FMM<Image, Set, DomainPredicate<Domain>, Distance > FMM;
614  Distance distance(map, set);
615  FMM fmm( map, set, dp, area, dist, distance );
616  fmm.compute();
617  trace.info() << fmm << std::endl;
618 
619  trace.endBlock();
620 
621  trace.beginBlock ( " DT computation " );
623  Predicate aPredicate(image,0);
625  typedef DistanceTransformation<SpaceND<dimension,int>, Predicate, LNorm> DT;
626  LNorm lnorm;
627  DT dt(&d,&aPredicate, &lnorm);
628  trace.info() << dt << std::endl;
629 
630  trace.endBlock();
631 
632  bool flagIsOk = true;
633 
634  trace.beginBlock ( " Comparison " );
635  //all points of result must be in map and have the same distance
636  typename Domain::ConstIterator it = d.begin();
637  typename Domain::ConstIterator itEnd = d.end();
638  for ( ; ( (it != itEnd)&&(flagIsOk) ); ++it)
639  {
640  if (set.find(*it) == set.end())
641  flagIsOk = false;
642  else
643  {
644  if (dt(*it) != map(*it))
645  flagIsOk = false;
646  }
647  }
648  trace.endBlock();
649 
650  return flagIsOk;
651 
652 }
653 
654 
655 
657 // Standard services - public :
658 
659 int main ( int argc, char** argv )
660 {
661  trace.beginBlock ( "Testing FMM" );
662  trace.info() << "Args:";
663  for ( int i = 0; i < argc; ++i )
664  trace.info() << " " << argv[ i ];
665  trace.info() << endl;
666 
667  //2d L2 tests
668  int size = 50;
669  int area = int( std::pow(double(2*size+1),2) )+1;
670  bool res
671  = testDisplayDT2d( size, area, std::sqrt(2*size*size) )
672  && testDisplayDT2d( size, area, size )
673  && testDisplayDT2d( size, 2*area, std::sqrt(2*size*size) )
674  && testDisplayDTFromCircle(size)
675  && accuracyTest(size)
676  ;
677 
678  size = 25;
679  area = 4*int( std::pow(double(size),3) );
680  //3d L2 test
681  res = res && testDisplayDT3d( size, area, std::sqrt(size*size*size) )
682  ;
683 
684  //3d L1 and comparison
685  size = 20;
686  area = int( std::pow(double(2*size+1),3) )+1;
687  res = res
688  && testComparison<3,1>( size, area, 3*size+1 )
689  ;
690  size = 5;
691  area = int( std::pow(double(2*size+1),4) ) + 1;
692  res = res
693  && testComparison<4,1>( size, area, 4*size+1 )
694  ;
695 
696  //&& ... other tests
697  trace.emphase() << ( res ? "Passed." : "Error." ) << endl;
698  trace.endBlock();
699  return res ? 0 : 1;
700 }
701 // //
void beginBlock(const std::string &keyword="")
Aim: The predicate returns true when the point predicate given at construction return false...
InnerPointsRange getInnerPointsRange() const
Definition: GridCurve.h:462
Aim: An adapter for viewing an associative image container like ImageContainerBySTLMap as a simple di...
Aim: Implementation of the linear in time distance transformation for separable metrics.
Aim: implements association bewteen points lying in a digital domain and values.
Definition: Image.h:69
Trace trace
Definition: Common.h:130
Aim: SpaceND is a utility class that defines the fundamental structure of a Digital Space in ND...
Definition: SpaceND.h:95
Custom style class redefining the fill color. You may use Board2D::Color::None for transparent color...
Definition: Board2D.h:342
Aim: A utility class for constructing surfaces (i.e. set of (n-1)-cells).
Definition: Surfaces.h:78
DGtal::uint32_t Dimension
Definition: Common.h:113
KhalimskySpaceND< 2, Integer > KSpace
Definition: StdDefs.h:77
Aim: Parallelepidec region of a digital space, model of a 'CDomain'.
STL namespace.
double endBlock()
Aim: This class template may be used to (linearly) convert scalar values in a given range into a colo...
functors namespace gathers all DGtal functors.
Aim: implements separable l_p metrics with exact predicates.
Aim: Fast Marching Method (FMM) for nd distance transforms.
Definition: FMM.h:150
bool initFromVector(const std::vector< Point > &aVectorOfPoints)
Aim: Class for the computation of the Euclidean distance at some point p, from the available distance...
IncidentPointsRange getIncidentPointsRange() const
Definition: GridCurve.h:486
Aim: Define a simple Foreground predicate thresholding image values given a single thresold...
std::ostream & emphase()
bool init(const Point &lower, const Point &upper, bool isClosed)
ConstIterator begin() const
DGtal is the top-level namespace which contains all DGtal functions and types.
Aim: model of CConstBidirectionalRange that adapts any range of elements bounded by two iterators [it...
Aim: Class for the computation of the L1-distance at some point p, from the available distance values...
std::ostream & info()
ConstIterator end() const
void saveEPS(const char *filename, PageSize size=Board::BoundingBox, double margin=10.0) const
Definition: Board.cpp:805
void setUnit(Unit unit)
Definition: Board.cpp:240
Aim: The predicate returns true when the given binary functor returns true for the two PointPredicate...
Aim: describes, in a cellular space of dimension n, a closed or open sequence of signed d-cells (or d...
Definition: GridCurve.h:172
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex...
OuterPointsRange getOuterPointsRange() const
Definition: GridCurve.h:474
Aim: A container class for storing sets of digital points within some given domain.
Aim: This class specializes a 'Board' class so as to display DGtal objects more naturally (with <<)...
Definition: Board2D.h:70
Aim: Class for the computation of the LInf-distance at some point p, from the available distance valu...