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"
42 using namespace DGtal;
49 template <
typename Integer>
53 return ( r % (after_last - first) ) + first;
59 template <
typename Integer,
typename NaivePlaneComputer>
62 int diameter,
unsigned int nbtries )
65 typedef typename Point::Component PointInteger;
72 if ( ( absA >= absB ) && ( absA >= absC ) )
74 else if ( ( absB >= absA ) && ( absB >= absC ) )
80 plane.
init( axis, 1, 1 );
83 unsigned int nbok = 0;
84 while ( nb != nbtries )
86 p[ 0 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
87 p[ 1 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
88 p[ 2 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
98 bool ok = plane.
extend( p );
99 ++nb; nbok += ok_ext ? 1 : 0;
100 ++nb; nbok += ok ? 1 : 0;
103 std::cerr <<
"[ERROR] p=" << p <<
" NOT IN plane=" << plane << std::endl;
106 std::cerr <<
" " << *it;
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])
118 std::cerr <<
"[ERROR] p=" << p <<
" was NOT extendable IN plane=" << plane << std::endl;
126 while ( nb != (nbtries * 11 ) / 10 )
128 p[ 0 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
129 p[ 1 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
130 p[ 2 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
139 PointInteger tmp = getRandomInteger<PointInteger>( 2, 5 )
140 * (2*getRandomInteger<PointInteger>( 0, 2 ) - 1 );
143 bool ok = ! plane.
extend( p );
144 ++nb; nbok += ok ? 1 : 0;
145 ++nb; nbok += ok_ext ? 1 : 0;
148 std::cerr <<
"[ERROR] p=" << p <<
" IN plane=" << plane << std::endl;
153 std::cerr <<
"[ERROR] p=" << p <<
" was extendable IN plane=" << plane << std::endl;
165 template <
typename Integer,
typename NaivePlaneComputer>
168 int diameter,
unsigned int nbtries )
171 typedef typename Point::Component PointInteger;
178 if ( ( absA >= absB ) && ( absA >= absC ) )
180 else if ( ( absB >= absA ) && ( absB >= absC ) )
186 plane.
init( axis, 1, 1 );
189 unsigned int nbok = 0;
190 while ( nb < nbtries )
192 std::vector<Point> points( 5 );
193 for (
unsigned int i = 0; i < 5; ++i )
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 );
204 ( ic.
ceilDiv( d - b * y - c * z, a ) );
break;
206 ( ic.
ceilDiv( d - a * x - c * z, b ) );
break;
208 ( ic.
ceilDiv( d - a * x - b * y, c ) );
break;
211 bool ok_ext = plane.
isExtendable( points.begin(), points.end() );
212 bool ok = plane.
extend( points.begin(), points.end() );
213 ++nb; nbok += ok_ext ? 1 : 0;
214 ++nb; nbok += ok ? 1 : 0;
217 std::cerr <<
"[ERROR] p=" << points[ 0 ] <<
" NOT IN plane=" << plane << std::endl;
220 std::cerr <<
" " << *it;
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])
232 std::cerr <<
"[ERROR] p=" << p <<
" was NOT extendable IN plane=" << plane << std::endl;
240 while ( nb < (nbtries * 11 ) / 10 )
242 p[ 0 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
243 p[ 1 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
244 p[ 2 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
253 PointInteger tmp = getRandomInteger<PointInteger>( 2, 5 )
254 * (2*getRandomInteger<PointInteger>( 0, 2 ) - 1 );
257 bool ok = ! plane.
extend( p );
258 ++nb; nbok += ok ? 1 : 0;
259 ++nb; nbok += ok_ext ? 1 : 0;
262 std::cerr <<
"[ERROR] p=" << p <<
" IN plane=" << plane << std::endl;
267 std::cerr <<
"[ERROR] p=" << p <<
" was extendable IN plane=" << plane << std::endl;
281 template <
typename Integer,
typename GenericNaivePlaneComputer>
284 int diameter,
unsigned int nbtries )
287 typedef typename Point::Component PointInteger;
294 if ( ( absA >= absB ) && ( absA >= absC ) )
296 else if ( ( absB >= absA ) && ( absB >= absC ) )
301 GenericNaivePlaneComputer plane;
305 unsigned int nbok = 0;
306 while ( nb != nbtries )
308 p[ 0 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
309 p[ 1 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
310 p[ 2 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
319 bool ok_ext = plane.isExtendable( p );
320 bool ok = plane.extend( p );
321 ++nb; nbok += ok_ext ? 1 : 0;
322 ++nb; nbok += ok ? 1 : 0;
325 std::cerr <<
"[ERROR] p=" << p <<
" NOT IN plane=" << plane << std::endl;
328 std::cerr <<
" " << *it;
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])
340 std::cerr <<
"[ERROR] p=" << p <<
" was NOT extendable IN plane=" << plane << std::endl;
348 while ( nb != (nbtries * 11 ) / 10 )
350 p[ 0 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
351 p[ 1 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
352 p[ 2 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
361 PointInteger tmp = getRandomInteger<PointInteger>( 2, 5 )
362 * (2*getRandomInteger<PointInteger>( 0, 2 ) - 1 );
364 bool ok_ext = ! plane.isExtendable( p );
365 bool ok = ! plane.extend( p );
366 ++nb; nbok += ok ? 1 : 0;
367 ++nb; nbok += ok_ext ? 1 : 0;
370 std::cerr <<
"[ERROR] p=" << p <<
" IN plane=" << plane << std::endl;
375 std::cerr <<
"[ERROR] p=" << p <<
" was extendable IN plane=" << plane << std::endl;
381 std::cerr <<
"plane = " << plane << std::endl;
386 template <
typename Integer,
typename NaivePlaneComputer>
388 checkPlanes(
unsigned int nbplanes,
int diameter,
unsigned int nbtries )
393 unsigned int nbok = 0;
394 for (
unsigned int nbp = 0; nbp < nbplanes; ++nbp )
400 if ( ( a != 0 ) || ( b != 0 ) || ( c != 0 ) )
402 ++nb; nbok += checkPlane<Integer, NaivePlaneComputer>( a, b, c, d, diameter, nbtries ) ? 1 : 0;
405 std::cerr <<
"[ERROR] (Simple extension) for plane " << a <<
" * x + "
406 << b <<
" * y + " << c <<
" * z = " << d << std::endl;
409 ++nb; nbok += checkPlaneGroupExtension<Integer, NaivePlaneComputer>( a, b, c, d, diameter, nbtries ) ? 1 : 0;
412 std::cerr <<
"[ERROR] (Group extension) for plane " << a <<
" * x + "
413 << b <<
" * y + " << c <<
" * z = " << d << std::endl;
424 template <
typename Integer,
typename NaivePlaneComputer>
427 int diameter,
unsigned int nbtries )
430 typedef typename NaivePlaneComputer::InternalScalar InternalScalar;
437 if ( ( absA >= absB ) && ( absA >= absC ) )
439 else if ( ( absB >= absA ) && ( absB >= absC ) )
445 unsigned int nbok = 0;
446 std::vector<Point> points( nbtries );
447 for (
unsigned int i = 0; i != nbtries; ++i )
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 );
464 << d <<
" <= " << a <<
"*x"
467 <<
" <= d + max(|a|,|b|,|c|)"
469 trace.
info() <<
"- " << points.size() <<
" points tested in diameter " << diameter
472 for (
unsigned int i = 0; i < 3; ++i )
474 std::pair<InternalScalar, InternalScalar> width
475 = NaivePlaneComputer::computeAxisWidth( i, points.begin(), points.end() );
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;
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;
492 template <
typename Integer,
typename NaivePlaneComputer>
494 checkWidths(
unsigned int nbplanes,
int diameter,
unsigned int nbtries )
499 unsigned int nbok = 0;
500 for (
unsigned int nbp = 0; nbp < nbplanes; ++nbp )
506 if ( ( a != 0 ) || ( b != 0 ) || ( c != 0 ) )
508 ++nb; nbok += checkWidth<Integer, NaivePlaneComputer>( a, b, c, d, diameter, nbtries ) ? 1 : 0;
511 std::cerr <<
"[ERROR] (checkWidth) for plane " << a <<
" * x + "
512 << b <<
" * y + " << c <<
" * z = " << d << std::endl;
527 unsigned int nbok = 0;
542 trace.
beginBlock (
"Testing block: ChordNaivePlaneComputer instantiation." );
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
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;
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;
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;
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;
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;
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;
586 plane2.
init( 2, 1, 1 );
590 trace.
info() <<
"(" << nbok <<
"/" << nb <<
") "
591 <<
" Plane2=" << plane2 << std::endl;
593 ++nb; nbok += checkPlane<Integer,NaivePlaneComputer>( 11, 5, 19, 20, 100, 100 ) ? 1 : 0;
595 <<
") checkPlane<Integer,NaivePlaneComputer>( 11, 5, 19, 20, 100, 100 )"
598 ++nb; nbok += checkGenericPlane<Integer,GenericNaivePlaneComputer>( 11, 5, 19, 20, 100, 100 ) ? 1 : 0;
600 <<
") checkGenericPlane<Integer,GenericNaivePlaneComputer>( 11, 5, 19, 20, 100, 100 )"
602 ++nb; nbok += checkGenericPlane<Integer,GenericNaivePlaneComputer>( 17, 33, 7, 10, 100, 100 ) ? 1 : 0;
604 <<
") checkGenericPlane<Integer,GenericNaivePlaneComputer>( 17, 33, 7, 10, 100, 100 )"
606 ++nb; nbok += checkPlane<Integer,NaivePlaneComputer>( 15, 8, 13, 15, 100, 100 ) ? 1 : 0;
608 <<
") checkPlane<Integer,NaivePlaneComputer>( 15, 8, 13, 15, 100, 100 )"
610 ++nb; nbok += checkGenericPlane<Integer,GenericNaivePlaneComputer>( 15, 8, 13, 15, 100, 100 ) ? 1 : 0;
612 <<
") checkGenericPlane<Integer,GenericNaivePlaneComputer>( 15, 8, 13, 15, 100, 100 )"
617 trace.
beginBlock (
"Testing block: ChordNaivePlaneComputer vertical instantiation." );
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
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
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
636 bool pppt3_inside = ppplane.
extend( pppt3 );
637 ++nb; nbok += pppt3_inside ==
true ? 1 : 0;
638 trace.
info() <<
"(" << nbok <<
"/" << nb <<
") Plane=" << ppplane
641 bool pppt4_inside = ppplane.
extend( pt4 );
642 ++nb; nbok += pppt4_inside ==
true ? 1 : 0;
643 trace.
info() <<
"(" << nbok <<
"/" << nb <<
") Plane=" << ppplane
649 trace.
beginBlock (
"Testing block: ChordNaivePlaneComputer vertical instantiation 2." );
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
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
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
673 template <
typename NaivePlaneComputer>
676 unsigned int nbplanes,
677 unsigned int nbpoints )
679 unsigned int nbok = 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.";
685 ++nb; nbok += checkPlanes<Scalar,NaivePlaneComputer>( nbplanes, diameter, nbpoints ) ? 1 : 0;
687 <<
") checkPlanes<Scalar,NaivePlaneComputer>()"
693 template <
typename GenericNaivePlaneComputer>
696 unsigned int nbplanes,
697 unsigned int nbpoints )
699 unsigned int nbok = 0;
701 typedef typename GenericNaivePlaneComputer::InternalScalar
Integer;
703 typedef typename Point::Coordinate PointInteger;
707 for (
unsigned int j = 0; j < nbplanes; ++j )
713 GenericNaivePlaneComputer plane;
715 if ( ( a >= b ) && ( a >= c ) ) axis = 0;
716 else if ( ( b >= a ) && ( b >= c ) ) axis = 1;
720 std::vector<Point> pts;
721 for (
unsigned int i = 0; i < nbpoints; ++i )
724 p[ 0 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
725 p[ 1 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
726 p[ 2 ] = getRandomInteger<PointInteger>( -diameter+1, diameter );
737 ++nb; nbok += plane.isExtendable( pts.begin(), pts.end() );
739 <<
") plane.isExtendable( pts.begin(), pts.end() )"
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() );
748 ++nb; nbok += check ? 1 : 0;
750 <<
") ! plane.isExtendable( pts.begin(), pts.end() )"
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 );
757 <<
") plane.extend( pts.begin(), pts.end() - 3)"
759 ++nb; nbok += ! plane.extend( pts.end() - 3, pts.end() );
761 <<
") ! plane.extend( pts.end() - 3, pts.end() )"
780 && checkManyPlanes<ChordNaivePlaneComputer<Z3i::Space, Z3i::Point, DGtal::int32_t> >( 4, 100, 200 )
782 && checkManyPlanes<ChordNaivePlaneComputer<Z3i::Space, Z3i::Point, DGtal::int32_t> >( 20, 100, 200 )
784 && checkManyPlanes<ChordNaivePlaneComputer<Z3i::Space, Z3i::Point, DGtal::int64_t> >( 2000, 100, 200 )
786 && checkExtendWithManyPoints<ChordGenericNaivePlaneComputer<Z3i::Space, Z3i::Point, DGtal::int64_t> >( 100, 100, 200 );
788 trace.
emphase() << ( res ?
"Passed." :
"Error." ) << endl;