DGtal 1.3.0
Loading...
Searching...
No Matches
testFMM.cpp
Go to the documentation of this file.
1
30#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
65using namespace std;
66using namespace DGtal;
67using namespace DGtal::functors;
68
70//
71template <typename TImage, typename TSet, int norm>
72struct DistanceTraits
73{
74 typedef LInfLocalDistance<TImage, TSet> Distance;
75};
76//partial specialization
77template <typename TImage, typename TSet>
78struct DistanceTraits<TImage, TSet, 1>
79{
80 typedef L1LocalDistance<TImage, TSet> Distance;
81};
82
84// digital circle generator
85template <typename TPoint>
86class BallPredicate
87{
88public:
89 typedef TPoint Point;
90
91public:
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 };
104private:
105 double myCx, myCy, myR;
106};
107
108template <typename TPoint>
109class BallFunctor
110{
111public:
112 typedef TPoint Point;
113 typedef double Value;
114public:
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 };
126private:
127 double myCx, myCy, myR;
128};
129
130
131template<typename TKSpace>
132void
133ballGenerator(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
167template< typename TIterator >
168void 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
195bool 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
233bool 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
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
274
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
309
310 typedef BallFunctor<Point> Functor;
311 Functor functor( 0, 0, radius );
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
347
348 typedef BallFunctor<Point> Functor;
349 Functor functor( 0, 0, radius );
350
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
382bool 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
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
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
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;
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;
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;
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
578template<Dimension dim, int norm>
579bool 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
659int 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) )
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// //
Aim: This class specializes a 'Board' class so as to display DGtal objects more naturally (with <<)....
Definition: Board2D.h:71
Aim: model of CConstBidirectionalRange that adapts any range of elements bounded by two iterators [it...
Aim: A container class for storing sets of digital points within some given domain.
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 separable l_p metrics with exact predicates.
Aim: Fast Marching Method (FMM) for nd distance transforms.
Definition: FMM.h:151
static void initFromIncidentPointsRange(const TIteratorOnPairs &itb, const TIteratorOnPairs &ite, Image &aImg, AcceptedPointSet &aSet, const Value &aValue, bool aFlagIsPositive=true)
void compute()
Value min() const
bool isValid() const
static void initFromBelsRange(const KSpace &aK, const TIteratorOnBels &itb, const TIteratorOnBels &ite, Image &aImg, AcceptedPointSet &aSet, const Value &aValue, bool aFlagIsPositive=true)
Value getMin() const
Value getMax() const
Value max() const
static void initFromPointsRange(const TIteratorOnPoints &itb, const TIteratorOnPoints &ite, Image &aImg, AcceptedPointSet &aSet, const Value &aValue)
Aim: describes, in a cellular space of dimension n, a closed or open sequence of signed d-cells (or d...
Definition: GridCurve.h:173
IncidentPointsRange getIncidentPointsRange() const
Definition: GridCurve.h:486
OuterPointsRange getOuterPointsRange() const
Definition: GridCurve.h:474
InnerPointsRange getInnerPointsRange() const
Definition: GridCurve.h:462
bool initFromVector(const std::vector< Point > &aVectorOfPoints)
Aim: This class template may be used to (linearly) convert scalar values in a given range into a colo...
Iterator for HyperRectDomain.
Aim: Parallelepidec region of a digital space, model of a 'CDomain'.
const ConstIterator & begin() const
const ConstIterator & end() const
Aim: implements association bewteen points lying in a digital domain and values.
Definition: Image.h:70
void setValue(const Point &aPoint, const Value &aValue)
Definition: Image.h:247
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
bool init(const Point &lower, const Point &upper, bool isClosed)
Specifies the upper and lower bounds for the maximal cells in this space.
Aim: Class for the computation of the L1-distance at some point p, from the available distance values...
Aim: Class for the computation of the Euclidean distance at some point p, from the available distance...
Aim: Class for the computation of the LInf-distance at some point p, from the available distance valu...
std::string className() const
static void track2DBoundaryPoints(std::vector< Point > &aVectorOfPoints, const KSpace &K, const SurfelAdjacency< KSpace::dimension > &surfel_adj, const PointPredicate &pp, const SCell &start_surfel)
static SCell findABel(const KSpace &K, const PointPredicate &pp, unsigned int nbtries=1000)
static void track2DBoundary(std::vector< SCell > &aSCellContour2D, const KSpace &K, const SurfelAdjacency< KSpace::dimension > &surfel_adj, const PointPredicate &pp, const SCell &start_surfel)
Aim: Represent adjacencies between surfel elements, telling if it follows an interior to exterior ord...
void beginBlock(const std::string &keyword="")
std::ostream & emphase()
std::ostream & info()
double endBlock()
Aim: Define a simple Foreground predicate thresholding image values given a single thresold....
void setUnit(Unit unit)
Definition: Board.cpp:240
Z3i::SCell SCell
functors namespace gathers all DGtal functors.
static const BoolFunctor2 andBF2
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::uint32_t Dimension
Definition: Common.h:137
Trace trace
Definition: Common.h:154
STL namespace.
Custom style class redefining the fill color. You may use Board2D::Color::None for transparent color.
Definition: Board2D.h:343
Represents a signed cell in a cellular grid space by its Khalimsky coordinates and a boolean value.
Aim: The predicate returns true when the given binary functor returns true for the two PointPredicate...
Aim: The predicate returning true iff the point is in the domain given at construction....
Aim: The predicate returns true when the point predicate given at construction return false....
int main()
Definition: testBits.cpp:56
MyPointD Point
Definition: testClone2.cpp:383
InHalfPlaneBySimple3x3Matrix< Point, double > Functor
KSpace K
bool testComparison(int size, int area, double dist)
Definition: testFMM.cpp:579
bool accuracyTest(int size)
Definition: testFMM.cpp:233
bool testDisplayDT3d(int size, int area, double distance)
Definition: testFMM.cpp:382
bool testDisplayDT2d(int size, int area, double distance)
Definition: testFMM.cpp:195
void draw(const TIterator &itb, const TIterator &ite, const int &size, std::string basename)
Definition: testFMM.cpp:168
bool testDisplayDTFromCircle(int size)
Definition: testFMM.cpp:438
void ballGenerator(const int &size, double aCx, double aCy, double aR, GridCurve< TKSpace > &gc)
Definition: testFMM.cpp:133
const Point aPoint(3, 4)
HyperRectDomain< Space > Domain