[gnuastro-commits] master d13c715 065/113: Imported work in master, no conflicts

[Mohammad Akhlaghi]

branch: master
commit d13c715ee4e7dca63b7e487a1a306a97e78afaa8
Merge: 2491c91 81252e2
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    Imported work in master, no conflicts
    
    There weren't any conflicts in this merge.
---
 NEWS                          |  4 +++
 bin/noisechisel/args.h        | 13 ++++++++
 bin/noisechisel/main.h        |  1 +
 bin/noisechisel/noisechisel.c |  8 ++---
 bin/noisechisel/ui.h          |  1 +
 bin/statistics/args.h         | 13 ++++++++
 bin/statistics/main.h         |  1 +
 bin/statistics/sky.c          | 17 ++++++++--
 bin/statistics/ui.h           |  1 +
 doc/gnuastro.texi             | 76 +++++++++++++++++++++++++++++++++++--------
 lib/label.c                   | 51 ++++++++++++++++-------------
 11 files changed, 142 insertions(+), 44 deletions(-)

diff --git a/NEWS b/NEWS
index ae73509..ba39e32 100644
--- a/NEWS
+++ b/NEWS
@@ -55,6 +55,9 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
 
   NoiseChisel:
     --rawoutput: only output the detection labels and Sky and its STD.
+    --ignoreblankinsky: don't set the pixels that are blank in the input to
+      blank in the Sky and Sky standard deviation outputs (when
+      `--oneelempertile' is not called).
     --label: label the connected detections. Until now this was the default
       behavior. However, from this release, NoiseChisel is only in charge
       of detection. Segmentation is done by a new program (Segment). Since
@@ -65,6 +68,7 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
 
   Statistics:
     --manualbinrange: histogram or CFP range can be outside of distribution.
+    --ignoreblankinsky: similar to same option in NoiseChisel.
 
   Libraries:
     gal_array_read: read array from any of known formats (FITS,TIFF,JPEG,...).
diff --git a/bin/noisechisel/args.h b/bin/noisechisel/args.h
index 3e6e159..e96b73d 100644
--- a/bin/noisechisel/args.h
+++ b/bin/noisechisel/args.h
@@ -137,6 +137,19 @@ struct argp_option program_options[] =
 
     /* Output options. */
     {
+      "ignoreblankinsky",
+      UI_KEY_IGNOREBLANKINSKY,
+      0,
+      0,
+      "Don't write input's blanks in the Sky output.",
+      GAL_OPTIONS_GROUP_OUTPUT,
+      &p->ignoreblankinsky,
+      GAL_OPTIONS_NO_ARG_TYPE,
+      GAL_OPTIONS_RANGE_0_OR_1,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
+    },
+    {
       "rawoutput",
       UI_KEY_RAWOUTPUT,
       0,
diff --git a/bin/noisechisel/main.h b/bin/noisechisel/main.h
index 05eb729..a2590bc 100644
--- a/bin/noisechisel/main.h
+++ b/bin/noisechisel/main.h
@@ -52,6 +52,7 @@ struct noisechiselparams
   char                  *whdu;  /* Wide kernel HDU.                       */
 
   uint8_t  continueaftercheck;  /* Don't abort after the check steps.     */
+  uint8_t    ignoreblankinsky;  /* Ignore input's blank values.           */
   uint8_t           rawoutput;  /* Only detection & 1 elem/tile output.   */
   uint8_t               label;  /* Label detections that are connected.   */
 
diff --git a/bin/noisechisel/noisechisel.c b/bin/noisechisel/noisechisel.c
index 04f6835..060b865 100644
--- a/bin/noisechisel/noisechisel.c
+++ b/bin/noisechisel/noisechisel.c
@@ -182,8 +182,8 @@ noisechisel_output(struct noisechiselparams *p)
   /* Write the Sky image into the output */
   if(p->sky->name) free(p->sky->name);
   p->sky->name="SKY";
-  gal_tile_full_values_write(p->sky, &p->cp.tl, 1, p->cp.output,
-                             NULL, PROGRAM_NAME);
+  gal_tile_full_values_write(p->sky, &p->cp.tl, !p->ignoreblankinsky,
+                             p->cp.output, NULL, PROGRAM_NAME);
   p->sky->name=NULL;
 
 
@@ -198,8 +198,8 @@ noisechisel_output(struct noisechiselparams *p)
   gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "MEDSTD", 0, &p->medstd, 0,
                         "Median raw tile standard deviation", 0,
                         p->input->unit);
-  gal_tile_full_values_write(p->std, &p->cp.tl, 1, p->cp.output, keys,
-                             PROGRAM_NAME);
+  gal_tile_full_values_write(p->std, &p->cp.tl, !p->ignoreblankinsky,
+                             p->cp.output, keys, PROGRAM_NAME);
   p->std->name=NULL;
 
   /* Let the user know that the output is written. */
diff --git a/bin/noisechisel/ui.h b/bin/noisechisel/ui.h
index 0c5e7b2..c3e72b2 100644
--- a/bin/noisechisel/ui.h
+++ b/bin/noisechisel/ui.h
@@ -94,6 +94,7 @@ enum option_keys_enum
   UI_KEY_CHECKDETECTION,
   UI_KEY_CHECKSKY,
   UI_KEY_RAWOUTPUT,
+  UI_KEY_IGNOREBLANKINSKY,
 };
 
 
diff --git a/bin/statistics/args.h b/bin/statistics/args.h
index 6b0fbb9..47beddf 100644
--- a/bin/statistics/args.h
+++ b/bin/statistics/args.h
@@ -520,6 +520,19 @@ struct argp_option program_options[] =
       GAL_OPTIONS_NOT_MANDATORY,
       GAL_OPTIONS_NOT_SET
     },
+    {
+      "ignoreblankinsky",
+      UI_KEY_IGNOREBLANKINSKY,
+      0,
+      0,
+      "Don't write input's blanks in the Sky output.",
+      UI_GROUP_SKY,
+      &p->ignoreblankinsky,
+      GAL_OPTIONS_NO_ARG_TYPE,
+      GAL_OPTIONS_RANGE_0_OR_1,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
+    },
 
 
 
diff --git a/bin/statistics/main.h b/bin/statistics/main.h
index 01c437b..fe58f90 100644
--- a/bin/statistics/main.h
+++ b/bin/statistics/main.h
@@ -90,6 +90,7 @@ struct statisticsparams
   size_t       smoothwidth;  /* Width of flat kernel to smooth interpd.  */
   uint8_t         checksky;  /* Save the steps for deriving the Sky.     */
   double    sclipparams[2];  /* Muliple and parameter of sigma clipping. */
+  uint8_t ignoreblankinsky;  /* Ignore input's blank values.             */
 
 
   /* Internal */
diff --git a/bin/statistics/sky.c b/bin/statistics/sky.c
index ffd56c1..f50359d 100644
--- a/bin/statistics/sky.c
+++ b/bin/statistics/sky.c
@@ -131,6 +131,7 @@ sky(struct statisticsparams *p)
   char *msg, *outname;
   struct timeval t0, t1;
   gal_data_t *num, *tmp;
+  uint8_t keepinputdir=p->cp.keepinputdir;
   struct gal_options_common_params *cp=&p->cp;
   struct gal_tile_two_layer_params *tl=&cp->tl;
 
@@ -246,13 +247,23 @@ sky(struct statisticsparams *p)
     }
 
 
-  /* Save the Sky and its standard deviation */
+  /* Save the Sky and its standard deviation. We want the output to have a
+     `_sky.fits' suffix. So we'll temporarily re-set `p->cp.keepinputdir'
+     if the user asked for a specific name. Note that we copied the actual
+     value in the `keepinputdir' above (in the definition). */
+  p->cp.keepinputdir = p->cp.output ? 1 : keepinputdir;
   outname=gal_checkset_automatic_output(&p->cp,
                                         ( p->cp.output
                                           ? p->cp.output
                                           : p->inputname ), "_sky.fits");
-  gal_tile_full_values_write(p->sky_t, tl, 1, outname, NULL, PROGRAM_NAME);
-  gal_tile_full_values_write(p->std_t, tl, 1, outname, NULL, PROGRAM_NAME);
+  p->sky_t->name="SKY";
+  p->std_t->name="SKY_STD";
+  p->cp.keepinputdir=keepinputdir;
+  gal_tile_full_values_write(p->sky_t, tl, !p->ignoreblankinsky, outname,
+                             NULL, PROGRAM_NAME);
+  gal_tile_full_values_write(p->std_t, tl, !p->ignoreblankinsky, outname,
+                             NULL, PROGRAM_NAME);
+  p->sky_t->name = p->std_t->name = NULL;
   if(!cp->quiet)
     printf("  - Written to `%s'.\n", outname);
 
diff --git a/bin/statistics/ui.h b/bin/statistics/ui.h
index e8d4bef..13efee9 100644
--- a/bin/statistics/ui.h
+++ b/bin/statistics/ui.h
@@ -95,6 +95,7 @@ enum option_keys_enum
   UI_KEY_MODMEDQDIFF,
   UI_KEY_SMOOTHWIDTH,
   UI_KEY_CHECKSKY,
+  UI_KEY_IGNOREBLANKINSKY,
   UI_KEY_SCLIPPARAMS,
 };
 
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index f400017..b5da0f1 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -7956,14 +7956,20 @@ save them into another plain text or FITS table.
 @section Fits
 
 The ``Flexible Image Transport System'', or FITS, is by far the most common
-data container format in astronomy. Although the full name of the standard
-invokes the idea that it is only for images, it also contains very complete
-and robust features for tables. It started off in the 1970s and was
-formally published as a standard in 1981, it was adopted by the
-International Astronomical Union (IAU) in 1982 and an IAU working group to
-maintain its future was defined in 1988. The FITS 2.0 and 3.0 standards
-were approved in 2000 and 2008 respectively, and the 4.0 draft has also
-been released recently, please see the
+data container format in astronomy and in constant use since the
+1970s. Archivability (future usage, simplicity) has been one of the primary
+design principles of this format. In the last few decades it has proved so
+useful and robust that the Vatican Library has also chosen FITS for its
+``long-term digital preservation''
+project@footnote{@url{https://www.vaticanlibrary.va/home.php?pag=progettodigit}}.
+
+Although the full name of the standard invokes the idea that it is only for
+images, it also contains very complete and robust features for tables. It
+started off in the 1970s and was formally published as a standard in 1981,
+it was adopted by the International Astronomical Union (IAU) in 1982 and an
+IAU working group to maintain its future was defined in 1988. The FITS 2.0
+and 3.0 standards were approved in 2000 and 2008 respectively, and the 4.0
+draft has also been released recently, please see the
 @url{https://fits.gsfc.nasa.gov/fits_standard.html, FITS standard document
 webpage} for the full text of all versions. Also see the
 @url{https://doi.org/10.1051/0004-6361/201015362, FITS 3.0 standard paper}
@@ -13722,6 +13728,14 @@ discontinuities do not show up in the final Sky 
values. The smoothing is
 done through convolution (see @ref{Convolution process}) with a flat
 kernel, so the value to this option must be an odd number.
 
+@item --ignoreblankinsky
+Don't set the input's blank pixels to blank in the output Sky and Sky
+standard deviation datasets. This is only applicable when the output has
+the same size as the input, in other words, when @option{--oneelempertile}
+isn't called. By default, blank values in the input (commonly on the edges
+which are outside the survey/field area) will be set to blank in the output
+Sky and Sky standard deviation also.
+
 @item --checksky
 Create a multi-extension FITS file showing the steps that were used to
 estimate the Sky value over the input, see @ref{Quantifying signal in a
@@ -14686,6 +14700,14 @@ NoiseChisel will abort once its desired extensions 
have been written. With
 NoiseChisel to continue with the rest of the processing, even after the
 requested check files are complete.
 
+@item --ignoreblankinsky
+Don't set the input's blank pixels to blank in the output Sky and Sky
+standard deviation datasets. This is only applicable when the output has
+the same size as the input, in other words, when @option{--oneelempertile}
+isn't called. By default, blank values in the input (commonly on the edges
+which are outside the survey/field area) will be set to blank in the output
+Sky and Sky standard deviation also.
+
 @item -l
 @itemx --label
 Run a connected-components algorithm on the finally detected pixels to
@@ -21353,11 +21375,18 @@ these flags. The currently recognized bits are stored 
in these macros:
 
 @cindex Blank data
 @item GAL_DATA_FLAG_BLANK_CH
-Marking that the dataset has been checked for blank values. Therefore, the
-value of the bit in @code{GAL_DATA_FLAG_HASBLANK} is reliable. Without
-this bit, when a dataset doesn't have any blank values (and this has been
-checked), the @code{GAL_DATA_FLAG_HASBLANK} bit will be zero so a checker
-has no way to know if this zero is real or if no check has been done yet.
+Marking that the dataset has been checked for blank values or not. When a
+dataset doesn't have any blank values, the @code{GAL_DATA_FLAG_HASBLANK}
+bit will be zero. But upon initialization, all bits also get a value of
+zero. Therefore, a checker needs this flag to see if the value in
+@code{GAL_DATA_FLAG_HASBLANK} is reliable (dataset has actually been parsed
+for a blank value) or not.
+
+Also, if it is necessary to re-check the presence of flags, you just have
+to set this flag to zero and call @code{gal_blank_present} for example to
+parse the dataset and check for blank values. Note that for improved
+efficiency, when this flag is set, @code{gal_blank_present} will not
+actually parse the dataset, it will just use @code{GAL_DATA_FLAG_HASBLANK}.
 
 @item GAL_DATA_FLAG_HASBLANK
 This bit has a value of @code{1} when the given dataset has blank
@@ -21368,7 +21397,8 @@ values, so there is no more need for further checks.
 @item GAL_DATA_FLAG_SORT_CH
 Marking that the dataset is already checked for being sorted or not and
 thus that the possible @code{0} values in @code{GAL_DATA_FLAG_SORTED_I} and
-@code{GAL_DATA_FLAG_SORTED_D} are meaningful.
+@code{GAL_DATA_FLAG_SORTED_D} are meaningful. The logic behind this is
+similar to that in @code{GAL_DATA_FLAG_BLANK_CH}.
 
 @item GAL_DATA_FLAG_SORTED_I
 This bit has a value of @code{1} when the given dataset is sorted in an
@@ -25801,6 +25831,24 @@ must have a @code{GAL_TYPE_INT32} type and be the same 
size as
 @code{values}. All local minima (maxima), when @code{min0_max1} is @code{1}
 (@code{0}), are considered rivers (watersheds) and given a label of
 @code{GAL_LABEL_RIVER} (see above).
+
+Note that rivers/watersheds will also be formed on the edges of the labeled
+regions or when the labeled pixels touch a blank pixel. Therefore this
+function will need to check for the presence of blank values. To be most
+efficient, it is thus recommended to use @code{gal_blank_present} (with
+@code{updateflag=1}) prior to calling this function (see @ref{Library blank
+values}. Once the flag has been set, no other function (including this one)
+that needs special behavior for blank pixels will have to parse the dataset
+any more.
+
+If you are sure your dataset doesn't have blank values (by the design of
+your software), to avoid an extra parsing of the dataset and improve
+performance, you can set the two bits manually (see the description of
+@code{flags} in @ref{Generic data container}):
+@example
+input->flags |=  GAL_DATA_FLAG_BLANK_CH; /* Set bit to 1. */
+input->flags &= ~GAL_DATA_FLAG_HASBLANK; /* Set bit to 0. */
+@end example
 @end deftypefun
 
 @deftypefun void gal_label_grow_indexs (gal_data_t @code{*labels}, gal_data_t 
@code{*indexs}, int @code{withrivers}, int @code{connectivity})
diff --git a/lib/label.c b/lib/label.c
index 6e719cb..44f7b54 100644
--- a/lib/label.c
+++ b/lib/label.c
@@ -181,13 +181,13 @@ gal_label_oversegment(gal_data_t *values, gal_data_t 
*indexs,
 {
   size_t ndim=values->ndim;
 
+  int hasblank;
   float *arr=values->array;
   gal_list_sizet_t *Q=NULL, *cleanup=NULL;
   size_t *a, *af, ind, *dsize=values->dsize;
   size_t *dinc=gal_dimension_increment(ndim, dsize);
   int32_t n1, nlab, rlab, curlab=1, *labs=labels->array;
 
-
   /* Sanity checks */
   label_check_type(values, GAL_TYPE_FLOAT32, "values", __func__);
   label_check_type(indexs, GAL_TYPE_SIZE_T,  "indexs", __func__);
@@ -200,6 +200,10 @@ gal_label_oversegment(gal_data_t *values, gal_data_t 
*indexs,
           "%zuD", __func__, indexs->ndim);
 
 
+  /* See if there are blank values in the input dataset. */
+  hasblank=gal_blank_present(values, 0);
+
+
   /*********************************************
    For checks and debugging:*
   gal_data_t *crop;
@@ -328,18 +332,18 @@ gal_label_oversegment(gal_data_t *values, gal_data_t 
*indexs,
                                         ? (n1==nlab ? n1 : GAL_LABEL_RIVER)
                                         : nlab )
 
-                                    /* If the data has blank pixels (recall
-                                       that blank in int32 is negative),
-                                       see if the neighbor is blank and if
-                                       so, set the label to a river. Since
-                                       the flag checking can be done
-                                       outside this loop, for datasets with
-                                       no blank element this last step will
-                                       be completley ignored. */
-                                    : ( ( (values->flag
-                                           & GAL_DATA_FLAG_HASBLANK)
-                                          && nlab==GAL_BLANK_INT32 )
-                                        ? GAL_LABEL_RIVER : n1 ) );
+                                    /* If the data has blank pixels, see if
+                                       the neighbor is blank. If so, set
+                                       the label to a river. Checking for
+                                       the presence of blank values in the
+                                       dataset can be done outside this
+                                       loop (or even outside this function
+                                       if flags are set). So to help the
+                                       compiler optimize the program, we'll
+                                       first use the pre-checked value. */
+                                    : ( ( hasblank && isnan(arr[nind]) )
+                                        ? GAL_LABEL_RIVER
+                                        : n1 ) );
                            }
 
                          /* If this neigbour has a label of zero, then we
@@ -419,16 +423,17 @@ gal_label_oversegment(gal_data_t *values, gal_data_t 
*indexs,
                                     ? ( nlab==n1 ? n1 : GAL_LABEL_RIVER )
                                     : nlab )
 
-                                /* If the dataset has blank values and this
-                                   neighbor is blank, then the pixel should
-                                   be a river. Note that the blank checking
-                                   can be optimized out, so if the input
-                                   doesn't have blank values,
-                                   `nlab==GAL_BLANK_INT32' will never be
-                                   checked. */
-                                : ( (values->flag & GAL_DATA_FLAG_HASBLANK)
-                                    && nlab==GAL_BLANK_INT32
-                                    ? GAL_LABEL_RIVER : n1 ) )
+                                /* If the data has blank pixels, see if the
+                                   neighbor is blank. If so, set the label
+                                   to a river. Checking for the presence of
+                                   blank values in the dataset can be done
+                                   outside this loop (or even outside this
+                                   function if flags are set). So to help
+                                   the compiler optimize the program, we'll
+                                   first use the pre-checked value. */
+                                : ( ( hasblank && isnan(arr[nind]) )
+                                    ? GAL_LABEL_RIVER
+                                    : n1 ) )
 
                             /* `nlab==0' (the neighbor lies in the other
                                domain (sky or detections). To avoid the