[gnuastro-commits] master c58950c 2/2: New NoiseChisel option to consider blank pixels as foreground

[Mohammad Akhlaghi]

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

    New NoiseChisel option to consider blank pixels as foreground
    
    Until now blank pixels were treated as foreground in NoiseChisel's binary
    processing steps. This would cause many false positives when there were
    lots of blank pixels in the dataset. With this commit the default behavior
    is to consider blank pixels as background and there is a new option
    (`--blankasforeground') if someone wants to have the old behavior.
---
 NEWS                        | 28 +++++++++++++++++++---------
 bin/noisechisel/args.h      | 13 +++++++++++++
 bin/noisechisel/detection.c | 43 +++++++++++++++++++++++++++++++++++++++++++
 bin/noisechisel/main.h      |  1 +
 bin/noisechisel/ui.h        |  1 +
 doc/gnuastro.texi           | 22 ++++++++++++++++++++++
 6 files changed, 99 insertions(+), 9 deletions(-)

diff --git a/NEWS b/NEWS
index 8c1a970..e26fdfb 100644
--- a/NEWS
+++ b/NEWS
@@ -13,6 +13,7 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
     --onedasimage: write output as an image if it has one dimension, not table.
 
   NoiseChisel:
+
     - New outlier identification algorithm for quantile thresholds. This is
       very useful when there are extended and bright sources in the
       dataset: the tiles containing very faint signal that pass the general
@@ -25,15 +26,24 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
        --outliersclip: Sigma-clipping parameters for the process.
        --outliersigma: Multiple of sigma to define an outlier.
     --skyfracnoblank: Ignore blank pixels when estimating the fraction of
-      undetected pixels for Sky estimation. NoiseChisel only measures the
-      Sky over the tiles that have a sufficiently large fraction of
-      undetected pixels. This is done to decrease the bias caused by faint
-      un-detected wings of bright galaxies or stars, see description of
-      `--minskyfrac' for more. Until now the reference for this fraction
-      was the whole tile size (irrespective of how many blank elements it
-      contains). With this option, it is now possible to ask for ignoring
-      blank pixels when calculating the fraction. This is useful when
-      blank/masked pixels are distributed across the image.
+          undetected pixels for Sky estimation. NoiseChisel only measures
+          the Sky over the tiles that have a sufficiently large fraction of
+          undetected pixels. This is done to decrease the bias caused by
+          faint un-detected wings of bright galaxies or stars, see
+          description of `--minskyfrac' for more. Until now the reference
+          for this fraction was the whole tile size (irrespective of how
+          many blank elements it contains). With this option, it is now
+          possible to ask for ignoring blank pixels when calculating the
+          fraction. This is useful when blank/masked pixels are distributed
+          across the image.
+    --blankasforeground: Treat blank elements as foreground (regions above
+          the threshold) in the binary processing steps: initial erosion
+          and opening as well as the filling holes and opening of
+          pseudo-detections. From this version, by default, blank elements
+          in the dataset are considered to be background, so if a
+          foreground pixel is touching it, it will be eroded. This option
+          is irrelevant if the datasets contains no blank elements, but can
+          help remove false positives that are touching blank elements.
 
   Statistics:
     - Sky estimation: new outlier estimation algorithm similar to NoiseChisel.
diff --git a/bin/noisechisel/args.h b/bin/noisechisel/args.h
index 0b88de9..b17acaf 100644
--- a/bin/noisechisel/args.h
+++ b/bin/noisechisel/args.h
@@ -265,6 +265,19 @@ struct argp_option program_options[] =
       GAL_OPTIONS_NOT_SET
     },
     {
+      "blankasforeground",
+      UI_KEY_BLANKASFOREGROUND,
+      0,
+      0,
+      "Blanks are foreground in erosion and opening.",
+      UI_GROUP_DETECTION,
+      &p->blankasforeground,
+      GAL_OPTIONS_NO_ARG_TYPE,
+      GAL_OPTIONS_RANGE_0_OR_1,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
+    },
+    {
       "erode",
       UI_KEY_ERODE,
       "INT",
diff --git a/bin/noisechisel/detection.c b/bin/noisechisel/detection.c
index 194923d..9ad3bca 100644
--- a/bin/noisechisel/detection.c
+++ b/bin/noisechisel/detection.c
@@ -30,6 +30,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 #include <gnuastro/fits.h>
 #include <gnuastro/label.h>
+#include <gnuastro/blank.h>
 #include <gnuastro/binary.h>
 #include <gnuastro/threads.h>
 #include <gnuastro/pointer.h>
@@ -55,8 +56,10 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 void
 detection_initial(struct noisechiselparams *p)
 {
+  float *f;
   char *msg;
   uint8_t *b, *bf;
+  int resetblank=0;
   struct timeval t0, t1;
 
 
@@ -78,6 +81,15 @@ detection_initial(struct noisechiselparams *p)
     }
 
 
+  /* Remove any blank elements from the binary image if requested. */
+  if(p->blankasforeground==0 && gal_blank_present(p->binary,0))
+    {
+      resetblank=1;
+      bf=(b=p->binary->array)+p->binary->size;
+      do *b = *b==GAL_BLANK_UINT8 ? 0 : *b; while(++b<bf);
+    }
+
+
   /* Erode the image. */
   if(!p->cp.quiet) gettimeofday(&t1, NULL);
   gal_binary_erode(p->binary, p->erode, p->erodengb==4 ? 1 : 2, 1);
@@ -115,6 +127,15 @@ detection_initial(struct noisechiselparams *p)
     }
 
 
+  /* Reset the blank values (if requested). */
+  if(resetblank)
+    {
+      f=p->input->array;
+      bf=(b=p->binary->array)+p->binary->size;
+      do *b = isnan(*f++) ? GAL_BLANK_UINT8 : *b; while(++b<bf);
+    }
+
+
   /* Label the connected components. */
   p->numinitialdets=gal_binary_connected_components(p->binary, &p->olabel,
                                                     p->binary->ndim);
@@ -218,6 +239,7 @@ detection_fill_holes_open(void *in_prm)
   struct noisechiselparams *p=fho_prm->p;
 
   void *tarray;
+  uint8_t *b, *bf;
   gal_data_t *tile, *copy, *tblock;
   size_t i, dsize[]={1,1,1,1,1,1,1,1,1,1}; /* For upto 10-Dimensions! */
 
@@ -245,9 +267,20 @@ detection_fill_holes_open(void *in_prm)
       /* Copy the tile into the contiguous patch of memory to work on, but
          first reset the size element so `gal_data_copy_to_allocated' knows
          there is enough space. */
+      copy->flag=0;
       copy->size=p->maxltcontig;
       gal_data_copy_to_allocated(tile, copy);
 
+      /* Take blank values to the background (set them to zero) if
+         necsesary. */
+      if( p->blankasforeground==0
+          && gal_blank_present(p->input,0)
+          && gal_blank_present(copy, 1) )
+        {
+          bf=(b=copy->array)+copy->size;
+          do *b = *b==GAL_BLANK_UINT8 ? 0 : *b; while(++b<bf);
+        }
+
       /* Fill the holes in this tile: holes with maximal connectivity means
          that they are most strongly bounded. */
       gal_binary_holes_fill(copy, copy->ndim, -1);
@@ -291,6 +324,8 @@ static size_t
 detection_pseudo_find(struct noisechiselparams *p, gal_data_t *workbin,
                       gal_data_t *worklab, int s0d1)
 {
+  float *f;
+  uint8_t *b, *bf;
   gal_data_t *bin;
   struct fho_params fho_prm={0, NULL, workbin, worklab, p};
 
@@ -339,6 +374,14 @@ detection_pseudo_find(struct noisechiselparams *p, 
gal_data_t *workbin,
           gal_threads_spin_off(detection_fill_holes_open, &fho_prm,
                                p->ltl.tottiles, p->cp.numthreads);
 
+          /* Reset the blank values (if they were changed). */
+          if( p->blankasforeground==0 && gal_blank_present(p->input,0) )
+            {
+              f=p->input->array;
+              bf=(b=workbin->array)+workbin->size;
+              do *b = isnan(*f++) ? GAL_BLANK_UINT8 : *b; while(++b<bf);
+            }
+
           /* Set the extension name based on the step. */
           switch(fho_prm.step)
             {
diff --git a/bin/noisechisel/main.h b/bin/noisechisel/main.h
index b0aef31..bfe485c 100644
--- a/bin/noisechisel/main.h
+++ b/bin/noisechisel/main.h
@@ -62,6 +62,7 @@ struct noisechiselparams
   double      outliersclip[2];  /* Outlier Sigma-clipping params.         */
   size_t          smoothwidth;  /* Interpolation: flat kernel to smooth.  */
   uint8_t        checkqthresh;  /* Save the quantile threhsold steps.     */
+  uint8_t   blankasforeground;  /* Blank as foreg. in erosion and opening.*/
   size_t                erode;  /* Number of erosions after thresholding. */
   size_t             erodengb;  /* Connectivity for erosion.              */
   float          noerodequant;  /* Quantile for no erosion.               */
diff --git a/bin/noisechisel/ui.h b/bin/noisechisel/ui.h
index bdff519..30593bc 100644
--- a/bin/noisechisel/ui.h
+++ b/bin/noisechisel/ui.h
@@ -86,6 +86,7 @@ enum option_keys_enum
   UI_KEY_OUTLIERSIGMA,
   UI_KEY_OUTLIERSCLIP,
   UI_KEY_CHECKQTHRESH,
+  UI_KEY_BLANKASFOREGROUND,
   UI_KEY_ERODENGB,
   UI_KEY_NOERODEQUANT,
   UI_KEY_OPENINGNGB,
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 0b23cc1..d9a3d17 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -15838,6 +15838,16 @@ a more complete description, see the latter half of 
@ref{Quantifying signal
 in a tile}.
 
 @item
address@hidden: allows blank pixels to be treated as
+foreground in NoiseChisel's binary operations: the initial erosion
+(@option{--erode}) and opening (@option{--open}) as well as the filling
+holes and opening step for defining pseudo-detections
+(@option{--dthresh}). In the published paper, blank pixels were treated as
+foreground by default. To avoid too many false positive near blank/masked
+regions, blank pixels are now considered to be in the background. This
+option will create the old behavior.
+
address@hidden
 @option{--skyfracnoblank}: To reduce the bias caused by undetected wings of
 galaxies and stars in the Sky measurements, NoiseChisel only uses tiles
 that have a sufficiently large fraction of undetected pixels. Until now the
@@ -16211,6 +16221,18 @@ output will have the same pixel size as the input, but 
with the
 @option{--oneelempertile} option, only one pixel will be used for each tile
 (see @ref{Processing options}).
 
address@hidden --blankasforeground
+In the erosion and opening steps below, treat blank elements as foreground
+(regions above the threshold). By default, blank elements in the dataset
+are considered to be background, so if a foreground pixel is touching it,
+it will be eroded. This option is irrelevant if the datasets contains no
+blank elements.
+
+When there are many blank elements in the dataset, treating them as
+foreground will systematically erode their regions less, therefore
+systematically creating more false positives. So use this option (when
+blank values are present) with care.
+
 @item -e INT
 @itemx --erode=INT
 @cindex Erosion