[gnuastro-commits] master 7323d4b 1/2: polygon.h and qsort.h documented in the book

[Mohammad Akhlaghi]

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

    polygon.h and qsort.h documented in the book
    
    The functions, macros and variables of these headers has been documented in
    the book.
---
 doc/gnuastro.texi      |  180 +++++++++++++++++++++++++++++++++++++++++++++++-
 lib/gnuastro/polygon.h |   75 --------------------
 lib/polygon.c          |  124 +++++++++++++++++++++------------
 3 files changed, 259 insertions(+), 120 deletions(-)

diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 6d98d2f..66648e0 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -514,6 +514,8 @@ Gnuastro library
 * FITS files::                  Working with FITS datat.
 * Linked lists::                Various types of linked lists.
 * Mesh grid for an image::      Breaking an image into a grid.
+* Polygons::                    Working with the vertices of a polygon.
+* Qsort functions::             Helper functions for Qsort.
 
 FITS files (@file{fits.h})
 
@@ -14763,6 +14765,8 @@ problems. It will stabilize with the removal of this 
notice. Check the
 * FITS files::                  Working with FITS datat.
 * Linked lists::                Various types of linked lists.
 * Mesh grid for an image::      Breaking an image into a grid.
+* Polygons::                    Working with the vertices of a polygon.
+* Qsort functions::             Helper functions for Qsort.
 @end menu
 
 @node Overall package, Array manipulation, Gnuastro library, Gnuastro library
@@ -15887,7 +15891,7 @@ within the nodes will also be discarded.
 @end deftypefun
 
 
address@hidden Mesh grid for an image,  , Linked lists, Gnuastro library
address@hidden Mesh grid for an image, Polygons, Linked lists, Gnuastro library
 @subsection Mesh grid for an image (@file{mesh.h})
 
 The basic concepts behind the mesh and channel grids in Gnuastro were
@@ -16103,16 +16107,190 @@ changed by this function.
 
 
 
address@hidden Polygons, Qsort functions, Mesh grid for an image, Gnuastro 
library
address@hidden Polygons (@file{polygon.h})
 
+Polygons are commonly necessary in image processing. In Gnuastro, they are
+used in ImageCrop (see @ref{ImageCrop}) for cutting out non-rectangular
+regions of a image. ImageWarp (see @ref{ImageWarp}) uses them to warp the
+images into a new pixel grid. The polygon related Gnuastro library macros
+and functions are introduced here.
 
+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.
 
address@hidden Macro GAL_POLYGON_MAX_CORNERS
+The largest number of vertices a polygon can have in this library.
address@hidden deffn
+
address@hidden Macro GAL_POLYGON_ROUND_ERR
address@hidden 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
address@hidden when @code{A>B-GAL_POLYGON_ROUND_ERR}.
address@hidden deffn
+
address@hidden 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
address@hidden 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:
+
address@hidden
+1st vertice: in[ordinds[0]*2], in[ordinds[0]*2+1]
+2nd vertice: in[ordinds[1]*2], in[ordinds[1]*2+1]
+3rd vertice: in[ordinds[2]*2], in[ordinds[2]*2+1]
+4th vertice: in[ordinds[3]*2], in[ordinds[3]*2+1]
address@hidden example
+
address@hidden Convex Hull
address@hidden
+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!
+
address@hidden
+            Concave Polygon        Convex Polygon
+
+             D --------C          D------------- C
+              \        |        E /              |
+               \E      |          \              |
+               /       |           \             |
+              A--------B             A ----------B
address@hidden example
+
+This is because we are always going to be calculating the area of the
+overlap between a quadrilateral and the pixel grid or the quadrilaterral
+its self.
+
+The @code{GAL_POLYGON_MAX_CORNERS} macro is defined so there will be no
+need to allocate these temporary arrays seprately. Since we are dealing
+with pixels, the polygon can't really have too many vertices.
+
address@hidden deftypefun
+
address@hidden 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.
address@hidden deftypefun
+
address@hidden int gal_polygon_pin (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.
address@hidden deftypefun
+
address@hidden 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
+one of the edges of a polygon, this will return @code{0}.
address@hidden deftypefun
+
address@hidden 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
address@hidden://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.
address@hidden deftypefun
 
 
 
 
address@hidden Qsort functions,  , Polygons, Gnuastro library
address@hidden Qsort functions (@file{qsort.h})
 
+The C programming language comes with the @code{qsort} (Quick sort)
+function. By allowing you to define the sorting based on the data,
address@hidden is a generic function which allows you to sort any kind of
+data structure. The functions here are all meant to be passed onto
address@hidden for sorting the respective types of data in the respective
+manner. Because of this all these functions have the same two argument
+types.
 
address@hidden {Global variable} {float *gal_qsort_index_arr}
+Pointer to a floating point array to use as a reference in
address@hidden, see the explanation there for
+more.
address@hidden deffn
 
address@hidden int gal_qsort_int_decreasing (const void @code{*a}, const void 
@code{*b})
+When passed to @code{qsort}, this function will sort an integer array in
+decreasing order.
address@hidden deftypefun
+
address@hidden int gal_qsort_int_increasing (const void @code{*a}, const void 
@code{*b})
+When passed to @code{qsort}, this function will sort an integer array in
+increasing order.
address@hidden deftypefun
+
address@hidden int gal_qsort_float_decreasing (const void @code{*a}, const void 
@code{*b})
+When passed to @code{qsort}, this function will sort a single precision
+floating point array in decreasing order.
address@hidden deftypefun
+
address@hidden int gal_qsort_float_increasing (const void @code{*a}, const void 
@code{*b})
+When passed to @code{qsort}, this function will sort a single precision
+floating point array in increasing order.
address@hidden deftypefun
+
address@hidden int gal_qsort_double_decreasing (const void @code{*a}, const 
void @code{*b})
+When passed to @code{qsort}, this function will sort a double precision
+floating point array in decreasing order.
address@hidden deftypefun
+
address@hidden int gal_qsort_double_increasing (const void @code{*a}, const 
void @code{*b})
+When passed to @code{qsort}, this function will sort a double precision
+floating point array in increasing order.
address@hidden deftypefun
+
address@hidden int gal_qsort_index_float_decreasing (const void @code{*a}, 
const void @code{*b})
+When passed to @code{qsort}, this function will sort a @code{size_t} array
+based on decreasing values in the @code{gal_qsort_index_arr} single
+precision floating poiny array. The floating point array will not be
+changed, it is only read. For example, if we have the following source
+code:
+
address@hidden
+#include <stdio.h>
+#include <stdlib.h>           /* qsort is defined in stdlib.h. */
+#include <gnuastro/qsort.h>
+
+int
+main (void)
address@hidden
+  size_t address@hidden, 1, 2, address@hidden;
+  float address@hidden,0.2,1.8,address@hidden;
+  gal_qsort_index_arr=f;
+  qsort(s, 4, sizeof(size_t), gal_qsort_index_float_decreasing);
+  printf("%lu, %lu, %lu, %lu\n", s[0], s[1], s[2], s[3]);
+  return EXIT_SUCCESS;
address@hidden
address@hidden example
+
address@hidden
+The output will be: @code{2, 0, 1, 3}.
address@hidden deftypefun
 
 @node Developing, GNU Astronomy Utilities list, Libraries, Top
 @chapter Developing
diff --git a/lib/gnuastro/polygon.h b/lib/gnuastro/polygon.h
index 85375e0..2325eee 100644
--- a/lib/gnuastro/polygon.h
+++ b/lib/gnuastro/polygon.h
@@ -74,81 +74,6 @@ gal_polygon_clip(double *s, size_t n, double *c, size_t m,
                  double *o, size_t *numcrn);
 
 
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-/***************************************************************/
-/**************            MACROS             ******************/
-/***************************************************************/
-/* The cross product of two points from the center. */
-#define GAL_POLYGON_CROSS_PRODUCT(A, B) ( (A)[0]*(B)[1] - (B)[0]*(A)[1] )
-
-
-
-
-/* Find the cross product (2*area) between three points. Each point is
-   assumed to be a pointer that has atleast two values within it. */
-#define GAL_POLYGON_TRI_CROSS_PRODUCT(A, B, C)    \
-  ( ( (B)[0]-(A)[0] ) * ( (C)[1]-(A)[1] ) -       \
-    ( (C)[0]-(A)[0] ) * ( (B)[1]-(A)[1] ) )       \
-
-
-
-
-
-/* We have the line A-B. We want to see if C is to the left of this
-   line or to its right. This function will return 1 if it is to the
-   left. It uses the basic property of vector multiplication: If the
-   three points are anti-clockwise (the point is to the left), then
-   the vector multiplication is positive, if it is negative, then it
-   is clockwise (c is to the right).
-
-   Ofcourse it is very important that A be below or equal to B in both
-   the X and Y directions. The rounding error might give
-   -0.0000000000001 (I didn't count the number of zeros!!) instead of
-   zero for the area. Zero would indicate that they are on the same
-   line in this case this should give a true result.
-*/
-#define GAL_POLYGON_LEFT_OF_LINE(A, B, C)                               \
-  ( GAL_POLYGON_TRI_CROSS_PRODUCT((A), (B), (C)) > -GAL_POLYGON_ROUND_ERR ) /* 
>= 0 */
-
-
-
-
-/* See if the three points are collinear, similar to GAL_POLYGON_LEFT_OF_LINE
-   except that the result has to be exactly zero. */
-#define GAL_POLYGON_COLLINEAR_WITH_LINE(A, B, C)                        \
-  (GAL_POLYGON_TRI_CROSS_PRODUCT((A), (B), (C)) > -GAL_POLYGON_ROUND_ERR \
-   && GAL_POLYGON_TRI_CROSS_PRODUCT((A), (B), (C)) < GAL_POLYGON_ROUND_ERR) /* 
== 0 */
-
-
-
-
-/* Similar to GAL_POLYGON_LEFT_OF_LINE except that if they are on the same 
line,
-   this will return 0 (so that it is not on the left). Therefore the
-   name is "proper left". */
-#define GAL_POLYGON_PROP_LEFT_OF_LINE(A, B, C)                          \
-  ( GAL_POLYGON_TRI_CROSS_PRODUCT((A), (B), (C)) > GAL_POLYGON_ROUND_ERR ) /* 
> 0   */
-
-
-#define GAL_POLYGON_MIN_OF_TWO(A, B) ((A)<(B)+GAL_POLYGON_ROUND_ERR ? (A) : 
(B))
-#define GAL_POLYGON_MAX_OF_TWO(A, B) ((A)>(B)-GAL_POLYGON_ROUND_ERR ? (A) : 
(B))
-
-
-
 __END_C_DECLS    /* From C++ preparations */
 
 #endif           /* __GAL_POLYGON_H__ */
diff --git a/lib/polygon.c b/lib/polygon.c
index a8224fa..11b134d 100644
--- a/lib/polygon.c
+++ b/lib/polygon.c
@@ -33,54 +33,90 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/polygon.h>
 
 
+
+
+
 /***************************************************************/
-/**************       Basic operations        ******************/
+/**************            MACROS             ******************/
 /***************************************************************/
+/* The cross product of two points from the center. */
+#define GAL_POLYGON_CROSS_PRODUCT(A, B) ( (A)[0]*(B)[1] - (B)[0]*(A)[1] )
+
+
+
+
+/* Find the cross product (2*area) between three points. Each point is
+   assumed to be a pointer that has atleast two values within it. */
+#define GAL_POLYGON_TRI_CROSS_PRODUCT(A, B, C)    \
+  ( ( (B)[0]-(A)[0] ) * ( (C)[1]-(A)[1] ) -       \
+    ( (C)[0]-(A)[0] ) * ( (B)[1]-(A)[1] ) )       \
+
+
+
+
 
-/* 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. Depending on the
-   transformation, the corners can have practically any order even if
-   before the transformation, they were ordered.
-
-   The input is an array containing the coordinates (two values) of
-   each corner. `n' is the number of corners. So the length of the
-   input should be 2*n. The output is an array with `n' elements
-   specifying the indexs in order. The reason the indexes are output
-   is that for all the pixels in the image, in a homographic
-   transform, the order is the same. So the input is unchanged, only
-   `n' values will be put in the ordinds array. Such that calling the
-   input coordinates in the following fashion will give an
-   anti-clockwise order for 4 points for example:
-
-   1st vertice: in[ordinds[0]*2], in[ordinds[0]*2+1]
-   2nd vertice: in[ordinds[1]*2], in[ordinds[1]*2+1]
-   3rd vertice: in[ordinds[2]*2], in[ordinds[2]*2+1]
-   4th vertice: in[ordinds[3]*2], in[ordinds[3]*2+1]
-
-   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!
-
-                 Concave Polygon        Convex Polygon
-
-                  D --------C          D------------- C
-                    \      |         E /            |
-                     \E    |           \            |
-                     /     |            \           |
-                   A--------B             A ---------B
-
-   This is because we are always going to be calculating the area of
-   the overlap between a quadrilateral and the pixel grid or the
-   quadrilaterral its self.
-
-   GAL_POLYGON_MAX_CORNERS is defined so there will be no need to allocate
-   these temporary arrays seprately. Since we are dealing with pixels,
-   the polygon can't really have too many vertices.
+/* We have the line A-B. We want to see if C is to the left of this
+   line or to its right. This function will return 1 if it is to the
+   left. It uses the basic property of vector multiplication: If the
+   three points are anti-clockwise (the point is to the left), then
+   the vector multiplication is positive, if it is negative, then it
+   is clockwise (c is to the right).
+
+   Ofcourse it is very important that A be below or equal to B in both
+   the X and Y directions. The rounding error might give
+   -0.0000000000001 (I didn't count the number of zeros!!) instead of
+   zero for the area. Zero would indicate that they are on the same
+   line in this case this should give a true result.
 */
+#define GAL_POLYGON_LEFT_OF_LINE(A, B, C)                               \
+  ( GAL_POLYGON_TRI_CROSS_PRODUCT((A), (B), (C)) > -GAL_POLYGON_ROUND_ERR ) /* 
>= 0 */
+
+
+
+
+/* See if the three points are collinear, similar to GAL_POLYGON_LEFT_OF_LINE
+   except that the result has to be exactly zero. */
+#define GAL_POLYGON_COLLINEAR_WITH_LINE(A, B, C)                        \
+  (GAL_POLYGON_TRI_CROSS_PRODUCT((A), (B), (C)) > -GAL_POLYGON_ROUND_ERR \
+   && GAL_POLYGON_TRI_CROSS_PRODUCT((A), (B), (C)) < GAL_POLYGON_ROUND_ERR) /* 
== 0 */
+
+
+
+
+/* Similar to GAL_POLYGON_LEFT_OF_LINE except that if they are on the same 
line,
+   this will return 0 (so that it is not on the left). Therefore the
+   name is "proper left". */
+#define GAL_POLYGON_PROP_LEFT_OF_LINE(A, B, C)                          \
+  ( GAL_POLYGON_TRI_CROSS_PRODUCT((A), (B), (C)) > GAL_POLYGON_ROUND_ERR ) /* 
> 0   */
+
+
+#define GAL_POLYGON_MIN_OF_TWO(A, B) ((A)<(B)+GAL_POLYGON_ROUND_ERR ? (A) : 
(B))
+#define GAL_POLYGON_MAX_OF_TWO(A, B) ((A)>(B)-GAL_POLYGON_ROUND_ERR ? (A) : 
(B))
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/***************************************************************/
+/**************       Basic operations        ******************/
+/***************************************************************/
+
+/* Sort the pixels in anti clock-wise order.*/
 void
 gal_polygon_ordered_corners(double *in, size_t n, size_t *ordinds)
 {