DGtal  1.5.beta
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)
36 #define FrechetShortcut_RECURSES
37 
38 #if !defined FrechetShortcut_h
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 
56 #define PRECISION 0.00000001
57 
58 namespace DGtal
59 {
60 
62  // class FrechetShortcut
109  template <typename TIterator,typename TInteger = typename IteratorCirculatorTraits<TIterator>::Value::Coordinate>
111  {
112  // ----------------------- Standard services ------------------------------
113  public:
114 
115  //entier
116 
118  typedef TInteger Integer;
119 
120 
121 
122  //required types
123  typedef TIterator ConstIterator;
126 
127  //2D point and 2D vector
130  typedef typename Vector::Coordinate Coordinate;
131 
136  class Backpath
137  {
138  private:
143 
144  protected:
145 
150  typedef struct occulter_attributes{
151  double angle_min; //
152  double angle_max; //
154 
159  typedef std::map <ConstIterator,occulter_attributes > occulter_list;
160 
161  public:
162  friend class FrechetShortcut<ConstIterator,Integer>;
163 
164 
165  public:
166 
167  typedef boost::icl::interval_set<double> IntervalSet;
168 
172  int myQuad;
173 
178  bool myFlag;
179 
181 
187 
192 
197 
204 
209  Backpath(const Backpath & other);
210 
211 
217  Backpath& operator=(const Backpath & other);
218 
219 
224 
228  void reset();
229 
234 
239 
247 
252 
257 
258 
259  }; // End of class Backpath
260 
261 
266  class Cone{
267 
268  public:
272  double myMin;
273 
277  double myMax;
278 
282  bool myInf; // true if the cone is infinite (the whole plane)
283 
284  Cone();
285 
291  Cone(double a0, double a1);
292 
303  Cone(double x, double y, double x0, double y0, double x1, double y1);
304 
309  Cone(const Cone& c) = default;
310 
316  Cone& operator=(const Cone& c);
317 
322  bool isEmpty() const;
323 
329 
336 
342 
347  void selfDisplay ( std::ostream & out) const;
348 
349  };
350 
351 
352 
353 
354  struct Tools
355  {
365  static bool isBetween(double i, double a, double b, double n)
366  {
367  if(a<=b)
368  return (i>=a && i<=b);
369  else
370  return ((i>=a && i<=n) || (i>=0 && i<=b));
371  }
372 
373 
374  // Code by Tim Voght
375  // http://local.wasp.uwa.edu.au/~pbourke/geometry/2circle/tvoght.c
376 
377  /* circle_circle_intersection() *
378  * Determine the points where 2 circles in a common plane intersect.
379  *
380  * int circle_circle_intersection(
381  * // center and radius of 1st circle
382  * double x0, double y0, double r0,
383  * // center and radius of 2nd circle
384  * double x1, double y1, double r1,
385  * // 1st intersection point
386  * double *xi, double *yi,
387  * // 2nd intersection point
388  * double *xi_prime, double *yi_prime)
389  *
390  * This is a public domain work. 3/26/2005 Tim Voght
391  *
392  */
400  static int circle_circle_intersection(double x0, double y0, double r0,
401  double x1, double y1, double r1,
402  double *xi, double *yi,
403  double *xi_prime, double *yi_prime)
404  {
405  // I. Sivignon 05/2024: fix to handle the case where r0=0 or
406  // r1=0. Enables the case where error=0.
407  // Other special cases where circles are tangent are treated
408  // below.
409  if(r0==0)
410  {
411  *xi = x0; *yi = y0;
412  *xi_prime = x0 ; *yi_prime = y0;
413  return 1;
414  }
415  else
416  if(r1==0)
417  {
418  *xi = x1; *yi = y1;
419  *xi_prime = x1 ; *yi_prime = y1;
420  return 1;
421  }
422 
423  double a, dx, dy, d, h, rx, ry;
424  double x2, y2;
425 
426  /* dx and dy are the vertical and horizontal distances between
427  * the circle centers.
428  */
429  dx = x1 - x0;
430  dy = y1 - y0;
431 
432  /* Determine the straight-line distance between the centers. */
433  //d = sqrt((dy*dy) + (dx*dx));
434  d = hypot(dx,dy); // Suggested by Keith Briggs
435 
436  if((r0+r1)*(r0+r1)-((dy*dy)+(dx*dx)) < PRECISION)
437  { // I. Sivignon 05/2024: fix to handle the case where
438  // circles are tangent but radii are non zero.
439 
440  double alpha= r0/d;
441  *xi = x0 + dx*alpha;
442  *yi = y0 + dy*alpha;
443  *xi_prime = *xi; *yi_prime = *yi;
444  return 1;
445  }
446 
447  /* Check for solvability. */
448  if (d > (r0 + r1))
449  {
450  /* no solution. circles do not intersect. */
451  std::cerr << "Warning : the two circles do not intersect -> should never happen" << std::endl;
452  return 0; //I. Sivignon 05/2010 should never happen for our specific use.
453  }
454  if (d < fabs(r0 - r1))
455  {
456  /* no solution. one circle is contained in the other */
457  return 0;
458  }
459 
460  // }
461  /* 'point 2' is the point where the line through the circle
462  * intersection points crosses the line between the circle
463  * centers.
464  */
465 
466  /* Determine the distance from point 0 to point 2. */
467  a = ((r0*r0) - (r1*r1) + (d*d)) / (2.0 * d) ;
468 
469  /* Determine the coordinates of point 2. */
470  x2 = x0 + (dx * a/d);
471  y2 = y0 + (dy * a/d);
472 
473 
474  /* Determine the distance from point 2 to either of the
475  * intersection points.
476  */
477  h = sqrt((r0*r0) - (a*a));
478 
479  /* Now determine the offsets of the intersection points from
480  * point 2.
481  */
482  rx = -dy * (h/d);
483  ry = dx * (h/d);
484 
485  /* Determine the absolute intersection points. */
486  *xi = x2 + rx;
487  *xi_prime = x2 - rx;
488  *yi = y2 + ry;
489  *yi_prime = y2 - ry;
490 
491  return 1;
492  }
493 
494 
495 
515  static int circleTangentPoints(double x, double y, double x1, double y1, double r1, double *xi, double *yi,
516  double *xi_prime, double *yi_prime)
517  {
518  double x0 = (x+x1)/2;
519  double y0 = (y+y1)/2;
520  double r0 = sqrt((x-x1)*(x-x1) + (y-y1)*(y-y1))/2;
521 
522  int res =
523  circle_circle_intersection(x0,y0,r0,x1,y1,r1,xi,yi,xi_prime,yi_prime);
524 
525 
526  return res;
527 
528  }
529 
530 
540  static double computeAngle(double x0, double y0, double x1, double y1)
541  {
542  double x = x1-x0;
543  double y = y1-y0;
544 
545  if(x!=0)
546  {
547  double alpha = y/x;
548 
549  if(x>0 && y>=0)
550  return atan(alpha);
551  else
552  if(x>0 && y<0)
553  return atan(alpha)+2*M_PI;
554  else
555  if(x<0)
556  return atan(alpha)+M_PI;
557  }
558  else
559  {
560  if(y==0)
561  return 0;
562  else
563  if(y>0)
564  return M_PI_2;
565  else
566  return 3*M_PI_2;
567  }
568  return -1;
569  }
570 
571 
572 
578  static double angleVectVect(Vector u, Vector v)
579  {
581  return acos((double)ic.dotProduct(u,v)/(u.norm()*v.norm()));
582  }
583 
589  static int computeChainCode(Point p, Point q)
590  {
591  int d;
592  Coordinate x2 = q[0];
593  Coordinate y2 = q[1];
594  Coordinate x1 = p[0];
595  Coordinate y1 = p[1];
596 
597  if(x2-x1==0)
598  if(y2-y1==1)
599  d=2;
600  else
601  d=6;
602  else
603  if(x2-x1==1)
604  if(y2-y1==0)
605  d=0;
606  else
607  if(y2-y1==1)
608  d=1;
609  else
610  d=7;
611  else
612  if(y2-y1==0)
613  d=4;
614  else
615  if(y2-y1==1)
616  d=3;
617  else
618  d=5;
619  return d;
620  }
621 
627  static int computeOctant(Point p, Point q)
628  {
629  int d = 0;
630  Coordinate x = q[0]-p[0];
631  Coordinate y = q[1]-p[1];
632 
633  if(x>=0)
634  {
635  if(y>=0)
636  {
637  if(x>y)
638  d=0; // 0 <= y < x
639  else
640  if(x!=0)
641  d=1; // 0 <= x <= y
642  }
643  else
644  {
645  if(x>=abs(y))
646  d=7; // 0 < abs(y) <= x
647  else
648  d=6; // 0 <= x < abs(y)
649  }
650  }
651  if(x<=0)
652  {
653  if(y>0)
654  if(abs(x)>=y)
655  d=3; //
656  else
657  d=2;
658  else
659  {
660  if(abs(x)>abs(y))
661  d=4;
662  else
663  if(x!=0)
664  d=5;
665  }
666  }
667  return d;
668 
669  }
670 
679  static int rot(int d, int quad)
680  {
681  return (d-quad+8)%8;
682  }
683 
689  static Vector chainCode2Vect(int d)
690  {
691  Vector v;
692  switch(d){
693 
694  case 0:{
695  v[0] = 1;
696  v[1] = 0;
697  break;
698  }
699  case 1:{
700  v[0] = 1;
701  v[1] = 1;
702  break;
703  }
704  case 2:{
705  v[0] = 0;
706  v[1] = 1;
707  break;
708  }
709  case 3:{
710  v[0] = -1;
711  v[1] = 1;
712  break;
713  }
714  case 4:{
715  v[0] = -1;
716  v[1] = 0;
717  break;
718  }
719  case 5:{
720  v[0] = -1;
721  v[1] = -1;
722  break;
723  }
724  case 6:{
725  v[0] = 0;
726  v[1] = -1;
727  break;
728  }
729  case 7:{
730  v[0] = 1;
731  v[1] = -1;
732  break;
733  }
734 
735  }
736 
737  return v;
738  }
739 
740 
741  };
742 
743 
749 
754  FrechetShortcut(double error);
755 
761  void init(const ConstIterator& it);
762 
764 
769  FrechetShortcut ( const Self & other );
770 
776  FrechetShortcut & operator= ( const Self & other );
777 
782 
789  bool operator==( const Self & other ) const;
790 
797  bool operator!=( const Self & other ) const;
798 
799 
800 
805 
806  // ----------------------- Interface --------------------------------------
807 public:
808 
815 
821  bool extendFront();
822 
823 
824  // ---------------------------- Accessors ----------------------------------
825 
826 
836 
837 
838  public:
839 
843  std::string className() const;
844 
845 
850  bool isValid() const;
851 
852  // ------------------------- Protected Datas ------------------------------
853  protected:
854 
858  double myError;
859 
864  std::vector <Backpath> myBackpath;
865 
870 
879 
880 
881 
882  // ------------------------- Private Datas --------------------------------
883  private:
884 
885 
886 
887  // ------------------------- Hidden services ------------------------------
888 
889  public:
890 
897  bool updateBackpath();
898 
905 
910  bool isBackpathOk();
911 
916 
920  void resetCone();
921 
927 
934 
939  bool updateWidth();
940 
941 
942 
943 
948  void selfDisplay ( std::ostream & out ) const ;
949 
950 
951 
952  // ------------------------- Internals ------------------------------------
953  private:
954 
955  }; // end of class FrechetShortcut
956 
957 
958  // Utils
959 
960 
967  template <typename TIterator,typename TInteger>
968  std::ostream&
969  operator<< ( std::ostream & out, const FrechetShortcut<TIterator,TInteger> & object );
970 
971 
972 
973 
974 
975 
976 } // namespace DGtal
977 
978 
980 // Includes inline functions.
981 #include "DGtal/geometry/curves/FrechetShortcut.ih"
982 // //
984 
985 #endif // !defined FrechetShortcut_h
986 
987 #undef FrechetShortcut_RECURSES
988 #endif // else defined(FrechetShortcut_RECURSES)
const FrechetShortcut< ConstIterator, Integer > * myS
struct DGtal::FrechetShortcut::Backpath::occulter_attributes occulter_attributes
boost::icl::interval_set< double > IntervalSet
void updateBackPathFirstQuad(int d, const ConstIterator &)
Backpath & operator=(const Backpath &other)
std::map< ConstIterator, occulter_attributes > occulter_list
Backpath(const Backpath &other)
Backpath(const FrechetShortcut< ConstIterator, Integer > *s, int q)
void selfDisplay(std::ostream &out) const
Cone(const Cone &c)=default
Cone & operator=(const Cone &c)
Cone(double a0, double a1)
Cone(double x, double y, double x0, double y0, double x1, double y1)
Aim: On-line computation Computation of the longest shortcut according to the Fr├ęchet distance for a ...
std::string className() const
FrechetShortcut< std::reverse_iterator< ConstIterator >, Integer > Reverse
ConstIterator begin() const
IteratorCirculatorTraits< ConstIterator >::Value Point
bool operator!=(const Self &other) const
Vector::Coordinate Coordinate
Reverse getReverse() const
IteratorCirculatorTraits< ConstIterator >::Value Vector
FrechetShortcut & operator=(const Self &other)
FrechetShortcut(const Self &other)
void selfDisplay(std::ostream &out) const
void init(const ConstIterator &it)
FrechetShortcut(double error)
FrechetShortcut getSelf()
BOOST_CONCEPT_ASSERT((concepts::CInteger< TInteger >))
bool operator==(const Self &other) const
ConstIterator end() const
std::vector< Backpath > myBackpath
FrechetShortcut< ConstIterator, Integer > Self
Integer dotProduct(const Vector2I &u, const Vector2I &v) const
DGtal is the top-level namespace which contains all DGtal functions and types.
std::ostream & operator<<(std::ostream &out, const ATu0v1< TKSpace, TLinearAlgebra > &object)
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)
static double angleVectVect(Vector u, Vector v)
static Vector chainCode2Vect(int d)
static int circleTangentPoints(double x, double y, double x1, double y1, double r1, double *xi, double *yi, double *xi_prime, double *yi_prime)
static bool isBetween(double i, double a, double b, double n)
static double computeAngle(double x0, double y0, double x1, double y1)
static int rot(int d, int quad)
static int computeChainCode(Point p, Point q)
static int computeOctant(Point p, Point q)
Aim: Concept checking for Integer Numbers. More precisely, this concept is a refinement of both CEucl...
Definition: CInteger.h:88