[gnuastro-commits] master 729d723 4/5: Crop: polygon vertices can now be concave (inner angles >180 deg)

[Mohammad Akhlaghi]

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

    Crop: polygon vertices can now be concave (inner angles >180 deg)
    
    Until now, the `--polygon' option of the Crop program could only work on
    convex polygons (where all inner points have an angle less than
    180). However, the deep fields of many astronomical surveys
    (stacked/coadded images in general), are a concave polygon and it was not
    possible to extract them accurately.
    
    With this commit, thanks to the two new functions that Sachin added a few
    commits earlier in the polygon library, it is now possible to also extract
    convex polygons from an image also. But Concave polygons are not unique, so
    it was also necessary to add the new `--polygonnosort' option to crop. With
    this option, the polygon vertices are no longer sorted, allowing users to
    extract the exact concave polygon that they want.
---
 NEWS                                              | 14 +++++++
 bin/crop/args.h                                   | 19 ++++++++--
 bin/crop/main.h                                   |  3 +-
 bin/crop/onecrop.c                                | 30 +++++++++------
 bin/crop/ui.c                                     |  8 ++--
 bin/crop/ui.h                                     |  3 +-
 doc/gnuastro.texi                                 | 46 ++++++++++++++++++-----
 tests/Makefile.am                                 |  4 +-
 tests/crop/{imgoutpolygon.sh => imgpolygonout.sh} |  2 +-
 9 files changed, 95 insertions(+), 34 deletions(-)

diff --git a/NEWS b/NEWS
index 5d97a7a..d4a660b 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
@@ -55,6 +64,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
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 54cd9ad..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((*isinside)(ipolygon, point, nvertices)==outpolygon)           \
+        if((*isinside)(ipolygon, point, nvertices)==polygonout)         \
           ba[i]=*bb;                                                    \
       }                                                                 \
     free(bb);                                                           \
@@ -347,7 +347,7 @@ 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;
 
@@ -366,9 +366,15 @@ 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];
@@ -379,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(crp->ipolygon, crp->p->nvertices)
-               ? &gal_polygon_isinside_convex
-               : &gal_polygon_isinside );
+  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.*/
@@ -529,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. */
@@ -545,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
@@ -572,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];
@@ -845,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 3bd5d4f..3921374 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -9483,7 +9483,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
 
@@ -9545,7 +9545,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
@@ -9715,14 +9717,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{:}).
@@ -9737,9 +9759,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           \
@@ -9767,7 +9788,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.
@@ -9776,6 +9797,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.
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