[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastrocommits] master 7323d4b 1/2: polygon.h and qsort.h documented
From: 
Mohammad Akhlaghi 
Subject: 
[gnuastrocommits] master 7323d4b 1/2: polygon.h and qsort.h documented in the book 
Date: 
Wed, 21 Sep 2016 23:48:30 +0000 (UTC) 
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 nonrectangular
+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 Roundoff error
+We have to consider floating point roundoff errors when dealing with
+polygons. For example we will take @code{A} as the maximum of @code{A} and
address@hidden when @code{A>BGAL_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 anticlockwise 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  \ 
+ /  \ 
+ AB 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 anticlockwise 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,
+SutherlandHodgman} polygon clipping algorithm. Note that the vertices of
+both polygons have to be sorted in an anticlockwise 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 AB. 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 anticlockwise (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
 anticlockwise 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  \ 
 /  \ 
 AB 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 AB. 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 anticlockwise (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 clockwise order.*/
void
gal_polygon_ordered_corners(double *in, size_t n, size_t *ordinds)
{