[gnuastro-commits] master 1236154: NoiseChisel: good error message when quantile threshold not found

[Mohammad Akhlaghi]

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

    NoiseChisel: good error message when quantile threshold not found
    
    Until now, when no tiles could be used to find the quantile threshold,
    NoiseChisel would crash with an error saying it can't allocate a zero-sized
    array, not saying exactly why this was happening!
    
    We had a very complete and descriptive error message *after* removing
    outliers, but not before it. So unless the user ran NoiseChisel again with
    `--checkqthresh' and inspected the extensions, they couldn't see if this
    error occurred before the outlier rejection or after it.
    
    With this commit, a function is defined to print the elaborate error
    message while mentioning at which point the error happens. Its suggestions
    are also increased if the error message happened after the outlier
    rejection step (since there are two more parameters).
    
    The particular crash of this bug was actually happening because
    `gal_statistics_no_blank_sorted' was trying to copy a zero-sized array. So
    with this commit, in such situations (when `size==0', it will just make a
    blank array when `inplace==0' instead of attempting to copy the input.
    
    Also, to implement this feature, the new `gal_blank_number' function was
    defined (as a small wrapper around `gal_statistics_number') in the library.
    
    This bug was reported by Raúl Infante Sainz.
    
    This fixes bug #55491.
---
 NEWS                        |  2 +
 bin/noisechisel/threshold.c | 91 ++++++++++++++++++++++++++++++++-------------
 doc/gnuastro.texi           | 13 ++++++-
 lib/blank.c                 | 28 ++++++++++++++
 lib/gnuastro/blank.h        |  3 ++
 lib/statistics.c            | 13 +++++--
 6 files changed, 119 insertions(+), 31 deletions(-)

diff --git a/NEWS b/NEWS
index 0dd41c5..5389881 100644
--- a/NEWS
+++ b/NEWS
@@ -38,6 +38,7 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
   Library:
     GAL_BLANK_LONG: new macro for the `long' type blank value.
     GAL_BLANK_ULONG: new macro for the `unsigned long' type blank value.
+    gal_blank_number: Return the number of blank elements in a dataset.
     gal_statistics_outlier_cumulative: Uses flatness of the cumulative
        distribution to find outliers.
 
@@ -57,6 +58,7 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
   bug #55333: Fits program crash when --write or --update have no value.
   bug #55439: Arithmetic segmentation fault when reusing dataset name.
   bug #55478: Memory mapping crashes when .gnuastro is not writable.
+  bug #55491: NoiseChisel crash when no tiles good for quantile thresholding.
 
 
 
diff --git a/bin/noisechisel/threshold.c b/bin/noisechisel/threshold.c
index 1132742..11a56a6 100644
--- a/bin/noisechisel/threshold.c
+++ b/bin/noisechisel/threshold.c
@@ -518,13 +518,72 @@ qthresh_on_tile(void *in_prm)
 
 
 
+static void
+threshold_good_error(size_t number, int before0_after1, size_t interpnumngb)
+{
+  before0_after1=1;
+
+  /* Set the differing strings. */
+  char *in1 = ( before0_after1
+                ? "after removing outliers"
+                : "for defining a quantile threshold" );
+  char *in2 = ( before0_after1
+                ? ""
+                : "NOTE that this is happening *BEFORE* outlier rejection "
+                  "(where the number may decrease even further).");
+  char *in3 = ( before0_after1
+                ? "\n"
+                  "  - (slightly) Increase `--outliersclip' to reject less "
+                  "as outliers.\n"
+                  "  - (slightly) Increase `--outliersigma' to reject less "
+                  "as outliers.\n"
+                : "\n");
+
+  /* Print the error message and abort. */
+  error(EXIT_FAILURE, 0, "%zu tiles usable %s!\n\n"
+
+        "This is smaller than the requested minimum value of %zu (value to "
+        "the `--interpnumngb' option). %s\n\n"
+
+        "There are several ways to address the problem. The best and most "
+        "highly recommended is to use a larger input if possible (when the "
+        "input is a crop from a larger dataset). If this is not the case, "
+        "or it doesn't solve the problem, you need to loosen the "
+        "parameters mentioned below in the respective order (and therefore "
+        "cause scatter/inaccuracy in the final result). Hence its best to "
+        "not loosen them too much (recall that you can see all the option "
+        "values to Gnuastro's programs by appending `-P' to the end of your "
+        "command).\n"
+        "  - (slightly) Decrease `--tilesize' so your tile-grid has more "
+        "tiles.\n"
+        "  - (slightly) Increase `--meanmedqdiff' to accept more tiles.%s"
+        "  - (slightly) Decrease `--interpnumngb' to be less than %zu.\n\n"
+
+        "---- Tip ----\n"
+        "Append your command with `--checkqthresh' to see the "
+        "successful tiles in relation with this dataset's contents "
+        "before this crash. A visual inspection will greatly help in "
+        "finding the cause/solution for this particular dataset (note "
+        "that the output of `--checkqthresh' is a multi-extension FITS "
+        "file).\n\n"
+        "To better understand this important step, please run the "
+        "following command (press `SPACE'/arrow-keys to navigate and "
+        "`Q' to return back to the command-line):\n\n"
+        "    $ info gnuastro \"Quantifying signal in a tile\"\n", number,
+        in1, interpnumngb, in2, in3, number);
+}
+
+
+
+
+
 void
 threshold_quantile_find_apply(struct noisechiselparams *p)
 {
   char *msg;
-  size_t nval;
   gal_data_t *num;
   struct timeval t1;
+  size_t nval, nblank;
   struct qthreshparams qprm;
   struct gal_options_common_params *cp=&p->cp;
   struct gal_tile_two_layer_params *tl=&cp->tl;
@@ -601,6 +660,10 @@ threshold_quantile_find_apply(struct noisechiselparams *p)
         }
     }
 
+  /* A small sanity check. */
+  nblank=gal_blank_number(qprm.erode_th, 1);
+  if( nblank > qprm.erode_th->size-cp->interpnumngb )
+    threshold_good_error(qprm.erode_th->size-nblank, 0, cp->interpnumngb);
 
   /* Remove outliers if requested. */
   if(p->outliersigma!=0.0)
@@ -616,31 +679,7 @@ threshold_quantile_find_apply(struct noisechiselparams *p)
   num=gal_statistics_number(qprm.erode_th);
   nval=((size_t *)(num->array))[0];
   if( nval < cp->interpnumngb )
-    error(EXIT_FAILURE, 0, "%zu tile(s) can be used for interpolation of the "
-          "quantile threshold values over the full dataset. This is smaller "
-          "than the requested minimum value of %zu (value to the "
-          "`--interpnumngb' option).\n\n"
-          "There are several ways to address the problem. The best and most "
-          "highly recommended is to use a larger input if possible (when the "
-          "input is a crop from a larger dataset). If that is not the case, "
-          "or it doesn't solve the problem, you need to loosen the "
-          "parameters (and therefore cause scatter in the final result). "
-          "Thus don't loosen them too much. Recall that you can see all the "
-          "option values to Gnuastro's programs by appending `-P' to the "
-          "end of your command.\n\n"
-          "  * Slightly decrease `--tilesize' to have more tiles.\n"
-          "  * Slightly increase `--meanmedqdiff' to accept more tiles.\n"
-          "  * Decrease `--outliersigma' to reject less tiles as outliers.\n"
-          "  * Decrease `--interpnumngb' to be smaller than %zu.\n\n"
-          "Append the previous command with `--checkqthresh' to see the "
-          "successful tiles and which were discarded as outliers. This will "
-          "help you find the cause/solution. Note that the output is a "
-          "multi-extension FITS file).\n\n"
-          "To better understand this important step, please run the "
-          "following command (press `SPACE'/arrow-keys to navigate and "
-          "`Q' to return back to the command-line):\n\n"
-          "    $ info gnuastro \"Quantifying signal in a tile\"\n",
-          nval, cp->interpnumngb, cp->interpnumngb);
+    threshold_good_error(nval, 1, cp->interpnumngb);
 
 
   /* Interpolate and smooth the derived values. */
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 1225654..5af6ba7 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -23741,6 +23741,14 @@ input->flags &= ~GAL_DATA_FLAG_BLANK_CH;
 @end example
 @end deftypefun
 
address@hidden size_t gal_blank_number (gal_data_t @code{*input}, int 
@code{updateflag})
+Return the number of blank elements in @code{input}. If
address@hidden, then the dataset blank keyword flags will be
+updated. See the description of @code{gal_blank_present} (above) for more
+on these flags. If @code{input==NULL}, then this function will return
address@hidden
address@hidden deftypefun
+
 
 @deftypefun {gal_data_t *} gal_blank_flag (gal_data_t @code{*input})
 Create a dataset of the the same size as the input, but with an
@@ -24091,11 +24099,12 @@ necessary allocations. If they aren't @code{NULL}, 
all input arrays
 separately copied (allocated) by this function for usage in @code{data}, so
 you can safely use one value to initialize many datasets or use statically
 allocated variables in this function call. Once you are done with the
-dataset, you can clean all the allocated spaces with
+dataset, you can free all the allocated spaces with
 @code{gal_data_free_contents}.
 
 If @code{array} is not @code{NULL}, it will be directly copied into
address@hidden>array} and no new space will be allocated for the array of this
address@hidden>array} (based on the total number of elements calculated from
address@hidden) and no new space will be allocated for the array of this
 dataset, this has many low-level advantages and can be used to work on
 regions of a dataset instead of the whole allocated array (see the
 description under @code{block} in @ref{Generic data container} for one
diff --git a/lib/blank.c b/lib/blank.c
index cbf5658..4445d5d 100644
--- a/lib/blank.c
+++ b/lib/blank.c
@@ -33,6 +33,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/tile.h>
 #include <gnuastro/blank.h>
 #include <gnuastro/pointer.h>
+#include <gnuastro/statistics.h>
 
 #include <gnuastro-internal/checkset.h>
 
@@ -462,6 +463,33 @@ gal_blank_present(gal_data_t *input, int updateflag)
 
 
 
+/* Return the number of blank elements in the dataset. */
+size_t
+gal_blank_number(gal_data_t *input, int updateflag)
+{
+  gal_data_t *number;
+  size_t num_not_blank;
+
+  if(input)
+    {
+      if( gal_blank_present(input, updateflag) )
+        {
+          number=gal_statistics_number(input);
+          num_not_blank=((size_t *)(number->array))[0];
+          gal_data_free(number);
+          return input->size - num_not_blank;
+        }
+      else
+        return 0;
+    }
+  else
+    return GAL_BLANK_SIZE_T;
+}
+
+
+
+
+
 
 
 /* Create a dataset of the the same size as the input, but with an uint8_t
diff --git a/lib/gnuastro/blank.h b/lib/gnuastro/blank.h
index 3ccc206..2e2490d 100644
--- a/lib/gnuastro/blank.h
+++ b/lib/gnuastro/blank.h
@@ -109,6 +109,9 @@ gal_blank_is(void *pointer, uint8_t type);
 int
 gal_blank_present(gal_data_t *input, int updateflag);
 
+size_t
+gal_blank_number(gal_data_t *input, int updateflag);
+
 gal_data_t *
 gal_blank_flag(gal_data_t *data);
 
diff --git a/lib/statistics.c b/lib/statistics.c
index 1b45453..5ba3c19 100644
--- a/lib/statistics.c
+++ b/lib/statistics.c
@@ -1421,9 +1421,14 @@ gal_statistics_no_blank_sorted(gal_data_t *input, int 
inplace)
       else
         sorted=noblank;
     }
-  /* When the input's size is zero, just return the actual input. */
+
+  /* Input's size was zero. Note that we cannot simply copy the zero-sized
+     input dataset, we'll have to allocate it here. */
   else
-    sorted = inplace ? input : gal_data_copy(input);
+    sorted = ( inplace
+               ? input
+               : gal_data_alloc(NULL, input->type, 0, NULL, input->wcs, 0,
+                                input->minmapsize, NULL, NULL, NULL) );
 
   /* Set the blank and sorted flags if the dataset has zero-elements. Even
      if having blank values or being sorted is not defined on a
@@ -1439,7 +1444,6 @@ gal_statistics_no_blank_sorted(gal_data_t *input, int 
inplace)
       sorted->flag &= ~GAL_DATA_FLAG_SORTED_D;
     }
 
-
   /* Return final array. */
   return sorted;
 }
@@ -2145,6 +2149,9 @@ gal_statistics_outlier_positive(gal_data_t *input, size_t 
window_size,
   /* Remove all blanks and sort the dataset. */
   nbs=gal_statistics_no_blank_sorted(input, inplace);
 
+  /* If all elements are blank, simply return the default (NULL) output. */
+  if(nbs->size==0) return out;
+
   /* Only continue if the window size is more than 2 elements (out
      "outlier" is hard to define on smaller datasets). */
   if(window_size>2)