gnuastro-commits
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[gnuastro-commits] master d2e2c22 5/5: Imported recent work on enabling


From: Mohammad Akhlaghi
Subject: [gnuastro-commits] master d2e2c22 5/5: Imported recent work on enabling concave polygons
Date: Sat, 21 Mar 2020 17:59:05 -0400 (EDT)

branch: master
commit d2e2c22912324fcf1a2134e2a0ea0bd3ebf7a034
Merge: 51ee136 729d723
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>

    Imported recent work on enabling concave polygons
    
    Sachin started his branch with his two new functions, then I concluded it
    by enabling concave polygon functionalty in Crop.
---
 NEWS                                              |  21 ++-
 bin/crop/args.h                                   |  19 ++-
 bin/crop/main.h                                   |   3 +-
 bin/crop/onecrop.c                                |  32 ++--
 bin/crop/ui.c                                     |   8 +-
 bin/crop/ui.h                                     |   3 +-
 doc/gnuastro.texi                                 | 178 ++++++++++++----------
 lib/gnuastro/polygon.h                            |   9 +-
 lib/polygon.c                                     |  98 +++++++++++-
 tests/Makefile.am                                 |   4 +-
 tests/crop/{imgoutpolygon.sh => imgpolygonout.sh} |   2 +-
 11 files changed, 274 insertions(+), 103 deletions(-)

diff --git a/NEWS b/NEWS
index 2c11e68..846761a 100644
--- a/NEWS
+++ b/NEWS
@@ -17,6 +17,15 @@ See the end of the file for license conditions.
      is convenient when you forget the specific name of the spectral line
      used within Gnuastro.
 
+  Crop:
+   --polygon: can now also crop concave polygons (when atleast one inner
+     angle is more than 180 degrees).
+   --polygonnosort: Don't sort the given set of vertices to the `--polygon'
+     option. For a concave polyton, this is redundant, but for a convex
+     polygon, this option may be necessary because there is no unique
+     solution and if automatic sorting is not disabled, the output polygon
+     may not correspond to what you wanted (the vertices will the correct).
+
   Fits:
    --datasum: Calculate and print the given HDU's "datasum" to stdout.
    --datetosec: Can also account for `Z' in the end of the date-time
@@ -40,8 +49,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
 
@@ -63,6 +74,11 @@ See the end of the file for license conditions.
      NaN. Until now, the corresponding program/library would abort,
      printing the problematic string and its location.
 
+  Crop:
+   - `--polygonout' is the new name for `--outpolygon'. Having `polygon' at
+     the start of the option name, makes it easier to find in the help list
+     and also to understand generally.
+
   NoiseChisel:
    - Until now, when NoiseChisel didn't detect any pixels, it just printed
      a message and wouldn't not make any output file. This was very
@@ -71,6 +87,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/bin/crop/args.h b/bin/crop/args.h
index 441fe57..1105af4 100644
--- a/bin/crop/args.h
+++ b/bin/crop/args.h
@@ -268,13 +268,26 @@ struct argp_option program_options[] =
       GAL_OPTIONS_NOT_SET
     },
     {
-      "outpolygon",
-      UI_KEY_OUTPOLYGON,
+      "polygonout",
+      UI_KEY_POLYGONOUT,
       0,
       0,
       "Keep the polygon's outside, mask the inside.",
       UI_GROUP_REGION,
-      &p->outpolygon,
+      &p->polygonout,
+      GAL_OPTIONS_NO_ARG_TYPE,
+      GAL_OPTIONS_RANGE_0_OR_1,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
+    },
+    {
+      "polygonnosort",
+      UI_KEY_POLYGONNOSORT,
+      0,
+      0,
+      "Don't sort polygon vertices (are anticlockwise).",
+      UI_GROUP_REGION,
+      &p->polygonnosort,
       GAL_OPTIONS_NO_ARG_TYPE,
       GAL_OPTIONS_RANGE_0_OR_1,
       GAL_OPTIONS_NOT_MANDATORY,
diff --git a/bin/crop/main.h b/bin/crop/main.h
index 924d731..6620098 100644
--- a/bin/crop/main.h
+++ b/bin/crop/main.h
@@ -96,7 +96,8 @@ struct cropparams
   gal_list_str_t     *coordcol;  /* Column in cat containing coordinates. */
   char                *section;  /* Section string.                       */
   char                *polygon;  /* Input string of polygon vertices.     */
-  uint8_t           outpolygon;  /* ==1: Keep the inner polygon region.   */
+  uint8_t           polygonout;  /* ==1: Keep the inner polygon region.   */
+  uint8_t        polygonnosort;  /* Don't sort polygon vertices.          */
 
   /* Internal */
   size_t                 numin;  /* Number of input images.               */
diff --git a/bin/crop/onecrop.c b/bin/crop/onecrop.c
index bfeb568..31c5855 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((*isinside)(ipolygon, point, nvertices)==polygonout)         \
           ba[i]=*bb;                                                    \
       }                                                                 \
     free(bb);                                                           \
@@ -347,7 +347,8 @@ polygonmask(struct onecropparams *crp, void *array, long 
*fpixel_i,
 {
   int type=crp->p->type;
   double *ipolygon, point[2];
-  int outpolygon=crp->p->outpolygon;
+  int polygonout=crp->p->polygonout;
+  int (*isinside)(double *, double *, size_t);
   size_t i, *ordinds, size=s0*s1, nvertices=crp->p->nvertices;
 
 
@@ -365,15 +366,28 @@ polygonmask(struct onecropparams *crp, void *array, long 
*fpixel_i,
 
 
   /* Find the order of the polygons and put the elements in the proper
-     order. Also subtract the fpixel_i coordinates from all the
-     vertices to bring them into the crop image coordinates.*/
-  gal_polygon_ordered_corners(crp->ipolygon, crp->p->nvertices, ordinds);
+     order.*/
+  if(crp->p->polygonnosort)
+    for(i=0;i<crp->p->nvertices;++i) ordinds[i]=i;
+  else
+    gal_polygon_ordered_corners(crp->ipolygon, crp->p->nvertices, ordinds);
+
+  /* Fill the final polygon vertice positions within `ipolygon' and also
+     the fpixel_i coordinates from all the vertices to bring them into the
+     crop image coordinates. */
   for(i=0;i<crp->p->nvertices;++i)
     {
       ipolygon[i*2  ] = crp->ipolygon[ordinds[i]*2]   - fpixel_i[0];
       ipolygon[i*2+1] = crp->ipolygon[ordinds[i]*2+1] - fpixel_i[1];
     }
 
+  /* Set the function for checking if a point is inside the polygon. For
+     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 );
 
   /* Go over all the pixels in the image and if they are within the
      polygon keep them if the user has asked for it.*/
@@ -521,7 +535,7 @@ onecrop_flpixel(struct onecropparams *crp)
         onecrop_parse_section(p, dsize, fpixel, lpixel);
       else if(p->polygon)       /* Defined by a polygon.  */
         {
-          if(p->outpolygon==0)
+          if(p->polygonout==0)
             onecrop_ipolygon_fl(p->ipolygon, p->nvertices, fpixel, lpixel);
         }
       else                      /* Defined by its center. */
@@ -537,7 +551,7 @@ onecrop_flpixel(struct onecropparams *crp)
         {
           /* Fill crp->ipolygon in wcspolygonpixel, then set flpixel*/
           fillcrpipolygon(crp);
-          if(p->outpolygon==0)
+          if(p->polygonout==0)
             onecrop_ipolygon_fl(crp->ipolygon, p->nvertices, fpixel, lpixel);
         }
       else
@@ -564,7 +578,7 @@ onecrop_flpixel(struct onecropparams *crp)
 
   /* If the user only wants regions outside to the polygon, then set
      the fpixel and lpixel to cover the full input image. */
-  if(p->polygon && p->outpolygon)
+  if(p->polygon && p->polygonout)
     {
       crp->fpixel[0]=crp->fpixel[1]=1;
       crp->lpixel[0]=dsize[1];
@@ -837,7 +851,7 @@ onecrop(struct onecropparams *crp)
       free(array);
     }
   else
-    if(p->polygon && p->outpolygon==0 && p->mode==IMGCROP_MODE_WCS)
+    if(p->polygon && p->polygonout==0 && p->mode==IMGCROP_MODE_WCS)
       free(crp->ipolygon);
 
   /* The crop is complete. */
diff --git a/bin/crop/ui.c b/bin/crop/ui.c
index 86cd039..5806708 100644
--- a/bin/crop/ui.c
+++ b/bin/crop/ui.c
@@ -377,10 +377,10 @@ ui_read_check_only_options(struct cropparams *p)
       if(p->nvertices<3)
         error(EXIT_FAILURE, 0, "a polygon has to have 3 or more vertices, "
               "you have only given %zu (%s)", p->nvertices, p->polygon);
-      if(p->outpolygon && p->numin>1)
-        error(EXIT_FAILURE, 0, "currently in WCS mode, outpolygon can only "
-              "be set to zero when there is one image, you have given %zu "
-              "images. For multiple images the region will be very large. "
+      if(p->polygonout && p->numin>1)
+        error(EXIT_FAILURE, 0, "currently in WCS mode, `--polygonout' can "
+              "only be set to zero when there is one image, you have given "
+              "%zu images. For multiple images the region will be very large. "
               "It is best if you first crop out the larger region you want "
               "into one image, then mask the polygon", p->numin);
     }
diff --git a/bin/crop/ui.h b/bin/crop/ui.h
index 58138e2..728a15b 100644
--- a/bin/crop/ui.h
+++ b/bin/crop/ui.h
@@ -68,7 +68,8 @@ enum option_keys_enum
   UI_KEY_CATHDU         = 1000,
   UI_KEY_HSTARTWCS,
   UI_KEY_HENDWCS,
-  UI_KEY_OUTPOLYGON,
+  UI_KEY_POLYGONOUT,
+  UI_KEY_POLYGONNOSORT,
   UI_KEY_CHECKCENTER,
 };
 
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index f388b9d..5cc057e 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -9482,7 +9482,7 @@ In Image mode there are two options to define the 
vertices of a region to crop:
 The former is lower-level (doesn't accept floating point vertices, and only a 
rectangular region can be defined), it is also only available in Image mode.
 Please see @ref{Crop section syntax} for a full description of this method.
 
-The latter option (@option{--polygon}) is a higher-level method to define any 
convex polygon (with any number of vertices) with floating point values.
+The latter option (@option{--polygon}) is a higher-level method to define any 
polygon (with any number of vertices) with floating point values.
 Please see the description of this option in @ref{Invoking astcrop} for its 
syntax.
 @end table
 
@@ -9544,7 +9544,9 @@ This string is given to the @option{--section} option.
 Therefore, the pixels along the first axis that are @mymath{\geq}@command{X1} 
and @mymath{\leq}@command{X2} will be included in the cropped image.
 The same goes for the second axis.
 Note that each different term will be read as an integer, not a float.
-This is a low-level option, for a higher-level way to specify region (any 
polygon, not just a box), please see the @option{--polygon} option in @ref{Crop 
options}.
+
+The reason it only accepts integers is that @option{--section} is a low-level 
option (which is also very fast!).
+For a higher-level way to specify region (any polygon, not just a box), please 
see the @option{--polygon} option in @ref{Crop options}.
 Also note that in the FITS standard, pixel indexes along each axis start from 
unity(1) not zero(0).
 
 @cindex Crop section format
@@ -9714,14 +9716,34 @@ If you want an even sided crop, you can run Crop 
afterwards with @option{--secti
 
 @item -l STR
 @itemx --polygon=STR
-String of crop polygon vertices.
-Note that currently only convex polygons should be used.
-In the future we will make it work for all kinds of polygons.
-Convex polygons are polygons that do not have an internal angle more than 180 
degrees.
+String of vertices to define a polygon to crop.
+For a convex polygon, the vertices can be in any order: Crop will 
automatically sort them in an anti-clockwise fashion and use it.
+It also attempts to sort the vertices of a concave polygon, but in a concave 
polygon, there is no unique way to sort the vertices, so the output may not be 
what you want.
+In such scenarios, it may be better to ensure that the list of vertices is 
already sorted the way you like in an anti-clockwise manner and use the 
@option{--polygonnosort} option to disable the automatic sorting.
+
+@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
+
 This option can be used both in the image and WCS modes, see @ref{Crop modes}.
 The cropped image will be the size of the rectangular region that completely 
encompasses the polygon.
 By default all the pixels that are outside of the polygon will be set as blank 
values (see @ref{Blank pixels}).
-However, if @option{--outpolygon} is called all pixels internal to the 
vertices will be set to blank.
+However, if @option{--polygonout} is called all pixels internal to the 
vertices will be set to blank.
+In WCS-mode, you may provide many FITS images/tiles: Crop will stitch them to 
produce this cropped region, then apply the polygon.
 
 The syntax for the polygon vertices is similar to, and simpler than, that for 
@option{--section}.
 In short, the dimensions of each coordinate are separated by a comma (@key{,}) 
and each vertex is separated by a colon (@key{:}).
@@ -9736,9 +9758,8 @@ For example, let's assume you want to work on the deepest 
part of the WFC3/IR im
   (53.134517,-27.787144), (53.161906,-27.807208) ]
 @end example
 
-They have provided mask images with only these pixels in the WFC3/IR images, 
but what if you also need to work on the same region in the full resolution ACS 
images? Also what if you want to use the CANDELS data for the shallow region? 
Running Crop with @option{--polygon} will easily pull out this region of the 
image for you irrespective of the resolution.
+They have provided mask images with only these pixels in the WFC3/IR images, 
but what if you also need to work on the same region in the full resolution ACS 
images? Also what if you want to use the CANDELS data for the shallow region? 
Running Crop with @option{--polygon} will easily pull out this region of the 
image for you, irrespective of the resolution.
 If you have set the operating mode to WCS mode in your nearest configuration 
file (see @ref{Configuration files}), there is no need to call 
@option{--mode=wcs} on the command line.
-You may also provide many FITS images/tiles and Crop will stitch them to 
produce this cropped region:
 
 @example
 $ astcrop --mode=wcs desired-filter-image(s).fits           \
@@ -9766,7 +9787,7 @@ $ v=$(awk 'NR==4' ds9.reg | sed -e's/polygon(//'        \
 $ astcrop --mode=wcs image.fits --polygon=$v
 @end example
 
-@item --outpolygon
+@item --polygonout
 Keep all the regions outside the polygon and mask the inner ones with blank 
pixels (see @ref{Blank pixels}).
 This is practically the inverse of the default mode of treating polygons.
 Note that this option only works when you have only provided one input image.
@@ -9775,6 +9796,11 @@ This can lead to a very large area if large surveys like 
COSMOS are used.
 So Crop will abort and notify you.
 In such cases, it is best to crop out the larger region you want, then mask 
the smaller region with this option.
 
+@item --polygonnosort
+Don't sort the given set of vertices to the @option{--polygon} option.
+For a concave polyton, this is redundant, but for a convex polygon, this 
option may be necessary because there is no unique for a concave polygon.
+If automatic sorting is not disabled, the output polygon may not correspond to 
what you had in mind, the vertices will be correct though.
+
 @item -s STR
 @itemx --section=STR
 Section of the input image which you want to be cropped.
@@ -23083,16 +23109,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.
@@ -23100,32 +23139,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]
@@ -23136,55 +23169,46 @@ 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
-
-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 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.
 
-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.
+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
 
+@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_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.
+@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}.
+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})
-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 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 d4f24d0..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,11 +61,17 @@ __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_pin(double *v, double *p, size_t n);
+gal_polygon_isinside(double *v, double *p, 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);
diff --git a/lib/polygon.c b/lib/polygon.c
index 4c51223..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,6 +252,59 @@ gal_polygon_area(double *v, 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--;
+      }
+
+      /* Increment `j' */
+      j=i++;
+
+      /* For a check:
+         printf("winding number: %ld, %.3f\n", wn);
+      */
+    }
+  return wn;
+}
+
+
+
+
+
 /* 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
@@ -225,7 +317,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 +334,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)
 {
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 9a0d857..6dd1e3d 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -98,7 +98,7 @@ endif
 if COND_CROP
   MAYBE_CROP_TESTS = crop/imgcat.sh crop/wcscat.sh crop/imgcenter.sh    \
   crop/imgcenternoblank.sh crop/section.sh crop/wcscenter.sh            \
-  crop/imgpolygon.sh crop/imgoutpolygon.sh crop/wcspolygon.sh
+  crop/imgpolygon.sh crop/imgpolygonout.sh crop/wcspolygon.sh
 
   crop/imgcat.sh: mkprof/mosaic1.sh.log
   crop/wcscat.sh: mkprof/mosaic1.sh.log mkprof/mosaic2.sh.log     \
@@ -109,7 +109,7 @@ if COND_CROP
   crop/wcscenter.sh: mkprof/mosaic1.sh.log mkprof/mosaic2.sh.log      \
                      mkprof/mosaic3.sh.log mkprof/mosaic4.sh.log
   crop/imgpolygon.sh: mkprof/mosaic1.sh.log
-  crop/imgoutpolygon.sh: mkprof/mosaic1.sh.log
+  crop/imgpolygonout.sh: mkprof/mosaic1.sh.log
   crop/wcspolygon.sh: mkprof/mosaic1.sh.log mkprof/mosaic2.sh.log \
                       mkprof/mosaic3.sh.log mkprof/mosaic4.sh.log
 endif
diff --git a/tests/crop/imgoutpolygon.sh b/tests/crop/imgpolygonout.sh
similarity index 96%
rename from tests/crop/imgoutpolygon.sh
rename to tests/crop/imgpolygonout.sh
index a6d04ed..722a4bf 100755
--- a/tests/crop/imgoutpolygon.sh
+++ b/tests/crop/imgpolygonout.sh
@@ -59,5 +59,5 @@ if [ ! -f $img      ]; then echo "$img does not exist.";   
exit 77; fi
 # string. Such programs will execute the command if present and help in
 # debugging when the developer doesn't have access to the user's system.
 $check_with_program $execname $img $cat --mode=img --zeroisnotblank     \
-                              --outpolygon --output=imgoutpolygon.fits  \
+                              --polygonout --output=imgpolygonout.fits  \
                               
--polygon=209,50:436.76,151:475.64,438.2:210.6,454.04:121.4,289.88



reply via email to

[Prev in Thread] Current Thread [Next in Thread]