DGtal  0.9.2
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:
159  friend class FrechetShortcut<ConstIterator,Integer>;
160 
161 
162  public:
163 
164  typedef boost::icl::interval_set<double> IntervalSet;
165 
169  int myQuad;
170 
175  bool myFlag;
176 
177  occulter_list myOcculters;
178 
183  IntervalSet myForbiddenIntervals;
184 
188  ConstIterator myIt;
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 
346  struct Tools
347  {
357  static bool isBetween(double i, double a, double b, double n)
358  {
359  if(a<=b)
360  if(i>=a && i<=b)
361  return true;
362  else
363  return false;
364  else
365  {
366  if((i>=a && i<=n) || (i>=0 && i<=b))
367  return true;
368  else
369  return false;
370 
371  }
372  }
373 
374 
375  // Code by Tim Voght
376  // http://local.wasp.uwa.edu.au/~pbourke/geometry/2circle/tvoght.c
377 
378  /* circle_circle_intersection() *
379  * Determine the points where 2 circles in a common plane intersect.
380  *
381  * int circle_circle_intersection(
382  * // center and radius of 1st circle
383  * double x0, double y0, double r0,
384  * // center and radius of 2nd circle
385  * double x1, double y1, double r1,
386  * // 1st intersection point
387  * double *xi, double *yi,
388  * // 2nd intersection point
389  * double *xi_prime, double *yi_prime)
390  *
391  * This is a public domain work. 3/26/2005 Tim Voght
392  *
393  */
401  static int circle_circle_intersection(double x0, double y0, double r0,
402  double x1, double y1, double r1,
403  double *xi, double *yi,
404  double *xi_prime, double *yi_prime)
405  {
406  double a, dx, dy, d, h, rx, ry;
407  double x2, y2;
408 
409  /* dx and dy are the vertical and horizontal distances between
410  * the circle centers.
411  */
412  dx = x1 - x0;
413  dy = y1 - y0;
414 
415  /* Determine the straight-line distance between the centers. */
416  //d = sqrt((dy*dy) + (dx*dx));
417  d = hypot(dx,dy); // Suggested by Keith Briggs
418 
419  /* Check for solvability. */
420  if (d > (r0 + r1))
421  {
422  /* no solution. circles do not intersect. */
423  std::cerr << "Warning : the two circles do not intersect -> should never happen" << std::endl;
424  return 0; //I. Sivignon 05/2010 should never happen for our specific use.
425  }
426  if (d < fabs(r0 - r1))
427  {
428  /* no solution. one circle is contained in the other */
429  return 0;
430  }
431 
432  /* 'point 2' is the point where the line through the circle
433  * intersection points crosses the line between the circle
434  * centers.
435  */
436 
437  /* Determine the distance from point 0 to point 2. */
438  a = ((r0*r0) - (r1*r1) + (d*d)) / (2.0 * d) ;
439 
440  /* Determine the coordinates of point 2. */
441  x2 = x0 + (dx * a/d);
442  y2 = y0 + (dy * a/d);
443 
444  /* Determine the distance from point 2 to either of the
445  * intersection points.
446  */
447  h = sqrt((r0*r0) - (a*a));
448 
449  /* Now determine the offsets of the intersection points from
450  * point 2.
451  */
452  rx = -dy * (h/d);
453  ry = dx * (h/d);
454 
455  /* Determine the absolute intersection points. */
456  *xi = x2 + rx;
457  *xi_prime = x2 - rx;
458  *yi = y2 + ry;
459  *yi_prime = y2 - ry;
460 
461  return 1;
462  }
463 
464 
465 
485  static int circleTangentPoints(double x, double y, double x1, double y1, double r1, double *xi, double *yi,
486  double *xi_prime, double *yi_prime)
487  {
488  double x0 = (x+x1)/2;
489  double y0 = (y+y1)/2;
490  double r0 = sqrt((x-x1)*(x-x1) + (y-y1)*(y-y1))/2;
491 
492  int res =
493  circle_circle_intersection(x0,y0,r0,x1,y1,r1,xi,yi,xi_prime,yi_prime);
494 
495  return res;
496 
497  }
498 
499 
509  static double computeAngle(double x0, double y0, double x1, double y1)
510  {
511  double x = x1-x0;
512  double y = y1-y0;
513 
514  if(x!=0)
515  {
516  double alpha = y/x;
517 
518  if(x>0 && y>=0)
519  return atan(alpha);
520  else
521  if(x>0 && y<0)
522  return atan(alpha)+2*M_PI;
523  else
524  if(x<0)
525  return atan(alpha)+M_PI;
526  }
527  else
528  {
529  if(y>0)
530  return M_PI_2;
531  else
532  return 3*M_PI_2;
533  }
534  return -1;
535  }
536 
537 
538 
544  static double angleVectVect(Vector u, Vector v)
545  {
547  return acos((double)ic.dotProduct(u,v)/(u.norm()*v.norm()));
548  }
549 
555  static int computeChainCode(Point p, Point q)
556  {
557  int d;
558  Coordinate x2 = q[0];
559  Coordinate y2 = q[1];
560  Coordinate x1 = p[0];
561  Coordinate y1 = p[1];
562 
563  if(x2-x1==0)
564  if(y2-y1==1)
565  d=2;
566  else
567  d=6;
568  else
569  if(x2-x1==1)
570  if(y2-y1==0)
571  d=0;
572  else
573  if(y2-y1==1)
574  d=1;
575  else
576  d=7;
577  else
578  if(y2-y1==0)
579  d=4;
580  else
581  if(y2-y1==1)
582  d=3;
583  else
584  d=5;
585  return d;
586  }
587 
593  static int computeOctant(Point p, Point q)
594  {
595  int d = 0;
596  Coordinate x = q[0]-p[0];
597  Coordinate y = q[1]-p[1];
598 
599  if(x>=0)
600  {
601  if(y>=0)
602  {
603  if(x>y)
604  d=0; // 0 <= y < x
605  else
606  if(x!=0)
607  d=1; // 0 <= x <= y
608  }
609  else
610  {
611  if(x>=abs(y))
612  d=7; // 0 < abs(y) <= x
613  else
614  d=6; // 0 <= x < abs(y)
615  }
616  }
617  if(x<=0)
618  {
619  if(y>0)
620  if(abs(x)>=y)
621  d=3; //
622  else
623  d=2;
624  else
625  {
626  if(abs(x)>abs(y))
627  d=4;
628  else
629  if(x!=0)
630  d=5;
631  }
632  }
633  return d;
634 
635  }
636 
645  static int rot(int d, int quad)
646  {
647  return (d-quad+8)%8;
648  }
649 
655  static Vector chainCode2Vect(int d)
656  {
657  Vector v;
658  switch(d){
659 
660  case 0:{
661  v[0] = 1;
662  v[1] = 0;
663  break;
664  }
665  case 1:{
666  v[0] = 1;
667  v[1] = 1;
668  break;
669  }
670  case 2:{
671  v[0] = 0;
672  v[1] = 1;
673  break;
674  }
675  case 3:{
676  v[0] = -1;
677  v[1] = 1;
678  break;
679  }
680  case 4:{
681  v[0] = -1;
682  v[1] = 0;
683  break;
684  }
685  case 5:{
686  v[0] = -1;
687  v[1] = -1;
688  break;
689  }
690  case 6:{
691  v[0] = 0;
692  v[1] = -1;
693  break;
694  }
695  case 7:{
696  v[0] = 1;
697  v[1] = -1;
698  break;
699  }
700 
701  }
702 
703  return v;
704  }
705 
706 
707  };
708 
709 
714  FrechetShortcut();
715 
720  FrechetShortcut(double error);
721 
727  void init(const ConstIterator& it);
728 
730 
735  FrechetShortcut ( const Self & other );
736 
742  FrechetShortcut & operator= ( const Self & other );
743 
747  Reverse getReverse() const;
748 
755  bool operator==( const Self & other ) const;
756 
763  bool operator!=( const Self & other ) const;
764 
765 
766 
771 
772  // ----------------------- Interface --------------------------------------
773 public:
774 
780  bool isExtendableFront();
781 
787  bool extendFront();
788 
789 
790  // ---------------------------- Accessors ----------------------------------
791 
792 
797  ConstIterator begin() const;
801  ConstIterator end() const;
802 
803 
804 
805 
806  public:
807 
811  std::string className() const;
812 
813 
818  bool isValid() const;
819 
820  // ------------------------- Protected Datas ------------------------------
821  protected:
822 
826  double myError;
827 
832  std::vector <Backpath> myBackpath;
833 
837  Cone myCone;
838 
842  ConstIterator myBegin;
846  ConstIterator myEnd;
847 
848 
849 
850  // ------------------------- Private Datas --------------------------------
851  private:
852 
853 
854 
855  // ------------------------- Hidden services ------------------------------
856 
857  public:
858 
865  bool updateBackpath();
866 
872  bool testUpdateBackpath();
873 
878  bool isBackpathOk();
879 
883  void resetBackpath();
884 
888  void resetCone();
889 
894  bool testUpdateWidth();
895 
901  Cone computeNewCone();
902 
907  bool updateWidth();
908 
909 
910 
911 
916  void selfDisplay ( std::ostream & out ) const ;
917 
918 
919 
920  // ------------------------- Internals ------------------------------------
921  private:
922 
923  }; // end of class FrechetShortcut
924 
925 
926  // Utils
927 
928 
935  template <typename TIterator,typename TInteger>
936  std::ostream&
937  operator<< ( std::ostream & out, const FrechetShortcut<TIterator,TInteger> & object );
938 
939 
940 
941 
942 
943 
944 } // namespace DGtal
945 
946 
948 // Includes inline functions.
949 #include "DGtal/geometry/curves/FrechetShortcut.ih"
950 // //
952 
953 #endif // !defined FrechetShortcut_h
954 
955 #undef FrechetShortcut_RECURSES
956 #endif // else defined(FrechetShortcut_RECURSES)
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
static double angleVectVect(Vector u, Vector v)
void selfDisplay(std::ostream &out) const
bool operator!=(const Self &other) const
void init(const ConstIterator &it)
Vector::Coordinate Coordinate
void updateBackPathFirstQuad(int d, const ConstIterator &)
BOOST_CONCEPT_ASSERT((concepts::CInteger< TInteger >))
Reverse getReverse() const
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()
ConstIterator begin() const
IteratorCirculatorTraits< ConstIterator >::Value Vector
std::string className() 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)
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)
Integer dotProduct(const Vector2I &u, const Vector2I &v) const
ConstIterator end() const
bool operator==(const Self &other) const
static int rot(int d, int quad)
FrechetShortcut< std::reverse_iterator< ConstIterator >, Integer > Reverse
bool isValid() const
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
void selfDisplay(std::ostream &out) const
std::map< ConstIterator, occulter_attributes > occulter_list