[gnuastro-commits] master 0dbf34e: Library (interpolate): Ability to define the nearest-ngb metric

[Mohammad Akhlaghi]

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

    Library (interpolate): Ability to define the nearest-ngb metric
    
    Until now, the nearest neighbor interpolation library would only use the
    Manhattan/taxicab distance between two points. But this could cause small
    `x'-like artifacts in the final result.
    
    With this commit, its also possible to ask for a radial metric. As a
    result, the programs that use this function (Arithmetic, NoiseChisel and
    Statistics) also needed a new option and thus the `interpmetric' was added
    as a tessellation option.
    
    This solution was implemented after a discussion with Ignacio Trujillo.
---
 NEWS                               | 15 +++++++++++-
 bin/arithmetic/arithmetic.c        |  5 ++--
 bin/gnuastro.conf                  |  1 +
 bin/noisechisel/threshold.c        |  5 ++--
 bin/noisechisel/ui.c               |  1 +
 bin/statistics/sky.c               |  5 ++--
 bin/statistics/statistics.c        |  1 +
 bin/statistics/ui.c                |  8 +++---
 doc/gnuastro.texi                  | 21 ++++++++++++++++
 lib/dimension.c                    | 14 ++++++++++-
 lib/gnuastro-internal/commonopts.h | 14 +++++++++++
 lib/gnuastro-internal/options.h    |  6 +++++
 lib/gnuastro/dimension.h           |  5 ++--
 lib/gnuastro/interpolate.h         | 17 +++++++++++--
 lib/interpolate.c                  | 28 +++++++++++++++++----
 lib/options.c                      | 50 ++++++++++++++++++++++++++++++++++++++
 16 files changed, 175 insertions(+), 21 deletions(-)

diff --git a/NEWS b/NEWS
index 92021bd..398d8fb 100644
--- a/NEWS
+++ b/NEWS
@@ -26,6 +26,9 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
      different file for reading the WCS of the output. This is useful when
      the default (theh WCS of the first dataset that is read) is not the
      required one.
+   --interpmetric: new option that is necessary with the
+     `interpolate-medianngb' operator. For more, see the description of
+     this option in NoiseChisel.
 
   Fits:
    - Add "title" to group FITS keywords with `--write=/,"title name". This
@@ -51,7 +54,13 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
      columns into one.
 
   NoiseChisel:
-
+   --interpmetric: Set the metric to use to identify the nearest neighbors
+     for tile interpolation (quantile threshold, initial Sky, and final
+     Sky). Until now only the manhattan/taxicab metric was available, which
+     is fast, but could cause 45-degree lines in the interpolation. From
+     this version, with this option, its also possible to use the radial
+     distance (which is now the default). Just note that if there are many
+     tiles over the image, the radial distance will be slower.
    --snthresh: Manually set the signal-to-noise ratio of true
      pseudo-detections. With this option, NoiseChisel will not attempt to
      find pseudo-detections over the noisy regions of the dataset, but will
@@ -67,10 +76,14 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
      --dopening: The number of openings to do after applying `--dthresh'.
      --dopeningngb: The connectivity (4 or 8) to use for `--dopening'.
 
+  Statistics:
+   --interpmetric: Similar to NoiseChisel.
+
   Library:
     GAL_BLANK_LONG: new macro for the `long' type blank value.
     GAL_BLANK_ULONG: new macro for the `unsigned long' type blank value.
     gal_blank_number: Return the number of blank elements in a dataset.
+    gal_dimension_dist_radial: Radial distance between two coordinates.
     gal_statistics_outlier_cumulative: Uses flatness of the cumulative
        distribution to find outliers.
 
diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index 68e3912..46ca143 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -728,8 +728,9 @@ arithmetic_interpolate(struct arithmeticparams *p, char 
*token)
   num_int = *((int32_t *)(num->array));
 
   /* Call the interpolation function. */
-  interpolated=gal_interpolate_close_neighbors(in, NULL, num_int,
-                                               p->cp.numthreads, 1, 0);
+  interpolated=gal_interpolate_close_neighbors(in, NULL, p->cp.interpmetric,
+                                               num_int, p->cp.numthreads,
+                                               1, 0);
 
   /* Clean up and push the interpolated array onto the stack. */
   gal_data_free(in);
diff --git a/bin/gnuastro.conf b/bin/gnuastro.conf
index 3652e6a..a901b99 100644
--- a/bin/gnuastro.conf
+++ b/bin/gnuastro.conf
@@ -28,6 +28,7 @@
  numchannels      1,1
  remainderfrac    0.1
  workoverch       0
+ interpmetric     radial
  interpnumngb     9
  interponlyblank  0
 
diff --git a/bin/noisechisel/threshold.c b/bin/noisechisel/threshold.c
index 47dcd19..4fe408a 100644
--- a/bin/noisechisel/threshold.c
+++ b/bin/noisechisel/threshold.c
@@ -282,8 +282,9 @@ 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->interpnumngb,
-                                      cp->numthreads, cp->interponlyblank, 1);
+  tmp=gal_interpolate_close_neighbors(*first, tl, cp->interpmetric,
+                                      cp->interpnumngb, cp->numthreads,
+                                      cp->interponlyblank, 1);
   gal_data_free(*first);
   gal_data_free(*second);
   if(third) gal_data_free(*third);
diff --git a/bin/noisechisel/ui.c b/bin/noisechisel/ui.c
index 318963a..6248748 100644
--- a/bin/noisechisel/ui.c
+++ b/bin/noisechisel/ui.c
@@ -135,6 +135,7 @@ ui_initialize_options(struct noisechiselparams *p,
         case GAL_OPTIONS_KEY_TILESIZE:
         case GAL_OPTIONS_KEY_MINMAPSIZE:
         case GAL_OPTIONS_KEY_NUMCHANNELS:
+        case GAL_OPTIONS_KEY_INTERPMETRIC:
         case GAL_OPTIONS_KEY_INTERPNUMNGB:
         case GAL_OPTIONS_KEY_REMAINDERFRAC:
           cp->coptions[i].mandatory=GAL_OPTIONS_MANDATORY;
diff --git a/bin/statistics/sky.c b/bin/statistics/sky.c
index 6d6a3fe..865cfc4 100644
--- a/bin/statistics/sky.c
+++ b/bin/statistics/sky.c
@@ -220,8 +220,9 @@ 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->interpnumngb,
-                                      cp->numthreads, cp->interponlyblank, 1);
+  tmp=gal_interpolate_close_neighbors(p->sky_t, tl, cp->interpmetric,
+                                      cp->interpnumngb, cp->numthreads,
+                                      cp->interponlyblank, 1);
   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 c92ca8b..333b432 100644
--- a/bin/statistics/statistics.c
+++ b/bin/statistics/statistics.c
@@ -251,6 +251,7 @@ statistics_interpolate_and_write(struct statisticsparams *p,
       && !(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);
diff --git a/bin/statistics/ui.c b/bin/statistics/ui.c
index dfa3c21..4e1c7d4 100644
--- a/bin/statistics/ui.c
+++ b/bin/statistics/ui.c
@@ -418,10 +418,10 @@ ui_read_check_only_options(struct statisticsparams *p)
     {
       /* Mandatory options. */
       if( isnan(p->meanmedqdiff) || isnan(p->sclipparams[0])
-          || p->cp.interpnumngb==0 )
-        error(EXIT_FAILURE, 0, "`--meanmedqdiff', `--sclipparams' and "
-              "`--interpnumngb' are mandatory when requesting Sky "
-              "measurement (`--sky')");
+          || p->cp.interpmetric==0 || p->cp.interpnumngb==0 )
+        error(EXIT_FAILURE, 0, "`--meanmedqdiff', `--sclipparams', "
+              "`--interpmetric' and `--interpnumngb' are mandatory when "
+              "requesting Sky measurement (`--sky')");
 
       /* If mode and median distance is a reasonable value. */
       if(p->meanmedqdiff>0.5)
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 1469ede..4c90c2a 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -7816,6 +7816,22 @@ is desirable).
 When values are to be interpolated, only change the values of the blank
 elements, keep the non-blank elements untouched.
 
address@hidden --interpmetric=STR
address@hidden Radial metric
address@hidden Taxicab metric
address@hidden Manhattan metric
address@hidden Metric: Manhattan, Taxicab, Radial
+The metric to use for finding nearest neighbors. Currently it only accepts
+the Manhattan (or taxicab) metric with @code{manhattan}, or the radial
+metric with @code{radial}.
+
+The Manhattan distance between two points is defined with
address@hidden|\Delta{x}|+|\Delta{y}|}. Thus the Manhattan metric has the
+advantage of being fast, but at the expense of being less accurate. The
+radial distance is the standard definition of distance in a Euclidean
+space: @mymath{\sqrt{\Delta{x}^2+\Delta{y}^2}}. It is accurate, but the
+multiplication and square root can slow down the processing.
+
 @item --interpnumngb=INT
 The number of nearby non-blank neighbors to use for interpolation.
 @end table
@@ -24497,6 +24513,11 @@ the two coordinates @code{a} and @code{b} (each an 
array of @code{ndim}
 elements).
 @end deftypefun
 
address@hidden float gal_dimension_dist_radial (size_t @code{*a}, size_t 
@code{*b}, size_t @code{ndim})
+Return the radial distance between the two coordinates @code{a} and
address@hidden (each an array of @code{ndim} elements).
address@hidden deftypefun
+
 @deftypefun {gal_data_t *} gal_dimension_collapse_sum (gal_data_t @code{*in}, 
size_t @code{c_dim}, gal_data_t @code{*weight})
 Collapse the input dataset (@code{in}) along the given dimension
 (@code{c_dim}, in C definition: starting from zero, from the slowest
diff --git a/lib/dimension.c b/lib/dimension.c
index 070265e..310b5cb 100644
--- a/lib/dimension.c
+++ b/lib/dimension.c
@@ -263,7 +263,7 @@ gal_dimension_index_to_coord(size_t index, size_t ndim, 
size_t *dsize,
 /************************************************************************/
 /********************           Distances          **********************/
 /************************************************************************/
-size_t
+float
 gal_dimension_dist_manhattan(size_t *a, size_t *b, size_t ndim)
 {
   size_t i, out=0;
@@ -275,6 +275,18 @@ gal_dimension_dist_manhattan(size_t *a, size_t *b, size_t 
ndim)
 
 
 
+float
+gal_dimension_dist_radial(size_t *a, size_t *b, size_t ndim)
+{
+  size_t i, out=0;
+  for(i=0;i<ndim;++i) out += (a[i]-b[i])*(a[i]-b[i]);
+  return sqrt(out);
+}
+
+
+
+
+
 
 
 
diff --git a/lib/gnuastro-internal/commonopts.h 
b/lib/gnuastro-internal/commonopts.h
index 4420947..96046c2 100644
--- a/lib/gnuastro-internal/commonopts.h
+++ b/lib/gnuastro-internal/commonopts.h
@@ -199,6 +199,20 @@ struct argp_option gal_commonopts_options[] =
       GAL_OPTIONS_NOT_SET
     },
     {
+      "interpmetric",
+      GAL_OPTIONS_KEY_INTERPMETRIC,
+      "INT",
+      0,
+      "Interpolation metric (radial, manhattan).",
+      GAL_OPTIONS_GROUP_TESSELLATION,
+      &cp->interpmetric,
+      GAL_TYPE_STRING,
+      GAL_OPTIONS_RANGE_ANY,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET,
+      gal_options_read_interpmetric
+    },
+    {
       "interpnumngb",
       GAL_OPTIONS_KEY_INTERPNUMNGB,
       "INT",
diff --git a/lib/gnuastro-internal/options.h b/lib/gnuastro-internal/options.h
index ba22bd3..2584ea9 100644
--- a/lib/gnuastro-internal/options.h
+++ b/lib/gnuastro-internal/options.h
@@ -121,6 +121,7 @@ enum options_common_keys
   GAL_OPTIONS_KEY_CHECKTILES,
   GAL_OPTIONS_KEY_ONEELEMPERTILE,
   GAL_OPTIONS_KEY_INTERPONLYBLANK,
+  GAL_OPTIONS_KEY_INTERPMETRIC,
   GAL_OPTIONS_KEY_INTERPNUMNGB,
 };
 
@@ -179,6 +180,7 @@ struct gal_options_common_params
   /* Tessellation. */
   struct gal_tile_two_layer_params tl; /* Two layer tessellation params.  */
   uint8_t      interponlyblank; /* Only interpolate over blank values.    */
+  uint8_t         interpmetric; /* Metric to use for nearest-ngb interp.  */
   size_t          interpnumngb; /* Number of neighbors for interpolation. */
 
   /* Input. */
@@ -274,6 +276,10 @@ void *
 gal_options_read_tableformat(struct argp_option *option, char *arg,
                              char *filename, size_t lineno, void *junk);
 
+void *
+gal_options_read_interpmetric(struct argp_option *option, char *arg,
+                              char *filename, size_t lineno, void *junk);
+
 gal_data_t *
 gal_options_parse_list_of_numbers(char *string, char *filename,
                                   size_t lineno);
diff --git a/lib/gnuastro/dimension.h b/lib/gnuastro/dimension.h
index 205820a..2bdc9c4 100644
--- a/lib/gnuastro/dimension.h
+++ b/lib/gnuastro/dimension.h
@@ -89,10 +89,11 @@ gal_dimension_index_to_coord(size_t index, size_t ndim, 
size_t *dsize,
 /************************************************************************/
 /********************           Distances          **********************/
 /************************************************************************/
-size_t
+float
 gal_dimension_dist_manhattan(size_t *a, size_t *b, size_t ndim);
 
-
+float
+gal_dimension_dist_radial(size_t *a, size_t *b, size_t ndim);
 
 
 
diff --git a/lib/gnuastro/interpolate.h b/lib/gnuastro/interpolate.h
index 1c72730..5f1ec26 100644
--- a/lib/gnuastro/interpolate.h
+++ b/lib/gnuastro/interpolate.h
@@ -46,6 +46,18 @@ __BEGIN_C_DECLS  /* From C++ preparations */
 
 
 
+
+/* Metrics to use for nearest-neighbor  */
+enum gal_interpolate_close_metric
+{
+ GAL_INTERPOLATE_CLOSE_METRIC_INVALID,
+
+ GAL_INTERPOLATE_CLOSE_METRIC_RADIAL,
+ GAL_INTERPOLATE_CLOSE_METRIC_MANHATTAN,
+};
+
+
+
 /* Types of interpolation. */
 enum gal_interpolate_1D_types
 {
@@ -65,8 +77,9 @@ enum gal_interpolate_1D_types
 gal_data_t *
 gal_interpolate_close_neighbors(gal_data_t *input,
                                 struct gal_tile_two_layer_params *tl,
-                                size_t numneighbors, size_t numthreads,
-                                int onlyblank, int aslinkedlist);
+                                uint8_t metric, size_t numneighbors,
+                                size_t numthreads, int onlyblank,
+                                int aslinkedlist);
 
 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 8277afc..bdc4f50 100644
--- a/lib/interpolate.c
+++ b/lib/interpolate.c
@@ -67,6 +67,8 @@ struct interpolate_params
   uint8_t                *thread_flags;
   int                        onlyblank;
   gal_list_void_t            *ngb_vals;
+  float (*metric)(size_t *, size_t *, size_t );
+
   struct gal_tile_two_layer_params *tl;
 };
 
@@ -87,11 +89,11 @@ interpolate_close_neighbors_on_thread(void *in_prm)
 
   /* Rest of variables. */
   void *nv;
-  float pdist;
+  float dist, pdist;
   uint8_t *b, *bf, *bb;
   gal_list_void_t *tvll;
   gal_list_dosizet_t *lQ, *sQ;
-  size_t ngb_counter, dist, pind, *dinc;
+  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;
   size_t size = (correct_index ? tl->tottilesinch : input->size);
@@ -233,7 +235,7 @@ interpolate_close_neighbors_on_thread(void *in_prm)
                  gal_dimension_index_to_coord(nind, ndim, dsize, ncoord);
 
                  /* Distance of this neighbor to the one to be filled. */
-                 dist=gal_dimension_dist_manhattan(icoord, ncoord, ndim);
+                 dist=prm->metric(icoord, ncoord, ndim);
 
                  /* Add this neighbor to the list. */
                  gal_list_dosizet_add(&lQ, &sQ, nind, dist);
@@ -328,8 +330,9 @@ interpolate_copy_input(gal_data_t *input, int aslinkedlist)
 gal_data_t *
 gal_interpolate_close_neighbors(gal_data_t *input,
                                 struct gal_tile_two_layer_params *tl,
-                                size_t numneighbors, size_t numthreads,
-                                int onlyblank, int aslinkedlist)
+                                uint8_t metric, size_t numneighbors,
+                                size_t numthreads, int onlyblank,
+                                int aslinkedlist)
 {
   gal_data_t *tin, *tout;
   struct interpolate_params prm;
@@ -354,6 +357,21 @@ gal_interpolate_close_neighbors(gal_data_t *input,
   prm.num          = aslinkedlist ? gal_list_data_number(input) : 1;
 
 
+  /* Set the metric. */
+  switch(metric)
+    {
+    case GAL_INTERPOLATE_CLOSE_METRIC_RADIAL:
+      prm.metric=gal_dimension_dist_radial;
+      break;
+    case GAL_INTERPOLATE_CLOSE_METRIC_MANHATTAN:
+      prm.metric=gal_dimension_dist_manhattan;
+      break;
+    default:
+      error(EXIT_FAILURE, 0, "%s: %d is not a valid metric identifier",
+            __func__, metric);
+    }
+
+
   /* Flag the blank values. */
   prm.blanks=gal_blank_flag(input);
 
diff --git a/lib/options.c b/lib/options.c
index db728ff..2e176ed 100644
--- a/lib/options.c
+++ b/lib/options.c
@@ -38,6 +38,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/threads.h>
 #include <gnuastro/pointer.h>
 #include <gnuastro/arithmetic.h>
+#include <gnuastro/interpolate.h>
 
 #include <gnuastro-internal/timing.h>
 #include <gnuastro-internal/options.h>
@@ -511,6 +512,55 @@ gal_options_read_tableformat(struct argp_option *option, 
char *arg,
 
 
 
+void *
+gal_options_read_interpmetric(struct argp_option *option, char *arg,
+                              char *filename, size_t lineno, void *junk)
+{
+  char *str;
+  if(lineno==-1)
+    {
+      switch(*(uint8_t *)(option->value))
+        {
+        case GAL_INTERPOLATE_CLOSE_METRIC_RADIAL:
+          gal_checkset_allocate_copy("radial", &str);
+          break;
+        case GAL_INTERPOLATE_CLOSE_METRIC_MANHATTAN:
+          gal_checkset_allocate_copy("manhattan", &str);
+          break;
+        default:
+          error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix 
the "
+                "problem. The code %u is not recognized as a nearest-neighbor "
+                "interpolation metric", __func__, PACKAGE_BUGREPORT,
+                *(uint8_t *)(option->value));
+        }
+      return str;
+    }
+  else
+    {
+      /* If the option is already set, just return. */
+      if(option->set) return NULL;
+
+      /* Set the value. */
+      if(       !strcmp(arg, "radial") )
+        *(uint8_t *)(option->value) = GAL_INTERPOLATE_CLOSE_METRIC_RADIAL;
+      else if ( !strcmp(arg, "manhattan") )
+        *(uint8_t *)(option->value) = GAL_INTERPOLATE_CLOSE_METRIC_MANHATTAN;
+      else
+        error_at_line(EXIT_FAILURE, 0, filename, lineno, "`%s' (value to `%s' "
+                      "option) isn't valid. Currently only `radial' and "
+                      "`manhattan' metrics are recognized for nearest neighbor 
"
+                      "interpolation", arg, option->name);
+
+      /* For no un-used variable warning. This function doesn't need the
+         pointer. */
+      return junk=NULL;
+    }
+}
+
+
+
+
+
 /* The input to this function is a string of any number of numbers
    separated by a comma (`,') and possibly containing fractions, for
    example: `1,2/3, 4.95'. The output `gal_data_t' contains the array of