DGtal  1.5.beta
testDistanceTransformation.cpp
Go to the documentation of this file.
1 
31 #include <iostream>
32 #include <iomanip>
33 #include "DGtal/base/Common.h"
34 
35 #include "DGtal/kernel/SpaceND.h"
36 #include "DGtal/kernel/domains/HyperRectDomain.h"
37 #include "DGtal/images/ImageSelector.h"
38 #include "DGtal/geometry/volumes/distance/DistanceTransformation.h"
39 #include "DGtal/geometry/volumes/distance/VoronoiMap.h"
40 #include "DGtal/geometry/volumes/distance/ExactPredicateLpSeparableMetric.h"
41 #include "DGtal/geometry/volumes/distance/InexactPredicateLpSeparableMetric.h"
42 #include "DGtal/io/colormaps/HueShadeColorMap.h"
43 #include "DGtal/io/colormaps/GrayscaleColorMap.h"
44 #include "DGtal/shapes/Shapes.h"
45 #include "DGtal/helpers/StdDefs.h"
46 #include "DGtal/shapes/ShapeFactory.h"
47 #include "DGtal/io/boards/Board2D.h"
48 #include "DGtal/images/SimpleThresholdForegroundPredicate.h"
50 
51 using namespace std;
52 using namespace DGtal;
53 using namespace DGtal::functors;
54 
56 // Functions for testing class DistanceTransformation.
58 
59 template<typename Image>
60 void randomSeeds(Image &input, const unsigned int nb, const int value)
61 {
62  typename Image::Point p, low = input.domain().lowerBound(), up = input.domain().upperBound();
63  typename Image::Vector ext;
64 
65  for (Dimension i = 0; i < Image::Domain::dimension; i++)
66  ext[i] = up[i] - low[i] + 1;
67 
68 
69  for (unsigned int k = 0 ; k < nb; k++)
70  {
71  for (unsigned int dim = 0; dim < Image::dimension; dim++)
72  {
73  p[dim] = rand() % (ext[dim]) + low[dim];
74  }
75  input.setValue(p, value);
76  }
77 }
78 
79 
80 
86 {
87  bool allfine;
88 
89  trace.beginBlock ( "Testing the whole DT computation" );
90 
91  typedef SpaceND<2> TSpace;
92  typedef TSpace::Point Point;
95  Point a ( 2, 2 );
96  Point b ( 15, 15 );
98  Image image ( Domain(a, b ));
99 
100  for ( unsigned k = 0; k < 49; k++ )
101  {
102  a[0] = ( k / 7 ) + 5;
103  a[1] = ( k % 7 ) + 5;
104  image.setValue ( a, 128 );
105  }
106  a= Point(2,2);
107 
109  Predicate aPredicate(image,0);
110 
112  L2Metric l2;
113  Domain dom(a,b);
115  VoronoiMap<Z2i::Space, Predicate, L2Metric> voro(&dom,&aPredicate,&l2);
116 
117  Board2D board;
119  Display2DFactory::drawImage<Gray>(board, image, (unsigned int)0, (unsigned int)255);
120  board.saveSVG ( "image-preDT.svg" );
121  //We just iterate on the Domain points and print out the point coordinates.
122  std::copy ( image.begin(),
123  image.end(),
124  std::ostream_iterator<unsigned int> ( std::cout, " " ) );
125 
126  trace.info()<<std::endl;
127  for(int i=2;i<=15;++i)
128  {
129  for(int j=2;j<=15;++j)
130  trace.info()<< image(Point(i,j))<<" ";
131  trace.info()<<std::endl;
132  }
133 
134 
135  trace.warning() << dt << endl;
136  trace.info() <<std::endl;
137 
140 
141  for (; it != itend; ++it)
142  {
143  std::cout << (*it) << " ";
144  }
145  std::cout << std::endl;
146 
147  trace.info()<<std::endl;
148  for(int i=2;i<=15;++i)
149  {
150  for(int j=2;j<=15;++j)
151  trace.info()<< dt(Point(i,j))<<" ";
152  trace.info()<<std::endl;
153  }
154 
155 
156  trace.info()<<std::endl;
157  for(int i=2;i<=15;++i)
158  {
159  for(int j=2;j<=15;++j)
160  {
161  Point p= dt.getVoronoiSite(Point(i,j));
162  if (p==Point(i,j))
163  trace.info()<<"-,- ";
164  else
165  trace.info()<< p[0]<<","<<p[1]<<" ";
166  }
167  trace.info()<<std::endl;
168  }
169 
170 
171  allfine = true;
172 
173  trace.info()<<std::endl;
174  for(int i=2;i<=15;++i)
175  {
176  for(int j=2;j<=15;++j)
177  {
178  Point p= voro(Point(i,j));
179  if (p != dt.getVoronoiSite(Point(i,j)))
180  allfine = false;
181  if (p==Point(i,j))
182  trace.info()<<"-,- ";
183  else
184  trace.info()<< p[0]<<","<<p[1]<<" ";
185  }
186  trace.info()<<std::endl;
187  }
188 
189  board.clear();
190  Display2DFactory::drawImage<Gray>(board, dt, (DGtal::int64_t)0, (DGtal::int64_t)16);
191  board.saveSVG ( "image-postDT.svg" );
192  trace.info() << dt << endl;
193  trace.endBlock();
194 
195  return allfine;
196 }
202 {
203  unsigned int nbok = 0;
204  unsigned int nb = 0;
205 
206  trace.beginBlock ( "Testing the Neg DT computation" );
207 
208  typedef Z2i::Space TSpace;
209  typedef TSpace::Point Point;
210  typedef Z2i::Domain Domain;
212  Point a ( -10, -10 );
213  Point b ( 10, 10 );
215  Image image ( Domain( a, b ));
216 
217  for(int y=-10; y<=10;y++)
218  for(int x=-10; x<=10;x++)
219  {
220  if ((abs(x)<7) && (abs(y)<5))
221  image.setValue(Point(x,y),1);
222  else
223  image.setValue(Point(x,y),0);
224  }
225 
227  Predicate aPredicate(image,0);
228 
230  L2Metric l2;
231  Board2D board;
233  Display2DFactory::drawImage<Gray>(board, image, (unsigned int)0, (unsigned int)1);
234  board.saveSVG ( "image-preDT-neg.svg" );
235 
236 
237  for(int y=-10; y<=10;y++)
238  {
239  for(int x=-10; x<=10;x++)
240  {
241  std::cout<<image(Point(x,y))<<" ";
242  }
243  std::cout<<std::endl;
244  }
245 
246  trace.info()<<"Domain "<<Domain(a,b)<<std::endl;
247  Domain dom(a,b);
249 
252  itend = dt.constRange().end();
253  it != itend ; ++it)
254  if ((*it) > maxv)
255  maxv = (*it);
256 
257  for(int y=-10; y<=10;y++)
258  {
259  for(int x=-10; x<=10;x++)
260  {
261  std::cout<<dt(Point(x,y))<<" ";
262  }
263  std::cout<<std::endl;
264  }
265 
266 
267 
268  trace.warning() << dt << endl;
269  trace.warning() << dt.domain() << endl;
270 
271 
272  trace.info() << "Exporting..." << endl;
273  board.clear();
274  Display2DFactory::drawImage<Gray>(board, dt, 0, maxv);
275  board.saveSVG ( "image-postDT-neg.svg" );
276 
277  trace.info() << "Done..." << endl;
278 
279  trace.info() << dt << endl;
280 
281  trace.endBlock();
282 
283  return nbok == nb;
284 }
285 
286 
288 {
289  unsigned int nbok = 0;
290  unsigned int nb = 0;
291 
292  trace.beginBlock ( "Testing the whole DT computation from a Set" );
293 
294  typedef SpaceND<2> TSpace;
296 
297 
298  Board2D board;
299 
300  AccFlower2D<Z2i::Space> flower(Z2i::Point(0,0), 30, 5,2,0);
301  Z2i::Domain domain(flower.getLowerBound(), flower.getUpperBound());
302  Z2i::DigitalSet aSet(domain);
303 
305 
306  // Since 0.6, models of CDigitalSet are models of CPointPredicate.
307  // SetPredicate<Z2i::DigitalSet> aPredicate(aSet);
310  L2Metric l2;
311  L2DT dt(&domain,&aSet, &l2);
314  L1Metric l1;
315  L1DT dt1(&domain,&aSet, &l1);
316 
317  L2DT::Value maxv = 0;
318  for ( L2DT::ConstRange::ConstIterator it = dt.constRange().begin(), itend = dt.constRange().end();
319  it != itend; ++it)
320  if ( (*it) > maxv)
321  maxv = (*it);
322  trace.error() << "MaxV="<<maxv<<std::endl;
323  Display2DFactory::drawImage<Hue>(board, dt, 0, maxv+1);
324  board.saveSVG ( "image-DTSet.svg" );
325 
326  board.clear();
327  maxv = 0;
328  for ( L1DT::ConstRange::ConstIterator it = dt1.constRange().begin(), itend = dt1.constRange().end();
329  it != itend; ++it)
330  if ( (*it) > maxv)
331  maxv = (*it);
332  trace.error() << "MaxV="<<maxv<<std::endl;
333  Display2DFactory::drawImage<Hue>(board, dt1, 0, maxv+1);
334  board.saveSVG ( "image-DTSet-l1.svg" );
335  trace.endBlock();
336 
337  return nbok == nb;
338 }
339 
345 {
346  unsigned int nbok = 0;
347  unsigned int nb = 0;
348 
349  trace.beginBlock ( "Testing DT computation with Infinity values at the first step" );
350 
351  typedef SpaceND<2> TSpace;
352  typedef TSpace::Point Point;
355 
356  Point a (0, 0 );
357  Point b ( 128, 128 );
358 
360  Image image ( Domain(a, b ));
361 
362  for ( Image::Iterator it = image.begin(), itend = image.end();it != itend; ++it)
363  *it = 128 ;
364 
365 
366  randomSeeds(image, 19, 0);
367 
369  Predicate aPredicate(image,0);
370 
372  L2Metric l2;
373  Domain dom(a,b);
375 
376  Board2D board;
378  Display2DFactory::drawImage<Hue>(board, image, (unsigned int)0, (unsigned int)150);
379  board.saveSVG ( "image-preDT-border.svg" );
380 
382  for ( DistanceTransformation<TSpace, Predicate, L2Metric>::ConstRange::ConstIterator it = dt.constRange().begin(), itend = dt.constRange().end();it != itend; ++it)
383  if ( (*it) > maxv)
384  maxv = (*it);
385 
387  for (unsigned int y = 0; y < 33; y++)
388  {
389  for (unsigned int x = 0; x < 33; x++)
390  {
391  std::cout << std::setw(4) << (*it) << " ";
392  ++it;
393  }
394  std::cout << std::endl;
395  }
396 
397  trace.warning() << dt << "MaxV = " << maxv << endl;
398 
399 
400  board.clear();
401  Display2DFactory::drawImage<Hue>(board, dt, (DGtal::int64_t)0, (DGtal::int64_t)maxv+1);
402  board.saveSVG ( "image-postDT-border.svg" );
403 
404 
405  trace.info() << dt << endl;
406 
407  trace.endBlock();
408 
409  return nbok == nb;
410 }
411 
412 
418 {
419  unsigned int nbok = 0;
420  unsigned int nb = 0;
421 
422  trace.beginBlock ( "Testing 3D DT computation" );
423 
424  typedef SpaceND<3> TSpace;
425  typedef TSpace::Point Point;
427  Point a ( 0, 0, 0 );
428  Point b ( 15, 15, 15 );
430  Image image ( Domain(a, b ));
431  Point c(8, 8, 8);
432  Domain dom(a, b);
433 
434  for (Domain::ConstIterator it = dom.begin(),
435  itend = dom.end(); it != itend; ++it)
436  {
437  if ( ((*it) - c).norm() < 7)
438  image.setValue ( *it, 128 );
439  }
440 
442  Predicate aPredicate(image,0);
443 
445  L2Metric l2;
447 
448  //We display the values on a 2D slice
449  for (unsigned int y = 0; y < 16; y++)
450  {
451  for (unsigned int x = 0; x < 16; x++)
452  {
453  Point p(x, y, 8);
454  std::cout << dt(p) << " ";
455  }
456  std::cout << std::endl;
457  }
458 
459 
460  trace.warning() << dt << endl;
461 
462  trace.endBlock();
463 
464  return nbok == nb;
465 }
466 
467 
469 {
470  unsigned int nbok = 0;
471  unsigned int nb = 0;
472 
473  trace.beginBlock ( "Testing DT computation with Infinity values at the first step" );
474 
475  typedef SpaceND<2> TSpace;
476  typedef TSpace::Point Point;
479 
480  Point a (0, 0 );
481  Point b ( 128, 128 );
482  Domain dom(a,b);
483 
485  Image image ( dom );
486 
487  for ( Image::Iterator it = image.begin(), itend = image.end();it != itend; ++it)
488  (*it) = 128;
489 
490 
491  randomSeeds(image, 19, 0);
492 
494  Predicate aPredicate(image,0);
495 
496 
497  // L_euc metric
499  L2Metric l2;
501  DT2 dt2(&dom, &aPredicate, &l2);
502  DT2 dt2_periodic( &dom, &aPredicate, &l2, {true, true} );
503 
504  //L_infinity metric
505  //typedef DistanceTransformation<TSpace,Predicate, 0> DT;
506  //DT dt(Domain(a,b), aPredicate);;
507 
508  //L_1 metric
510  L1Metric l1;
512  DT1 dt1(&dom,&aPredicate,&l1);
513  DT1 dt1_periodic( &dom, &aPredicate, &l1, {true, true} );
514 
515  DGtal::int64_t maxv = 0;
516  for ( DistanceTransformation<TSpace,Predicate, L2Metric>::ConstRange::ConstIterator it = dt2.constRange().begin(), itend = dt2.constRange().end();it != itend; ++it)
517  if ( (*it) > maxv)
518  maxv = (*it);
519 
520  trace.warning() << dt2 << "MaxV = " << maxv << endl;
521  //We display the values on a 2D slice
522  for (unsigned int y = 0; y < 16; y++)
523  {
524  for (unsigned int x = 0; x < 16; x++)
525  {
526  Point p(x, y);
527  std::cout << std::setw(4) << dt2(p) << " ";
528  }
529  std::cout << std::endl;
530  }
531 
532  trace.info()<< "Exporting to SVG"<<endl;
533 
534  Board2D board;
536  Display2DFactory::drawImage<Hue>(board, dt2, (DGtal::int64_t)0, (DGtal::int64_t)maxv+1);
537  board.saveSVG ( "image-DT-l2.svg" );
538 
539  board.clear();
540  Display2DFactory::drawImage<Hue>(board, dt2_periodic, (DGtal::int64_t)0, (DGtal::int64_t)maxv+1);
541  board.saveSVG ( "image-DT-l2-periodic.svg" );
542 
543  trace.info()<< "done"<<endl;
544 
545 
546 
547  trace.info()<< "max L1"<<endl;
548  maxv = 0;
550  itend = dt1.constRange().end();
551  it2 != itend; ++it2)
552  {
553  if ( *it2 > maxv)
554  maxv = (*it2);
555  }
556 
557  trace.info()<< "Exporting to SVG L1"<<endl;
558  board.clear();
559  Display2DFactory::drawImage<Hue>(board, dt1, (DGtal::int64_t)0, (DGtal::int64_t)maxv+1);
560  board.saveSVG ( "image-DT-l1.svg" );
561 
562  board.clear();
563  Display2DFactory::drawImage<Hue>(board, dt1_periodic, (DGtal::int64_t)0, (DGtal::int64_t)maxv+1);
564  board.saveSVG ( "image-DT-l1-periodic.svg" );
565 
566  trace.info()<< "done"<<endl;
567 
568 
569  trace.endBlock();
570 
571  return nbok == nb;
572 }
573 
574 template <typename Space, int norm>
575 bool testCompareExactInexact(unsigned int size, unsigned int nb)
576 {
577  trace.beginBlock("Checking Exact/Inexct predicate metrics");
579  typedef InexactPredicateLpSeparableMetric<Space> MetricInex;
581  typedef typename Space::Point Point;
582  typedef DigitalSetBySTLSet<Domain> Set;
583  // typedef NotPointPredicate<SetPredicate<Set> > NegPredicate;
584  typedef functors::NotPointPredicate<Set> NegPredicate;
585 
586  Point low=Point::diagonal(0),
587  up=Point::diagonal(size);
588 
589  Domain domain(low,up);
590  Set set(domain);
591 
592  for(unsigned int i = 0; i<nb; ++i)
593  {
594  Point p;
595  for(unsigned int dim=0; dim<Space::dimension;++dim)
596  p[dim] = rand() % size;
597  set.insert(p);
598  }
599 
600  trace.info()<< "Testing metrics "<<MetricEx()<<" "<<MetricInex(norm)<<std::endl;
601  trace.info()<< "Testing space dimension "<<Space::dimension<<std::endl;
602  trace.info()<< "Inserting "<<set.size() << " points."<<std::endl;
603 
604  // SetPredicate<Set> setPred(set);
605  NegPredicate negPred(set);
606 
609  MetricEx metricEx;
610  MetricInex metricInex(norm);
611  DTEx dtex(&domain, &negPred, &metricEx);
612  DTIn dtinex(&domain, &negPred, &metricInex);
613 
614  double MSE=0.0;
615  typename DTEx::ConstRange::ConstIterator it=dtex.constRange().begin(), itend=dtex.constRange().end();
616  typename DTIn::ConstRange::ConstIterator it2 = dtinex.constRange().begin();
617  for( ; it != itend; ++it, ++it2)
618  MSE += ((*it) - (*it2))*((*it) - (*it2));
619 
620  trace.warning()<<"Resulting MSE= "<<MSE;
621  trace.endBlock();
622  return true;
623 }
624 
626 // Standard services - public :
627 
628 int main ( int argc, char** argv ){
629 
630  trace.beginBlock ( "Testing class DistanceTransformation" );
631  trace.info() << "Args:";
632  for ( int i = 0; i < argc; ++i )
633  trace.info() << " " << argv[ i ];
634  trace.info() << endl;
635 
637  && testDTFromSet()
640  && testChessboard()
641  && testDTFromSet()
642  && testCompareExactInexact<Z2i::Space, 2>(50, 50)
643  && testCompareExactInexact<Z3i::Space, 2>(50, 50)
644  && testCompareExactInexact<Z2i::Space, 4>(50, 50)
645  && testCompareExactInexact<Z3i::Space, 4>(50, 50)
646  ;
647  //&& ... other tests
648  trace.emphase() << ( res ? "Passed." : "Error." ) << endl;
649  trace.endBlock();
650  return res ? 0 : 1;
651 }
652 // //
Aim: Model of the concept StarShaped represents any accelerated flower in the plane.
Definition: AccFlower2D.h:65
RealPoint getLowerBound() const
Definition: AccFlower2D.h:134
RealPoint getUpperBound() const
Definition: AccFlower2D.h:143
Aim: This class specializes a 'Board' class so as to display DGtal objects more naturally (with <<)....
Definition: Board2D.h:71
This class adapts any iterator so that operator* returns another element than the one pointed to by t...
Aim: A wrapper class around a STL associative container for storing sets of digital points within som...
Aim: A container class for storing sets of digital points within some given domain.
Aim: Implementation of the linear in time distance transformation for separable metrics.
SeparableMetric::Value Value
Definition of the image value type.
Aim: implements separable l_p metrics with exact predicates.
Aim: This class template may be used to (linearly) convert scalar values in a given range into gray l...
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 & end() const
const ConstIterator & begin() const
Aim: implements association bewteen points lying in a digital domain and values.
Definition: Image.h:70
Aim: implements separable l_p metrics with approximated predicates.
Aim: A utility class for constructing different shapes (balls, diamonds, and others).
std::ostream & error()
void beginBlock(const std::string &keyword="")
std::ostream & emphase()
std::ostream & info()
std::ostream & warning()
double endBlock()
Aim: Implementation of the linear in time Voronoi map construction.
Definition: VoronoiMap.h:127
Aim: Define a simple Foreground predicate thresholding image values given a single thresold....
void clear(const DGtal::Color &color=DGtal::Color::None)
Definition: Board.cpp:151
void saveSVG(const char *filename, PageSize size=Board::BoundingBox, double margin=10.0) const
Definition: Board.cpp:1011
void setUnit(Unit unit)
Definition: Board.cpp:239
MyDigitalSurface::ConstIterator ConstIterator
functors namespace gathers all DGtal functors.
DGtal is the top-level namespace which contains all DGtal functions and types.
boost::int64_t int64_t
signed 94-bit integer.
Definition: BasicTypes.h:74
DGtal::uint32_t Dimension
Definition: Common.h:136
Trace trace
Definition: Common.h:153
Aim: The predicate returns true when the point predicate given at construction return false....
MyPointD Point
Definition: testClone2.cpp:383
bool testChessboard()
bool testDistanceTransformationBorder()
int main(int argc, char **argv)
bool testCompareExactInexact(unsigned int size, unsigned int nb)
bool testDistanceTransformation()
void randomSeeds(Image &input, const unsigned int nb, const int value)
bool testDistanceTransformationNeg()
bool testDTFromSet()
bool testDistanceTransformation3D()
Domain domain
Image image(domain)
ImageContainerBySTLVector< Domain, Value > Image
HyperRectDomain< Space > Domain
unsigned int dim(const Vector &z)