DGtal 1.3.0
Loading...
Searching...
No Matches
testVoronoiMap.cpp
Go to the documentation of this file.
1
31#include <iostream>
32#include <array>
33#include <algorithm>
34
35#include "DGtal/base/Common.h"
36#include "DGtal/helpers/StdDefs.h"
37#include "DGtal/images/CConstImage.h"
38#include "DGtal/geometry/volumes/distance/VoronoiMap.h"
39#include "DGtal/geometry/volumes/distance/ExactPredicateLpSeparableMetric.h"
40#include "DGtal/geometry/volumes/distance/InexactPredicateLpSeparableMetric.h"
41#include "DGtal/geometry/volumes/distance/DistanceTransformation.h"
42#include "DGtal/kernel/BasicPointPredicates.h"
43#include "DGtal/io/boards/Board2D.h"
44#include "DGtal/io/colormaps/HueShadeColorMap.h"
46
47using namespace std;
48using namespace DGtal;
50// Functions for testing class VoronoiMap.
52
53
54template<typename Point>
55double mynorm(const Point &point, const double p)
56{
57 double res=0.0;
58 for(unsigned int i=0; i< Point::dimension; i++)
59 res += std::pow ( (double)abs(point[i]) , p);
60
61 return std::pow(res, 1.0/p);
62}
63
64template <typename VoroMap>
65void saveVoroMap(const std::string &filename,const VoroMap &output,const double p)
66{
68
69 const double maxdt = mynorm( output.domain().upperBound() - output.domain().lowerBound(), p );
70
71 Board2D board;
72 Hue hue(0, maxdt);
73
74 for(typename VoroMap::Domain::ConstIterator it = output.domain().begin(),
75 itend = output.domain().end();
76 it != itend; ++it)
77 {
78 typename VoroMap::Value point = output(*it);
79 board << CustomStyle( (*it).className(), new CustomColors( hue(mynorm(point- (*it),p)),
80 hue(mynorm(point- (*it),p))))
81 << (*it);
82 }
83
84 board.saveSVG(filename.c_str());
85}
86
95template < typename Point, typename Domain >
97{
98 auto const & lowerBound = aDomain.lowerBound();
99 auto const & upperBound = aDomain.upperBound();
100
101 for ( std::size_t i = 0; i < Domain::dimension; ++i )
102 aPoint[i] = ( aPoint[i] - 2*lowerBound[i] + upperBound[i] + 1 ) % ( upperBound[i] - lowerBound[i] + 1 ) + lowerBound[i];
103
104 return aPoint;
105}
106
117template < std::size_t N >
118std::array<bool, N> getPeriodicityFromInteger( std::size_t anInteger )
119{
120 std::array<bool, N> periodicity;
121 for ( std::size_t i = 0, mask = 1 ; i < N ; ++i, mask *= 2 )
122 periodicity[i] = anInteger & mask;
123
124 return periodicity;
125}
126
133template < std::size_t N >
134std::string formatPeriodicity( std::array<bool, N> const& aPeriodicity )
135{
136 std::string str;
137 for ( std::size_t i = 0; i < N; ++i )
138 str += aPeriodicity[i] ? '1' : '0';
139
140 return str;
141}
142
145template < typename Set, typename Voro>
146bool checkVoronoi(const Set &aSet, const Voro & voro)
147{
148 typedef typename Voro::Point Point;
149
150 // Domain extent
151 Point const extent = voro.domain().upperBound() - voro.domain().lowerBound() + Point::diagonal();
152
153 // Site shifting depending on the periodicity.
154 std::vector< typename Voro::PeriodicitySpec > periodicShift;
155 for ( std::size_t i = 0; i < ( 1u << Voro::Space::dimension ); ++i )
156 {
157 const auto periodicity = getPeriodicityFromInteger< Voro::Space::dimension >( i );
158
159 // Checking if this periodicity possibility is valid.
160 bool isValid = true;
161 for ( size_t j = 0; j < periodicity.size(); ++j )
162 if ( periodicity[j] && ! voro.isPeriodic(j) )
163 {
164 isValid = false;
165 break;
166 }
167
168 if ( isValid )
169 periodicShift.push_back( periodicity );
170 }
171
172 // Checking site for all domain points.
173 for ( auto const& pt : voro.domain() )
174 {
175 // Calculating reference (raw-)distance.
176 const Point psite = voro(pt);
177 const auto dist = voro.metric()->rawDistance( pt, psite );
178
179 // Checking that the reference site is actually a given site.
180 // Also testing projectPoint method.
181 if ( std::find( aSet.begin(), aSet.end(), voro.projectPoint( psite ) ) == aSet.end() )
182 {
183 trace.error() << "The Voro point " << psite
184 << " projected to " << voro.projectPoint( psite )
185 << " is not a valid site." << std::endl;
186 return false;
187 }
188
189 // Checking if we can found a better site.
190 for ( auto site : aSet )
191 {
192 // Trying all shifting possibilities of the site depending on the domain periodicity.
193 for ( auto const & periodicity : periodicShift )
194 {
195 auto currSite = site;
196
197 // Shifting site.
198 for ( std::size_t dim = 0; dim < Voro::Space::dimension ; ++dim )
199 if ( periodicity[dim] )
200 currSite[dim] += ( pt[dim] < currSite[dim] ? -1 : 1 ) * extent[dim];
201
202 // Checking raw-distance.
203 const auto dbis = voro.metric()->rawDistance( pt, currSite );
204 if ( dbis < dist )
205 {
206 trace.error() << "DT Error at " << pt
207 << " Voro:" << psite << " (" << dist << ")"
208 << " from set:" << site << "(" << dbis << ")"
209 << " projected from " << currSite << "." << std::endl;
210 return false;
211 }
212 }
213 }
214 }
215
216 return true;
217}
218
219
221{
224
225 return true;
226}
227
232bool testVoronoiMap( std::array<bool, 2> const& periodicity = { {false, false} } )
233{
234 unsigned int nbok = 0;
235 unsigned int nb = 0;
236
237 trace.beginBlock ( "Testing VoronoiMap2D with periodicity " + formatPeriodicity(periodicity) );
238
239 Z2i::Point a(-10,-10);
240 Z2i::Point b(10,10);
241 Z2i::Domain domain(a,b);
242
243 Z2i::DigitalSet mySet(domain);
244
245 for(Z2i::Domain::ConstIterator it = domain.begin(), itend = domain.end();
246 it != itend;
247 ++it)
248 mySet.insertNew( *it );
249
250
251 Z2i::DigitalSet sites(domain);
252
253 sites.insertNew( Z2i::Point(3,-6));
254 sites.insertNew( Z2i::Point(9,0));
255 sites.insertNew( Z2i::Point(-2,0));
256
257 for(Z2i::DigitalSet::ConstIterator it = sites.begin(), itend = sites.end();
258 it != itend; ++it)
259 mySet.erase (*it);
260
261
262
265 L2Metric l2;
266 Voro2 voro(&domain, &mySet, &l2, periodicity);
267
268 for(int j=-10; j <= 10; j++)
269 {
270 for(int i=-10; i<=10; i++)
271 trace.info() << "("<<voro( Z2i::Point(i,j))[0]<<","<< voro( Z2i::Point(i,j))[1]<<") ";
272 trace.info()<<std::endl;
273 }
274
275 trace.info()<<"Exporting to SVG"<<std::endl;
276
277 Board2D board;
278 for(Voro2::OutputImage::Domain::ConstIterator it = voro.domain().begin(),
279 itend = voro.domain().end();
280 it != itend; ++it)
281 {
282 const auto p = calcPointModuloDomain( voro(*it), voro.domain() );
283 const unsigned char c = ( p[1]*13 + p[0] * 7) % 256;
284 board << CustomStyle( (*it).className(), new CustomColors(Color(c,c,c),Color(c,c,c)))
285 << (*it);
286 }
287
288 board.saveSVG( ( "Voromap." + formatPeriodicity(periodicity) + ".svg" ).c_str() );
289
290 trace.beginBlock("Validating the Voronoi Map");
291 nbok += checkVoronoi(sites,voro) ? 1 : 0;
292 nb++;
293 trace.endBlock();
294
295 trace.endBlock();
296
297 return nbok == nb;
298}
299
300
301
306template<typename Set>
307bool testVoronoiMapFromSites2D( const Set &aSet,
308 const std::string &name,
309 std::array<bool, 2> const& periodicity = { {false, false} } )
310{
311 unsigned int nbok = 0;
312 unsigned int nb = 0;
313
314 Set mySet(aSet.domain());
315
316 for(typename Set::Domain::ConstIterator it = aSet.domain().begin(), itend = aSet.domain().end();
317 it != itend;
318 ++it)
319 mySet.insertNew( *it );
320
321
322 for(typename Set::ConstIterator it = aSet.begin(), itend = aSet.end();
323 it != itend; ++it)
324 mySet.erase (*it);
325
326
327 trace.beginBlock(" Voro computation");
330 L2Metric l2;
331 Voro2 voro(aSet.domain(), mySet, l2, periodicity );
332
333 trace.endBlock();
334
335 trace.beginBlock(" Voronoi computation l_6");
337 L6Metric l6;
339 Voro6 voro6( aSet.domain(), mySet, l6, periodicity );
340 trace.endBlock();
341
342
343
344 trace.beginBlock(" DT computation");
346 DT dt( aSet.domain(), mySet, l2, periodicity );
347 trace.endBlock();
348
349
350 if ( (aSet.domain().upperBound()[1] - aSet.domain().lowerBound()[1]) <20)
351 {
352 for(int j= aSet.domain().lowerBound()[1]; j <= aSet.domain().upperBound()[1]; j++)
353 {
354 for(int i=aSet.domain().lowerBound()[0]; i<=aSet.domain().upperBound()[0]; i++)
355 if ( aSet.find( Z2i::Point(i,j) ) != aSet.end() )
356 std::cout <<"X ";
357 else
358 std::cout<<"0 ";
359 trace.info()<<std::endl;
360 }
361
362 trace.info() << std::endl;
363
364 for(int j= aSet.domain().lowerBound()[1]; j <= aSet.domain().upperBound()[1]; j++)
365 {
366 for(int i=aSet.domain().lowerBound()[0]; i<=aSet.domain().upperBound()[0]; i++)
367 trace.info() << "("<<voro( Z2i::Point(i,j))[0]<<","<< voro( Z2i::Point(i,j))[1]<<") ";
368 trace.info()<<std::endl;
369 }
370 }
371
372 Board2D board;
373 board << voro.domain();
374 for(typename Voro2::OutputImage::Domain::ConstIterator it = voro.domain().begin(), itend = voro.domain().end();
375 it != itend; ++it)
376 {
377 if (!mySet(*it))
378 board << (*it);
379 }
380 std::string filename = "Voromap-" + name + "-orig." + formatPeriodicity(periodicity) + ".svg";
381 board.saveSVG(filename.c_str());
382
383 board.clear();
384 board << voro.domain();
385 board.setPenColor(Color(0,0,0));
386 for(typename Voro2::OutputImage::Domain::ConstIterator it = voro.domain().begin(), itend = voro.domain().end();
387 it != itend; ++it)
388 {
389 Z2i::Point p = voro(*it);
390 if ((p != (*it)) && (p != voro.domain().upperBound() + Z2i::Point::diagonal(1))
391 && (p != voro.domain().lowerBound()))
392 Display2DFactory::draw( board, p - (*it), (*it));
393 }
394
395 filename = "Voromap-" + name + "-diag." + formatPeriodicity(periodicity) + ".svg";
396 board.saveSVG(filename.c_str());
397
398
399 board.clear();
400 board << voro.domain();
401 for(typename Voro2::OutputImage::Domain::ConstIterator it = voro.domain().begin(), itend = voro.domain().end();
402 it != itend; ++it)
403 {
404 const auto p = calcPointModuloDomain( voro(*it), voro.domain() );
405 unsigned char c = (p[1]*13 + p[0] * 7) % 256;
406 board << CustomStyle( (*it).className(), new CustomColors(Color(c,c,c),Color(c,c,c)))
407 << (*it);
408 }
409
410 filename = "Voromap-" + name + "." + formatPeriodicity(periodicity) + ".svg";
411 board.saveSVG(filename.c_str());
412
413 filename = "Voromap-" + name + "-hue." + formatPeriodicity(periodicity) + ".svg";
414 saveVoroMap(filename.c_str(),voro,2);
415
416
417 board.clear();
418 for(typename Voro6::OutputImage::Domain::ConstIterator it = voro6.domain().begin(),
419 itend = voro6.domain().end();
420 it != itend; ++it)
421 {
422 Z2i::Point p = voro6(*it);
423 if (p != (*it))
424 Display2DFactory::draw( board, p - (*it), (*it));
425 }
426
427 filename = "Voromap-diag-l6-" + name + "." + formatPeriodicity(periodicity) + ".svg";
428 board.saveSVG(filename.c_str());
429
430 board.clear();
431 for(typename Voro6::OutputImage::Domain::ConstIterator it = voro6.domain().begin(), itend = voro6.domain().end();
432 it != itend; ++it)
433 {
434 const auto p = calcPointModuloDomain( voro6(*it), voro6.domain() );
435 unsigned char c = (p[1]*13 + p[0] * 7) % 256;
436 board << CustomStyle( (*it).className(), new CustomColors(Color(c,c,c),Color(c,c,c)))
437 << (*it);;
438 }
439
440 filename = "Voromap-l6" + name + "." + formatPeriodicity(periodicity) + ".svg";
441 board.saveSVG(filename.c_str());
442
443 filename = "Voromap-hue-l6-" + name + "." + formatPeriodicity(periodicity) + ".svg";
444 saveVoroMap(filename.c_str(),voro6,3);
445
446
447 trace.beginBlock("Validating the Voronoi Map");
448 nbok += checkVoronoi(aSet, voro) ? 1 : 0;
449 nb++;
450 trace.endBlock();
451
452 trace.beginBlock("Validating the Voronoi Map l_6");
453 nbok += checkVoronoi(aSet, voro6) ? 1 : 0;
454 nb++;
455 trace.endBlock();
456
457 return nbok == nb;
458}
459
464template<typename Set>
465bool testVoronoiMapFromSites( const Set &aSet )
466{
467 std::array<bool, Set::dimension> periodicity;
468 periodicity.fill( false );
469 return testVoronoiMapFromSites( aSet, periodicity );
470}
471
472template<typename Set>
473bool testVoronoiMapFromSites( const Set &aSet, std::array<bool, Set::Space::dimension> const & periodicity )
474{
475 unsigned int nbok = 0;
476 unsigned int nb = 0;
477
478 Set mySet(aSet.domain());
479
480 for(typename Set::Domain::ConstIterator it = aSet.domain().begin(),
481 itend = aSet.domain().end();
482 it != itend;
483 ++it)
484 mySet.insertNew( *it );
485
486
487 for(typename Set::ConstIterator it = aSet.begin(), itend = aSet.end();
488 it != itend; ++it)
489 mySet.erase (*it);
490
491
492
493 trace.beginBlock(" Voronoi computation l_2");
496 L2Metric l2;
497 Voro2 voro(aSet.domain(), mySet, l2, periodicity);
498 trace.endBlock();
499
500 trace.beginBlock("Validating the Voronoi Map");
501 nbok += checkVoronoi(aSet, voro) ? 1 : 0;
502 nb++;
503 trace.endBlock();
504
505 trace.beginBlock(" Voronoi computation l_3");
508 L3Metric l3;
509 Voro3 voro3(aSet.domain(), mySet, l3, periodicity);
510 trace.endBlock();
511
512 trace.beginBlock("Validating the Voronoi Map l_3");
513 nbok += checkVoronoi(aSet, voro3) ? 1 : 0;
514 nb++;
515 trace.endBlock();
516
517 trace.beginBlock(" DT computation");
519 DT dt(aSet.domain(), mySet, l2, periodicity);
520 trace.endBlock();
521
522 return nbok == nb;
523}
524
525
527{
528
529 Z2i::Point a(-10,-10);
530 Z2i::Point b(10,10);
531 Z2i::Domain domain(a,b);
532
533 Z2i::DigitalSet sites(domain);
534 bool ok = true;
535
536 sites.insertNew( Z2i::Point(3,-6));
537 sites.insertNew( Z2i::Point(9,0));
538 sites.insertNew( Z2i::Point(-3,0));
539
540 for ( std::size_t i = 0; i < 4; ++i )
541 {
542 auto const periodicity = getPeriodicityFromInteger<2>(i);
543 trace.beginBlock( "Simple2D with periodicity " + formatPeriodicity(periodicity) );
544 ok = ok && testVoronoiMapFromSites2D( sites, "simple", periodicity );
545 trace.endBlock();
546 }
547
548 return ok;
549
550}
551
553{
554 unsigned size=16;
555 Z2i::Point a(0,0);
556 Z2i::Point b(size,size);
557 Z2i::Domain domain(a,b);
558
559 Z2i::DigitalSet sites(domain);
560 bool ok = true;
561
562 for(unsigned int i = 0 ; i < size; ++i)
563 {
564 Z2i::Point p( rand() % (b[0]) - a[0], rand() % (b[1]) + a[1] );
565 sites.insert( p );
566 }
567
568 for ( std::size_t i = 0; i < 4; ++i )
569 {
570 auto const periodicity = getPeriodicityFromInteger<2>(i);
571 trace.beginBlock( "Random 2D with periodicity " + formatPeriodicity(periodicity) );
572 ok = ok && testVoronoiMapFromSites2D( sites, "random", periodicity );
573 trace.endBlock();
574 }
575
576 for(unsigned int i = 0 ; i < size*size-size; ++i)
577 {
578 Z2i::Point p( rand() % (b[0]) - a[0], rand() % (b[1]) + a[1] );
579 sites.insert( p );
580 }
581
582 for ( std::size_t i = 0; i < 4; ++i )
583 {
584 auto const periodicity = getPeriodicityFromInteger<2>(i);
585 trace.beginBlock( "Random 2D (dense) with periodicity " + formatPeriodicity(periodicity) );
586 ok = ok && testVoronoiMapFromSites2D( sites, "random-dense", periodicity );
587 trace.endBlock();
588 }
589
590 return ok;
591}
592
593
595{
596
597 Z3i::Point a(-10,-10,-10);
598 Z3i::Point b(10,10,10);
599 Z3i::Domain domain(a,b);
600
601 Z3i::DigitalSet sites(domain);
602 bool ok = true;
603
604 sites.insertNew( Z3i::Point(0,0,-6));
605 sites.insertNew( Z3i::Point(6,0,0));
606 sites.insertNew( Z3i::Point(-6,0,3));
607
608 for ( auto i = 0; i < 8; ++i )
609 {
610 auto const periodicity = getPeriodicityFromInteger<3>(i);
611 trace.beginBlock( "Simple3D with periodicity " + formatPeriodicity(periodicity) );
612 ok = ok && testVoronoiMapFromSites( sites, periodicity );
613 trace.endBlock();
614 }
615
616 return ok;
617}
618
620{
621 std::size_t const N = 16;
622
623 Z3i::Point a(0, 0, 0);
624 Z3i::Point b(N, N, N);
625 Z3i::Domain domain(a,b);
626
627 Z3i::DigitalSet sites(domain);
628 bool ok = true;
629
630 for(unsigned int i = 0 ; i < N; ++i)
631 {
632 Z3i::Point p( rand() % (b[0]) - a[0],
633 rand() % (b[1]) + a[1],
634 rand() % (b[2]) + a[2] );
635 sites.insert( p );
636 }
637
638 for ( std::size_t i = 0; i < 8; ++i )
639 {
640 auto const periodicity = getPeriodicityFromInteger<3>(i);
641 trace.beginBlock( "Random 3D with periodicity " + formatPeriodicity(periodicity) );
642 ok = ok && testVoronoiMapFromSites( sites, periodicity );
643 trace.endBlock();
644 }
645
646 return ok;
647}
648
649
650
652{
653
654 typedef SpaceND<4> Space4;
655 Space4::Point a(0,0,0,0);
656 Space4::Point b(5,5,5,5);
658
660 bool ok = true;
661
662 sites.insertNew( Space4::Point(1,4,1,1));
663 sites.insertNew( Space4::Point(3,1,3,1));
664 sites.insertNew( Space4::Point(0,0,0,0));
665
666 for ( std::size_t i = 0; i < 16; ++i )
667 {
668 auto const periodicity = getPeriodicityFromInteger<4>(i);
669 trace.beginBlock( "Simple4D with periodicity " + formatPeriodicity(periodicity) );
670 ok = ok && testVoronoiMapFromSites( sites, periodicity );
671 trace.endBlock();
672 }
673
674 return ok;
675}
676
677
679// Standard services - public :
680
681int main( int argc, char** argv )
682{
683 trace.beginBlock ( "Testing class VoronoiMap" );
684 trace.info() << "Args:";
685 for ( int i = 0; i < argc; ++i )
686 trace.info() << " " << argv[ i ];
687 trace.info() << endl;
688
689 bool res = testCheckConcept()
690 && testVoronoiMap()
691 && testVoronoiMap( { {true, false} } )
692 && testVoronoiMap( { {false, true} } )
693 && testVoronoiMap( { {true, true} } )
694 && testSimple2D()
696 && testSimple3D()
698 && testSimple4D()
699 ; // && ... other tests
700
701 trace.emphase() << ( res ? "Passed." : "Error." ) << endl;
702 trace.endBlock();
703 return res ? 0 : 1;
704}
705// //
Aim: This class specializes a 'Board' class so as to display DGtal objects more naturally (with <<)....
Definition: Board2D.h:71
Structure representing an RGB triple with alpha component.
Definition: Color.h:68
Aim: A wrapper class around a STL associative container for storing sets of digital points within som...
Container::const_iterator ConstIterator
ConstIterator type of the container;.
Aim: A container class for storing sets of digital points within some given domain.
void insertNew(const Point &p)
Aim: Implementation of the linear in time distance transformation for separable metrics.
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 a colo...
Iterator for HyperRectDomain.
const ConstIterator & begin() const
const Point & lowerBound() const
const Point & upperBound() const
const ConstIterator & end() const
static Self diagonal(Component val=1)
void beginBlock(const std::string &keyword="")
std::ostream & emphase()
std::ostream & error()
std::ostream & info()
double endBlock()
Aim: Implementation of the linear in time Voronoi map construction.
Definition: VoronoiMap.h:127
Board & setPenColor(const DGtal::Color &color)
Definition: Board.cpp:298
void clear(const DGtal::Color &color=DGtal::Color::None)
Definition: Board.cpp:152
void saveSVG(const char *filename, PageSize size=Board::BoundingBox, double margin=10.0) const
Definition: Board.cpp:1012
ExactPredicateLpSeparableMetric< Space, 2 > L2Metric
Definition: StdDefs.h:118
DGtal is the top-level namespace which contains all DGtal functions and types.
Trace trace
Definition: Common.h:154
STL namespace.
Custom style class redefining the pen color and the fill color. You may use Board2D::Color::None for ...
Definition: Board2D.h:279
static void draw(DGtal::Board2D &board, const DGtal::DiscreteExteriorCalculus< dimEmbedded, dimAmbient, TLinearAlgebraBackend, TInteger > &calculus)
Aim: Defines the concept describing a read-only image, which is a refinement of CPointFunctor.
Definition: CConstImage.h:95
int main()
Definition: testBits.cpp:56
MyPointD Point
Definition: testClone2.cpp:383
Domain domain
const Point aPoint(3, 4)
double mynorm(const Point &point, const double p)
bool testSimple2D()
Point calcPointModuloDomain(Point aPoint, Domain const &aDomain)
bool testSimple3D()
bool testCheckConcept()
bool testVoronoiMap(std::array< bool, 2 > const &periodicity={ {false, false} })
bool testSimpleRandom3D()
std::string formatPeriodicity(std::array< bool, N > const &aPeriodicity)
std::array< bool, N > getPeriodicityFromInteger(std::size_t anInteger)
void saveVoroMap(const std::string &filename, const VoroMap &output, const double p)
bool testVoronoiMapFromSites(const Set &aSet)
bool testSimple4D()
bool testVoronoiMapFromSites2D(const Set &aSet, const std::string &name, std::array< bool, 2 > const &periodicity={ {false, false} })
bool testSimpleRandom2D()
bool checkVoronoi(const Set &aSet, const Voro &voro)