[gnuastro-commits] master aea662d: Mean and Median filtering added to Arithmetic program

[Mohammad Akhlaghi]

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

    Mean and Median filtering added to Arithmetic program
    
    The Arithmetic program (not library currently) now has two new operators:
    `filter-median' and `filter-mean'. Currently they are implemented as a
    function directly in `bin/arithmetic/arithmetic.c', but we will later move
    it to the libraries (either as a separate library in `gnuastro/filter.h' or
    as part of an existing header.
---
 NEWS                        |   4 +
 bin/arithmetic/arithmetic.c | 428 +++++++++++++++++++++++++++++++++++++++-----
 bin/arithmetic/arithmetic.h |  15 ++
 doc/gnuastro.texi           |  45 ++++-
 lib/arithmetic.c            |   2 +
 lib/gnuastro/arithmetic.h   |   2 +
 lib/statistics.c            |   2 +-
 7 files changed, 453 insertions(+), 45 deletions(-)

diff --git a/NEWS b/NEWS
index 6ea4787..0fa2662 100644
--- a/NEWS
+++ b/NEWS
@@ -8,6 +8,10 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
   All programs: a value of `0' to the `--numthreads' option will use the
   number of threads available to the system at run time.
 
+  Arithmetic: The new operators `filter-median' and `filter-mean' can be
+  used to filter (smooth) the input. The size of the filter can be set as
+  the other operands to these operators.
+
   BuildProgram: The new `--la' option allows the identification of a
   different Libtool `.la' file for Libtool linking information.
 
diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index b213e7d..36ed627 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -31,6 +31,9 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 #include <gnuastro/wcs.h>
 #include <gnuastro/fits.h>
+#include <gnuastro/threads.h>
+#include <gnuastro/dimension.h>
+#include <gnuastro/statistics.h>
 #include <gnuastro/arithmetic.h>
 
 #include <gnuastro-internal/checkset.h>
@@ -60,7 +63,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
     CTYPE a=*(CTYPE *)(data->array); if(a>0) return a;    }
 
 static size_t
-set_number_of_operands(struct arithmeticparams *p, gal_data_t *data,
+pop_number_of_operands(struct arithmeticparams *p, gal_data_t *data,
                        char *token_string)
 {
   /* Check if its a number. */
@@ -121,6 +124,313 @@ set_number_of_operands(struct arithmeticparams *p, 
gal_data_t *data,
 
 
 
+/**********************************************************************/
+/****************         Filtering operators         *****************/
+/**********************************************************************/
+#define ARITHMETIC_FILTER_DIM 10
+
+struct arithmetic_filter_p
+{
+  int      operator;            /* The type of filtering.          */
+  size_t     *fsize;            /* Filter size.                    */
+  size_t   *hpfsize;            /* Positive Half-filter size.      */
+  size_t   *hnfsize;            /* Negative Half-filter size.      */
+  gal_data_t *input;            /* Input dataset.                  */
+  gal_data_t   *out;            /* Output dataset.                 */
+
+  int      hasblank;            /* If the dataset has blank values.*/
+};
+
+
+
+
+
+/* Main filtering work function. */
+static void *
+arithmetic_filter(void *in_prm)
+{
+  struct gal_threads_params *tprm=(struct gal_threads_params *)in_prm;
+  struct arithmetic_filter_p *afp=(struct arithmetic_filter_p *)tprm->params;
+  gal_data_t *input=afp->input;
+
+  size_t ind,index;
+  int out_type_checked=0;
+  gal_data_t *result=NULL;
+  size_t *hpfsize=afp->hpfsize, *hnfsize=afp->hnfsize;
+  size_t *tsize, *dsize=input->dsize, *fsize=afp->fsize;
+  size_t i, j, coord[ARITHMETIC_FILTER_DIM], ndim=input->ndim;
+  size_t start[ARITHMETIC_FILTER_DIM], end[ARITHMETIC_FILTER_DIM];
+  gal_data_t *tile=gal_data_alloc(NULL, input->type, ndim, afp->fsize, NULL,
+                                  0, -1, NULL, NULL, NULL);
+
+  /* Prepare the tile. */
+  free(tile->array);
+  tsize=tile->dsize;
+  tile->block=input;
+
+
+  /* Go over all the pixels that were assigned to this thread. */
+  for(i=0; tprm->indexs[i] != GAL_BLANK_SIZE_T; ++i)
+    {
+      /* Get the coordinate of the pixel. */
+      ind=tprm->indexs[i];
+      gal_dimension_index_to_coord(ind, ndim, dsize, coord);
+
+      /* See which dimensions need trimming. */
+      tile->size=1;
+      for(j=0;j<ndim;++j)
+        {
+          /* Estimate the coordinate of the filter's starting point. Note
+             that we are dealing with size_t (unsigned int) type here, so
+             there are no negatives. A negative result will produce an
+             extremely large number, so instead of checking for negative,
+             we can just see if the result of a subtraction is less than
+             the width of the input. */
+          if( (coord[j] - hnfsize[j] > dsize[j])
+              || (coord[j] + hpfsize[j] >= dsize[j]) )
+            {
+              start[j] = ( (coord[j] - hnfsize[j] > dsize[j])
+                           ? 0 : coord[j] - hnfsize[j] );
+              end[j]   = ( (coord[j] + hpfsize[j] >= dsize[j])
+                           ? dsize[j]
+                           : coord[j] + hpfsize[j] + 1);
+              tsize[j] = end[j] - start[j];
+            }
+          else  /* We are NOT on the edge (given requested filter width). */
+            {
+              tsize[j] = fsize[j];
+              start[j] = coord[j] - hnfsize[j];
+            }
+          tile->size *= tsize[j];
+        }
+
+      /* Set the tile's starting pointer. */
+      index=gal_dimension_coord_to_index(ndim, dsize, start);
+      tile->array=gal_data_ptr_increment(input->array, index, input->type);
+
+      /* Do the necessary calculation. */
+      switch(afp->operator)
+        {
+        case ARITHMETIC_OP_FILTER_MEDIAN:
+          result=gal_statistics_median(tile, 0);
+          break;
+
+        case ARITHMETIC_OP_FILTER_MEAN:
+          result=gal_statistics_mean(tile);
+          break;
+
+        default:
+          error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
+                "fix the problem. `afp->operator' code %d is not "
+                "recognized", PACKAGE_BUGREPORT, __func__, afp->operator);
+        }
+
+      /* Make sure the output array type and result's type are the same. We
+         only need to do this once, but we'll suffice to once for each
+         thread for simplicify of the code, it is too negligible to have
+         any real affect. */
+      if( out_type_checked==0)
+        {
+          if(result->type!=afp->out->type )
+            error(EXIT_FAILURE, 0, "%s: a bug! please contact us at %s so "
+                  "we can address the problem. The tyes of `result' and "
+                  "`out' aren't the same, they are respectively: `%s' and "
+                  "`%s'", __func__, PACKAGE_BUGREPORT,
+                  gal_type_name(result->type, 1),
+                  gal_type_name(afp->out->type, 1));
+          out_type_checked=1;
+        }
+
+      /* Copy the result into the output array. */
+      memcpy(gal_data_ptr_increment(afp->out->array, ind, afp->out->type),
+             result->array, gal_type_sizeof(afp->out->type));
+
+      /* Clean up. */
+      gal_data_free(result);
+    }
+
+
+  /* Clean up. */
+  tile->array=NULL;
+  tile->block=NULL;
+  gal_data_free(tile);
+
+
+  /* Wait for all the other threads to finish, then return. */
+  if(tprm->b) pthread_barrier_wait(tprm->b);
+  return NULL;
+}
+
+
+
+
+
+static void
+wrapper_for_filter(struct arithmeticparams *p, char *token, int operator)
+{
+  size_t i=0, ndim, one=1;
+  int type=GAL_TYPE_INVALID;
+  struct arithmetic_filter_p afp={0};
+  size_t fsize[ARITHMETIC_FILTER_DIM];
+  gal_data_t *tmp, *tmp2, *zero, *comp, *fsize_list=NULL;
+  size_t hnfsize[ARITHMETIC_FILTER_DIM], hpfsize[ARITHMETIC_FILTER_DIM];
+
+
+  /* Get the input's number of dimensions. */
+  afp.input=operands_pop(p, token);
+  afp.operator=operator;
+  ndim=afp.input->ndim;
+  afp.hnfsize=hnfsize;
+  afp.hpfsize=hpfsize;
+  afp.fsize=fsize;
+
+
+  /* A small sanity check. */
+  if(ndim>ARITHMETIC_FILTER_DIM)
+    error(EXIT_FAILURE, 0, "%s: currently only datasets with less than "
+          "%d dimensions are acceptable. The input has %zu dimensions",
+          __func__, ARITHMETIC_FILTER_DIM, ndim);
+
+
+  /* A zero value for checking the value of input widths. */
+  zero=gal_data_alloc(NULL, GAL_TYPE_INT32, 1, &one, NULL, 1, -1, NULL,
+                      NULL, NULL);
+
+
+  /* Based on the dimensions of the popped operand, pop the necessary
+     number of operands. */
+  for(i=0;i<ndim;++i)
+    gal_list_data_add(&fsize_list, operands_pop(p, token));
+
+
+  /* Make sure the filter size only has single values. */
+  i=0;
+  for(tmp=fsize_list; tmp!=NULL; tmp=tmp->next)
+    {
+      ++i;
+      if(tmp->size!=1)
+        error(EXIT_FAILURE, 0, "the filter length values given to the "
+              "filter operators can only be numbers. Value number %zu has "
+              "%zu elements, so its an array", ndim-i-1, tmp->size);
+    }
+
+
+  /* If the input only has one element, filtering makes no sense, so don't
+     waste time, just add the input onto the stack. */
+  if(afp.input->size==1) afp.out=afp.input;
+  else
+    {
+      /* Allocate an array for the size of the filter and fill it in. The
+         values must be written in the inverse order since the user gives
+         dimensions with the FITS standard. */
+      i=ndim-1;
+      for(tmp=fsize_list; tmp!=NULL; tmp=tmp->next)
+        {
+          /* Make sure the user has given an integer type. */
+          if(tmp->type==GAL_TYPE_FLOAT32 || tmp->type==GAL_TYPE_FLOAT64)
+            error(EXIT_FAILURE, 0, "lengths of filter along dimensions "
+                  "must be integer values, not floats. The given length "
+                  "along dimension %zu is a float", ndim-i);
+
+          /* Make sure it isn't negative. */
+          comp=gal_arithmetic(GAL_ARITHMETIC_OP_GT, 0, tmp, zero);
+          if( *(uint8_t *)(comp->array) == 0 )
+            error(EXIT_FAILURE, 0, "lengths of filter along dimensions "
+                  "must be positive. The given length in dimension %zu"
+                  "is either zero or negative", ndim-i);
+
+          /* Convert the input into size_t and put it into the array that
+             keeps the filter size. */
+          tmp2=gal_data_copy_to_new_type(tmp, GAL_TYPE_SIZE_T);
+          fsize[ i ] = *(size_t *)(tmp2->array);
+          gal_data_free(tmp2);
+
+          /* If the width is larger than the input's size, change the width
+             to the input's size. */
+          if( fsize[i] > afp.input->dsize[i] )
+            error(EXIT_FAILURE, 0, "%s: the filter size along dimension %zu "
+                  "(%zu) is greater than the input's length in that "
+                  "dimension (%zu)", __func__, i, fsize[i],
+                  afp.input->dsize[i]);
+
+          /* Go onto the previous dimension. */
+          --i;
+        }
+
+
+      /* Set the half filter sizes. Note that when the size is an odd
+         number, the number of pixels before and after the actual pixel are
+         equal, but for an even number, we will look into one element more
+         when looking before than the ones after. */
+      for(i=0;i<ndim;++i)
+        {
+          if( fsize[i]%2 )
+            hnfsize[i]=hpfsize[i]=fsize[i]/2;
+          else
+            { hnfsize[i]=fsize[i]/2; hpfsize[i]=fsize[i]/2-1; }
+        }
+
+
+      /* See if the input has blank pixels. */
+      afp.hasblank=gal_blank_present(afp.input, 1);
+
+      /* Set the type of the output dataset. */
+      switch(operator)
+        {
+        case ARITHMETIC_OP_FILTER_MEDIAN:
+          type=afp.input->type;
+          break;
+
+        case ARITHMETIC_OP_FILTER_MEAN:
+          type=GAL_TYPE_FLOAT64;
+          break;
+
+        default:
+          error(EXIT_FAILURE, 0, "%s: a bug! please contact us at %s to fix "
+                "the problem. The `operator' code %d is not recognized",
+                PACKAGE_BUGREPORT, __func__, operator);
+        }
+
+      /* Allocate the output dataset. Note that filtering doesn't change
+         the units of the dataset. */
+      afp.out=gal_data_alloc(NULL, type, ndim, afp.input->dsize,
+                             afp.input->wcs, 0, afp.input->minmapsize,
+                             NULL, afp.input->unit, NULL);
+
+      /* Spin off threads for each pixel. */
+      gal_threads_spin_off(arithmetic_filter, &afp, afp.input->size,
+                           p->cp.numthreads);
+    }
+
+
+  /* Add the output to the top of the stack. */
+  operands_add(p, NULL, afp.out);
+
+
+  /* Clean up and add the output on top of the stack */
+  gal_data_free(afp.input);
+  gal_list_data_free(fsize_list);
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
 /***************************************************************/
 /*************      Reverse Polish algorithm       *************/
 /***************************************************************/
@@ -274,6 +584,13 @@ reversepolish(struct arithmeticparams *p)
           else if (!strcmp(token->v, "float64"))
             { op=GAL_ARITHMETIC_OP_TO_FLOAT64;    nop=1;  }
 
+          /* Filters. */
+          else if (!strcmp(token->v, "filter-median"))
+            { op=ARITHMETIC_OP_FILTER_MEDIAN;     nop=0;  }
+          else if (!strcmp(token->v, "filter-mean"))
+            { op=ARITHMETIC_OP_FILTER_MEAN;       nop=0;  }
+
+
           /* Finished checks with known operators */
           else
             error(EXIT_FAILURE, 0, "the argument \"%s\" could not be "
@@ -281,51 +598,76 @@ reversepolish(struct arithmeticparams *p)
                   token->v);
 
 
-          /* Pop the necessary number of operators. Note that the operators
-             are poped from a linked list (which is last-in-first-out). So
-             for the operators which need a specific order, the first poped
-             operand is actally the last (right most, in in-fix notation)
-             input operand.*/
-          switch(nop)
+          /* See if the arithmetic library must be called or not. */
+          if(nop)
             {
-            case 1:
-              d1=operands_pop(p, token->v);
-              break;
-
-            case 2:
-              d2=operands_pop(p, token->v);
-              d1=operands_pop(p, token->v);
-              break;
-
-            case 3:
-              d3=operands_pop(p, token->v);
-              d2=operands_pop(p, token->v);
-              d1=operands_pop(p, token->v);
-              break;
-
-            case -1:
-              /* This case is when the number of operands is itself an
-                 operand. So the first popped operand must be an integer
-                 number, we will use that to construct a linked list of any
-                 number of operands within the single `d1' pointer. */
-              d1=NULL;
-              numop=set_number_of_operands(p, operands_pop(p, token->v),
-                                           token->v);
-              for(i=0;i<numop;++i)
-                gal_list_data_add(&d1, operands_pop(p, token->v));
-              break;
-
-            default:
-              error(EXIT_FAILURE, 0, "No operators need %d operands", nop);
+              /* Pop the necessary number of operators. Note that the
+                 operators are poped from a linked list (which is
+                 last-in-first-out). So for the operators which need a
+                 specific order, the first poped operand is actally the
+                 last (right most, in in-fix notation) input operand.*/
+              switch(nop)
+                {
+                case 1:
+                  d1=operands_pop(p, token->v);
+                  break;
+
+                case 2:
+                  d2=operands_pop(p, token->v);
+                  d1=operands_pop(p, token->v);
+                  break;
+
+                case 3:
+                  d3=operands_pop(p, token->v);
+                  d2=operands_pop(p, token->v);
+                  d1=operands_pop(p, token->v);
+                  break;
+
+                case -1:
+                  /* This case is when the number of operands is itself an
+                     operand. So the first popped operand must be an
+                     integer number, we will use that to construct a linked
+                     list of any number of operands within the single `d1'
+                     pointer. */
+                  d1=NULL;
+                  numop=pop_number_of_operands(p, operands_pop(p, token->v),
+                                               token->v);
+                  for(i=0;i<numop;++i)
+                    gal_list_data_add(&d1, operands_pop(p, token->v));
+                  break;
+
+                default:
+                  error(EXIT_FAILURE, 0, "no operators, `%s' needs %d "
+                        "operand(s)", token->v, nop);
+                }
+
+
+              /* Run the arithmetic operation. Note that `gal_arithmetic'
+                 is a variable argument function (like printf). So the
+                 number of arguments it uses depend on the operator. So
+                 when the operator doesn't need three operands, the extra
+                 arguments will be ignored. */
+              operands_add(p, NULL, gal_arithmetic(op, flags, d1, d2, d3));
             }
 
-
-          /* Run the arithmetic operation. Note that `gal_arithmetic' is a
-             variable argument function (like printf). So the number of
-             arguments it uses depend on the operator. So when the operator
-             doesn't need three operands, the extra arguments will be
-             ignored. */
-          operands_add(p, NULL, gal_arithmetic(op, flags, d1, d2, d3));
+          /* No need to call the arithmetic library, call the proper
+             wrappers directly. */
+          else
+            {
+              switch(op)
+                {
+                case ARITHMETIC_OP_FILTER_MEAN:
+                case ARITHMETIC_OP_FILTER_MEDIAN:
+                  wrapper_for_filter(p, token->v, op);
+                  break;
+
+                default:
+                  error(EXIT_FAILURE, 0, "%s: a bug! please contact us at "
+                        "%s to fix the problem. The code %d is not "
+                        "recognized for `op'", __func__, PACKAGE_BUGREPORT,
+                        op);
+                }
+            }
         }
     }
 
diff --git a/bin/arithmetic/arithmetic.h b/bin/arithmetic/arithmetic.h
index e7f2ea0..833c5bf 100644
--- a/bin/arithmetic/arithmetic.h
+++ b/bin/arithmetic/arithmetic.h
@@ -23,6 +23,21 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #ifndef IMGARITH_H
 #define IMGARITH_H
 
+
+#include <gnuastro/arithmetic.h>
+
+
+/* These are operator codes for functions that aren't in the arithmetic
+   library. */
+enum arithmetic_prog_operators
+{
+  ARITHMETIC_OP_FILTER_MEDIAN=GAL_ARITHMETIC_OP_LAST_CODE,
+  ARITHMETIC_OP_FILTER_MEAN,
+};
+
+
+
+
 void
 imgarith(struct arithmeticparams *p);
 
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 305dcfc..daeb86e 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -8694,6 +8694,49 @@ standard deviation of the respective pixels in all 
operands in the stack.
 Similar to @command{min}, but the pixels of the output will contain
 the median of the respective pixels in all operands in the stack.
 
address@hidden filter-mean
+Apply mean filtering (or @url{https://en.wikipedia.org/wiki/Moving_average,
+moving average}) on the input dataset. During mean filtering, each pixel
+(data element) is replaced by the mean value of all its surrounding pixels
+(excluding blank values). The number of surrounding pixels in each
+dimension (to calculate the mean) is determined through the earlier
+operands that have been pushed onto the stack prior to the input
+dataset. The number of necessary operands is determined by the
+dimensionality of the input dataset (first popped operand). The order of
+the dimensions on the command-line is the order in FITS format. Here is one
+example:
+
address@hidden
+$ astarithmetic 5 4 image.fits filter-mean
address@hidden example
+
address@hidden
+In this example, each pixel is replaced by the mean of a 5 by 4 box around
+it. The box is 5 pixels along the first FITS dimension (horizontal when
+viewed in ds9) and 4 pixels along the second FITS dimension (vertical).
+
+Each pixel will be placed in the center of the box that the mean is
+calculated on. If the given width along a dimension is even, then the
+center is assumed to be between the pixels (not in the center of a
+pixel). When the pixel is close to the center, the pixels of the box that
+fall outside the image are ignored. Therefore, on the edge, less points
+will be used in calculating the mean.
+
+The final effect of mean filtering is to smooth the input image, it is
+essentially a convolution with a kernel that has identical values for all
+its pixels (is flat), see @ref{Convolution process}.
+
address@hidden filter-median
+Apply @url{https://en.wikipedia.org/wiki/Median_filter, median filtering}
+on the input dataset. This is very similar to @command{filter-mean}, except
+that instead of the mean value of the box pixels, the median value is used
+to replace a pixel value. For more on how to use this operator, please see
address@hidden
+
+The median is less susceptible to outliers compared to the mean. As a
+result, after median filtering, the pixel values will be more discontinuous
+than mean filtering.
+
 @item lt
 Less than: If the second popped (or left operand in infix notation, see
 @ref{Reverse polish notation}) value is smaller than the first popped
@@ -17481,7 +17524,7 @@ number of threads on the system and spin-off threads. 
You can see a
 demonstration of using these functions in @ref{Library demo -
 multi-threaded operation}.
 
address@hidden {C @code{struct})} gal_threads_params
address@hidden {C @code{struct}} gal_threads_params
 Structure keeping the parameters of each thread. When each thread is
 created, a pointer to this structure is passed to it. The @code{params}
 element can be the pointer to a structure defined by the user which
diff --git a/lib/arithmetic.c b/lib/arithmetic.c
index 6e20539..db61e9c 100644
--- a/lib/arithmetic.c
+++ b/lib/arithmetic.c
@@ -1484,6 +1484,8 @@ gal_arithmetic(int operator, unsigned char flags, ...)
       out=arithmetic_abs(flags, d1);
       break;
 
+
+    /* Multi-operand operators */
     case GAL_ARITHMETIC_OP_MIN:
     case GAL_ARITHMETIC_OP_MAX:
     case GAL_ARITHMETIC_OP_NUM:
diff --git a/lib/gnuastro/arithmetic.h b/lib/gnuastro/arithmetic.h
index 5d50518..ba20c6a 100644
--- a/lib/gnuastro/arithmetic.h
+++ b/lib/gnuastro/arithmetic.h
@@ -132,6 +132,8 @@ enum gal_arithmetic_operators
   GAL_ARITHMETIC_OP_TO_INT64,     /* Convert to int64_t.                   */
   GAL_ARITHMETIC_OP_TO_FLOAT32,   /* Convert to float32.                   */
   GAL_ARITHMETIC_OP_TO_FLOAT64,   /* Convert to float64.                   */
+
+  GAL_ARITHMETIC_OP_LAST_CODE,    /* Last code of the library operands.    */
 };
 
 
diff --git a/lib/statistics.c b/lib/statistics.c
index 6d5e013..a85ce88 100644
--- a/lib/statistics.c
+++ b/lib/statistics.c
@@ -1805,7 +1805,7 @@ gal_statistics_sigma_clip(gal_data_t *input, float 
multip, float param,
       meanstd=gal_statistics_mean_std(nbs);
 
       /* Put the three final values in usable (with a type) pointers. */
-      med = median_d->array;
+      med  = median_d->array;
       mean = meanstd->array;
       std  = &((double *)(meanstd->array))[1];