DGtal  1.2.0
testArithmeticalDSSConvexHull.cpp
Go to the documentation of this file.
1
31 #include <iostream>
32 #include "DGtal/base/Common.h"
33 #include "DGtal/geometry/curves/ArithmeticalDSL.h"
34 #include "DGtal/geometry/curves/ArithmeticalDSSConvexHull.h"
35 #include "DGtal/geometry/curves/ArithmeticalDSS.h"
36
37 #include "DGtal/kernel/SpaceND.h"
38 #include "DGtal/kernel/NumberTraits.h"
39 #include "DGtal/arithmetic/LatticePolytope2D.h"
40
42
43 using namespace std;
44 using namespace DGtal;
45 using namespace functions;
46
48 // Functions for testing functions smartCH and reversedSmartCH.
50
52 // smartCH
60 template <typename DSL>
62 {
63  typedef typename DSL::Point Point;
64  typedef typename DSL::Vector Vector;
65  typedef typename DSL::Integer Integer;
66  typedef typename DSL::Position Position;
67
68  unsigned int nbok = 0;
69  unsigned int nb = 0;
70
71  trace.beginBlock ( "One simple test..." );
72
73  trace.info() << aDSL << std::endl;
74
75  //forward test
77  std::vector<Point> lch, uch;
78  Vector v = smartCH( aDSL, Point(0,0), l,
79  std::back_inserter(uch), std::back_inserter(lch) );
80
82  nbok++;
83  nb++;
84  trace.info() << "(" << nbok << "/" << nb << ") " << std::endl;
86  nbok++;
87  nb++;
88  trace.info() << "(" << nbok << "/" << nb << ") " << std::endl; //eq. 8
90  nbok++;
91  nb++;
92  trace.info() << "(" << nbok << "/" << nb << ") " << std::endl;
93
94  //test on the size of the convex hulls
96  unsigned int threshold = (unsigned int) std::ceil( std::log(bound) / std::log(1.618) );
97  if ( (lch.size()+uch.size()-1) <= threshold )
98  nbok++;
99  nb++;
100  trace.info() << "(" << nbok << "/" << nb << ") " << std::endl;
101
102  trace.endBlock();
103
104  return (nbok == nb);
105 }
106
116 template <typename DSL>
117 bool basicTest(typename DSL::Coordinate a,
118  typename DSL::Coordinate b)
119 {
120  unsigned int nbok = 0;
121  unsigned int nb = 0;
122
124  for (typename DSL::Integer mu = 0; ( (-mu < aDSL.omega())&&(nbok == nb) ); --mu)
125  {
126  if (basicTest(DSL(a,b,mu)))
127  nbok++;
128  nb++;
129  }
130
131  return (nbok == nb);
132 }
133
138 {
139  trace.beginBlock ( "Testing block (without length constraint)..." );
140
141  bool res = true;
142  res = res
143  && basicTest<NaiveDSL<DGtal::int32_t> >(5,8)
145  && basicTest<NaiveDSL<DGtal::int32_t> >(12,29)
147  && basicTest<NaiveDSL<DGtal::int32_t> >(70,29)
149  && basicTest<NaiveDSL<DGtal::int32_t> >(-70,29)
151  && basicTest<NaiveDSL<DGtal::int32_t> >(70,-29)
152  && basicTest<NaiveDSL<DGtal::int32_t> >(-29,-70)
153  && basicTest<NaiveDSL<DGtal::int32_t> >(-70,-29)
154
156  && basicTest<StandardDSL<DGtal::int32_t> >(8,13)
158  && basicTest<StandardDSL<DGtal::int32_t> >(29,70)
160  && basicTest<StandardDSL<DGtal::int32_t> >(-29,70)
162  && basicTest<StandardDSL<DGtal::int32_t> >(29,-70)
164  && basicTest<StandardDSL<DGtal::int32_t> >(-29,-70)
166
167
168  && basicTest<NaiveDSL<DGtal::int32_t> >(33, 109)
169  && basicTest<NaiveDSL<DGtal::int32_t> >(109, 360)
170  && basicTest<NaiveDSL<DGtal::int32_t> >(8, 73)
172 #ifdef WITH_BIGINTEGER
173  && basicTest<NaiveDSL<DGtal::int32_t,DGtal::BigInteger> >(57, 520)
175 #endif
176
177  ;
178
179  trace.endBlock();
180
181  return res;
182 }
183
184
195 template <typename DSL>
196 bool comparisonLeftHull(typename DSL::Coordinate a, typename DSL::Coordinate b)
197 {
198  ASSERT(a >= 0);
199  ASSERT(b >= 0);
200  ASSERT(a < b);
201
202  typedef typename DSL::Point Point;
203  typedef typename DSL::Vector Vector;
204  typedef typename DSL::Coordinate Coordinate;
205
206  unsigned int nbok = 0;
207  unsigned int nb = 0;
208
209  trace.beginBlock ( "Comparison ..." );
210  DSL inputDSL(a,b,0);
211
212  trace.info() << a << " " << b << std::endl;
213  trace.info() << "testing every mu between 0 and -" << inputDSL.omega() << std::endl;
214
215  for (typename DSL::Integer mu = 0; ( (mu-1 >= -inputDSL.omega())&&(nbok == nb) ); --mu)
216  {
217  trace.info() << "mu=" << mu << ", testing every length between 1 and 2*" << inputDSL.omega() << std::endl;
218  inputDSL = DSL(a,b,mu);
219
220  for (typename DSL::Position l = 1; ( (l <= 2*inputDSL.patternLength())&&(nbok == nb) ); ++l)
221  {
222  //trace.info() << "l=" << l << std::endl;
223
224  //smartCH
225  std::vector<Point> lch, uch;
226  Vector v = smartCH( inputDSL, Point(0,0), l,
227  std::back_inserter(uch), std::back_inserter(lch) );
228
229  // std::copy( lch.begin(), lch.end(), std::ostream_iterator<Point>(std::cout, " ") );
230  // std::cout << std::endl;
231  // std::copy( uch.begin(), uch.end(), std::ostream_iterator<Point>(std::cout, " ") );
232  // std::cout << std::endl;
233
234  Vector shift = -inputDSL.shift();
235  //algorithm of Charrier and Buzer
236  typedef SpaceND<2, Coordinate> Space2;
237  typedef LatticePolytope2D<Space2> CIP;
238  CIP cip;
239  typename CIP::HalfSpace line( typename CIP::Vector(a,-b), mu );
240  typename CIP::HalfSpace line2( typename CIP::Vector(a,-b), mu-1 );
241  //NB: since only closed half-space are used,
242  // we must run the procedure twice, for mu and mu-1
243  // in order to get the lower left hull (not included)
244  // and the upper left hull (included)
245  typename CIP::HalfSpace constraint( typename CIP::Vector(shift[1],-shift[0]), l );
246  std::vector<typename CIP::Point> inch, outch, inch2, outch2;
247  inch.push_back( typename CIP::Point(shift[0],shift[1]) );
248  inch2.push_back( typename CIP::Point(shift[0],shift[1]) );
249  ASSERT( line(inch[0]) );
250  ASSERT( constraint(inch[0]) );
251  outch.push_back( typename CIP::Point(0,0) );
252  outch2.push_back( typename CIP::Point(0,0) );
253  ASSERT( (!line2(outch[0])) );
254  ASSERT( constraint(outch[0]) );
255  typename CIP::Vector vBezout(1,0);
256  cip.getAllPointsOfHull(inch, outch, vBezout, line, constraint);
257  cip.getAllPointsOfHull(inch2, outch2, vBezout, line2, constraint);
258
259  // std::copy( inch2.begin(), inch2.end(), std::ostream_iterator<typename CIP::Point>(std::cout, " ") );
260  // std::cout << std::endl;
261  // std::copy( outch.begin(), outch.end(), std::ostream_iterator<typename CIP::Point>(std::cout, " ") );
262  // std::cout << std::endl;
263
264  //comparisons
265  std::unique(inch2.begin(), inch2.end());
266  if (std::equal(lch.begin(), lch.end(), inch2.begin()))
267  nbok++;
268  nb++;
269
270  std::unique(outch.begin(), outch.end());
271  if (std::equal(uch.begin(), uch.end(), outch.begin()))
272  nbok++;
273  nb++;
274
275  }
276  }
277
278  trace.endBlock();
279
280  return (nbok == nb);
281 }
282
287 {
288  trace.beginBlock ( "Testing block (with length constraint)..." );
289
290  bool res = true;
291  res = res
292  && comparisonLeftHull<NaiveDSL<DGtal::int32_t> >(5,8)
294  && comparisonLeftHull<NaiveDSL<DGtal::int32_t> >(12,29)
296
297  && comparisonLeftHull<StandardDSL<DGtal::int32_t> >(5,8)
299  && comparisonLeftHull<StandardDSL<DGtal::int32_t> >(12,29)
301  ;
302
303  trace.endBlock();
304
305  return res;
306 }
307
317 template <typename DSL>
319  typename DSL::Position x, typename DSL::Position y)
320 {
321  ASSERT((y-x) > 0);
322
323  typedef typename DSL::Point Point;
324  typedef typename DSL::Vector Vector;
325  typedef typename DSL::Integer Integer;
326
328
329  //running smartCH
330  std::vector<Point> lch, uch;
331  Vector v = smartCH( aDSL, startingPoint, (y-x),
332  std::back_inserter(uch), std::back_inserter(lch) );
333
334  //computing the slope and the remainder of the last upper leaning point
335  Point upperLeaningPoint = uch.back();
336  Integer intercept = (static_cast<Integer>(upperLeaningPoint[0])*static_cast<Integer>(v[1])
337  -static_cast<Integer>(upperLeaningPoint[1])*static_cast<Integer>(v[0]));
338
339  // trace.info() << "(" << v[1] << ", " << v[0] << ", " << intercept << ")" << std::endl;
340  return DSL( v[1],v[0],intercept );
341 }
342
343
355 template <typename DSL>
357  typename DSL::Position x, typename DSL::Position y)
358 {
359  ASSERT((y-x) > 0);
360
361  typedef typename DSL::Point Point;
362  typedef typename DSL::Coordinate Coordinate;
363  typedef typename DSL::Integer Integer;
364  typedef typename DSL::ConstIterator ConstIterator;
366
371
374  ASSERT (it != ite);
375
376  DSS dss = DSS(*it);
377  for (++it; (it != ite); ++it)
378  {
379  dss.extendFront(*it);
380  }
381
382  //trace.info() << "(" << dss.a() << ", " << dss.b() << ", " << dss.mu() << ")" << std::endl;
383  return DSL(dss.a(), dss.b(), dss.mu());
384 }
385
386
396 template <typename DSL>
398  typename DSL::Position x, typename DSL::Position y)
399 {
400  DSL dsl1 = smartCHSubsegment(aDSL, x, y);
401  DSL dsl2 = trivialSubsegment(aDSL, x, y);
402
403  return (dsl1 == dsl2);
404 }
405
415 template <typename DSL>
416 bool comparisonSubsegment(typename DSL::Coordinate a, typename DSL::Coordinate b)
417 {
418  unsigned int nbok = 0;
419  unsigned int nb = 0;
420
421  trace.beginBlock ( "Subsegment comparison ..." );
422
424  for (typename DSL::Integer mu = 0; ( (mu-1 >= -aDSL.omega())&&(nbok == nb) ); --mu)
425  {
426  //trace.info() << "mu=" << mu << std::endl;
427
428  typename DSL::Position f = -aDSL.patternLength();
429  for (typename DSL::Position l = 1; ( (l <= 2*aDSL.patternLength())&&(nbok == nb) ); ++l)
430  {
431  //trace.info() << "f=" << f << " l=" << l << std::endl;
432
433  if (comparisonSubsegment(DSL(a, b, mu), f, f+l))
434  nbok++;
435  nb++;
436  }
437
438  }
439
440  trace.endBlock();
441
442  return (nb == nbok);
443 }
444
449 {
450  trace.beginBlock ( "Testing block for the subsegment problem..." );
451
452  bool res = true;
453  res = res
454  && comparisonSubsegment<NaiveDSL<DGtal::int32_t> >(5,8)
456  && comparisonSubsegment<NaiveDSL<DGtal::int32_t> >(12,29)
458  && comparisonSubsegment<NaiveDSL<DGtal::int32_t> >(8,5)
460  && comparisonSubsegment<NaiveDSL<DGtal::int32_t> >(-8,5)
462  && comparisonSubsegment<NaiveDSL<DGtal::int32_t> >(8,-5)
464  && comparisonSubsegment<NaiveDSL<DGtal::int32_t> >(-8,-5)
465
467  && comparisonSubsegment<StandardDSL<DGtal::int32_t> >(8,13)
469  && comparisonSubsegment<StandardDSL<DGtal::int32_t> >(29,70)
471  && comparisonSubsegment<StandardDSL<DGtal::int32_t> >(-5,8)
473  && comparisonSubsegment<StandardDSL<DGtal::int32_t> >(5,-8)
475  && comparisonSubsegment<StandardDSL<DGtal::int32_t> >(-5,-8)
477
478 #ifdef WITH_BIGINTEGER
479  && comparisonSubsegment<NaiveDSL<DGtal::int32_t,DGtal::BigInteger> >(5,8)
481 #endif
482  ;
483
484  trace.endBlock();
485
486  return res;
487 }
488
490 // reversedSmartCH
498 template <typename DSL>
500 {
501  typedef typename DSL::Point Point;
502  typedef typename DSL::Vector Vector;
503  typedef typename DSL::Coordinate Coordinate;
504  typedef typename DSL::Integer Integer;
505  typedef typename DSL::Position Position;
506
507  unsigned int nbok = 0;
508  unsigned int nb = 0;
509
510  trace.beginBlock ( "One simple test..." );
511
512  //bounding DSS
514  Point A(0,0);
518
519  trace.info() << dss << std::endl;
520
521  //computation
522  std::vector<Point> lch, uch;
524  std::back_inserter(uch), std::back_inserter(lch) );
525
526  trace.info() << v << lch.back() << uch.back() << std::endl;
527
528  if ( (uch.back() == A) && (lch.back() == A - aDSL.shift()) )
529  nbok++;
530  nb++;
531  trace.info() << "(" << nbok << "/" << nb << ") " << std::endl;
532
533  Vector LU = lch.back() - uch.back();
534  if ( (v[0]*LU[1] - v[1]*LU[0]) == NumberTraits<Coordinate>::ONE )
535  nbok++;
536  nb++;
537  trace.info() << "(" << nbok << "/" << nb << ") " << std::endl;
538
539  trace.endBlock();
540
541  return (nbok == nb);
542 }
543
553 template <typename DSL>
554 bool basicTest2(typename DSL::Coordinate a,
555  typename DSL::Coordinate b)
556 {
557  unsigned int nbok = 0;
558  unsigned int nb = 0;
559
561  for (typename DSL::Integer mu = 0; ( (-mu < aDSL.omega())&&(nbok == nb) ); --mu)
562  {
563  if (basicTest2(DSL(a,b,mu)))
564  nbok++;
565  nb++;
566  }
567
568  return (nbok == nb);
569 }
570
575 {
576  trace.beginBlock ( "Testing block (without length constraint)..." );
577
578  bool res = true;
579  res = res
580  && basicTest2<NaiveDSL<DGtal::int32_t> >(5,8)
582  && basicTest2<NaiveDSL<DGtal::int32_t> >(12,29)
584  && basicTest2<NaiveDSL<DGtal::int32_t> >(70,29)
586  && basicTest2<NaiveDSL<DGtal::int32_t> >(-70,29)
588  && basicTest2<NaiveDSL<DGtal::int32_t> >(70,-29)
590  && basicTest2<NaiveDSL<DGtal::int32_t> >(-70,-29)
591
593  && basicTest2<StandardDSL<DGtal::int32_t> >(8,13)
595  && basicTest2<StandardDSL<DGtal::int32_t> >(29,70)
597  && basicTest2<StandardDSL<DGtal::int32_t> >(-29,70)
599  && basicTest2<StandardDSL<DGtal::int32_t> >(29,-70)
601  && basicTest2<StandardDSL<DGtal::int32_t> >(-29,-70)
603
604
605  && basicTest2<NaiveDSL<DGtal::int32_t> >(33, 109)
606  && basicTest2<NaiveDSL<DGtal::int32_t> >(109, 360)
607  && basicTest2<NaiveDSL<DGtal::int32_t> >(8, 73)
609 #ifdef WITH_BIGINTEGER
610  && basicTest2<NaiveDSL<DGtal::int32_t,DGtal::BigInteger> >(57, 520)
612 #endif
613
614  ;
615
616  trace.endBlock();
617
618  return res;
619 }
620
628 template <typename DSS>
629 typename DSS::DSL reversedSmartCHSubsegment(const DSS& aDSS,
630  typename DSS::Position aBound)
631 {
633
634  typedef typename DSS::Point Point;
635  typedef typename DSS::Vector Vector;
636  typedef typename DSS::Integer Integer;
637  typedef typename DSS::DSL DSL;
638
640
641  //running reversedSmartCH
642  std::vector<Point> lch, uch;
643  Vector v = reversedSmartCH( aDSS, aBound,
644  std::back_inserter(uch), std::back_inserter(lch) );
645
646  //computing the slope and the remainder of the last upper leaning point
647  Point upperLeaningPoint = uch.back();
648  Integer intercept = (static_cast<Integer>(upperLeaningPoint[0])*static_cast<Integer>(v[1])
649  -static_cast<Integer>(upperLeaningPoint[1])*static_cast<Integer>(v[0]));
650
651  //trace.info() << "(" << v[1] << ", " << v[0] << ", " << intercept << ")" << std::endl;
652  return DSL( v[1],v[0],intercept );
653 }
654
663 template <typename DSS>
664 bool comparisonSubsegment2(const DSS& aDSS, typename DSS::Position aBound)
665 {
666  typedef typename DSS::DSL DSL;
667  DSL dsl1 = reversedSmartCHSubsegment(aDSS, aBound);
669
670  return (dsl1 == dsl2);
671 }
672
673
674
684 template <typename DSL>
685 bool comparisonSubsegment2(typename DSL::Coordinate a, typename DSL::Coordinate b)
686 {
687  unsigned int nbok = 0;
688  unsigned int nb = 0;
689
690  trace.beginBlock ( "Subsegment comparison ..." );
691
693  for (typename DSL::Integer mu = 0; ( (mu-1 >= -aDSL.omega())&&(nbok == nb) ); --mu)
694  {
695  trace.info() << "mu=" << mu << std::endl;
696
697  //computation of a bounding DSS
698  typedef typename DSL::Point Point;
699  typedef typename DSL::Coordinate Coordinate;
700  typedef typename DSL::Integer Integer;
702
707
709
710  //test for a left subsegment
711  for (typename DSL::Position l = 1; ( (l <= 2*aDSL.patternLength())&&(nbok == nb) ); ++l)
712  {
713  trace.info() << "l=" << l << std::endl;
714
715  if (comparisonSubsegment2(dss, l))
716  nbok++;
717  nb++;
718  }
719
720  }
721
722  trace.endBlock();
723
724  return (nb == nbok);
725 }
726
731 {
732  trace.beginBlock ( "Testing block for the subsegment problem..." );
733
734  bool res = true;
735  res = res
736  && comparisonSubsegment2<NaiveDSL<DGtal::int32_t> >(5,8)
738  && comparisonSubsegment2<NaiveDSL<DGtal::int32_t> >(12,29)
740  && comparisonSubsegment2<NaiveDSL<DGtal::int32_t> >(8,5)
742  && comparisonSubsegment2<NaiveDSL<DGtal::int32_t> >(-8,5)
744  && comparisonSubsegment2<NaiveDSL<DGtal::int32_t> >(8,-5)
746  && comparisonSubsegment2<NaiveDSL<DGtal::int32_t> >(-8,-5)
747
749  && comparisonSubsegment2<StandardDSL<DGtal::int32_t> >(8,13)
751  && comparisonSubsegment2<StandardDSL<DGtal::int32_t> >(29,70)
753  && comparisonSubsegment2<StandardDSL<DGtal::int32_t> >(-5,8)
755  && comparisonSubsegment2<StandardDSL<DGtal::int32_t> >(5,-8)
757  && comparisonSubsegment2<StandardDSL<DGtal::int32_t> >(-5,-8)
759
760 #ifdef WITH_BIGINTEGER
761  && comparisonSubsegment2<NaiveDSL<DGtal::int32_t,DGtal::BigInteger> >(5,8)
763 #endif
764  ;
765
766  trace.endBlock();
767
768  return res;
769 }
770
772 // Standard services - public :
773 int main( int argc, char** argv )
774 {
775  trace.beginBlock ( "Testing class ArithmeticalDSSConvexHull" );
776  trace.info() << "Args:";
777  for ( int i = 0; i < argc; ++i )
778  trace.info() << " " << argv[ i ];
779  trace.info() << endl;
780
781  if (argc >= 4)
782  { //basic test with a specified DSL
783  std::istringstream issb(argv[1]);
784  int b;
785  issb >> b;
786  std::istringstream issa(argv[2]);
787  int a;
788  issa >> a;
789  if ( (a <= 0)||(b <= 0) )
790  {
791  std::cerr << " a and b should be strictly positive " << std::endl;
792  return 1;
793  }
794  std::istringstream issmu(argv[3]);
795  int mu;
796  issmu >> mu;
797  if ( (mu > 0)||(mu <= -(std::max(std::abs(a),std::abs(b)))) )
798  {
799  std::cerr << " mu should be within the range ]-max(|a|,|b|); 0] " << std::endl;
800  return 1;
801  }
802
803  bool res = basicTest( NaiveDSL<DGtal::int32_t>(a,b,mu) )
805  trace.emphase() << ( res ? "Passed." : "Error." ) << endl;
806  return res ? 0 : 1;
807  }
808
809  //all automatic tests
810  bool res = true;
811  res = res
814  && testSubsegment()
816  && testSubsegment2()
817  ;
818
819  trace.emphase() << ( res ? "Passed." : "Error." ) << endl;
820  trace.endBlock();
821  return res ? 0 : 1;
822 }
823 // //
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...
void beginBlock(const std::string &keyword="")
std::ostream & emphase()
std::ostream & info()
double endBlock()
Point::Coordinate Integer
MyDigitalSurface::ConstIterator ConstIterator
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
Definition: Common.h:154
Aim: The traits class for all models of Cinteger.
Definition: NumberTraits.h:533
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()
int main(int argc, char **argv)