Loading [MathJax]/extensions/TeX/AMSsymbols.js
DGtal 2.0.0
testArithmeticalDSSConvexHull.cpp
Go to the documentation of this file.
1
16
29
31#include <iostream>
32#include <sstream>
33#include "DGtal/base/Common.h"
34#include "DGtal/geometry/curves/ArithmeticalDSL.h"
35#include "DGtal/geometry/curves/ArithmeticalDSSConvexHull.h"
36#include "DGtal/geometry/curves/ArithmeticalDSS.h"
37
38#include "DGtal/kernel/SpaceND.h"
39#include "DGtal/kernel/NumberTraits.h"
40#include "DGtal/arithmetic/LatticePolytope2D.h"
41
43
44using namespace std;
45using namespace DGtal;
46using namespace functions;
47
49// Functions for testing functions smartCH and reversedSmartCH.
51
53// smartCH
61template <typename DSL>
62bool basicTest(const DSL& aDSL)
63{
64 typedef typename DSL::Point Point;
65 typedef typename DSL::Vector Vector;
66 typedef typename DSL::Integer Integer;
67 typedef typename DSL::Position Position;
68
69 unsigned int nbok = 0;
70 unsigned int nb = 0;
71
72 trace.beginBlock ( "One simple test..." );
73
74 trace.info() << aDSL << std::endl;
75
76 //forward test
77 Position l = (2*aDSL.patternLength());
78 std::vector<Point> lch, uch;
79 Vector v = smartCH( aDSL, Point(0,0), l,
80 std::back_inserter(uch), std::back_inserter(lch) );
81
82 if (v == Vector(aDSL.b(),aDSL.a())) //eq. 7
83 nbok++;
84 nb++;
85 trace.info() << "(" << nbok << "/" << nb << ") " << std::endl;
86 if ( aDSL.remainder(lch.back()) == aDSL.mu()-1 ) //eq. 8
87 nbok++;
88 nb++;
89 trace.info() << "(" << nbok << "/" << nb << ") " << std::endl; //eq. 8
90 if ( aDSL.remainder(uch.back()) == aDSL.mu() )
91 nbok++;
92 nb++;
93 trace.info() << "(" << nbok << "/" << nb << ") " << std::endl;
94
95 //test on the size of the convex hulls
96 double bound = NumberTraits<Integer>::castToDouble(aDSL.patternLength());
97 unsigned int threshold = (unsigned int) std::ceil( std::log(bound) / std::log(1.618) );
98 if ( (lch.size()+uch.size()-1) <= threshold )
99 nbok++;
100 nb++;
101 trace.info() << "(" << nbok << "/" << nb << ") " << std::endl;
102
103 trace.endBlock();
104
105 return (nbok == nb);
106}
107
117template <typename DSL>
118bool basicTest(typename DSL::Coordinate a,
119 typename DSL::Coordinate b)
120{
121 unsigned int nbok = 0;
122 unsigned int nb = 0;
123
124 DSL aDSL(a,b,0);
125 for (typename DSL::Integer mu = 0; ( (-mu < aDSL.omega())&&(nbok == nb) ); --mu)
126 {
127 if (basicTest(DSL(a,b,mu)))
128 nbok++;
129 nb++;
130 }
131
132 return (nbok == nb);
133}
134
139{
140 trace.beginBlock ( "Testing block (without length constraint)..." );
141
142 bool res = true;
143 res = res
155
167
168
175
176 ;
177
178 trace.endBlock();
179
180 return res;
181}
182
183
194template <typename DSL>
195bool comparisonLeftHull(typename DSL::Coordinate a, typename DSL::Coordinate b)
196{
197 ASSERT(a >= 0);
198 ASSERT(b >= 0);
199 ASSERT(a < b);
200
201 typedef typename DSL::Point Point;
202 typedef typename DSL::Vector Vector;
203 typedef typename DSL::Coordinate Coordinate;
204
205 unsigned int nbok = 0;
206 unsigned int nb = 0;
207
208 trace.beginBlock ( "Comparison ..." );
209 DSL inputDSL(a,b,0);
210
211 trace.info() << a << " " << b << std::endl;
212 trace.info() << "testing every mu between 0 and -" << inputDSL.omega() << std::endl;
213
214 for (typename DSL::Integer mu = 0; ( (mu-1 >= -inputDSL.omega())&&(nbok == nb) ); --mu)
215 {
216 trace.info() << "mu=" << mu << ", testing every length between 1 and 2*" << inputDSL.omega() << std::endl;
217 inputDSL = DSL(a,b,mu);
218
219 for (typename DSL::Position l = 1; ( (l <= 2*inputDSL.patternLength())&&(nbok == nb) ); ++l)
220 {
221 //trace.info() << "l=" << l << std::endl;
222
223 //smartCH
224 std::vector<Point> lch, uch;
225 Vector v = smartCH( inputDSL, Point(0,0), l,
226 std::back_inserter(uch), std::back_inserter(lch) );
227
228 // std::copy( lch.begin(), lch.end(), std::ostream_iterator<Point>(std::cout, " ") );
229 // std::cout << std::endl;
230 // std::copy( uch.begin(), uch.end(), std::ostream_iterator<Point>(std::cout, " ") );
231 // std::cout << std::endl;
232
233 Vector shift = -inputDSL.shift();
234 //algorithm of Charrier and Buzer
235 typedef SpaceND<2, Coordinate> Space2;
236 typedef LatticePolytope2D<Space2> CIP;
237 CIP cip;
238 typename CIP::HalfSpace line( typename CIP::Vector(a,-b), mu );
239 typename CIP::HalfSpace line2( typename CIP::Vector(a,-b), mu-1 );
240 //NB: since only closed half-space are used,
241 // we must run the procedure twice, for mu and mu-1
242 // in order to get the lower left hull (not included)
243 // and the upper left hull (included)
244 typename CIP::HalfSpace constraint( typename CIP::Vector(shift[1],-shift[0]), l );
245 std::vector<typename CIP::Point> inch, outch, inch2, outch2;
246 inch.push_back( typename CIP::Point(shift[0],shift[1]) );
247 inch2.push_back( typename CIP::Point(shift[0],shift[1]) );
248 ASSERT( line(inch[0]) );
249 ASSERT( constraint(inch[0]) );
250 outch.push_back( typename CIP::Point(0,0) );
251 outch2.push_back( typename CIP::Point(0,0) );
252 ASSERT( (!line2(outch[0])) );
253 ASSERT( constraint(outch[0]) );
254 typename CIP::Vector vBezout(1,0);
255 cip.getAllPointsOfHull(inch, outch, vBezout, line, constraint);
256 cip.getAllPointsOfHull(inch2, outch2, vBezout, line2, constraint);
257
258 // std::copy( inch2.begin(), inch2.end(), std::ostream_iterator<typename CIP::Point>(std::cout, " ") );
259 // std::cout << std::endl;
260 // std::copy( outch.begin(), outch.end(), std::ostream_iterator<typename CIP::Point>(std::cout, " ") );
261 // std::cout << std::endl;
262
263 //comparisons
264 std::unique(inch2.begin(), inch2.end());
265 if (std::equal(lch.begin(), lch.end(), inch2.begin()))
266 nbok++;
267 nb++;
268
269 std::unique(outch.begin(), outch.end());
270 if (std::equal(uch.begin(), uch.end(), outch.begin()))
271 nbok++;
272 nb++;
273
274 }
275 }
276
277 trace.endBlock();
278
279 return (nbok == nb);
280}
281
286{
287 trace.beginBlock ( "Testing block (with length constraint)..." );
288
289 bool res = true;
290 res = res
295
300 ;
301
302 trace.endBlock();
303
304 return res;
305}
306
316template <typename DSL>
317DSL smartCHSubsegment(const DSL& aDSL,
318 typename DSL::Position x, typename DSL::Position y)
319{
320 ASSERT((y-x) > 0);
321
322 typedef typename DSL::Point Point;
323 typedef typename DSL::Vector Vector;
324 typedef typename DSL::Integer Integer;
325
326 Point startingPoint = aDSL.getPoint(x);
327
328 //running smartCH
329 std::vector<Point> lch, uch;
330 Vector v = smartCH( aDSL, startingPoint, (y-x),
331 std::back_inserter(uch), std::back_inserter(lch) );
332
333 //computing the slope and the remainder of the last upper leaning point
334 Point upperLeaningPoint = uch.back();
335 Integer intercept = (static_cast<Integer>(upperLeaningPoint[0])*static_cast<Integer>(v[1])
336 -static_cast<Integer>(upperLeaningPoint[1])*static_cast<Integer>(v[0]));
337
338 // trace.info() << "(" << v[1] << ", " << v[0] << ", " << intercept << ")" << std::endl;
339 return DSL( v[1],v[0],intercept );
340}
341
342
354template <typename DSL>
355DSL trivialSubsegment(const DSL& aDSL,
356 typename DSL::Position x, typename DSL::Position y)
357{
358 ASSERT((y-x) > 0);
359
360 typedef typename DSL::Point Point;
361 typedef typename DSL::Coordinate Coordinate;
362 typedef typename DSL::Integer Integer;
363 typedef typename DSL::ConstIterator ConstIterator;
365
366 Point startingPoint = aDSL.getPoint(x);
367 ASSERT( aDSL(startingPoint) );
368 Point endingPoint = aDSL.getPoint(y);
369 ASSERT( aDSL(endingPoint) );
370
371 ConstIterator it = aDSL.begin(startingPoint);
372 ConstIterator ite = aDSL.end(endingPoint);
373 ASSERT (it != ite);
374
375 DSS dss = DSS(*it);
376 for (++it; (it != ite); ++it)
377 {
378 dss.extendFront(*it);
379 }
380
381 //trace.info() << "(" << dss.a() << ", " << dss.b() << ", " << dss.mu() << ")" << std::endl;
382 return DSL(dss.a(), dss.b(), dss.mu());
383}
384
385
395template <typename DSL>
396bool comparisonSubsegment(const DSL& aDSL,
397 typename DSL::Position x, typename DSL::Position y)
398{
399 DSL dsl1 = smartCHSubsegment(aDSL, x, y);
400 DSL dsl2 = trivialSubsegment(aDSL, x, y);
401
402 return (dsl1 == dsl2);
403}
404
414template <typename DSL>
415bool comparisonSubsegment(typename DSL::Coordinate a, typename DSL::Coordinate b)
416{
417 unsigned int nbok = 0;
418 unsigned int nb = 0;
419
420 trace.beginBlock ( "Subsegment comparison ..." );
421
422 DSL aDSL(a, b, 0);
423 for (typename DSL::Integer mu = 0; ( (mu-1 >= -aDSL.omega())&&(nbok == nb) ); --mu)
424 {
425 //trace.info() << "mu=" << mu << std::endl;
426
427 typename DSL::Position f = -aDSL.patternLength();
428 for (typename DSL::Position l = 1; ( (l <= 2*aDSL.patternLength())&&(nbok == nb) ); ++l)
429 {
430 //trace.info() << "f=" << f << " l=" << l << std::endl;
431
432 if (comparisonSubsegment(DSL(a, b, mu), f, f+l))
433 nbok++;
434 nb++;
435 }
436
437 }
438
439 trace.endBlock();
440
441 return (nb == nbok);
442}
443
448{
449 trace.beginBlock ( "Testing block for the subsegment problem..." );
450
451 bool res = true;
452 res = res
464
478 ;
479
480 trace.endBlock();
481
482 return res;
483}
484
486// reversedSmartCH
494template <typename DSL>
495bool basicTest2(const DSL& aDSL)
496{
497 typedef typename DSL::Point Point;
498 typedef typename DSL::Vector Vector;
499 typedef typename DSL::Coordinate Coordinate;
500 typedef typename DSL::Integer Integer;
501 typedef typename DSL::Position Position;
502
503 unsigned int nbok = 0;
504 unsigned int nb = 0;
505
506 trace.beginBlock ( "One simple test..." );
507
508 //bounding DSS
510 Point A(0,0);
511 Position l = (2*aDSL.patternLength());
512 Point B = aDSL.getPoint( aDSL.position(A) + l + 1 );
513 DSS dss(aDSL.begin(A), aDSL.begin(B) );
514
515 trace.info() << dss << std::endl;
516
517 //computation
518 std::vector<Point> lch, uch;
520 std::back_inserter(uch), std::back_inserter(lch) );
521
522 trace.info() << v << lch.back() << uch.back() << std::endl;
523
524 if ( (uch.back() == A) && (lch.back() == A - aDSL.shift()) )
525 nbok++;
526 nb++;
527 trace.info() << "(" << nbok << "/" << nb << ") " << std::endl;
528
529 Vector LU = lch.back() - uch.back();
530 if ( (v[0]*LU[1] - v[1]*LU[0]) == NumberTraits<Coordinate>::ONE )
531 nbok++;
532 nb++;
533 trace.info() << "(" << nbok << "/" << nb << ") " << std::endl;
534
535 trace.endBlock();
536
537 return (nbok == nb);
538}
539
549template <typename DSL>
550bool basicTest2(typename DSL::Coordinate a,
551 typename DSL::Coordinate b)
552{
553 unsigned int nbok = 0;
554 unsigned int nb = 0;
555
556 DSL aDSL(a,b,0);
557 for (typename DSL::Integer mu = 0; ( (-mu < aDSL.omega())&&(nbok == nb) ); --mu)
558 {
559 if (basicTest2(DSL(a,b,mu)))
560 nbok++;
561 nb++;
562 }
563
564 return (nbok == nb);
565}
566
571{
572 trace.beginBlock ( "Testing block (without length constraint)..." );
573
574 bool res = true;
575 res = res
587
599
600
607 ;
608
609 trace.endBlock();
610
611 return res;
612}
613
621template <typename DSS>
622typename DSS::DSL reversedSmartCHSubsegment(const DSS& aDSS,
623 typename DSS::Position aBound)
624{
625 ASSERT( (aBound - aDSS.position(aDSS.back())) > 0 );
626
627 typedef typename DSS::Point Point;
628 typedef typename DSS::Vector Vector;
629 typedef typename DSS::Integer Integer;
630 typedef typename DSS::DSL DSL;
631
632 Point startingPoint = aDSS.dsl().getPoint(aBound);
633
634 //running reversedSmartCH
635 std::vector<Point> lch, uch;
636 Vector v = reversedSmartCH( aDSS, aBound,
637 std::back_inserter(uch), std::back_inserter(lch) );
638
639 //computing the slope and the remainder of the last upper leaning point
640 Point upperLeaningPoint = uch.back();
641 Integer intercept = (static_cast<Integer>(upperLeaningPoint[0])*static_cast<Integer>(v[1])
642 -static_cast<Integer>(upperLeaningPoint[1])*static_cast<Integer>(v[0]));
643
644 //trace.info() << "(" << v[1] << ", " << v[0] << ", " << intercept << ")" << std::endl;
645 return DSL( v[1],v[0],intercept );
646}
647
656template <typename DSS>
657bool comparisonSubsegment2(const DSS& aDSS, typename DSS::Position aBound)
658{
659 typedef typename DSS::DSL DSL;
660 DSL dsl1 = reversedSmartCHSubsegment(aDSS, aBound);
661 DSL dsl2 = trivialSubsegment(aDSS.dsl(), aDSS.position(aDSS.back()), aBound);
662
663 return (dsl1 == dsl2);
664}
665
666
667
677template <typename DSL>
678bool comparisonSubsegment2(typename DSL::Coordinate a, typename DSL::Coordinate b)
679{
680 unsigned int nbok = 0;
681 unsigned int nb = 0;
682
683 trace.beginBlock ( "Subsegment comparison ..." );
684
685 DSL aDSL(a, b, 0);
686 for (typename DSL::Integer mu = 0; ( (mu-1 >= -aDSL.omega())&&(nbok == nb) ); --mu)
687 {
688 trace.info() << "mu=" << mu << std::endl;
689
690 //computation of a bounding DSS
691 typedef typename DSL::Point Point;
692 typedef typename DSL::Coordinate Coordinate;
693 typedef typename DSL::Integer Integer;
695
696 Point startingPoint = aDSL.getPoint(0);
697 ASSERT( aDSL(startingPoint) );
698 Point endingPoint = aDSL.getPoint(2*aDSL.patternLength()+1);
699 ASSERT( aDSL(endingPoint) );
700
701 DSS dss = DSS(aDSL.begin(startingPoint), aDSL.begin(endingPoint));
702
703 //test for a left subsegment
704 for (typename DSL::Position l = 1; ( (l <= 2*aDSL.patternLength())&&(nbok == nb) ); ++l)
705 {
706 trace.info() << "l=" << l << std::endl;
707
708 if (comparisonSubsegment2(dss, l))
709 nbok++;
710 nb++;
711 }
712
713 }
714
715 trace.endBlock();
716
717 return (nb == nbok);
718}
719
724{
725 trace.beginBlock ( "Testing block for the subsegment problem..." );
726
727 bool res = true;
728 res = res
740
754 ;
755
756 trace.endBlock();
757
758 return res;
759}
760
762// Standard services - public :
763int main( int argc, char** argv )
764{
765 trace.beginBlock ( "Testing class ArithmeticalDSSConvexHull" );
766 trace.info() << "Args:";
767 for ( int i = 0; i < argc; ++i )
768 trace.info() << " " << argv[ i ];
769 trace.info() << endl;
770
771 if (argc >= 4)
772 { //basic test with a specified DSL
773 std::istringstream issb(argv[1]);
774 int b;
775 issb >> b;
776 std::istringstream issa(argv[2]);
777 int a;
778 issa >> a;
779 if ( (a <= 0)||(b <= 0) )
780 {
781 std::cerr << " a and b should be strictly positive " << std::endl;
782 return 1;
783 }
784 std::istringstream issmu(argv[3]);
785 int mu;
786 issmu >> mu;
787 if ( (mu > 0)||(mu <= -(std::max(std::abs(a),std::abs(b)))) )
788 {
789 std::cerr << " mu should be within the range ]-max(|a|,|b|); 0] " << std::endl;
790 return 1;
791 }
792
793 bool res = basicTest( NaiveDSL<DGtal::int32_t>(a,b,mu) )
795 trace.emphase() << ( res ? "Passed." : "Error." ) << endl;
796 return res ? 0 : 1;
797 }
798
799 //all automatic tests
800 bool res = true;
801 res = res
804 && testSubsegment()
806 && testSubsegment2()
807 ;
808
809 trace.emphase() << ( res ? "Passed." : "Error." ) << endl;
810 trace.endBlock();
811 return res ? 0 : 1;
812}
813// //
Aim: This class represents a naive (resp. standard) digital straight segment (DSS),...
Aim: Represents a 2D polytope, i.e. a convex polygon, in the two-dimensional digital plane....
Aim: This class is an alias of ArithmeticalDSS for naive DSL. It represents a naive digital straight ...
Aim: This class is an alias of ArithmeticalDSS for standard DSL. It represents a standard digital str...
DigitalPlane::Point Vector
functions namespace gathers all DGtal functionsxs.
PointVector smartCH(const PointVector &aFirstPoint, const Coordinate &aRemainderBound, const Position &aPositionBound, const PointVector &aStep, const Coordinate &aRStep, const PointVector &aShift, const Coordinate &aRShift, const PositionFunctor &aPositionFunctor, OutputIterator uIto, OutputIterator lIto)
Procedure that computes the lower and upper left hull of a DSS of first point aFirstPoint,...
PointVector reversedSmartCH(PointVector U, PointVector L, PointVector V, const Position &aFirstPosition, const Position &aLastPosition, const PositionFunctor &aPositionFunctor, OutputIterator uIto, OutputIterator lIto)
Procedure that computes the lower and upper left hull of the left subsegment of a greater DSS charact...
DGtal is the top-level namespace which contains all DGtal functions and types.
Trace trace
STL namespace.
Aim: The traits class for all models of Cinteger.
DSL smartCHSubsegment(const DSL &aDSL, typename DSL::Position x, typename DSL::Position y)
bool comparisonLeftHull(typename DSL::Coordinate a, typename DSL::Coordinate b)
bool testWithoutLengthConstraint2()
bool basicTest2(const DSL &aDSL)
bool testWithoutLengthConstraint()
bool comparisonSubsegment2(const DSS &aDSS, typename DSS::Position aBound)
bool basicTest(const DSL &aDSL)
DSL trivialSubsegment(const DSL &aDSL, typename DSL::Position x, typename DSL::Position y)
bool testWithLengthConstraint()
DSS::DSL reversedSmartCHSubsegment(const DSS &aDSS, typename DSS::Position aBound)
bool comparisonSubsegment(const DSL &aDSL, typename DSL::Position x, typename DSL::Position y)
bool comparisonSubsegment(const DSS &aDSS, typename DSS::Position x, typename DSS::Position y)
int main()
Definition testBits.cpp:56
MyPointD Point