FrechetShortcut.h

 
#pragma once
 
#if defined(FrechetShortcut_RECURSES)
#error Recursive header files inclusion detected in FrechetShortcut.h
#else // defined(FrechetShortcut_RECURSES)
#define FrechetShortcut_RECURSES
 
#if !defined FrechetShortcut_h
#define FrechetShortcut_h
 
// Inclusions
#include "DGtal/base/Common.h"
#include "DGtal/base/ConceptUtils.h"
#include "DGtal/kernel/PointVector.h"
#include "DGtal/arithmetic/IntegerComputer.h"
#include <boost/icl/interval_set.hpp>
#include <map>
 
 
#include "DGtal/geometry/curves/SegmentComputerUtils.h"
 
 
#define PRECISION 0.00000001
 
namespace DGtal
{
  
  // class FrechetShortcut
  template <typename TIterator,typename TInteger = typename IteratorCirculatorTraits<TIterator>::Value::Coordinate>
    class FrechetShortcut
    {
      // ----------------------- Standard services ------------------------------
    public:
    
    //entier
    
    BOOST_CONCEPT_ASSERT(( concepts::CInteger<TInteger> ) );
    typedef TInteger Integer;
    
    
    
    //required types
    typedef TIterator ConstIterator;
    typedef FrechetShortcut<ConstIterator,Integer> Self; 
    typedef FrechetShortcut<std::reverse_iterator<ConstIterator>,Integer> Reverse;
    
    //2D point and 2D vector
    typedef typename IteratorCirculatorTraits<ConstIterator>::Value Point; 
    typedef typename IteratorCirculatorTraits<ConstIterator>::Value Vector; 
    typedef typename Vector::Coordinate Coordinate;
    
  class Backpath
  {
  private:
    const FrechetShortcut<ConstIterator,Integer> *myS;
    
  protected: 
    
    typedef struct occulter_attributes{
      double angle_min; // 
      double angle_max; // 
    } occulter_attributes;
    
    typedef std::map <ConstIterator,occulter_attributes > occulter_list;
    
  public:
    friend class FrechetShortcut<ConstIterator,Integer>;
    
    
  public:
    
    typedef boost::icl::interval_set<double> IntervalSet;
    
      int myQuad; 
      
      bool myFlag;   
      
      occulter_list myOcculters;
      
      IntervalSet myForbiddenIntervals;
      
      ConstIterator myIt; 
      
      Backpath();
      
      Backpath(const FrechetShortcut<ConstIterator,Integer> *s ,int q);
 
      Backpath(const Backpath & other);
 
      
      Backpath& operator=(const Backpath & other);
    
 
      ~Backpath();
      
      void reset();
      
      void addPositivePoint();
 
      void addNegativePoint();
 
      void updateBackPathFirstQuad(int d, const ConstIterator&);
      
      void updateOcculters();
 
      void updateIntervals();
      
   
    }; // End of class Backpath
 
 
  class Cone{
    
  public:
    double myMin;
    
    double myMax;
    
    bool myInf; // true if the cone is infinite (the whole plane)
    
    Cone();
    
    Cone(double a0, double a1);
 
    Cone(double x, double y, double x0, double y0, double x1, double y1);
 
    Cone(const Cone& c) = default;
    
    Cone& operator=(const Cone& c);
    
    bool isEmpty() const;
 
    void intersectCones(Cone c);
    
    Cone intersectConesSimple(Cone c);
    
    Cone symmetricalCone();
    
    void selfDisplay ( std::ostream & out) const;
    
    };
    
    
 
    
    struct Tools
    {
      static  bool isBetween(double i, double a, double b, double n)
      {
        if(a<=b)
          return (i>=a && i<=b);
        else
          return ((i>=a && i<=n) || (i>=0 && i<=b));
      }
      
      
      // Code by Tim Voght
      // http://local.wasp.uwa.edu.au/~pbourke/geometry/2circle/tvoght.c
      
      /* circle_circle_intersection() *
       * Determine the points where 2 circles in a common plane intersect.
       *
       * int circle_circle_intersection(
       *                                // center and radius of 1st circle
       *                                double x0, double y0, double r0,
       *                                // center and radius of 2nd circle
       *                                double x1, double y1, double r1,
       *                                // 1st intersection point
       *                                double *xi, double *yi,              
       *                                // 2nd intersection point
       *                                double *xi_prime, double *yi_prime)
       *
       * This is a public domain work. 3/26/2005 Tim Voght
       *
       */
      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)
      {
        // I. Sivignon 05/2024: fix to handle the case where r0=0 or
        // r1=0. Enables the case where error=0. 
        // Other special cases where circles are tangent are treated
        // below. 
        if(r0==0)
          {
            *xi = x0; *yi = y0;
            *xi_prime = x0 ; *yi_prime = y0;
            return 1;
          }
        else
          if(r1==0)
            {
              *xi = x1; *yi = y1;
              *xi_prime = x1 ; *yi_prime = y1;
              return 1;
            }
        
        double a, dx, dy, d, h, rx, ry;
        double x2, y2;
        
        /* dx and dy are the vertical and horizontal distances between
         * the circle centers.
         */
        dx = x1 - x0;
        dy = y1 - y0;
 
        /* Determine the straight-line distance between the centers. */
        //d = sqrt((dy*dy) + (dx*dx));
        d = hypot(dx,dy); // Suggested by Keith Briggs
 
        if((r0+r1)*(r0+r1)-((dy*dy)+(dx*dx)) < PRECISION)   
          {  // I. Sivignon 05/2024: fix to handle the case where
             // circles are tangent but radii are non zero.
            
            double alpha= r0/d;
            *xi = x0 + dx*alpha;
            *yi = y0 + dy*alpha;
            *xi_prime = *xi; *yi_prime = *yi; 
            return 1;
          }
        
        /* Check for solvability. */
        if (d > (r0 + r1))
        {
          /* no solution. circles do not intersect. */
          std::cerr  << "Warning : the two circles do not intersect -> should never happen" << std::endl;
          return 0; //I. Sivignon 05/2010 should never happen for our specific use.
        }
        if (d < fabs(r0 - r1))
          {
            /* no solution. one circle is contained in the other */
            return 0;
          }
 
        //   }
        /* 'point 2' is the point where the line through the circle
         * intersection points crosses the line between the circle
         * centers.  
         */
        
        /* Determine the distance from point 0 to point 2. */
        a = ((r0*r0) - (r1*r1) + (d*d)) / (2.0 * d) ;
 
        /* Determine the coordinates of point 2. */
        x2 = x0 + (dx * a/d);
        y2 = y0 + (dy * a/d);
 
 
        /* Determine the distance from point 2 to either of the
       * intersection points.
       */
        h = sqrt((r0*r0) - (a*a));
 
        /* Now determine the offsets of the intersection points from
         * point 2.
         */
        rx = -dy * (h/d);
        ry = dx * (h/d);
        
        /* Determine the absolute intersection points. */
        *xi = x2 + rx;
        *xi_prime = x2 - rx;
        *yi = y2 + ry;
        *yi_prime = y2 - ry;
        
        return 1;
      }
      
 
 
      static int circleTangentPoints(double x, double y, double x1, double y1, double r1, double *xi, double *yi, 
                              double *xi_prime, double *yi_prime)
      {
        double x0 = (x+x1)/2;
        double y0 = (y+y1)/2;
        double r0 = sqrt((x-x1)*(x-x1) + (y-y1)*(y-y1))/2;
        
        int res =
          circle_circle_intersection(x0,y0,r0,x1,y1,r1,xi,yi,xi_prime,yi_prime);
 
        
        return res;
        
      }
    
 
      static double computeAngle(double x0, double y0, double x1, double y1)
      {
        double x = x1-x0;
        double y = y1-y0;
        
        if(x!=0)
          {
            double alpha = y/x;
          
            if(x>0 && y>=0)
              return atan(alpha);
            else
              if(x>0 && y<0)
                return atan(alpha)+2*M_PI;
              else
                if(x<0)
                  return atan(alpha)+M_PI;
          }
        else
          {
            if(y==0)
              return 0;
            else 
              if(y>0)
                return M_PI_2;
              else
                return 3*M_PI_2;
          }
        return -1;
      }      
      
      
      
      static double angleVectVect(Vector u, Vector v)
      {
        IntegerComputer<Integer> ic;
        return acos((double)ic.dotProduct(u,v)/(u.norm()*v.norm()));
      }
        
      static  int computeChainCode(Point p, Point q)
      {
        int d;
        Coordinate x2 = q[0];
        Coordinate y2 = q[1];
        Coordinate x1 = p[0];
        Coordinate y1 = p[1];
        
        if(x2-x1==0)
          if(y2-y1==1)
            d=2;
          else
            d=6;
        else
          if(x2-x1==1)
            if(y2-y1==0)
              d=0;
            else
              if(y2-y1==1)
                d=1;
              else
                d=7;
          else
            if(y2-y1==0)
              d=4;
            else
              if(y2-y1==1)
                d=3;
              else
                d=5;
        return d;
      }
    
    static int computeOctant(Point p, Point q)
    {
      int d = 0;
      Coordinate x = q[0]-p[0];
      Coordinate y = q[1]-p[1];
      
      if(x>=0)
        {
          if(y>=0)
            {
              if(x>y)
                d=0; // 0 <= y < x  
              else
                if(x!=0)
                  d=1; // 0 <= x <= y
            }
          else
            {
              if(x>=abs(y)) 
                d=7; // 0 < abs(y) <= x 
              else
                d=6; // 0 <= x < abs(y)
            }
        }
      if(x<=0)
        {
          if(y>0)
            if(abs(x)>=y)
              d=3; // 
            else
              d=2;
          else
            {
              if(abs(x)>abs(y))
                d=4;
              else
                if(x!=0)
                  d=5;
            }
        }
      return d;
      
    }
    
      static int rot(int d, int quad)
      {
        return (d-quad+8)%8;
      }
      
    static Vector chainCode2Vect(int d)
    {
      Vector v;
      switch(d){
        
      case 0:{
        v[0] = 1;
        v[1] = 0;
        break;
      }
      case 1:{
        v[0] = 1;
        v[1] = 1;
        break;
      }
      case 2:{
        v[0] = 0;
        v[1] = 1;
        break;
      }
      case 3:{
        v[0] = -1;
        v[1] = 1;
        break;
      }
      case 4:{
        v[0] = -1;
        v[1] = 0;
        break;
      }
      case 5:{
        v[0] = -1;
        v[1] = -1;
        break;
      }
      case 6:{
        v[0] = 0;
        v[1] = -1;
        break;
      }
      case 7:{
        v[0] = 1;
        v[1] = -1;
        break;
      } 
        
      }
      
      return v;
    }
    
    
  };
  
    
  FrechetShortcut();
  
  FrechetShortcut(double error);
 
  void init(const ConstIterator& it);
  
  FrechetShortcut getSelf();
 
  FrechetShortcut ( const Self & other );
  
  FrechetShortcut & operator= ( const Self & other );
 
  Reverse getReverse() const;
  
  bool operator==( const Self & other ) const;
    
  bool operator!=( const Self & other ) const;
  
 
  
  ~FrechetShortcut(){};
 
    // ----------------------- Interface --------------------------------------
public:
 
  bool isExtendableFront();
  
  bool extendFront();
    
  
  // ---------------------------- Accessors ----------------------------------
  
  
  ConstIterator begin() const;
  ConstIterator end() const;
  
 
  public:  
  
  std::string className() const;
  
  
  bool isValid() const;
  
  // ------------------------- Protected Datas ------------------------------
 protected:
  
    double myError;
    
    std::vector <Backpath> myBackpath;
  
    Cone myCone;
  
  ConstIterator myBegin;
  ConstIterator myEnd;
  
 
  
  // ------------------------- Private Datas --------------------------------
 private:
 
    
  
  // ------------------------- Hidden services ------------------------------
 
 public:
 
  bool updateBackpath(); 
 
  bool testUpdateBackpath();
  
  bool isBackpathOk();
  
  void resetBackpath();
  
  void resetCone();
  
  bool testUpdateWidth();
  
  Cone computeNewCone();
  
  bool updateWidth();
  
    
 
    
  void selfDisplay ( std::ostream & out ) const ; 
  
  
  
  // ------------------------- Internals ------------------------------------
 private:
  
  }; // end of class FrechetShortcut
  
 
  // Utils
  
  
  template <typename TIterator,typename TInteger>
  std::ostream&
    operator<< ( std::ostream & out, const FrechetShortcut<TIterator,TInteger> & object );
 
 
  
  
  
  
} // namespace DGtal
 
 
// Includes inline functions.
#include "DGtal/geometry/curves/FrechetShortcut.ih"
//                                                                           //
 
#endif // !defined FrechetShortcut_h
 
#undef FrechetShortcut_RECURSES
#endif // else defined(FrechetShortcut_RECURSES)