[gnuastro-commits] master adeb1c6 016/113: NoiseChisel's convolution step in 3D completed

[Mohammad Akhlaghi]

branch: master
commit adeb1c677c1a2441e8bfb7d433665b751c97f6e0
Author: Mohammad Akhlaghi <akhlaghi@gnu.org>
Commit: Mohammad Akhlaghi <akhlaghi@gnu.org>

    NoiseChisel's convolution step in 3D completed
    
    Work on adapting NoiseChisel to 3D is ongoing and with this commit, a major
    block has been removed: 3D spatial convolution is now tested and and
    working on 3D data. So NoiseChisel is now ready to start doing its
    thresholding on the images.
    
    To do the job, a default 3D kernel has been defined with this commit as a
    separate header that is included into the `ui_prepare_kernel' function (it
    was too large to be put manually into the full source and would make
    parsing by eye hard).
    
    In the process, the following issues have also been synchronized with the
    master branch:
    
     - A bug found by Vladimir in the Warp program is now implemented in this
       branch also.
    
     - The corrections and bug fixes in the tessellation and convolution
       libraries of this commit were also implemented in the master branch.
---
 NEWS                                   |   2 +
 bin/noisechisel/Makefile.am            |   2 +-
 bin/noisechisel/astnoisechisel-3d.conf |  23 ++++
 bin/noisechisel/kernel-3d.h            |  98 +++++++++++++++
 bin/noisechisel/noisechisel.c          |   6 +
 bin/noisechisel/ui.c                   | 115 +++++++++++------
 bin/warp/ui.c                          |   2 +-
 lib/convolve.c                         | 223 ++++++++++++++-------------------
 lib/gnuastro/convolve.h                |   2 +-
 lib/gnuastro/tile.h                    |   6 +-
 lib/tile.c                             |   6 +-
 11 files changed, 311 insertions(+), 174 deletions(-)

diff --git a/NEWS b/NEWS
index c59f194..964206a 100644
--- a/NEWS
+++ b/NEWS
@@ -99,6 +99,8 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
 
   Crashes on 32-bit and big-endian systems (bug #51476).
 
+  Warp's align matrix when second dimention must be reversed (bug #51536).
+
   Reading BZERO for unsigned 64-bit integers (bug #51555).
 
   Arithmetic with one file and no operators (bug #51559).
diff --git a/bin/noisechisel/Makefile.am b/bin/noisechisel/Makefile.am
index 64023c4..a0dd0cb 100644
--- a/bin/noisechisel/Makefile.am
+++ b/bin/noisechisel/Makefile.am
@@ -34,7 +34,7 @@ astnoisechisel_SOURCES = main.c ui.c clumps.c detection.c 
noisechisel.c \
   sky.c segmentation.c threshold.c
 
 EXTRA_DIST = main.h authors-cite.h args.h ui.h clumps.h detection.h     \
-  noisechisel.h segmentation.h sky.h threshold.h
+  kernel-3d.h noisechisel.h segmentation.h sky.h threshold.h
 
 
 
diff --git a/bin/noisechisel/astnoisechisel-3d.conf 
b/bin/noisechisel/astnoisechisel-3d.conf
new file mode 100644
index 0000000..f7cb3a2
--- /dev/null
+++ b/bin/noisechisel/astnoisechisel-3d.conf
@@ -0,0 +1,23 @@
+# Default parameters (System) for NoiseChisel.
+# NoiseChisel is part of GNU Astronomy Utitlies.
+#
+# Use the long option name of each parameter followed by a value. The name
+# and value should be separated by atleast one white-space character (for
+# example ` '[space], or tab). Lines starting with `#' are ignored.
+#
+# For more information, please run these commands:
+#
+#  $ astnoisechisel --help               # Full list of options, short doc.
+#  $ astnoisechisel -P                   # Print all options and used values.
+#  $ info astnoisechisel                 # All options and input/output.
+#  $ info gnuastro "Configuration files" # How to use configuration files.
+#
+# Copying and distribution of this file, with or without modification, are
+# permitted in any medium without royalty provided the copyright notice and
+# this notice are preserved.  This file is offered as-is, without any
+# warranty.
+
+# Tessellation
+ numchannels    1,1,1
+ tilesize       20,20,20
+ largetilesize  50,50,50
diff --git a/bin/noisechisel/kernel-3d.h b/bin/noisechisel/kernel-3d.h
new file mode 100644
index 0000000..6a3dbed
--- /dev/null
+++ b/bin/noisechisel/kernel-3d.h
@@ -0,0 +1,98 @@
+size_t kernel_3d_dsize[3]={7, 13, 13};
+float kernel_3d[1183]={0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 1.485799e-07, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 2.746858e-07, 1.807235e-06, 3.431401e-06, 1.838885e-06, 
2.798509e-07, 0, 0, 0, 0,
+0, 0, 0, 0, 1.859928e-06, 1.146277e-05, 2.261949e-05, 1.195174e-05, 
1.743884e-06, 0, 0, 0, 0,
+0, 0, 0, 1.457355e-07, 3.317235e-06, 2.253982e-05, 4.076775e-05, 2.208206e-05, 
3.396871e-06, 1.474514e-07, 0, 0, 0,
+0, 0, 0, 0, 1.8064e-06, 1.195441e-05, 2.113626e-05, 1.198961e-05, 1.72606e-06, 
0, 0, 0, 0,
+0, 0, 0, 0, 2.834052e-07, 1.820066e-06, 3.449281e-06, 1.89521e-06, 
2.803539e-07, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 1.44856e-07, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 1.234796e-08, 9.878367e-08, 1.020775e-06, 9.878367e-08, 
1.234796e-08, 0, 0, 0, 0,
+0, 0, 0, 4.939184e-08, 7.326431e-06, 4.693479e-05, 8.266516e-05, 4.69633e-05, 
6.999201e-06, 4.939184e-08, 0, 0, 0,
+0, 0, 1.234796e-08, 6.85008e-06, 0.0001643393, 0.001056819, 0.001962634, 
0.001032768, 0.0001649413, 6.859012e-06, 1.234796e-08, 0, 0,
+0, 0, 9.878367e-08, 4.576423e-05, 0.001047376, 0.00684491, 0.01280988, 
0.006805513, 0.001082551, 4.551054e-05, 9.878367e-08, 0, 0,
+0, 0, 9.986172e-07, 8.921415e-05, 0.001990294, 0.01254342, 0.02374495, 
0.0126314, 0.001925766, 8.45066e-05, 1.011709e-06, 0, 0,
+0, 0, 9.878367e-08, 4.649093e-05, 0.001060791, 0.006822467, 0.01247382, 
0.007007569, 0.001065096, 4.461473e-05, 9.878367e-08, 0, 0,
+0, 0, 1.234796e-08, 7.275932e-06, 0.0001640168, 0.00106862, 0.001938186, 
0.001079752, 0.0001645134, 7.191654e-06, 1.234796e-08, 0, 0,
+0, 0, 0, 4.939184e-08, 6.87645e-06, 4.4285e-05, 8.624209e-05, 4.65561e-05, 
6.938224e-06, 4.939184e-08, 0, 0, 0,
+0, 0, 0, 0, 1.234796e-08, 9.878367e-08, 1.048286e-06, 9.878367e-08, 
1.234796e-08, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 6.173979e-09, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 6.173979e-09, 1.975673e-07, 3.458353e-06, 6.472969e-06, 3.47083e-06, 
1.975673e-07, 6.173979e-09, 0, 0, 0,
+0, 0, 6.173979e-09, 7.902694e-07, 4.513072e-05, 0.0002927972, 0.0005428808, 
0.0002944311, 4.433255e-05, 7.902694e-07, 6.173979e-09, 0, 0,
+0, 0, 1.975673e-07, 4.542675e-05, 0.001045313, 0.006685175, 0.01239065, 
0.00679263, 0.001044551, 4.52873e-05, 1.975673e-07, 0, 0,
+0, 0, 3.391631e-06, 0.0002917077, 0.006763012, 0.0440967, 0.08055502, 
0.04359836, 0.006787267, 0.0002969475, 3.295347e-06, 0, 0,
+0, 6.173979e-09, 6.250698e-06, 0.0005457952, 0.01255168, 0.08104536, 0.150408, 
0.08023416, 0.01253388, 0.000535345, 6.40636e-06, 6.173979e-09, 0,
+0, 0, 1.580539e-06, 0.0002882469, 0.006713877, 0.04338823, 0.08082015, 
0.04357146, 0.006704987, 0.0002913436, 1.580539e-06, 0, 0,
+0, 0, 1.975673e-07, 4.445576e-05, 0.00105628, 0.006745212, 0.0125202, 
0.006786238, 0.001058999, 4.54068e-05, 1.975673e-07, 0, 0,
+0, 0, 6.173979e-09, 7.902694e-07, 4.506374e-05, 0.00029103, 0.0005455695, 
0.0002928559, 4.558354e-05, 7.902694e-07, 6.173979e-09, 0, 0,
+0, 0, 0, 6.173979e-09, 1.975673e-07, 3.41195e-06, 6.359724e-06, 1.580539e-06, 
1.975673e-07, 6.173979e-09, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 6.173979e-09, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 1.234796e-08, 9.878367e-08, 9.757257e-07, 9.878367e-08, 
1.234796e-08, 0, 0, 0, 0,
+0, 0, 0, 4.939184e-08, 6.978876e-06, 4.588304e-05, 8.4791e-05, 4.664944e-05, 
7.174984e-06, 4.939184e-08, 0, 0, 0,
+0, 0, 1.234796e-08, 7.018961e-06, 0.0001637559, 0.001060732, 0.001981034, 
0.0010486, 0.0001597786, 7.174912e-06, 1.234796e-08, 0, 0,
+0, 0, 9.878367e-08, 4.663123e-05, 0.001044997, 0.006716454, 0.01270765, 
0.006864154, 0.001047177, 4.60964e-05, 9.878367e-08, 0, 0,
+0, 0, 1.012701e-06, 8.26145e-05, 0.002013659, 0.0127906, 0.023479, 0.01250077, 
0.001970029, 8.642203e-05, 9.922311e-07, 0, 0,
+0, 0, 9.878367e-08, 4.552128e-05, 0.001030169, 0.006799623, 0.01275112, 
0.007000111, 0.001071859, 4.659725e-05, 9.878367e-08, 0, 0,
+0, 0, 1.234796e-08, 7.110504e-06, 0.0001573054, 0.001044423, 0.001932773, 
0.001080642, 0.0001619574, 6.927679e-06, 1.234796e-08, 0, 0,
+0, 0, 0, 4.939184e-08, 6.894978e-06, 4.457595e-05, 8.273318e-05, 4.60869e-05, 
7.289148e-06, 4.939184e-08, 0, 0, 0,
+0, 0, 0, 0, 1.234796e-08, 9.878367e-08, 9.978808e-07, 9.878367e-08, 
1.234796e-08, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 1.520727e-07, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 2.777842e-07, 1.819656e-06, 3.497153e-06, 1.862236e-06, 
2.800668e-07, 0, 0, 0, 0,
+0, 0, 0, 0, 1.814004e-06, 1.22698e-05, 2.187492e-05, 1.239061e-05, 
1.843132e-06, 0, 0, 0, 0,
+0, 0, 0, 1.385601e-07, 3.35185e-06, 2.165249e-05, 4.046146e-05, 2.121244e-05, 
3.392139e-06, 1.534694e-07, 0, 0, 0,
+0, 0, 0, 0, 1.827893e-06, 1.199963e-05, 2.301149e-05, 1.166906e-05, 
1.823535e-06, 0, 0, 0, 0,
+0, 0, 0, 0, 2.809536e-07, 1.730059e-06, 3.307363e-06, 1.8598e-06, 
2.946499e-07, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 1.439844e-07, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
diff --git a/bin/noisechisel/noisechisel.c b/bin/noisechisel/noisechisel.c
index 4a71130..b6ba506 100644
--- a/bin/noisechisel/noisechisel.c
+++ b/bin/noisechisel/noisechisel.c
@@ -71,6 +71,12 @@ noisechisel_convolve(struct noisechiselparams *p)
       gal_fits_img_write(p->input, p->detectionname, NULL, PROGRAM_STRING);
       gal_fits_img_write(p->conv, p->detectionname, NULL, PROGRAM_STRING);
     }
+
+  if(p->conv->ndim==3)
+    {
+      printf("\n... end of %s ...\n", __func__);
+      exit(0);
+    }
 }
 
 
diff --git a/bin/noisechisel/ui.c b/bin/noisechisel/ui.c
index 3165f33..7c53eee 100644
--- a/bin/noisechisel/ui.c
+++ b/bin/noisechisel/ui.c
@@ -391,29 +391,78 @@ ui_set_output_names(struct noisechiselparams *p)
 
 
 
+/* Prepare the kernel, either from a file, or from the default arrays
+   (depending on the dimensionality). The default kernels were created as
+   follows:
+
+   2D: The following commands were put in a script and run. It will create
+   a plain text array along with a FITS image.
+
+       set -o errexit           # Stop if a program returns false.
+       echo "0 0.0 0.0 3 2 0 0 1 1 5" > tmp.txt
+       export GSL_RNG_TYPE=ranlxs2
+       export GSL_RNG_SEED=1
+       astmkprof tmp.txt --oversample=1 --envseed --numrandom=10000 \
+                 --tolerance=0.01 --nomerged
+       astcrop 0_tmp.fits --section=2:*-1,2:*-1 --zeroisnotblank    \
+               --output=fwhm2.fits
+       astconvertt fwhm2.fits --output=fwhm2.txt
+       rm 0_tmp.fits tmp.txt
+
+   3D: MakeProfiles was initially run with these commands:
+
+       $ export GSL_RNG_SEED=1
+       $ export GSL_RNG_TYPE=ranlxs2
+       $ astmkprof --kernel=gaussian-3d,2,5,0.5 --oversample=1 --envseed
+
+   The resulting fits file was converted to text with the this C program to
+   generate the `kernel-3d.h' header file (until ConvertType supports 3D
+   datasets) which is then "include"d into this function.
+
+       #include <stdio.h>
+       #include <stdlib.h>
+       #include <gnuastro/fits.h>
+
+       int
+       main(void)
+       {
+         size_t i;
+         float *arr;
+         gal_data_t *img=gal_fits_img_read_to_type("kernel.fits", "1",
+                                                   GAL_TYPE_FLOAT32, -1);
+
+         arr=img->array;
+
+         printf("size_t kernel_3d_dsize[3]={%zu, %zu, %zu};\n",
+                img->dsize[0], img->dsize[1], img->dsize[2]);
+         printf("float kernel_3d[%zu]={", img->size);
+         for(i=0;i<img->size;++i)
+           {
+             if(i>0)
+               {
+                 if(i % img->dsize[2]                 == 0 ) printf("\n");
+                if(i % (img->dsize[2]*img->dsize[1]) == 0 ) printf("\n");
+              }
+             printf(i==(img->size-1) ? "%.7g": "%.7g, ", arr[i]);
+           }
+         printf("};\n");
+
+         return EXIT_SUCCESS;
+       }
+
+   Assuming this C program is in a file named `kernel.c', it can be
+   compiled, run and saved into `kernel-3d.h' with the command below:
+
+       $ astbuildprog -q kernel.c > kernel-3d.h
+*/
 static void
 ui_prepare_kernel(struct noisechiselparams *p)
 {
-  /* The default kernel. It was created by saving the following commands in a
-     script and running it. It will create a plain text array along with a
-     FITS image. The crop is because the first and last rows and columns of
-     all PSFs made by MakeProfiles are blank (zero) when you run with
-     oversample=1. You can keep the spaces when copying and pasting ;-). Just
-     make it executable and run it.
-
-     set -o errexit           # Stop if a program returns false.
-     echo "0 0.0 0.0 3 2 0 0 1 1 5" > tmp.txt
-     export GSL_RNG_TYPE=ranlxs2
-     export GSL_RNG_SEED=1
-     astmkprof tmp.txt --oversample=1 --envseed --numrandom=10000 \
-     --tolerance=0.01 --nomerged
-     astcrop 0_tmp.fits --section=2:*-1,2:*-1 --zeroisnotblank    \
-     --output=fwhm2.fits
-     astconvertt fwhm2.fits --output=fwhm2.txt
-     rm 0_tmp.fits tmp.txt
-  */
+#include <kernel-3d.h>
+  float *f, *ff, *k;
+  size_t ndim=p->input->ndim;
   size_t kernel_2d_dsize[2]={11,11};
-  float *f, *ff, *k, kernel_2d[121]=
+  float  kernel_2d[121]=
     {
       0, 0, 0, 0, 0, 6.57699e-09, 0, 0, 0, 0, 0,
 
@@ -456,14 +505,14 @@ ui_prepare_kernel(struct noisechiselparams *p)
     {
       /* Allocate space for the kernel (we don't want to use the statically
          allocated array. */
-      p->kernel=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, 2,
-                               kernel_2d_dsize, NULL, 0, p->cp.minmapsize,
-                               NULL, NULL, NULL);
-
-      /* Now copy the staticly allocated array into it. */
-      k=p->kernel->array;
-      ff=(f=kernel_2d)+gal_dimension_total_size(2, p->kernel->dsize);
-      do *k++=*f; while(++f<ff);
+      p->kernel=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, p->input->ndim,
+                               ndim==2 ? kernel_2d_dsize : kernel_3d_dsize,
+                               NULL, 0, p->cp.minmapsize, NULL, NULL, NULL);
+
+      /* Copy the staticly allocated default array into `p->kernel'. */
+      k = ndim==2 ? kernel_2d : kernel_3d;
+      ff = (f=p->kernel->array) + p->kernel->size;
+      do *f=*k++; while(++f<ff);
     }
 }
 
@@ -566,13 +615,6 @@ ui_preparations(struct noisechiselparams *p)
     gal_checkset_allocate_copy("INPUT", &p->input->name);
 
 
-  /* NoiseChisel currently only works on 2D datasets (images). */
-  if(p->input->ndim!=2)
-    error(EXIT_FAILURE, 0, "%s (hdu: %s) has %zu dimensions but NoiseChisel "
-          "can only operate on 2D datasets (images)", p->inputname, p->cp.hdu,
-          p->input->ndim);
-
-
   /* Check for blank values to help later processing. AFTERWARDS, set the
      USE_ZERO flag, so the zero-bit (if the input doesn't have any blank
      value) will be meaningful. */
@@ -681,7 +723,10 @@ ui_read_check_inputs_setup(int argc, char *argv[],
       if(p->kernelname)
         printf("  - Kernel: %s (hdu: %s)\n", p->kernelname, p->khdu);
       else
-        printf("  - Kernel: FWHM=2 pixel Gaussian.\n");
+        printf(p->kernel->ndim==2
+               ? "  - Kernel: FWHM=2 pixel Gaussian.\n"
+               : "  - Kernel: FWHM=2 pixel Gaussian (half extent in "
+               "3rd dimension).\n");
     }
 }
 
diff --git a/bin/warp/ui.c b/bin/warp/ui.c
index 0e2df12..3ecbf1a 100644
--- a/bin/warp/ui.c
+++ b/bin/warp/ui.c
@@ -518,7 +518,7 @@ ui_matrix_make_align(struct warpparams *p, double *tmatrix)
       x[0] = w[0]<0 ? 1.0f : -1.0f;  /* Has to be negative. */
       x[1] = 0.0f;
       x[2] = 0.0f;
-      x[3] = w[3]>0 ? 1.0f : 1.0f;   /* Has to be positive. */
+      x[3] = w[3]>0 ? 1.0f : -1.0f;  /* Has to be positive. */
     }
   else
     {
diff --git a/lib/convolve.c b/lib/convolve.c
index 847671b..48287ff 100644
--- a/lib/convolve.c
+++ b/lib/convolve.c
@@ -1,5 +1,5 @@
 /*********************************************************************
-convolve -- Convolve a dataset with a given kernel.
+Convolve -- Convolve a dataset with a given kernel.
 This is part of GNU Astronomy Utilities (Gnuastro) package.
 
 Original author:
@@ -91,8 +91,10 @@ convolve_tile_is_on_edge(size_t *h, size_t *start_end_coord, 
size_t *k,
 struct per_thread_spatial_prm
 {
   /* Internally stored/used values. */
+  size_t             id;     /* ID of tile being operatred on.           */
   gal_data_t      *tile;     /* Tile this thread is working on.          */
-  gal_data_t   *overlap;     /* Overlap pointer and starting point.      */
+  gal_data_t *i_overlap;     /* Overlap tile over input dataset.         */
+  gal_data_t *k_overlap;     /* Overlap tile over kernel dataset.        */
   size_t *overlap_start;     /* Starting coordinate of kernel overlap.   */
   size_t  *kernel_start;     /* Kernel starting point.                   */
   size_t    *host_start;     /* Starting coordinate of host.             */
@@ -100,7 +102,6 @@ struct per_thread_spatial_prm
                                 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.*/
 };
 
@@ -124,23 +125,24 @@ 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).  */
 static int
-convolve_spatial_overlap(struct per_thread_spatial_prm *pprm,
-                         int tocorrect)
+convolve_spatial_overlap(struct per_thread_spatial_prm *pprm, int tocorrect)
 {
   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 ndim=block->ndim;
 
-  size_t *pp, *ppf, *hs;
-  size_t overlap_inc, ndim=block->ndim;
+  size_t *kd=pprm->k_overlap->dsize;
+  size_t *pp, *ppf, *hs, increment, size=1;
   size_t *h=dsize, *os=pprm->overlap_start;
+  size_t *p, *pf, *od=pprm->i_overlap->dsize;
   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
@@ -153,14 +155,15 @@ convolve_spatial_overlap(struct per_thread_spatial_prm 
*pprm,
 
 
   /* Coordinate to start convolution for this pixel. */
+  pf = (p=pprm->pix) + ndim;
   do
     {
       /* Initialize the overlap for this dimension (we'll assume it
          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
+      /* When the tile is on the edge, 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. 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. */
@@ -206,6 +209,12 @@ convolve_spatial_overlap(struct per_thread_spatial_prm 
*pprm,
               if(is_start) *od -= *k/2 - *p;
               if(is_end)   *od -= *p + *k/2 - *h + 1;
 
+              /* Put the overlap size into the kernel's overlap `dsize'
+                 also and then use it to update the total size of the
+                 overlap. */
+              *kd++ = *od;
+              size *= *od;
+
               /* Increment and finalize. */
               ++h;
               ++k;
@@ -221,38 +230,22 @@ convolve_spatial_overlap(struct per_thread_spatial_prm 
*pprm,
         {
           /* Set the values. */
           *ks++ = 0;
-          *od++ = *k;
+          size *= *k;
+          *kd++ = *od = *k;
           *os++ = *p - *k/2;
 
           /* Increment. */
           ++h;
           ++k;
+          ++od;
         }
     }
   while(++p<pf);
 
 
-  /* 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("\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);
-    }
-  */
+  /* Update the `size' element of both overlap datasets. */
+  pprm->i_overlap->size = pprm->k_overlap->size = size;
 
-  /* Set the increment to start working on the kernel. */
-  pprm->k_start_inc = ( full_overlap
-                        ? 0
-                        : gal_dimension_coord_to_index(ndim,
-                                                       kernel->dsize,
-                                                       pprm->kernel_start) );
 
   /* Make correction.
 
@@ -270,17 +263,26 @@ convolve_spatial_overlap(struct per_thread_spatial_prm 
*pprm,
          the channel). */
   hs=pprm->host_start;
   if(tocorrect)
-    { ppf=(pp=pprm->pix)+ndim; do *pp -= *hs++; while(++pp<ppf); }
+    { ppf=(pp=pprm->pix)           + ndim; do *pp -= *hs++; while(++pp<ppf); }
   else
-    { ppf=(pp=pprm->overlap_start)+ndim; do *pp += *hs++; while(++pp<ppf); }
+    { ppf=(pp=pprm->overlap_start) + ndim; do *pp += *hs++; while(++pp<ppf); }
+
+
+  /* Set the starting point of the dataset overlap tile. */
+  increment=gal_dimension_coord_to_index(ndim, block->dsize,
+                                         pprm->overlap_start);
+  pprm->i_overlap->array=gal_data_ptr_increment(block->array, increment,
+                                                block->type);
 
 
-  /* 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,
-                                           pprm->overlap_start);
-  pprm->overlap->array=gal_data_ptr_increment(block->array, overlap_inc,
-                                              block->type);
+  /* Set the starting point of the kernel overlap tile. */
+  increment = ( full_overlap
+                ? 0
+                : gal_dimension_coord_to_index(ndim,
+                                               kernel->dsize,
+                                               pprm->kernel_start) );
+  pprm->k_overlap->array=gal_data_ptr_increment(kernel->array, increment,
+                                                kernel->type);
   return full_overlap;
 }
 
@@ -292,35 +294,37 @@ convolve_spatial_overlap(struct per_thread_spatial_prm 
*pprm,
 static void
 convolve_spatial_tile(struct per_thread_spatial_prm *pprm)
 {
+  gal_data_t *tile=pprm->tile;
 
-  double sum, ksum;
   int full_overlap;
+  double sum, ksum;
   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];
+  size_t j, ndim=block->ndim, csize=tile->dsize[ndim-1];
+  gal_data_t *i_overlap=pprm->i_overlap, *k_overlap=pprm->k_overlap;
 
   /* Variables for scanning a tile (`i_*') and the region around every
      pixel of a tile (`o_*'). */
   size_t start_fastdim;
-  size_t k_inc, i_inc, i_ninc, i_st_en[2], o_inc, o_ninc, o_st_en[2];
+  size_t i_inc, i_ninc, i_st_en[2];
 
   /* These variables depend on the type of the input. */
-  float *kv, *iv, *ivf, *i_start, *o_start, *k_start;
+  float *i_start;
   float *in_v, *in=block->array, *out=cprm->out->array;
 
 
-  /* Starting pixel for this tile. Note that when we are in `tocorrect'
-     mode, this position has to be  */
+  /* Starting pixel for the host of this tile. Note that when we are in
+     `convoverch' mode, `host' refers to the fully allocated block of
+     memory. */
   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
-     start and end are the first two allocated spaces in
-     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'. */
+     space for the start and end coordinates is stored in `p->pix'). 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, pprm->pix, cprm->convoverch);
   start_fastdim = pprm->pix[ndim-1];
 
@@ -334,15 +338,8 @@ convolve_spatial_tile(struct per_thread_spatial_prm *pprm)
      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", pprm->pix[0], pprm->pix[1]);
-      exit(0);
-    }
-  */
-  /* Go over the tile. */
+
+  /* Parse over all the tile elements. */
   i_inc=0; i_ninc=1;
   i_start=gal_tile_start_end_ind_inclusive(tile, block, i_st_en);
   while( i_st_en[0] + i_inc <= i_st_en[1] )
@@ -352,10 +349,10 @@ convolve_spatial_tile(struct per_thread_spatial_prm *pprm)
       pprm->pix[ndim-1]=start_fastdim;
 
       /* Go over each pixel to convolve. */
-      for(i=0;i<csize;++i)
+      for(j=0;j<csize;++j)
         {
           /* Pointer to the pixel under consideration. */
-          in_v = i_start + i_inc + i;
+          in_v = i_start + i_inc + j;
 
           /* If the input on this pixel is a NaN, then just set the output
              to NaN too and go onto the next pixel. `in_v' is the pointer
@@ -378,54 +375,21 @@ convolve_spatial_tile(struct per_thread_spatial_prm *pprm)
                   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] )
-                    {
-                      /* 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
+                  /* Initialize the necessary values. */
+                  sum  = 0.0L;
+                  ksum = cprm->edgecorrection ? 0.0L : 1.0L;
+
+                  /* Parse over both the overlap tiles. */
+                  GAL_TILE_PO_OISET(float, float, i_overlap, k_overlap, 1, 0, {
+                      if( !isnan(*i) )
                         {
-                          if( !isnan(*iv) )
-                            {
-                              sum += *iv * *kv;
-                              if(cprm->edgecorrection) ksum += *kv;
-                            }
-                          ++kv;
+                          sum += *i * *o;
+                          if(cprm->edgecorrection) ksum += *o;
                         }
-                      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);
-                    }
+                    });
 
                   /* Set the output value. */
-                  out[ in_v - in ] = ( ksum==0.0f
-                                       ? NAN
-                                       : sum/ksum );
+                  out[ in_v - in ] = ksum==0.0L ? NAN : sum/ksum;
                 }
             }
 
@@ -439,7 +403,7 @@ convolve_spatial_tile(struct per_thread_spatial_prm *pprm)
                                         pprm->pix);
     }
   /*
-  if(tile_ind==2053)
+  if(pprm->id==2053)
     printf("... done.\n");
   */
 }
@@ -459,7 +423,7 @@ convolve_spatial_on_thread(void *inparam)
   size_t i;
   size_t ndim=block->ndim;
   struct per_thread_spatial_prm *pprm=&cprm->pprm[tprm->id];
-  size_t *dsize=gal_data_malloc_array(GAL_TYPE_SIZE_T,ndim, __func__,
+  size_t *dsize=gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim, __func__,
                                       "dsize");
 
 
@@ -469,47 +433,47 @@ convolve_spatial_on_thread(void *inparam)
 
 
   /* Initialize/Allocate necessary items for this thread. */
-  pprm->cprm           = cprm;
-  pprm->pix            = gal_data_malloc_array(GAL_TYPE_SIZE_T, 2*ndim,
-                                               __func__, "pprm->pix");
-  pprm->host_start     = gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim,
-                                               __func__, "pprm->host_start");
-  pprm->kernel_start   = gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim,
-                                               __func__,
-                                               "pprm->kernel_start");
-  pprm->overlap_start  = gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim,
+  pprm->cprm          = cprm;
+  pprm->pix           = gal_data_malloc_array(GAL_TYPE_SIZE_T, 2*ndim,
+                                              __func__, "pprm->pix");
+  pprm->host_start    = gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim,
+                                              __func__, "pprm->host_start");
+  pprm->kernel_start  = gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim,
+                                              __func__, "pprm->kernel_start");
+  pprm->overlap_start = gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim,
                                                __func__,
-                                               "pprm->overlap_start");
-  pprm->overlap        = gal_data_alloc(NULL, block->type, ndim, dsize,
-                                        NULL, 0, -1, NULL, NULL, NULL);
+                                              "pprm->overlap_start");
+  pprm->i_overlap     = gal_data_alloc(NULL, block->type, ndim, dsize,
+                                       NULL, 0, -1, NULL, NULL, NULL);
+  pprm->k_overlap     = gal_data_alloc(NULL, cprm->kernel->type, ndim, dsize,
+                                       NULL, 0, -1, NULL, NULL, NULL);
   free(dsize);
-  free(pprm->overlap->array);
-  pprm->overlap->block = cprm->block;
+  free(pprm->i_overlap->array);
+  free(pprm->k_overlap->array);
+  pprm->i_overlap->block = cprm->block;
+  pprm->k_overlap->block = cprm->kernel;
 
 
   /* Go over all the tiles given to this thread. */
   for(i=0; tprm->indexs[i] != GAL_BLANK_SIZE_T; ++i)
     {
       /* Set this tile's pointer into this thread's parameters. */
-      pprm->tile = &cprm->tiles[ tprm->indexs[i] ];
+      pprm->id   = tprm->indexs[i];
+      pprm->tile = &cprm->tiles[ pprm->id ];
 
       /* 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(pprm->pix);
   free(pprm->host_start);
   free(pprm->kernel_start);
   free(pprm->overlap_start);
-  gal_data_free(pprm->overlap);
+  gal_data_free(pprm->i_overlap);
+  gal_data_free(pprm->k_overlap);
   if(tprm->b) pthread_barrier_wait(tprm->b);
   return NULL;
 }
@@ -533,8 +497,7 @@ gal_convolve_spatial_general(gal_data_t *tiles, gal_data_t 
*kernel,
   if(tiles->ndim!=kernel->ndim)
     error(EXIT_FAILURE, 0, "%s: The number of dimensions between the kernel "
           "and input should be the same", __func__);
-  if( block->type!=GAL_TYPE_FLOAT32
-      || kernel->type!=GAL_TYPE_FLOAT32 )
+  if( block->type!=GAL_TYPE_FLOAT32 || kernel->type!=GAL_TYPE_FLOAT32 )
     error(EXIT_FAILURE, 0, "%s: only accepts `float32' type input and "
           "kernel currently", __func__);
 
diff --git a/lib/gnuastro/convolve.h b/lib/gnuastro/convolve.h
index 563baed..07ba554 100644
--- a/lib/gnuastro/convolve.h
+++ b/lib/gnuastro/convolve.h
@@ -1,5 +1,5 @@
 /*********************************************************************
-convolve -- Convolve the dataset with a given kernel.
+Convolve -- Convolve the dataset with a given kernel.
 This is part of GNU Astronomy Utilities (Gnuastro) package.
 
 Original author:
diff --git a/lib/gnuastro/tile.h b/lib/gnuastro/tile.h
index 2809b23..e33f1d8 100644
--- a/lib/gnuastro/tile.h
+++ b/lib/gnuastro/tile.h
@@ -218,9 +218,9 @@ gal_tile_full_free_contents(struct 
gal_tile_two_layer_params *tl);
         else                                                            \
           if(gal_data_dsize_is_different(IN, OTHER) )                   \
             {                                                           \
-              /* The `error' function, is a GNU extension and this is  */ \
-              /* a header, not a library which the user has to compile */ \
-              /* every time (on their own system).                     */ \
+              /* The `error' function, is a GNU extension and this */   \
+              /* is a header, not a library which the user has to  */   \
+              /* compile every time (on their own system).         */   \
               fprintf(stderr, "GAL_TILE_PO_OISET: when "                \
                       "`PARSE_OTHER' is non-zero, the sizes of `IN' "   \
                       "and `OTHER' must be equal (in all "              \
diff --git a/lib/tile.c b/lib/tile.c
index 5457375..cbff41e 100644
--- a/lib/tile.c
+++ b/lib/tile.c
@@ -329,8 +329,8 @@ gal_tile_block(gal_data_t *tile)
 
   num_increment      coord         increment
   -------------      -----         ---------
-         0          (...0,0,0)     b[n-1]: fastest dimension of the block.
-         1          (...0,1,0)     Similar to previous
+         1          (...0,0,0)     b[n-1]: fastest dimension of the block.
+         2          (...0,1,0)     Similar to previous
          .              .               .
          .              .               .
        t[n-2]       (...1,0,0)     (b[n-2] * b[n-1]) - ( (t[n-2]-1) * b[n-1] )
@@ -385,7 +385,7 @@ gal_tile_block_increment(gal_data_t *block, size_t *tsize,
       else
         {
           increment=(b[1] * b[2]) - ( (t[1]-1) * b[2] );
-          if(coord) { ++coord[0]; coord[1]=coord[2]=0; }
+          if(coord) { ++coord[0]; coord[1] -= t[1]-1; coord[2]=0; }
         }
       break;
     }