gnuastro-commits
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[gnuastro-commits] master 8ba1d77: Library (interpolte.h): Min and Max i


From: Mohammad Akhlaghi
Subject: [gnuastro-commits] master 8ba1d77: Library (interpolte.h): Min and Max in nearest neighbor
Date: Sun, 26 Jul 2020 20:57:25 -0400 (EDT)

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

    Library (interpolte.h): Min and Max in nearest neighbor
    
    Until now, the only operator for nearest neighbor interpolation was the
    median. But the low-level function was easily expandable to any type of
    statistic!
    
    So with this commit, we are now adding minimum and maximum also (for
    example using maximum of nearest neighbors). This was based on the need to
    interpolate between the blank pixels of saturated stars to avoid
    segmentation complications while working on a project with Zahra Sharbaf.
---
 NEWS                        |  19 ++++++--
 bin/arithmetic/arithmetic.c |  36 ++++++++++++---
 bin/arithmetic/arithmetic.h |   2 +
 bin/noisechisel/threshold.c |   7 +--
 bin/statistics/sky.c        |   7 +--
 bin/statistics/statistics.c |  11 ++---
 doc/gnuastro.texi           | 104 ++++++++++++++++++++++++--------------------
 lib/gnuastro/interpolate.h  |  30 +++++++++----
 lib/interpolate.c           |  60 ++++++++++++++++---------
 lib/options.c               |   8 ++--
 10 files changed, 185 insertions(+), 99 deletions(-)

diff --git a/NEWS b/NEWS
index 7d53e9a..22d813c 100644
--- a/NEWS
+++ b/NEWS
@@ -8,6 +8,11 @@ See the end of the file for license conditions.
 ** New features
 
   Arithmetic:
+   - New operators:
+     - interpolate-minngb: Fill blanks with minimum of nearest neighbors.
+     - interpolate-maxngb: Fill blanks with maximum of nearest neighbors.
+          This can be useful to fill the blank values of saturated stars
+          for example.
    - To force integers to floats, you can also put a '.' or '.0' after
      them. Until now, it was only possibly by putting an 'f' after
      them. Hence while '5' will be read as an integer, '5.', '5.0' or '5f'
@@ -47,6 +52,13 @@ See the end of the file for license conditions.
 
   Library:
    - Spectral lines library: SiIII, OIII, CIV, NV and rest of Lyman series.
+   - GAL_INTERPOLATE_NEIGHBORS_METRIC_RADIAL: Radial metric for interpolation.
+   - GAL_INTERPOLATE_NEIGHBORS_METRIC_MANHATTAN: Mahattan distance.
+   - GAL_INTERPOLATE_NEIGHBORS_METRIC_INVALID: For error-handling/completeness.
+   - GAL_INTERPOLATE_NEIGHBORS_FUNC_MIN: Use minimum for interpolation.
+   - GAL_INTERPOLATE_NEIGHBORS_FUNC_MAX: Use maximum for interpolation.
+   - GAL_INTERPOLATE_NEIGHBORS_FUNC_MEDIAN: Use median for interpolation.
+   - GAL_INTERPOLATE_NEIGHBORS_FUNC_INVALID: for error-handling/completeness.
    - GAL_CONFIG_HAVE_WCSLIB_DIS_H: if the host's WCSLIB supports distortions.
    - GAL_CONFIG_HAVE_WCSLIB_MJDREF: if the host's WCSLIB reads MJDREF keyword.
    - GAL_WCS_FLTERROR: Limit to identify a floating point error for WCS.
@@ -72,12 +84,13 @@ See the end of the file for license conditions.
      changed).
 
   Library:
-   - gal_type_string_to_number: Numbers ending in '.' or '.0' will be
-     parsed as floating point. Until now, it would only parse numbers as
-     floating point if they had non-zero decimals.
+   - gal_interpolate_neighbors: new name for gal_interpolate_close_neighbors.
    - gal_statistics_outlier_bydistance: new name for the old
      'gal_statistics_outlier_positive'. It can now use the same algorithm
      for negative outliers with a new argument.
+   - gal_type_string_to_number: Numbers ending in '.' or '.0' will be
+     parsed as floating point. Until now, it would only parse numbers as
+     floating point if they had non-zero decimals.
 
 ** Bugs fixed
   bug #58434: MakeCatalog crash when ordering is required and no usable pixels.
diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index 316ea6d..080516c 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -703,9 +703,9 @@ arithmetic_invert(struct arithmeticparams *p, char *token)
 
 
 static void
-arithmetic_interpolate(struct arithmeticparams *p, char *token)
+arithmetic_interpolate(struct arithmeticparams *p, int operator, char *token)
 {
-  int num_int;
+  int num_int, interpop;
   gal_data_t *interpolated;
 
   /* First pop the number of nearby neighbors.*/
@@ -731,10 +731,28 @@ arithmetic_interpolate(struct arithmeticparams *p, char 
*token)
   num=gal_data_copy_to_new_type_free(num, GAL_TYPE_INT32);
   num_int = *((int32_t *)(num->array));
 
+  /* Set the interpolation operator. */
+  switch(operator)
+    {
+    case ARITHMETIC_OP_INTERPOLATE_MINNGB:
+      interpop=GAL_INTERPOLATE_NEIGHBORS_FUNC_MIN;
+      break;
+    case ARITHMETIC_OP_INTERPOLATE_MAXNGB:
+      interpop=GAL_INTERPOLATE_NEIGHBORS_FUNC_MAX;
+      break;
+    case ARITHMETIC_OP_INTERPOLATE_MEDIANNGB:
+      interpop=GAL_INTERPOLATE_NEIGHBORS_FUNC_MEDIAN;
+      break;
+    default:
+      error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix "
+            "the problem. The operator code %d isn't recognized",
+            __func__, PACKAGE_BUGREPORT, operator);
+    }
+
   /* Call the interpolation function. */
-  interpolated=gal_interpolate_close_neighbors(in, NULL, p->cp.interpmetric,
-                                               num_int, p->cp.numthreads,
-                                               1, 0);
+  interpolated=gal_interpolate_neighbors(in, NULL, p->cp.interpmetric,
+                                         num_int, p->cp.numthreads,
+                                         1, 0, interpop);
 
   /* Clean up and push the interpolated array onto the stack. */
   gal_data_free(in);
@@ -1019,6 +1037,10 @@ arithmetic_set_operator(char *string, size_t 
*num_operands)
         { op=ARITHMETIC_OP_FILL_HOLES;            *num_operands=0; }
       else if (!strcmp(string, "invert"))
         { op=ARITHMETIC_OP_INVERT;                *num_operands=0; }
+      else if (!strcmp(string, "interpolate-minngb"))
+        { op=ARITHMETIC_OP_INTERPOLATE_MINNGB;    *num_operands=0; }
+      else if (!strcmp(string, "interpolate-maxngb"))
+        { op=ARITHMETIC_OP_INTERPOLATE_MAXNGB;    *num_operands=0; }
       else if (!strcmp(string, "interpolate-medianngb"))
         { op=ARITHMETIC_OP_INTERPOLATE_MEDIANNGB; *num_operands=0; }
       else if (!strcmp(string, "collapse-sum"))
@@ -1142,8 +1164,10 @@ arithmetic_operator_run(struct arithmeticparams *p, int 
operator,
           arithmetic_invert(p, operator_string);
           break;
 
+        case ARITHMETIC_OP_INTERPOLATE_MINNGB:
+        case ARITHMETIC_OP_INTERPOLATE_MAXNGB:
         case ARITHMETIC_OP_INTERPOLATE_MEDIANNGB:
-          arithmetic_interpolate(p, operator_string);
+          arithmetic_interpolate(p, operator, operator_string);
           break;
 
         case ARITHMETIC_OP_COLLAPSE_SUM:
diff --git a/bin/arithmetic/arithmetic.h b/bin/arithmetic/arithmetic.h
index 41d78a2..a18f70e 100644
--- a/bin/arithmetic/arithmetic.h
+++ b/bin/arithmetic/arithmetic.h
@@ -40,6 +40,8 @@ enum arithmetic_prog_operators
   ARITHMETIC_OP_CONNECTED_COMPONENTS,
   ARITHMETIC_OP_FILL_HOLES,
   ARITHMETIC_OP_INVERT,
+  ARITHMETIC_OP_INTERPOLATE_MINNGB,
+  ARITHMETIC_OP_INTERPOLATE_MAXNGB,
   ARITHMETIC_OP_INTERPOLATE_MEDIANNGB,
   ARITHMETIC_OP_COLLAPSE_SUM,
   ARITHMETIC_OP_COLLAPSE_MIN,
diff --git a/bin/noisechisel/threshold.c b/bin/noisechisel/threshold.c
index 698dd45..868cf60 100644
--- a/bin/noisechisel/threshold.c
+++ b/bin/noisechisel/threshold.c
@@ -276,9 +276,10 @@ threshold_interp_smooth(struct noisechiselparams *p, 
gal_data_t **first,
   /* Do the interpolation of both arrays. */
   (*first)->next = *second;
   if(third) (*second)->next = *third;
-  tmp=gal_interpolate_close_neighbors(*first, tl, cp->interpmetric,
-                                      cp->interpnumngb, cp->numthreads,
-                                      cp->interponlyblank, 1);
+  tmp=gal_interpolate_neighbors(*first, tl, cp->interpmetric,
+                                cp->interpnumngb, cp->numthreads,
+                                cp->interponlyblank, 1,
+                                GAL_INTERPOLATE_NEIGHBORS_FUNC_MEDIAN);
   gal_data_free(*first);
   gal_data_free(*second);
   if(third) gal_data_free(*third);
diff --git a/bin/statistics/sky.c b/bin/statistics/sky.c
index 757d2c0..a307aec 100644
--- a/bin/statistics/sky.c
+++ b/bin/statistics/sky.c
@@ -220,9 +220,10 @@ sky(struct statisticsparams *p)
   /* Interpolate the Sky and its standard deviation. */
   if(!cp->quiet) gettimeofday(&t1, NULL);
   p->sky_t->next=p->std_t;
-  tmp=gal_interpolate_close_neighbors(p->sky_t, tl, cp->interpmetric,
-                                      cp->interpnumngb, cp->numthreads,
-                                      cp->interponlyblank, 1);
+  tmp=gal_interpolate_neighbors(p->sky_t, tl, cp->interpmetric,
+                                cp->interpnumngb, cp->numthreads,
+                                cp->interponlyblank, 1,
+                                GAL_INTERPOLATE_NEIGHBORS_FUNC_MEDIAN);
   gal_data_free(p->sky_t);
   gal_data_free(p->std_t);
   p->sky_t=tmp;
diff --git a/bin/statistics/statistics.c b/bin/statistics/statistics.c
index 70ed6c8..9e5f6b7 100644
--- a/bin/statistics/statistics.c
+++ b/bin/statistics/statistics.c
@@ -271,11 +271,12 @@ statistics_interpolate_and_write(struct statisticsparams 
*p,
   if( p->interpolate
       && !(p->cp.interponlyblank && gal_blank_present(values, 1)==0) )
     {
-      interpd=gal_interpolate_close_neighbors(values, &cp->tl,
-                                              cp->interpmetric,
-                                              cp->interpnumngb,
-                                              cp->numthreads,
-                                              cp->interponlyblank, 0);
+      interpd=gal_interpolate_neighbors(values, &cp->tl,
+                                        cp->interpmetric,
+                                        cp->interpnumngb,
+                                        cp->numthreads,
+                                        cp->interponlyblank, 0,
+                                        GAL_INTERPOLATE_NEIGHBORS_FUNC_MEDIAN);
       gal_data_free(values);
       values=interpd;
     }
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 4b5fece..34f3c7c 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -10567,6 +10567,13 @@ Interpolate all the blank elements of the second 
popped operand with the median
 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 interpolate-minngb
+Similar to @code{interpolate-medianngb}, but will fill the blank values of the 
dataset with the minimum value of the nearest neighbors.
+
+@item interpolate-maxngb
+Similar to @code{interpolate-medianngb}, but will fill the blank values of the 
dataset with the maximum value of the nearest neighbors.
+One useful implementation of this operator is to fill the saturated pixels of 
stars in images.
+
 @item collapse-sum
 Collapse the given dataset (second popped operand), by summing all elements 
along the first popped operand (a dimension in FITS standard: counting from 
one, from fastest dimension).
 The returned dataset has one dimension less compared to the input.
@@ -24992,54 +24999,59 @@ is much faster.
 
 @cindex Sky line
 @cindex Interpolation
-During data analysis, it happens that parts of the data cannot be given a
-value, but one is necessary for the higher-level analysis. For example a
-very bright star saturated part of your image and you need to fill in the
-saturated pixels with some values. Another common usage case are masked
-sky-lines in 1D spectra that similarly need to be assigned a value for
-higher-level analysis. In other situations, you might want a value in an
-arbitrary point: between the elements/pixels where you have data. The
-functions described in this section are for such operations.
+During data analysis, it happens that parts of the data cannot be given a 
value, but one is necessary for the higher-level analysis.
+For example a very bright star saturated part of your image and you need to 
fill in the saturated pixels with some values.
+Another common usage case are masked sky-lines in 1D spectra that similarly 
need to be assigned a value for higher-level analysis.
+In other situations, you might want a value in an arbitrary point: between the 
elements/pixels where you have data.
+The functions described in this section are for such operations.
 
 @cindex GNU Scientific Library
-The parametric interpolations discussed below are wrappers around the
-interpolation functions of the GNU Scientific Library (or GSL, see @ref{GNU
-Scientific Library}). To identify the different GSL interpolation types,
-Gnuastro's @file{gnuastro/interpolate.h} header file contains macros that
-are discussed below. The GSL wrappers provided here are not yet complete
-because we are too busy. If you need them, please consider helping us in
-adding them to Gnuastro's library. Your would be very welcome and
-appreciated.
-
-@deftypefun {gal_data_t *} gal_interpolate_close_neighbors (gal_data_t 
@code{*input}, struct gal_tile_two_layer_params @code{*tl}, size_t 
@code{numneighbors}, size_t @code{numthreads}, int @code{onlyblank}, int 
@code{aslinkedlist})
-Interpolate the values in the image using the median value of their
-@code{numneighbors} closest neighbors. This function is non-parametric and
-thus agnostic to the input's number of dimension. If @code{onlyblank} is
-non-zero, then only blank elements will be interpolated and pixels that
-already have a value will be left untouched. This function is
-multi-threaded and will run on @code{numthreads} threads (see
-@code{gal_threads_number} in @ref{Multithreaded programming}).
-
-@code{tl} is Gnuastro's two later tessellation structure used to define
-tiles over an image and is fully described in @ref{Tile grid}. When
-@code{tl!=NULL}, then it is assumed that the @code{input->array} contains
-one value per tile and interpolation will respect certain tessellation
-properties, for example to not interpolate over channel borders.
-
-If several datasets have the same set of blank values, you don't need to
-call this function multiple times. When @code{aslinkedlist} is non-zero,
-then @code{input} will be seen as a @ref{List of gal_data_t}. In this case,
-the same neighbors will be used for all the datasets in the list. Of
-course, the values for each dataset will be different, so a different value
-will be written in the each dataset, but the neighbor checking that is the
-most CPU intensive part will only be done once.
-
-This is a non-parametric and robust function for interpolation. The
-interpolated values are also always within the range of the non-blank
-values and strong outliers do not get created. However, this type of
-interpolation must be used with care when there are gradients. This is
-because it is non-parametric and if there aren't enough neighbors,
-step-like features can be created.
+The parametric interpolations discussed below are wrappers around the 
interpolation functions of the GNU Scientific Library (or GSL, see @ref{GNU 
Scientific Library}).
+To identify the different GSL interpolation types, Gnuastro's 
@file{gnuastro/interpolate.h} header file contains macros that are discussed 
below.
+The GSL wrappers provided here are not yet complete because we are too busy.
+If you need them, please consider helping us in adding them to Gnuastro's 
library.
+Your would be very welcome and appreciated.
+
+@deffn Macro GAL_INTERPOLATE_NEIGHBORS_METRIC_RADIAL
+@deffnx Macro GAL_INTERPOLATE_NEIGHBORS_METRIC_MANHATTAN
+@deffnx Macro GAL_INTERPOLATE_NEIGHBORS_METRIC_INVALID
+The metric used to find distance for nearest neighbor interpolation.
+A radial metric uses the simple Euclidean function to find the distance 
between two pixels.
+A manhattan metric will always be an integer and is like steps (but is also 
much faster to calculate than radial metric because it doesn't need a square 
root calculation).
+@end deffn
+
+@deffn Macro GAL_INTERPOLATE_NEIGHBORS_FUNC_MIN
+@deffnx Macro GAL_INTERPOLATE_NEIGHBORS_FUNC_MAX
+@deffnx Macro GAL_INTERPOLATE_NEIGHBORS_FUNC_MEDIAN
+@deffnx Macro GAL_INTERPOLATE_NEIGHBORS_FUNC_INVALID
+@cindex Saturated stars
+The various types of nearest-neighbor interporation functions for 
@code{gal_interpolate_neighbors}.
+The names are descriptive for the operation they do, so we won't go into much 
more detail here.
+The median operator will be one of the most used, but operators like the 
maximum are good to fill the center of saturated stars.
+@end deffn
+
+@deftypefun {gal_data_t *} gal_interpolate_neighbors (gal_data_t 
@code{*input}, struct gal_tile_two_layer_params @code{*tl}, uint8_t 
@code{metric}, size_t @code{numneighbors}, size_t @code{numthreads}, int 
@code{onlyblank}, int @code{aslinkedlist}, int @code{function})
+
+Interpolate the values in the input dataset using a calculated statistics from 
the distribution of their @code{numneighbors} closest neighbors.
+The desired statistics is determined from the @code{func} argument, which 
takes any of the @code{GAL_INTERPOLATE_NEIGHBORS_FUNC_} macros (see above).
+This function is non-parametric and thus agnostic to the input's number of 
dimension or shape of the distribution.
+
+Distance can be defined on different metrics that are dentified through 
@code{metric} (taking values determined by the 
@code{GAL_INTERPOLATE_NEIGHBORS_METRIC_} macros described above).
+If @code{onlyblank} is non-zero, then only blank elements will be interpolated 
and pixels that already have a value will be left untouched.
+This function is multi-threaded and will run on @code{numthreads} threads (see 
@code{gal_threads_number} in @ref{Multithreaded programming}).
+
+@code{tl} is Gnuastro's tessellation structure used to define tiles over an 
image and is fully described in @ref{Tile grid}.
+When @code{tl!=NULL}, then it is assumed that the @code{input->array} contains 
one value per tile and interpolation will respect certain tessellation 
properties, for example to not interpolate over channel borders.
+
+If several datasets have the same set of blank values, you don't need to call 
this function multiple times.
+When @code{aslinkedlist} is non-zero, then @code{input} will be seen as a 
@ref{List of gal_data_t}.
+In this case, the same neighbors will be used for all the datasets in the list.
+Of course, the values for each dataset will be different, so a different value 
will be written in the each dataset, but the neighbor checking that is the most 
CPU intensive part will only be done once.
+
+This is a non-parametric and robust function for interpolation.
+The interpolated values are also always within the range of the non-blank 
values and strong outliers do not get created.
+However, this type of interpolation must be used with care when there are 
gradients.
+This is because it is non-parametric and if there aren't enough neighbors, 
step-like features can be created.
 @end deftypefun
 
 @deffn Macro GAL_INTERPOLATE_1D_INVALID
diff --git a/lib/gnuastro/interpolate.h b/lib/gnuastro/interpolate.h
index 65fc48e..2b84d0d 100644
--- a/lib/gnuastro/interpolate.h
+++ b/lib/gnuastro/interpolate.h
@@ -48,12 +48,24 @@ __BEGIN_C_DECLS  /* From C++ preparations */
 
 
 /* Metrics to use for nearest-neighbor  */
-enum gal_interpolate_close_metric
+enum gal_interpolate_neighbors_metric
 {
- GAL_INTERPOLATE_CLOSE_METRIC_INVALID,
+ GAL_INTERPOLATE_NEIGHBORS_METRIC_INVALID,
 
- GAL_INTERPOLATE_CLOSE_METRIC_RADIAL,
- GAL_INTERPOLATE_CLOSE_METRIC_MANHATTAN,
+ GAL_INTERPOLATE_NEIGHBORS_METRIC_RADIAL,
+ GAL_INTERPOLATE_NEIGHBORS_METRIC_MANHATTAN,
+};
+
+
+
+/* Types of interpolation for nearest-neighbor. */
+enum gal_interpolate_neighbors_func
+{
+ GAL_INTERPOLATE_NEIGHBORS_FUNC_INVALID,
+
+ GAL_INTERPOLATE_NEIGHBORS_FUNC_MIN,
+ GAL_INTERPOLATE_NEIGHBORS_FUNC_MAX,
+ GAL_INTERPOLATE_NEIGHBORS_FUNC_MEDIAN,
 };
 
 
@@ -75,11 +87,11 @@ enum gal_interpolate_1D_types
 
 
 gal_data_t *
-gal_interpolate_close_neighbors(gal_data_t *input,
-                                struct gal_tile_two_layer_params *tl,
-                                uint8_t metric, size_t numneighbors,
-                                size_t numthreads, int onlyblank,
-                                int aslinkedlist);
+gal_interpolate_neighbors(gal_data_t *input,
+                             struct gal_tile_two_layer_params *tl,
+                             uint8_t metric, size_t numneighbors,
+                             size_t numthreads, int onlyblank,
+                             int aslinkedlist, int function);
 
 gsl_spline *
 gal_interpolate_1d_make_gsl_spline(gal_data_t *X, gal_data_t *Y, int type_1d);
diff --git a/lib/interpolate.c b/lib/interpolate.c
index c57dd96..341ddcf 100644
--- a/lib/interpolate.c
+++ b/lib/interpolate.c
@@ -57,8 +57,9 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 
 /* Parameters for interpolation on threads. */
-struct interpolate_params
+struct interpolate_ngb_params
 {
+  int                         function;
   gal_data_t                    *input;
   size_t                           num;
   gal_data_t                      *out;
@@ -78,11 +79,14 @@ struct interpolate_params
 
 /* Run the interpolation on many threads. */
 static void *
-interpolate_close_neighbors_on_thread(void *in_prm)
+interpolate_neighbors_on_thread(void *in_prm)
 {
-  /* Variables that others depend on. */
+  /* Low-level variables that others depend on. */
   struct gal_threads_params *tprm=(struct gal_threads_params *)in_prm;
-  struct interpolate_params *prm=(struct interpolate_params *)(tprm->params);
+  struct interpolate_ngb_params *prm=
+    (struct interpolate_ngb_params *)(tprm->params);
+
+  /* Higher-level variables. */
   struct gal_tile_two_layer_params *tl=prm->tl;
   int correct_index=(tl && tl->totchannels>1 && !tl->workoverch);
   gal_data_t *input=prm->input;
@@ -95,7 +99,7 @@ interpolate_close_neighbors_on_thread(void *in_prm)
   gal_list_dosizet_t *lQ, *sQ;
   size_t ngb_counter, pind, *dinc;
   size_t i, index, fullind, chstart=0, ndim=input->ndim;
-  gal_data_t *median, *tin, *tout, *tnear, *nearest=NULL;
+  gal_data_t *value, *tin, *tout, *tnear, *nearest=NULL;
   size_t size = (correct_index ? tl->tottilesinch : input->size);
   size_t *dsize = (correct_index ? tl->numtilesinch : input->dsize);
   size_t *icoord=gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__,
@@ -255,19 +259,34 @@ interpolate_close_neighbors_on_thread(void *in_prm)
                   "interpolation", __func__, ngb_counter, prm->numneighbors);
         }
 
-      /* Calculate the median of the values and write it in the output. */
+      /* Calculate the desired statistic, and write it in the output. */
       tout=prm->out;
       for(tnear=nearest; tnear!=NULL; tnear=tnear->next)
         {
-          /* Find the median and copy it, but first, reset the flags (which
-             remain from the last time). */
+          /* Find the desired statistic and copy it, but first, reset the
+             flags (which remain from the last time). */
           tnear->flag &= ~(GAL_DATA_FLAG_SORT_CH | GAL_DATA_FLAG_BLANK_CH);
-          median=gal_statistics_median(tnear, 1);
+          switch(prm->function)
+            {
+            case GAL_INTERPOLATE_NEIGHBORS_FUNC_MIN:
+              value=gal_statistics_minimum(tnear); break;
+              break;
+            case GAL_INTERPOLATE_NEIGHBORS_FUNC_MAX:
+              value=gal_statistics_maximum(tnear); break;
+              break;
+            case GAL_INTERPOLATE_NEIGHBORS_FUNC_MEDIAN:
+              value=gal_statistics_median(tnear, 1); break;
+            default:
+              error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s "
+                    "to fix the problem. The value %d is not a recognized "
+                    "interpolation function identifier", __func__,
+                    PACKAGE_BUGREPORT, prm->function);
+            }
           memcpy(gal_pointer_increment(tout->array, fullind, tout->type),
-                 median->array, gal_type_sizeof(tout->type));
+                 value->array, gal_type_sizeof(tout->type));
 
           /* Clean up and go to next array. */
-          gal_data_free(median);
+          gal_data_free(value);
           tout=tout->next;
         }
     }
@@ -328,14 +347,14 @@ interpolate_copy_input(gal_data_t *input, int 
aslinkedlist)
    assumed that the tile values correspond to given tessellation. Such that
    'input[i]' corresponds to 'tiles[i]' in the tessellation. */
 gal_data_t *
-gal_interpolate_close_neighbors(gal_data_t *input,
-                                struct gal_tile_two_layer_params *tl,
-                                uint8_t metric, size_t numneighbors,
-                                size_t numthreads, int onlyblank,
-                                int aslinkedlist)
+gal_interpolate_neighbors(gal_data_t *input,
+                             struct gal_tile_two_layer_params *tl,
+                             uint8_t metric, size_t numneighbors,
+                             size_t numthreads, int onlyblank,
+                             int aslinkedlist, int function)
 {
   gal_data_t *tin, *tout;
-  struct interpolate_params prm;
+  struct interpolate_ngb_params prm;
   size_t ngbvnum=numthreads*numneighbors;
   int permute=(tl && tl->totchannels>1 && tl->workoverch);
 
@@ -352,6 +371,7 @@ gal_interpolate_close_neighbors(gal_data_t *input,
   prm.tl           = tl;
   prm.ngb_vals     = NULL;
   prm.input        = input;
+  prm.function     = function;
   prm.onlyblank    = onlyblank;
   prm.numneighbors = numneighbors;
   prm.num          = aslinkedlist ? gal_list_data_number(input) : 1;
@@ -360,10 +380,10 @@ gal_interpolate_close_neighbors(gal_data_t *input,
   /* Set the metric. */
   switch(metric)
     {
-    case GAL_INTERPOLATE_CLOSE_METRIC_RADIAL:
+    case GAL_INTERPOLATE_NEIGHBORS_METRIC_RADIAL:
       prm.metric=gal_dimension_dist_radial;
       break;
-    case GAL_INTERPOLATE_CLOSE_METRIC_MANHATTAN:
+    case GAL_INTERPOLATE_NEIGHBORS_METRIC_MANHATTAN:
       prm.metric=gal_dimension_dist_manhattan;
       break;
     default:
@@ -438,7 +458,7 @@ gal_interpolate_close_neighbors(gal_data_t *input,
 
 
   /* Spin off the threads. */
-  gal_threads_spin_off(interpolate_close_neighbors_on_thread, &prm,
+  gal_threads_spin_off(interpolate_neighbors_on_thread, &prm,
                        input->size, numthreads);
 
 
diff --git a/lib/options.c b/lib/options.c
index 11a7f48..a22a5ec 100644
--- a/lib/options.c
+++ b/lib/options.c
@@ -521,10 +521,10 @@ gal_options_read_interpmetric(struct argp_option *option, 
char *arg,
     {
       switch(*(uint8_t *)(option->value))
         {
-        case GAL_INTERPOLATE_CLOSE_METRIC_RADIAL:
+        case GAL_INTERPOLATE_NEIGHBORS_METRIC_RADIAL:
           gal_checkset_allocate_copy("radial", &str);
           break;
-        case GAL_INTERPOLATE_CLOSE_METRIC_MANHATTAN:
+        case GAL_INTERPOLATE_NEIGHBORS_METRIC_MANHATTAN:
           gal_checkset_allocate_copy("manhattan", &str);
           break;
         default:
@@ -542,9 +542,9 @@ gal_options_read_interpmetric(struct argp_option *option, 
char *arg,
 
       /* Set the value. */
       if(       !strcmp(arg, "radial") )
-        *(uint8_t *)(option->value) = GAL_INTERPOLATE_CLOSE_METRIC_RADIAL;
+        *(uint8_t *)(option->value) = GAL_INTERPOLATE_NEIGHBORS_METRIC_RADIAL;
       else if ( !strcmp(arg, "manhattan") )
-        *(uint8_t *)(option->value) = GAL_INTERPOLATE_CLOSE_METRIC_MANHATTAN;
+        *(uint8_t *)(option->value) = 
GAL_INTERPOLATE_NEIGHBORS_METRIC_MANHATTAN;
       else
         error_at_line(EXIT_FAILURE, 0, filename, lineno, "'%s' (value to "
                       "'%s' option) isn't valid. Currently only 'radial' "



reply via email to

[Prev in Thread] Current Thread [Next in Thread]