eic-opticks/ana/ellipse_intersect_circle.hh

0001 #pragma once
0002
0003 /**
0004 ellipse_intersect_circle.hh
0005 -----------------------------
0006
0007 Numerical intersection of an ellipse and a circle
0008 returning 0,1,2,3 or 4 intersection points.
0009
0010 Developed within opticks/ana, see::
0011
0012     ana/ellipse_intersect_circle.hh
0013     ana/ellipse_intersect_circle.cc
0014     ana/shape.py
0015     analytic/gdml.py
0016     ana/gplt.py
0017     ana/gpmt.py
0018
0019 Started from::
0020
0021     def points_inside_circle(points, center, radius):
0022         """
0023         :param points: (n,2) array of points
0024         :param center: (2,) coordinates of circle center
0025         :param radius:
0026         :return mask: boolean array of dimension (n,2) indicating if points are within the circle
0027         """
0028         return np.sqrt(np.sum(np.square(points-center),1)) - radius < 0.
0029
0030 **/
0031
0032 #include <stdio.h>
0033 #include <math.h>
0034 #include <assert.h>
0035
0036 //#define WITH_MASK 1
0037
0038 struct Point
0039 {
0040     double x ;
0041     double y ;
0042     int    i ;
0043     int    n ;
0044
0045     void dump(const char* msg) const
0046     {
0047         printf("Point %s : i %d n %d  (%10.4f, %10.4f)  \n", msg, i, n, x, y );
0048     }
0049
0050     void copy(const Point& other)
0051     {
0052         x = other.x ;
0053         y = other.y ;
0054         i = other.i ;
0055         n = other.n ;
0056     }
0057
0058     void set(double x_, double y_, int i_, int n_ )
0059     {
0060         x = x_ ;
0061         y = y_ ;
0062         i = i_ ;
0063         n = n_ ;
0064     }
0065
0066     static Point make(double px, double py, int i=-1, int n=-1)
0067     {
0068         Point p ;
0069         p.x = px;
0070         p.y = py;
0071         p.i = i ;
0072         p.n = n ;
0073         return p ;
0074     }
0075
0076     double distance(const Point& other) const
0077     {
0078         double dx = other.x - x ;
0079         double dy = other.y - y ;
0080         return sqrt(dx*dx + dy*dy) ;
0081     }
0082 };
0083
0084 struct Circle
0085 {
0086     Point center ;
0087     double radius ;
0088
0089     static Circle make(double cx, double cy, double r )
0090     {
0091         Circle c ;
0092         c.center = Point::make(cx, cy);
0093         c.radius = r ;
0094         return c ;
0095     }
0096
0097     bool is_point_inside( const Point& point)
0098     {
0099         double d = center.distance(point) ;
0100         bool inside = d - radius < 0. ;
0101         return inside ;
0102     }
0103 };
0104
0105
0106 struct Mask
0107 {
0108     bool* a ;
0109     int n ;
0110
0111     Mask(int n_)
0112        :
0113        a(new bool[n_]),
0114        n(n_)
0115     {
0116     }
0117
0118
0119     int first() const
0120     {
0121         for(int i=0 ; i < n ; i++) if(a[i]) return i ;
0122         return -1 ;
0123     }
0124
0125     int last() const
0126     {
0127         for(int i=0 ; i < n ; i++) if(a[n-1-i]) return n-1-i ;
0128         return -1 ;
0129     }
0130
0131     int count() const
0132     {
0133         int tot = 0 ;
0134         for(int i=0 ; i < n ; i++) if(a[i]) tot += 1 ;
0135         return tot ;
0136     }
0137
0138     void dump(const char* msg) const
0139     {
0140         printf("Mask %s n %d first %d last %d count %d \n", msg, n, first(), last(), count() );
0141     }
0142
0143
0144 };
0145
0146
0147 struct Ellipse
0148 {
0149     Point center ;
0150     double ax ;
0151     double ay ;
0152
0153     static Ellipse make( double px, double py, double ax_, double ay_ )
0154     {
0155          Ellipse e ;
0156          e.center = Point::make( px, py);
0157          e.ax = ax_ ;
0158          e.ay = ay_ ;
0159          return e ;
0160     }
0161
0162     void get_point( Point& p, int i, int n ) const
0163     {
0164          /**
0165          With typical ranges, do not get to fraction 1 : so does not reach 2pi and the initial point
0166
0167             range of i        :  0 -> n-1
0168             range of fraction :  0/n -> (n-1)/n  : 0 -> 1-1/n
0169
0170          **/
0171          assert( i > -1 && i < n );
0172
0173          double fraction = double(i)/double(n) ;
0174          double t = M_PI*2.*fraction ;
0175
0176          p.x = ax*cos(t) + center.x ;
0177          p.y = ay*sin(t) + center.y ;
0178          p.i = i ;
0179          p.n = n ;
0180     }
0181
0182     void dump_point( const char* msg, int i, int n ) const
0183     {
0184         Point point ;
0185         get_point(point, i, n );
0186         printf("%s  i %3d/%d   point %10.4f %10.4f \n", msg, i, n, point.x, point.y );
0187     }
0188
0189     void dump_points( const Mask& mask ) const
0190     {
0191          for(int i=0 ; i < mask.n ; i++)
0192          {
0193              if(!mask.a[i]) continue ;
0194              dump_point("", i, mask.n );
0195          }
0196     }
0197
0198 };
0199
0200
0201 struct Intersect
0202 {
0203     Point* p ;
0204     int mx ;
0205     int count ;
0206
0207     static Intersect make(int mx_)
0208     {
0209          Intersect isect ;
0210          isect.p = new Point[mx_] ;
0211          isect.mx = mx_ ;
0212          isect.count = 0 ;
0213          return isect ;
0214     }
0215
0216     void add(const Point& point)
0217     {
0218         assert( count < mx );
0219         p[count].copy(point) ;
0220         count += 1 ;
0221     }
0222
0223     void dump(const char* msg)
0224     {
0225         printf("Intersect %s : count %d crossings  \n", msg, count );
0226         for(int i=0 ; i < count ; i++)
0227         {
0228             p[i].dump("p");
0229         }
0230
0231     }
0232
0233 };
0234
0235
0236
0237 struct Ellipse_Intersect_Circle
0238 {
0239     Ellipse ellipse ;
0240     Circle circle ;
0241     int n ;
0242 #ifdef WITH_MASK
0243     Mask* mask ;
0244 #endif
0245     Intersect intersect ;
0246     bool verbose ;
0247
0248     enum {
0249        INSIDE,
0250        OUTSIDE
0251     };
0252
0253
0254     static Ellipse_Intersect_Circle make(double e_cx, double e_cy, double e_ax, double e_ay,  double c_cx, double c_cy, double c_r, int n, bool verbose )
0255     {
0256         Ellipse_Intersect_Circle ec ;
0257         ec.ellipse = Ellipse::make( e_cx, e_cy, e_ax, e_ay );
0258         ec.circle = Circle::make( c_cx, c_cy, c_r );
0259         ec.n = n ;
0260         ec.intersect = Intersect::make(4) ;
0261         ec.verbose = verbose ;
0262
0263         int ni = ec.find_intersects();
0264         if(verbose)
0265         {
0266             printf("Ellipse_Intersect_Circle found %d intersects\n", ni );
0267             ec.intersect.dump("EC") ;
0268         }
0269
0270
0271 #ifdef WITH_MASK
0272         ec.mask = new Mask(n) ;
0273         ec.find_mask();
0274 #endif
0275
0276         return ec ;
0277     }
0278
0279
0280     int find_intersects()
0281     {
0282          /*
0283          A direct comparison of states between the first
0284          and last is done in order to "close" the ellipse.
0285
0286          Crossings are checked for all points around the ellipse
0287          traversed in an anti-clockwise direction starting at right hand extremity
0288          and are collected in that order.
0289
0290          */
0291
0292          Point e ;
0293          ellipse.get_point(e, 0, n );
0294          int s0 = circle.is_point_inside(e) ? INSIDE : OUTSIDE ;
0295          int s_prev = s0 ;
0296
0297          for(int i=1 ; i < n ; i++)
0298          {
0299              ellipse.get_point(e, i, n );
0300              int s = circle.is_point_inside(e) ? INSIDE : OUTSIDE ;
0301              if( s != s_prev ) intersect.add(e) ;
0302              s_prev = s ;
0303          }
0304
0305          if( s_prev != s0 ) intersect.add(e) ;
0306
0307          return intersect.count ;
0308     }
0309
0310 #ifdef WITH_MASK
0311     int find_mask()
0312     {
0313          Point e ;
0314          int num_inside = 0 ;
0315          for(int i=0 ; i < mask->n ; i++)
0316          {
0317              ellipse.get_point(e, i, mask->n );
0318              bool inside = circle.is_point_inside(e);
0319              mask->a[i] = inside ;
0320              if(inside) num_inside += 1 ;
0321          }
0322          if(verbose) mask->dump("ellipse points inside circle");
0323          return num_inside ;
0324     }
0325 #endif
0326
0327 };
0328
0329