DGtal  1.2.0
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 
45 using namespace std;
46 using namespace DGtal;
47 
49 // Functions for testing class PowerMap.
51 
62 template < std::size_t N >
63 std::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 
78 template < std::size_t N >
79 std::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 
90 template < typename PowerMap>
91 bool 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 
183 bool 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 
288 template < typename Set >
289 bool testPowerMapFromSites( const Set &aSet )
290 {
291  std::array<bool, Set::dimension> periodicity;
292  periodicity.fill( false );
293  return testPowerMapFromSites( aSet, periodicity );
294 }
295 
296 template < typename Set >
297 bool 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");
318  typedef PowerMap< Image, L2PowerMetric > Power2;
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");
330  typedef PowerMap< Image, L3PowerMetric > Power3;
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 
428 int 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:271
Point projectPoint(Point aPoint) const
bool isPeriodic(const Dimension n) const
Definition: PowerMap.h:291
const PowerSeparableMetric * metricPtr() const
Definition: PowerMap.h:263
const Domain & domain() const
Definition: PowerMap.h:235
Aim: Implementation of the linear in time reverse distance transformation for separable metrics.
std::ostream & error()
void beginBlock(const std::string &keyword="")
std::ostream & emphase()
std::ostream & info()
double endBlock()
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
Aim: Defines the concept describing a read-only image, which is a refinement of CPointFunctor.
Definition: CConstImage.h:95
MyPointD Point
Definition: testClone2.cpp:383
bool testSimple2D()
int main(int argc, char **argv)
bool testSimple3D()
bool testCheckConcept()
std::string formatPeriodicity(std::array< bool, N > const &aPeriodicity)
bool checkPowerMap(const PowerMap &aPowerMap)
bool testSimple4D()
std::array< bool, N > getPeriodicityFromInteger(std::size_t anInteger)
bool testPowerMapFromSites(const Set &aSet)
bool testPowerMap(std::array< bool, 2 > const &aPeriodicity={{ false, false }})
Domain domain
Image image(domain)
unsigned int dim(const Vector &z)