[gnuastro-commits] master 3479c92: NoiseChisel can now ignore some qthresh tiles

[Mohammad Akhlaghi]

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

    NoiseChisel can now ignore some qthresh tiles
    
    Until now, NoiseChisel would use all the tiles with a successful
    measurement. As a result, if we had a region with a large diffuse signal,
    the large diffuse signal would be treated like an elevated background and
    effectively detected as sky.
    
    From this commit onwards, NoiseChisel has a `--qthreshtilequant' option to
    discard tiles above a given quantile.
---
 NEWS                                |  11 +++
 bin/noisechisel/args.h              |  13 +++
 bin/noisechisel/astnoisechisel.conf |   1 +
 bin/noisechisel/main.h              |   1 +
 bin/noisechisel/threshold.c         | 152 +++++++++++++++++++++++++++++++++++-
 bin/noisechisel/ui.h                |   1 +
 doc/gnuastro.texi                   |  97 ++++++++++++++++++++---
 lib/gnuastro-internal/commonopts.h  |   2 +-
 8 files changed, 263 insertions(+), 15 deletions(-)

diff --git a/NEWS b/NEWS
index d065e6a..a5b711b 100644
--- a/NEWS
+++ b/NEWS
@@ -16,6 +16,12 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
   are done on the dataset convolved with `--kernel' as they were
   before. Since it is time consuming, this is an optional feature.
 
+  NoiseChisel: with the new `--qthreshtilequant' option, it is now possible
+  to discard high-valued (outlier) tiles before estimating qthresh over the
+  whole image. This can be useful in detecting very large diffuse/flat
+  regions that would otherwise be detected as background (and effectively
+  removed).
+
   NoiseChisel: the finally selected true detections are now grown based on
   signal contiguity, not by blind dilation. The growth process is the same
   as the growing of clumps to define objects. Only for true detections, the
@@ -29,6 +35,11 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
 
 ** Removed features
 
+  NoiseChisel: The `--dilate' and `--dilatengb' options have been
+  removed. Growing of true detections is no longer done through dilation
+  but through the `--detgrowquant' and `--detgrowmaxholesize' options (see
+  above).
+
 ** Changed features
 
   `gal_binary_fill_holes' now accepts a `connectivity' and `maxsize'
diff --git a/bin/noisechisel/args.h b/bin/noisechisel/args.h
index 80de2ce..cf9934b 100644
--- a/bin/noisechisel/args.h
+++ b/bin/noisechisel/args.h
@@ -221,6 +221,19 @@ struct argp_option program_options[] =
       GAL_OPTIONS_NOT_SET
     },
     {
+      "qthreshtilequant",
+      UI_KEY_QTHRESHTILEQUANT,
+      "FLT",
+      0,
+      "Remove tiles at higher quantiles.",
+      ARGS_GROUP_DETECTION,
+      &p->qthreshtilequant,
+      GAL_TYPE_FLOAT32,
+      GAL_OPTIONS_RANGE_GE_0_LE_1,
+      GAL_OPTIONS_MANDATORY,
+      GAL_OPTIONS_NOT_SET
+    },
+    {
       "smoothwidth",
       UI_KEY_SMOOTHWIDTH,
       "INT",
diff --git a/bin/noisechisel/astnoisechisel.conf 
b/bin/noisechisel/astnoisechisel.conf
index a037355..b3ab828 100644
--- a/bin/noisechisel/astnoisechisel.conf
+++ b/bin/noisechisel/astnoisechisel.conf
@@ -30,6 +30,7 @@
  mirrordist         1.5
  modmedqdiff       0.01
  qthresh            0.3
+ qthreshtilequant   1.0
  smoothwidth          3
  erode                2
  erodengb             4
diff --git a/bin/noisechisel/main.h b/bin/noisechisel/main.h
index 5d08e59..c994c4f 100644
--- a/bin/noisechisel/main.h
+++ b/bin/noisechisel/main.h
@@ -61,6 +61,7 @@ struct noisechiselparams
   float            mirrordist;  /* Maximum distance to check mode sym.    */
   float           modmedqdiff;  /* Difference between mode and median.    */
   float               qthresh;  /* Quantile threshold on convolved image. */
+  float      qthreshtilequant;  /* Remove tiles with lower quantile.      */
   size_t          smoothwidth;  /* Width of flat kernel to smooth.        */
   uint8_t        checkqthresh;  /* Save the quantile threhsold steps.     */
   size_t                erode;  /* Number of erosions after thresholding. */
diff --git a/bin/noisechisel/threshold.c b/bin/noisechisel/threshold.c
index 34ae8fd..9ab5dd5 100644
--- a/bin/noisechisel/threshold.c
+++ b/bin/noisechisel/threshold.c
@@ -499,6 +499,149 @@ qthresh_on_tile(void *in_prm)
 
 
 
+/* The main working function for `threshold_qthresh_clean'. The main
+   purpose/problem is this: when we have channels, the qthresh values for
+   each channel should be treated independently. */
+static void
+threshold_qthresh_clean_work(struct noisechiselparams *p, gal_data_t *first,
+                             gal_data_t *second, gal_data_t *third,
+                             size_t start, size_t number)
+{
+  gal_data_t *quantile;
+  size_t i, osize=first->size;
+  float *oa1=NULL, *oa2=NULL, *oa3=NULL;
+  float q, *arr1=NULL, *arr2=NULL, *arr3=NULL;
+
+  /* A small sanity check. */
+  if(first->type!=GAL_TYPE_FLOAT32)
+    error(EXIT_FAILURE, 0, "%s: datatype has to be float32", __func__);
+
+  /* Correct the arrays (if necessary). IMPORTANT: The datasets are
+     multi-dimensional. However, when estimating the quantile, their
+     dimensionality doesn't matter (only the `size' element is checked by
+     `gal_statistics_quantile', not `ndim' or `dsize'). So we just need to
+     correct `size' if channels are to be considered. */
+  if(start || number!=first->size)
+    {
+      /* Keep the original values for re-setting later. */
+      oa1=first->array;
+      oa2=second->array;
+      if(third) oa3=third->array;
+
+      /* Increment the array pointers. */
+      first->array=gal_data_ptr_increment(first->array, start, first->type);
+      second->array=gal_data_ptr_increment(second->array, start,
+                                           second->type);
+      if(third)
+        third->array=gal_data_ptr_increment(third->array, start, third->type);
+
+      /* Correct their sizes. */
+      first->size=number;
+      second->size=number;
+      if(third) third->size=number;
+    }
+
+  /* Find the quantile and remove all tiles that are more than it in the
+     first array. */
+  arr1=first->array;
+  quantile=gal_statistics_quantile(first, p->qthreshtilequant, 0);
+  q=*((float *)(quantile->array));
+  for(i=0;i<first->size;++i)
+    /* Just note that we have blank (NaN) values, so to avoid doing a
+       NaN check with `isnan', we will check if the value is below the
+       quantile, if it succeeds (isn't NaN and is below the quantile),
+       then we'll put it's actual value, otherwise, a NaN. */
+    arr1[i] = arr1[i]<q ? arr1[i] : NAN;
+  gal_data_free(quantile);
+
+  /* Second quantile threshold. */
+  arr2=second->array;
+  quantile=gal_statistics_quantile(second, p->qthreshtilequant, 0);
+  q=*((float *)(quantile->array));
+  for(i=0;i<second->size;++i)
+    arr2[i] = arr2[i]<q ? arr2[i] : NAN;
+  gal_data_free(quantile);
+
+  /* The third (if it exists). */
+  if(third)
+    {
+      arr3=third->array;
+      quantile=gal_statistics_quantile(third, p->qthreshtilequant, 0);
+      q=*((float *)(quantile->array));
+      for(i=0;i<third->size;++i)
+        arr3[i] = arr3[i]<q ? arr3[i] : NAN;
+      gal_data_free(quantile);
+    }
+
+  /* Make sure all three have the same NaN pixels. */
+  for(i=0;i<first->size;++i)
+    if( isnan(arr1[i]) || isnan(arr2[i]) || isnan(arr3[i]) )
+      {
+        arr1[i] = arr2[i] = NAN;
+        if(third) arr3[i] = NAN;
+      }
+
+  /* Correct the values if they were changed. */
+  if(start || number!=osize)
+    {
+      first->array=oa1;
+      second->array=oa2;
+      first->size = second->size = osize;
+      if(third) { third->array=oa3; third->size=osize; }
+    }
+}
+
+
+
+
+
+/* Clean higher valued quantile thresholds: useful when the diffuse (almost
+   flat) structures are much larger than the tile size. */
+static void
+threshold_qthresh_clean(struct noisechiselparams *p, gal_data_t *first,
+                        gal_data_t *second, gal_data_t *third,
+                        char *filename)
+{
+  size_t i;
+  struct gal_tile_two_layer_params *tl=&p->cp.tl;
+
+  /* A small sanity check: */
+  if(first->size!=tl->tottiles)
+    error(EXIT_FAILURE, 0, "%s: `first->size' and `tl->tottiles' must have "
+          "the same value, but they don't: %zu, %zu", __func__, first->size,
+          tl->tottiles);
+
+  /* If the input is from a tile structure and the user has asked to ignore
+     channels, then re-order the values. */
+  for(i=0;i<tl->totchannels;++i)
+    threshold_qthresh_clean_work(p, first, second, third,
+                                 i*tl->tottilesinch, tl->tottilesinch);
+
+  /* If the user wants to see the steps. */
+  if(p->qthreshname)
+    {
+      first->name="QTHRESH_ERODE_CLEAN";
+      second->name="QTHRESH_NOERODE_CLEAN";
+      gal_tile_full_values_write(first, tl, 1, p->qthreshname, NULL,
+                                 PROGRAM_NAME);
+      gal_tile_full_values_write(second, tl, 1, p->qthreshname, NULL,
+                                 PROGRAM_NAME);
+      first->name=second->name=NULL;
+
+      if(third)
+        {
+          third->name="QTHRESH_EXPAND_CLEAN";
+          gal_tile_full_values_write(third, tl, 1, p->qthreshname,
+                                     NULL, PROGRAM_NAME);
+          third->name=NULL;
+        }
+    }
+}
+
+
+
+
+
 void
 threshold_quantile_find_apply(struct noisechiselparams *p)
 {
@@ -582,9 +725,16 @@ threshold_quantile_find_apply(struct noisechiselparams *p)
     }
 
 
+  /* Remove higher thresholds if requested. */
+  if(p->qthreshtilequant!=1.0)
+    threshold_qthresh_clean(p, qprm.erode_th, qprm.noerode_th,
+                            qprm.expand_th ? qprm.expand_th : NULL,
+                            p->qthreshname);
+
+
   /* Interpolate and smooth the derived values. */
   threshold_interp_smooth(p, &qprm.erode_th, &qprm.noerode_th,
-                          qprm.expand_th ? & qprm.expand_th : NULL,
+                          qprm.expand_th ? &qprm.expand_th : NULL,
                           p->qthreshname);
 
 
diff --git a/bin/noisechisel/ui.h b/bin/noisechisel/ui.h
index bfe3b22..c8209ae 100644
--- a/bin/noisechisel/ui.h
+++ b/bin/noisechisel/ui.h
@@ -67,6 +67,7 @@ enum option_keys_enum
   UI_KEY_ONLYDETECTION,
   UI_KEY_GROWNCLUMPS,
   UI_KEY_SMOOTHWIDTH,
+  UI_KEY_QTHRESHTILEQUANT,
   UI_KEY_CHECKQTHRESH,
   UI_KEY_ERODENGB,
   UI_KEY_NOERODEQUANT,
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index a567e7f..172b369 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -12225,15 +12225,43 @@ successfully.
 Before using NoiseChisel it is strongly recommended to read
 @url{https://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa [2015]} to
 gain a good understanding of what it does and how each parameter influences
-the output. Thanks to that paper, there is no more need to continue this
-introduction any further and we can just dive into the details of running
-NoiseChisel in @ref{Invoking astnoisechisel}. However, the paper cannot
-undergo any further updates, but NoiseChisel will evolve: better algorithms
-or steps will be found, thus options will be added or removed. So this book
-is the final and definitive guide. To make the transition form the paper to
-this book easier (and encourage reading the paper), below you can see the
-major changes since that paper was published.
+the output. Thanks to that paper, there is no need to go into the details
+of the major processing steps. Hence we can just dive into the details of
+running NoiseChisel in @ref{Invoking astnoisechisel}.
+
+However, the paper cannot undergo any further updates, but NoiseChisel will
+evolve: better algorithms or steps will be found, thus options will be
+added or removed. So this book is the final and definitive guide. For a
+more detailed list of changes in each release, please follow the
address@hidden file. The @file{NEWS} file is in the released Gnuastro tarball
+(see @ref{Release tarball}). You can also view the most recent @file{NEWS}
+file @url{http://git.savannah.gnu.org/cgit/gnuastro.git/plain/NEWS,
address@hidden online version of the @file{NEWS} file here may
+contain features that have been implemented, but not yet officially
+released as a tarball (see @ref{Downloading the source}). Therefore, please
+be sure to look at the dates and versions above each group of changed
+features and make sure it corresponds to your installed version. It is
+hence recommended to always stay up to date, see @ref{Announcements}).}.
+
+To make the transition form the paper to this book easier (and encourage
+reading the paper), below you can see the major changes since the paper was
+published. First, the options that have been removed are discussed,
+followed by those that have been added.
+
address@hidden
+Removed options:
address@hidden
address@hidden
address@hidden: In the paper, true detections were dilated for a final
+dig into the noise. However, simple 8-connected dilation can produce boxy
+results which are not realistic and could miss diffuse flux. The finalq dig
+into the noise is now done by ``grow''ing the true detections, similar to
+how true clumps were grown, see the description of @option{--detgrowquant}
+below and in @ref{Detection options} for more on the new alternative.
address@hidden itemize
 
address@hidden
+Added options:
 @itemize
 
 @item
@@ -12261,6 +12289,11 @@ give this option a value of 1 (only the largest valued 
pixel in the input
 will not be eroded).
 
 @item
address@hidden: to manually remove the measured qthresh from
+some tiles. This feature helps in detecting large and extended diffuse
+(almost flat) signal when necessary, see @ref{Detection options}.
+
address@hidden
 @option{--detgrowquant}: is used to grow the final true detections until a
 given quantile in the same way that clumps are grown during segmentation
 (compare columns 2 and 3 in Figure 10 of the paper). It replaces the old
@@ -12274,6 +12307,11 @@ noise properties and how cleanly it was Sky 
subtracted. The new
 size to fill as part of this growth, see the description in @ref{Detection
 options} for more details.
 
+This new growth process can be much more successful in detecting diffuse
+flux around true detections compared to dilation and give more realistic
+results, but it can also increase the NoiseChisel run time (depending on
+the given value and input size).
+
 @item
 @option{--cleangrowndet}: A process to further clean/remove the possibility
 of false detections, see the descriptions under this option in
@@ -12281,11 +12319,6 @@ of false detections, see the descriptions under this 
option in
 
 @end itemize
 
-For a more detailed list of updates in each release, please follow the
address@hidden file. The @file{NEWS} file is in the released Gnuastro tarball
-(see @ref{Release tarball}). You can also see it online at
address@hidden://git.savannah.gnu.org/cgit/gnuastro.git/plain/NEWS}.
-
 @node Invoking astnoisechisel,  , NoiseChisel changes after publication, 
NoiseChisel
 @subsection Invoking NoiseChisel
 
@@ -12590,6 +12623,44 @@ is complete, we will have a binary (two valued) image. 
The pixels above the
 threshold are known as foreground pixels (have a value of 1) while those
 which lie below the threshold are known as background (have a value of 0).
 
address@hidden --qthreshtilequant=FLT
+Only keep tiles which have a q-thresh value above the given quantile of the
+all successful tiles. Hence, when given a value of @code{1}, this option
+will be ignored. When there is more than one channel (and
address@hidden is not called), quantile calculation and application
+will be done on each channel independently.
+
+This option is useful when a large and diffuse (almost flat within each
+tile) signal exists with very small regions of Sky. The flat-ness of the
+profile will cause it to successfully pass the tests of @ref{Quantifying
+signal in a tile}. As a result, without this option the flat and diffuse
+signal will be interpretted as sky. In such cases, you can see the status
+of the tiles with the @option{--checkqthresh} option (first image extension
+is enough) and select a quantile through this option to ignore the measured
+values in the higher-valued tiles.
+
+This option can also be useful when there are large bright objects in the
+image with large flat wings which can also pass the tests and give outlier
+values. When there is a sky gradient over the image (mainly due to
+post-processing issues like bad flat fielding), this option must be set to
address@hidden so it is completely ignored and the sky gradient is accurately
+measured and subtracted.
+
+To get an estimate of the measured qthresh distribution, you can run the
+following commands. The first will create the qthresh check image with one
+value/pixel per tile (see @ref{Processing options}). Open the image in a
+FITS viewer and inspect it. The second command below will print the basic
+information about the distribution of values and the third will print the
+value at the 0.4 quantile. Recall that Gnuastro's Statistics program
+ignores blank values (in this case: tiles with significant signal), see
address@hidden
+
address@hidden
+$ astnoisechisel image.fits --checkqthresh --oneelempertile
+$ aststatistics image_qthresh.fits
+$ aststatistics image_qthresh.fits --quantile=0.4
address@hidden example
+
 @item --smoothwidth=INT
 Width of flat kernel used to smooth the interpolated quantile thresholds,
 see @option{--qthresh} for more.
diff --git a/lib/gnuastro-internal/commonopts.h 
b/lib/gnuastro-internal/commonopts.h
index 78ad0d1..be69dee 100644
--- a/lib/gnuastro-internal/commonopts.h
+++ b/lib/gnuastro-internal/commonopts.h
@@ -164,7 +164,7 @@ struct argp_option gal_commonopts_options[] =
       GAL_OPTIONS_KEY_ONEELEMPERTILE,
       0,
       0,
-      "Only display 1 element/tile, not full input res.",
+      "Display 1 element/tile, not full input res.",
       GAL_OPTIONS_GROUP_TESSELLATION,
       &cp->tl.oneelempertile,
       GAL_OPTIONS_NO_ARG_TYPE,