[gnuastro-commits] master a22b8d0 3/5: Library (polygon): improved documentation, new functions in NEWS

[Mohammad Akhlaghi]

branch: master
commit a22b8d021c86c6cceb2e07666eeb88d1165ce8bd
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>

    Library (polygon): improved documentation, new functions in NEWS
    
    With the new set of polygon functions defined in the previous commit, it
    was necessary to improve the relevant part of the documentation so they
    blend in better with the others.
    
    For example the discussion on the definitions/differences of a concave and
    convex polygon is now brough to the start of the `polygon.h' subsection and
    the new functions were placed in a better logical order
    (`gal_polygon_isinside' is now placed before
    `gal_polygon_isinside_convex'). This re-ordering was also reflected in
    `polygon.c' and `polygon.h'.
    
    Also, Sachin's name was added as a contributing author in the two polygon
    library files.
---
 NEWS                   |   7 ++-
 doc/gnuastro.texi      | 137 +++++++++++++++++++++-----------------------
 lib/gnuastro/polygon.h |  13 +++--
 lib/polygon.c          | 151 +++++++++++++++++++++++++------------------------
 4 files changed, 152 insertions(+), 156 deletions(-)

diff --git a/NEWS b/NEWS
index b6e84ac..5d97a7a 100644
--- a/NEWS
+++ b/NEWS
@@ -40,8 +40,10 @@ See the end of the file for license conditions.
 
   Library:
    - GAL_SPECLINES_INVALID_MAX: Total number of spectral lines, plus 1.
-   - gal_txt_trim_space: trim white space before and after a string.
    - GAL_ARITHMETIC_OP_QUANTILE: operator for `gal_arithmetic'.
+   - gal_txt_trim_space: trim white space before and after a string.
+   - gal_polygon_isconvex: identify if a polygon is convex or concave.
+   - gal_polygon_isinside: if point is inside polygon (convex or concave).
 
 ** Removed features
 
@@ -61,6 +63,9 @@ See the end of the file for license conditions.
      value of zero. As a result, the Sky and Sky standard deviation
      extensions will be measured over all the tiles.
 
+  Library:
+   - gal_polygon_isinside_convex: new name for `gal_polygon_pin'.
+
 ** Bugs fixed
   bug #57300: MakeCatalog memory crash when input dataset has units.
   bug #57301: MakeCatalog using river sum instead of mean times by clump area.
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 64844a5..3bd5d4f 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -23086,16 +23086,29 @@ pixels of the overlap will be put in @code{fpixel_o} 
and @code{lpixel_o}.
 @node Polygons, Qsort functions, Bounding box, Gnuastro library
 @subsection Polygons (@file{polygon.h})
 
-Polygons are commonly necessary in image processing. In Gnuastro, they are
-used in Crop (see @ref{Crop}) for cutting out non-rectangular regions of a
-image. Warp (see @ref{Warp}) uses them to warp the images into a new pixel
-grid. The polygon related Gnuastro library macros and functions are
-introduced here.
+Polygons are commonly necessary in image processing.
+For example in Crop they are used for cutting out non-rectangular regions of a 
image (see @ref{Crop}), and in Warp, for mapping different pixel grids over 
each other (see @ref{Warp}).
 
-In all the functions here the vertices (and points) are defined as an
-array. So a polygon with 4 vertices will be identified with an array of 8
-elements with the first two elements keeping the 2D coordinates of the
-first vertice and so on.
+@cindex Convex polygons
+@cindex Concave polygons
+@cindex Polygons, Convex
+@cindex Polygons, Concave
+Polygons come in two classes: convex and concave (or generally, non-convex!), 
see below for a demonstration.
+Convex polygons are those where all inner angles are less than 180 degrees.
+By contrast, a convex polygon is one where an inner angle may be more than 180 
degress.
+
+@example
+            Concave Polygon        Convex Polygon
+
+             D --------C          D------------- C
+              \        |        E /              |
+               \E      |          \              |
+               /       |           \             |
+              A--------B             A ----------B
+@end example
+
+In all the functions here the vertices (and points) are defined as an array.
+So a polygon with 4 vertices will be identified with an array of 8 elements 
with the first two elements keeping the 2D coordinates of the first vertice and 
so on.
 
 @deffn Macro GAL_POLYGON_MAX_CORNERS
 The largest number of vertices a polygon can have in this library.
@@ -23103,32 +23116,26 @@ The largest number of vertices a polygon can have in 
this library.
 
 @deffn Macro GAL_POLYGON_ROUND_ERR
 @cindex Round-off error
-We have to consider floating point round-off errors when dealing with
-polygons. For example we will take @code{A} as the maximum of @code{A} and
-@code{B} when @code{A>B-GAL_POLYGON_ROUND_ERR}.
+We have to consider floating point round-off errors when dealing with polygons.
+For example we will take @code{A} as the maximum of @code{A} and @code{B} when 
@code{A>B-GAL_POLYGON_ROUND_ERR}.
 @end deffn
 
 @deftypefun void gal_polygon_ordered_corners (double @code{*in}, size_t 
@code{n}, size_t @code{*ordinds})
-We have a simple polygon (that can result from projection, so its edges
-don't collide or it doesn't have holes) and we want to order its corners in
-an anticlockwise fashion. This is necessary for clipping it and finding its
-area later. The input vertices can have practically any order.
-
-The input (@code{in}) is an array containing the coordinates (two values)
-of each vertice. @code{n} is the number of corners. So @code{in} should
-have @code{2*n} elements. The output (@code{ordinds}) is an array with
-@code{n} elements specifying the indexs in order. This array must have been
-allocated before calling this function. The indexes are output for more
-generic usage, for example in a homographic transform (necessary in warping
-an image, see @ref{Warping basics}), the necessary order of vertices is the
-same for all the pixels. In other words, only the positions of the vertices
-change, not the way they need to be ordered. Therefore, this function would
-only be necessary once.
-
-As a summary, the input is unchanged, only @code{n} values will be put in
-the @code{ordinds} array. Such that calling the input coordinates in the
-following fashion will give an anti-clockwise order when there are 4
-vertices:
+We have a simple polygon (that can result from projection, so its edges don't 
collide or it doesn't have holes) and we want to order its corners in an 
anticlockwise fashion.
+This is necessary for clipping it and finding its area later.
+The input vertices can have practically any order.
+
+The input (@code{in}) is an array containing the coordinates (two values) of 
each vertice.
+@code{n} is the number of corners.
+So @code{in} should have @code{2*n} elements.
+The output (@code{ordinds}) is an array with @code{n} elements specifying the 
indexs in order.
+This array must have been allocated before calling this function.
+The indexes are output for more generic usage, for example in a homographic 
transform (necessary in warping an image, see @ref{Warping basics}), the 
necessary order of vertices is the same for all the pixels.
+In other words, only the positions of the vertices change, not the way they 
need to be ordered.
+Therefore, this function would only be necessary once.
+
+As a summary, the input is unchanged, only @code{n} values will be put in the 
@code{ordinds} array.
+Such that calling the input coordinates in the following fashion will give an 
anti-clockwise order when there are 4 vertices:
 
 @example
 1st vertice: in[ordinds[0]*2], in[ordinds[0]*2+1]
@@ -23139,42 +23146,35 @@ vertices:
 
 @cindex Convex Hull
 @noindent
-The implementation of this is very similar to the Graham scan in finding
-the Convex Hull. However, in projection we will never have a concave
-polygon (the left condition below, where this algorithm will get to E
-before D), we will always have a convex polygon (right case) or E won't
-exist!
-
-@example
-            Concave Polygon        Convex Polygon
-
-             D --------C          D------------- C
-              \        |        E /              |
-               \E      |          \              |
-               /       |           \             |
-              A--------B             A ----------B
-@end example
+The implementation of this is very similar to the Graham scan in finding the 
Convex Hull.
+However, in projection we will never have a concave polygon (the left 
condition below, where this algorithm will get to E before D), we will always 
have a convex polygon (right case) or E won't exist!
+This is because we are always going to be calculating the area of the overlap 
between a quadrilateral and the pixel grid or the quadrilateral its self.
 
-This is because we are always going to be calculating the area of the
-overlap between a quadrilateral and the pixel grid or the quadrilateral
-its self.
+The @code{GAL_POLYGON_MAX_CORNERS} macro is defined so there will be no need 
to allocate these temporary arrays separately.
+Since we are dealing with pixels, the polygon can't really have too many 
vertices.
+@end deftypefun
 
-The @code{GAL_POLYGON_MAX_CORNERS} macro is defined so there will be no
-need to allocate these temporary arrays separately. Since we are dealing
-with pixels, the polygon can't really have too many vertices.
+@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.
 @end deftypefun
 
 @deftypefun double gal_polygon_area (double @code{*v}, size_t @code{n})
-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.
+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})
+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}.
 @end deftypefun
 
 @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.
+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}.
+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})
@@ -23182,21 +23182,10 @@ Similar to @code{gal_polygon_isinside_convex}, except 
that if the point @code{p}
 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.
-@end deftypefun
-
-@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.
-@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. Note that the vertices of
-both polygons have to be sorted in an anti-clock-wise manner.
+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.
+Note that the vertices of both polygons have to be sorted in an 
anti-clock-wise manner.
 @end deftypefun
 
 
diff --git a/lib/gnuastro/polygon.h b/lib/gnuastro/polygon.h
index 2838bb1..8d4c1e2 100644
--- a/lib/gnuastro/polygon.h
+++ b/lib/gnuastro/polygon.h
@@ -5,6 +5,7 @@ This is part of GNU Astronomy Utilities (Gnuastro) package.
 Original author:
      Mohammad Akhlaghi <address@hidden>
 Contributing author(s):
+     Sachin Kumar Singh <address@hidden>
 Copyright (C) 2015-2020, Free Software Foundation, Inc.
 
 Gnuastro is free software: you can redistribute it and/or modify it
@@ -60,20 +61,20 @@ __BEGIN_C_DECLS  /* From C++ preparations */
 void
 gal_polygon_ordered_corners(double *in, size_t n, size_t *ordinds);
 
+int
+gal_polygon_isconvex(double *v, size_t n);
+
 double
 gal_polygon_area(double *v, size_t n);
 
 int
-gal_polygon_isinside_convex(double *v, double *p, size_t n);
-
-int
-gal_polygon_ppropin(double *v, double *p, size_t n);
+gal_polygon_isinside(double *v, double *p, size_t n);
 
 int
-gal_polygon_isinside(double *v, double *p, size_t n);
+gal_polygon_isinside_convex(double *v, double *p, size_t n);
 
 int
-gal_polygon_isconvex(double *v, size_t n);
+gal_polygon_ppropin(double *v, double *p, size_t n);
 
 void
 gal_polygon_clip(double *s, size_t n, double *c, size_t m,
diff --git a/lib/polygon.c b/lib/polygon.c
index 2d794ff..dec8698 100644
--- a/lib/polygon.c
+++ b/lib/polygon.c
@@ -5,6 +5,7 @@ This is part of GNU Astronomy Utilities (Gnuastro) package.
 Original author:
      Mohammad Akhlaghi <address@hidden>
 Contributing author(s):
+     Sachin Kumar Singh <address@hidden>
 Copyright (C) 2015-2020, Free Software Foundation, Inc.
 
 Gnuastro is free software: you can redistribute it and/or modify it
@@ -183,6 +184,44 @@ gal_polygon_ordered_corners(double *in, size_t n, size_t 
*ordinds)
 
 
 
+/* 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;
+}
+
+
+
+
+
 /* The area of a polygon is the sum of the vector products of all the
    vertices in a counterclockwise order. See the Wikipedia page for
    Polygon for more information.
@@ -213,61 +252,6 @@ gal_polygon_area(double *v, size_t n)
 
 
 
-/* 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 p[0], p[1]) and we want to see if it is inside the
-   polygon or not.
-
-   If the point is inside the polygon, it will always be to the left
-   of the edge connecting the two vertices when the vertices are
-   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)
-{
-  size_t i=0, j=n-1;
-
-  while(i<n)
-    {
-      if( GAL_POLYGON_LEFT_OF_LINE(&v[j*2], &v[i*2], p) )
-        j=i++;
-      else return 0;
-    }
-  return 1;
-}
-
-
-
-
-
-/* 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)
-{
-  size_t i=0, j=n-1;
-
-  while(i<n)
-    {
-      /* For a check:
-      printf("(%.2f, %.2f), (%.2f, %.2f), (%.2f, %.2f)=%f\n",
-             v[j*2], v[j*2+1], v[i*2], v[i*2+1], p[0], p[1],
-             GAL_POLYGON_TRI_CROSS_PRODUCT(&v[j*2], &v[i*2], p));
-      */
-      if( GAL_POLYGON_PROP_LEFT_OF_LINE(&v[j*2], &v[i*2], p) )
-        j=i++;
-      else
-        return 0;
-    }
-  return 1;
-}
-
-
-
-
-
 /* 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
@@ -321,37 +305,54 @@ gal_polygon_isinside(double *v, double *p, size_t n)
 
 
 
-/* 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.
+/* 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 p[0], p[1]) and we want to see if it is inside the
+   polygon or not.
 
-   return 1: convex polygon
-   return 0: concave polygon
-   */
+   If the point is inside the polygon, it will always be to the left
+   of the edge connecting the two vertices when the vertices are
+   traversed in order. See the comments above `gal_polygon_area' for an
+   explanation about i and j and the loop.*/
 int
-gal_polygon_isconvex(double *v, size_t n)
+gal_polygon_isinside_convex(double *v, double *p, size_t n)
 {
-  size_t i;
-  int flag=1;
+  size_t i=0, j=n-1;
 
-  /* Check the first n-1 edges made by n points. */
-  for(i=0; i<n-2; i++)
+  while(i<n)
     {
-      if( GAL_POLYGON_LEFT_OF_LINE(&v[i*2], &v[(i+1)*2], &v[(i+2)*2]) )
-        continue;
-      else
-        return 0;
+      if( GAL_POLYGON_LEFT_OF_LINE(&v[j*2], &v[i*2], p) )
+        j=i++;
+      else return 0;
     }
+  return 1;
+}
 
-  /* Check the edge between nth and 1st point */
-  if(flag)
+
+
+
+
+/* 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)
+{
+  size_t i=0, j=n-1;
+
+  while(i<n)
     {
-      if( GAL_POLYGON_LEFT_OF_LINE(&v[(n-2)*2], &v[(n-1)*2], &v[0]) )
-        return 1;
+      /* For a check:
+      printf("(%.2f, %.2f), (%.2f, %.2f), (%.2f, %.2f)=%f\n",
+             v[j*2], v[j*2+1], v[i*2], v[i*2+1], p[0], p[1],
+             GAL_POLYGON_TRI_CROSS_PRODUCT(&v[j*2], &v[i*2], p));
+      */
+      if( GAL_POLYGON_PROP_LEFT_OF_LINE(&v[j*2], &v[i*2], p) )
+        j=i++;
       else
         return 0;
     }
-
   return 1;
 }