DGtal  0.9.3
testChordNaivePlaneComputer.cpp
Go to the documentation of this file.
1 
30 #include <cstdlib>
32 #include <iostream>
33 #include "DGtal/base/Common.h"
34 #include "DGtal/helpers/StdDefs.h"
35 #include "DGtal/kernel/CPointPredicate.h"
36 #include "DGtal/geometry/surfaces/CAdditivePrimitiveComputer.h"
37 #include "DGtal/geometry/surfaces/ChordNaivePlaneComputer.h"
38 #include "DGtal/geometry/surfaces/ChordGenericNaivePlaneComputer.h"
40 
41 using namespace std;
42 using namespace DGtal;
43 using namespace DGtal::concepts;
44 
46 // Functions for testing class ChordNaivePlaneComputer.
48 
49 template <typename Integer>
50 Integer getRandomInteger( Integer first, Integer after_last )
51 {
52  Integer r = (Integer) rand();
53  return ( r % (after_last - first) ) + first;
54 }
55 
59 template <typename Integer, typename NaivePlaneComputer>
60 bool
61 checkPlane( Integer a, Integer b, Integer c, Integer d,
62  int diameter, unsigned int nbtries )
63 {
64  typedef typename NaivePlaneComputer::Point Point;
65  typedef typename Point::Component PointInteger;
67  Integer absA = ic.abs( a );
68  Integer absB = ic.abs( b );
69  Integer absC = ic.abs( c );
70  Integer x, y, z;
71  Dimension axis;
72  if ( ( absA >= absB ) && ( absA >= absC ) )
73  axis = 0;
74  else if ( ( absB >= absA ) && ( absB >= absC ) )
75  axis = 1;
76  else
77  axis = 2;
78  Point p;
79  NaivePlaneComputer plane;
80  plane.init( axis, 1, 1 );
81  // Checks that points within the naive plane are correctly recognized.
82  unsigned int nb = 0;
83  unsigned int nbok = 0;
84  while ( nb != nbtries )
85  {
86  p[ 0 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
87  p[ 1 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
88  p[ 2 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
89  x = (Integer) p[ 0 ];
90  y = (Integer) p[ 1 ];
91  z = (Integer) p[ 2 ];
92  switch ( axis ) {
93  case 0: p[ 0 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - b * y - c * z, a ) ); break;
94  case 1: p[ 1 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - c * z, b ) ); break;
95  case 2: p[ 2 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - b * y, c ) ); break;
96  }
97  bool ok_ext = plane.isExtendable( p ); // should be ok
98  bool ok = plane.extend( p ); // should be ok
99  ++nb, nbok += ok_ext ? 1 : 0;
100  ++nb, nbok += ok ? 1 : 0;
101  if ( ! ok )
102  {
103  std::cerr << "[ERROR] p=" << p << " NOT IN plane=" << plane << std::endl;
104  for ( typename NaivePlaneComputer::ConstIterator it = plane.begin(), itE = plane.end();
105  it != itE; ++it )
106  std::cerr << " " << *it;
107  std::cerr << endl;
108  std::cerr << "d <= a*x+b*y+c*z <= d+max(a,b,c)"
109  << d << " <= " << a << "*" << p[0]
110  << "+" << b << "*" << p[1]
111  << "+" << c << "*" << p[2]
112  << " = " << (a*p[0]+b*p[1]+c*p[2])
113  << std::endl;
114  break;
115  }
116  if ( ! ok_ext )
117  {
118  std::cerr << "[ERROR] p=" << p << " was NOT extendable IN plane=" << plane << std::endl;
119  break;
120  }
121  // else
122  // std::cerr << "[OK] p=" << p << " IN plane=" << plane << std::endl;
123  }
124 
125  // Checks that points outside the naive plane are correctly recognized as outliers.
126  while ( nb != (nbtries * 11 ) / 10 )
127  {
128  p[ 0 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
129  p[ 1 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
130  p[ 2 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
131  x = (Integer) p[ 0 ];
132  y = (Integer) p[ 1 ];
133  z = (Integer) p[ 2 ];
134  switch ( axis ) {
135  case 0: p[ 0 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - b * y - c * z, a ) ); break;
136  case 1: p[ 1 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - c * z, b ) ); break;
137  case 2: p[ 2 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - b * y, c ) ); break;
138  }
139  PointInteger tmp = getRandomInteger<PointInteger>( 2, 5 )
140  * (2*getRandomInteger<PointInteger>( 0, 2 ) - 1 );
141  p[ axis ] += tmp;
142  bool ok_ext = ! plane.isExtendable( p ); // should *not* be ok
143  bool ok = ! plane.extend( p ); // should *not* be ok
144  ++nb, nbok += ok ? 1 : 0;
145  ++nb, nbok += ok_ext ? 1 : 0;
146  if ( ! ok )
147  {
148  std::cerr << "[ERROR] p=" << p << " IN plane=" << plane << std::endl;
149  break;
150  }
151  if ( ! ok_ext )
152  {
153  std::cerr << "[ERROR] p=" << p << " was extendable IN plane=" << plane << std::endl;
154  break;
155  }
156  // else
157  // std::cerr << "[OK] p=" << p << " IN plane=" << plane << std::endl;
158  }
159  return nb == nbok;
160 }
161 
165 template <typename Integer, typename NaivePlaneComputer>
166 bool
167 checkPlaneGroupExtension( Integer a, Integer b, Integer c, Integer d,
168  int diameter, unsigned int nbtries )
169 {
170  typedef typename NaivePlaneComputer::Point Point;
171  typedef typename Point::Component PointInteger;
173  Integer absA = ic.abs( a );
174  Integer absB = ic.abs( b );
175  Integer absC = ic.abs( c );
176  Integer x, y, z;
177  Dimension axis;
178  if ( ( absA >= absB ) && ( absA >= absC ) )
179  axis = 0;
180  else if ( ( absB >= absA ) && ( absB >= absC ) )
181  axis = 1;
182  else
183  axis = 2;
184  Point p;
185  NaivePlaneComputer plane;
186  plane.init( axis, 1, 1 );
187  // Checks that points within the naive plane are correctly recognized.
188  unsigned int nb = 0;
189  unsigned int nbok = 0;
190  while ( nb < nbtries )
191  {
192  std::vector<Point> points( 5 );
193  for ( unsigned int i = 0; i < 5; ++i )
194  {
195  Point & pp = points[ i ];
196  pp[ 0 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
197  pp[ 1 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
198  pp[ 2 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
199  x = (Integer) pp[ 0 ];
200  y = (Integer) pp[ 1 ];
201  z = (Integer) pp[ 2 ];
202  switch ( axis ) {
203  case 0: pp[ 0 ] = NumberTraits<Integer>::castToInt64_t
204  ( ic.ceilDiv( d - b * y - c * z, a ) ); break;
205  case 1: pp[ 1 ] = NumberTraits<Integer>::castToInt64_t
206  ( ic.ceilDiv( d - a * x - c * z, b ) ); break;
207  case 2: pp[ 2 ] = NumberTraits<Integer>::castToInt64_t
208  ( ic.ceilDiv( d - a * x - b * y, c ) ); break;
209  }
210  }
211  bool ok_ext = plane.isExtendable( points.begin(), points.end() ); // should be ok
212  bool ok = plane.extend( points.begin(), points.end() ); // should be ok
213  ++nb, nbok += ok_ext ? 1 : 0;
214  ++nb, nbok += ok ? 1 : 0;
215  if ( ! ok )
216  {
217  std::cerr << "[ERROR] p=" << points[ 0 ] << " NOT IN plane=" << plane << std::endl;
218  for ( typename NaivePlaneComputer::ConstIterator it = plane.begin(), itE = plane.end();
219  it != itE; ++it )
220  std::cerr << " " << *it;
221  std::cerr << endl;
222  std::cerr << "d <= a*x+b*y+c*z <= d+max(a,b,c)"
223  << d << " <= " << a << "*" << p[0]
224  << "+" << b << "*" << p[1]
225  << "+" << c << "*" << p[2]
226  << " = " << (a*p[0]+b*p[1]+c*p[2])
227  << std::endl;
228  break;
229  }
230  if ( ! ok_ext )
231  {
232  std::cerr << "[ERROR] p=" << p << " was NOT extendable IN plane=" << plane << std::endl;
233  break;
234  }
235  // else
236  // std::cerr << "[OK] p=" << p << " IN plane=" << plane << std::endl;
237  }
238 
239  // Checks that points outside the naive plane are correctly recognized as outliers.
240  while ( nb < (nbtries * 11 ) / 10 )
241  {
242  p[ 0 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
243  p[ 1 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
244  p[ 2 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
245  x = (Integer) p[ 0 ];
246  y = (Integer) p[ 1 ];
247  z = (Integer) p[ 2 ];
248  switch ( axis ) {
249  case 0: p[ 0 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - b * y - c * z, a ) ); break;
250  case 1: p[ 1 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - c * z, b ) ); break;
251  case 2: p[ 2 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - b * y, c ) ); break;
252  }
253  PointInteger tmp = getRandomInteger<PointInteger>( 2, 5 )
254  * (2*getRandomInteger<PointInteger>( 0, 2 ) - 1 );
255  p[ axis ] += tmp;
256  bool ok_ext = ! plane.isExtendable( p ); // should *not* be ok
257  bool ok = ! plane.extend( p ); // should *not* be ok
258  ++nb, nbok += ok ? 1 : 0;
259  ++nb, nbok += ok_ext ? 1 : 0;
260  if ( ! ok )
261  {
262  std::cerr << "[ERROR] p=" << p << " IN plane=" << plane << std::endl;
263  break;
264  }
265  if ( ! ok_ext )
266  {
267  std::cerr << "[ERROR] p=" << p << " was extendable IN plane=" << plane << std::endl;
268  break;
269  }
270  // else
271  // std::cerr << "[OK] p=" << p << " IN plane=" << plane << std::endl;
272  }
273  return nb == nbok;
274 }
275 
276 
277 
281 template <typename Integer, typename GenericNaivePlaneComputer>
282 bool
283 checkGenericPlane( Integer a, Integer b, Integer c, Integer d,
284  int diameter, unsigned int nbtries )
285 {
286  typedef typename GenericNaivePlaneComputer::Point Point;
287  typedef typename Point::Component PointInteger;
289  Integer absA = ic.abs( a );
290  Integer absB = ic.abs( b );
291  Integer absC = ic.abs( c );
292  Integer x, y, z;
293  Dimension axis;
294  if ( ( absA >= absB ) && ( absA >= absC ) )
295  axis = 0;
296  else if ( ( absB >= absA ) && ( absB >= absC ) )
297  axis = 1;
298  else
299  axis = 2;
300  Point p;
301  GenericNaivePlaneComputer plane;
302  plane.init( 1, 1 );
303  // Checks that points within the naive plane are correctly recognized.
304  unsigned int nb = 0;
305  unsigned int nbok = 0;
306  while ( nb != nbtries )
307  {
308  p[ 0 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
309  p[ 1 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
310  p[ 2 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
311  x = (Integer) p[ 0 ];
312  y = (Integer) p[ 1 ];
313  z = (Integer) p[ 2 ];
314  switch ( axis ) {
315  case 0: p[ 0 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - b * y - c * z, a ) ); break;
316  case 1: p[ 1 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - c * z, b ) ); break;
317  case 2: p[ 2 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - b * y, c ) ); break;
318  }
319  bool ok_ext = plane.isExtendable( p ); // should be ok
320  bool ok = plane.extend( p ); // should be ok
321  ++nb, nbok += ok_ext ? 1 : 0;
322  ++nb, nbok += ok ? 1 : 0;
323  if ( ! ok )
324  {
325  std::cerr << "[ERROR] p=" << p << " NOT IN plane=" << plane << std::endl;
326  for ( typename GenericNaivePlaneComputer::ConstIterator it = plane.begin(), itE = plane.end();
327  it != itE; ++it )
328  std::cerr << " " << *it;
329  std::cerr << endl;
330  std::cerr << "d <= a*x+b*y+c*z <= d+max(a,b,c)"
331  << d << " <= " << a << "*" << p[0]
332  << "+" << b << "*" << p[1]
333  << "+" << c << "*" << p[2]
334  << " = " << (a*p[0]+b*p[1]+c*p[2])
335  << std::endl;
336  break;
337  }
338  if ( ! ok_ext )
339  {
340  std::cerr << "[ERROR] p=" << p << " was NOT extendable IN plane=" << plane << std::endl;
341  break;
342  }
343  // else
344  // std::cerr << "[OK] p=" << p << " IN plane=" << plane << std::endl;
345  }
346 
347  // Checks that points outside the naive plane are correctly recognized as outliers.
348  while ( nb != (nbtries * 11 ) / 10 )
349  {
350  p[ 0 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
351  p[ 1 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
352  p[ 2 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
353  x = (Integer) p[ 0 ];
354  y = (Integer) p[ 1 ];
355  z = (Integer) p[ 2 ];
356  switch ( axis ) {
357  case 0: p[ 0 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - b * y - c * z, a ) ); break;
358  case 1: p[ 1 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - c * z, b ) ); break;
359  case 2: p[ 2 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - b * y, c ) ); break;
360  }
361  PointInteger tmp = getRandomInteger<PointInteger>( 2, 5 )
362  * (2*getRandomInteger<PointInteger>( 0, 2 ) - 1 );
363  p[ axis ] += tmp;
364  bool ok_ext = ! plane.isExtendable( p ); // should *not* be ok
365  bool ok = ! plane.extend( p ); // should *not* be ok
366  ++nb, nbok += ok ? 1 : 0;
367  ++nb, nbok += ok_ext ? 1 : 0;
368  if ( ! ok )
369  {
370  std::cerr << "[ERROR] p=" << p << " IN plane=" << plane << std::endl;
371  break;
372  }
373  if ( ! ok_ext )
374  {
375  std::cerr << "[ERROR] p=" << p << " was extendable IN plane=" << plane << std::endl;
376  break;
377  }
378  // else
379  // std::cerr << "[OK] p=" << p << " IN plane=" << plane << std::endl;
380  }
381  std::cerr << "plane = " << plane << std::endl;
382  return nb == nbok;
383 }
384 
385 
386 template <typename Integer, typename NaivePlaneComputer>
387 bool
388 checkPlanes( unsigned int nbplanes, int diameter, unsigned int nbtries )
389 {
390  //using namespace Z3i;
391  //typedef ChordNaivePlaneComputer<Z3, Integer> NaivePlaneComputer;
392  unsigned int nb = 0;
393  unsigned int nbok = 0;
394  for ( unsigned int nbp = 0; nbp < nbplanes; ++nbp )
395  {
396  Integer a = getRandomInteger<Integer>( (Integer) 0, (Integer) diameter / 2 );
397  Integer b = getRandomInteger<Integer>( (Integer) 0, (Integer) diameter / 2 );
398  Integer c = getRandomInteger<Integer>( (Integer) 0, (Integer) diameter / 2 );
399  Integer d = getRandomInteger<Integer>( (Integer) 0, (Integer) diameter / 2 );
400  if ( ( a != 0 ) || ( b != 0 ) || ( c != 0 ) )
401  {
402  ++nb, nbok += checkPlane<Integer, NaivePlaneComputer>( a, b, c, d, diameter, nbtries ) ? 1 : 0;
403  if ( nb != nbok )
404  {
405  std::cerr << "[ERROR] (Simple extension) for plane " << a << " * x + "
406  << b << " * y + " << c << " * z = " << d << std::endl;
407  break;
408  }
409  ++nb, nbok += checkPlaneGroupExtension<Integer, NaivePlaneComputer>( a, b, c, d, diameter, nbtries ) ? 1 : 0;
410  if ( nb != nbok )
411  {
412  std::cerr << "[ERROR] (Group extension) for plane " << a << " * x + "
413  << b << " * y + " << c << " * z = " << d << std::endl;
414  break;
415  }
416  }
417  }
418  return nb == nbok;
419 }
420 
424 template <typename Integer, typename NaivePlaneComputer>
425 bool
426 checkWidth( Integer a, Integer b, Integer c, Integer d,
427  int diameter, unsigned int nbtries )
428 {
429  typedef typename NaivePlaneComputer::Point Point;
430  typedef typename NaivePlaneComputer::InternalScalar InternalScalar;
432  Integer absA = ic.abs( a );
433  Integer absB = ic.abs( b );
434  Integer absC = ic.abs( c );
435  Integer x, y, z;
436  Dimension axis;
437  if ( ( absA >= absB ) && ( absA >= absC ) )
438  axis = 0;
439  else if ( ( absB >= absA ) && ( absB >= absC ) )
440  axis = 1;
441  else
442  axis = 2;
443  // Checks that points within the naive plane are correctly recognized.
444  unsigned int nb = 0;
445  unsigned int nbok = 0;
446  std::vector<Point> points( nbtries );
447  for ( unsigned int i = 0; i != nbtries; ++i )
448  {
449  Point & p = points[ i ];
450  p[ 0 ] = getRandomInteger<Integer>( -diameter+1, diameter );
451  p[ 1 ] = getRandomInteger<Integer>( -diameter+1, diameter );
452  p[ 2 ] = getRandomInteger<Integer>( -diameter+1, diameter );
453  x = (Integer) p[ 0 ];
454  y = (Integer) p[ 1 ];
455  z = (Integer) p[ 2 ];
456  switch ( axis ) {
457  case 0: p[ 0 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - b * y - c * z, a ) ); break;
458  case 1: p[ 1 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - c * z, b ) ); break;
459  case 2: p[ 2 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - b * y, c ) ); break;
460  }
461  }
462  trace.beginBlock( "Computing axis width." );
463  trace.info() << "- plane is "
464  << d << " <= " << a << "*x"
465  << "+" << b << "*y"
466  << "+" << c << "*z"
467  << " <= d + max(|a|,|b|,|c|)"
468  << std::endl;
469  trace.info() << "- " << points.size() << " points tested in diameter " << diameter
470  << std::endl;
471  double min = -1.0;
472  for ( unsigned int i = 0; i < 3; ++i )
473  {
474  std::pair<InternalScalar, InternalScalar> width
475  = NaivePlaneComputer::computeAxisWidth( i, points.begin(), points.end() );
476  double wn = NumberTraits<InternalScalar>::castToDouble( width.first );
477  double wd = NumberTraits<InternalScalar>::castToDouble( width.second );
478  trace.info() << " (" << i << ") width=" << (wn/wd) << std::endl;
479  if ( min < 0.0 ) min = wn/wd;
480  else if ( wn/wd < min ) min = wn/wd;
481  }
482  ++nb, nbok += (min < 1.0 ) ? 1 : 0;
483  trace.info() << "(" << nbok << "/" << nb << ") min width = " << min
484  << " < 1.0" << std::endl;
485  ++nb, nbok += (0.9 < min ) ? 1 : 0;
486  trace.info() << "(" << nbok << "/" << nb << ") min width = " << min
487  << " > 0.9" << std::endl;
488  trace.endBlock();
489  return nb == nbok;
490 }
491 
492 template <typename Integer, typename NaivePlaneComputer>
493 bool
494 checkWidths( unsigned int nbplanes, int diameter, unsigned int nbtries )
495 {
496  //using namespace Z3i;
497  //typedef ChordNaivePlaneComputer<Z3, Integer> NaivePlaneComputer;
498  unsigned int nb = 0;
499  unsigned int nbok = 0;
500  for ( unsigned int nbp = 0; nbp < nbplanes; ++nbp )
501  {
502  Integer a = getRandomInteger<Integer>( (Integer) 0, (Integer) diameter / 2 );
503  Integer b = getRandomInteger<Integer>( (Integer) 0, (Integer) diameter / 2 );
504  Integer c = getRandomInteger<Integer>( (Integer) 0, (Integer) diameter / 2 );
505  Integer d = getRandomInteger<Integer>( (Integer) 0, (Integer) diameter / 2 );
506  if ( ( a != 0 ) || ( b != 0 ) || ( c != 0 ) )
507  {
508  ++nb, nbok += checkWidth<Integer, NaivePlaneComputer>( a, b, c, d, diameter, nbtries ) ? 1 : 0;
509  if ( nb != nbok )
510  {
511  std::cerr << "[ERROR] (checkWidth) for plane " << a << " * x + "
512  << b << " * y + " << c << " * z = " << d << std::endl;
513  break;
514  }
515  }
516  }
517  return nb == nbok;
518 }
519 
520 
526 {
527  unsigned int nbok = 0;
528  unsigned int nb = 0;
529  typedef DGtal::int64_t Integer;
530  typedef DGtal::Z3i::Z3 Space;
531  typedef DGtal::Z3i::Point Point;
533  typedef ChordGenericNaivePlaneComputer<Space, Point, Integer> GenericNaivePlaneComputer;
534 
535  BOOST_CONCEPT_ASSERT(( CAdditivePrimitiveComputer< NaivePlaneComputer > ));
537  BOOST_CONCEPT_ASSERT(( boost::ForwardContainer< NaivePlaneComputer > ));
539  BOOST_CONCEPT_ASSERT(( CPointPredicate< NaivePlaneComputer::Primitive > ));
541 
542  trace.beginBlock ( "Testing block: ChordNaivePlaneComputer instantiation." );
543  NaivePlaneComputer plane;
544  Point pt0( 0, 0, 0 );
545  plane.init( 2, 1, 1 );
546  bool pt0_inside = plane.extend( pt0 );
547  ++nb, nbok += pt0_inside == true ? 1 : 0;
548  trace.info() << "(" << nbok << "/" << nb << ") Plane=" << plane
549  << std::endl;
550  Point pt1( Point( 2, 0, 0 ) );
551  bool pt1_inside = plane.extend( pt1 );
552  ++nb, nbok += pt1_inside == true ? 1 : 0;
553  trace.info() << "(" << nbok << "/" << nb << ") add " << pt1
554  << " Plane=" << plane << std::endl;
555  Point pt2( Point( 0, 2, 2 ) );
556  bool pt2_inside = plane.extend( pt2 );
557  ++nb, nbok += pt2_inside == true ? 1 : 0;
558  trace.info() << "(" << nbok << "/" << nb << ") add " << pt2
559  << " Plane=" << plane << std::endl;
560 
561  Point pt3( Point( 1, 1, 1 ) );
562  bool pt3_inside = plane.extend( pt3 );
563  ++nb, nbok += pt3_inside == true ? 1 : 0;
564  trace.info() << "(" << nbok << "/" << nb << ") add " << pt3
565  << " Plane=" << plane << std::endl;
566 
567  Point pt4( Point( -10, -10, 10 ) );
568  bool pt4_inside = plane.extend( pt4 );
569  ++nb, nbok += pt4_inside == false ? 1 : 0;
570  trace.info() << "(" << nbok << "/" << nb << ") impossible add " << pt4
571  << " Plane=" << plane << std::endl;
572 
573  Point pt5 = pt2 + Point( 1, 0, 1 );
574  bool pt5_inside = plane.extend( pt5 );
575  ++nb, nbok += pt5_inside == true ? 1 : 0;
576  trace.info() << "(" << nbok << "/" << nb << ") add " << pt5
577  << " Plane=" << plane << std::endl;
578 
579  Point pt6 = pt5 + Point( 6, 0, 2 );
580  bool pt6_inside = plane.extend( pt6 );
581  ++nb, nbok += pt6_inside == true ? 1 : 0;
582  trace.info() << "(" << nbok << "/" << nb << ") add " << pt6
583  << " Plane=" << plane << std::endl;
584 
585  NaivePlaneComputer plane2;
586  plane2.init( 2, 1, 1 );
587  plane2.extend( Point( 10, 0, 0 ) );
588  plane2.extend( Point( 0, 8, 0 ) );
589  plane2.extend( Point( 0, 0, 6 ) );
590  trace.info() << "(" << nbok << "/" << nb << ") "
591  << " Plane2=" << plane2 << std::endl;
592 
593  ++nb, nbok += checkPlane<Integer,NaivePlaneComputer>( 11, 5, 19, 20, 100, 100 ) ? 1 : 0;
594  trace.info() << "(" << nbok << "/" << nb
595  << ") checkPlane<Integer,NaivePlaneComputer>( 11, 5, 19, 20, 100, 100 )"
596  << std::endl;
597 
598  ++nb, nbok += checkGenericPlane<Integer,GenericNaivePlaneComputer>( 11, 5, 19, 20, 100, 100 ) ? 1 : 0;
599  trace.info() << "(" << nbok << "/" << nb
600  << ") checkGenericPlane<Integer,GenericNaivePlaneComputer>( 11, 5, 19, 20, 100, 100 )"
601  << std::endl;
602  ++nb, nbok += checkGenericPlane<Integer,GenericNaivePlaneComputer>( 17, 33, 7, 10, 100, 100 ) ? 1 : 0;
603  trace.info() << "(" << nbok << "/" << nb
604  << ") checkGenericPlane<Integer,GenericNaivePlaneComputer>( 17, 33, 7, 10, 100, 100 )"
605  << std::endl;
606  ++nb, nbok += checkPlane<Integer,NaivePlaneComputer>( 15, 8, 13, 15, 100, 100 ) ? 1 : 0;
607  trace.info() << "(" << nbok << "/" << nb
608  << ") checkPlane<Integer,NaivePlaneComputer>( 15, 8, 13, 15, 100, 100 )"
609  << std::endl;
610  ++nb, nbok += checkGenericPlane<Integer,GenericNaivePlaneComputer>( 15, 8, 13, 15, 100, 100 ) ? 1 : 0;
611  trace.info() << "(" << nbok << "/" << nb
612  << ") checkGenericPlane<Integer,GenericNaivePlaneComputer>( 15, 8, 13, 15, 100, 100 )"
613  << std::endl;
614  trace.endBlock();
615 
616  {
617  trace.beginBlock ( "Testing block: ChordNaivePlaneComputer vertical instantiation." );
618  NaivePlaneComputer ppplane;
619  Point pppt0( 0, 0, 0 );
620  ppplane.init( 2, 5, 2 );
621  bool pppt0_inside = ppplane.extend( pppt0 );
622  ++nb, nbok += pppt0_inside == true ? 1 : 0;
623  trace.info() << "(" << nbok << "/" << nb << ") Plane=" << ppplane
624  << std::endl;
625  Point pppt1( 3, 2, 2 );
626  bool pppt1_inside = ppplane.extend( pppt1 );
627  ++nb, nbok += pppt1_inside == true ? 1 : 0;
628  trace.info() << "(" << nbok << "/" << nb << ") Plane=" << ppplane
629  << std::endl;
630  Point pppt2( 0, 0, 1 );
631  bool pppt2_inside = ppplane.extend( pppt2 );
632  ++nb, nbok += pppt2_inside == true ? 1 : 0;
633  trace.info() << "(" << nbok << "/" << nb << ") Plane=" << ppplane
634  << std::endl;
635  Point pppt3 = pppt1 + Point( 0, 0, 1 );
636  bool pppt3_inside = ppplane.extend( pppt3 );
637  ++nb, nbok += pppt3_inside == true ? 1 : 0;
638  trace.info() << "(" << nbok << "/" << nb << ") Plane=" << ppplane
639  << std::endl;
640  Point pppt4 = pppt3 + Point( 0, 0, 1 );
641  bool pppt4_inside = ppplane.extend( pt4 );
642  ++nb, nbok += pppt4_inside == true ? 1 : 0;
643  trace.info() << "(" << nbok << "/" << nb << ") Plane=" << ppplane
644  << std::endl;
645  trace.endBlock();
646  }
647 
648  {
649  trace.beginBlock ( "Testing block: ChordNaivePlaneComputer vertical instantiation 2." );
650  NaivePlaneComputer pplane;
651  pplane.init( 1, 1, 1 );
652  Point ppt0( -6, -3, 5 );
653  bool ppt0_inside = pplane.extend( ppt0 );
654  ++nb, nbok += ppt0_inside == true ? 1 : 0;
655  trace.info() << "(" << nbok << "/" << nb << ") Plane=" << pplane
656  << std::endl;
657  Point ppt1( 4, 4, -5 );
658  bool ppt1_inside = pplane.extend( ppt1 );
659  ++nb, nbok += ppt1_inside == true ? 1 : 0;
660  trace.info() << "(" << nbok << "/" << nb << ") Plane=" << pplane
661  << std::endl;
662  Point ppt2( -5, -2, 4 );
663  bool ppt2_inside = pplane.extend( ppt2 );
664  ++nb, nbok += ppt2_inside == true ? 1 : 0;
665  trace.info() << "(" << nbok << "/" << nb << ") Plane=" << pplane
666  << std::endl;
667  trace.endBlock();
668  }
669 
670  return nbok == nb;
671 }
672 
673 template <typename NaivePlaneComputer>
674 bool
675 checkManyPlanes( unsigned int diameter,
676  unsigned int nbplanes,
677  unsigned int nbpoints )
678 {
679  unsigned int nbok = 0;
680  unsigned int nb = 0;
681  typedef typename NaivePlaneComputer::InternalScalar Scalar;
682  stringstream ss (stringstream::out);
683  ss << "Testing block: Diameter is " << diameter << ". Check " << nbplanes << " planes with " << nbpoints << " points each.";
684  trace.beginBlock ( ss.str() );
685  ++nb, nbok += checkPlanes<Scalar,NaivePlaneComputer>( nbplanes, diameter, nbpoints ) ? 1 : 0;
686  trace.info() << "(" << nbok << "/" << nb
687  << ") checkPlanes<Scalar,NaivePlaneComputer>()"
688  << std::endl;
689  trace.endBlock();
690  return nbok == nb;
691 }
692 
693 template <typename GenericNaivePlaneComputer>
694 bool
695 checkExtendWithManyPoints( unsigned int diameter,
696  unsigned int nbplanes,
697  unsigned int nbpoints )
698 {
699  unsigned int nbok = 0;
700  unsigned int nb = 0;
701  typedef typename GenericNaivePlaneComputer::InternalScalar Integer;
702  typedef typename GenericNaivePlaneComputer::Point Point;
703  typedef typename Point::Coordinate PointInteger;
705 
706  trace.beginBlock( "checkExtendWithManyPoints" );
707  for ( unsigned int j = 0; j < nbplanes; ++j )
708  {
709  Integer a = getRandomInteger<Integer>( (Integer) 0, (Integer) diameter / 2 );
710  Integer b = getRandomInteger<Integer>( (Integer) 0, (Integer) diameter / 2 );
711  Integer c = getRandomInteger<Integer>( (Integer) 1, (Integer) diameter / 2 );
712  Integer d = getRandomInteger<Integer>( (Integer) 0, (Integer) diameter / 2 );
713  GenericNaivePlaneComputer plane;
714  Dimension axis;
715  if ( ( a >= b ) && ( a >= c ) ) axis = 0;
716  else if ( ( b >= a ) && ( b >= c ) ) axis = 1;
717  else axis = 2;
718  plane.init( 1, 1 );
719 
720  std::vector<Point> pts;
721  for ( unsigned int i = 0; i < nbpoints; ++i )
722  {
723  Point p;
724  p[ 0 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
725  p[ 1 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
726  p[ 2 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
727  Integer x = (Integer) p[ 0 ];
728  Integer y = (Integer) p[ 1 ];
729  Integer z = (Integer) p[ 2 ];
730  switch( axis ) {
731  case 0: p[ 0 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - b * y - c * z, a ) ); break;
732  case 1: p[ 1 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - c * z, b ) ); break;
733  case 2: p[ 2 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - b * y, c ) ); break;
734  }
735  pts.push_back( p );
736  }
737  ++nb, nbok += plane.isExtendable( pts.begin(), pts.end() ); // should be ok
738  trace.info() << "(" << nbok << "/" << nb
739  << ") plane.isExtendable( pts.begin(), pts.end() )"
740  << std::endl;
741  Point & any0 = pts[ getRandomInteger<int>( 0, pts.size() ) ];
742  pts.push_back( any0 + Point(1,0,0) );
743  Point & any1 = pts[ getRandomInteger<int>( 0, pts.size() ) ];
744  pts.push_back( any1 + Point(0,1,0) );
745  Point & any2 = pts[ getRandomInteger<int>( 0, pts.size() ) ];
746  pts.push_back( any2 + Point(0,0,1) );
747  bool check = ! plane.isExtendable( pts.begin(), pts.end() ); // should not be ok
748  ++nb, nbok += check ? 1 : 0;
749  trace.info() << "(" << nbok << "/" << nb
750  << ") ! plane.isExtendable( pts.begin(), pts.end() )"
751  << std::endl;
752  if ( ! check )
753  trace.warning() << plane << " last=" << pts.back() << std::endl
754  << "a=" << a << " b=" << b << " c=" << c << " d=" << d << std::endl;
755  ++nb, nbok += plane.extend( pts.begin(), pts.end() - 3 ); // should be ok
756  trace.info() << "(" << nbok << "/" << nb
757  << ") plane.extend( pts.begin(), pts.end() - 3)"
758  << std::endl;
759  ++nb, nbok += ! plane.extend( pts.end() - 3, pts.end() ); // should not be ok
760  trace.info() << "(" << nbok << "/" << nb
761  << ") ! plane.extend( pts.end() - 3, pts.end() )"
762  << std::endl;
763  }
764  trace.endBlock();
765  return nb == nbok;
766 }
767 
768 
769 
771 // Standard services - public :
772 
773 int main( int /*argc*/, char**/* argv */)
774 {
775  using namespace Z3i;
776 
777  trace.beginBlock ( "Testing class ChordNaivePlaneComputer" );
778  bool res = true
780  && checkManyPlanes<ChordNaivePlaneComputer<Z3i::Space, Z3i::Point, DGtal::int32_t> >( 4, 100, 200 )
782  && checkManyPlanes<ChordNaivePlaneComputer<Z3i::Space, Z3i::Point, DGtal::int32_t> >( 20, 100, 200 )
784  && checkManyPlanes<ChordNaivePlaneComputer<Z3i::Space, Z3i::Point, DGtal::int64_t> >( 2000, 100, 200 )
786  && checkExtendWithManyPoints<ChordGenericNaivePlaneComputer<Z3i::Space, Z3i::Point, DGtal::int64_t> >( 100, 100, 200 );
787 
788  trace.emphase() << ( res ? "Passed." : "Error." ) << endl;
789  trace.endBlock();
790  return res ? 0 : 1;
791 }
792 // //
void beginBlock(const std::string &keyword="")
void init(Dimension axis, InternalInteger diameter, InternalInteger widthNumerator=NumberTraits< InternalInteger >::ONE, InternalInteger widthDenominator=NumberTraits< InternalInteger >::ONE)
MyDigitalSurface::ConstIterator ConstIterator
Trace trace
Definition: Common.h:137
COBANaivePlaneComputer< Z3, InternalInteger > NaivePlaneComputer
ConstIterator begin() const
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:120
bool checkGenericPlane(Integer a, Integer b, Integer c, Integer d, int diameter, unsigned int nbtries)
Aim: A class that contains the chord-based algorithm for recognizing pieces of digital planes of give...
bool checkExtendWithManyPoints(unsigned int diameter, unsigned int nbplanes, unsigned int nbpoints)
STL namespace.
double endBlock()
bool isExtendable(const Point &p) const
ConstIterator end() const
Aim: Gathers several functions useful for concept checks.
Aim: Defines the concept describing an object that computes some primitive from input points given gr...
static Integer abs(IntegerParamType a)
Aim: Defines a predicate on a point.
std::ostream & emphase()
bool extend(const Point &p)
Aim: The traits class for all models of Cinteger.
Definition: NumberTraits.h:69
Integer getRandomInteger(Integer first, Integer after_last)
bool checkManyPlanes(unsigned int diameter, unsigned int nbplanes, unsigned int nbpoints)
bool checkWidth(Integer a, Integer b, Integer c, Integer d, int diameter, unsigned int nbtries)
bool checkPlane(Integer a, Integer b, Integer c, Integer d, int diameter, unsigned int nbtries)
bool checkWidths(unsigned int nbplanes, int diameter, unsigned int nbtries)
DGtal is the top-level namespace which contains all DGtal functions and types.
MyPointD Point
Definition: testClone2.cpp:383
bool checkPlanes(unsigned int nbplanes, int diameter, unsigned int nbtries)
int main(int, char **)
Aim: A class that recognizes pieces of digital planes of given axis width. When the width is 1...
bool testChordNaivePlaneComputer()
std::ostream & warning()
std::ostream & info()
Go to http://www.sgi.com/tech/stl/ForwardContainer.html.
Definition: Boost.dox:110
boost::int64_t int64_t
signed 94-bit integer.
Definition: BasicTypes.h:74
Integer ceilDiv(IntegerParamType na, IntegerParamType nb) const
bool checkPlaneGroupExtension(Integer a, Integer b, Integer c, Integer d, int diameter, unsigned int nbtries)