DGtal  0.9.2
testDistanceTransformation.cpp
1 
30 #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 
85 bool testDistanceTransformation()
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.getVoronoiVector(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.getVoronoiVector(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 }
201 bool testDistanceTransformationNeg()
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);
248  DistanceTransformation<TSpace, Predicate , L2Metric> dt(&dom, &aPredicate, &l2);
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 
287 bool testDTFromSet()
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 
344 bool testDistanceTransformationBorder()
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;
353  typedef HyperRectDomain<TSpace> Domain;
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);
374  DistanceTransformation<TSpace, Predicate, L2Metric> dt(&dom, &aPredicate, &l2);
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 
417 bool testDistanceTransformation3D()
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;
426  typedef HyperRectDomain<TSpace> Domain;
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;
446  DistanceTransformation<TSpace, Predicate, L2Metric> dt(&dom, &aPredicate,&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 
468 bool testChessboard()
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;
477  typedef HyperRectDomain<TSpace> Domain;
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 
503  //L_infinity metric
504  //typedef DistanceTransformation<TSpace,Predicate, 0> DT;
505  //DT dt(Domain(a,b), aPredicate);;
506 
507  //L_1 metric
509  L1Metric l1;
511  DT1 dt1(&dom,&aPredicate,&l1);
512 
513  DGtal::int64_t maxv = 0;
514  for ( DistanceTransformation<TSpace,Predicate, L2Metric>::ConstRange::ConstIterator it = dt2.constRange().begin(), itend = dt2.constRange().end();it != itend; ++it)
515  if ( (*it) > maxv)
516  maxv = (*it);
517 
518  trace.warning() << dt2 << "MaxV = " << maxv << endl;
519  //We display the values on a 2D slice
520  for (unsigned int y = 0; y < 16; y++)
521  {
522  for (unsigned int x = 0; x < 16; x++)
523  {
524  Point p(x, y);
525  std::cout << std::setw(4) << dt2(p) << " ";
526  }
527  std::cout << std::endl;
528  }
529 
530  trace.info()<< "Exporting to SVG"<<endl;
531 
532  Board2D board;
534  Display2DFactory::drawImage<Hue>(board, dt2, (DGtal::int64_t)0, (DGtal::int64_t)maxv+1);
535  board.saveSVG ( "image-DT-l2.svg" );
536  trace.info()<< "done"<<endl;
537 
538 
539 
540  trace.info()<< "max L1"<<endl;
541  maxv = 0;
543  itend = dt1.constRange().end();
544  it2 != itend; ++it2)
545  {
546  if ( *it2 > maxv)
547  maxv = (*it2);
548  }
549 
550  trace.info()<< "Exporting to SVG L1"<<endl;
551  board.clear();
552  Display2DFactory::drawImage<Hue>(board, dt1, (DGtal::int64_t)0, (DGtal::int64_t)maxv+1);
553  board.saveSVG ( "image-DT-l1.svg" );
554  trace.info()<< "done"<<endl;
555 
556 
557  trace.endBlock();
558 
559  return nbok == nb;
560 }
561 
562 template <typename Space, int norm>
563 bool testCompareExactInexact(unsigned int size, unsigned int nb)
564 {
565  trace.beginBlock("Checking Exact/Inexct predicate metrics");
567  typedef InexactPredicateLpSeparableMetric<Space> MetricInex;
568  typedef HyperRectDomain<Space> Domain;
569  typedef typename Space::Point Point;
570  typedef DigitalSetBySTLSet<Domain> Set;
571  // typedef NotPointPredicate<SetPredicate<Set> > NegPredicate;
572  typedef functors::NotPointPredicate<Set> NegPredicate;
573 
574  Point low=Point::diagonal(0),
575  up=Point::diagonal(size);
576 
577  Domain domain(low,up);
578  Set set(domain);
579 
580  for(unsigned int i = 0; i<nb; ++i)
581  {
582  Point p;
583  for(unsigned int dim=0; dim<Space::dimension;++dim)
584  p[dim] = rand() % size;
585  set.insert(p);
586  }
587 
588  trace.info()<< "Testing metrics "<<MetricEx()<<" "<<MetricInex(norm)<<std::endl;
589  trace.info()<< "Testing space dimension "<<Space::dimension<<std::endl;
590  trace.info()<< "Inserting "<<set.size() << " points."<<std::endl;
591 
592  // SetPredicate<Set> setPred(set);
593  NegPredicate negPred(set);
594 
597  MetricEx metricEx;
598  MetricInex metricInex(norm);
599  DTEx dtex(&domain, &negPred, &metricEx);
600  DTIn dtinex(&domain, &negPred, &metricInex);
601 
602  double MSE=0.0;
603  typename DTEx::ConstRange::ConstIterator it=dtex.constRange().begin(), itend=dtex.constRange().end();
604  typename DTIn::ConstRange::ConstIterator it2 = dtinex.constRange().begin();
605  for( ; it != itend; ++it, ++it2)
606  MSE += ((*it) - (*it2))*((*it) - (*it2));
607 
608  trace.warning()<<"Resulting MSE= "<<MSE;
609  trace.endBlock();
610  return true;
611 }
612 
614 // Standard services - public :
615 
616 int main ( int argc, char** argv ){
617 
618  trace.beginBlock ( "Testing class DistanceTransformation" );
619  trace.info() << "Args:";
620  for ( int i = 0; i < argc; ++i )
621  trace.info() << " " << argv[ i ];
622  trace.info() << endl;
623 
624  bool res = testDistanceTransformation() && testDistanceTransformationNeg()
625  && testDTFromSet()
626  && testDistanceTransformationBorder()
627  && testDistanceTransformation3D()
628  && testChessboard()
629  && testDTFromSet()
630  && testCompareExactInexact<Z2i::Space, 2>(50, 50)
631  && testCompareExactInexact<Z3i::Space, 2>(50, 50)
632  && testCompareExactInexact<Z2i::Space, 4>(50, 50)
633  && testCompareExactInexact<Z3i::Space, 4>(50, 50)
634  ;
635  //&& ... other tests
636  trace.emphase() << ( res ? "Passed." : "Error." ) << endl;
637  trace.endBlock();
638  return res ? 0 : 1;
639 }
640 // //
void beginBlock(const std::string &keyword="")
Aim: The predicate returns true when the point predicate given at construction return false...
Aim: Implementation of the linear in time Voronoi map construction.
Definition: VoronoiMap.h:113
ExactPredicateLpSeparableMetric< Space, 1 > L1Metric
Definition: StdDefs.h:119
static const Dimension dimension
static constants to store the dimension.
Definition: SpaceND.h:132
This class adapts any iterator so that operator* returns another element than the one pointed to by t...
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
ExactPredicateLpSeparableMetric< Space, 2 > L2Metric
Definition: StdDefs.h:118
Aim: SpaceND is a utility class that defines the fundamental structure of a Digital Space in ND...
Definition: SpaceND.h:95
DGtal::uint32_t Dimension
Definition: Common.h:113
Aim: implements separable l_p metrics with approximated predicates.
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.
void setValue(const Point &aPoint, const Value &aValue)
Definition: Image.h:247
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: Define a simple Foreground predicate thresholding image values given a single thresold...
std::ostream & emphase()
Aim: Model of the concept StarShaped represents any accelerated flower in the plane.
Definition: AccFlower2D.h:64
DGtal is the top-level namespace which contains all DGtal functions and types.
Aim: A wrapper class around a STL associative container for storing sets of digital points within som...
std::ostream & warning()
std::ostream & info()
const Domain & domain() const
Definition: Image.h:192
void setUnit(Unit unit)
Definition: Board.cpp:240
Aim: A utility class for constructing different shapes (balls, diamonds, and others).
std::ostream & error()
boost::int64_t int64_t
signed 94-bit integer.
Definition: BasicTypes.h:74
Aim: A container class for storing sets of digital points within some given domain.
SeparableMetric::Value Value
Definition of the image value type.
Aim: This class specializes a 'Board' class so as to display DGtal objects more naturally (with <<)...
Definition: Board2D.h:70