DGtal  0.9.2
testVoronoiMap.cpp
1 
30 #include <iostream>
32 #include "DGtal/base/Common.h"
33 #include "DGtal/helpers/StdDefs.h"
34 #include "DGtal/images/CConstImage.h"
35 #include "DGtal/geometry/volumes/distance/VoronoiMap.h"
36 #include "DGtal/geometry/volumes/distance/ExactPredicateLpSeparableMetric.h"
37 #include "DGtal/geometry/volumes/distance/InexactPredicateLpSeparableMetric.h"
38 #include "DGtal/geometry/volumes/distance/DistanceTransformation.h"
39 #include "DGtal/kernel/BasicPointPredicates.h"
40 #include "DGtal/io/boards/Board2D.h"
41 #include "DGtal/io/colormaps/HueShadeColorMap.h"
43 
44 using namespace std;
45 using namespace DGtal;
47 // Functions for testing class VoronoiMap.
49 
50 
51 template<typename Point>
52 double mynorm(const Point &point, const double p)
53 {
54  double res=0.0;
55  for(unsigned int i=0; i< Point::dimension; i++)
56  res += std::pow ( (double)abs(point[i]) , p);
57 
58  return std::pow(res, 1.0/p);
59 }
60 
61 template <typename VoroMap>
62 void saveVoroMap(const std::string &filename,const VoroMap &output,const double p)
63 {
64  typedef HueShadeColorMap<double,2> Hue;
65  double maxdt=0.0;
66 
67  for ( typename VoroMap::Domain::ConstIterator it = output.domain().begin(), itend = output.domain().end();
68  it != itend; ++it)
69  {
70  typename VoroMap::Value point = output(*it);
71  if ( mynorm(point-(*it),p) > maxdt)
72  maxdt = mynorm(point-(*it),p);
73  }
74  trace.error() << "MaxDT="<<maxdt<<std::endl;
75 
76  Board2D board;
77  Hue hue(0,maxdt);
78 
79  for(typename VoroMap::Domain::ConstIterator it = output.domain().begin(),
80  itend = output.domain().end();
81  it != itend; ++it)
82  {
83  typename VoroMap::Value point = output(*it);
84  board << CustomStyle( (*it).className(), new CustomColors( hue(mynorm(point- (*it),p)),
85  hue(mynorm(point- (*it),p))))
86  << (*it);
87  }
88 
89  board.saveSVG(filename.c_str());
90 }
91 
92 
93 /* Is Validate the VoronoiMap
94  */
95 template < typename Set, typename Image>
96 bool checkVoronoiL2(const Set &aSet, const Image & voro)
97 {
98  typedef typename Image::Point Point;
99 
100  for(typename Image::Domain::ConstIterator it = voro.domain().begin(), itend = voro.domain().end();
101  it != itend; ++it)
102  {
103  Point psite = voro(*it);
104  Point p = (*it);
105  DGtal::int64_t d=0;
106  for(DGtal::Dimension i=0; i<Point::dimension;++i)
107  d+= (p[i]-psite[i])*(p[i]-psite[i]);
108 
109  for(typename Set::ConstIterator itset = aSet.begin(), itendSet = aSet.end();
110  itset != itendSet;
111  ++itset)
112  {
113  DGtal::int64_t dbis=0;
114  for(DGtal::Dimension i=0; i<Point::dimension;++i)
115  dbis+= (p[i]-(*itset)[i])*(p[i]-(*itset)[i]);
116  if ( dbis < d)
117  {
118  trace.error() << "DT Error at "<<p<<" Voro:"<<psite<<" ("<<d<<") from set:"
119  << (*itset) << "("<<dbis<<")"<<std::endl;
120  return false;
121  }
122  }
123  }
124  return true;
125 }
126 
127 
128 bool testCheckConcept()
129 {
132 
133  return true;
134 }
135 
140 bool testVoronoiMap()
141 {
142  unsigned int nbok = 0;
143  unsigned int nb = 0;
144 
145  trace.beginBlock ( "Testing VoronoiMap2D ..." );
146 
147  Z2i::Point a(-10,-10);
148  Z2i::Point b(10,10);
149  Z2i::Domain domain(a,b);
150 
151  Z2i::DigitalSet mySet(domain);
152 
153  for(Z2i::Domain::ConstIterator it = domain.begin(), itend = domain.end();
154  it != itend;
155  ++it)
156  mySet.insertNew( *it );
157 
158 
159  Z2i::DigitalSet sites(domain);
160 
161  sites.insertNew( Z2i::Point(0,-6));
162  sites.insertNew( Z2i::Point(6,0));
163  sites.insertNew( Z2i::Point(-6,0));
164 
165  for(Z2i::DigitalSet::ConstIterator it = sites.begin(), itend = sites.end();
166  it != itend; ++it)
167  mySet.erase (*it);
168 
169 
170 
173  L2Metric l2;
174  Voro2 voro(&domain, &mySet,&l2);
175 
176  for(int j=-10; j <= 10; j++)
177  {
178  for(int i=-10; i<=10; i++)
179  trace.info() << "("<<voro( Z2i::Point(i,j))[0]<<","<< voro( Z2i::Point(i,j))[1]<<") ";
180  trace.info()<<std::endl;
181  }
182 
183  trace.info()<<"Exporting o SVG"<<std::endl;
184 
185  Board2D board;
186  for(Voro2::OutputImage::Domain::ConstIterator it = voro.domain().begin(),
187  itend = voro.domain().end();
188  it != itend; ++it)
189  {
190  Z2i::Point p = voro(*it);
191  unsigned char c = (p[1]*13 + p[0] * 7) % 256;
192  board << CustomStyle( (*it).className(), new CustomColors(Color(c,c,c),Color(c,c,c)))
193  << (*it);
194  }
195 
196  board.saveSVG("Voromap.svg");
197 
198  nbok += checkVoronoiL2(sites,voro) ? 1 : 0;
199  nb++;
200  trace.info() << "(" << nbok << "/" << nb << ") "
201  << "Voronoi diagram is valid !" << std::endl;
202  trace.endBlock();
203 
204 
205 
206 
207  return nbok == nb;
208 }
209 
210 
211 
216 template<typename Set>
217 bool testVoronoiMapFromSites2D(const Set &aSet, const std::string &name)
218 {
219  unsigned int nbok = 0;
220  unsigned int nb = 0;
221 
222  Set mySet(aSet.domain());
223 
224  for(typename Set::Domain::ConstIterator it = aSet.domain().begin(), itend = aSet.domain().end();
225  it != itend;
226  ++it)
227  mySet.insertNew( *it );
228 
229 
230  for(typename Set::ConstIterator it = aSet.begin(), itend = aSet.end();
231  it != itend; ++it)
232  mySet.erase (*it);
233 
234 
235  trace.beginBlock(" Voro computation");
238  L2Metric l2;
239  Voro2 voro(aSet.domain(), mySet, l2 );
240 
241  trace.endBlock();
242 
243  trace.beginBlock(" Voronoi computation l_6");
245  L6Metric l6;
247  Voro6 voro6(aSet.domain(), mySet, l6 );
248  trace.endBlock();
249 
250 
251 
252  trace.beginBlock(" DT computation");
254  DT dt(aSet.domain(), mySet, l2);
255  trace.endBlock();
256 
257 
258  if ( (aSet.domain().upperBound()[1] - aSet.domain().lowerBound()[1]) <20)
259  {
260  for(int j= aSet.domain().lowerBound()[1]; j <= aSet.domain().upperBound()[1]; j++)
261  {
262  for(int i=aSet.domain().lowerBound()[0]; i<=aSet.domain().upperBound()[0]; i++)
263  if ( aSet.find( Z2i::Point(i,j) ) != aSet.end() )
264  std::cout <<"X ";
265  else
266  std::cout<<"0 ";
267  trace.info()<<std::endl;
268  }
269 
270  trace.info() << std::endl;
271 
272  for(int j= aSet.domain().lowerBound()[1]; j <= aSet.domain().upperBound()[1]; j++)
273  {
274  for(int i=aSet.domain().lowerBound()[0]; i<=aSet.domain().upperBound()[0]; i++)
275  trace.info() << "("<<voro( Z2i::Point(i,j))[0]<<","<< voro( Z2i::Point(i,j))[1]<<") ";
276  trace.info()<<std::endl;
277  }
278  }
279 
280  Board2D board;
281  board << voro.domain();
282  for(typename Voro2::OutputImage::Domain::ConstIterator it = voro.domain().begin(), itend = voro.domain().end();
283  it != itend; ++it)
284  {
285  if (!mySet(*it))
286  board << (*it);
287  }
288  std::string filename= "Voromap-"+name+"-orig.svg";
289  board.saveSVG(filename.c_str());
290 
291  board.clear();
292  board << voro.domain();
293  board.setPenColor(Color(0,0,0));
294  for(typename Voro2::OutputImage::Domain::ConstIterator it = voro.domain().begin(), itend = voro.domain().end();
295  it != itend; ++it)
296  {
297  Z2i::Point p = voro(*it);
298  if ((p != (*it)) && (p != voro.domain().upperBound() + Z2i::Point::diagonal(1))
299  && (p != voro.domain().lowerBound()))
300  Display2DFactory::draw( board, p - (*it), (*it));
301  }
302 
303  filename= "Voromap-"+name+"-diag.svg";
304  board.saveSVG(filename.c_str());
305 
306 
307  board.clear();
308  board << voro.domain();
309  for(typename Voro2::OutputImage::Domain::ConstIterator it = voro.domain().begin(), itend = voro.domain().end();
310  it != itend; ++it)
311  {
312  Z2i::Point p = voro(*it);
313  unsigned char c = (p[1]*13 + p[0] * 7) % 256;
314  // if ((p != (*it)) && (p != voro.domain().upperBound() + Z2i::Point::diagonal(1))
315  // && (p != voro.domain().lowerBound()))
316  board << CustomStyle( (*it).className(), new CustomColors(Color(c,c,c),Color(c,c,c)))
317  << (*it);;
318  }
319 
320  filename= "Voromap-"+name+".svg";
321  board.saveSVG(filename.c_str());
322  filename= "Voromap-"+name+"-hue.svg";
323  saveVoroMap(filename.c_str(),voro,2);
324 
325 
326 
327  board.clear();
328  for(typename Voro6::OutputImage::Domain::ConstIterator it = voro6.domain().begin(),
329  itend = voro6.domain().end();
330  it != itend; ++it)
331  {
332  Z2i::Point p = voro6(*it);
333  if (p != (*it))
334  Display2DFactory::draw( board, p - (*it), (*it));
335  }
336 
337  filename= "Voromap-diag-l6-"+name+".svg";
338  board.saveSVG(filename.c_str());
339 
340  board.clear();
341  for(typename Voro6::OutputImage::Domain::ConstIterator it = voro6.domain().begin(), itend = voro6.domain().end();
342  it != itend; ++it)
343  {
344  Z2i::Point p = voro6(*it);
345  unsigned char c = (p[1]*13 + p[0] * 7) % 256;
346  board << CustomStyle( (*it).className(), new CustomColors(Color(c,c,c),Color(c,c,c)))
347  << (*it);;
348  }
349 
350  filename= "Voromap-l6"+name+".svg";
351  board.saveSVG(filename.c_str());
352  filename= "Voromap-hue-l6-"+name+".svg";
353  saveVoroMap(filename.c_str(),voro6,3);
354 
355 
356  nbok += checkVoronoiL2(aSet,voro) ? 1 : 0;
357  nb++;
358  trace.info() << "(" << nbok << "/" << nb << ") "
359  << "Voronoi diagram is valid !" << std::endl;
360 
361  return nbok == nb;
362 }
363 
368 template<typename Set>
369 bool testVoronoiMapFromSites(const Set &aSet)
370 {
371  unsigned int nbok = 0;
372  unsigned int nb = 0;
373 
374  Set mySet(aSet.domain());
375 
376  for(typename Set::Domain::ConstIterator it = aSet.domain().begin(),
377  itend = aSet.domain().end();
378  it != itend;
379  ++it)
380  mySet.insertNew( *it );
381 
382 
383  for(typename Set::ConstIterator it = aSet.begin(), itend = aSet.end();
384  it != itend; ++it)
385  mySet.erase (*it);
386 
387 
388 
389  trace.beginBlock(" Voronoi computation");
392  L2Metric l2;
393  Voro2 voro(aSet.domain(), mySet, l2);
394  trace.endBlock();
395 
396 
397  trace.beginBlock(" Voronoi computation l_3");
400  L3Metric l3;
401  Voro3 voro3(aSet.domain(), mySet, l3);
402  trace.endBlock();
403 
404 
405  trace.beginBlock(" DT computation");
407  DT dt(aSet.domain(), mySet, l2);
408 
409  trace.endBlock();
410 
411 
412  trace.beginBlock("Validating the Voronoi Map");
413  nbok += (checkVoronoiL2(aSet,voro) )? 1 : 0;
414  trace.endBlock();
415  nb++;
416  trace.info() << "(" << nbok << "/" << nb << ") "
417  << "Voronoi diagram is valid !" << std::endl;
418 
419  return nbok == nb;
420 }
421 
422 
423 bool testSimple2D()
424 {
425 
426  Z2i::Point a(-10,-10);
427  Z2i::Point b(10,10);
428  Z2i::Domain domain(a,b);
429 
430  Z2i::DigitalSet sites(domain);
431  bool ok;
432 
433  trace.beginBlock("Simple2D");
434  sites.insertNew( Z2i::Point(0,-6));
435  sites.insertNew( Z2i::Point(6,0));
436  sites.insertNew( Z2i::Point(-6,0));
437 
438  ok = testVoronoiMapFromSites2D<Z2i::DigitalSet>(sites,"simple");
439  trace.endBlock();
440 
441  return ok;
442 
443 }
444 
445 bool testSimpleRandom2D()
446 {
447  unsigned size=16;
448  Z2i::Point a(0,0);
449  Z2i::Point b(size,size);
450  Z2i::Domain domain(a,b);
451 
452  Z2i::DigitalSet sites(domain);
453  bool ok;
454 
455  trace.beginBlock("Random 2D");
456  for(unsigned int i = 0 ; i < size; ++i)
457  {
458  Z2i::Point p( rand() % (b[0]) - a[0], rand() % (b[1]) + a[1] );
459  sites.insert( p );
460  }
461  ok = testVoronoiMapFromSites2D<Z2i::DigitalSet>(sites,"random");
462  trace.endBlock();
463 
464  trace.beginBlock("Random 2D (dense)");
465  for(unsigned int i = 0 ; i < size*size-size; ++i)
466  {
467  Z2i::Point p( rand() % (b[0]) - a[0], rand() % (b[1]) + a[1] );
468  sites.insert( p );
469  }
470  ok = ok && testVoronoiMapFromSites2D<Z2i::DigitalSet>(sites,"random-dense");
471  trace.endBlock();
472 
473  return ok;
474 
475 }
476 
477 
478 bool testSimple3D()
479 {
480 
481  Z3i::Point a(-10,-10,-10);
482  Z3i::Point b(10,10,10);
483  Z3i::Domain domain(a,b);
484 
485  Z3i::DigitalSet sites(domain);
486  bool ok;
487 
488  trace.beginBlock("Simple3D");
489  sites.insertNew( Z3i::Point(0,0,-6));
490  sites.insertNew( Z3i::Point(6,0,0));
491  sites.insertNew( Z3i::Point(-6,0,3));
492 
493  ok = testVoronoiMapFromSites<Z3i::DigitalSet>(sites);
494  trace.endBlock();
495 
496  return ok;
497 
498 }
499 
500 bool testSimpleRandom3D()
501 {
502 
503  Z3i::Point a(0,0,0);
504  Z3i::Point b(64,64,64);
505  Z3i::Domain domain(a,b);
506 
507  Z3i::DigitalSet sites(domain);
508  bool ok;
509 
510  trace.beginBlock("Random 3D");
511  for(unsigned int i = 0 ; i < 64; ++i)
512  {
513  Z3i::Point p( rand() % (b[0]) - a[0],
514  rand() % (b[1]) + a[1],
515  rand() % (b[2]) + a[2] );
516  sites.insert( p );
517  }
518  ok = testVoronoiMapFromSites<Z3i::DigitalSet>(sites);
519  trace.endBlock();
520 
521  return ok;
522 
523 }
524 
525 
526 
527 bool testSimple4D()
528 {
529 
530  typedef SpaceND<4> Space4;
531  Space4::Point a(0,0,0,0);
532  Space4::Point b(5,5,5,5);
533  HyperRectDomain<Space4> domain(a,b);
534 
536  bool ok;
537 
538  trace.beginBlock("Simple4D");
539  sites.insertNew( Space4::Point(1,4,1,1));
540  sites.insertNew( Space4::Point(3,1,3,1));
541  sites.insertNew( Space4::Point(0,0,0,0));
542 
543  ok = testVoronoiMapFromSites< DigitalSetBySTLSet< HyperRectDomain<Space4> > >(sites);
544  trace.endBlock();
545 
546  return ok;
547 
548 }
549 
550 
552 // Standard services - public :
553 
554 int main( int argc, char** argv )
555 {
556  trace.beginBlock ( "Testing class VoronoiMap" );
557  trace.info() << "Args:";
558  for ( int i = 0; i < argc; ++i )
559  trace.info() << " " << argv[ i ];
560  trace.info() << endl;
561 
562  bool res = testCheckConcept()
563  && testVoronoiMap()
564  && testSimple2D()
565  && testSimpleRandom2D()
566  && testSimple3D()
567  && testSimpleRandom3D()
568  && testSimple4D()
569  ; // && ... other tests
570  trace.emphase() << ( res ? "Passed." : "Error." ) << endl;
571  trace.endBlock();
572  return res ? 0 : 1;
573 }
574 // //
void beginBlock(const std::string &keyword="")
Aim: Implementation of the linear in time Voronoi map construction.
Definition: VoronoiMap.h:113
const ConstIterator & begin() const
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
const ConstIterator & end() const
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
STL namespace.
double endBlock()
Custom style class redefining the pen color and the fill color. You may use Board2D::Color::None for ...
Definition: Board2D.h:278
Aim: This class template may be used to (linearly) convert scalar values in a given range into a colo...
Aim: implements separable l_p metrics with exact predicates.
std::ostream & emphase()
void clear(const DGtal::Color &color=DGtal::Color::None)
Definition: Board.cpp:152
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...
void saveSVG(const char *filename, PageSize size=Board::BoundingBox, double margin=10.0) const
Definition: Board.cpp:1012
std::ostream & info()
Structure representing an RGB triple with alpha component.
Definition: Color.h:66
Board & setPenColor(const DGtal::Color &color)
Definition: Board.cpp:298
const Domain & domain() const
Definition: Image.h:192
Aim: Defines the concept describing a read-only image, which is a refinement of CPointFunctor.
Definition: CConstImage.h:94
std::ostream & error()
Container::const_iterator ConstIterator
ConstIterator type of the container;.
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.
Aim: This class specializes a 'Board' class so as to display DGtal objects more naturally (with <<)...
Definition: Board2D.h:70