[gnuastro-commits] master e6d7c9a 1/2: gal_label_oversegment now checks for blank values

[Mohammad Akhlaghi]

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

    gal_label_oversegment now checks for blank values
    
    Until now, `gal_label_oversegment' wouldn't check for blank values,
    therefore, if the user hadn't already set the blank value flag, it wouldn't
    put rivers on pixels touching blank values. So a call to
    `gal_blank_present' was added to this function and a single variable is now
    set to see if blank values are present or not (making it easier to read
    than bit-flag checking deep into the loop).
    
    Also, the description in the bit flag `GAL_DATA_FLAG_BLANK_CH' was edited
    to be more clear.
---
 doc/gnuastro.texi | 38 ++++++++++++++++++++++++++++++++------
 lib/label.c       | 51 ++++++++++++++++++++++++++++-----------------------
 2 files changed, 60 insertions(+), 29 deletions(-)

diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index fac52c5..cf67d25 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -21051,11 +21051,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
address@hidden 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
@@ -21066,7 +21073,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
address@hidden are meaningful.
address@hidden 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
@@ -25466,6 +25474,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
address@hidden) 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
address@hidden in @ref{Generic data container}):
address@hidden
+input->flags |=  GAL_DATA_FLAG_BLANK_CH; /* Set bit to 1. */
+input->flags &= ~GAL_DATA_FLAG_HASBLANK; /* Set bit to 0. */
address@hidden 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