DGtal  0.9.4.1
FrechetShortcut.h
1 
18 #pragma once
19 
32 #if defined(FrechetShortcut_RECURSES)
33 #error Recursive header files inclusion detected in FrechetShortcut.h
34 #else // defined(FrechetShortcut_RECURSES)
35 
36 #define FrechetShortcut_RECURSES
37 
38 #if !defined FrechetShortcut_h
39 
40 #define FrechetShortcut_h
41 
43 // Inclusions
44 #include "DGtal/base/Common.h"
45 #include "DGtal/base/ConceptUtils.h"
46 #include "DGtal/kernel/PointVector.h"
47 #include "DGtal/arithmetic/IntegerComputer.h"
48 #include <boost/icl/interval_set.hpp>
49 #include <map>
50 
52 
53 #include "DGtal/geometry/curves/SegmentComputerUtils.h"
54 
55 namespace DGtal
56 {
57 
59  // class FrechetShortcut
106  template <typename TIterator,typename TInteger = typename IteratorCirculatorTraits<TIterator>::Value::Coordinate>
108  {
109  // ----------------------- Standard services ------------------------------
110  public:
111 
112  //entier
113 
115  typedef TInteger Integer;
116 
117 
118 
119  //required types
120  typedef TIterator ConstIterator;
123 
124  //2D point and 2D vector
127  typedef typename Vector::Coordinate Coordinate;
128 
133  class Backpath
134  {
135  private:
140 
141  protected:
142 
147  typedef struct occulter_attributes{
148  double angle_min; //
149  double angle_max; //
151 
156  typedef std::map <ConstIterator,occulter_attributes > occulter_list;
157 
158  public:
160 
161 
162  public:
163 
164  typedef boost::icl::interval_set<double> IntervalSet;
165 
169  int myQuad;
170 
175  bool myFlag;
176 
178 
184 
189 
193  Backpath();
194 
201 
206  Backpath(const Backpath & other);
207 
208 
214  Backpath& operator=(const Backpath & other);
215 
216 
220  ~Backpath();
221 
225  void reset();
226 
230  void addPositivePoint();
231 
235  void addNegativePoint();
236 
243  void updateBackPathFirstQuad(int d, const ConstIterator&);
244 
248  void updateOcculters();
249 
253  void updateIntervals();
254 
255 
256  }; // End of class Backpath
257 
258 
263  class Cone{
264 
265  public:
269  double myMin;
270 
274  double myMax;
275 
279  bool myInf; // true if the cone is infinite (the whole plane)
280 
281  Cone();
282 
288  Cone(double a0, double a1);
289 
300  Cone(double x, double y, double x0, double y0, double x1, double y1);
301 
306  bool isEmpty();
307 
313  Cone& operator=(const Cone& c);
314 
319  void intersectCones(Cone c);
320 
327 
333 
338  void selfDisplay ( std::ostream & out) const;
339 
340  };
341 
342 
343 
344 
345  struct Tools
346  {
356  static bool isBetween(double i, double a, double b, double n)
357  {
358  if(a<=b)
359  return (i>=a && i<=b);
360  else
361  return ((i>=a && i<=n) || (i>=0 && i<=b));
362  }
363 
364 
365  // Code by Tim Voght
366  // http://local.wasp.uwa.edu.au/~pbourke/geometry/2circle/tvoght.c
367 
368  /* circle_circle_intersection() *
369  * Determine the points where 2 circles in a common plane intersect.
370  *
371  * int circle_circle_intersection(
372  * // center and radius of 1st circle
373  * double x0, double y0, double r0,
374  * // center and radius of 2nd circle
375  * double x1, double y1, double r1,
376  * // 1st intersection point
377  * double *xi, double *yi,
378  * // 2nd intersection point
379  * double *xi_prime, double *yi_prime)
380  *
381  * This is a public domain work. 3/26/2005 Tim Voght
382  *
383  */
391  static int circle_circle_intersection(double x0, double y0, double r0,
392  double x1, double y1, double r1,
393  double *xi, double *yi,
394  double *xi_prime, double *yi_prime)
395  {
396  double a, dx, dy, d, h, rx, ry;
397  double x2, y2;
398 
399  /* dx and dy are the vertical and horizontal distances between
400  * the circle centers.
401  */
402  dx = x1 - x0;
403  dy = y1 - y0;
404 
405  /* Determine the straight-line distance between the centers. */
406  //d = sqrt((dy*dy) + (dx*dx));
407  d = hypot(dx,dy); // Suggested by Keith Briggs
408 
409  /* Check for solvability. */
410  if (d > (r0 + r1))
411  {
412  /* no solution. circles do not intersect. */
413  std::cerr << "Warning : the two circles do not intersect -> should never happen" << std::endl;
414  return 0; //I. Sivignon 05/2010 should never happen for our specific use.
415  }
416  if (d < fabs(r0 - r1))
417  {
418  /* no solution. one circle is contained in the other */
419  return 0;
420  }
421 
422  /* 'point 2' is the point where the line through the circle
423  * intersection points crosses the line between the circle
424  * centers.
425  */
426 
427  /* Determine the distance from point 0 to point 2. */
428  a = ((r0*r0) - (r1*r1) + (d*d)) / (2.0 * d) ;
429 
430  /* Determine the coordinates of point 2. */
431  x2 = x0 + (dx * a/d);
432  y2 = y0 + (dy * a/d);
433 
434  /* Determine the distance from point 2 to either of the
435  * intersection points.
436  */
437  h = sqrt((r0*r0) - (a*a));
438 
439  /* Now determine the offsets of the intersection points from
440  * point 2.
441  */
442  rx = -dy * (h/d);
443  ry = dx * (h/d);
444 
445  /* Determine the absolute intersection points. */
446  *xi = x2 + rx;
447  *xi_prime = x2 - rx;
448  *yi = y2 + ry;
449  *yi_prime = y2 - ry;
450 
451  return 1;
452  }
453 
454 
455 
475  static int circleTangentPoints(double x, double y, double x1, double y1, double r1, double *xi, double *yi,
476  double *xi_prime, double *yi_prime)
477  {
478  double x0 = (x+x1)/2;
479  double y0 = (y+y1)/2;
480  double r0 = sqrt((x-x1)*(x-x1) + (y-y1)*(y-y1))/2;
481 
482  int res =
483  circle_circle_intersection(x0,y0,r0,x1,y1,r1,xi,yi,xi_prime,yi_prime);
484 
485  return res;
486 
487  }
488 
489 
499  static double computeAngle(double x0, double y0, double x1, double y1)
500  {
501  double x = x1-x0;
502  double y = y1-y0;
503 
504  if(x!=0)
505  {
506  double alpha = y/x;
507 
508  if(x>0 && y>=0)
509  return atan(alpha);
510  else
511  if(x>0 && y<0)
512  return atan(alpha)+2*M_PI;
513  else
514  if(x<0)
515  return atan(alpha)+M_PI;
516  }
517  else
518  {
519  if(y>0)
520  return M_PI_2;
521  else
522  return 3*M_PI_2;
523  }
524  return -1;
525  }
526 
527 
528 
534  static double angleVectVect(Vector u, Vector v)
535  {
537  return acos((double)ic.dotProduct(u,v)/(u.norm()*v.norm()));
538  }
539 
545  static int computeChainCode(Point p, Point q)
546  {
547  int d;
548  Coordinate x2 = q[0];
549  Coordinate y2 = q[1];
550  Coordinate x1 = p[0];
551  Coordinate y1 = p[1];
552 
553  if(x2-x1==0)
554  if(y2-y1==1)
555  d=2;
556  else
557  d=6;
558  else
559  if(x2-x1==1)
560  if(y2-y1==0)
561  d=0;
562  else
563  if(y2-y1==1)
564  d=1;
565  else
566  d=7;
567  else
568  if(y2-y1==0)
569  d=4;
570  else
571  if(y2-y1==1)
572  d=3;
573  else
574  d=5;
575  return d;
576  }
577 
583  static int computeOctant(Point p, Point q)
584  {
585  int d = 0;
586  Coordinate x = q[0]-p[0];
587  Coordinate y = q[1]-p[1];
588 
589  if(x>=0)
590  {
591  if(y>=0)
592  {
593  if(x>y)
594  d=0; // 0 <= y < x
595  else
596  if(x!=0)
597  d=1; // 0 <= x <= y
598  }
599  else
600  {
601  if(x>=abs(y))
602  d=7; // 0 < abs(y) <= x
603  else
604  d=6; // 0 <= x < abs(y)
605  }
606  }
607  if(x<=0)
608  {
609  if(y>0)
610  if(abs(x)>=y)
611  d=3; //
612  else
613  d=2;
614  else
615  {
616  if(abs(x)>abs(y))
617  d=4;
618  else
619  if(x!=0)
620  d=5;
621  }
622  }
623  return d;
624 
625  }
626 
635  static int rot(int d, int quad)
636  {
637  return (d-quad+8)%8;
638  }
639 
645  static Vector chainCode2Vect(int d)
646  {
647  Vector v;
648  switch(d){
649 
650  case 0:{
651  v[0] = 1;
652  v[1] = 0;
653  break;
654  }
655  case 1:{
656  v[0] = 1;
657  v[1] = 1;
658  break;
659  }
660  case 2:{
661  v[0] = 0;
662  v[1] = 1;
663  break;
664  }
665  case 3:{
666  v[0] = -1;
667  v[1] = 1;
668  break;
669  }
670  case 4:{
671  v[0] = -1;
672  v[1] = 0;
673  break;
674  }
675  case 5:{
676  v[0] = -1;
677  v[1] = -1;
678  break;
679  }
680  case 6:{
681  v[0] = 0;
682  v[1] = -1;
683  break;
684  }
685  case 7:{
686  v[0] = 1;
687  v[1] = -1;
688  break;
689  }
690 
691  }
692 
693  return v;
694  }
695 
696 
697  };
698 
699 
704  FrechetShortcut();
705 
710  FrechetShortcut(double error);
711 
717  void init(const ConstIterator& it);
718 
720 
725  FrechetShortcut ( const Self & other );
726 
732  FrechetShortcut & operator= ( const Self & other );
733 
737  Reverse getReverse() const;
738 
745  bool operator==( const Self & other ) const;
746 
753  bool operator!=( const Self & other ) const;
754 
755 
756 
761 
762  // ----------------------- Interface --------------------------------------
763 public:
764 
770  bool isExtendableFront();
771 
777  bool extendFront();
778 
779 
780  // ---------------------------- Accessors ----------------------------------
781 
782 
787  ConstIterator begin() const;
791  ConstIterator end() const;
792 
793 
794 
795 
796  public:
797 
801  std::string className() const;
802 
803 
808  bool isValid() const;
809 
810  // ------------------------- Protected Datas ------------------------------
811  protected:
812 
816  double myError;
817 
822  std::vector <Backpath> myBackpath;
823 
827  Cone myCone;
828 
837 
838 
839 
840  // ------------------------- Private Datas --------------------------------
841  private:
842 
843 
844 
845  // ------------------------- Hidden services ------------------------------
846 
847  public:
848 
855  bool updateBackpath();
856 
862  bool testUpdateBackpath();
863 
868  bool isBackpathOk();
869 
873  void resetBackpath();
874 
878  void resetCone();
879 
884  bool testUpdateWidth();
885 
891  Cone computeNewCone();
892 
897  bool updateWidth();
898 
899 
900 
901 
906  void selfDisplay ( std::ostream & out ) const ;
907 
908 
909 
910  // ------------------------- Internals ------------------------------------
911  private:
912 
913  }; // end of class FrechetShortcut
914 
915 
916  // Utils
917 
918 
925  template <typename TIterator,typename TInteger>
926  std::ostream&
927  operator<< ( std::ostream & out, const FrechetShortcut<TIterator,TInteger> & object );
928 
929 
930 
931 
932 
933 
934 } // namespace DGtal
935 
936 
938 // Includes inline functions.
939 #include "DGtal/geometry/curves/FrechetShortcut.ih"
940 // //
942 
943 #endif // !defined FrechetShortcut_h
944 
945 #undef FrechetShortcut_RECURSES
946 #endif // else defined(FrechetShortcut_RECURSES)
Integer dotProduct(const Vector2I &u, const Vector2I &v) const
static Vector chainCode2Vect(int d)
Aim: On-line computation Computation of the longest shortcut according to the Fr├ęchet distance for a ...
static bool isBetween(double i, double a, double b, double n)
IteratorCirculatorTraits< ConstIterator >::Value Point
MyDigitalSurface::ConstIterator ConstIterator
static double angleVectVect(Vector u, Vector v)
bool operator==(const Self &other) const
ConstIterator end() const
void init(const ConstIterator &it)
Vector::Coordinate Coordinate
void updateBackPathFirstQuad(int d, const ConstIterator &)
BOOST_CONCEPT_ASSERT((concepts::CInteger< TInteger >))
FrechetShortcut< ConstIterator, Integer > Self
Aim: Concept checking for Integer Numbers. More precisely, this concept is a refinement of both CEucl...
Definition: CInteger.h:87
static int computeOctant(Point p, Point q)
struct DGtal::FrechetShortcut::Backpath::occulter_attributes occulter_attributes
FrechetShortcut & operator=(const Self &other)
static int computeChainCode(Point p, Point q)
FrechetShortcut getSelf()
IteratorCirculatorTraits< ConstIterator >::Value Vector
void selfDisplay(std::ostream &out) const
ConstIterator begin() const
Backpath & operator=(const Backpath &other)
static int circleTangentPoints(double x, double y, double x1, double y1, double r1, double *xi, double *yi, double *xi_prime, double *yi_prime)
std::string className() const
bool operator!=(const Self &other) const
DGtal is the top-level namespace which contains all DGtal functions and types.
Cone & operator=(const Cone &c)
static int circle_circle_intersection(double x0, double y0, double r0, double x1, double y1, double r1, double *xi, double *yi, double *xi_prime, double *yi_prime)
Reverse getReverse() const
void selfDisplay(std::ostream &out) const
static int rot(int d, int quad)
FrechetShortcut< std::reverse_iterator< ConstIterator >, Integer > Reverse
const FrechetShortcut< ConstIterator, Integer > * myS
boost::icl::interval_set< double > IntervalSet
static double computeAngle(double x0, double y0, double x1, double y1)
std::vector< Backpath > myBackpath
std::map< ConstIterator, occulter_attributes > occulter_list