DGtal  1.4.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 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 
178 
184 
189 
194 
201 
206  Backpath(const Backpath & other);
207 
208 
214  Backpath& operator=(const Backpath & other);
215 
216 
221 
225  void reset();
226 
231 
236 
244 
249 
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  Cone(const Cone& c) = default;
307 
313  Cone& operator=(const Cone& c);
314 
319  bool isEmpty() const;
320 
326 
333 
339 
344  void selfDisplay ( std::ostream & out) const;
345 
346  };
347 
348 
349 
350 
351  struct Tools
352  {
362  static bool isBetween(double i, double a, double b, double n)
363  {
364  if(a<=b)
365  return (i>=a && i<=b);
366  else
367  return ((i>=a && i<=n) || (i>=0 && i<=b));
368  }
369 
370 
371  // Code by Tim Voght
372  // http://local.wasp.uwa.edu.au/~pbourke/geometry/2circle/tvoght.c
373 
374  /* circle_circle_intersection() *
375  * Determine the points where 2 circles in a common plane intersect.
376  *
377  * int circle_circle_intersection(
378  * // center and radius of 1st circle
379  * double x0, double y0, double r0,
380  * // center and radius of 2nd circle
381  * double x1, double y1, double r1,
382  * // 1st intersection point
383  * double *xi, double *yi,
384  * // 2nd intersection point
385  * double *xi_prime, double *yi_prime)
386  *
387  * This is a public domain work. 3/26/2005 Tim Voght
388  *
389  */
397  static int circle_circle_intersection(double x0, double y0, double r0,
398  double x1, double y1, double r1,
399  double *xi, double *yi,
400  double *xi_prime, double *yi_prime)
401  {
402  double a, dx, dy, d, h, rx, ry;
403  double x2, y2;
404 
405  /* dx and dy are the vertical and horizontal distances between
406  * the circle centers.
407  */
408  dx = x1 - x0;
409  dy = y1 - y0;
410 
411  /* Determine the straight-line distance between the centers. */
412  //d = sqrt((dy*dy) + (dx*dx));
413  d = hypot(dx,dy); // Suggested by Keith Briggs
414 
415  /* Check for solvability. */
416  if (d > (r0 + r1))
417  {
418  /* no solution. circles do not intersect. */
419  std::cerr << "Warning : the two circles do not intersect -> should never happen" << std::endl;
420  return 0; //I. Sivignon 05/2010 should never happen for our specific use.
421  }
422  if (d < fabs(r0 - r1))
423  {
424  /* no solution. one circle is contained in the other */
425  return 0;
426  }
427 
428  /* 'point 2' is the point where the line through the circle
429  * intersection points crosses the line between the circle
430  * centers.
431  */
432 
433  /* Determine the distance from point 0 to point 2. */
434  a = ((r0*r0) - (r1*r1) + (d*d)) / (2.0 * d) ;
435 
436  /* Determine the coordinates of point 2. */
437  x2 = x0 + (dx * a/d);
438  y2 = y0 + (dy * a/d);
439 
440  /* Determine the distance from point 2 to either of the
441  * intersection points.
442  */
443  h = sqrt((r0*r0) - (a*a));
444 
445  /* Now determine the offsets of the intersection points from
446  * point 2.
447  */
448  rx = -dy * (h/d);
449  ry = dx * (h/d);
450 
451  /* Determine the absolute intersection points. */
452  *xi = x2 + rx;
453  *xi_prime = x2 - rx;
454  *yi = y2 + ry;
455  *yi_prime = y2 - ry;
456 
457  return 1;
458  }
459 
460 
461 
481  static int circleTangentPoints(double x, double y, double x1, double y1, double r1, double *xi, double *yi,
482  double *xi_prime, double *yi_prime)
483  {
484  double x0 = (x+x1)/2;
485  double y0 = (y+y1)/2;
486  double r0 = sqrt((x-x1)*(x-x1) + (y-y1)*(y-y1))/2;
487 
488  int res =
489  circle_circle_intersection(x0,y0,r0,x1,y1,r1,xi,yi,xi_prime,yi_prime);
490 
491  return res;
492 
493  }
494 
495 
505  static double computeAngle(double x0, double y0, double x1, double y1)
506  {
507  double x = x1-x0;
508  double y = y1-y0;
509 
510  if(x!=0)
511  {
512  double alpha = y/x;
513 
514  if(x>0 && y>=0)
515  return atan(alpha);
516  else
517  if(x>0 && y<0)
518  return atan(alpha)+2*M_PI;
519  else
520  if(x<0)
521  return atan(alpha)+M_PI;
522  }
523  else
524  {
525  if(y>0)
526  return M_PI_2;
527  else
528  return 3*M_PI_2;
529  }
530  return -1;
531  }
532 
533 
534 
540  static double angleVectVect(Vector u, Vector v)
541  {
543  return acos((double)ic.dotProduct(u,v)/(u.norm()*v.norm()));
544  }
545 
551  static int computeChainCode(Point p, Point q)
552  {
553  int d;
554  Coordinate x2 = q[0];
555  Coordinate y2 = q[1];
556  Coordinate x1 = p[0];
557  Coordinate y1 = p[1];
558 
559  if(x2-x1==0)
560  if(y2-y1==1)
561  d=2;
562  else
563  d=6;
564  else
565  if(x2-x1==1)
566  if(y2-y1==0)
567  d=0;
568  else
569  if(y2-y1==1)
570  d=1;
571  else
572  d=7;
573  else
574  if(y2-y1==0)
575  d=4;
576  else
577  if(y2-y1==1)
578  d=3;
579  else
580  d=5;
581  return d;
582  }
583 
589  static int computeOctant(Point p, Point q)
590  {
591  int d = 0;
592  Coordinate x = q[0]-p[0];
593  Coordinate y = q[1]-p[1];
594 
595  if(x>=0)
596  {
597  if(y>=0)
598  {
599  if(x>y)
600  d=0; // 0 <= y < x
601  else
602  if(x!=0)
603  d=1; // 0 <= x <= y
604  }
605  else
606  {
607  if(x>=abs(y))
608  d=7; // 0 < abs(y) <= x
609  else
610  d=6; // 0 <= x < abs(y)
611  }
612  }
613  if(x<=0)
614  {
615  if(y>0)
616  if(abs(x)>=y)
617  d=3; //
618  else
619  d=2;
620  else
621  {
622  if(abs(x)>abs(y))
623  d=4;
624  else
625  if(x!=0)
626  d=5;
627  }
628  }
629  return d;
630 
631  }
632 
641  static int rot(int d, int quad)
642  {
643  return (d-quad+8)%8;
644  }
645 
651  static Vector chainCode2Vect(int d)
652  {
653  Vector v;
654  switch(d){
655 
656  case 0:{
657  v[0] = 1;
658  v[1] = 0;
659  break;
660  }
661  case 1:{
662  v[0] = 1;
663  v[1] = 1;
664  break;
665  }
666  case 2:{
667  v[0] = 0;
668  v[1] = 1;
669  break;
670  }
671  case 3:{
672  v[0] = -1;
673  v[1] = 1;
674  break;
675  }
676  case 4:{
677  v[0] = -1;
678  v[1] = 0;
679  break;
680  }
681  case 5:{
682  v[0] = -1;
683  v[1] = -1;
684  break;
685  }
686  case 6:{
687  v[0] = 0;
688  v[1] = -1;
689  break;
690  }
691  case 7:{
692  v[0] = 1;
693  v[1] = -1;
694  break;
695  }
696 
697  }
698 
699  return v;
700  }
701 
702 
703  };
704 
705 
711 
716  FrechetShortcut(double error);
717 
723  void init(const ConstIterator& it);
724 
726 
731  FrechetShortcut ( const Self & other );
732 
738  FrechetShortcut & operator= ( const Self & other );
739 
744 
751  bool operator==( const Self & other ) const;
752 
759  bool operator!=( const Self & other ) const;
760 
761 
762 
767 
768  // ----------------------- Interface --------------------------------------
769 public:
770 
777 
783  bool extendFront();
784 
785 
786  // ---------------------------- Accessors ----------------------------------
787 
788 
798 
799 
800 
801 
802  public:
803 
807  std::string className() const;
808 
809 
814  bool isValid() const;
815 
816  // ------------------------- Protected Datas ------------------------------
817  protected:
818 
822  double myError;
823 
828  std::vector <Backpath> myBackpath;
829 
834 
843 
844 
845 
846  // ------------------------- Private Datas --------------------------------
847  private:
848 
849 
850 
851  // ------------------------- Hidden services ------------------------------
852 
853  public:
854 
861  bool updateBackpath();
862 
869 
874  bool isBackpathOk();
875 
880 
884  void resetCone();
885 
891 
898 
903  bool updateWidth();
904 
905 
906 
907 
912  void selfDisplay ( std::ostream & out ) const ;
913 
914 
915 
916  // ------------------------- Internals ------------------------------------
917  private:
918 
919  }; // end of class FrechetShortcut
920 
921 
922  // Utils
923 
924 
931  template <typename TIterator,typename TInteger>
932  std::ostream&
933  operator<< ( std::ostream & out, const FrechetShortcut<TIterator,TInteger> & object );
934 
935 
936 
937 
938 
939 
940 } // namespace DGtal
941 
942 
944 // Includes inline functions.
945 #include "DGtal/geometry/curves/FrechetShortcut.ih"
946 // //
948 
949 #endif // !defined FrechetShortcut_h
950 
951 #undef FrechetShortcut_RECURSES
952 #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