[gnuastro-commits] master 80cdd4f 118/125: Correction of convolution over channel borders

[Mohammad Akhlaghi]

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

    Correction of convolution over channel borders
    
    When convolution is done in each channel separately, the boundaries of the
    channels cause a discontinuity. This is ofcourse the desired effect. But
    afterwards, when the effect of the channels has been subtracted, we want to
    remove those edges in an already convolved image (this is what happens in
    NoiseChisel). Since convolution is the most time consuming step of
    NoiseChisel (or many other tools), is a significant over-kill to convolve
    the whole image again just to remove the edges!
    
    To fix this issue, the new `gal_convolve_spatial_correct_ch_edge' function
    was defined. It will only convolve the pixels on the edges of channels. So
    all the other pixels in the image will be left untouched and it is possible
    to correct an already convolved image really fast and easy. Since the
    implementation of convolution was effectively the same as before, only some
    corrections were made to the main spatial convolution function (the new
    `gal_convolve_spatial_general'). As a result the `gal_convolve_spatial' and
    `gal_convolve_spatial_correct_ch_edge' are now just wrappers for that more
    general function.
    
    This was done as part of work on NoiseChisel, and in the process (to make
    things more tidy in NoiseChisel), the Sky and its STD calculation and
    subtractky functions were moved to a new `sky.c' source file.
---
 bin/noisechisel/Makefile.am            |   5 +-
 bin/noisechisel/args.h                 |   2 +-
 bin/noisechisel/detection.c            |  23 +-
 bin/noisechisel/noisechisel.c          |  92 ++++++-
 bin/noisechisel/sky.c                  | 273 ++++++++++++++++++++
 bin/noisechisel/{threshold.h => sky.h} |  29 +--
 bin/noisechisel/threshold.c            | 180 +------------
 bin/noisechisel/threshold.h            |  10 +-
 lib/convolve.c                         | 449 ++++++++++++++++++++++-----------
 lib/gnuastro/convolve.h                |   6 +
 10 files changed, 703 insertions(+), 366 deletions(-)

diff --git a/bin/noisechisel/Makefile.am b/bin/noisechisel/Makefile.am
index bb9793c..b5fe436 100644
--- a/bin/noisechisel/Makefile.am
+++ b/bin/noisechisel/Makefile.am
@@ -30,10 +30,11 @@ bin_PROGRAMS = astnoisechisel
 
 astnoisechisel_LDADD = -lgnuastro
 
-astnoisechisel_SOURCES = main.c ui.c detection.c noisechisel.c threshold.c
+astnoisechisel_SOURCES = main.c ui.c detection.c noisechisel.c sky.c     \
+  threshold.c
 
 EXTRA_DIST = main.h authors-cite.h args.h ui.h detection.h noisechisel.h \
-  threshold.h
+  sky.h threshold.h
 
 
 
diff --git a/bin/noisechisel/args.h b/bin/noisechisel/args.h
index ee5f503..9c713ae 100644
--- a/bin/noisechisel/args.h
+++ b/bin/noisechisel/args.h
@@ -397,7 +397,7 @@ struct argp_option program_options[] =
       0,
       "Final sky and its STD steps in a file.",
       ARGS_GROUP_DETECTION,
-      &p->checkdetection,
+      &p->checksky,
       GAL_OPTIONS_NO_ARG_TYPE,
       GAL_OPTIONS_RANGE_0_OR_1,
       GAL_OPTIONS_NOT_MANDATORY,
diff --git a/bin/noisechisel/detection.c b/bin/noisechisel/detection.c
index d67fd8b..7df73e5 100644
--- a/bin/noisechisel/detection.c
+++ b/bin/noisechisel/detection.c
@@ -37,6 +37,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include "main.h"
 
 #include "ui.h"
+#include "sky.h"
 #include "threshold.h"
 
 
@@ -586,10 +587,13 @@ detection_remove_false_initial(struct noisechiselparams 
*p,
         ++b;
       }
     while(++l<lf);
-  else
-    do
-      if(*l!=GAL_BLANK_UINT32)
-        *l = newlabels[ *l ];
+  else                               /* We need the binary array even when */
+    do                               /* there is no dilation: the binary   */
+      {                              /* array is used for estimating the   */
+        if(*l!=GAL_BLANK_UINT32)     /* Sky and its STD. */
+          *b = ( *l = newlabels[ *l ] ) > 0;
+        ++b;
+      }
     while(++l<lf);
 
 
@@ -622,7 +626,7 @@ detection(struct noisechiselparams *p)
 
   /* Find the Sky and its Standard Deviation from the initial detectios. */
   if(!p->cp.quiet) gettimeofday(&t1, NULL);
-  threshold_sky_and_std(p);
+  sky_and_std(p, p->detskyname);
   if(!p->cp.quiet)
     gal_timing_report(&t1, "Initial (crude) Sky and its STD found.", 2);
 
@@ -655,7 +659,6 @@ detection(struct noisechiselparams *p)
     }
 
 
-
   /* If the user asked for dilation, then apply it. */
   if(p->dilate)
     {
@@ -678,8 +681,12 @@ detection(struct noisechiselparams *p)
     }
 
 
-
-  /*  */
+  /* p->binary was used to keep the initial pseudo-detection threshold. But
+     we don't need it any more, so we'll just free it and put the `workbin'
+     array in its place. Note that `workbin' has a map of all the detected
+     objects, which is still necessary during NoiseChisel. */
+  gal_data_free(p->binary);
+  p->binary=workbin;
 
 
   /* If the user wanted to check the threshold and hasn't called
diff --git a/bin/noisechisel/noisechisel.c b/bin/noisechisel/noisechisel.c
index 2244460..1fb90d1 100644
--- a/bin/noisechisel/noisechisel.c
+++ b/bin/noisechisel/noisechisel.c
@@ -29,20 +29,34 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/convolve.h>
 
 #include <gnuastro-internal/timing.h>
+#include <gnuastro-internal/checkset.h>
 
 #include "main.h"
 
+#include "ui.h"
+#include "sky.h"
 #include "detection.h"
+#include "threshold.h"
 
 
 
 
+
+
+
+
+
+
+/***********************************************************************/
+/*************  Wrapper functions (for clean high-level) ***************/
+/***********************************************************************/
 static void
 noisechisel_convolve(struct noisechiselparams *p)
 {
   struct timeval t1;
   struct gal_tile_two_layer_params *tl=&p->cp.tl;
 
+  /* Do the convolution. */
   if(!p->cp.quiet) gettimeofday(&t1, NULL);
   p->conv=gal_convolve_spatial(tl->tiles, p->kernel, p->cp.numthreads,
                                1, tl->workoverch);
@@ -59,15 +73,89 @@ noisechisel_convolve(struct noisechiselparams *p)
 
 
 
+static void
+noisechisel_find_sky_subtract(struct noisechiselparams *p)
+{
+  /* Free the initial Sky and its STD estimates. */
+  gal_data_free(p->sky);
+  gal_data_free(p->std);
+
+  /* Find the final Sky value. */
+  sky_and_std(p, p->skyname);
+
+  /* Abort if the user only wanted to see until this point.*/
+  if(p->skyname && !p->continueaftercheck)
+    ui_abort_after_check(p, p->skyname, "derivation of final Sky (and its "
+                         "STD) value");
+
+  /* Subtract the Sky from the Input and Convolved (necessary for
+     segmentation) images. */
+  sky_subtract(p);
+}
+
+
+
+
+
+/* If convolution was not done while respecting channel edges (when there
+   is more than one channel, pixels outside the edge weren't used in the
+   convolution), then correct it. */
+static void
+noisechisel_convolve_correct_ch_edges(struct noisechiselparams *p)
+{
+  struct gal_tile_two_layer_params *tl=&p->cp.tl;
+
+  gal_checkset_check_remove_file("conv-correct.fits", 0, 0);
+  gal_fits_img_write(p->conv, "conv-correct.fits", NULL, PROGRAM_STRING);
+
+  if( tl->totchannels>1 && tl->workoverch==0 )
+    gal_convolve_spatial_correct_ch_edge(tl->tiles, p->kernel,
+                                         p->cp.numthreads, 1, p->conv);
+
+  gal_fits_img_write(p->conv, "conv-correct.fits", NULL, PROGRAM_STRING);
+
+  if(!p->cp.quiet)
+    gal_timing_report(NULL, "Corrected convolution of touching channel "
+                      "edges", 1);
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/***********************************************************************/
+/*************             High level function           ***************/
+/***********************************************************************/
 void
 noisechisel(struct noisechiselparams *p)
 {
   /* Convolve the image. */
   noisechisel_convolve(p);
 
-  /* Do the initial detection: */
+  /* Do the initial detection. */
   detection_initial(p);
 
-  /* Remove false detections */
+  /* Remove false detections. */
   detection(p);
+
+  /* Find the Sky value and subtract it. */
+  noisechisel_find_sky_subtract(p);
+
+  /* Correct the convolved image channel edges if necessary. */
+  noisechisel_convolve_correct_ch_edges(p);
 }
diff --git a/bin/noisechisel/sky.c b/bin/noisechisel/sky.c
new file mode 100644
index 0000000..92dac5e
--- /dev/null
+++ b/bin/noisechisel/sky.c
@@ -0,0 +1,273 @@
+/*********************************************************************
+NoiseChisel - Detect and segment signal in a noisy dataset.
+NoiseChisel is part of GNU Astronomy Utilities (Gnuastro) package.
+
+Original author:
+     Mohammad Akhlaghi <address@hidden>
+Contributing author(s):
+Copyright (C) 2015, Free Software Foundation, Inc.
+
+Gnuastro is free software: you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation, either version 3 of the License, or (at your
+option) any later version.
+
+Gnuastro is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
+**********************************************************************/
+#include <config.h>
+
+#include <stdio.h>
+#include <errno.h>
+#include <error.h>
+#include <string.h>
+#include <stdlib.h>
+
+#include <gnuastro/tile.h>
+#include <gnuastro/fits.h>
+#include <gnuastro/blank.h>
+#include <gnuastro/threads.h>
+#include <gnuastro/statistics.h>
+
+#include "main.h"
+
+#include "ui.h"
+#include "threshold.h"
+
+
+
+
+
+
+
+
+/****************************************************************
+ ************            Estimate the Sky            ************
+ ****************************************************************/
+static void *
+sky_mean_std_undetected(void *in_prm)
+{
+  struct gal_threads_params *tprm=(struct gal_threads_params *)in_prm;
+  struct noisechiselparams *p=(struct noisechiselparams *)tprm->params;
+
+  double *darr, s, s2;
+  int type=p->sky->type;
+  size_t i, tind, numsky, dsize=2;
+  gal_data_t *tile, *meanstd_d, *meanstd, *bintile;
+
+
+  /* A dataset to keep the mean and STD in double type. */
+  meanstd_d=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &dsize,
+                           NULL, 0, -1, NULL, NULL, NULL);
+  darr=meanstd_d->array;
+
+
+  /* An empty dataset to replicate a tile on the binary array. */
+  bintile=gal_data_alloc(NULL, GAL_TYPE_UINT8, 1, &dsize,
+                         NULL, 0, -1, NULL, NULL, NULL);
+  free(bintile->array);
+  free(bintile->dsize);
+  bintile->block=p->binary;
+  bintile->ndim=p->binary->ndim;
+
+
+  /* Go over all the tiles given to this thread. */
+  for(i=0; tprm->indexs[i] != GAL_THREADS_NON_THRD_INDEX; ++i)
+    {
+      /* Basic definitions */
+      numsky=0;
+      tind = tprm->indexs[i];
+      tile = &p->cp.tl.tiles[tind];
+
+      /* Correct the fake binary tile's properties to be the same as this
+         one, then count the number of zero valued elements in it. */
+      bintile->size=tile->size;
+      bintile->dsize=tile->dsize;
+      bintile->array=gal_tile_block_relative_to_other(tile, p->binary);
+      GAL_TILE_PARSE_OPERATE({if(!*o) numsky++;}, tile, bintile, 1, 1);
+
+      /* Only continue, if the fraction of Sky values are less than the
+         requested fraction. */
+      if( (float)(numsky)/(float)(tile->size) > p->minbfrac)
+        {
+          /* Calculate the mean and STD over this tile. */
+          s=s2=0.0f;
+          GAL_TILE_PARSE_OPERATE(
+                                 {
+                                   if(!*o)
+                                     {
+                                       s  += *i;
+                                       s2 += *i * *i;
+                                     }
+                                 }, tile, bintile, 1, 1);
+          darr[0]=s/numsky;
+          darr[1]=sqrt( (s2-s*s/numsky)/numsky );
+
+          /* Convert the mean and std into the same type as the sky and std
+             arrays. */
+          meanstd=gal_data_copy_to_new_type(meanstd_d, type);
+
+          /* Copy the mean and STD to their respective places in the tile
+             arrays. */
+          memcpy(gal_data_ptr_increment(p->sky->array, tind, type),
+                 meanstd->array, gal_type_sizeof(type));
+          memcpy(gal_data_ptr_increment(p->std->array, tind, type),
+                 gal_data_ptr_increment(meanstd->array, 1, type),
+                 gal_type_sizeof(type));
+        }
+      else
+        {
+          gal_blank_write(gal_data_ptr_increment(p->sky->array, tind, type),
+                          type);
+          gal_blank_write(gal_data_ptr_increment(p->std->array, tind, type),
+                          type);
+        }
+    }
+
+  /* Clean up and wait for other threads to finish and abort. */
+  bintile->array=NULL;
+  bintile->dsize=NULL;
+  gal_data_free(bintile);
+  gal_data_free(meanstd_d);
+  if(tprm->b) pthread_barrier_wait(tprm->b);
+  return NULL;
+}
+
+
+
+
+
+void
+sky_and_std(struct noisechiselparams *p, char *checkname)
+{
+  gal_data_t *tmp;
+  struct gal_options_common_params *cp=&p->cp;
+  struct gal_tile_two_layer_params *tl=&cp->tl;
+
+
+  /* When the check image has the same resolution as the input, write the
+     binary array as a reference to help in the comparison. */
+  if(checkname && !tl->oneelempertile)
+    gal_fits_img_write(p->binary, checkname, NULL, PROGRAM_STRING);
+
+
+  /* Allocate space for the mean and standard deviation. */
+  p->sky=gal_data_alloc(NULL, p->input->type, p->input->ndim, tl->numtiles,
+                        NULL, 0, cp->minmapsize, "SKY", p->input->unit, NULL);
+  p->std=gal_data_alloc(NULL, p->input->type, p->input->ndim, tl->numtiles,
+                        NULL, 0, cp->minmapsize, "STD", p->input->unit, NULL);
+
+
+  /* Find the Sky and its STD on proper tiles. */
+  gal_threads_spin_off(sky_mean_std_undetected, p, tl->tottiles,
+                       cp->numthreads);
+  if(checkname)
+    {
+      gal_tile_full_values_write(p->sky, tl, checkname, PROGRAM_STRING);
+      gal_tile_full_values_write(p->std, tl, checkname, PROGRAM_STRING);
+    }
+
+  /* Get the basic information about the standard deviation
+     distribution. */
+  tmp=gal_statistics_median(p->std, 0);
+  tmp=gal_data_copy_to_new_type(tmp, GAL_TYPE_FLOAT32);
+  memcpy(&p->medstd, tmp->array, sizeof p->medstd);
+  free(tmp);
+
+  tmp=gal_statistics_minimum(p->std);
+  tmp=gal_data_copy_to_new_type(tmp, GAL_TYPE_FLOAT32);
+  memcpy(&p->minstd, tmp->array, sizeof p->minstd);
+  free(tmp);
+
+  tmp=gal_statistics_maximum(p->std);
+  tmp=gal_data_copy_to_new_type(tmp, GAL_TYPE_FLOAT32);
+  memcpy(&p->maxstd, tmp->array, sizeof p->maxstd);
+  free(tmp);
+
+  /* In case the image is in electrons or counts per second, the standard
+     deviation of the noise will become smaller than unity, so we need to
+     correct it in the S/N calculation. So, we'll calculate the correction
+     factor here. */
+  p->cpscorr = p->minstd>1 ? 1.0f : p->minstd;
+
+  /* Interpolate and smooth the derived values. */
+  threshold_interp_smooth(p, &p->sky, &p->std, checkname);
+
+
+  /* If a check was requested, abort NoiseChisel. */
+  if(checkname && !p->continueaftercheck)
+    ui_abort_after_check(p, checkname, "showing derivation of Sky value"
+                         "and its standard deviation, or STD");
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/****************************************************************
+ ************            Subtract the Sky            ************
+ ****************************************************************/
+void
+sky_subtract(struct noisechiselparams *p)
+{
+  size_t tid;
+  void *tarray=NULL;
+  float *sky=p->sky->array;
+  gal_data_t *tile, *tblock=NULL;
+
+  /* A small sanity check. */
+  if(p->sky->type!=GAL_TYPE_FLOAT32)
+    error(EXIT_FAILURE, 0, "`sky_subtract' currently only works on float32 "
+          "type sky values.");
+
+  /* Go over all the tiles. */
+  for(tid=0; tid<p->cp.tl.tottiles; ++tid)
+    {
+      /* For easy reading. */
+      tile=&p->cp.tl.tiles[tid];
+
+      /* First subtract the Sky value from the input image. */
+      GAL_TILE_PARSE_OPERATE({*i-=sky[tid];}, tile, NULL, 0, 0);
+
+      /* Change to the convolved image. */
+      tarray=tile->array;
+      tblock=tile->block;
+      tile->array=gal_tile_block_relative_to_other(tile, p->conv);
+      tile->block=p->conv;
+
+      /* The threshold is always low. So for the majority of non-NaN
+         pixels in the image, the condition above will be true. If we
+         come over a NaN pixel, then by definition of NaN, all
+         conditionals will fail.
+
+         If an image doesn't have any NaN pixels, only the pixels below
+         the threshold have to be checked for a NaN which are by
+         definition a very small fraction of the total pixels. And if
+         there are NaN pixels in the image. */
+      GAL_TILE_PARSE_OPERATE({*i-=sky[tid];}, tile, NULL, 0, 0);
+
+      /* Revert back to the original block. */
+      tile->array=tarray;
+      tile->block=tblock;
+    }
+}
diff --git a/bin/noisechisel/threshold.h b/bin/noisechisel/sky.h
similarity index 63%
copy from bin/noisechisel/threshold.h
copy to bin/noisechisel/sky.h
index 159699f..08a7d36 100644
--- a/bin/noisechisel/threshold.h
+++ b/bin/noisechisel/sky.h
@@ -20,34 +20,13 @@ General Public License for more details.
 You should have received a copy of the GNU General Public License
 along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
 **********************************************************************/
-#ifndef THRESHOLD_H
-#define THRESHOLD_H
-
-#define THRESHOLD_NO_ERODE_VALUE 2
-
-
-enum threshold_type
-  {
-    THRESHOLD_QUANTILES,
-    THRESHOLD_SKY_STD,
-  };
-
-
-void
-threshold_quantile_find_apply(struct noisechiselparams *p);
+#ifndef SKY_H
+#define SKY_H
 
 void
-threshold_write_sn_table(struct noisechiselparams *p, gal_data_t *sntable,
-                         gal_data_t *snind, char *filename,
-                         struct gal_linkedlist_stll *comments);
+sky_and_std(struct noisechiselparams *p, char *checkname);
 
 void
-threshold_sky_and_std(struct noisechiselparams *p);
-
-
-void
-threshold_apply(struct noisechiselparams *p, float *value1, float *value2,
-                int type);
-
+sky_subtract(struct noisechiselparams *p);
 
 #endif
diff --git a/bin/noisechisel/threshold.c b/bin/noisechisel/threshold.c
index e47a794..9f9fc87 100644
--- a/bin/noisechisel/threshold.c
+++ b/bin/noisechisel/threshold.c
@@ -210,7 +210,7 @@ threshold_write_sn_table(struct noisechiselparams *p, 
gal_data_t *insn,
 /***************     Interpolation and smoothing      *****************/
 /**********************************************************************/
 /* Interpolate and smooth the values for each tile over the whole image. */
-static void
+void
 threshold_interp_smooth(struct noisechiselparams *p, gal_data_t **first,
                         gal_data_t **second, char *filename)
 {
@@ -475,181 +475,3 @@ threshold_quantile_find_apply(struct noisechiselparams *p)
   if(p->qthreshname && !p->continueaftercheck)
     ui_abort_after_check(p, p->qthreshname, "quantile threshold check");
 }
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-/****************************************************************
- ************       Mean and STD of undetected       ************
- ****************************************************************/
-static void *
-threshold_mean_std_undetected(void *in_prm)
-{
-  struct gal_threads_params *tprm=(struct gal_threads_params *)in_prm;
-  struct noisechiselparams *p=(struct noisechiselparams *)tprm->params;
-
-  double *darr, s, s2;
-  int type=p->sky->type;
-  size_t i, tind, numsky, dsize=2;
-  gal_data_t *tile, *meanstd_d, *meanstd, *bintile;
-
-
-  /* A dataset to keep the mean and STD in double type. */
-  meanstd_d=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &dsize,
-                           NULL, 0, -1, NULL, NULL, NULL);
-  darr=meanstd_d->array;
-
-
-  /* An empty dataset to replicate a tile on the binary array. */
-  bintile=gal_data_alloc(NULL, GAL_TYPE_UINT8, 1, &dsize,
-                         NULL, 0, -1, NULL, NULL, NULL);
-  free(bintile->array);
-  free(bintile->dsize);
-  bintile->block=p->binary;
-  bintile->ndim=p->binary->ndim;
-
-
-  /* Go over all the tiles given to this thread. */
-  for(i=0; tprm->indexs[i] != GAL_THREADS_NON_THRD_INDEX; ++i)
-    {
-      /* Basic definitions */
-      numsky=0;
-      tind = tprm->indexs[i];
-      tile = &p->cp.tl.tiles[tind];
-
-      /* Correct the fake binary tile's properties to be the same as this
-         one, then count the number of zero valued elements in it. */
-      bintile->size=tile->size;
-      bintile->dsize=tile->dsize;
-      bintile->array=gal_tile_block_relative_to_other(tile, p->binary);
-      GAL_TILE_PARSE_OPERATE({if(!*o) numsky++;}, tile, bintile, 1, 1);
-
-      /* Only continue, if the fraction of Sky values are less than the
-         requested fraction. */
-      if( (float)(numsky)/(float)(tile->size) > p->minbfrac)
-        {
-          /* Calculate the mean and STD over this tile. */
-          s=s2=0.0f;
-          GAL_TILE_PARSE_OPERATE(
-                                 {
-                                   if(!*o)
-                                     {
-                                       s  += *i;
-                                       s2 += *i * *i;
-                                     }
-                                 }, tile, bintile, 1, 1);
-          darr[0]=s/numsky;
-          darr[1]=sqrt( (s2-s*s/numsky)/numsky );
-
-          /* Convert the mean and std into the same type as the sky and std
-             arrays. */
-          meanstd=gal_data_copy_to_new_type(meanstd_d, type);
-
-          /* Copy the mean and STD to their respective places in the tile
-             arrays. */
-          memcpy(gal_data_ptr_increment(p->sky->array, tind, type),
-                 meanstd->array, gal_type_sizeof(type));
-          memcpy(gal_data_ptr_increment(p->std->array, tind, type),
-                 gal_data_ptr_increment(meanstd->array, 1, type),
-                 gal_type_sizeof(type));
-        }
-      else
-        {
-          gal_blank_write(gal_data_ptr_increment(p->sky->array, tind, type),
-                          type);
-          gal_blank_write(gal_data_ptr_increment(p->std->array, tind, type),
-                          type);
-        }
-    }
-
-  /* Clean up and wait for other threads to finish and abort. */
-  bintile->array=NULL;
-  bintile->dsize=NULL;
-  gal_data_free(bintile);
-  gal_data_free(meanstd_d);
-  if(tprm->b) pthread_barrier_wait(tprm->b);
-  return NULL;
-}
-
-
-
-
-
-void
-threshold_sky_and_std(struct noisechiselparams *p)
-{
-  gal_data_t *tmp;
-  struct gal_options_common_params *cp=&p->cp;
-  struct gal_tile_two_layer_params *tl=&cp->tl;
-
-
-  /* When the check image has the same resolution as the input, write the
-     binary array as a reference to help in the comparison. */
-  if(p->detskyname && !tl->oneelempertile)
-    gal_fits_img_write(p->binary, p->detskyname, NULL, PROGRAM_STRING);
-
-
-  /* Allocate space for the mean and standard deviation. */
-  p->sky=gal_data_alloc(NULL, p->input->type, p->input->ndim, tl->numtiles,
-                        NULL, 0, cp->minmapsize, "SKY", p->input->unit, NULL);
-  p->std=gal_data_alloc(NULL, p->input->type, p->input->ndim, tl->numtiles,
-                        NULL, 0, cp->minmapsize, "STD", p->input->unit, NULL);
-
-
-  /* Find the Sky and its STD on proper tiles. */
-  gal_threads_spin_off(threshold_mean_std_undetected, p, tl->tottiles,
-                       cp->numthreads);
-  if(p->detskyname)
-    {
-      gal_tile_full_values_write(p->sky, tl, p->detskyname, PROGRAM_STRING);
-      gal_tile_full_values_write(p->std, tl, p->detskyname, PROGRAM_STRING);
-    }
-
-  /* Get the basic information about the standard deviation
-     distribution. */
-  tmp=gal_statistics_median(p->std, 0);
-  tmp=gal_data_copy_to_new_type(tmp, GAL_TYPE_FLOAT32);
-  memcpy(&p->medstd, tmp->array, sizeof p->medstd);
-  free(tmp);
-
-  tmp=gal_statistics_minimum(p->std);
-  tmp=gal_data_copy_to_new_type(tmp, GAL_TYPE_FLOAT32);
-  memcpy(&p->minstd, tmp->array, sizeof p->minstd);
-  free(tmp);
-
-  tmp=gal_statistics_maximum(p->std);
-  tmp=gal_data_copy_to_new_type(tmp, GAL_TYPE_FLOAT32);
-  memcpy(&p->maxstd, tmp->array, sizeof p->maxstd);
-  free(tmp);
-
-  /* In case the image is in electrons or counts per second, the standard
-     deviation of the noise will become smaller than unity, so we need to
-     correct it in the S/N calculation. So, we'll calculate the correction
-     factor here. */
-  p->cpscorr = p->minstd>1 ? 1.0f : p->minstd;
-
-  /* Interpolate and smooth the derived values. */
-  threshold_interp_smooth(p, &p->sky, &p->std, p->detskyname);
-
-
-  /* If a check was requested, abort NoiseChisel. */
-  if(p->detskyname && !p->continueaftercheck)
-    ui_abort_after_check(p, p->detskyname, "showing derivation of Sky value"
-                         "and its standard deviation, or STD");
-}
diff --git a/bin/noisechisel/threshold.h b/bin/noisechisel/threshold.h
index 159699f..606dddd 100644
--- a/bin/noisechisel/threshold.h
+++ b/bin/noisechisel/threshold.h
@@ -34,7 +34,8 @@ enum threshold_type
 
 
 void
-threshold_quantile_find_apply(struct noisechiselparams *p);
+threshold_apply(struct noisechiselparams *p, float *value1, float *value2,
+                int type);
 
 void
 threshold_write_sn_table(struct noisechiselparams *p, gal_data_t *sntable,
@@ -42,12 +43,11 @@ threshold_write_sn_table(struct noisechiselparams *p, 
gal_data_t *sntable,
                          struct gal_linkedlist_stll *comments);
 
 void
-threshold_sky_and_std(struct noisechiselparams *p);
-
+threshold_interp_smooth(struct noisechiselparams *p, gal_data_t **first,
+                        gal_data_t **second, char *filename);
 
 void
-threshold_apply(struct noisechiselparams *p, float *value1, float *value2,
-                int type);
+threshold_quantile_find_apply(struct noisechiselparams *p);
 
 
 #endif
diff --git a/lib/convolve.c b/lib/convolve.c
index 87e00f9..ca0a685 100644
--- a/lib/convolve.c
+++ b/lib/convolve.c
@@ -86,13 +86,38 @@ convolve_tile_is_on_edge(size_t *h, size_t 
*start_end_coord, size_t *k,
 /*********************************************************************/
 /********************     Spatial convolution     ********************/
 /*********************************************************************/
+/* Each thread needs one of these structures. */
+struct per_thread_spatial_prm
+{
+  /* Internally stored/used values. */
+  gal_data_t      *tile;     /* Tile this thread is working on.          */
+  gal_data_t   *overlap;     /* Overlap pointer and starting point.      */
+  size_t *overlap_start;     /* Starting coordinate of kernel overlap.   */
+  size_t  *kernel_start;     /* Kernel starting point.                   */
+  size_t    *host_start;     /* Starting coordinate of host.             */
+  size_t           *pix;     /* 2*ndim: starting and ending of tile,
+                                Later, just the pixel being convolved.   */
+  int           on_edge;     /* If the tile is on the edge or not.       */
+  gal_data_t      *host;     /* Size of host (channel or block).         */
+  size_t    k_start_inc;     /* Increment for kernel.                    */
+  struct spatial_params *cprm; /* Link to main structure for all threads.*/
+};
+
+
+
+
+/* This is for the full list. */
 struct spatial_params
 {
-  gal_data_t    *out;
-  gal_data_t  *tiles;
-  gal_data_t *kernel;
-  int     convoverch;
-  int edgecorrection;
+  /* Main input/output parameters. */
+  gal_data_t       *out;     /* Output data structure.                   */
+  gal_data_t     *tiles;     /* Tiles over the input image.              */
+  gal_data_t     *block;     /* Pointer to block for this tile.          */
+  gal_data_t    *kernel;     /* Kernel to convolve with input.           */
+  gal_data_t *tocorrect;     /* (possible) convolved image to correct.   */
+  int        convoverch;     /* Ignore channel edges in convolution.     */
+  int    edgecorrection;     /* Correct convolution's edge effects.      */
+  struct per_thread_spatial_prm *pprm; /* Array of per-thread parameters.*/
 };
 
 
@@ -100,28 +125,32 @@ struct spatial_params
 
 /* Define the overlap of the kernel and image over this part of the image,
    the necessary input image parameters are stored in `overlap' (its
-   `array' and `dsize' elements). The pointer to the kernel array that the
-   overlap starts will be returned. */
+   `array' and `dsize' elements).  */
 static int
-convolve_spatial_overlap(int on_edge, size_t *parse_coords, size_t *hsize,
-                         gal_data_t *block, gal_data_t *kernel,
-                         gal_data_t *overlap, size_t *k_start_inc,
-                         size_t tile_ind)
+convolve_spatial_overlap(struct per_thread_spatial_prm *pprm,
+                         int tocorrect)
 {
-  size_t overlap_inc, ndim=kernel->ndim;
+  struct spatial_params *cprm=pprm->cprm;
+  gal_data_t *block=cprm->block, *kernel=cprm->kernel;
+  size_t *dsize = tocorrect ? block->dsize : pprm->host->dsize;
+
+  size_t *pp, *ppf, *hs;
+  size_t overlap_inc, ndim=block->ndim;
+  size_t *h=dsize, *os=pprm->overlap_start;
+  size_t *k=kernel->dsize, *ks=pprm->kernel_start;
   int full_overlap=1, dim_full_overlap, is_start, is_end;
+  size_t *p=pprm->pix, *pf=pprm->pix+ndim, *od=pprm->overlap->dsize;
+
+
+  /* In to-correct mode, the pix position needs to be relative to the
+     block. */
+  if(tocorrect)
+    {
+      hs=pprm->host_start;
+      ppf=(pp=pprm->pix)+ndim; do *pp += *hs++; while(++pp<ppf);
+    }
 
-  size_t *pix           = parse_coords; /* needs 2*ndim, also for end. */
-  size_t *overlap_start = parse_coords + ndim * 3;
-  size_t *kernel_start  = parse_coords + ndim * 4;
-  size_t *host_start    = parse_coords + ndim * 5;
 
-  size_t *p=pix, *pf=pix+ndim, *os=overlap_start, *h=hsize;
-  size_t *k=kernel->dsize, *ks=kernel_start, *od=overlap->dsize;
-  /*
-  if(tile_ind==2053)
-    printf("pix: %zu, %zu\n", pix[0], pix[1]);
-  */
   /* Coordinate to start convolution for this pixel. */
   do
     {
@@ -129,16 +158,17 @@ convolve_spatial_overlap(int on_edge, size_t 
*parse_coords, size_t *hsize,
          overlaps because this is the most common case usually). */
       dim_full_overlap=1;
 
-
       /* When the tile is on the edge, still some pixels in it can have
          full overlap. So using the `dim_full_overlap', we will do the same
-         thing we do for the tiles that don't overlap for them. */
-      if(on_edge)
+         thing we do for the tiles that don't overlap for them. When
+         `tocorrect!=0', then only pixels that are on the edge of the tile
+         will get to this point, so it must always be checked. */
+      if( tocorrect ? 1 : pprm->on_edge )
         {
           /* See if this pixel is on the start and/or end of the dimension
              relative to the kernel. */
-          is_start = *p < *k/2;
-          is_end   = *p + *k/2 >= *h;
+          is_start = *p          <  *k/2;
+          is_end   = (*p + *k/2) >= *h;
           if( is_start || is_end )
             {
               /* Overlapping with the outside.
@@ -199,36 +229,57 @@ convolve_spatial_overlap(int on_edge, size_t 
*parse_coords, size_t *hsize,
         }
     }
   while(++p<pf);
-  /*
-  if(tile_ind==2053)
+
+
+  /* To check, add an `int check' argument to the function.
+  if(check)
     {
+      printf("pix (within %s): %zu, %zu\n", tocorrect ? "full" : "channel",
+             pprm->pix[0], pprm->pix[1]);
       printf("\tk/2: %zu, %zu\n", kernel->dsize[0]/2, kernel->dsize[1]/2);
-      printf("\toverlap_start: %zu, %zu\n", overlap_start[0],
-             overlap_start[1]);
-      printf("\toverlap->dsize: %zu, %zu\n", overlap->dsize[0],
-             overlap->dsize[1]);
-      exit(0);
+      printf("\th: %zu, %zu\n", dsize[0], dsize[1]);
+      printf("\toverlap_start: %zu, %zu\n", pprm->overlap_start[0],
+             pprm->overlap_start[1]);
+      printf("\toverlap->dsize: %zu, %zu\n", pprm->overlap->dsize[0],
+             pprm->overlap->dsize[1]);
+      printf("\tfulloverlap: %d\n", full_overlap);
     }
   */
+
   /* Set the increment to start working on the kernel. */
-  *k_start_inc = ( full_overlap
-                   ? 0
-                   : gal_dimension_coord_to_index(ndim, kernel->dsize,
-                                                  kernel_start) );
-
-  /* Add the host's starting location (necessary when convolution over the
-     host/channel is treated independently). Until now we worked as if the
-     the host/channel is the full image so the edges don't get mixed. But
-     from now on we will be working over the allocated block to look at
-     pixel values, so we need to convert the location to the proper place
-     within the allocated array. */
-  gal_dimension_add_coords(overlap_start, host_start, overlap_start, ndim);
+  pprm->k_start_inc = ( full_overlap
+                        ? 0
+                        : gal_dimension_coord_to_index(ndim,
+                                                       kernel->dsize,
+                                                       pprm->kernel_start) );
+
+  /* Make correction.
+
+      Normal mode (when `tocorrect==0'): add the host's starting location
+         (necessary when convolution over the host/channel is treated
+         independently). In this mode, until now we were working as if the
+         the host/channel is the full image so the edges don't get
+         mixed. But from now on we will be working over the allocated block
+         to look at pixel values, so we need to convert the location to the
+         proper place within the allocated array.
+
+      To-correct mode: The boundaries were calculated with respect to the
+         block, so we don't need to correct `overlap_start'. But we need to
+         correct the pixel position back to its original state (relative to
+         the channel). */
+  hs=pprm->host_start;
+  if(tocorrect)
+    { ppf=(pp=pprm->pix)+ndim; do *pp -= *hs++; while(++pp<ppf); }
+  else
+    { ppf=(pp=pprm->overlap_start)+ndim; do *pp += *hs++; while(++pp<ppf); }
+
 
   /* Set the increment to start working on the overlap region and use that
      to set the starting pointer of the overlap region. */
-  overlap_inc=gal_dimension_coord_to_index(ndim, block->dsize, overlap_start);
-  overlap->array=gal_data_ptr_increment(block->array, overlap_inc,
-                                        block->type);
+  overlap_inc=gal_dimension_coord_to_index(ndim, block->dsize,
+                                           pprm->overlap_start);
+  pprm->overlap->array=gal_data_ptr_increment(block->array, overlap_inc,
+                                              block->type);
   return full_overlap;
 }
 
@@ -238,34 +289,30 @@ convolve_spatial_overlap(int on_edge, size_t 
*parse_coords, size_t *hsize,
 
 /* Convolve over one tile that is not touching the edge. */
 static void
-convolve_spatial_tile(struct spatial_params *cprm, size_t tile_ind,
-                      size_t *parse_coords, gal_data_t *overlap)
+convolve_spatial_tile(struct per_thread_spatial_prm *pprm)
 {
-  gal_data_t *tile=&cprm->tiles[tile_ind];
-  gal_data_t *block=gal_tile_block(tile), *kernel=cprm->kernel;
 
   double sum, ksum;
-  int on_edge, full_overlap;
-  size_t i, ndim=tile->ndim, csize=tile->dsize[ndim-1];
-  gal_data_t *host=cprm->convoverch ? block : tile->block;
+  int full_overlap;
+  struct spatial_params *cprm=pprm->cprm;
+  gal_data_t *tile=pprm->tile, *overlap=pprm->overlap;
+  gal_data_t *block=cprm->block, *kernel=cprm->kernel;
+  size_t i, ndim=block->ndim, csize=tile->dsize[ndim-1];
 
   /* Variables for scanning a tile (`i_*') and the region around every
      pixel of a tile (`o_*'). */
-  size_t start_fastdim, k_start_inc;
+  size_t start_fastdim;
   size_t k_inc, i_inc, i_ninc, i_st_en[2], o_inc, o_ninc, o_st_en[2];
 
   /* These variables depend on the type of the input. */
   float *kv, *iv, *ivf, *i_start, *o_start, *k_start;
   float *in_v, *in=block->array, *out=cprm->out->array;
 
-  /* The `parse_coords' array was allocated once before this function for
-     all the tiles that are given to a thread. It has the space for all the
-     necessary coordinates. At first, `pix' will be the starting element of
-     the tile, but then it will be incremented as we carpet the tile. So we
-     aren't calling it `start'. */
-  size_t *pix           = parse_coords; /* needs 2*ndim (also for end). */
-  size_t *o_c           = parse_coords + ndim * 2;
-  size_t *host_start    = parse_coords + ndim * 5;
+
+  /* Starting pixel for this tile. Note that when we are in `tocorrect'
+     mode, this position has to be  */
+  pprm->host=cprm->convoverch ? block : tile->block;
+  gal_tile_start_coord(pprm->host, pprm->host_start);
 
 
   /* Set the starting and ending coordinates of this tile (recall that the
@@ -273,20 +320,24 @@ convolve_spatial_tile(struct spatial_params *cprm, size_t 
tile_ind,
      parse_coords). When `convoverch' is set, we want to convolve over the
      whole allocated block, not just one channel. So in effect, it is the
      same as `rel_block' in `gal_tile_start_end_coord'. */
-  gal_tile_start_end_coord(tile, parse_coords, cprm->convoverch);
-  start_fastdim = pix[ndim-1];
+  gal_tile_start_end_coord(tile, pprm->pix, cprm->convoverch);
+  start_fastdim = pprm->pix[ndim-1];
 
-  /* Set the starting coordinate of the host for this tile. */
-  gal_tile_start_coord(host, host_start);
 
   /* See if this tile is on the edge or not. */
-  on_edge=convolve_tile_is_on_edge(host->dsize, parse_coords, kernel->dsize,
-                                   ndim);
+  pprm->on_edge=convolve_tile_is_on_edge(pprm->host->dsize, pprm->pix,
+                                         kernel->dsize, ndim);
+
+
+  /* If it isn't on the edge and we are correcting an already convolved
+     image (`tocorrect!=NULL'), then this tile can be ignored. */
+  if(cprm->tocorrect && pprm->on_edge==0) return;
+
   /*
   if(tile_ind==2053)
     {
       printf("\ntile %zu...\n", tile_ind);
-      printf("\tpix: %zu, %zu\n", pix[0], pix[1]);
+      printf("\tpix: %zu, %zu\n", pprm->pix[0], pprm->pix[1]);
       exit(0);
     }
   */
@@ -297,7 +348,7 @@ convolve_spatial_tile(struct spatial_params *cprm, size_t 
tile_ind,
     {
       /* Initialize the value along the fastest dimension (it is not
          incremented during `gal_tile_block_increment'). */
-      pix[ndim-1]=start_fastdim;
+      pprm->pix[ndim-1]=start_fastdim;
 
       /* Go over each pixel to convolve. */
       for(i=0;i<csize;++i)
@@ -313,66 +364,78 @@ convolve_spatial_tile(struct spatial_params *cprm, size_t 
tile_ind,
           else
             {
               /* Define the overlap region. */
-              full_overlap=convolve_spatial_overlap(on_edge, parse_coords,
-                                                    host->dsize, block,
-                                                    kernel, overlap,
-                                                    &k_start_inc, tile_ind);
-
-              /* Set the starting pixel over the image (`o_start'). */
-              o_start=gal_tile_start_end_ind_inclusive(overlap, block,
-                                                       o_st_en);
-
-              /* Set the starting kernel pixel. Note that `kernel_array' is
-                 `void *' (pointer arithmetic is not defined on it). So we
-                 will first put it in `k_start, and then increment that. */
-              k_start=kernel->array; k_start+=k_start_inc;
-
-              /* Go over the kernel-overlap region. */
-              ksum = cprm->edgecorrection ? 0.0f : 1.0f;
-              sum=0.0f; k_inc=0; o_inc=0; o_ninc=1; kv=k_start;
-              while( o_st_en[0] + o_inc <= o_st_en[1] )
+              full_overlap=convolve_spatial_overlap(pprm, 0);
+
+              /* If tocorrect has been given and we have full overlap, then
+                 just ignore this pixel. */
+              if( !(cprm->tocorrect && full_overlap) )
                 {
-                  /* Go over the contiguous region. When there is full
-                     overlap, we don't need to calculate incrementation
-                     over the kernel, it is always a single
-                     incrementation. But when we have partial overlap,
-                     we'll need to calculate a different incrementation. */
-                  ivf = ( iv = o_start + o_inc ) + overlap->dsize[ndim-1];
-                  if(full_overlap==0) kv = k_start + k_inc;
-                  do
+
+                  /* If we are in correct mode, then re-calculate the
+                     full-overlap and all the other necessary paramters as
+                     if the channels didn't exist. */
+                  if(cprm->tocorrect)
+                    full_overlap=convolve_spatial_overlap(pprm, 1);
+
+                  /* Set the starting pixel over the image (`o_start'). */
+                  o_start=gal_tile_start_end_ind_inclusive(overlap, block,
+                                                           o_st_en);
+
+                  /* Set the starting kernel pixel. Note that
+                     `kernel_array' is `void *' (pointer arithmetic is not
+                     defined on it). So we will first put it in `k_start,
+                     and then increment that. */
+                  k_start=kernel->array; k_start += pprm->k_start_inc;
+
+                  /* Go over the kernel-overlap region. */
+                  ksum = cprm->edgecorrection ? 0.0f : 1.0f;
+                  sum=0.0f; k_inc=0; o_inc=0; o_ninc=1; kv=k_start;
+                  while( o_st_en[0] + o_inc <= o_st_en[1] )
                     {
-                      if( !isnan(*iv) )
+                      /* Go over the contiguous region. When there is full
+                         overlap, we don't need to calculate incrementation
+                         over the kernel, it is always a single
+                         incrementation. But when we have partial overlap,
+                         we'll need to calculate a different
+                         incrementation. */
+                      ivf = ( iv = o_start + o_inc ) + overlap->dsize[ndim-1];
+                      if(full_overlap==0) kv = k_start + k_inc;
+                      do
                         {
-                          sum += *iv * *kv;
-                          if(cprm->edgecorrection) ksum += *kv;
+                          if( !isnan(*iv) )
+                            {
+                              sum += *iv * *kv;
+                              if(cprm->edgecorrection) ksum += *kv;
+                            }
+                          ++kv;
                         }
-                      ++kv;
+                      while(++iv<ivf);
+
+                      /* Update the incrementation to the next contiguous
+                         region of memory over this tile. */
+                      o_inc += gal_tile_block_increment(block, overlap->dsize,
+                                                        o_ninc++, NULL);
+                      if(full_overlap==0)
+                        k_inc += gal_tile_block_increment(kernel,
+                                                          overlap->dsize,
+                                                          o_ninc-1, NULL);
                     }
-                  while(++iv<ivf);
-
-                  /* Update the incrementation to the next contiguous
-                     region of memory over this tile. Note that the
-                     contents of `o_c' are irrelevant here.*/
-                  o_inc += gal_tile_block_increment(block, overlap->dsize,
-                                                    o_ninc++, o_c);
-                  if(full_overlap==0)
-                    k_inc += gal_tile_block_increment(kernel, overlap->dsize,
-                                                      o_ninc-1, NULL);
-                }
 
-              /* Set the output value. */
-              out[ in_v - in ] = ( ksum==0.0f
-                                   ? NAN
-                                   : sum/ksum );
+                  /* Set the output value. */
+                  out[ in_v - in ] = ( ksum==0.0f
+                                       ? NAN
+                                       : sum/ksum );
+                }
             }
 
           /* Increment the last coordinate. */
-          pix[ndim-1]++;
+          pprm->pix[ndim-1]++;
         }
 
       /* Increase the increment from the start of the tile for the next
          contiguous patch. */
-      i_inc += gal_tile_block_increment(block, tile->dsize, i_ninc++, pix);
+      i_inc += gal_tile_block_increment(block, tile->dsize, i_ninc++,
+                                        pprm->pix);
     }
   /*
   if(tile_ind==2053)
@@ -390,26 +453,54 @@ convolve_spatial_on_thread(void *inparam)
 {
   struct gal_threads_params *tprm=(struct gal_threads_params *)inparam;
   struct spatial_params *cprm=(struct spatial_params *)(tprm->params);
+  gal_data_t *block=cprm->block;
 
   size_t i;
-  gal_data_t overlap, *block=gal_tile_block(cprm->tiles);
-  size_t *parse_coords=gal_data_malloc_array(GAL_TYPE_SIZE_T, 6*block->ndim);
+  struct per_thread_spatial_prm *pprm=&cprm->pprm[tprm->id];
+  size_t ndim=block->ndim, *dsize=gal_data_malloc_array(GAL_TYPE_SIZE_T,ndim);
 
-  /* Initialize the necessary overlap parameters. */
-  overlap.block=block;
-  overlap.ndim=block->ndim;
-  overlap.dsize=gal_data_malloc_array(GAL_TYPE_SIZE_T, block->ndim);
+
+  /* Set all dsize values to 1 (the values within `overlap->dsize' will be
+     changed during convolution). */
+  for(i=0;i<ndim;++i) dsize[i]=1;
+
+
+  /* Initialize/Allocate necessary items for this thread. */
+  pprm->cprm           = cprm;
+  pprm->pix            = gal_data_malloc_array(GAL_TYPE_SIZE_T, 2*ndim);
+  pprm->host_start     = gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim);
+  pprm->kernel_start   = gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim);
+  pprm->overlap_start  = gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim);
+  pprm->overlap        = gal_data_alloc(NULL, block->type, ndim, dsize,
+                                        NULL, 0, -1, NULL, NULL, NULL);
+  free(dsize);
+  free(pprm->overlap->array);
+  pprm->overlap->block = cprm->block;
 
 
   /* Go over all the tiles given to this thread. */
   for(i=0; tprm->indexs[i] != GAL_THREADS_NON_THRD_INDEX; ++i)
-    convolve_spatial_tile(cprm, tprm->indexs[i], parse_coords, &overlap);
+    {
+      /* Set this tile's pointer into this thread's parameters. */
+      pprm->tile = &cprm->tiles[ tprm->indexs[i] ];
+
+      /* Do the convolution on this tile. */
+      convolve_spatial_tile(pprm);
+    }
+
+
+  /* Set the overlap dataset's array to NULL, it was used to point to
+     different parts of the image during convolution. */
+  pprm->overlap->array=NULL;
 
 
   /* Clean up, wait until all other threads finish, then return. In a
      single thread situation, `tprm->b==NULL'. */
-  free(parse_coords);
-  free(overlap.dsize);
+  free(pprm->pix);
+  free(pprm->host_start);
+  free(pprm->kernel_start);
+  free(pprm->overlap_start);
+  gal_data_free(pprm->overlap);
   if(tprm->b) pthread_barrier_wait(tprm->b);
   return NULL;
 }
@@ -418,19 +509,18 @@ convolve_spatial_on_thread(void *inparam)
 
 
 
-/* Convolve a dataset with a given kernel in the spatial domain. Spatial
-   convolution can be greatly sped up if it is done on separate tiles over
-   the image (on multiple threads). So as input, you can either give tile
-   values or one full array. Just note that if you give a single array as
-   input, the `next' element has to be `NULL'.*/
-gal_data_t *
-gal_convolve_spatial(gal_data_t *tiles, gal_data_t *kernel,
-                     size_t numthreads, int edgecorrection, int convoverch)
+/* General spatial convolve function. This function is called by both
+   `gal_convolve_spatial' and */
+static gal_data_t *
+gal_convolve_spatial_general(gal_data_t *tiles, gal_data_t *kernel,
+                             size_t numthreads, int edgecorrection,
+                             int convoverch, gal_data_t *tocorrect)
 {
   char *name;
   struct spatial_params params;
   gal_data_t *out, *block=gal_tile_block(tiles);
 
+
   /* Small sanity checks. */
   if(tiles->ndim!=kernel->ndim)
     error(EXIT_FAILURE, 0, "The number of dimensions between the kernel and "
@@ -440,25 +530,96 @@ gal_convolve_spatial(gal_data_t *tiles, gal_data_t 
*kernel,
     error(EXIT_FAILURE, 0, "`gal_convolve_spatial' currently only works on "
           "`float32' type input and kernel");
 
-  /* Allocate the output convolved dataset. */
-  name = ( block->name
-           ? gal_checkset_malloc_cat("CONVL_", block->name) : NULL );
-  out=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, block->ndim, block->dsize,
-                     block->wcs, 0, block->minmapsize, name,
-                     block->unit, NULL);
+
+  /* Set the output datastructure.  */
+  if(tocorrect) out=tocorrect;
+  else
+    {
+      name = ( block->name
+               ? gal_checkset_malloc_cat("CONVL_", block->name) : NULL );
+      out=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, block->ndim, block->dsize,
+                         block->wcs, 0, block->minmapsize, name,
+                         block->unit, NULL);
+      if(name) free(name);
+    }
+
 
   /* Set the pointers in the parameters structure. */
   params.out=out;
   params.tiles=tiles;
+  params.block=block;
   params.kernel=kernel;
+  params.tocorrect=tocorrect;
   params.convoverch=convoverch;
   params.edgecorrection=edgecorrection;
 
+
+  /* Allocate the per-thread parameters. */
+  errno=0;
+  params.pprm=malloc(numthreads * sizeof *params.pprm);
+  if(params.pprm==NULL)
+    error(EXIT_FAILURE, 0, "%zu bytes for `params.p' in "
+          "`gal_convolve_spatial_general'", numthreads * sizeof *params.pprm);
+
+
   /* Do the spatial convolution on threads. */
   gal_threads_spin_off(convolve_spatial_on_thread, &params,
                        gal_data_num_in_ll(tiles), numthreads);
 
+
   /* Clean up and return the output array. */
-  if(name) free(name);
+  free(params.pprm);
   return out;
 }
+
+
+
+
+
+/* Convolve a dataset with a given kernel in the spatial domain. Spatial
+   convolution can be greatly sped up if it is done on separate tiles over
+   the image (on multiple threads). So as input, you can either give tile
+   values or one full array. Just note that if you give a single array as
+   input, the `next' element has to be `NULL'.*/
+gal_data_t *
+gal_convolve_spatial(gal_data_t *tiles, gal_data_t *kernel,
+                     size_t numthreads, int edgecorrection, int convoverch)
+{
+  /* Call the general function. */
+  return gal_convolve_spatial_general(tiles, kernel, numthreads,
+                                      edgecorrection, convoverch, NULL);
+}
+
+
+
+
+
+/* Correct the edges of channels in an already convolved image when it was
+   initially convolved with `gal_convolve_spatial' with `convoverch==0'. In
+   that case, strong boundaries exist on the tile edges. So if you later
+   need to remove those boundaries, you can call this function, it will
+   only do convolution on the tiles that are near the edge, not the full
+   area, so it is much faster. */
+void
+gal_convolve_spatial_correct_ch_edge(gal_data_t *tiles, gal_data_t *kernel,
+                                     size_t numthreads, int edgecorrection,
+                                     gal_data_t *tocorrect)
+{
+  gal_data_t *block=gal_tile_block(tiles);
+
+  /* Some small sanity checks. */
+  if( gal_data_dsize_is_different(block, tocorrect) )
+    error(EXIT_FAILURE, 0, "the `tocorrect' dataset has to have the same "
+          "dimensions/size as the block of the `tiles' input in "
+          "`gal_convolve_spatial_correct_ch_edge'");
+  if( block->type != tocorrect->type )
+    error(EXIT_FAILURE, 0, "the `tocorrect' dataset has to have the same "
+          "type as the block of the `tiles' input in "
+          "`gal_convolve_spatial_correct_ch_edge'. The given types are `%s' "
+          "and `%s' respectively", gal_type_to_string(tocorrect->type, 1),
+          gal_type_to_string(block->type, 1));
+
+  /* Call the general function, which will do the correction. */
+  gal_convolve_spatial_general(tiles, kernel, numthreads,
+                               edgecorrection, 0, tocorrect);
+}
diff --git a/lib/gnuastro/convolve.h b/lib/gnuastro/convolve.h
index a8bbfd0..563baed 100644
--- a/lib/gnuastro/convolve.h
+++ b/lib/gnuastro/convolve.h
@@ -49,6 +49,12 @@ gal_convolve_spatial(gal_data_t *tiles, gal_data_t *kernel,
                      size_t numthreads, int edgecorrection, int convoverch);
 
 
+void
+gal_convolve_spatial_correct_ch_edge(gal_data_t *tiles, gal_data_t *kernel,
+                                     size_t numthreads, int edgecorrection,
+                                     gal_data_t *tocorrect);
+
+
 
 __END_C_DECLS    /* From C++ preparations */