[gnuastro-commits] master 09ffaa7 1/5: Library(polygon): Added new functions inside polygon.c

[Mohammad Akhlaghi]

branch: master
commit 09ffaa7086d5446ac9731ab95521d1c20bc64f7c
Author: Sachin Kumar Singh <address@hidden>
Commit: Sachin Kumar Singh <address@hidden>

    Library(polygon): Added new functions inside polygon.c
    
    Two new functions are added in polygon.c. `gal_polygon_isconvex` is used
    to check whether the polygon is concave or convex.
    `gal_polygon_isinside` is used to check whether a given point
    lies inside the polygon or not. The function `gal_polygon_pin` is
    changed to `gal_polygon_isinside_convex`.
---
 bin/crop/onecrop.c     |   4 +-
 doc/gnuastro.texi      |  12 +++++-
 lib/gnuastro/polygon.h |   8 +++-
 lib/polygon.c          | 105 +++++++++++++++++++++++++++++++++++++++++++++++--
 4 files changed, 122 insertions(+), 7 deletions(-)

diff --git a/bin/crop/onecrop.c b/bin/crop/onecrop.c
index bfeb568..9f5f6a4 100644
--- a/bin/crop/onecrop.c
+++ b/bin/crop/onecrop.c
@@ -334,7 +334,7 @@ onecrop_ipolygon_fl(double *ipolygon, size_t nvertices, 
long *fpixel,
     for(i=0;i<size;++i)                                                 \
       {                                                                 \
         point[0]=i%s1+1; point[1]=i/s1+1;                               \
-        if(gal_polygon_pin(ipolygon, point, nvertices)==outpolygon)     \
+        if((*func_ptr)(ipolygon, point, nvertices)==outpolygon)           \
           ba[i]=*bb;                                                    \
       }                                                                 \
     free(bb);                                                           \
@@ -348,6 +348,7 @@ polygonmask(struct onecropparams *crp, void *array, long 
*fpixel_i,
   int type=crp->p->type;
   double *ipolygon, point[2];
   int outpolygon=crp->p->outpolygon;
+  int (*func_ptr)(double *, double *, size_t);
   size_t i, *ordinds, size=s0*s1, nvertices=crp->p->nvertices;
 
 
@@ -374,6 +375,7 @@ polygonmask(struct onecropparams *crp, void *array, long 
*fpixel_i,
       ipolygon[i*2+1] = crp->ipolygon[ordinds[i]*2+1] - fpixel_i[1];
     }
 
+  func_ptr=&gal_polygon_isinside;
 
   /* Go over all the pixels in the image and if they are within the
      polygon keep them if the user has asked for it.*/
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index e6414dd..f98afdc 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -23172,17 +23172,25 @@ that @code{v[0]} and @code{v[1]} are the positions of 
the first vertice to
 be considered.
 @end deftypefun
 
-@deftypefun int gal_polygon_pin (double @code{*v}, double @code{*p}, size_t 
@code{n})
+@deftypefun int gal_polygon_isinside_convex (double @code{*v}, double 
@code{*p}, size_t @code{n})
 Return @code{1} if the point @code{p} is within the polygon whose vertices
 are defined by @code{v} and @code{0} otherwise. Note that the vertices of
 the polygon have to be sorted in an anti-clock-wise manner.
 @end deftypefun
 
 @deftypefun int gal_polygon_ppropin (double @code{*v}, double @code{*p}, 
size_t @code{n})
-Similar to @code{gal_polygon_pin}, except that if the point @code{p} is on
+Similar to @code{gal_polygon_isinside_convex}, except that if the point 
@code{p} is on
 one of the edges of a polygon, this will return @code{0}.
 @end deftypefun
 
+@deftypefun int gal_polygon_isinside(double @code{*v}, double @code{*p}, 
size_t @code{n})
+Returns @code{0} if point @code{p} in inside a polygon, either convex or 
concave, whose vertices are defined by @code{v} and @code{0} otherwise.
+This function uses the 
@url{https://en.wikipedia.org/wiki/Point_in_polygon#Winding_number_algorithm, 
winding number algorithm}, to check the points.
+
+@deftypefun int gal_polygon_isconvex(double @code{*v}, size_t @code{n})
+Returns @code{1} if the polygon is convex with vertices defined by @code{v} 
and @code{0} if it is a concave polygon.
+Note that the vertices of the polygon should be sorted in an anti-clockwise 
manner.
+
 @deftypefun void gal_polygon_clip (double @code{*s}, size_t @code{n}, double 
@code{*c}, size_t @code{m}, double @code{*o}, size_t @code{*numcrn})
 Clip (find the overlap of) two polygons. This function uses the
 @url{https://en.wikipedia.org/wiki/Sutherland%E2%80%93Hodgman_algorithm,
diff --git a/lib/gnuastro/polygon.h b/lib/gnuastro/polygon.h
index d4f24d0..2838bb1 100644
--- a/lib/gnuastro/polygon.h
+++ b/lib/gnuastro/polygon.h
@@ -64,11 +64,17 @@ double
 gal_polygon_area(double *v, size_t n);
 
 int
-gal_polygon_pin(double *v, double *p, size_t n);
+gal_polygon_isinside_convex(double *v, double *p, size_t n);
 
 int
 gal_polygon_ppropin(double *v, double *p, size_t n);
 
+int
+gal_polygon_isinside(double *v, double *p, size_t n);
+
+int
+gal_polygon_isconvex(double *v, size_t n);
+
 void
 gal_polygon_clip(double *s, size_t n, double *c, size_t m,
                  double *o, size_t *numcrn);
diff --git a/lib/polygon.c b/lib/polygon.c
index 4c51223..33cf85f 100644
--- a/lib/polygon.c
+++ b/lib/polygon.c
@@ -225,7 +225,7 @@ gal_polygon_area(double *v, size_t n)
    traversed in order. See the comments above `gal_polygon_area' for an
    explanation about i and j and the loop.*/
 int
-gal_polygon_pin(double *v, double *p, size_t n)
+gal_polygon_isinside_convex(double *v, double *p, size_t n)
 {
   size_t i=0, j=n-1;
 
@@ -242,8 +242,8 @@ gal_polygon_pin(double *v, double *p, size_t n)
 
 
 
-/* Similar to gal_polygon_pin, except that if the point is on one of the
-   edges of a polygon, this will return 0. */
+/* Similar to gal_polygon_isinside_convex, except that if the point
+   is on one of the edges of a polygon, this will return 0. */
 int
 gal_polygon_ppropin(double *v, double *p, size_t n)
 {
@@ -268,6 +268,105 @@ gal_polygon_ppropin(double *v, double *p, size_t n)
 
 
 
+/* This fuction test if a point is inside the polygon using
+  winding number algorithm.
+
+  See its wiki here:
+  https://en.wikipedia.org/wiki/Point_in_polygon#Winding_number_algorithm
+
+  We have a polygon with `n' vertices whose vertices are in the array
+  `v' (with 2*n elements). Such that v[0], v[1] are the two
+  coordinates of the first vertice. The vertices also have to be
+  sorted in a counter clockwise fashion. We also have a point (with
+  coordinates (x, y) == (p[0], p[1]) and we want to see if it is inside the
+  polygon or not.
+
+  If point is outside the polygon, return 0.
+  If point is inside the polygon, return a non-zero number.
+  */
+int
+gal_polygon_isinside(double *v, double *p, size_t n)
+{
+  /* The winding number(wn) keeping the count of number of times the ray
+     crosses the polygon. */
+  size_t wn=0, i=0, j=n-1;
+
+  /* Loop through all the edges of the polygon*/
+  while(i<n)
+    {
+    /* edge from v[i] to v[i+1] in upward direction */
+    if(v[j*2+1] <= p[1]){
+      if(v[i*2+1] > p[1]){
+        /* p left of edge is an upward intersection, increase wn */
+        if( GAL_POLYGON_TRI_CROSS_PRODUCT(&v[j*2], &v[i*2], p) > 0 ){
+          wn++;
+        }
+      }
+    }
+    else{
+      /* edge from v[i] to v[i+1] in downward direction */
+      if(v[i*2+1] <= p[1]){
+        /* p right of edge is a downward intersection, decrease wn */
+        if( GAL_POLYGON_TRI_CROSS_PRODUCT(&v[j*2], &v[i*2], p) < 0 ){
+          wn--;
+        }
+      }
+    }
+    j=i++;
+  /* For a check:
+    printf("winding number: %ld, %.3f\n", wn);
+    */
+    }
+  return wn;
+
+}
+
+
+
+
+
+/* This function checks if the polygon is convex or concave by
+   testing all 3 consecutive points of the sorted polygon.
+   If any of the test returns false, the the polygon is concave
+   else it is convex.
+   return 1: convex polygon
+   return 0: concave polygon
+   */
+int
+gal_polygon_isconvex(double *v, size_t n)
+{
+  size_t i;
+  int flag=1;
+
+  /* check the first n-1 edges made by n points. */
+  for(i=0; i<n-2; i++)
+    {
+    if( GAL_POLYGON_LEFT_OF_LINE(&v[i*2], &v[(i+1)*2], &v[(i+2)*2]) )
+      {
+        continue;
+      }
+    else
+      {
+        return 0;
+      }
+    }
+
+  /* check the edge between nth and 1st point */
+  if(flag)
+    {
+      if( GAL_POLYGON_LEFT_OF_LINE(&v[(n-2)*2], &v[(n-1)*2], &v[0]) )
+        return 1;
+      else
+        return 0;
+  }
+
+  return 1;
+}
+
+
+
+
+
 /* Find the intersection of a line segment (Aa--Ab) and an infinite
    line (Ba--Bb) and put the intersection point in the output `o'. All
    the points are assumed to be two element double arrays already