[gnuastro-commits] master 698df79 1/2: Library(polygon): Functions to check and convert to counter-clokwise direction

[Mohammad Akhlaghi]

branch: master
commit 698df79ab7c2a2c428d74d1abc930be6122793ac
Author: Sachin Kumar Singh <address@hidden>
Commit: Sachin Kumar Singh <address@hidden>

    Library(polygon): Functions to check and convert to counter-clokwise 
direction
    
    New functions, `gal_polygon_is_counterclockwise` and
    `gal_polygon_to_counterclockwise`, are added to check if the given vertices
    are present in counter-clockwise direction and if not then arrange these
    points to counter-clockwise direction.
    
    Names of few functions were stylystic changed to be more intutive and
    readable. Documentations were added for all these changes.
---
 bin/crop/onecrop.c     |   6 +--
 doc/gnuastro.texi      |  24 ++++++++---
 lib/gnuastro/polygon.h |  12 ++++--
 lib/polygon.c          | 113 +++++++++++++++++++++++++++++++++++++++++++++++--
 4 files changed, 139 insertions(+), 16 deletions(-)

diff --git a/bin/crop/onecrop.c b/bin/crop/onecrop.c
index 31c5855..c984f26 100644
--- a/bin/crop/onecrop.c
+++ b/bin/crop/onecrop.c
@@ -385,9 +385,9 @@ polygonmask(struct onecropparams *crp, void *array, long 
*fpixel_i,
      concave polygons the process is more complex and thus
      slower. Therefore when the polygon is convex, its better to use the
      simpler/faster function. */
-  isinside = ( gal_polygon_isconvex(ipolygon, crp->p->nvertices)
-               ? gal_polygon_isinside_convex
-               : gal_polygon_isinside );
+  isinside = ( gal_polygon_is_convex(ipolygon, crp->p->nvertices)
+               ? gal_polygon_is_inside_convex
+               : gal_polygon_is_inside );
 
   /* 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 db027f4..1bf8882 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -23192,7 +23192,7 @@ The @code{GAL_POLYGON_MAX_CORNERS} macro is defined so 
there will be no need to
 Since we are dealing with pixels, the polygon can't really have too many 
vertices.
 @end deftypefun
 
-@deftypefun int gal_polygon_isconvex(double @code{*v}, size_t @code{n})
+@deftypefun int gal_polygon_is_convex(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.
 @end deftypefun
@@ -23202,24 +23202,36 @@ Find the area of a polygon with vertices defined in 
@code{v}.
 @code{v} points to an array of doubles which keep the positions of the 
vertices such that @code{v[0]} and @code{v[1]} are the positions of the first 
vertice to be considered.
 @end deftypefun
 
-@deftypefun int gal_polygon_isinside(double @code{*v}, double @code{*p}, 
size_t @code{n})
+@deftypefun int gal_polygon_is_inside(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.
 The vertices of the polygon are defined by @code{v} and @code{0} otherwise, 
they have to be ordered in an anti-clockwise manner.
 This function uses the 
@url{https://en.wikipedia.org/wiki/Point_in_polygon#Winding_number_algorithm, 
winding number algorithm}, to check the points.
-Note that this is a generic function (working on both concave and convex 
polygons, so if you know before-hand that your polygon is convex, it is much 
more efficient to use @code{gal_polygon_isinside_convex}.
+Note that this is a generic function (working on both concave and convex 
polygons, so if you know before-hand that your polygon is convex, it is much 
more efficient to use @code{gal_polygon_is_inside_convex}.
 @end deftypefun
 
-@deftypefun int gal_polygon_isinside_convex (double @code{*v}, double 
@code{*p}, size_t @code{n})
+@deftypefun int gal_polygon_is_inside_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}.
-The polygon is assumed to be convex, for a more generic function that deals 
with concave and convex polygons, see @code{gal_polygon_isinside}.
+The polygon is assumed to be convex, for a more generic function that deals 
with concave and convex polygons, see @code{gal_polygon_is_inside}.
 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_isinside_convex}, except that if the point 
@code{p} is on
+Similar to @code{gal_polygon_is_inside_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_is_counterclockwise(double @code{*v}, size_t 
@code{n})
+Returns  @code{1} if polygon is sorted in anti-clockwise and @code{0} 
otherwise.
+It uses winding which defines the relative order in which the vertices of a 
polygon are listed to determine the orientation of vertices.
+If orientation is positive vertices are in clockwise direction, else is 
negative for anti-clockwise direction. Zero sum implies a figure like @code{8}, 
with equal orientation in both direction.
+@end deftypefun
+
+@deftypefun int gal_polygon_to_counterclockwise(double @code{*v}, size_t 
@code{n})
+Arrange the vertices of polygon in counter-clockwise direction if input is in 
clockwise direction, otherwise do nothing.
+It uses @code{gal_polygon_is_counterclockwise} to determine the orientation of 
polygons.
+Return @code{1} on successful execution, @code{0} otherwise.
+@end deftypefun
+
 @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, 
Sutherland-Hodgman} polygon clipping algorithm.
diff --git a/lib/gnuastro/polygon.h b/lib/gnuastro/polygon.h
index 8d4c1e2..e45b4be 100644
--- a/lib/gnuastro/polygon.h
+++ b/lib/gnuastro/polygon.h
@@ -62,20 +62,26 @@ void
 gal_polygon_ordered_corners(double *in, size_t n, size_t *ordinds);
 
 int
-gal_polygon_isconvex(double *v, size_t n);
+gal_polygon_is_convex(double *v, size_t n);
 
 double
 gal_polygon_area(double *v, size_t n);
 
 int
-gal_polygon_isinside(double *v, double *p, size_t n);
+gal_polygon_is_inside(double *v, double *p, size_t n);
 
 int
-gal_polygon_isinside_convex(double *v, double *p, size_t n);
+gal_polygon_is_inside_convex(double *v, double *p, size_t n);
 
 int
 gal_polygon_ppropin(double *v, double *p, size_t n);
 
+int
+gal_polygon_is_counterclockwise(double *v, size_t n);
+
+int
+gal_polygon_to_counterclockwise(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 dec8698..c0fe4b1 100644
--- a/lib/polygon.c
+++ b/lib/polygon.c
@@ -192,7 +192,7 @@ gal_polygon_ordered_corners(double *in, size_t n, size_t 
*ordinds)
    return 0: concave polygon
    */
 int
-gal_polygon_isconvex(double *v, size_t n)
+gal_polygon_is_convex(double *v, size_t n)
 {
   size_t i;
   int flag=1;
@@ -266,7 +266,7 @@ gal_polygon_area(double *v, size_t n)
   If point is inside the polygon, return a non-zero number.
   */
 int
-gal_polygon_isinside(double *v, double *p, size_t n)
+gal_polygon_is_inside(double *v, double *p, size_t n)
 {
   /* The winding number(wn) keeping the count of number of times the ray
      crosses the polygon. */
@@ -317,7 +317,7 @@ gal_polygon_isinside(double *v, double *p, 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_isinside_convex(double *v, double *p, size_t n)
+gal_polygon_is_inside_convex(double *v, double *p, size_t n)
 {
   size_t i=0, j=n-1;
 
@@ -334,7 +334,7 @@ gal_polygon_isinside_convex(double *v, double *p, size_t n)
 
 
 
-/* Similar to gal_polygon_isinside_convex, except that if the point
+/* Similar to gal_polygon_is_inside_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)
@@ -360,6 +360,111 @@ gal_polygon_ppropin(double *v, double *p, size_t n)
 
 
 
+/* This function uses the concept of winding, which defines the
+   relative order in which the vertices of a polygon are
+   listed to determine the orientation of vertices. If orientation
+   is positive vertices are in clockwise direction, else is negative
+   for counter-clockwise direction. Zero sum implies a figure like 8, with
+   equal orientation in both direction.
+
+   See the link below for a detailed description:
+   "https://www.element84.com/blog/determining-the-winding-of-a-
+    polygon-given-as-a-set-of-ordered-points"
+
+   return 1: sorted in counter-clockwise order.
+   return 0: sorted clockwise order or weight of orientation is equal.
+   */
+int
+gal_polygon_is_counterclockwise(double *v, size_t n)
+{
+  double sum=0.0f;
+  size_t i=0, j=n-1;
+
+  while(i<n)
+    {
+      sum+=( (v[i*2]-v[j*2])*(v[i*2+1]+v[j*2+1]) );
+      j=i++;
+    }
+
+  if(sum>0) return 0;
+  else      return 1;
+
+}
+
+
+
+
+
+/* This function checks if the vertices are actually sorted in
+   the counterclockwise. If they are do nothing, otherwise if
+   they are clockwise, convert them to counter-clockwise direction
+   return 1: success
+   return 0: error
+   */
+int
+gal_polygon_to_counterclockwise(double *v, size_t n)
+{
+  size_t i, j=0;
+  size_t *permutation;
+  gal_data_t *temp=NULL;
+  int orientation=gal_polygon_is_counterclockwise(v, n);
+
+  switch(orientation)
+    {
+      /* Polygon is already sorted in counterclockwise order, do nothing */
+      case 1:
+        return 1;
+
+      /* Polygon is sorted in clockwise order, convert to counterclockwise
+        direction. */
+      case 0:
+        /* Allocate space for permutation array, which stores the order of
+          index in which the vertices are to be ordered for a
+          counter-clockwise direction. */
+        permutation=gal_pointer_allocate(GAL_TYPE_SIZE_T, 2*n, 0,
+                                          __func__, "permutation");
+        for(i=0; i<=2*n-1;)
+          {
+            /* Below sequence of steps ensures that the permutation has
+               indexes reversed but in order (x,y). Simple reversing would
+               have turned it in (y,x) format. */
+            j++;
+            permutation[j]  =(2*n-1-i++);
+            permutation[j-1]=(2*n-1-i++);
+            j++;
+          }
+
+        /* Put the vertices in the `gal_data_t' object */
+        temp=gal_data_alloc(v, GAL_TYPE_FLOAT64, 1, &n, NULL, 0,
+                          -1, 0, NULL, NULL, NULL);
+
+        /* Apply permutations to just reverse the order of clockwise to
+          counter-clockwise. */
+        gal_permutation_apply(temp, permutation);
+
+        temp->array=NULL;
+
+        /* Free allocated spaces. */
+        gal_data_free(temp);
+        free(permutation);
+        return 1;
+
+
+      default:
+        error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s"
+          "polygon orientation code %d not recognized",
+          __func__, PACKAGE_BUGREPORT, orientation);
+        return 0;
+    }
+
+  return 0;
+
+}
+
+
+
+
+
 /* 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