DGtal 1.3.0
Loading...
Searching...
No Matches
testPowerMap.cpp
Go to the documentation of this file.
1
31#include <iostream>
32#include <array>
33
34#include "DGtal/base/Common.h"
35#include "DGtal/helpers/StdDefs.h"
36#include "DGtal/geometry/volumes/distance/PowerMap.h"
37#include "DGtal/geometry/volumes/distance/ExactPredicateLpPowerSeparableMetric.h"
38#include "DGtal/geometry/volumes/distance/ReverseDistanceTransformation.h"
39#include "DGtal/kernel/sets/DigitalSetDomain.h"
40#include "DGtal/images/ImageContainerBySTLMap.h"
41#include "DGtal/images/ImageContainerBySTLVector.h"
42#include "DGtal/kernel/sets/DigitalSetBySTLSet.h"
44
45using namespace std;
46using namespace DGtal;
47
49// Functions for testing class PowerMap.
51
62template < std::size_t N >
63std::array<bool, N> getPeriodicityFromInteger( std::size_t anInteger )
64{
65 std::array<bool, N> periodicity;
66 for ( std::size_t i = 0, mask = 1 ; i < N ; ++i, mask *= 2 )
67 periodicity[i] = anInteger & mask;
68
69 return periodicity;
70}
71
78template < std::size_t N >
79std::string formatPeriodicity( std::array<bool, N> const& aPeriodicity )
80{
81 std::string str;
82 for ( std::size_t i = 0; i < N; ++i )
83 str += aPeriodicity[i] ? '1' : '0';
84
85 return str;
86}
87
90template < typename PowerMap>
91bool checkPowerMap( const PowerMap & aPowerMap)
92{
93 using Point = typename PowerMap::Point;
94
95 // Domain extent
96 Point const extent = aPowerMap.domain().upperBound() - aPowerMap.domain().lowerBound() + Point::diagonal();
97
98 // Site shifting depending on the periodicity.
99 std::vector< typename PowerMap::PeriodicitySpec > periodicShift;
100 for ( std::size_t i = 0; i < ( 1u << PowerMap::Space::dimension ); ++i )
101 {
102 const auto periodicity = getPeriodicityFromInteger< PowerMap::Space::dimension >( i );
103
104 // Checking if this periodicity possibility is valid.
105 bool isValid = true;
106 for ( std::size_t j = 0; j < periodicity.size(); ++j )
107 if ( periodicity[j] && ! aPowerMap.isPeriodic(j) )
108 {
109 isValid = false;
110 break;
111 }
112
113 if ( isValid )
114 periodicShift.push_back( periodicity );
115 }
116
117 // Weight image
118 auto const& image = *(aPowerMap.weightImagePtr());
119
120 // Power metric
121 auto const& metric = *(aPowerMap.metricPtr());
122
123 // Pruning sites.
125 for ( auto const& pt : image.domain() )
126 if ( image(pt) > 0 )
127 aSet.insertNew( pt );
128
129 // Checking site for all domain points.
130 for ( auto const& pt : aPowerMap.domain() )
131 {
132 // Getting reference size.
133 const Point psite = aPowerMap(pt);
134 const Point boundedPSite = aPowerMap.projectPoint( psite );
135
136 // Checking that the reference site is actually a given site.
137 // Also testing projectPoint method.
138 if ( std::find( aSet.begin(), aSet.end(), boundedPSite ) == aSet.end() )
139 {
140 trace.error() << "The PowerMap point " << psite
141 << " projected to " << aPowerMap.projectPoint( psite )
142 << " is not a valid site." << std::endl;
143 return false;
144 }
145
146 // Calculating reference power distance.
147 const auto dist = metric.powerDistance( pt, psite, image( boundedPSite ) );
148
149 // Checking if we can found a better site.
150 for ( auto site : aSet )
151 {
152 // Trying all shifting possibilities of the site depending on the domain periodicity.
153 for ( auto const & periodicity : periodicShift )
154 {
155 auto currSite = site;
156
157 // Shifting site.
158 for ( std::size_t dim = 0; dim < PowerMap::Space::dimension ; ++dim )
159 if ( periodicity[dim] )
160 currSite[dim] += ( pt[dim] < currSite[dim] ? -1 : 1 ) * extent[dim];
161
162 // Checking raw-distance.
163 const auto dbis = metric.powerDistance( pt, currSite, image( site ) );
164 if ( dbis < dist )
165 {
166 trace.error() << "DT Error at " << pt
167 << " PowerMap:" << psite << "@" << image(boundedPSite) << " (" << dist << ")"
168 << " from set:" << site << "@" << image(site) << " (" << dbis << ")"
169 << " projected from " << currSite << "." << std::endl;
170 return false;
171 }
172 }
173 }
174 }
175
176 return true;
177}
178
183bool testPowerMap( std::array<bool, 2> const& aPeriodicity = {{ false, false }} )
184{
185 unsigned int nbok = 0;
186 unsigned int nb = 0;
187
188 trace.beginBlock ( "Testing PowerMap2D with periodicity " + formatPeriodicity(aPeriodicity) );
189
191 Z2i::Domain domainLarge(Z2i::Point(0,0),Z2i::Point(10,10));
192
194 set.insertNew(Z2i::Point(3,3));
195 //set.insertNew(Z2i::Point(3,7));
196 set.insertNew(Z2i::Point(7,7));
197
200
201 Image image( new SetDomain( set ) );
202
203 //Setting some values
204 image.setValue(Z2i::Point(3,3), 9);
205 // image.setValue(Z2i::Point(3,7), 0);
206 image.setValue(Z2i::Point(7,7), 16);
207
208 Z2i::L2PowerMetric l2power;
209 PowerMap<Image, Z2i::L2PowerMetric> power(&domainLarge, &image, &l2power, aPeriodicity);
210 for(unsigned int i=0; i<11; i++)
211 {
212 for(unsigned int j=0; j<11; j++)
213 if (image.domain().isInside(Z2i::Point(i,j)))
214 trace.info()<< image(Z2i::Point(i,j))<<" ";
215 else
216 trace.info()<< "0 ";
217 trace.info()<<std::endl;
218 }
219 trace.info()<<std::endl;
220
221 //Power Map
222 for(unsigned int i=0; i<11; i++)
223 {
224 for(unsigned int j=0; j<11; j++)
225 trace.info()<< power(Z2i::Point(i,j))[0]<<","<<power(Z2i::Point(i,j))[1]<<" ";
226 trace.info()<<std::endl;
227 }
228 trace.info() << "REDT" << std::endl;
229 trace.info() << std::endl;
230
231 // Checking PowerMap validity.
232 trace.beginBlock("Validating the Power Map");
233 nbok += checkPowerMap( power ) ? 1 : 0;
234 nb++;
235 trace.endBlock();
236 trace.info() << std::endl;
237
238 //Reconstruction
239 for(unsigned int i=0; i<11; i++)
240 {
241 for(unsigned int j=0; j<11; j++)
242 {
243 Z2i::Point p(i,j);
244 DGtal::int64_t dist = (i-power(p)[0])*(i-power(p)[0]) +
245 ( j-power(p)[1])*(j-power(p)[1]) - image( power.projectPoint(power(p)) );
246 trace.info()<< dist;
247 }
248 std::cerr<<std::endl;
249 }
250 trace.info()<<std::endl;
251
252 //Reconstruction
253 for(unsigned int i=0; i<11; i++)
254 {
255 for(unsigned int j=0; j<11; j++)
256 {
257 Z2i::Point p(i,j);
258 DGtal::int32_t dist = (i-power(p)[0])*(i-power(p)[0]) +
259 ( j-power(p)[1])*(j-power(p)[1]) - image( power.projectPoint(power(p)) );
260 if (dist>=0)
261 std::cerr<< "0 ";
262 else
263 std::cerr<< "X ";
264 }
265 std::cerr<<std::endl;
266 }
267
268 nbok ++;
269 nb++;
270 trace.info() << "(" << nbok << "/" << nb << ") "
271 << "true == true" << std::endl;
272 trace.endBlock();
273
274 return nbok == nb;
275}
276
278{
281 using PowerMetric = Z2i::L2PowerMetric;
282
283 BOOST_CONCEPT_ASSERT(( concepts::CConstImage< PowerMap< Image, PowerMetric > >));
284
285 return true;
286}
287
288template < typename Set >
289bool testPowerMapFromSites( const Set &aSet )
290{
291 std::array<bool, Set::dimension> periodicity;
292 periodicity.fill( false );
293 return testPowerMapFromSites( aSet, periodicity );
294}
295
296template < typename Set >
297bool testPowerMapFromSites( const Set &aSet, std::array<bool, Set::Space::dimension> const & periodicity )
298{
299 unsigned int nbok = 0;
300 unsigned int nb = 0;
301
302 using SetDomain = DigitalSetDomain< Set >;
304
305 // Weight image.
306 Image image( new SetDomain( aSet ) );
307
308 // Setting weights
309 typename Image::Value weight = 2;
310 for ( auto const & pt : image.domain() )
311 {
312 image.setValue( pt, weight );
313 weight *= 2;
314 }
315
316 trace.beginBlock(" Power Map computation l_2");
319 L2PowerMetric l2;
320 Power2 power2( aSet.domain(), image, l2, periodicity);
321 trace.endBlock();
322
323 trace.beginBlock("Validating the Power Map");
324 nbok += checkPowerMap( power2 ) ? 1 : 0;
325 nb++;
326 trace.endBlock();
327
328 trace.beginBlock(" Power Map computation l_3");
331 L3PowerMetric l3;
332 Power3 power3( aSet.domain(), image, l3, periodicity);
333 trace.endBlock();
334
335 trace.beginBlock("Validating the Power Map l_3");
336 nbok += checkPowerMap( power3 ) ? 1 : 0;
337 nb++;
338 trace.endBlock();
339
340 trace.beginBlock("Reverse DT computation");
342 RDT rdt(aSet.domain(), power2.weightImagePtr(), l2, periodicity);
343 trace.endBlock();
344
345 return nbok == nb;
346}
347
349{
350
351 Z2i::Point a(-10,-10);
352 Z2i::Point b(10,10);
353 Z2i::Domain domain(a,b);
354
356 bool ok = true;
357
358 sites.insertNew( Z2i::Point(3,-6));
359 sites.insertNew( Z2i::Point(9,0));
360 sites.insertNew( Z2i::Point(-3,0));
361
362 for ( std::size_t i = 0; i < 4; ++i )
363 {
364 auto const periodicity = getPeriodicityFromInteger<2>(i);
365 trace.beginBlock( "Simple2D with periodicity " + formatPeriodicity(periodicity) );
366 ok = ok && testPowerMapFromSites( sites, periodicity );
367 trace.endBlock();
368 }
369
370 return ok;
371
372}
373
375{
376
377 Z3i::Point a(-10,-10,-10);
378 Z3i::Point b(10,10,10);
379 Z3i::Domain domain(a,b);
380
382 bool ok = true;
383
384 sites.insertNew( Z3i::Point(0,0,-6));
385 sites.insertNew( Z3i::Point(6,0,0));
386 sites.insertNew( Z3i::Point(-6,0,3));
387
388 for ( std::size_t i = 0; i < 8; ++i )
389 {
390 auto const periodicity = getPeriodicityFromInteger<3>(i);
391 trace.beginBlock( "Simple3D with periodicity " + formatPeriodicity(periodicity) );
392 ok = ok && testPowerMapFromSites( sites, periodicity );
393 trace.endBlock();
394 }
395
396 return ok;
397}
398
400{
401
402 typedef SpaceND<4> Space4;
403 Space4::Point a(0,0,0,0);
404 Space4::Point b(5,5,5,5);
406
408 bool ok = true;
409
410 sites.insertNew( Space4::Point(1,4,1,1));
411 sites.insertNew( Space4::Point(3,1,3,1));
412 sites.insertNew( Space4::Point(0,0,0,0));
413
414 for ( std::size_t i = 0; i < 16; ++i )
415 {
416 auto const periodicity = getPeriodicityFromInteger<4>(i);
417 trace.beginBlock( "Simple4D with periodicity " + formatPeriodicity(periodicity) );
418 ok = ok && testPowerMapFromSites( sites, periodicity );
419 trace.endBlock();
420 }
421
422 return ok;
423}
424
426// Standard services - public :
427
428int main( int argc, char** argv )
429{
430 trace.beginBlock ( "Testing class PowerMap" );
431 trace.info() << "Args:";
432 for ( int i = 0; i < argc; ++i )
433 trace.info() << " " << argv[ i ];
434 trace.info() << endl;
435
436 bool res = testCheckConcept()
437 && testPowerMap()
438 && testPowerMap( {{ true, false }} )
439 && testPowerMap( {{ false, true }} )
440 && testPowerMap( {{ true, true }} )
441 && testSimple2D()
442 && testSimple3D()
443 && testSimple4D()
444 ; // && ... other tests
445
446 trace.emphase() << ( res ? "Passed." : "Error." ) << endl;
447 trace.endBlock();
448 return res ? 0 : 1;
449}
450// //
Aim: A container class for storing sets of digital points within some given domain.
ConstIterator end() const
void insertNew(const Point &p)
ConstIterator begin() const
Aim: Constructs a domain limited to the given digital set.
Aim: implements weighted separable l_p metrics with exact predicates.
Aim: implements association bewteen points lying in a digital domain and values.
Definition: Image.h:70
Aim: Implementation of the linear in time Power map construction.
Definition: PowerMap.h:111
const WeightImage * weightImagePtr() const
Definition: PowerMap.h:275
Point projectPoint(Point aPoint) const
bool isPeriodic(const Dimension n) const
Definition: PowerMap.h:295
const Domain & domain() const
Definition: PowerMap.h:239
const PowerSeparableMetric * metricPtr() const
Definition: PowerMap.h:267
Space::Point Point
Definition: PowerMap.h:124
Aim: Implementation of the linear in time reverse distance transformation for separable metrics.
void beginBlock(const std::string &keyword="")
std::ostream & emphase()
std::ostream & error()
std::ostream & info()
double endBlock()
ExactPredicateLpPowerSeparableMetric< Space, 2 > L2PowerMetric
Definition: StdDefs.h:120
T power(const T &aVal, const unsigned int exponent)
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
Trace trace
Definition: Common.h:154
boost::int32_t int32_t
signed 32-bit integer.
Definition: BasicTypes.h:72
STL namespace.
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
bool testSimple2D()
bool testSimple3D()
bool testCheckConcept()
std::string formatPeriodicity(std::array< bool, N > const &aPeriodicity)
std::array< bool, N > getPeriodicityFromInteger(std::size_t anInteger)
bool checkPowerMap(const PowerMap &aPowerMap)
bool testSimple4D()
bool testPowerMapFromSites(const Set &aSet)
bool testPowerMap(std::array< bool, 2 > const &aPeriodicity={{ false, false }})
Domain domain