[gnuastro-commits] master 63b989f 2/2: Nearest neighbor interpolation in Arithmetic

[Mohammad Akhlaghi]

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

    Nearest neighbor interpolation in Arithmetic
    
    Using the (currently single function) interpolation library of Gnuastro,
    Arithmetic can now fill-in blank values using the nearest neighbors with
    the `interpolate-medianngb' operator.
---
 NEWS                        |  1 +
 bin/arithmetic/arithmetic.c | 50 +++++++++++++++++++++++++++++++++++++++++++++
 bin/arithmetic/arithmetic.h |  1 +
 doc/gnuastro.texi           |  7 +++++++
 lib/interpolate.c           |  4 ++--
 5 files changed, 61 insertions(+), 2 deletions(-)

diff --git a/NEWS b/NEWS
index 45356df..3ab0703 100644
--- a/NEWS
+++ b/NEWS
@@ -22,6 +22,7 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
     - filter-sigclip-median: sigma-clipped, median filter operator.
     - connected-components: label the connected elements of the input.
     - invert: subtract the maximum of unsigned types (absorption to emission).
+    - interpolate-medianngb: Interpolate (only blanks) with nearest neighbors.
 
   ConvertType:
     - TIFF images can also be used as input.
diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index 63487bc..ac1a183 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -37,6 +37,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/dimension.h>
 #include <gnuastro/statistics.h>
 #include <gnuastro/arithmetic.h>
+#include <gnuastro/interpolate.h>
 
 #include <gnuastro-internal/checkset.h>
 
@@ -642,6 +643,49 @@ arithmetic_invert(struct arithmeticparams *p, char *token)
 
 
 
+static void
+arithmetic_interpolate(struct arithmeticparams *p, char *token)
+{
+  int num_int;
+  gal_data_t *interpolated;
+
+  /* First pop the number of nearby neighbors.*/
+  gal_data_t *num = operands_pop(p, token);
+
+  /* Then pop the actual dataset to interpolate. */
+  gal_data_t *in = operands_pop(p, token);
+
+  /* Do proper sanity checks on `num'. */
+  if(num->size!=1)
+    error(EXIT_FAILURE, 0, "the first popped operand to "
+          "`interpolate-medianngb' must be a single number. However, "
+          "it has %zu elements", num->size);
+  if(num->type==GAL_TYPE_FLOAT32 || num->type==GAL_TYPE_FLOAT64)
+    error(EXIT_FAILURE, 0, "the first popped operand to "
+          "`interpolate-medianngb' is the number of nearby neighbors (a "
+          "counter, an integer). It must NOT be a floating point.\n\n"
+          "If its already an integer, but in a floating point container, "
+          "you can use the `int32' operator to convert it to a 32-bit "
+          "integer for example");
+
+  /* Convert the given number to a 32-bit integer and read it in. */
+  num=gal_data_copy_to_new_type_free(num, GAL_TYPE_INT32);
+  num_int = *((int32_t *)(num->array));
+
+  /* Call the interpolation function. */
+  interpolated=gal_interpolate_close_neighbors(in, NULL, num_int,
+                                               p->cp.numthreads, 1, 0);
+
+  /* Clean up and push the interpolated array onto the stack. */
+  gal_data_free(in);
+  gal_data_free(num);
+  operands_add(p, NULL, interpolated);
+}
+
+
+
+
+
 
 
 
@@ -828,6 +872,8 @@ reversepolish(struct arithmeticparams *p)
             { op=ARITHMETIC_OP_CONNECTED_COMPONENTS;  nop=0;  }
           else if (!strcmp(token->v, "invert"))
             { op=ARITHMETIC_OP_INVERT;                nop=0;  }
+          else if (!strcmp(token->v, "interpolate-medianngb"))
+            { op=ARITHMETIC_OP_INTERPOLATE_MEDIANNGB; nop=0;  }
 
 
           /* Finished checks with known operators */
@@ -915,6 +961,10 @@ reversepolish(struct arithmeticparams *p)
                   arithmetic_invert(p, token->v);
                   break;
 
+                case ARITHMETIC_OP_INTERPOLATE_MEDIANNGB:
+                  arithmetic_interpolate(p, token->v);
+                  break;
+
                 default:
                   error(EXIT_FAILURE, 0, "%s: a bug! please contact us at "
                         "%s to fix the problem. The code %d is not "
diff --git a/bin/arithmetic/arithmetic.h b/bin/arithmetic/arithmetic.h
index fcf48d1..e5eca69 100644
--- a/bin/arithmetic/arithmetic.h
+++ b/bin/arithmetic/arithmetic.h
@@ -39,6 +39,7 @@ enum arithmetic_prog_operators
   ARITHMETIC_OP_DILATE,
   ARITHMETIC_OP_CONNECTED_COMPONENTS,
   ARITHMETIC_OP_INVERT,
+  ARITHMETIC_OP_INTERPOLATE_MEDIANNGB,
 };
 
 
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 4f67c8a..6b8759c 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -10169,6 +10169,13 @@ dataset. This operator and its necessary operands are 
almost identical to
 median value (which is less affected by outliers than the mean) is added
 back to the stack.
 
address@hidden interpolate-medianngb
+Interpolate all the blank elements of the second popped operand with the
+median of its nearest non-blank neighbors. The number of the nearest
+non-blank neighbors used to calculate the median is given by the first
+popped operand. Note that the distance of the nearest non-blank neighbors
+is irrelevant in this interpolation.
+
 @item erode
 @cindex Erosion
 Erode the foreground pixels (with value @code{1}) of the input dataset
diff --git a/lib/interpolate.c b/lib/interpolate.c
index 1276aaa..8d66deb 100644
--- a/lib/interpolate.c
+++ b/lib/interpolate.c
@@ -333,8 +333,8 @@ gal_interpolate_close_neighbors(gal_data_t *input,
   int permute=(tl && tl->totchannels>1 && tl->workoverch);
 
 
-  /* If there are no blank values in the array we should only fill blank
-     values, then simply copy the input and abort. */
+  /* If there are no blank values in the array, AND we should only fill
+     blank values, then simply copy the input and abort. */
   if( (input->flag | GAL_DATA_FLAG_BLANK_CH)     /* Zero bit is meaningful.*/
       && !(input->flag | GAL_DATA_FLAG_HASBLANK) /* There are no blanks.   */
       && onlyblank )                             /* Only interpolate blank.*/