[gnuastro-commits] master f4741ec 2/2: NoiseChisel ignoring tiles with zero standard deviation

[Mohammad Akhlaghi]

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

    NoiseChisel ignoring tiles with zero standard deviation
    
    It might happen that the input dataset has zero values for regions with no
    data (essentially using zero as a mask value). In other detection methods
    (that have above-Sky thresholds) this isn't an issue. But NoiseChisel's
    threshold is below the Sky value, so the zero valued pixels will actually
    be considered as detections. We therefore encourge users to set zero valued
    pixels to NaN/blank.
    
    When there are tiles with fully zero-valued pixels, they may thus be used
    for estimating the pseudo-detection Sky and Sky STD and result in
    zero-valued Sky (and its STD) values. This will later cause the
    counts-per-second correction (`cpscorr') to be zero and produce `inf' S/N
    values.
    
    With this commit, a check is added after finding the Sky and its Standard
    deviation on a dataset. If the Standard deviation is zero, the tile will be
    given a blank value (to be interpolated later).
    
    This issue was found based on a dataset provided by Brandon Kelly.
    
    This fixes bug #53825.
---
 NEWS                         |  1 +
 THANKS                       |  1 +
 bin/noisechisel/detection.c  |  1 -
 bin/noisechisel/sky.c        | 47 ++++++++++++++++++++++++++++++--------------
 doc/announce-acknowledge.txt |  1 +
 5 files changed, 35 insertions(+), 16 deletions(-)

diff --git a/NEWS b/NEWS
index 560f7af..a870f12 100644
--- a/NEWS
+++ b/NEWS
@@ -198,6 +198,7 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
   bug #53407: Instrumental noise in MakeNoise should be squared.
   bug #53424: Sigma-clipping seg-faults when there are no more elements.
   bug #53580: Warp crash when no WCS present.
+  bug #53825: NoiseChisel not accounting for fully zero-valued tiles.
 
 
 
diff --git a/THANKS b/THANKS
index 5103035..22fafa8 100644
--- a/THANKS
+++ b/THANKS
@@ -34,6 +34,7 @@ support in Gnuastro. The list is ordered alphabetically (by 
family name).
     Brandon Invergo                      address@hidden
     Aurélien Jarno                       address@hidden
     Lee Kelvin                           address@hidden
+    Brandon Kelly                        address@hidden
     Mohammad-Reza Khellat                address@hidden
     Floriane Leclercq                    address@hidden
     Alan Lefor                           address@hidden
diff --git a/bin/noisechisel/detection.c b/bin/noisechisel/detection.c
index 095b450..26ed03a 100644
--- a/bin/noisechisel/detection.c
+++ b/bin/noisechisel/detection.c
@@ -612,7 +612,6 @@ detection_sn(struct noisechiselparams *p, gal_data_t 
*worklab, size_t num,
           }
     }
 
-
   /* If we are in Sky mode, the sizes have to be corrected */
   if(s0d1D2==0)
     {
diff --git a/bin/noisechisel/sky.c b/bin/noisechisel/sky.c
index 815bc0f..2d003c9 100644
--- a/bin/noisechisel/sky.c
+++ b/bin/noisechisel/sky.c
@@ -57,7 +57,7 @@ sky_mean_std_undetected(void *in_prm)
   struct noisechiselparams *p=(struct noisechiselparams *)tprm->params;
 
   double *darr, s, s2;
-  int type=p->sky->type;
+  int setblank, type=p->sky->type;
   size_t i, tind, numsky, dsize=2;
   gal_data_t *tile, *meanstd_d, *meanstd, *bintile;
 
@@ -94,6 +94,7 @@ sky_mean_std_undetected(void *in_prm)
 
       /* Only continue, if the fraction of Sky values are less than the
          requested fraction. */
+      setblank=0;
       if( (float)(numsky)/(float)(tile->size) > p->minskyfrac)
         {
           /* Calculate the mean and STD over this tile. */
@@ -108,22 +109,37 @@ sky_mean_std_undetected(void *in_prm)
           darr[0]=s/numsky;
           darr[1]=sqrt( (s2-s*s/numsky)/numsky );
 
-          /* Convert the mean and std into the same type as the sky and std
-             arrays. */
-          meanstd=gal_data_copy_to_new_type(meanstd_d, type);
-
-          /* Copy the mean and STD to their respective places in the tile
-             arrays. */
-          memcpy(gal_pointer_increment(p->sky->array, tind, type),
-                 meanstd->array, gal_type_sizeof(type));
-          memcpy(gal_pointer_increment(p->std->array, tind, type),
-                 gal_pointer_increment(meanstd->array, 1, type),
-                 gal_type_sizeof(type));
-
-          /* Clean up. */
-          gal_data_free(meanstd);
+          /* When there are zero-valued pixels on the edges of the dataset
+             (that have not been set to NaN/blank), given special
+             conditions, the whole zero-valued region can get a binary
+             value of 1 and so the Sky and its standard deviation can
+             become zero. So, we need ignore such tiles. */
+          if(darr[1]==0.0f)
+            setblank=1;
+          else
+            {
+              /* Convert the mean and std into the same type as the sky and
+                 std arrays. */
+              meanstd=gal_data_copy_to_new_type(meanstd_d, type);
+
+              /* Copy the mean and STD to their respective places in the
+                 tile arrays. */
+              memcpy(gal_pointer_increment(p->sky->array, tind, type),
+                     meanstd->array, gal_type_sizeof(type));
+              memcpy(gal_pointer_increment(p->std->array, tind, type),
+                     gal_pointer_increment(meanstd->array, 1, type),
+                     gal_type_sizeof(type));
+
+              /* Clean up. */
+              gal_data_free(meanstd);
+            }
         }
       else
+        setblank=1;
+
+      /* If the tile is marked as being blank, write blank values into
+         it. */
+      if(setblank==1)
         {
           gal_blank_write(gal_pointer_increment(p->sky->array, tind, type),
                           type);
@@ -204,6 +220,7 @@ sky_and_std(struct noisechiselparams *p, char *checkname)
      factor here. */
   p->cpscorr = p->minstd>1 ? 1.0f : p->minstd;
 
+
   /* Interpolate and smooth the derived values. */
   threshold_interp_smooth(p, &p->sky, &p->std, NULL, checkname);
 
diff --git a/doc/announce-acknowledge.txt b/doc/announce-acknowledge.txt
index c962921..c74b8fa 100644
--- a/doc/announce-acknowledge.txt
+++ b/doc/announce-acknowledge.txt
@@ -4,6 +4,7 @@ Leindert Boogaard
 Nima Dehdilani
 Antonio Diaz Diaz
 Lee Kelvin
+Brandon Kelly
 Guillaume Mahler
 Bertrand Pain
 Ole Streicher