DGtal 1.4.0
Loading...
Searching...
No Matches
testCOBANaivePlaneComputer.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/COBANaivePlaneComputer.h"
38#include "DGtal/geometry/surfaces/COBAGenericNaivePlaneComputer.h"
40
41using namespace std;
42using namespace DGtal;
43using namespace DGtal::concepts;
44
46// Functions for testing class COBANaivePlaneComputer.
48
49template <typename Integer>
50Integer getRandomInteger( const Integer & first, const Integer & after_last )
51{
52 Integer r = (Integer) rand();
53 return ( r % (after_last - first) ) + first;
54}
55
59template <typename Integer, typename NaivePlaneComputer>
60bool
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;
80 plane.init( axis, diameter, 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 break;
105 }
106 if ( ! ok_ext )
107 {
108 std::cerr << "[ERROR] p=" << p << " was NOT extendable IN plane=" << plane << std::endl;
109 break;
110 }
111 // else
112 // std::cerr << "[OK] p=" << p << " IN plane=" << plane << std::endl;
113 }
114
115 // Checks that points outside the naive plane are correctly recognized as outliers.
116 while ( nb != (nbtries * 11 ) / 10 )
117 {
118 p[ 0 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
119 p[ 1 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
120 p[ 2 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
121 x = (Integer) p[ 0 ];
122 y = (Integer) p[ 1 ];
123 z = (Integer) p[ 2 ];
124 switch ( axis ) {
125 case 0: p[ 0 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - b * y - c * z, a ) ); break;
126 case 1: p[ 1 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - c * z, b ) ); break;
127 case 2: p[ 2 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - b * y, c ) ); break;
128 }
129 PointInteger tmp = getRandomInteger<PointInteger>( 2, 5 )
130 * (2*getRandomInteger<PointInteger>( 0, 2 ) - 1 );
131 p[ axis ] += tmp;
132 bool ok_ext = ! plane.isExtendable( p ); // should *not* be ok
133 bool ok = ! plane.extend( p ); // should *not* be ok
134 ++nb; nbok += ok ? 1 : 0;
135 ++nb; nbok += ok_ext ? 1 : 0;
136 if ( ! ok )
137 {
138 std::cerr << "[ERROR] p=" << p << " IN plane=" << plane << std::endl;
139 break;
140 }
141 if ( ! ok_ext )
142 {
143 std::cerr << "[ERROR] p=" << p << " was extendable IN plane=" << plane << std::endl;
144 break;
145 }
146 // else
147 // std::cerr << "[OK] p=" << p << " IN plane=" << plane << std::endl;
148 }
149 return nb == nbok;
150}
151
152
156template <typename Integer, typename GenericNaivePlaneComputer>
157bool
159 int diameter, unsigned int nbtries )
160{
161 typedef typename GenericNaivePlaneComputer::Point Point;
162 typedef typename Point::Component PointInteger;
164 Integer absA = ic.abs( a );
165 Integer absB = ic.abs( b );
166 Integer absC = ic.abs( c );
167 Integer x, y, z;
168 Dimension axis;
169 if ( ( absA >= absB ) && ( absA >= absC ) )
170 axis = 0;
171 else if ( ( absB >= absA ) && ( absB >= absC ) )
172 axis = 1;
173 else
174 axis = 2;
175 Point p;
176 GenericNaivePlaneComputer plane;
177 plane.init( diameter, 1, 1 );
178 // Checks that points within the naive plane are correctly recognized.
179 unsigned int nb = 0;
180 unsigned int nbok = 0;
181 while ( nb != nbtries )
182 {
183 p[ 0 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
184 p[ 1 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
185 p[ 2 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
186 x = (Integer) p[ 0 ];
187 y = (Integer) p[ 1 ];
188 z = (Integer) p[ 2 ];
189 switch ( axis ) {
190 case 0: p[ 0 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - b * y - c * z, a ) ); break;
191 case 1: p[ 1 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - c * z, b ) ); break;
192 case 2: p[ 2 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - b * y, c ) ); break;
193 }
194 bool ok_ext = plane.isExtendable( p ); // should be ok
195 bool ok = plane.extend( p ); // should be ok
196 ++nb; nbok += ok_ext ? 1 : 0;
197 ++nb; nbok += ok ? 1 : 0;
198 if ( ! ok )
199 {
200 std::cerr << "[ERROR] p=" << p << " NOT IN plane=" << plane << std::endl;
201 break;
202 }
203 if ( ! ok_ext )
204 {
205 std::cerr << "[ERROR] p=" << p << " was NOT extendable IN plane=" << plane << std::endl;
206 break;
207 }
208 // else
209 // std::cerr << "[OK] p=" << p << " IN plane=" << plane << std::endl;
210 }
211
212 // Checks that points outside the naive plane are correctly recognized as outliers.
213 while ( nb != (nbtries * 11 ) / 10 )
214 {
215 p[ 0 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
216 p[ 1 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
217 p[ 2 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
218 x = (Integer) p[ 0 ];
219 y = (Integer) p[ 1 ];
220 z = (Integer) p[ 2 ];
221 switch ( axis ) {
222 case 0: p[ 0 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - b * y - c * z, a ) ); break;
223 case 1: p[ 1 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - c * z, b ) ); break;
224 case 2: p[ 2 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - b * y, c ) ); break;
225 }
226 PointInteger tmp = getRandomInteger<PointInteger>( 2, 5 )
227 * (2*getRandomInteger<PointInteger>( 0, 2 ) - 1 );
228 p[ axis ] += tmp;
229 bool ok_ext = ! plane.isExtendable( p ); // should *not* be ok
230 bool ok = ! plane.extend( p ); // should *not* be ok
231 ++nb; nbok += ok ? 1 : 0;
232 ++nb; nbok += ok_ext ? 1 : 0;
233 if ( ! ok )
234 {
235 std::cerr << "[ERROR] p=" << p << " IN plane=" << plane << std::endl;
236 break;
237 }
238 if ( ! ok_ext )
239 {
240 std::cerr << "[ERROR] p=" << p << " was extendable IN plane=" << plane << std::endl;
241 break;
242 }
243 // else
244 // std::cerr << "[OK] p=" << p << " IN plane=" << plane << std::endl;
245 }
246 std::cerr << "plane = " << plane << std::endl;
247 return nb == nbok;
248}
249
250
251template <typename Integer, typename NaivePlaneComputer>
252bool
253checkPlanes( unsigned int nbplanes, int diameter, unsigned int nbtries )
254{
255 //using namespace Z3i;
256 //typedef COBANaivePlaneComputer<Z3, Integer> NaivePlaneComputer;
257 unsigned int nb = 0;
258 unsigned int nbok = 0;
259 for ( unsigned int nbp = 0; nbp < nbplanes; ++nbp )
260 {
261 Integer a = getRandomInteger<Integer>( (Integer) 0, (Integer) diameter / 2 );
262 Integer b = getRandomInteger<Integer>( (Integer) 0, (Integer) diameter / 2 );
263 Integer c = getRandomInteger<Integer>( (Integer) 0, (Integer) diameter / 2 );
264 Integer d = getRandomInteger<Integer>( (Integer) 0, (Integer) diameter / 2 );
265 if ( ( a != 0 ) || ( b != 0 ) || ( c != 0 ) )
266 {
267 ++nb; nbok += checkPlane<Integer, NaivePlaneComputer>( a, b, c, d, diameter, nbtries ) ? 1 : 0;
268 if ( nb != nbok )
269 {
270 std::cerr << "[ERROR] for plane " << a << " * x + "
271 << b << " * y + " << c << " * z = " << d << std::endl;
272 break;
273 }
274 }
275 }
276 return nb == nbok;
277}
278
284{
285 unsigned int nbok = 0;
286 unsigned int nb = 0;
287 using namespace Z3i;
289 typedef COBAGenericNaivePlaneComputer<Z3, BigInteger> GenericNaivePlaneComputer;
290
291 BOOST_CONCEPT_ASSERT(( CAdditivePrimitiveComputer< NaivePlaneComputer > ));
293 BOOST_CONCEPT_ASSERT(( boost::ForwardContainer< NaivePlaneComputer > ));
295 BOOST_CONCEPT_ASSERT(( CPointPredicate< NaivePlaneComputer::Primitive > ));
297
298 trace.beginBlock ( "Testing block: COBANaivePlaneComputer instantiation." );
299 NaivePlaneComputer plane;
300 Point pt0( 0, 0, 0 );
301 plane.init( 2, 100, 3, 2 );
302 bool pt0_inside = plane.extend( pt0 );
303 FATAL_ERROR(pt0_inside);
304
305 trace.info() << "(" << nbok << "/" << nb << ") Plane=" << plane
306 << std::endl;
307 Point pt1( Point( 8, 1, 3 ) );
308 bool pt1_inside = plane.extend( pt1 );
309 ++nb; nbok += pt1_inside == true ? 1 : 0;
310 trace.info() << "(" << nbok << "/" << nb << ") add " << pt1
311 << " Plane=" << plane << std::endl;
312 Point pt2( Point( 2, 7, 1 ) );
313 bool pt2_inside = plane.extend( pt2 );
314 ++nb; nbok += pt2_inside == true ? 1 : 0;
315 trace.info() << "(" << nbok << "/" << nb << ") add " << pt2
316 << " Plane=" << plane << std::endl;
317
318 Point pt3( Point( 0, 5, 17 ) );
319 bool pt3_inside = plane.extend( pt3 );
320 ++nb; nbok += pt3_inside == false ? 1 : 0;
321 trace.info() << "(" << nbok << "/" << nb << ") add " << pt3
322 << " Plane=" << plane << std::endl;
323
324 Point pt4( Point( -10, -10, 10 ) );
325 bool pt4_inside = plane.extend( pt4 );
326 ++nb; nbok += pt4_inside == false ? 1 : 0;
327 trace.info() << "(" << nbok << "/" << nb << ") add " << pt4
328 << " Plane=" << plane << std::endl;
329
330 Point pt5 = pt0 + pt1 + pt2 + Point( 0, 0, 2 );
331 bool pt5_inside = plane.extend( pt5 );
332 ++nb; nbok += pt5_inside == true ? 1 : 0;
333 trace.info() << "(" << nbok << "/" << nb << ") add " << pt5
334 << " Plane=" << plane << std::endl;
335
336 NaivePlaneComputer plane2;
337 plane2.init( 2, 100, 1, 1 );
338 plane2.extend( Point( 10, 0, 0 ) );
339 plane2.extend( Point( 0, 8, 0 ) );
340 plane2.extend( Point( 0, 0, 6 ) );
341 trace.info() << "(" << nbok << "/" << nb << ") "
342 << " Plane2=" << plane2 << std::endl;
343
344 ++nb; nbok += checkPlane<Integer,NaivePlaneComputer>( 11, 5, 19, 20, 100, 100 ) ? 1 : 0;
345 trace.info() << "(" << nbok << "/" << nb
346 << ") checkPlane<Integer,NaivePlaneComputer>( 11, 5, 19, 20, 100, 100 )"
347 << std::endl;
348
349 ++nb; nbok += checkGenericPlane<Integer,GenericNaivePlaneComputer>( 11, 5, 19, 20, 100, 100 ) ? 1 : 0;
350 trace.info() << "(" << nbok << "/" << nb
351 << ") checkGenericPlane<Integer,GenericNaivePlaneComputer>( 11, 5, 19, 20, 100, 100 )"
352 << std::endl;
353 ++nb; nbok += checkGenericPlane<Integer,GenericNaivePlaneComputer>( 17, 33, 7, 10, 100, 100 ) ? 1 : 0;
354 trace.info() << "(" << nbok << "/" << nb
355 << ") checkGenericPlane<Integer,GenericNaivePlaneComputer>( 17, 33, 7, 10, 100, 100 )"
356 << std::endl;
357 ++nb; nbok += checkPlane<Integer,NaivePlaneComputer>( 15, 8, 13, 15, 100, 100 ) ? 1 : 0;
358 trace.info() << "(" << nbok << "/" << nb
359 << ") checkPlane<Integer,NaivePlaneComputer>( 15, 8, 13, 15, 100, 100 )"
360 << std::endl;
361 ++nb; nbok += checkGenericPlane<Integer,GenericNaivePlaneComputer>( 15, 8, 13, 15, 100, 100 ) ? 1 : 0;
362 trace.info() << "(" << nbok << "/" << nb
363 << ") checkGenericPlane<Integer,GenericNaivePlaneComputer>( 15, 8, 13, 15, 100, 100 )"
364 << std::endl;
365 trace.endBlock();
366 return nbok == nb;
367}
368
369template <typename NaivePlaneComputer>
370bool
371checkManyPlanes( unsigned int diameter,
372 unsigned int nbplanes,
373 unsigned int nbpoints )
374{
375 unsigned int nbok = 0;
376 unsigned int nb = 0;
378 stringstream ss (stringstream::out);
379 ss << "Testing block: Diameter is " << diameter << ". Check " << nbplanes << " planes with " << nbpoints << " points each.";
380 trace.beginBlock ( ss.str() );
381 ++nb; nbok += checkPlanes<Integer,NaivePlaneComputer>( nbplanes, diameter, nbpoints ) ? 1 : 0;
382 trace.info() << "(" << nbok << "/" << nb
383 << ") checkPlanes<Integer,NaivePlaneComputer>()"
384 << std::endl;
385 trace.endBlock();
386 return nbok == nb;
387}
388
392template <typename NaivePlaneComputer>
393unsigned int maxDiameter( unsigned int min, unsigned int max )
394{
395 while ( min < max )
396 {
397 unsigned int middle = (min+max)/2;
398 bool ok = checkManyPlanes<NaivePlaneComputer>( middle, 2, 2000 );
399 if ( ok ) min = middle+1;
400 else max = middle;
401 }
402 return min-1;
403}
404
405template <typename GenericNaivePlaneComputer>
406bool
407checkExtendWithManyPoints( unsigned int diameter,
408 unsigned int nbplanes,
409 unsigned int nbpoints )
410{
411 unsigned int nbok = 0;
412 unsigned int nb = 0;
413 typedef typename GenericNaivePlaneComputer::InternalInteger Integer;
414 typedef typename GenericNaivePlaneComputer::Point Point;
415 typedef typename Point::Coordinate PointInteger;
417
418 trace.beginBlock( "checkExtendWithManyPoints" );
419 for ( unsigned int j = 0; j < nbplanes; ++j )
420 {
421 Integer a = getRandomInteger<Integer>( (Integer) 0, (Integer) diameter / 2 );
422 Integer b = getRandomInteger<Integer>( (Integer) 0, (Integer) diameter / 2 );
423 Integer c = getRandomInteger<Integer>( (Integer) 1, (Integer) diameter / 2 );
424 Integer d = getRandomInteger<Integer>( (Integer) 0, (Integer) diameter / 2 );
425 GenericNaivePlaneComputer plane;
426 Dimension axis;
427 if ( ( a >= b ) && ( a >= c ) ) axis = 0;
428 else if ( ( b >= a ) && ( b >= c ) ) axis = 1;
429 else axis = 2;
430 plane.init( diameter, 1, 1 );
431
432 std::vector<Point> pts;
433 for ( unsigned int i = 0; i < nbpoints; ++i )
434 {
435 Point p;
436 p[ 0 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
437 p[ 1 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
438 p[ 2 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
439 Integer x = (Integer) p[ 0 ];
440 Integer y = (Integer) p[ 1 ];
441 Integer z = (Integer) p[ 2 ];
442 switch( axis ) {
443 case 0: p[ 0 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - b * y - c * z, a ) ); break;
444 case 1: p[ 1 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - c * z, b ) ); break;
445 case 2: p[ 2 ] = NumberTraits<Integer>::castToInt64_t( ic.ceilDiv( d - a * x - b * y, c ) ); break;
446 }
447 pts.push_back( p );
448 }
449 ++nb; nbok += plane.isExtendable( pts.begin(), pts.end() ); // should be ok
450 trace.info() << "(" << nbok << "/" << nb
451 << ") plane.isExtendable( pts.begin(), pts.end() )"
452 << std::endl;
453 Point & any0 = pts[ getRandomInteger<int>( 0, pts.size() ) ];
454 pts.push_back( any0 + Point(1,0,0) );
455 Point & any1 = pts[ getRandomInteger<int>( 0, pts.size() ) ];
456 pts.push_back( any1 + Point(0,1,0) );
457 Point & any2 = pts[ getRandomInteger<int>( 0, pts.size() ) ];
458 pts.push_back( any2 + Point(0,0,1) );
459 bool check = ! plane.isExtendable( pts.begin(), pts.end() ); // should not be ok
460 ++nb; nbok += check ? 1 : 0;
461 trace.info() << "(" << nbok << "/" << nb
462 << ") ! plane.isExtendable( pts.begin(), pts.end() )"
463 << std::endl;
464 if ( ! check )
465 trace.warning() << plane << " last=" << pts.back() << std::endl
466 << "a=" << a << " b=" << b << " c=" << c << " d=" << d << std::endl;
467 ++nb; nbok += plane.extend( pts.begin(), pts.end() - 3 ); // should be ok
468 trace.info() << "(" << nbok << "/" << nb
469 << ") plane.extend( pts.begin(), pts.end() - 3)"
470 << std::endl;
471 ++nb; nbok += ! plane.extend( pts.end() - 3, pts.end() ); // should not be ok
472 trace.info() << "(" << nbok << "/" << nb
473 << ") ! plane.extend( pts.end() - 3, pts.end() )"
474 << std::endl;
475 }
476 trace.endBlock();
477 return nb == nbok;
478}
479
481// Standard services - public :
482
483int main( int /*argc*/, char** /*argv*/ )
484{
485 using namespace Z3i;
486
487 // Max diameter is ~20 for int32_t, ~500 for int64_t, any with BigInteger.
488 trace.beginBlock ( "Testing class COBANaivePlaneComputer" );
489 bool res = true
495
496 trace.emphase() << ( res ? "Passed." : "Error." ) << endl;
497 trace.endBlock();
498 // trace.beginBlock ( "Max diameter for COBANaivePlaneComputer<Z3, int32_t>" );
499 // unsigned int maxd = maxDiameter<COBANaivePlaneComputer<Z3, DGtal::int32_t> >( 10, 1000 );
500 // trace.emphase() << maxd << endl;
501 // trace.endBlock();
502 // trace.beginBlock ( "Max diameter for COBANaivePlaneComputer<Z3, int64_t>" );
503 // unsigned int maxd2 = maxDiameter<COBANaivePlaneComputer<Z3, DGtal::int32_t> >( 100, 100000 );
504 // trace.emphase() << maxd2 << endl;
505 // trace.endBlock();
506 return res ? 0 : 1;
507}
508// //
Aim: A class that recognizes pieces of digital planes of given axis width. When the width is 1,...
void init(Dimension axis, InternalInteger diameter, InternalInteger widthNumerator=NumberTraits< InternalInteger >::ONE, InternalInteger widthDenominator=NumberTraits< InternalInteger >::ONE)
bool isExtendable(const Point &p) const
bool extend(const Point &p)
Aim: This class gathers several types and methods to make computation with integers.
static Integer abs(IntegerParamType a)
Integer ceilDiv(IntegerParamType na, IntegerParamType nb) const
void beginBlock(const std::string &keyword="")
std::ostream & warning()
std::ostream & emphase()
std::ostream & info()
double endBlock()
COBANaivePlaneComputer< Z3, InternalInteger > NaivePlaneComputer
Aim: Gathers several functions useful for concept checks.
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::uint32_t Dimension
Definition Common.h:136
Trace trace
Definition Common.h:153
STL namespace.
Aim: The traits class for all models of Cinteger.
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
int max(int a, int b)
int main()
Definition testBits.cpp:56
unsigned int maxDiameter(unsigned int min, unsigned int max)
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)
Integer getRandomInteger(const Integer &first, const Integer &after_last)
bool checkExtendWithManyPoints(unsigned int diameter, unsigned int nbplanes, unsigned int nbpoints)
bool testCOBANaivePlaneComputer()
bool checkManyPlanes(unsigned int diameter, unsigned int nbplanes, unsigned int nbpoints)
MyPointD Point