34 #include "DGtal/base/Common.h"
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"
45 #include "DGtal/images/ImageSelector.h"
46 #include "DGtal/geometry/volumes/distance/DistanceTransformation.h"
47 #include "DGtal/geometry/volumes/distance/ExactPredicateLpSeparableMetric.h"
50 #include "DGtal/geometry/volumes/distance/FMM.h"
53 #include "DGtal/io/colormaps/HueShadeColorMap.h"
54 #include "DGtal/io/boards/Board2D.h"
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"
66 using namespace DGtal;
71 template <
typename TImage,
typename TSet,
int norm>
77 template <
typename TImage,
typename TSet>
78 struct DistanceTraits<TImage, TSet, 1>
85 template <
typename TPo
int>
93 BallPredicate(
double aCx,
double aCy,
double aR ):
94 myCx(aCx), myCy(aCy), myR(aR)
97 bool operator()(
const TPoint&
aPoint)
const
99 double d = std::sqrt( std::pow( (myCx-
aPoint[0] ), 2)
100 + std::pow( (myCy-
aPoint[1] ), 2) );
101 if (d <= myR)
return true;
105 double myCx, myCy, myR;
108 template <
typename TPo
int>
112 typedef TPoint
Point;
113 typedef double Value;
116 BallFunctor(
double aCx,
double aCy,
double aR ):
117 myCx(aCx), myCy(aCy), myR(aR)
118 { ASSERT(myR > 0); };
122 double d = std::sqrt( std::pow( (myCx-
aPoint[0] ), 2)
123 + std::pow( (myCy-
aPoint[1] ), 2) );
127 double myCx, myCy, myR;
131 template<
typename TKSpace>
136 ASSERT( aR < (
double) size );
148 std::cerr <<
" error in creating KSpace." << std::endl;
152 BallPredicate<Point> dig(aCx, aCy, aR);
157 std::vector<Point> points;
163 std::cerr <<
" error in finding a bel." << std::endl;
167 template<
typename TIterator >
168 void draw(
const TIterator& itb,
const TIterator& ite,
const int& size, std::string basename)
170 typedef typename std::iterator_traits<TIterator>::value_type Pair;
171 typedef typename Pair::first_type
Point;
178 for ( ; it != ite; ++it)
186 s << basename <<
".eps";
187 b.saveEPS(s.str().c_str());
203 Domain d(Point::diagonal(-size), Point::diagonal(size));
208 Image map( d, (size*size) );
209 map.
setValue( Point::diagonal(0), 0.0 );
217 FMM fmm(map, set, dp, area, distance);
225 s <<
"DTFrom2dPt-" << size <<
"-" << area <<
"-" << distance;
226 draw(map.begin(), map.end(), size, s.str());
241 Domain d(Point::diagonal(-size), Point::diagonal(size));
242 double h = 1.0/(double)size;
245 int radius = (size/2);
246 typedef BallPredicate<Point> Predicate;
247 Predicate predicate( 0, 0, radius );
250 trace.
info() <<
" # circle of radius 0.5 "
251 <<
"digitized in a square of size 1 "
252 <<
"at step h=" << h << endl;
256 KSpace K;
K.
init( Point::diagonal(-size), Point::diagonal(size),
true);
259 std::vector<KSpace::SCell> vSCells;
262 double diff1, diff2, diff3 = 0.0;
276 FMM::initFromBelsRange(
K,
277 vSCells.begin(), vSCells.end(),
283 FMM fmm(map, set, predicate);
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;
310 typedef BallFunctor<Point>
Functor;
311 Functor functor( 0, 0, radius );
313 FMM::initFromBelsRange(
K,
314 vSCells.begin(), vSCells.end(),
319 FMM fmm(map, set, predicate);
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;
348 typedef BallFunctor<Point>
Functor;
349 Functor functor( 0, 0, radius );
351 FMM::initFromBelsRange(
K,
352 vSCells.begin(), vSCells.end(),
356 Distance distance(map, set);
357 FMM fmm(map, set, predicate, distance);
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;
375 return ( (diff1 >= diff2)&&(diff2 >= diff3) );
390 Domain d(Point::diagonal(-size), Point::diagonal(size));
396 map.
setValue( Point::diagonal(0), 0.0 );
404 FMM fmm(map, set, dp, area, distance);
417 for ( ; it != d.
end(); ++it)
430 s <<
"DTFrom3dPt-" << size <<
"-" << area <<
"-" << distance
432 b.saveEPS(s.str().c_str());
446 Domain d(Point::diagonal(-size), Point::diagonal(size));
456 double radius = (rand()%size);
457 trace.
info() <<
" #ball c(" << 0 <<
"," << 0 <<
") r=" << radius << endl;
458 ballGenerator<KSpace>( size, 0, 0, radius, gc );
461 unsigned int nbok = 0;
467 typedef BallPredicate<Point> Predicate;
474 FMM::initFromPointsRange(r.begin(), r.end(), map, set, 0.5);
477 Predicate bp(0,0,radius);
478 FMM fmm(map, set, bp);
482 trace.
info() << nbok <<
"/" << ++nb << std::endl;
489 s <<
"DTInCircle-" << size;
490 draw(map.begin(), map.end(), size, s.str());
507 FMM::initFromPointsRange(r.begin(), r.end(), map, set, 0.5);
510 BallPredicate<Point> bp(0,0,radius);
511 PointPredicate nbp( bp );
512 Predicate pred( nbp, dp, andBF2 );
513 FMM fmm(map, set, pred);
517 trace.
info() << nbok <<
"/" << ++nb << std::endl;
524 s <<
"DTOutCircle-" << size;
525 draw(map.begin(), map.end(), size, s.str());
540 FMM::initFromIncidentPointsRange(r.begin(), r.end(), map, set, 0.5);
543 FMM fmm(map, set, dp);
547 trace.
info() << nbok <<
"/" << ++nb << std::endl;
555 s <<
"DTfromCircle-" << size;
556 draw(map.begin(), map.end(), size, s.str());
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;
578 template<Dimension dim,
int norm>
587 Domain d(Point::diagonal(-size), Point::diagonal(size));
593 map.
setValue( Point::diagonal(0), 0);
596 set.insert( Point::diagonal(0) );
599 typedef Image Image2;
603 for ( ; dit != ditEnd; ++dit)
605 image.setValue(*dit, 128);
607 image.setValue(Point::diagonal(0), 0);
614 Distance distance(map, set);
615 FMM fmm( map, set, dp, area, dist, distance );
623 Predicate aPredicate(image,0);
627 DT dt(&d,&aPredicate, &lnorm);
632 bool flagIsOk =
true;
638 for ( ; ( (it != itEnd)&&(flagIsOk) ); ++it)
640 if (set.find(*it) == set.end())
644 if (dt(*it) != map(*it))
659 int main (
int argc,
char** argv )
663 for (
int i = 0; i < argc; ++i )
669 int area = int( std::pow(
double(2*size+1),2) )+1;
679 area = 4*int( std::pow(
double(size),3) );
686 area = int( std::pow(
double(2*size+1),3) )+1;
688 && testComparison<3,1>( size, area, 3*size+1 )
691 area = int( std::pow(
double(2*size+1),4) ) + 1;
693 && testComparison<4,1>( size, area, 4*size+1 )
697 trace.
emphase() << ( res ?
"Passed." :
"Error." ) << endl;