[gnuastro-commits] master b03f7f9 113/125: New inplementation of interpolation complete

[Mohammad Akhlaghi]

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

    New inplementation of interpolation complete
    
    An initial working implementation of interpolation with the new approach is
    now complete. In the previous method, interpolation was heavily intertwined
    with the tile structure, but now it can be done completely independent of
    it. However, it is possible for the user to give a
    `gal_tile_two_layer_params' structure to benefit from (for example to only
    interpolate over a channel).
    
    Other updates in this commit:
    
     - Also in this commit, the old `gal_tile_function_on_threads' function was
       changed to the more generic `gal_threads_spin_off'. In practice, the old
       function just used the number of tiles, but the new version takes in the
       number of actions as one argument, so it is usable in any context.
    
     - The new `gal_dimension_dist_manhattan' function will calculate the
       manhattan distance of two coordinates.
---
 bin/statistics/args.h             |  26 ++++
 bin/statistics/aststatistics.conf |   4 +
 bin/statistics/main.h             |   2 +
 bin/statistics/statistics.c       |  13 +-
 bin/statistics/ui.c               |   2 +
 bin/statistics/ui.h               |   2 +
 lib/convolve.c                    |   6 +-
 lib/dimension.c                   |  30 ++++
 lib/gnuastro/dimension.h          |  10 ++
 lib/gnuastro/interpolate.h        |  16 +-
 lib/gnuastro/threads.h            |  17 +++
 lib/gnuastro/tile.h               |  22 ---
 lib/interpolate.c                 | 311 +++++++++++++++++++++++++++++++++++---
 lib/threads.c                     | 168 ++++++++++++++++++++
 lib/tile.c                        | 156 -------------------
 15 files changed, 577 insertions(+), 208 deletions(-)

diff --git a/bin/statistics/args.h b/bin/statistics/args.h
index 40eb667..01b8b72 100644
--- a/bin/statistics/args.h
+++ b/bin/statistics/args.h
@@ -114,6 +114,32 @@ struct argp_option program_options[] =
       GAL_OPTIONS_NOT_MANDATORY,
       GAL_OPTIONS_NOT_SET
     },
+    {
+      "interponlyblank",
+      ARGS_OPTION_KEY_INTERPONLYBLANK,
+      0,
+      0,
+      "Only interpolate over the blank values.",
+      GAL_OPTIONS_GROUP_TESSELLATION,
+      &p->interponlyblank,
+      GAL_OPTIONS_NO_ARG_TYPE,
+      GAL_OPTIONS_RANGE_0_OR_1,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
+    },
+    {
+      "interpnumngb",
+      ARGS_OPTION_KEY_INTERPNUMNGB,
+      0,
+      0,
+      "Number of neighbors to use for interpolation.",
+      GAL_OPTIONS_GROUP_TESSELLATION,
+      &p->interpnumngb,
+      GAL_OPTIONS_NO_ARG_TYPE,
+      GAL_OPTIONS_RANGE_GT_0,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
+    },
 
 
 
diff --git a/bin/statistics/aststatistics.conf 
b/bin/statistics/aststatistics.conf
index 25810e3..32db56a 100644
--- a/bin/statistics/aststatistics.conf
+++ b/bin/statistics/aststatistics.conf
@@ -19,6 +19,10 @@
 
 # Input image:
 
+# Tessellation
+ interponlyblank      1
+ interpnumngb         9
+
 # Histogram and CFP settings
  numasciibins        70
  asciiheight         10
diff --git a/bin/statistics/main.h b/bin/statistics/main.h
index 3c449c1..6e8cfd8 100644
--- a/bin/statistics/main.h
+++ b/bin/statistics/main.h
@@ -66,6 +66,8 @@ struct statisticsparams
   float           quantmax;  /* Quantile maximum.                        */
   uint8_t           ontile;  /* Do single value calculations on tiles.   */
   uint8_t      interpolate;  /* Use interpolation to fill blank tiles.   */
+  uint8_t  interponlyblank;  /* Only interpolate over blank values.      */
+  size_t      interpnumngb;  /* Number of neighbors for interpolation.   */
 
   uint8_t        asciihist;  /* Print an ASCII histogram.                */
   uint8_t         asciicfp;  /* Print an ASCII cumulative frequency plot.*/
diff --git a/bin/statistics/statistics.c b/bin/statistics/statistics.c
index 6a104f3..a31c6e7 100644
--- a/bin/statistics/statistics.c
+++ b/bin/statistics/statistics.c
@@ -237,12 +237,21 @@ static void
 statistics_interpolate_and_write(struct statisticsparams *p,
                                  gal_data_t *values)
 {
+  gal_data_t *interpd;
   char *output=p->cp.output;
   struct gal_tile_two_layer_params *tl=&p->cp.tl;
 
+
   /* Do the interpolation (if necessary). */
-  if(p->interpolate && gal_blank_present(values))
-    gal_tile_full_interpolate(values, tl);
+  if( p->interpolate
+      && !(p->interponlyblank && gal_blank_present(values)==0) )
+    {
+      interpd=gal_interpolate_close_neighbors(values, tl, p->interpnumngb,
+                                              p->cp.numthreads,
+                                              p->interponlyblank);
+      gal_data_free(values);
+      values=interpd;
+    }
 
   /* Write the values. */
   gal_tile_full_values_write(values, tl, output, PROGRAM_STRING);
diff --git a/bin/statistics/ui.c b/bin/statistics/ui.c
index 56f7199..282c56f 100644
--- a/bin/statistics/ui.c
+++ b/bin/statistics/ui.c
@@ -33,6 +33,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/qsort.h>
 #include <gnuastro/blank.h>
 #include <gnuastro/table.h>
+#include <gnuastro/threads.h>
 #include <gnuastro/arithmetic.h>
 #include <gnuastro/linkedlist.h>
 #include <gnuastro/statistics.h>
@@ -126,6 +127,7 @@ ui_initialize_options(struct statisticsparams *p,
   cp->program_bibtex     = PROGRAM_BIBTEX;
   cp->program_authors    = PROGRAM_AUTHORS;
   cp->coptions           = gal_commonopts_options;
+  cp->numthreads         = gal_threads_number();
 
   /* Program-specific initializers */
   p->lessthan            = NAN;
diff --git a/bin/statistics/ui.h b/bin/statistics/ui.h
index b681f4f..2ab02b8 100644
--- a/bin/statistics/ui.h
+++ b/bin/statistics/ui.h
@@ -72,6 +72,8 @@ enum option_keys_enum
   ARGS_OPTION_KEY_ONEBINSTART,
   ARGS_OPTION_KEY_MAXBINONE,
   ARGS_OPTION_KEY_MIRRORDIST,
+  ARGS_OPTION_KEY_INTERPONLYBLANK,
+  ARGS_OPTION_KEY_INTERPNUMNGB,
 };
 
 
diff --git a/lib/convolve.c b/lib/convolve.c
index e3b1e3d..880704a 100644
--- a/lib/convolve.c
+++ b/lib/convolve.c
@@ -388,7 +388,7 @@ convolve_spatial_tile(struct spatial_params *cprm, size_t 
tile_ind,
 static void *
 convolve_spatial_on_thread(void *inparam)
 {
-  struct gal_tile_thread_param *tprm=(struct gal_tile_thread_param *)inparam;
+  struct gal_threads_params *tprm=(struct gal_threads_params *)inparam;
   struct spatial_params *cprm=(struct spatial_params *)(tprm->params);
 
   size_t i;
@@ -450,8 +450,8 @@ gal_convolve_spatial(gal_data_t *tiles, gal_data_t *kernel,
   params.edgecorrection=edgecorrection;
 
   /* Do the spatial convolution on threads. */
-  gal_tile_function_on_threads(tiles, convolve_spatial_on_thread,
-                               numthreads, &params);
+  gal_threads_spin_off(convolve_spatial_on_thread, &params,
+                       gal_data_num_in_ll(tiles), numthreads);
 
   /* Clean up and return the output array. */
   return out;
diff --git a/lib/dimension.c b/lib/dimension.c
index edc6672..cced659 100644
--- a/lib/dimension.c
+++ b/lib/dimension.c
@@ -212,3 +212,33 @@ gal_dimension_index_to_coord(size_t ind, size_t ndim, 
size_t *dsize,
       free(dinc);
     }
 }
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/************************************************************************/
+/********************           Distances          **********************/
+/************************************************************************/
+size_t
+gal_dimension_dist_manhattan(size_t *a, size_t *b, size_t ndim)
+{
+  size_t i, out=0;
+  for(i=0;i<ndim;++i) out += labs(a[i]-b[i]);
+  return out;
+}
diff --git a/lib/gnuastro/dimension.h b/lib/gnuastro/dimension.h
index 60a0c8b..95e6908 100644
--- a/lib/gnuastro/dimension.h
+++ b/lib/gnuastro/dimension.h
@@ -80,6 +80,16 @@ gal_dimension_index_to_coord(size_t ind, size_t ndim, size_t 
*dsize,
 
 
 /************************************************************************/
+/********************           Distances          **********************/
+/************************************************************************/
+size_t
+gal_dimension_dist_manhattan(size_t *a, size_t *b, size_t ndim);
+
+
+
+
+
+/************************************************************************/
 /********************          Neighbors           **********************/
 /************************************************************************/
 /* Purpose
diff --git a/lib/gnuastro/interpolate.h b/lib/gnuastro/interpolate.h
index 3c2495a..30aab88 100644
--- a/lib/gnuastro/interpolate.h
+++ b/lib/gnuastro/interpolate.h
@@ -26,6 +26,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 /* Include other headers if necessary here. Note that other header files
    must be included before the C++ preparations below */
 #include <gnuastro/data.h>
+#include <gnuastro/tile.h>
 
 /* C++ Preparations */
 #undef __BEGIN_C_DECLS
@@ -42,8 +43,19 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 /* Actual header contants (the above were for the Pre-processor). */
 __BEGIN_C_DECLS  /* From C++ preparations */
 
-void
-gal_interpolate(gal_data_t *input);
+
+
+
+
+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);
+
+
+
+
 
 __END_C_DECLS    /* From C++ preparations */
 
diff --git a/lib/gnuastro/threads.h b/lib/gnuastro/threads.h
index 6c0b912..5943f01 100644
--- a/lib/gnuastro/threads.h
+++ b/lib/gnuastro/threads.h
@@ -114,6 +114,23 @@ gal_threads_attr_barrier_init(pthread_attr_t *attr, 
pthread_barrier_t *b,
 
 
 
+
+/*******************************************************************/
+/************     Run a function on multiple threads  **************/
+/*******************************************************************/
+struct gal_threads_params
+{
+  size_t            id; /* Id of this thread.                            */
+  void         *params; /* Input structure for higher-level settings.    */
+  size_t       *indexs; /* Indexes of actions to be done in this thread. */
+  pthread_barrier_t *b; /* Pointer the barrier for all threads.          */
+};
+
+void
+gal_threads_spin_off(void *(*function)(void *), void *caller_params,
+                     size_t numactions, size_t numthreads);
+
+
 __END_C_DECLS    /* From C++ preparations */
 
 #endif           /* __GAL_THREADS_H__ */
diff --git a/lib/gnuastro/tile.h b/lib/gnuastro/tile.h
index f4c3420..766e08a 100644
--- a/lib/gnuastro/tile.h
+++ b/lib/gnuastro/tile.h
@@ -137,28 +137,6 @@ gal_tile_full_values_write(gal_data_t *tilevalues,
                            struct gal_tile_two_layer_params *tl,
                            char *filename, char *program_string);
 
-void
-gal_tile_full_interpolate(gal_data_t *tilevalues,
-                          struct gal_tile_two_layer_params *tl);
-
-
-
-
-
-/*********************************************************************/
-/********************         On threads          ********************/
-/*********************************************************************/
-struct gal_tile_thread_param
-{
-  size_t            id; /* Id of this thread.                            */
-  void         *params; /* Input structure for higher-level settings.    */
-  size_t       *indexs; /* Indexes of actions to be done in this thread. */
-  pthread_barrier_t *b; /* Pointer the barrier for all threads.          */
-};
-
-void
-gal_tile_function_on_threads(gal_data_t *tiles, void *(*meshfunc)(void *),
-                             size_t numthreads, void *caller_params);
 
 
 
diff --git a/lib/interpolate.c b/lib/interpolate.c
index 0723d46..5ae15da 100644
--- a/lib/interpolate.c
+++ b/lib/interpolate.c
@@ -23,39 +23,304 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <config.h>
 
 #include <stdio.h>
+#include <errno.h>
+#include <error.h>
 #include <stdlib.h>
+#include <string.h>
 
 #include <gnuastro/data.h>
 #include <gnuastro/fits.h>
 #include <gnuastro/blank.h>
+#include <gnuastro/threads.h>
 #include <gnuastro/dimension.h>
+#include <gnuastro/statistics.h>
+#include <gnuastro/interpolate.h>
+#include <gnuastro/permutation.h>
 
 
 
 
 
-/* Interpolate blank values in an array. */
-void
-gal_interpolate(gal_data_t *input)
+
+/*********************************************************************/
+/********************      Nearest neighbor       ********************/
+/*********************************************************************/
+/* We want the flags to be powers of two so we can use bit-wise checks. */
+#define INTERPOLATE_FLAGS_NO       0
+#define INTERPOLATE_FLAGS_CHECKED  1<<0
+#define INTERPOLATE_FLAGS_BLANK    1<<1
+
+
+
+
+
+/* Paramters for interpolation on threads. */
+struct interpolate_params
 {
-  size_t F;
-  uint8_t *flagarr;
-  gal_data_t *flag;
-  size_t *dinc=gal_dimension_increment(input->ndim, input->dsize);
-
-  /* Set a value of 1 for every element that must be interpolated. */
-  flag=gal_blank_flag(input);
-
-  /* Find the proper value for each neighbor. */
-  flagarr=flag->array;
-  for(F=0; F<input->size; ++F)
-    if(flagarr[F])
-      {
-        printf("to be filled: %zu\n", F);
-        GAL_DIMENSION_NEIGHBOR_OP(F, input->ndim, input->dsize, 2,
-                                  dinc, {printf("\tneighbor: %zu\n", nind);});
-      }
-
-  /* Clean up. */
-  gal_data_free(flag);
+  gal_data_t                    *input;
+  gal_data_t                      *out;
+  gal_data_t                   *blanks;
+  size_t                  numneighbors;
+  uint8_t                *thread_flags;
+  void                       *ngb_vals;
+  int                        onlyblank;
+  struct gal_tile_two_layer_params *tl;
+};
+
+
+
+
+
+/* Run the interpolation on many threads. */
+static void *
+interpolate_close_neighbors_on_thread(void *in_prm)
+{
+  /* 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 gal_tile_two_layer_params *tl=prm->tl;
+  int correct_index=(tl && tl->totchannels>1 && tl->workoverch==0);
+  gal_data_t *input=prm->input;
+
+  /* Rest of variables. */
+  float pdist;
+  gal_data_t *nearest, *single;
+  size_t ngb_counter, dist, *dinc;
+  struct gal_linkedlist_tosll *lQ, *sQ;
+  uint8_t *b, *bf, *bb, type=input->type;
+  void *values=input->array, *out=prm->out->array;
+  size_t pind, twidth=gal_data_sizeof(input->type);
+  size_t i, index, fullind, chstart, ndim=input->ndim;
+  size_t size = (correct_index ? tl->tottilesinch : input->size);
+  size_t *icoord=gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, ndim);
+  size_t *ncoord=gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, ndim);
+  size_t *dsize = (correct_index ? tl->numtilesinch : input->dsize);
+  uint8_t *fullflag=&prm->thread_flags[tprm->id*input->size], *flag=fullflag;
+  void *nv=gal_data_ptr_increment(prm->ngb_vals, tprm->id*prm->numneighbors,
+                                  input->type);
+
+  /* Initializations. */
+  bb=prm->blanks->array;
+  bf=(b=fullflag)+input->size;
+  dinc=gal_dimension_increment(ndim, dsize);
+  do *b = *bb++ ? INTERPOLATE_FLAGS_BLANK : 0; while(++b<bf);
+
+
+  /* Put the nearest neighbor values into a structure for easy processing
+     later. */
+  nearest=gal_data_alloc(nv, input->type, 1, &prm->numneighbors, NULL, 0,
+                         0, NULL, NULL, NULL);
+
+
+  /* Go over all the points given to this thread. */
+  for(i=0; tprm->indexs[i] != GAL_THREADS_NON_THRD_INDEX; ++i)
+    {
+      /* For easy reading. */
+      fullind=tprm->indexs[i];
+
+
+      /* If the caller only wanted to interpolate over blank values and
+         this value is not blank (we know from the flags), then just set
+         the output value at this element to the input value and go to the
+         next element. */
+      if(prm->onlyblank && !(fullflag[fullind] & INTERPOLATE_FLAGS_BLANK) )
+        {
+          /* It should be `input->array', not `values'! Because when
+             channels are treated separately, `values' is going to
+             change. */
+          memcpy(gal_data_ptr_increment(out,          fullind, type),
+                 gal_data_ptr_increment(input->array, fullind, type),
+                 twidth);
+          continue;
+        }
+
+
+      /* Correct the index (if necessary). When the values come from a
+         tiled dataset, the caller might want to interpolate the values of
+         each channel separately (not mix values from different
+         channels). In such a case, the tiles of each channel (and their
+         values in `input' are contiguous. So we need to correct
+         `tprm->indexs[i]' (which is the index over the whole tessellation,
+         including all channels). */
+      if(correct_index)
+        {
+          /* Index of this tile in its channel. */
+          index = fullind % tl->tottilesinch;
+
+          /* Index of the first tile in this channel. */
+          chstart = (fullind / tl->tottilesinch) * tl->tottilesinch;
+
+          /* Correct the values and flag pointers so we can only work in
+             this channel. */
+          values = gal_data_ptr_increment(input->array, chstart, type);
+          flag = gal_data_ptr_increment(fullflag, chstart,
+                                        GAL_DATA_TYPE_UINT8);
+        }
+      else
+        index=fullind;
+
+
+      /* Reset all checked bits in the flags array to 0. */
+      ngb_counter=0;
+      bf=(b=flag)+size;
+      do *b &= ~(INTERPOLATE_FLAGS_CHECKED); while(++b<bf);
+
+
+      /* Get the coordinates of this pixel (to be interpolated). */
+      gal_dimension_index_to_coord(index, ndim, dsize, icoord);
+
+
+      /* Start parsing the neighbors. We will use a two-way ordered linked
+         list structure. To start from the nearest and go out to the
+         farthest. */
+      lQ=sQ=NULL;
+      gal_linkedlist_add_to_tosll_end(&lQ, &sQ, index, 0.0f);
+      while(sQ)
+        {
+          /* Pop-out (p) an index from the queue: */
+          gal_linkedlist_pop_from_tosll_start(&lQ, &sQ, &pind, &pdist);
+
+          /* If this isn't a blank value (recall that `index' was the first
+             element added to the list which might be blank). */
+          if( !(flag[pind] & INTERPOLATE_FLAGS_BLANK) )
+            {
+              /* Copy the value into the `nv' array. */
+              memcpy(gal_data_ptr_increment(nv, ngb_counter, type),
+                     gal_data_ptr_increment(values, pind, type),
+                     twidth);
+
+              /* If we have filled all the elements, break out. */
+              if(++ngb_counter>=prm->numneighbors) break;
+            }
+
+          /* Go over all the neighbors of this popped pixel and add them to
+             the list of neighbors to be checked. */
+          GAL_DIMENSION_NEIGHBOR_OP(pind, ndim, dsize, 1, dinc,
+           {
+             /* Only look at neighbors that have not been checked. VERY
+                IMPORTANT: we must not check for blank values here,
+                otherwise we won't be able to parse over extended blank
+                regions. */
+             if( !(flag[nind] & INTERPOLATE_FLAGS_CHECKED) )
+               {
+                 /* Get the coordinates of this neighbor. */
+                 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);
+
+                 /* Add this neighbor to the list. */
+                 gal_linkedlist_add_to_tosll_end(&lQ, &sQ, nind, dist);
+
+                 /* Flag this neighbor as checked. */
+                 flag[nind] |= INTERPOLATE_FLAGS_CHECKED;
+               }
+           } );
+
+          /* If there are no more meshes to add to the queue, then this
+             shows, there were not enough points for
+             interpolation. Normally, this loop should only be exited
+             through the `currentnum>=numnearest' check above. */
+          if(sQ==NULL)
+            error(EXIT_FAILURE, 0, "only %zu neighbors found while you had "
+                  "asked to use %zu neighbors for close neighbor "
+                  "interpolation", ngb_counter, prm->numneighbors);
+        }
+
+      /* Calculate the median of the values and write it in. */
+      single=gal_statistics_median(nearest, 1);
+      memcpy(gal_data_ptr_increment(out, fullind, type), single->array,
+             twidth);
+
+      /* Clean up. */
+      gal_data_free(single);
+    }
+
+
+  /* Clean up, wait for all the other threads to finish and return. */
+  free(icoord);
+  free(ncoord);
+  nearest->array=NULL;
+  gal_data_free(nearest);
+  if(tprm->b) pthread_barrier_wait(tprm->b);
+  return NULL;
+}
+
+
+
+
+
+/* Interpolate blank values in an array. If the `tl!=NULL', then it is
+   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,
+                                size_t numneighbors, size_t numthreads,
+                                int onlyblank)
+{
+  struct interpolate_params prm;
+  int permute=(tl && tl->totchannels>1 && tl->workoverch);
+
+
+  /* Initialize the constant parameters. */
+  prm.tl=tl;
+  prm.input=input;
+  prm.onlyblank=onlyblank;
+  prm.numneighbors=numneighbors;
+
+
+  /* Flag the blank values. */
+  prm.blanks=gal_blank_flag(input);
+
+
+  /* Allocate space for the output. */
+  prm.out=gal_data_alloc(NULL, input->type, input->ndim, input->dsize,
+                         input->wcs, 0, input->minmapsize, "INTERPOLATED",
+                         input->unit, NULL);
+
+
+  /* If the input is from a tile structure and the user has asked to ignore
+     channels, then re-order the values. */
+  if(permute)
+    {
+      /* Prepare the permutation (if necessary/not already defined). */
+      gal_tile_full_permutation(tl);
+
+      /* Re-order values to ignore channels (if necessary). */
+      gal_permutation_apply(input, tl->permutation);
+      gal_permutation_apply(prm.blanks, tl->permutation);
+    }
+
+
+  /* Allocate space for all the flag values of all the threads here (memory
+     in each thread is limited) and this is cleaner. */
+  prm.thread_flags=gal_data_malloc_array(GAL_DATA_TYPE_UINT8,
+                                         numthreads*input->size);
+
+
+  /* Allocate space for the list of neighbor values in each thread. */
+  prm.ngb_vals=gal_data_malloc_array(input->type, numthreads*numneighbors);
+
+
+  /* Spin off the threads. */
+  gal_threads_spin_off(interpolate_close_neighbors_on_thread, &prm,
+                       input->size, numthreads);
+
+
+  /* If the values were permuted for the interpolation, then re-order the
+     values back to their original location (so they correspond to their
+     tile indexs. */
+  if(permute)
+    {
+      gal_permutation_apply_inverse(input, tl->permutation);
+      gal_permutation_apply_inverse(prm.out, tl->permutation);
+    }
+
+
+  /* Clean up and return. */
+  gal_data_free(prm.blanks);
+  free(prm.thread_flags);
+  free(prm.ngb_vals);
+  return prm.out;
 }
diff --git a/lib/threads.c b/lib/threads.c
index 099312e..c1bfa46 100644
--- a/lib/threads.c
+++ b/lib/threads.c
@@ -236,3 +236,171 @@ gal_threads_attr_barrier_init(pthread_attr_t *attr, 
pthread_barrier_t *b,
   err=pthread_barrier_init(b, NULL, limit);
   if(err) error(EXIT_FAILURE, 0, "thread barrier not initialized");
 }
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/*******************************************************************/
+/************     Run a function on multiple threads  **************/
+/*******************************************************************/
+/* Run a given function on the given tiles. The function has to be
+   link-able with your final executable and has to have only one `void *'
+   argument and return a `void *' value. To have access to
+   variables/parameters in the function, you have to define a structure and
+   pass its pointer as `caller_params'.
+
+   Here is one simple example. At least two functions and one structure are
+   necessary to use this function.
+
+     --------- Parameters to keep values you need ---------
+     struct my_params
+     {
+       int    value1;
+       double value2;
+       float  *array;
+     }
+
+
+     ---------    Function to run on each thread  ---------
+     void *
+     run_on_thread(void *in_prm)
+     {
+       size_t i;
+       struct gal_threads_params *tprm=(struct gal_threads_params *)in_prm;
+       struct my_params *prm=(struct my_params *)(tprm->params);
+
+       for(i=0; tprm->indexs[i] != GAL_THREADS_NON_THRD_INDEX; ++i)
+       {
+
+           THE INDEX OF THE TARGET IS NOW AVAILABLE AS
+           `tprm->indexs[i]'. YOU CAN USE IT IN WHAT EVER MANNER YOU LIKE
+           ALONG WITH THE SET OF VARIABLES/ARRAYS in `prm'.
+
+       }
+
+       if(tprm->b) pthread_barrier_wait(tprm->b);
+       return NULL;
+     }
+
+
+     ---------         High-level function        ---------
+     int
+     higher_level_function(float *array, size_t num_in_array, int value1)
+     {
+        double value2;
+        struct my_params;
+        size_t numthreads;
+
+        my_params.value1=value1;
+        my_params.value2=value2;
+        my_params.arary=array;
+
+        gal_threads_spin_off(run_on_thread, &my_params, num_in_array,
+                             numthreads);
+
+        return 1;
+     }
+
+  For real world applications of this function, you can also inspect the
+  Gnuastro source. There are also many cases in Gnuastro where we benefit
+  from this function. Please run the following command from the top source
+  directory of Gnuastro to see where:
+
+      $ grep -r gal_threads_spin_off ./
+*/
+void
+gal_threads_spin_off(void *(*function)(void *), void *caller_params,
+                     size_t numactions, size_t numthreads)
+{
+  int err;
+  pthread_t t;          /* All thread ids saved in this, not used. */
+  pthread_attr_t attr;
+  pthread_barrier_t b;
+  struct gal_threads_params *prm;
+  size_t i, *indexs, thrdcols, numbarriers;
+
+  /* If there are no actions, then just return. */
+  if(numactions==0) return;
+
+  /* Sanity check. */
+  if(numthreads==0)
+    error(EXIT_FAILURE, 0, "the number of threads (`numthreads') in "
+          "`gal_threads_spin_off' cannot be zero");
+
+  /* Allocate the array of parameters structure structures. */
+  errno=0;
+  prm=malloc(numthreads*sizeof *prm);
+  if(prm==NULL)
+    {
+      fprintf(stderr, "%zu bytes could not be allocated for prm.",
+              numthreads*sizeof *prm);
+      exit(EXIT_FAILURE);
+    }
+
+  /* Distribute the actions into the threads: */
+  gal_threads_dist_in_threads(numactions, numthreads, &indexs, &thrdcols);
+
+  /* Do the job: when only one thread is necessary, there is no need to
+     spin off one thread, just call the function directly (spinning off
+     threads is expensive). This is for the generic thread spinner
+     function, not this simple function where `numthreads' is a
+     constant. */
+  if(numthreads==1)
+    {
+      prm[0].id=0;
+      prm[0].b=NULL;
+      prm[0].indexs=indexs;
+      prm[0].params=caller_params;
+      function(&prm[0]);
+    }
+  else
+    {
+      /* Initialize the attributes. Note that this running thread
+         (that spinns off the nt threads) is also a thread, so the
+         number the barriers should be one more than the number of
+         threads spinned off. */
+      numbarriers = (numactions<numthreads ? numactions : numthreads) + 1;
+      gal_threads_attr_barrier_init(&attr, &b, numbarriers);
+
+      /* Spin off the threads: */
+      for(i=0;i<numthreads;++i)
+        if(indexs[i*thrdcols]!=GAL_THREADS_NON_THRD_INDEX)
+          {
+            prm[i].id=i;
+            prm[i].b=&b;
+            prm[i].params=caller_params;
+            prm[i].indexs=&indexs[i*thrdcols];
+            err=pthread_create(&t, &attr, function, &prm[i]);
+            if(err)
+              {
+                fprintf(stderr, "can't create thread %zu", i);
+                exit(EXIT_FAILURE);
+              }
+          }
+
+      /* Wait for all threads to finish and free the spaces. */
+      pthread_barrier_wait(&b);
+      pthread_attr_destroy(&attr);
+      pthread_barrier_destroy(&b);
+    }
+
+  /* Clean up. */
+  free(prm);
+  free(indexs);
+}
diff --git a/lib/tile.c b/lib/tile.c
index 9dc2531..330fe02 100644
--- a/lib/tile.c
+++ b/lib/tile.c
@@ -970,68 +970,6 @@ gal_tile_full_values_write(gal_data_t *tilevalues,
 
 
 
-/* Interpolate the blank tile values while respecting the channels (if
-   requested). The output maybe a newly allocated array (if there are
-   channels), or the same array as input with interpolated values.
-
-   IMPORTANT: it is assumed that the values are in the same order as the
-   tiles.
-
-                      tile[i]  -->   tilevalues[i]                    */
-void
-gal_tile_full_interpolate(gal_data_t *tilevalues,
-                          struct gal_tile_two_layer_params *tl)
-{
-  size_t i;
-  void *chvstart;
-  gal_data_t *tointerp;
-
-  /* Ignore channels (if requested). */
-  if(tl->workoverch)
-    {
-      /* Prepare the permutation (if necessary/not already defined). */
-      gal_tile_full_permutation(tl);
-
-      /* Re-order values to ignore channels (if necessary). */
-      gal_permutation_apply(tilevalues, tl->permutation);
-
-      /* Interpolate the values. */
-      gal_interpolate(tilevalues);
-
-      /* Re-order values back to their original order. */
-      gal_permutation_apply_inverse(tilevalues, tl->permutation);
-    }
-  else
-    {
-      /* Interpolate each channel individually. */
-      for(i=0;i<tl->totchannels;++i)
-        {
-          /* Set the starting pointer for this channel IN THE VALUES
-             array. */
-          chvstart=gal_data_ptr_increment(tilevalues->array,
-                                          i*tl->tottilesinch,
-                                          tilevalues->type);
-
-          /* Make a cover-data-strcuture for the values in this
-             channel. */
-          tointerp=gal_data_alloc(chvstart, tilevalues->type,
-                                  tilevalues->ndim, tl->numtilesinch,
-                                  NULL, 0, 0, NULL, NULL, NULL);
-
-          /* Interpolate over the tiles in this channel. */
-          gal_interpolate(tointerp);
-
-          /* Clean up. */
-          tointerp->array=NULL;
-          gal_data_free(tointerp);
-        }
-    }
-}
-
-
-
-
-
 /* Clean up the allocated spaces in the parameters. */
 void
 gal_tile_full_free_contents(struct gal_tile_two_layer_params *tl)
@@ -1049,97 +987,3 @@ gal_tile_full_free_contents(struct 
gal_tile_two_layer_params *tl)
   if(tl->tiles)    gal_data_array_free(tl->tiles,    tl->tottiles,    0);
   if(tl->channels) gal_data_array_free(tl->channels, tl->totchannels, 0);
 }
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-/*********************************************************************/
-/********************         On threads          ********************/
-/*********************************************************************/
-/* Run a given function on the given tiles. */
-void
-gal_tile_function_on_threads(gal_data_t *tiles, void *(*function)(void *),
-                             size_t numthreads, void *caller_params)
-{
-  int err;
-  pthread_t t;          /* All thread ids saved in this, not used. */
-  pthread_attr_t attr;
-  pthread_barrier_t b;
-  struct gal_tile_thread_param *prm;
-  size_t i, *indexs, thrdcols, numbarriers;
-  size_t numtiles=gal_data_num_in_ll(tiles);
-
-  /* Allocate the array of parameters structure structures. */
-  prm=malloc(numthreads*sizeof *prm);
-  if(prm==NULL)
-    {
-      fprintf(stderr, "%zu bytes could not be allocated for prm.",
-              numthreads*sizeof *prm);
-      exit(EXIT_FAILURE);
-    }
-
-  /* Distribute the actions into the threads: */
-  gal_threads_dist_in_threads(numtiles, numthreads, &indexs, &thrdcols);
-
-  /* Do the job: when only one thread is necessary, there is no need to
-     spin off one thread, just call the function directly (spinning off
-     threads is expensive). This is for the generic thread spinner
-     function, not this simple function where `numthreads' is a
-     constant. */
-  if(numthreads==1)
-    {
-      prm[0].id=0;
-      prm[0].b=NULL;
-      prm[0].indexs=indexs;
-      prm[0].params=caller_params;
-      function(&prm[0]);
-    }
-  else
-    {
-      /* Initialize the attributes. Note that this running thread
-         (that spinns off the nt threads) is also a thread, so the
-         number the barriers should be one more than the number of
-         threads spinned off. */
-      numbarriers = (numtiles<numthreads ? numtiles : numthreads) + 1;
-      gal_threads_attr_barrier_init(&attr, &b, numbarriers);
-
-      /* Spin off the threads: */
-      for(i=0;i<numthreads;++i)
-        if(indexs[i*thrdcols]!=GAL_THREADS_NON_THRD_INDEX)
-          {
-            prm[i].id=i;
-            prm[i].b=&b;
-            prm[i].params=caller_params;
-            prm[i].indexs=&indexs[i*thrdcols];
-            err=pthread_create(&t, &attr, function, &prm[i]);
-            if(err)
-              {
-                fprintf(stderr, "can't create thread %zu", i);
-                exit(EXIT_FAILURE);
-              }
-          }
-
-      /* Wait for all threads to finish and free the spaces. */
-      pthread_barrier_wait(&b);
-      pthread_attr_destroy(&attr);
-      pthread_barrier_destroy(&b);
-    }
-
-  /* Clean up. */
-  free(indexs);
-}