[gnuastro-commits] master bbe151f2: Arithmetic: new MAD-clipping filter operators added

[Mohammad Akhlaghi]

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

    Arithmetic: new MAD-clipping filter operators added
    
    Until now, for robust filtering we only had the sigma-clipping
    solution. However, recently MAD-clipping has been added in Gnuastro which
    is more robust to outliers.
    
    With this commit, using those new libraries, two new operators have been
    added to the Arithmetic program for this.
    
    These operators have been added after a nice discussion with Paola Dimauro
    and Mathias Urbano.
---
 NEWS                         |  5 +++++
 THANKS                       |  1 +
 bin/arithmetic/arithmetic.c  | 33 ++++++++++++++++++++++++++++-----
 bin/arithmetic/arithmetic.h  |  2 ++
 doc/announce-acknowledge.txt |  2 ++
 doc/gnuastro.texi            | 13 +++++++++++++
 6 files changed, 51 insertions(+), 5 deletions(-)

diff --git a/NEWS b/NEWS
index bdd7f864..87566844 100644
--- a/NEWS
+++ b/NEWS
@@ -14,6 +14,11 @@ See the end of the file for license conditions.
 
   - New operators:
 
+    - filter-madclip-mean: filter/smooth the input using MAD-clipped mean.
+
+    - filter-madclip-median: filter/smooth the input using MAD-clipped
+      median.
+
     - madclip-maskfilled: mask (set to NaN) all input elements that are
       outliers (defined by MAD-clipping). Combined with the stacking
       operators this allows removing large contiugous outliers in your
diff --git a/THANKS b/THANKS
index a5d7ae79..e318febe 100644
--- a/THANKS
+++ b/THANKS
@@ -99,6 +99,7 @@ support in Gnuastro. The list is ordered alphabetically (by 
family name).
     Markus Schaney                       markus@riseup.net
     Martin Guerrero Roncel               mar@iaa.es
     Martin Kuemmel                       mkuemmel@usm.lmu.de
+    Mathias Urbano                       mathias.urbano@astro.unistra.fr
     Michael H.F. Wilkinson               m.h.f.wilkinson@gmail.com
     Michael Stein                        mstein@astro.rub.de
     Michel Tallon                        mtallon@obs.univ-lyon1.fr
diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index 794cfa0a..973e1564 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -288,12 +288,22 @@ arithmetic_filter(void *in_prm)
 
 
         case ARITHMETIC_OP_FILTER_SIGCLIP_MEAN:
+        case ARITHMETIC_OP_FILTER_MADCLIP_MEAN:
         case ARITHMETIC_OP_FILTER_SIGCLIP_MEDIAN:
+        case ARITHMETIC_OP_FILTER_MADCLIP_MEDIAN:
           /* The median is always available with a sigma-clip, but the mean
              needs to be explicitly requested. */
           if(afp->operator == ARITHMETIC_OP_FILTER_SIGCLIP_MEAN)
             clipflags = GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_MEAN;
-          sigclip=gal_statistics_clip_sigma(tile, afp->sclip_multip,
+
+          /* Do the main operation. */
+          if(    afp->operator==ARITHMETIC_OP_FILTER_SIGCLIP_MEAN
+              || afp->operator==ARITHMETIC_OP_FILTER_SIGCLIP_MEDIAN )
+            sigclip=gal_statistics_clip_sigma(tile, afp->sclip_multip,
+                                              afp->sclip_param, clipflags,
+                                              0, 1);
+          else
+            sigclip=gal_statistics_clip_mad(tile, afp->sclip_multip,
                                             afp->sclip_param, clipflags,
                                             0, 1);
 
@@ -301,8 +311,10 @@ arithmetic_filter(void *in_prm)
           switch(afp->operator)
             {
             case ARITHMETIC_OP_FILTER_SIGCLIP_MEAN:
+            case ARITHMETIC_OP_FILTER_MADCLIP_MEAN:
               sind = GAL_STATISTICS_CLIP_OUTCOL_MEAN; break;
             case ARITHMETIC_OP_FILTER_SIGCLIP_MEDIAN:
+            case ARITHMETIC_OP_FILTER_MADCLIP_MEDIAN:
               sind = GAL_STATISTICS_CLIP_OUTCOL_MEDIAN; break;
             default:
               error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at "
@@ -368,8 +380,10 @@ wrapper_for_filter(struct arithmeticparams *p, char 
*token, int operator)
   size_t fsize[ARITHMETIC_FILTER_DIM];
   gal_data_t *tmp, *tmp2, *zero, *comp, *params_list=NULL;
   size_t hnfsize[ARITHMETIC_FILTER_DIM], hpfsize[ARITHMETIC_FILTER_DIM];
-  int issigclip=(operator==ARITHMETIC_OP_FILTER_SIGCLIP_MEAN
-                 || operator==ARITHMETIC_OP_FILTER_SIGCLIP_MEDIAN);
+  int isclip=(   operator==ARITHMETIC_OP_FILTER_SIGCLIP_MEAN
+              || operator==ARITHMETIC_OP_FILTER_MADCLIP_MEAN
+              || operator==ARITHMETIC_OP_FILTER_SIGCLIP_MEDIAN
+              || operator==ARITHMETIC_OP_FILTER_MADCLIP_MEDIAN );
 
 
   /* Get the input's number of dimensions. */
@@ -395,7 +409,7 @@ wrapper_for_filter(struct arithmeticparams *p, char *token, 
int operator)
 
   /* Based on the first popped operand's dimensions and the operator, of
      pop the necessary number of operands. */
-  nparams = ndim + (issigclip ? 2 : 0 );
+  nparams = ndim + (isclip ? 2 : 0 );
   for(i=0;i<nparams;++i)
     {
       /* Add this to the list of parameters. */
@@ -412,7 +426,7 @@ wrapper_for_filter(struct arithmeticparams *p, char *token, 
int operator)
 
   /* If this is a sigma-clipping filter, the top two operands are the
      sigma-clipping parameters. */
-  if(issigclip)
+  if(isclip)
     {
       /* Read the sigma-clipping multiple (first element in the list). */
       tmp=gal_list_data_pop(&params_list);
@@ -497,6 +511,7 @@ wrapper_for_filter(struct arithmeticparams *p, char *token, 
int operator)
             { hnfsize[i]=fsize[i]/2; hpfsize[i]=fsize[i]/2-1; }
         }
 
+
       /* For a test.
       printf("fsize: %zu, %zu\n", fsize[0], fsize[1]);
       printf("hnfsize: %zu, %zu\n", hnfsize[0], hnfsize[1]);
@@ -509,11 +524,13 @@ wrapper_for_filter(struct arithmeticparams *p, char 
*token, int operator)
         {
         case ARITHMETIC_OP_FILTER_MEDIAN:
         case ARITHMETIC_OP_FILTER_SIGCLIP_MEDIAN:
+        case ARITHMETIC_OP_FILTER_MADCLIP_MEDIAN:
           type=afp.input->type;
           break;
 
         case ARITHMETIC_OP_FILTER_MEAN:
         case ARITHMETIC_OP_FILTER_SIGCLIP_MEAN:
+        case ARITHMETIC_OP_FILTER_MADCLIP_MEAN:
           type=GAL_TYPE_FLOAT64;
           break;
 
@@ -1457,8 +1474,12 @@ arithmetic_set_operator(char *string, size_t 
*num_operands, int *inlib)
         { op=ARITHMETIC_OP_FILTER_MEDIAN;         *num_operands=0; }
       else if (!strcmp(string, "filter-sigclip-mean"))
         { op=ARITHMETIC_OP_FILTER_SIGCLIP_MEAN;   *num_operands=0; }
+      else if (!strcmp(string, "filter-madclip-mean"))
+        { op=ARITHMETIC_OP_FILTER_MADCLIP_MEAN;   *num_operands=0; }
       else if (!strcmp(string, "filter-sigclip-median"))
         { op=ARITHMETIC_OP_FILTER_SIGCLIP_MEDIAN; *num_operands=0; }
+      else if (!strcmp(string, "filter-madclip-median"))
+        { op=ARITHMETIC_OP_FILTER_MADCLIP_MEDIAN; *num_operands=0; }
       else if (!strcmp(string, "erode"))
         { op=ARITHMETIC_OP_ERODE;                 *num_operands=0; }
       else if (!strcmp(string, "dilate"))
@@ -1733,7 +1754,9 @@ arithmetic_operator_run(struct arithmeticparams *p, int 
operator,
         case ARITHMETIC_OP_FILTER_MEAN:
         case ARITHMETIC_OP_FILTER_MEDIAN:
         case ARITHMETIC_OP_FILTER_SIGCLIP_MEAN:
+        case ARITHMETIC_OP_FILTER_MADCLIP_MEAN:
         case ARITHMETIC_OP_FILTER_SIGCLIP_MEDIAN:
+        case ARITHMETIC_OP_FILTER_MADCLIP_MEDIAN:
           wrapper_for_filter(p, operator_string, operator);
           break;
 
diff --git a/bin/arithmetic/arithmetic.h b/bin/arithmetic/arithmetic.h
index 96cfbfb9..943e29f4 100644
--- a/bin/arithmetic/arithmetic.h
+++ b/bin/arithmetic/arithmetic.h
@@ -35,6 +35,8 @@ enum arithmetic_prog_operators
   ARITHMETIC_OP_FILTER_MEAN,
   ARITHMETIC_OP_FILTER_SIGCLIP_MEAN,
   ARITHMETIC_OP_FILTER_SIGCLIP_MEDIAN,
+  ARITHMETIC_OP_FILTER_MADCLIP_MEAN,
+  ARITHMETIC_OP_FILTER_MADCLIP_MEDIAN,
   ARITHMETIC_OP_ERODE,
   ARITHMETIC_OP_DILATE,
   ARITHMETIC_OP_NUMBER_NEIGHBORS,
diff --git a/doc/announce-acknowledge.txt b/doc/announce-acknowledge.txt
index 2d0f3f99..580b5e01 100644
--- a/doc/announce-acknowledge.txt
+++ b/doc/announce-acknowledge.txt
@@ -6,7 +6,9 @@ Greg Wooledge
 Hamed Altafi
 Jesús Vega
 Juan Castillo Ramírez
+Mathias Urbano
 Ooldooz Kabood
+Paola Dimauro
 Phil Wyett
 Rahna Payyasseri Thanduparackal
 Raul Infante-Sainz
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index c794e9a6..72fccb1d 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -22287,10 +22287,23 @@ $ astarithmetic 3 0.2 5 4 image.fits 
filter-sigclip-mean
 
 The median (which needs a sorted dataset) is necessary for 
@mymath{\sigma}-clipping, therefore @code{filter-sigclip-mean} can be 
significantly slower than @code{filter-mean}.
 However, if there are strong outliers in the dataset that you want to ignore 
(for example, emission lines on a spectrum when finding the continuum), this is 
a much better solution.
+Note that @mymath{\sigma}-clipping is easily biased by outliers, so 
@code{filter-madclip-mean} is more robust.
 
 @item filter-sigclip-median
 Apply a @mymath{\sigma}-clipped median filtering onto the input dataset.
 This operator and its necessary operands are almost identical to 
@code{filter-sigclip-mean}, except that after @mymath{\sigma}-clipping, the 
median value (which is less affected by outliers than the mean) is added back 
to the stack.
+Note that @mymath{\sigma}-clipping is easily biased by outliers, so 
@code{filter-madclip-median} is more robust.
+
+@item filter-madclip-mean
+Apply a MAD-clipped mean filtering onto the input dataset.
+MAD stands for Median Absolute Deviation, this is much more robust to outliers 
than the @mymath{\sigma}-clipping.
+This operator is called in a similar way to @code{filter-madclip-mean}, for 
example in the command below, we use @mymath{4\times}MAD clipping with a 
termination criteria of 0.01 on an image:
+@example
+$ astarithmetic 4 0.01 5 5 image.fits filter-madclip-mean
+@end example
+
+@item filter-madclip-median
+Apply a MAD-clipped median filtering onto the input dataset; see the 
description of @code{filter-madclip-mean} for more.
 @end table
 
 @node Pooling operators, Interpolation operators, Filtering operators, 
Arithmetic operators