DGtal  1.5.beta
testChordNaivePlaneComputer.cpp
Go to the documentation of this file.
1 
31 #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>
51 {
52  Integer r = (Integer) rand();
53  return ( r % (after_last - first) ) + first;
54 }
55 
59 template <typename Integer, typename NaivePlaneComputer>
60 bool
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 ] = (PointInteger)NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - b * y - c * z, a ) ); break;
94  case 1: p[ 1 ] = (PointInteger)NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - c * z, b ) ); break;
95  case 2: p[ 2 ] = (PointInteger)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 ] = (PointInteger)NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - b * y - c * z, a ) ); break;
136  case 1: p[ 1 ] = (PointInteger)NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - c * z, b ) ); break;
137  case 2: p[ 2 ] = (PointInteger)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
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 ] = (PointInteger)NumberTraits<Integer>::castToInt64_t
204  ( ic.ceilDiv( d - b * y - c * z, a ) ); break;
205  case 1: pp[ 1 ] = (PointInteger)NumberTraits<Integer>::castToInt64_t
206  ( ic.ceilDiv( d - a * x - c * z, b ) ); break;
207  case 2: pp[ 2 ] = (PointInteger)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 ] = (PointInteger)NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - b * y - c * z, a ) ); break;
250  case 1: p[ 1 ] = (PointInteger)NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - c * z, b ) ); break;
251  case 2: p[ 2 ] = (PointInteger)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
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 ] = (PointInteger)NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - b * y - c * z, a ) ); break;
316  case 1: p[ 1 ] = (PointInteger)NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - c * z, b ) ); break;
317  case 2: p[ 2 ] = (PointInteger)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 ] = (PointInteger)NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - b * y - c * z, a ) ); break;
358  case 1: p[ 1 ] = (PointInteger)NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - c * z, b ) ); break;
359  case 2: p[ 2 ] = (PointInteger)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
427  int diameter, unsigned int nbtries )
428 {
429  typedef typename NaivePlaneComputer::Point Point;
430  typedef typename NaivePlaneComputer::InternalScalar InternalScalar;
431  typedef typename Point::Component PointInteger;
433  Integer absA = ic.abs( a );
434  Integer absB = ic.abs( b );
435  Integer absC = ic.abs( c );
436  Integer x, y, z;
437  Dimension axis;
438  if ( ( absA >= absB ) && ( absA >= absC ) )
439  axis = 0;
440  else if ( ( absB >= absA ) && ( absB >= absC ) )
441  axis = 1;
442  else
443  axis = 2;
444  // Checks that points within the naive plane are correctly recognized.
445  unsigned int nb = 0;
446  unsigned int nbok = 0;
447  std::vector<Point> points( nbtries );
448  for ( unsigned int i = 0; i != nbtries; ++i )
449  {
450  Point & p = points[ i ];
451  p[ 0 ] = getRandomInteger<Integer>( -diameter+1, diameter );
452  p[ 1 ] = getRandomInteger<Integer>( -diameter+1, diameter );
453  p[ 2 ] = getRandomInteger<Integer>( -diameter+1, diameter );
454  x = (Integer) p[ 0 ];
455  y = (Integer) p[ 1 ];
456  z = (Integer) p[ 2 ];
457  switch ( axis ) {
458  case 0: p[ 0 ] = (PointInteger)NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - b * y - c * z, a ) ); break;
459  case 1: p[ 1 ] = (PointInteger)NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - c * z, b ) ); break;
460  case 2: p[ 2 ] = (PointInteger)NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - b * y, c ) ); break;
461  }
462  }
463  trace.beginBlock( "Computing axis width." );
464  trace.info() << "- plane is "
465  << d << " <= " << a << "*x"
466  << "+" << b << "*y"
467  << "+" << c << "*z"
468  << " <= d + max(|a|,|b|,|c|)"
469  << std::endl;
470  trace.info() << "- " << points.size() << " points tested in diameter " << diameter
471  << std::endl;
472  double min = -1.0;
473  for ( unsigned int i = 0; i < 3; ++i )
474  {
475  std::pair<InternalScalar, InternalScalar> width
476  = NaivePlaneComputer::computeAxisWidth( i, points.begin(), points.end() );
477  double wn = NumberTraits<InternalScalar>::castToDouble( width.first );
478  double wd = NumberTraits<InternalScalar>::castToDouble( width.second );
479  trace.info() << " (" << i << ") width=" << (wn/wd) << std::endl;
480  if ( min < 0.0 ) min = wn/wd;
481  else if ( wn/wd < min ) min = wn/wd;
482  }
483  ++nb; nbok += (min < 1.0 ) ? 1 : 0;
484  trace.info() << "(" << nbok << "/" << nb << ") min width = " << min
485  << " < 1.0" << std::endl;
486  ++nb; nbok += (0.9 < min ) ? 1 : 0;
487  trace.info() << "(" << nbok << "/" << nb << ") min width = " << min
488  << " > 0.9" << std::endl;
489  trace.endBlock();
490  return nb == nbok;
491 }
492 
493 template <typename Integer, typename NaivePlaneComputer>
494 bool
495 checkWidths( unsigned int nbplanes, int diameter, unsigned int nbtries )
496 {
497  //using namespace Z3i;
498  //typedef ChordNaivePlaneComputer<Z3, Integer> NaivePlaneComputer;
499  unsigned int nb = 0;
500  unsigned int nbok = 0;
501  for ( unsigned int nbp = 0; nbp < nbplanes; ++nbp )
502  {
503  Integer a = getRandomInteger<Integer>( (Integer) 0, (Integer) diameter / 2 );
504  Integer b = getRandomInteger<Integer>( (Integer) 0, (Integer) diameter / 2 );
505  Integer c = getRandomInteger<Integer>( (Integer) 0, (Integer) diameter / 2 );
506  Integer d = getRandomInteger<Integer>( (Integer) 0, (Integer) diameter / 2 );
507  if ( ( a != 0 ) || ( b != 0 ) || ( c != 0 ) )
508  {
509  ++nb; nbok += checkWidth<Integer, NaivePlaneComputer>( a, b, c, d, diameter, nbtries ) ? 1 : 0;
510  if ( nb != nbok )
511  {
512  std::cerr << "[ERROR] (checkWidth) for plane " << a << " * x + "
513  << b << " * y + " << c << " * z = " << d << std::endl;
514  break;
515  }
516  }
517  }
518  return nb == nbok;
519 }
520 
521 
527 {
528  unsigned int nbok = 0;
529  unsigned int nb = 0;
530  typedef DGtal::int64_t Integer;
531  typedef DGtal::Z3i::Z3 Space;
532  typedef DGtal::Z3i::Point Point;
534  typedef ChordGenericNaivePlaneComputer<Space, Point, Integer> GenericNaivePlaneComputer;
535 
536  BOOST_CONCEPT_ASSERT(( CAdditivePrimitiveComputer< NaivePlaneComputer > ));
538  BOOST_CONCEPT_ASSERT(( boost::ForwardContainer< NaivePlaneComputer > ));
540  BOOST_CONCEPT_ASSERT(( CPointPredicate< NaivePlaneComputer::Primitive > ));
542 
543  trace.beginBlock ( "Testing block: ChordNaivePlaneComputer instantiation." );
544  NaivePlaneComputer plane;
545  Point pt0( 0, 0, 0 );
546  plane.init( 2, 1, 1 );
547  bool pt0_inside = plane.extend( pt0 );
548  ++nb; nbok += pt0_inside == true ? 1 : 0;
549  trace.info() << "(" << nbok << "/" << nb << ") Plane=" << plane
550  << std::endl;
551  Point pt1( Point( 2, 0, 0 ) );
552  bool pt1_inside = plane.extend( pt1 );
553  ++nb; nbok += pt1_inside == true ? 1 : 0;
554  trace.info() << "(" << nbok << "/" << nb << ") add " << pt1
555  << " Plane=" << plane << std::endl;
556  Point pt2( Point( 0, 2, 2 ) );
557  bool pt2_inside = plane.extend( pt2 );
558  ++nb; nbok += pt2_inside == true ? 1 : 0;
559  trace.info() << "(" << nbok << "/" << nb << ") add " << pt2
560  << " Plane=" << plane << std::endl;
561 
562  Point pt3( Point( 1, 1, 1 ) );
563  bool pt3_inside = plane.extend( pt3 );
564  ++nb; nbok += pt3_inside == true ? 1 : 0;
565  trace.info() << "(" << nbok << "/" << nb << ") add " << pt3
566  << " Plane=" << plane << std::endl;
567 
568  Point pt4( Point( -10, -10, 10 ) );
569  bool pt4_inside = plane.extend( pt4 );
570  ++nb; nbok += pt4_inside == false ? 1 : 0;
571  trace.info() << "(" << nbok << "/" << nb << ") impossible add " << pt4
572  << " Plane=" << plane << std::endl;
573 
574  Point pt5 = pt2 + Point( 1, 0, 1 );
575  bool pt5_inside = plane.extend( pt5 );
576  ++nb; nbok += pt5_inside == true ? 1 : 0;
577  trace.info() << "(" << nbok << "/" << nb << ") add " << pt5
578  << " Plane=" << plane << std::endl;
579 
580  Point pt6 = pt5 + Point( 6, 0, 2 );
581  bool pt6_inside = plane.extend( pt6 );
582  ++nb; nbok += pt6_inside == true ? 1 : 0;
583  trace.info() << "(" << nbok << "/" << nb << ") add " << pt6
584  << " Plane=" << plane << std::endl;
585 
586  NaivePlaneComputer plane2;
587  plane2.init( 2, 1, 1 );
588  plane2.extend( Point( 10, 0, 0 ) );
589  plane2.extend( Point( 0, 8, 0 ) );
590  plane2.extend( Point( 0, 0, 6 ) );
591  trace.info() << "(" << nbok << "/" << nb << ") "
592  << " Plane2=" << plane2 << std::endl;
593 
594  ++nb; nbok += checkPlane<Integer,NaivePlaneComputer>( 11, 5, 19, 20, 100, 100 ) ? 1 : 0;
595  trace.info() << "(" << nbok << "/" << nb
596  << ") checkPlane<Integer,NaivePlaneComputer>( 11, 5, 19, 20, 100, 100 )"
597  << std::endl;
598 
599  ++nb; nbok += checkGenericPlane<Integer,GenericNaivePlaneComputer>( 11, 5, 19, 20, 100, 100 ) ? 1 : 0;
600  trace.info() << "(" << nbok << "/" << nb
601  << ") checkGenericPlane<Integer,GenericNaivePlaneComputer>( 11, 5, 19, 20, 100, 100 )"
602  << std::endl;
603  ++nb; nbok += checkGenericPlane<Integer,GenericNaivePlaneComputer>( 17, 33, 7, 10, 100, 100 ) ? 1 : 0;
604  trace.info() << "(" << nbok << "/" << nb
605  << ") checkGenericPlane<Integer,GenericNaivePlaneComputer>( 17, 33, 7, 10, 100, 100 )"
606  << std::endl;
607  ++nb; nbok += checkPlane<Integer,NaivePlaneComputer>( 15, 8, 13, 15, 100, 100 ) ? 1 : 0;
608  trace.info() << "(" << nbok << "/" << nb
609  << ") checkPlane<Integer,NaivePlaneComputer>( 15, 8, 13, 15, 100, 100 )"
610  << std::endl;
611  ++nb; nbok += checkGenericPlane<Integer,GenericNaivePlaneComputer>( 15, 8, 13, 15, 100, 100 ) ? 1 : 0;
612  trace.info() << "(" << nbok << "/" << nb
613  << ") checkGenericPlane<Integer,GenericNaivePlaneComputer>( 15, 8, 13, 15, 100, 100 )"
614  << std::endl;
615  trace.endBlock();
616 
617  {
618  trace.beginBlock ( "Testing block: ChordNaivePlaneComputer vertical instantiation." );
619  NaivePlaneComputer ppplane;
620  Point pppt0( 0, 0, 0 );
621  ppplane.init( 2, 5, 2 );
622  bool pppt0_inside = ppplane.extend( pppt0 );
623  ++nb; nbok += pppt0_inside == true ? 1 : 0;
624  trace.info() << "(" << nbok << "/" << nb << ") Plane=" << ppplane
625  << std::endl;
626  Point pppt1( 3, 2, 2 );
627  bool pppt1_inside = ppplane.extend( pppt1 );
628  ++nb; nbok += pppt1_inside == true ? 1 : 0;
629  trace.info() << "(" << nbok << "/" << nb << ") Plane=" << ppplane
630  << std::endl;
631  Point pppt2( 0, 0, 1 );
632  bool pppt2_inside = ppplane.extend( pppt2 );
633  ++nb; nbok += pppt2_inside == true ? 1 : 0;
634  trace.info() << "(" << nbok << "/" << nb << ") Plane=" << ppplane
635  << std::endl;
636  Point pppt3 = pppt1 + Point( 0, 0, 1 );
637  bool pppt3_inside = ppplane.extend( pppt3 );
638  ++nb; nbok += pppt3_inside == true ? 1 : 0;
639  trace.info() << "(" << nbok << "/" << nb << ") Plane=" << ppplane
640  << std::endl;
641  Point pppt4 = pppt3 + Point( 0, 0, 1 );
642  bool pppt4_inside = ppplane.extend( pt4 );
643  ++nb; nbok += pppt4_inside == true ? 1 : 0;
644  trace.info() << "(" << nbok << "/" << nb << ") Plane=" << ppplane
645  << std::endl;
646  trace.endBlock();
647  }
648 
649  {
650  trace.beginBlock ( "Testing block: ChordNaivePlaneComputer vertical instantiation 2." );
651  NaivePlaneComputer pplane;
652  pplane.init( 1, 1, 1 );
653  Point ppt0( -6, -3, 5 );
654  bool ppt0_inside = pplane.extend( ppt0 );
655  ++nb; nbok += ppt0_inside == true ? 1 : 0;
656  trace.info() << "(" << nbok << "/" << nb << ") Plane=" << pplane
657  << std::endl;
658  Point ppt1( 4, 4, -5 );
659  bool ppt1_inside = pplane.extend( ppt1 );
660  ++nb; nbok += ppt1_inside == true ? 1 : 0;
661  trace.info() << "(" << nbok << "/" << nb << ") Plane=" << pplane
662  << std::endl;
663  Point ppt2( -5, -2, 4 );
664  bool ppt2_inside = pplane.extend( ppt2 );
665  ++nb; nbok += ppt2_inside == true ? 1 : 0;
666  trace.info() << "(" << nbok << "/" << nb << ") Plane=" << pplane
667  << std::endl;
668  trace.endBlock();
669  }
670 
671  return nbok == nb;
672 }
673 
674 template <typename NaivePlaneComputer>
675 bool
676 checkManyPlanes( unsigned int diameter,
677  unsigned int nbplanes,
678  unsigned int nbpoints )
679 {
680  unsigned int nbok = 0;
681  unsigned int nb = 0;
682  typedef typename NaivePlaneComputer::InternalScalar Scalar;
683  stringstream ss (stringstream::out);
684  ss << "Testing block: Diameter is " << diameter << ". Check " << nbplanes << " planes with " << nbpoints << " points each.";
685  trace.beginBlock ( ss.str() );
686  ++nb; nbok += checkPlanes<Scalar,NaivePlaneComputer>( nbplanes, diameter, nbpoints ) ? 1 : 0;
687  trace.info() << "(" << nbok << "/" << nb
688  << ") checkPlanes<Scalar,NaivePlaneComputer>()"
689  << std::endl;
690  trace.endBlock();
691  return nbok == nb;
692 }
693 
694 template <typename GenericNaivePlaneComputer>
695 bool
696 checkExtendWithManyPoints( unsigned int diameter,
697  unsigned int nbplanes,
698  unsigned int nbpoints )
699 {
700  unsigned int nbok = 0;
701  unsigned int nb = 0;
702  typedef typename GenericNaivePlaneComputer::InternalScalar Integer;
703  typedef typename GenericNaivePlaneComputer::Point Point;
704  typedef typename Point::Coordinate PointInteger;
706 
707  trace.beginBlock( "checkExtendWithManyPoints" );
708  for ( unsigned int j = 0; j < nbplanes; ++j )
709  {
710  Integer a = getRandomInteger<Integer>( (Integer) 0, (Integer) diameter / 2 );
711  Integer b = getRandomInteger<Integer>( (Integer) 0, (Integer) diameter / 2 );
712  Integer c = getRandomInteger<Integer>( (Integer) 1, (Integer) diameter / 2 );
713  Integer d = getRandomInteger<Integer>( (Integer) 0, (Integer) diameter / 2 );
714  GenericNaivePlaneComputer plane;
715  Dimension axis;
716  if ( ( a >= b ) && ( a >= c ) ) axis = 0;
717  else if ( ( b >= a ) && ( b >= c ) ) axis = 1;
718  else axis = 2;
719  plane.init( 1, 1 );
720 
721  std::vector<Point> pts;
722  for ( unsigned int i = 0; i < nbpoints; ++i )
723  {
724  Point p;
725  p[ 0 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
726  p[ 1 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
727  p[ 2 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
728  Integer x = (Integer) p[ 0 ];
729  Integer y = (Integer) p[ 1 ];
730  Integer z = (Integer) p[ 2 ];
731  switch( axis ) {
732  case 0: p[ 0 ] = (Integer)NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - b * y - c * z, a ) ); break;
733  case 1: p[ 1 ] = (Integer)NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - c * z, b ) ); break;
734  case 2: p[ 2 ] = (Integer)NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - b * y, c ) ); break;
735  }
736  pts.push_back( p );
737  }
738  ++nb; nbok += plane.isExtendable( pts.begin(), pts.end() ); // should be ok
739  trace.info() << "(" << nbok << "/" << nb
740  << ") plane.isExtendable( pts.begin(), pts.end() )"
741  << std::endl;
742  Point & any0 = pts[ getRandomInteger<int>( 0, pts.size() ) ];
743  pts.push_back( any0 + Point(1,0,0) );
744  Point & any1 = pts[ getRandomInteger<int>( 0, pts.size() ) ];
745  pts.push_back( any1 + Point(0,1,0) );
746  Point & any2 = pts[ getRandomInteger<int>( 0, pts.size() ) ];
747  pts.push_back( any2 + Point(0,0,1) );
748  bool check = ! plane.isExtendable( pts.begin(), pts.end() ); // should not be ok
749  ++nb; nbok += check ? 1 : 0;
750  trace.info() << "(" << nbok << "/" << nb
751  << ") ! plane.isExtendable( pts.begin(), pts.end() )"
752  << std::endl;
753  if ( ! check )
754  trace.warning() << plane << " last=" << pts.back() << std::endl
755  << "a=" << a << " b=" << b << " c=" << c << " d=" << d << std::endl;
756  ++nb; nbok += plane.extend( pts.begin(), pts.end() - 3 ); // should be ok
757  trace.info() << "(" << nbok << "/" << nb
758  << ") plane.extend( pts.begin(), pts.end() - 3)"
759  << std::endl;
760  ++nb; nbok += ! plane.extend( pts.end() - 3, pts.end() ); // should not be ok
761  trace.info() << "(" << nbok << "/" << nb
762  << ") ! plane.extend( pts.end() - 3, pts.end() )"
763  << std::endl;
764  }
765  trace.endBlock();
766  return nb == nbok;
767 }
768 
769 
770 
772 // Standard services - public :
773 
774 int main( int /*argc*/, char**/* argv */)
775 {
776  using namespace Z3i;
777 
778  trace.beginBlock ( "Testing class ChordNaivePlaneComputer" );
779  bool res = true
781  && checkManyPlanes<ChordNaivePlaneComputer<Z3i::Space, Z3i::Point, DGtal::int32_t> >( 4, 100, 200 )
783  && checkManyPlanes<ChordNaivePlaneComputer<Z3i::Space, Z3i::Point, DGtal::int32_t> >( 20, 100, 200 )
785  && checkManyPlanes<ChordNaivePlaneComputer<Z3i::Space, Z3i::Point, DGtal::int64_t> >( 2000, 100, 200 )
787  && checkExtendWithManyPoints<ChordGenericNaivePlaneComputer<Z3i::Space, Z3i::Point, DGtal::int64_t> >( 100, 100, 200 );
788 
789  trace.emphase() << ( res ? "Passed." : "Error." ) << endl;
790  trace.endBlock();
791  return res ? 0 : 1;
792 }
793 // //
void init(Dimension axis, InternalInteger diameter, InternalInteger widthNumerator=NumberTraits< InternalInteger >::ONE, InternalInteger widthDenominator=NumberTraits< InternalInteger >::ONE)
ConstIterator end() const
ConstIterator begin() const
bool isExtendable(const Point &p) const
bool extend(const Point &p)
Aim: A class that recognizes pieces of digital planes of given axis width. When the width is 1,...
Aim: A class that contains the chord-based algorithm for recognizing pieces of digital planes of give...
static Integer abs(IntegerParamType a)
Integer ceilDiv(IntegerParamType na, IntegerParamType nb) const
void beginBlock(const std::string &keyword="")
std::ostream & emphase()
std::ostream & info()
std::ostream & warning()
double endBlock()
COBANaivePlaneComputer< Z3, InternalInteger > NaivePlaneComputer
MyDigitalSurface::ConstIterator ConstIterator
Aim: Gathers several functions useful for concept checks.
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
DGtal::uint32_t Dimension
Definition: Common.h:136
Trace trace
Definition: Common.h:153
Aim: The traits class for all models of Cinteger.
Definition: NumberTraits.h:564
Aim: Defines the concept describing an object that computes some primitive from input points given gr...
Aim: Defines a predicate on a point.
Go to http://www.sgi.com/tech/stl/ForwardContainer.html.
Definition: Boost.dox:110
bool checkWidth(Integer a, Integer b, Integer c, Integer d, int diameter, unsigned int nbtries)
bool testChordNaivePlaneComputer()
int main(int, char **)
bool checkPlane(Integer a, Integer b, Integer c, Integer d, int diameter, unsigned int nbtries)
bool checkPlanes(unsigned int nbplanes, int diameter, unsigned int nbtries)
bool checkGenericPlane(Integer a, Integer b, Integer c, Integer d, int diameter, unsigned int nbtries)
bool checkExtendWithManyPoints(unsigned int diameter, unsigned int nbplanes, unsigned int nbpoints)
Integer getRandomInteger(Integer first, Integer after_last)
bool checkWidths(unsigned int nbplanes, int diameter, unsigned int nbtries)
bool checkPlaneGroupExtension(Integer a, Integer b, Integer c, Integer d, int diameter, unsigned int nbtries)
bool checkManyPlanes(unsigned int diameter, unsigned int nbplanes, unsigned int nbpoints)
MyPointD Point
Definition: testClone2.cpp:383