[gnuastro-commits] master e0512a5: NoiseChisel can ignore blanks in fraction of undetected pixels

[Mohammad Akhlaghi]

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

    NoiseChisel can ignore blanks in fraction of undetected pixels
    
    Until now, the reference for using a tile for estimating the Sky value was
    the number of undetected pixels divided by the total tile area. With this
    commit a new `--skyfracnoblank' option is added. When called, the
    denominator of the fraction will be the number of non-blank elements within
    the tile (not whole tile size). This helps in situtations where the lare
    large blank/masked regions scattered across the image, causing a failure in
    all tiles.
---
 NEWS                   | 10 ++++++++++
 bin/noisechisel/args.h | 13 +++++++++++++
 bin/noisechisel/main.h |  3 ++-
 bin/noisechisel/sky.c  | 12 ++++++++----
 bin/noisechisel/ui.h   |  1 +
 doc/gnuastro.texi      | 40 ++++++++++++++++++++++++++++++++--------
 6 files changed, 66 insertions(+), 13 deletions(-)

diff --git a/NEWS b/NEWS
index 0ada736..8c1a970 100644
--- a/NEWS
+++ b/NEWS
@@ -24,6 +24,16 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
       configured with these options.
        --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.
 
   Statistics:
     - Sky estimation: new outlier estimation algorithm similar to NoiseChisel.
diff --git a/bin/noisechisel/args.h b/bin/noisechisel/args.h
index 56de700..0b88de9 100644
--- a/bin/noisechisel/args.h
+++ b/bin/noisechisel/args.h
@@ -330,6 +330,19 @@ struct argp_option program_options[] =
       GAL_OPTIONS_NOT_SET
     },
     {
+      "skyfracnoblank",
+      UI_KEY_SKYFRACNOBLANK,
+      0,
+      0,
+      "No blanks in tile undetected frac. (minskyfrac).",
+      UI_GROUP_DETECTION,
+      &p->skyfracnoblank,
+      GAL_OPTIONS_NO_ARG_TYPE,
+      GAL_OPTIONS_RANGE_0_OR_1,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
+    },
+    {
       "minskyfrac",
       UI_KEY_MINSKYFRAC,
       "FLT",
diff --git a/bin/noisechisel/main.h b/bin/noisechisel/main.h
index c0d64fd..b0aef31 100644
--- a/bin/noisechisel/main.h
+++ b/bin/noisechisel/main.h
@@ -57,7 +57,6 @@ struct noisechiselparams
   uint8_t               label;  /* Label detections that are connected.   */
 
   float          meanmedqdiff;  /* Difference between mode and median.    */
-  float            minskyfrac;  /* Undetected area min. frac. in tile.    */
   float               qthresh;  /* Quantile threshold on convolved image. */
   float          outliersigma;  /* Multiple of sigma to define outlier.   */
   double      outliersclip[2];  /* Outlier Sigma-clipping params.         */
@@ -68,6 +67,8 @@ struct noisechiselparams
   float          noerodequant;  /* Quantile for no erosion.               */
   size_t              opening;  /* Depth of opening after erosion.        */
   size_t           openingngb;  /* Connectivity to use for opening.       */
+  uint8_t      skyfracnoblank;  /* No blanks in estimating non-det frac.  */
+  float            minskyfrac;  /* Undetected area min. frac. in tile.    */
   double         sigmaclip[2];  /* Sigma-clipping parameters.             */
   uint8_t         checkdetsky;  /* Check pseudo-detection sky value.      */
   float               dthresh;  /* Sigma threshold for Pseudo-detections. */
diff --git a/bin/noisechisel/sky.c b/bin/noisechisel/sky.c
index 582f0bb..1b6e9c6 100644
--- a/bin/noisechisel/sky.c
+++ b/bin/noisechisel/sky.c
@@ -59,8 +59,8 @@ sky_mean_std_undetected(void *in_prm)
   struct noisechiselparams *p=(struct noisechiselparams *)tprm->params;
 
   int setblank, type=GAL_TYPE_FLOAT32;
-  size_t twidth=gal_type_sizeof(GAL_TYPE_FLOAT32);
   size_t i, tind, numsky, bdsize=2, ndim=p->sky->ndim;
+  size_t refarea, twidth=gal_type_sizeof(GAL_TYPE_FLOAT32);
   gal_data_t *tile, *fusage, *busage, *bintile, *sigmaclip;
 
 
@@ -88,20 +88,24 @@ sky_mean_std_undetected(void *in_prm)
       numsky=0;
       tind = tprm->indexs[i];
       tile = &p->cp.tl.tiles[tind];
+      refarea = p->skyfracnoblank ? 0 : tile->size;
 
       /* Correct the fake binary tile's properties to be the same as this
          one, then count the number of zero valued elements in it. Note
          that the `CHECK_BLANK' flag of `GAL_TILE_PARSE_OPERATE' is set to
-         1. So blank values in the input array are not counted also. */
+         1. So blank values in the input array are not counted. */
       bintile->size=tile->size;
       bintile->dsize=tile->dsize;
       bintile->array=gal_tile_block_relative_to_other(tile, p->binary);
-      GAL_TILE_PARSE_OPERATE(tile, bintile, 1, 1, {if(!*o) numsky++;});
+      GAL_TILE_PARSE_OPERATE(tile, bintile, 1, 1, {
+          if(p->skyfracnoblank) ++refarea;
+          if(!*o)               ++numsky;
+        });
 
       /* Only continue, if the fraction of Sky values is less than the
          requested fraction. */
       setblank=0;
-      if( (float)(numsky)/(float)(tile->size) > p->minskyfrac)
+      if( (float)(numsky)/(float)(refarea) > p->minskyfrac)
         {
           /* Re-initialize the usage array's size information (will be
              corrected to this tile's size by
diff --git a/bin/noisechisel/ui.h b/bin/noisechisel/ui.h
index 4c66ef8..bdff519 100644
--- a/bin/noisechisel/ui.h
+++ b/bin/noisechisel/ui.h
@@ -89,6 +89,7 @@ enum option_keys_enum
   UI_KEY_ERODENGB,
   UI_KEY_NOERODEQUANT,
   UI_KEY_OPENINGNGB,
+  UI_KEY_SKYFRACNOBLANK,
   UI_KEY_CHECKDETSKY,
   UI_KEY_CHECKSN,
   UI_KEY_DETGROWMAXHOLESIZE,
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 05e4488..0b23cc1 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -1897,13 +1897,14 @@ should do all the work on a finer pixel grid. In the 
end he can
 re-sample the result to the initially desired grid size.
 
 @item
-Convolve the image with a PSF image that is over-sampled to the same
-value as the mock image. Since he wants to finish in a reasonable time
-and the PSF kernel will be very large due to oversampling, he has to
-use frequency domain convolution which has the side effect of dimming
-the edges of the image. So in the first step above he also has to
-build the image to be larger by at least half the width of the PSF
-convolution kernel on each edge.
address@hidden PSF
+Convolve the image with a point spread function (PSF, see @ref{PSF}) that
+is over-sampled to the same resolution as the mock image. Since he wants to
+finish in a reasonable time and the PSF kernel will be very large due to
+oversampling, he has to use frequency domain convolution which has the side
+effect of dimming the edges of the image. So in the first step above he
+also has to build the image to be larger by at least half the width of the
+PSF convolution kernel on each edge.
 
 @item
 With all the transformations complete, the image should be re-sampled
@@ -1960,7 +1961,7 @@ $ astmkprof -P
 @noindent
 In Gnuastro, column counting starts from 1, so the columns are ordered such
 that the first column (number 1) can be an ID he specifies for each object
-(and MakeProfiles ignores), each subsequent column is used used for another
+(and MakeProfiles ignores), each subsequent column is used for another
 property of the profile. It is also possible to use column names for the
 values of these options and change these defaults, but Sufi preferred to
 stick to the defaults. Fortunately MakeProfiles has the capability to also
@@ -15837,6 +15838,16 @@ a more complete description, see the latter half of 
@ref{Quantifying signal
 in a tile}.
 
 @item
address@hidden: 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
+reference for this fraction was the whole tile size. 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. For more, see the description of this option in @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
@@ -16262,6 +16273,19 @@ them. Once opening is complete, we have @emph{initial} 
detections.
 The structuring element used for opening, see @option{--erodengb} for more
 information about a structuring element.
 
address@hidden --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 (value given to
address@hidden). By default this fraction is found by dividing
+number of undetected pixels in a tile by the tile's area. But this default
+behavior ignores the possibility of blank pixels. In situations that
+blank/masked pixels are scattered across the image and if they are large
+enough, all the tiles can fail the @option{--minskyfrac} test, thus not
+allowing NoiseChisel to proceed. With this option, such scenarios can be
+fixed: the denominator of the fraction will be the number of non-blank
+elements in the tile, not the total tile area.
+
 @item -B FLT
 @itemx --minskyfrac=FLT
 Minimum fraction (value between 0 and 1) of Sky (undetected) areas in a