DGtal  0.6.devel
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
testCOBANaivePlane.cpp
1 
30 
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/COBANaivePlane.h"
37 #include "DGtal/geometry/surfaces/COBAGenericNaivePlane.h"
39 
40 using namespace std;
41 using namespace DGtal;
42 
44 // Functions for testing class COBANaivePlane.
46 
47 template <typename Integer>
48 Integer getRandomInteger( const Integer & first, const Integer & after_last )
49 {
50  Integer r = (Integer) random();
51  return ( r % (after_last - first) ) + first;
52 }
53 
57 template <typename Integer, typename NaivePlane>
58 bool
59 checkPlane( Integer a, Integer b, Integer c, Integer d,
60  int diameter, unsigned int nbtries )
61 {
62  typedef typename NaivePlane::Point Point;
63  typedef typename Point::Component PointInteger;
64  IntegerComputer<Integer> ic;
65  Integer absA = ic.abs( a );
66  Integer absB = ic.abs( b );
67  Integer absC = ic.abs( c );
68  Integer x, y, z;
69  Dimension axis;
70  if ( ( absA >= absB ) && ( absA >= absC ) )
71  axis = 0;
72  else if ( ( absB >= absA ) && ( absB >= absC ) )
73  axis = 1;
74  else
75  axis = 2;
76  Point p;
77  NaivePlane plane;
78  plane.init( axis, diameter, 1, 1 );
79  // Checks that points within the naive plane are correctly recognized.
80  unsigned int nb = 0;
81  unsigned int nbok = 0;
82  while ( nb != nbtries )
83  {
84  p[ 0 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
85  p[ 1 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
86  p[ 2 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
87  x = (Integer) p[ 0 ];
88  y = (Integer) p[ 1 ];
89  z = (Integer) p[ 2 ];
90  switch ( axis ) {
91  case 0: p[ 0 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - b * y - c * z, a ) ); break;
92  case 1: p[ 1 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - c * z, b ) ); break;
93  case 2: p[ 2 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - b * y, c ) ); break;
94  }
95  bool ok_ext = plane.isExtendable( p ); // should be ok
96  bool ok = plane.extend( p ); // should be ok
97  ++nb, nbok += ok_ext ? 1 : 0;
98  ++nb, nbok += ok ? 1 : 0;
99  if ( ! ok )
100  {
101  std::cerr << "[ERROR] p=" << p << " NOT IN plane=" << plane << std::endl;
102  break;
103  }
104  if ( ! ok_ext )
105  {
106  std::cerr << "[ERROR] p=" << p << " was NOT extendable IN plane=" << plane << std::endl;
107  break;
108  }
109  // else
110  // std::cerr << "[OK] p=" << p << " IN plane=" << plane << std::endl;
111  }
112 
113  // Checks that points outside the naive plane are correctly recognized as outliers.
114  while ( nb != (nbtries * 11 ) / 10 )
115  {
116  p[ 0 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
117  p[ 1 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
118  p[ 2 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
119  x = (Integer) p[ 0 ];
120  y = (Integer) p[ 1 ];
121  z = (Integer) p[ 2 ];
122  switch ( axis ) {
123  case 0: p[ 0 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - b * y - c * z, a ) ); break;
124  case 1: p[ 1 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - c * z, b ) ); break;
125  case 2: p[ 2 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - b * y, c ) ); break;
126  }
127  PointInteger tmp = getRandomInteger<PointInteger>( 2, 5 )
128  * (2*getRandomInteger<PointInteger>( 0, 2 ) - 1 );
129  p[ axis ] += tmp;
130  bool ok_ext = ! plane.isExtendable( p ); // should *not* be ok
131  bool ok = ! plane.extend( p ); // should *not* be ok
132  ++nb, nbok += ok ? 1 : 0;
133  ++nb, nbok += ok_ext ? 1 : 0;
134  if ( ! ok )
135  {
136  std::cerr << "[ERROR] p=" << p << " IN plane=" << plane << std::endl;
137  break;
138  }
139  if ( ! ok_ext )
140  {
141  std::cerr << "[ERROR] p=" << p << " was extendable IN plane=" << plane << std::endl;
142  break;
143  }
144  // else
145  // std::cerr << "[OK] p=" << p << " IN plane=" << plane << std::endl;
146  }
147  return nb == nbok;
148 }
149 
150 
154 template <typename Integer, typename GenericNaivePlane>
155 bool
156 checkGenericPlane( Integer a, Integer b, Integer c, Integer d,
157  int diameter, unsigned int nbtries )
158 {
159  typedef typename GenericNaivePlane::Point Point;
160  typedef typename Point::Component PointInteger;
161  IntegerComputer<Integer> ic;
162  Integer absA = ic.abs( a );
163  Integer absB = ic.abs( b );
164  Integer absC = ic.abs( c );
165  Integer x, y, z;
166  Dimension axis;
167  if ( ( absA >= absB ) && ( absA >= absC ) )
168  axis = 0;
169  else if ( ( absB >= absA ) && ( absB >= absC ) )
170  axis = 1;
171  else
172  axis = 2;
173  Point p;
174  GenericNaivePlane plane;
175  plane.init( diameter, 1, 1 );
176  // Checks that points within the naive plane are correctly recognized.
177  unsigned int nb = 0;
178  unsigned int nbok = 0;
179  while ( nb != nbtries )
180  {
181  p[ 0 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
182  p[ 1 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
183  p[ 2 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
184  x = (Integer) p[ 0 ];
185  y = (Integer) p[ 1 ];
186  z = (Integer) p[ 2 ];
187  switch ( axis ) {
188  case 0: p[ 0 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - b * y - c * z, a ) ); break;
189  case 1: p[ 1 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - c * z, b ) ); break;
190  case 2: p[ 2 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - b * y, c ) ); break;
191  }
192  bool ok_ext = plane.isExtendable( p ); // should be ok
193  bool ok = plane.extend( p ); // should be ok
194  ++nb, nbok += ok_ext ? 1 : 0;
195  ++nb, nbok += ok ? 1 : 0;
196  if ( ! ok )
197  {
198  std::cerr << "[ERROR] p=" << p << " NOT IN plane=" << plane << std::endl;
199  break;
200  }
201  if ( ! ok_ext )
202  {
203  std::cerr << "[ERROR] p=" << p << " was NOT extendable IN plane=" << plane << std::endl;
204  break;
205  }
206  // else
207  // std::cerr << "[OK] p=" << p << " IN plane=" << plane << std::endl;
208  }
209 
210  // Checks that points outside the naive plane are correctly recognized as outliers.
211  while ( nb != (nbtries * 11 ) / 10 )
212  {
213  p[ 0 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
214  p[ 1 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
215  p[ 2 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
216  x = (Integer) p[ 0 ];
217  y = (Integer) p[ 1 ];
218  z = (Integer) p[ 2 ];
219  switch ( axis ) {
220  case 0: p[ 0 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - b * y - c * z, a ) ); break;
221  case 1: p[ 1 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - c * z, b ) ); break;
222  case 2: p[ 2 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - b * y, c ) ); break;
223  }
224  PointInteger tmp = getRandomInteger<PointInteger>( 2, 5 )
225  * (2*getRandomInteger<PointInteger>( 0, 2 ) - 1 );
226  p[ axis ] += tmp;
227  bool ok_ext = ! plane.isExtendable( p ); // should *not* be ok
228  bool ok = ! plane.extend( p ); // should *not* be ok
229  ++nb, nbok += ok ? 1 : 0;
230  ++nb, nbok += ok_ext ? 1 : 0;
231  if ( ! ok )
232  {
233  std::cerr << "[ERROR] p=" << p << " IN plane=" << plane << std::endl;
234  break;
235  }
236  if ( ! ok_ext )
237  {
238  std::cerr << "[ERROR] p=" << p << " was extendable IN plane=" << plane << std::endl;
239  break;
240  }
241  // else
242  // std::cerr << "[OK] p=" << p << " IN plane=" << plane << std::endl;
243  }
244  std::cerr << "plane = " << plane << std::endl;
245  return nb == nbok;
246 }
247 
248 
249 template <typename Integer, typename NaivePlane>
250 bool
251 checkPlanes( unsigned int nbplanes, int diameter, unsigned int nbtries )
252 {
253  //using namespace Z3i;
254  //typedef COBANaivePlane<Z3, Integer> NaivePlane;
255  unsigned int nb = 0;
256  unsigned int nbok = 0;
257  for ( unsigned int nbp = 0; nbp < nbplanes; ++nbp )
258  {
259  Integer a = getRandomInteger<Integer>( (Integer) 0, (Integer) diameter / 2 );
260  Integer b = getRandomInteger<Integer>( (Integer) 0, (Integer) diameter / 2 );
261  Integer c = getRandomInteger<Integer>( (Integer) 0, (Integer) diameter / 2 );
262  Integer d = getRandomInteger<Integer>( (Integer) 0, (Integer) diameter / 2 );
263  if ( ( a != 0 ) || ( b != 0 ) || ( c != 0 ) )
264  {
265  ++nb, nbok += checkPlane<Integer, NaivePlane>( a, b, c, d, diameter, nbtries ) ? 1 : 0;
266  if ( nb != nbok )
267  {
268  std::cerr << "[ERROR] for plane " << a << " * x + "
269  << b << " * y + " << c << " * z = " << d << std::endl;
270  break;
271  }
272  }
273  }
274  return nb == nbok;
275 }
276 
281 bool testCOBANaivePlane()
282 {
283  unsigned int nbok = 0;
284  unsigned int nb = 0;
285  using namespace Z3i;
286  typedef BigInteger Integer;
287  typedef COBANaivePlane<Z3, BigInteger> NaivePlane;
288  typedef COBAGenericNaivePlane<Z3, BigInteger> GenericNaivePlane;
289 
290  BOOST_CONCEPT_ASSERT(( CPointPredicate< NaivePlane > ));
291  BOOST_CONCEPT_ASSERT(( boost::ForwardContainer< NaivePlane > ));
292  BOOST_CONCEPT_ASSERT(( CPointPredicate< GenericNaivePlane > ));
293  BOOST_CONCEPT_ASSERT(( boost::ForwardContainer< GenericNaivePlane > ));
294 
295  trace.beginBlock ( "Testing block: COBANaivePlane instantiation." );
296  NaivePlane plane;
297  Point pt0( 0, 0, 0 );
298  plane.init( 2, 100, 3, 2 );
299  bool pt0_inside = plane.extend( pt0 );
300  trace.info() << "(" << nbok << "/" << nb << ") Plane=" << plane
301  << std::endl;
302  Point pt1( Point( 8, 1, 3 ) );
303  bool pt1_inside = plane.extend( pt1 );
304  ++nb, nbok += pt1_inside == true ? 1 : 0;
305  trace.info() << "(" << nbok << "/" << nb << ") add " << pt1
306  << " Plane=" << plane << std::endl;
307  Point pt2( Point( 2, 7, 1 ) );
308  bool pt2_inside = plane.extend( pt2 );
309  ++nb, nbok += pt2_inside == true ? 1 : 0;
310  trace.info() << "(" << nbok << "/" << nb << ") add " << pt2
311  << " Plane=" << plane << std::endl;
312 
313  Point pt3( Point( 0, 5, 17 ) );
314  bool pt3_inside = plane.extend( pt3 );
315  ++nb, nbok += pt3_inside == false ? 1 : 0;
316  trace.info() << "(" << nbok << "/" << nb << ") add " << pt3
317  << " Plane=" << plane << std::endl;
318 
319  Point pt4( Point( -10, -10, 10 ) );
320  bool pt4_inside = plane.extend( pt4 );
321  ++nb, nbok += pt4_inside == false ? 1 : 0;
322  trace.info() << "(" << nbok << "/" << nb << ") add " << pt4
323  << " Plane=" << plane << std::endl;
324 
325  Point pt5 = pt0 + pt1 + pt2 + Point( 0, 0, 2 );
326  bool pt5_inside = plane.extend( pt5 );
327  ++nb, nbok += pt5_inside == true ? 1 : 0;
328  trace.info() << "(" << nbok << "/" << nb << ") add " << pt5
329  << " Plane=" << plane << std::endl;
330 
331  NaivePlane plane2;
332  plane2.init( 2, 100, 1, 1 );
333  plane2.extend( Point( 10, 0, 0 ) );
334  plane2.extend( Point( 0, 8, 0 ) );
335  plane2.extend( Point( 0, 0, 6 ) );
336  trace.info() << "(" << nbok << "/" << nb << ") "
337  << " Plane2=" << plane2 << std::endl;
338 
339  ++nb, nbok += checkPlane<Integer,NaivePlane>( 11, 5, 19, 20, 100, 100 ) ? 1 : 0;
340  trace.info() << "(" << nbok << "/" << nb
341  << ") checkPlane<Integer,NaivePlane>( 11, 5, 19, 20, 100, 100 )"
342  << std::endl;
343 
344  ++nb, nbok += checkGenericPlane<Integer,GenericNaivePlane>( 11, 5, 19, 20, 100, 100 ) ? 1 : 0;
345  trace.info() << "(" << nbok << "/" << nb
346  << ") checkGenericPlane<Integer,GenericNaivePlane>( 11, 5, 19, 20, 100, 100 )"
347  << std::endl;
348  ++nb, nbok += checkGenericPlane<Integer,GenericNaivePlane>( 17, 33, 7, 10, 100, 100 ) ? 1 : 0;
349  trace.info() << "(" << nbok << "/" << nb
350  << ") checkGenericPlane<Integer,GenericNaivePlane>( 17, 33, 7, 10, 100, 100 )"
351  << std::endl;
352  ++nb, nbok += checkPlane<Integer,NaivePlane>( 15, 8, 13, 15, 100, 100 ) ? 1 : 0;
353  trace.info() << "(" << nbok << "/" << nb
354  << ") checkPlane<Integer,NaivePlane>( 15, 8, 13, 15, 100, 100 )"
355  << std::endl;
356  ++nb, nbok += checkGenericPlane<Integer,GenericNaivePlane>( 15, 8, 13, 15, 100, 100 ) ? 1 : 0;
357  trace.info() << "(" << nbok << "/" << nb
358  << ") checkGenericPlane<Integer,GenericNaivePlane>( 15, 8, 13, 15, 100, 100 )"
359  << std::endl;
360  trace.endBlock();
361  return nbok == nb;
362 }
363 
364 template <typename NaivePlane>
365 bool
366 checkManyPlanes( unsigned int diameter,
367  unsigned int nbplanes,
368  unsigned int nbpoints )
369 {
370  unsigned int nbok = 0;
371  unsigned int nb = 0;
372  typedef typename NaivePlane::InternalInteger Integer;
373  stringstream ss (stringstream::out);
374  ss << "Testing block: Diameter is " << diameter << ". Check " << nbplanes << " planes with " << nbpoints << " points each.";
375  trace.beginBlock ( ss.str() );
376  ++nb, nbok += checkPlanes<Integer,NaivePlane>( nbplanes, diameter, nbpoints ) ? 1 : 0;
377  trace.info() << "(" << nbok << "/" << nb
378  << ") checkPlanes<Integer,NaivePlane>()"
379  << std::endl;
380  trace.endBlock();
381  return nbok == nb;
382 }
383 
387 template <typename NaivePlane>
388 unsigned int maxDiameter( unsigned int min, unsigned int max )
389 {
390  while ( min < max )
391  {
392  unsigned int middle = (min+max)/2;
393  bool ok = checkManyPlanes<NaivePlane>( middle, 2, 2000 );
394  if ( ok ) min = middle+1;
395  else max = middle;
396  }
397  return min-1;
398 }
400 // Standard services - public :
401 
402 int main( int /*argc*/, char** /*argv*/ )
403 {
404  using namespace Z3i;
405  // Max diameter is ~20 for int32_t, ~500 for int64_t, any with BigInteger.
406  trace.beginBlock ( "Testing class COBANaivePlane" );
407  bool res = true
408  && testCOBANaivePlane()
409  && checkManyPlanes<COBANaivePlane<Z3, DGtal::int32_t> >( 20, 100, 200 )
410  && checkManyPlanes<COBANaivePlane<Z3, DGtal::int64_t> >( 500, 100, 200 )
411  && checkManyPlanes<COBANaivePlane<Z3, DGtal::BigInteger> >( 10000, 10, 200 );
412  trace.emphase() << ( res ? "Passed." : "Error." ) << endl;
413  trace.endBlock();
414  // trace.beginBlock ( "Max diameter for COBANaivePlane<Z3, int32_t>" );
415  // unsigned int maxd = maxDiameter<COBANaivePlane<Z3, DGtal::int32_t> >( 10, 1000 );
416  // trace.emphase() << maxd << endl;
417  // trace.endBlock();
418  // trace.beginBlock ( "Max diameter for COBANaivePlane<Z3, int64_t>" );
419  // unsigned int maxd2 = maxDiameter<COBANaivePlane<Z3, DGtal::int32_t> >( 100, 100000 );
420  // trace.emphase() << maxd2 << endl;
421  // trace.endBlock();
422  return res ? 0 : 1;
423 }
424 // //
Generated on Wed Dec 19 2012 19:10:37 for DGtal by   doxygen 1.8.1.1