[gnuastro-commits] master eaab08ce: Library (arithmetic.h and dimension.h): maskfilled operators max area

[Mohammad Akhlaghi]

branch: master
commit eaab08cee634361d9b5b0e94e5ab51d296ee1a97
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    Library (arithmetic.h and dimension.h): maskfilled operators max area
    
    Until now, when we were using the newly added maskfilled operators and one
    of the input datasets was fully bad (for example a much higher noise level
    than the rest), it because of filling holes, the outer edges of the bad
    dataset would remain, while most of the inner region would be masked. This
    would make the thin outer edges have different statistics than the main
    area of the final stack!
    
    With this commit, we now check the area covered by the mask and if it
    covers more than 95% of the image's size, then the whole image is masked.
---
 doc/gnuastro.texi | 12 +++++++++++-
 lib/arithmetic.c  | 14 +++++++++++---
 lib/dimension.c   | 10 +++++++++-
 3 files changed, 31 insertions(+), 5 deletions(-)

diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 72fccb1d..ea589c2b 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -22016,7 +22016,12 @@ $ astarithmetic image.fits medianvalue -q
 @item sigclip-maskfilled
 Mask (set to blank/NaN) all the outlying elements (defined by @mymath{\sigma} 
or MAD clipping) in the inputs and put all the inputs back on the stack.
 The first popped operand is the termination criteria of the clipping, the 
second popped operand is the multiple of @mymath{\sigma} or MAD and the third 
is the number of input datasets that will be popped for the actual operation.
-If you are not yet familiar with these @mymath{\sigma} or MAD clipping, it is 
recommended to read this tutorial: @ref{Clipping outliers}.
+If you are not yet familiar with @mymath{\sigma} or MAD clipping, it is 
recommended to read this tutorial: @ref{Clipping outliers}.
+
+When more than 95@mymath{\%} of the area of an operand is masked, the full 
operand will be masked.
+This is necesasry in like this: one of your inputs has many outliers (for 
example it is much more noisy than the rest or its sky level has not been 
subtracted properly).
+Because this operator fills holes between outlying pixels, most of the area of 
the input will be masked, but the thin edges (where there are no ``holes'') 
will remain, causing different statistics in those thin edges of that input in 
your final stack.
+Through this mask coverage fraction (which is currently 
hard-coded@footnote{Please get in touch with us at @code{bug-gnuastro@@gnu.org} 
if you notice this problem and feel the fraction needs to be lowered (or 
generally to be set in each run).}), we ensure that such thin edges do not 
cause artifacts in the final stack.
 
 For example, with the second command below, we are masking the MAD clipped 
pixels of the 9 inputs (that are generated in @ref{Clipping outliers}) and 
writing them as separate HDUs of the output.
 The clipping is done with 5 times the MAD and the clipping starts when the 
relative difference between subsequent MADs is 0.01.
@@ -22722,6 +22727,11 @@ This is the most robust method to reject outliers; for 
more on filled re-clippin
 For a more general tutorial on rejecting outliers, see @ref{Clipping outliers}.
 If you have not done this tutorial yet, we recommend you to take an hour or so 
and go through that tutorial for optimal understanding and results.
 
+When more than 95@mymath{\%} of the area of an operand is masked, the full 
operand will be masked.
+This is necesasry in like this: one of your inputs has many outliers (for 
example it is much more noisy than the rest or its sky level has not been 
subtracted properly).
+Because this operator fills holes between outlying pixels, most of the area of 
the input will be masked, but the thin edges (where there are no ``holes'') 
will remain, causing different statistics in those thin edges of that input in 
your final stack.
+Through this mask coverage fraction (which is currently 
hard-coded@footnote{Please get in touch with us at @code{bug-gnuastro@@gnu.org} 
if you notice this problem and feel the fraction needs to be lowered (or 
generally to be set in each run).}), we ensure that such thin edges do not 
cause artifacts in the final stack.
+
 For example, with the command below, the pixels of the input 2 dimensional 
@file{image.fits} will be collapsed to a single dimension output.
 The first popped operand is @code{2}, so it will collapse all the pixels that 
are vertically on top of each other.
 Such that the output will have the same number of pixels as the horizontal 
axis of the input.
diff --git a/lib/arithmetic.c b/lib/arithmetic.c
index b07dc049..a7907a93 100644
--- a/lib/arithmetic.c
+++ b/lib/arithmetic.c
@@ -2509,12 +2509,13 @@ arithmetic_multioperand_clip_mask_worker(void *in_prm)
     (struct arithmetic_multioperand_clip_mask_params *)tprm->params;
 
   /* Subsequent definitions. */
+  uint8_t *u, *uf;
   gal_data_t *use, *tmp;
-  size_t i, index, counter, ndim=p->list->ndim;
+  size_t i, index, counter, sumflag, ndim=p->list->ndim;
   int aflags=GAL_ARITHMETIC_FLAG_NUMOK; /* Don't free the inputs. */
 
-  /* Go over all the actions (pixels in this case) that were assigned to
-     this thread. */
+  /* Go over all the actions (input datasets in this case) that were
+     assigned to this thread. */
   for(i=0; tprm->indexs[i] != GAL_BLANK_SIZE_T; ++i)
     {
       /* For easy reading. */
@@ -2542,6 +2543,13 @@ arithmetic_multioperand_clip_mask_worker(void *in_prm)
       tmp=gal_binary_erode(tmp, 2, ndim, 1);
       tmp=gal_binary_dilate(tmp, 2, ndim, 1);
 
+      /* If the covered area is larger than a certain fraction of the
+         area, flag the whole dataset. */
+      sumflag=0;
+      uf=(u=tmp->array)+tmp->size; do sumflag+=*u; while(++u<uf);
+      if(sumflag>0.95*tmp->size)
+        {uf=(u=tmp->array)+tmp->size; do *u=1; while(++u<uf);}
+
       /* Set all the 1-valued pixels in the binary image to NaN in the
          input. */
       gal_blank_flag_apply(use, tmp);
diff --git a/lib/dimension.c b/lib/dimension.c
index 46898146..1eb36862 100644
--- a/lib/dimension.c
+++ b/lib/dimension.c
@@ -989,8 +989,9 @@ dimension_collapse_sortbased_fill(struct 
dimension_sortbased_p *p,
                                   gal_data_t *conv, uint8_t clipflags,
                                   size_t wdsize, int check)
 {
-  size_t one=1;
   int std1_mad0=0;
+  uint8_t *u, *uf;
+  size_t one=1, sumflag;
   float *farr=fstat->array;
   gal_data_t *formask=conv?conv:work, *out=NULL;
   gal_data_t *tmp, *multip, *upper, *lower, *center, *spread;
@@ -1076,6 +1077,13 @@ dimension_collapse_sortbased_fill(struct 
dimension_sortbased_p *p,
   tmp=gal_binary_erode(tmp, 2, formask->ndim, 1);
   tmp=gal_binary_dilate(tmp, formask->ndim==1?4:2, formask->ndim, 1);
 
+  /* If the covered area is larger than a certain fraction of the
+     area, flag the whole dataset. */
+  sumflag=0;
+  uf=(u=tmp->array)+tmp->size; do sumflag+=*u; while(++u<uf);
+  if(sumflag>0.95*tmp->size)
+    {uf=(u=tmp->array)+tmp->size; do *u=1; while(++u<uf);}
+
   /* Apply the flag onto the input (to set the pixels to NaN). */
   gal_blank_flag_apply(work, tmp);