DGtal  1.3.beta
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;
66  IntegerComputer<Integer> ic;
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;
172  IntegerComputer<Integer> ic;
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;
288  IntegerComputer<Integer> ic;
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;
431  IntegerComputer<Integer> ic;
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 
525 bool testChordNaivePlaneComputer()
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;
532  typedef ChordNaivePlaneComputer<Space, Point, Integer> NaivePlaneComputer;
533  typedef ChordGenericNaivePlaneComputer<Space, Point, Integer> GenericNaivePlaneComputer;
534 
535  BOOST_CONCEPT_ASSERT(( CAdditivePrimitiveComputer< NaivePlaneComputer > ));
536  BOOST_CONCEPT_ASSERT(( CAdditivePrimitiveComputer< GenericNaivePlaneComputer > ));
537  BOOST_CONCEPT_ASSERT(( boost::ForwardContainer< NaivePlaneComputer > ));
538  BOOST_CONCEPT_ASSERT(( boost::ForwardContainer< GenericNaivePlaneComputer > ));
539  BOOST_CONCEPT_ASSERT(( CPointPredicate< NaivePlaneComputer::Primitive > ));
540  BOOST_CONCEPT_ASSERT(( CPointPredicate< GenericNaivePlaneComputer::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;
704  IntegerComputer<Integer> ic;
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
779  && testChordNaivePlaneComputer()
780  && checkManyPlanes<ChordNaivePlaneComputer<Z3i::Space, Z3i::Point, DGtal::int32_t> >( 4, 100, 200 )
781  && checkManyPlanes<ChordNaivePlaneComputer<Z3i::Space, Z3i::Point, DGtal::int32_t> >( 8, 100, 200 )
782  && checkManyPlanes<ChordNaivePlaneComputer<Z3i::Space, Z3i::Point, DGtal::int32_t> >( 20, 100, 200 )
783  && checkManyPlanes<ChordNaivePlaneComputer<Z3i::Space, Z3i::Point, DGtal::int32_t> >( 100, 100, 200 )
784  && checkManyPlanes<ChordNaivePlaneComputer<Z3i::Space, Z3i::Point, DGtal::int64_t> >( 2000, 100, 200 )
785  && checkWidths<DGtal::int64_t, ChordNaivePlaneComputer<Z3i::Space, Z3i::Point, DGtal::int64_t> >( 100, 1000000, 1000 )
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 // //
DGtal::ChordNaivePlaneComputer
Aim: A class that contains the chord-based algorithm for recognizing pieces of digital planes of give...
Definition: ChordNaivePlaneComputer.h:157
boost::ForwardContainer
Go to http://www.sgi.com/tech/stl/ForwardContainer.html.
Definition: Boost.dox:110
ConstIterator
MyDigitalSurface::ConstIterator ConstIterator
Definition: greedy-plane-segmentation-ex2.cpp:93
checkPlanes
bool checkPlanes(unsigned int nbplanes, int diameter, unsigned int nbtries)
Definition: testChordNaivePlaneComputer.cpp:388
main
int main(int, char **)
Definition: testChordNaivePlaneComputer.cpp:773
testChordNaivePlaneComputer
bool testChordNaivePlaneComputer()
Definition: testChordNaivePlaneComputer.cpp:525
NaivePlaneComputer
COBANaivePlaneComputer< Z3, InternalInteger > NaivePlaneComputer
Definition: greedy-plane-segmentation-ex2.cpp:88
DGtal::COBANaivePlaneComputer::init
void init(Dimension axis, InternalInteger diameter, InternalInteger widthNumerator=NumberTraits< InternalInteger >::ONE, InternalInteger widthDenominator=NumberTraits< InternalInteger >::ONE)
DGtal::Trace::endBlock
double endBlock()
checkPlaneGroupExtension
bool checkPlaneGroupExtension(Integer a, Integer b, Integer c, Integer d, int diameter, unsigned int nbtries)
Definition: testChordNaivePlaneComputer.cpp:167
checkPlane
bool checkPlane(Integer a, Integer b, Integer c, Integer d, int diameter, unsigned int nbtries)
Definition: testChordNaivePlaneComputer.cpp:61
checkManyPlanes
bool checkManyPlanes(unsigned int diameter, unsigned int nbplanes, unsigned int nbpoints)
Definition: testChordNaivePlaneComputer.cpp:675
DGtal::COBANaivePlaneComputer::extend
bool extend(const Point &p)
checkGenericPlane
bool checkGenericPlane(Integer a, Integer b, Integer c, Integer d, int diameter, unsigned int nbtries)
Definition: testChordNaivePlaneComputer.cpp:283
DGtal::IntegerComputer::ceilDiv
Integer ceilDiv(IntegerParamType na, IntegerParamType nb) const
DGtal::NumberTraits
Aim: The traits class for all models of Cinteger.
Definition: NumberTraits.h:562
DGtal::IntegerComputer< Integer >
DGtal::Trace::emphase
std::ostream & emphase()
DGtal::trace
Trace trace
Definition: Common.h:154
DGtal::COBANaivePlaneComputer::end
ConstIterator end() const
DGtal::Dimension
DGtal::uint32_t Dimension
Definition: Common.h:137
DGtal::Trace::beginBlock
void beginBlock(const std::string &keyword="")
DGtal::IntegerComputer::abs
static Integer abs(IntegerParamType a)
DGtal::SpaceND
Definition: SpaceND.h:95
DGtal::Trace::info
std::ostream & info()
DGtal::concepts
Aim: Gathers several functions useful for concept checks.
Definition: CPositiveIrreducibleFraction.h:54
DGtal
DGtal is the top-level namespace which contains all DGtal functions and types.
checkExtendWithManyPoints
bool checkExtendWithManyPoints(unsigned int diameter, unsigned int nbplanes, unsigned int nbpoints)
Definition: testChordNaivePlaneComputer.cpp:695
DGtal::ChordGenericNaivePlaneComputer
Aim: A class that recognizes pieces of digital planes of given axis width. When the width is 1,...
Definition: ChordGenericNaivePlaneComputer.h:118
checkWidths
bool checkWidths(unsigned int nbplanes, int diameter, unsigned int nbtries)
Definition: testChordNaivePlaneComputer.cpp:494
Integer
Point::Coordinate Integer
Definition: examplePlaneProbingParallelepipedEstimator.cpp:44
DGtal::COBANaivePlaneComputer::begin
ConstIterator begin() const
DGtal::PointVector< dim, Integer >
checkWidth
bool checkWidth(Integer a, Integer b, Integer c, Integer d, int diameter, unsigned int nbtries)
Definition: testChordNaivePlaneComputer.cpp:426
DGtal::concepts::CPointPredicate
Aim: Defines a predicate on a point.
Definition: CPointPredicate.h:80
Space
SpaceND< 2 > Space
Definition: testSimpleRandomAccessRangeFromPoint.cpp:42
DGtal::int64_t
boost::int64_t int64_t
signed 94-bit integer.
Definition: BasicTypes.h:74
DGtal::COBANaivePlaneComputer::isExtendable
bool isExtendable(const Point &p) const
DGtal::Trace::warning
std::ostream & warning()
DGtal::COBANaivePlaneComputer< Z3, InternalInteger >
Point
MyPointD Point
Definition: testClone2.cpp:383
DGtal::COBANaivePlaneComputer< Z3, InternalInteger >::ConstIterator
PointSet::const_iterator ConstIterator
Definition: COBANaivePlaneComputer.h:142
DGtal::concepts::CAdditivePrimitiveComputer
Aim: Defines the concept describing an object that computes some primitive from input points given gr...
Definition: CAdditivePrimitiveComputer.h:99
getRandomInteger
Integer getRandomInteger(Integer first, Integer after_last)
Definition: testChordNaivePlaneComputer.cpp:50
Generated on Thu May 19 2022 18:08:46 for DGtal by   doxygen 1.8.17