[gnuastro-commits] master 95ba655 100/125: Spatial convolution with new tessellation

[Mohammad Akhlaghi]

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

    Spatial convolution with new tessellation
    
    Spatial convolution has been fully implemented using the new tessellation
    system defined in `lib/tile.c'. It can also account for individual
    convolution on each channel as before, but not using any extra structure!
    It was fully done with `gal_data_t'.
---
 bin/convolve/args.h           |  12 +--
 bin/convolve/astconvolve.conf |   2 +-
 bin/convolve/convolve.c       |   5 +-
 bin/convolve/main.h           |   4 +-
 bin/convolve/ui.c             |  14 +--
 bin/convolve/ui.h             |   2 +-
 doc/gnuastro.texi             | 157 ++++++++++++++-------------
 lib/convolve.c                | 243 +++++++++++++++++++++++++++++++++---------
 lib/tile.c                    |  38 +++++--
 9 files changed, 323 insertions(+), 154 deletions(-)

diff --git a/bin/convolve/args.h b/bin/convolve/args.h
index ab310b4..98bd78e 100644
--- a/bin/convolve/args.h
+++ b/bin/convolve/args.h
@@ -122,17 +122,17 @@ struct argp_option program_options[] =
     /* Tile grid options. */
     {
       0, 0, 0, 0,
-      "Tile grid (only for spatial domain):",
+      "Tessellation (only for spatial domain):",
       ARGS_GROUP_MESH_GRID
     },
     {
-      "tile",
-      ARGS_OPTION_KEY_TILE,
+      "tilesize",
+      ARGS_OPTION_KEY_TILESIZE,
       "INT[,INT]",
       0,
-      "Size of tiles along each dim. (FITS order).",
+      "Size of regular tiles on each dim. (FITS order).",
       ARGS_GROUP_MESH_GRID,
-      &p->tile,
+      &p->tilesize,
       GAL_DATA_TYPE_SIZE_T,
       GAL_OPTIONS_RANGE_GT_0,
       GAL_OPTIONS_NOT_MANDATORY,
@@ -158,7 +158,7 @@ struct argp_option program_options[] =
       ARGS_OPTION_KEY_REMAINDERFRAC,
       "FLT",
       0,
-      "Fraction of remainers in each dim to chop.",
+      "Fraction of remainder to split last tile.",
       ARGS_GROUP_MESH_GRID,
       &p->remainderfrac,
       GAL_DATA_TYPE_FLOAT32,
diff --git a/bin/convolve/astconvolve.conf b/bin/convolve/astconvolve.conf
index 941e714..da3d0f5 100644
--- a/bin/convolve/astconvolve.conf
+++ b/bin/convolve/astconvolve.conf
@@ -25,7 +25,7 @@
  type               float32
 
 # Mesh grid:
- tile               30,30
+ tilesize           30,30
  numchannels        1,1
  remainderfrac      0.1
  convoverch         0
diff --git a/bin/convolve/convolve.c b/bin/convolve/convolve.c
index ff523d1..f179440 100644
--- a/bin/convolve/convolve.c
+++ b/bin/convolve/convolve.c
@@ -261,7 +261,7 @@ removepaddingcorrectroundoff(struct convolveparams *p)
   size_t ps1=p->ps1;
   size_t i, hi0, hi1;
   size_t *isize=p->input->dsize;
-  double *o, *input=p->input->array;
+  float *o, *input=p->input->array;
   double *d, *df, *start, *pimg=p->pimg;
 
   /* Set all the necessary parameters to crop the desired region. hi0
@@ -776,14 +776,13 @@ convolve(struct convolveparams *p)
   if(p->domain==CONVOLVE_DOMAIN_SPATIAL)
     {
       /* Prepare the mesh structure. */
-      gal_tile_all_position_two_layers(p->input, p->channel, p->tile,
+      gal_tile_all_position_two_layers(p->input, p->channelsize, p->tilesize,
                                        p->remainderfrac, &channels, &tiles);
 
       /* Save the tile IDs if they are requested. */
       if(p->tilesname)
         gal_tile_block_check_tiles(tiles, p->tilesname, PROGRAM_NAME);
 
-
       /* Do the spatial convolution. */
       out=gal_convolve_spatial(tiles, p->kernel, p->cp.numthreads,
                                p->convoverch);
diff --git a/bin/convolve/main.h b/bin/convolve/main.h
index ab8d157..4d08549 100644
--- a/bin/convolve/main.h
+++ b/bin/convolve/main.h
@@ -79,7 +79,7 @@ struct convolveparams
   uint8_t       nokernelnorm;  /* Do not normalize the kernel.            */
   double        minsharpspec;  /* Deconvolution: min spect. of sharp img. */
   uint8_t     checkfreqsteps;  /* View the frequency domain steps.        */
-  size_t               *tile;  /* Size of tiles along each dim. (C order).*/
+  size_t           *tilesize;  /* Size of tiles along each dim. (C order).*/
   size_t        *numchannels;  /* No. of tiles along each dim. (C order). */
   float        remainderfrac;  /* Frac. of remainers in each dim to cut.  */
   uint8_t         convoverch;  /* Convolve over channel borders.          */
@@ -88,7 +88,7 @@ struct convolveparams
   size_t          makekernel;  /* Make a kernel to create input.          */
 
   /* Internal */
-  size_t            *channel;  /* Size of channels along each dimension.  */
+  size_t        *channelsize;  /* Size of channels along each dimension.  */
   int                 domain;  /* Frequency or spatial domain conv.       */
   gal_data_t          *input;  /* Input image array.                      */
   gal_data_t         *kernel;  /* Input Kernel array.                     */
diff --git a/bin/convolve/ui.c b/bin/convolve/ui.c
index 513c3ee..eb51aff 100644
--- a/bin/convolve/ui.c
+++ b/bin/convolve/ui.c
@@ -233,16 +233,16 @@ ui_read_check_only_options(struct convolveparams *p)
   /* If we are in the spatial domain, make sure that the necessary
      parameters are set. */
   if( p->domain==CONVOLVE_DOMAIN_SPATIAL )
-    if( p->tile==NULL || p->numchannels==NULL )
+    if( p->tilesize==NULL || p->numchannels==NULL )
       {
-        if( p->tile==NULL && p->numchannels==NULL )
+        if( p->tilesize==NULL && p->numchannels==NULL )
           error(EXIT_FAILURE, 0, "in spatial convolution, `--numchannels' "
-                "and `--tile' are mandatory");
+                "and `--tilesize' are mandatory");
         else
           error(EXIT_FAILURE, 0, "in spatial convolution, `--%s' is "
                 "mandatory: you should use it to set the %s",
-                p->tile ? "numchannels" : "tile",
-                ( p->tile
+                p->tilesize ? "numchannels" : "tilesize",
+                ( p->tilesize
                   ? "number of channels along each dimension of the input"
                   : "size of tiles to cover the input along each "
                   "dimension" ) );
@@ -360,8 +360,8 @@ ui_preparations(struct convolveparams *p)
                 PROGRAM_NAME);
     }
   else
-    p->channel=gal_tile_all_sanity_check(p->filename, p->cp.hdu, p->input,
-                                         p->tile, p->numchannels);
+    p->channelsize=gal_tile_all_sanity_check(p->filename, p->cp.hdu, p->input,
+                                             p->tilesize, p->numchannels);
 
 
 
diff --git a/bin/convolve/ui.h b/bin/convolve/ui.h
index a4df9f5..53413ad 100644
--- a/bin/convolve/ui.h
+++ b/bin/convolve/ui.h
@@ -39,7 +39,7 @@ enum option_keys_enum
   ARGS_OPTION_KEY_KHDU           = 'u',
   ARGS_OPTION_KEY_MINSHARPSPEC   = 'H',
   ARGS_OPTION_KEY_CHECKFREQSTEPS = 'C',
-  ARGS_OPTION_KEY_TILE           = 't',
+  ARGS_OPTION_KEY_TILESIZE       = 't',
   ARGS_OPTION_KEY_NUMCHANNELS    = 'c',
   ARGS_OPTION_KEY_REMAINDERFRAC  = 'r',
   ARGS_OPTION_KEY_DOMAIN         = 'd',
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 2623302..9d01695 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -10793,26 +10793,57 @@ are ultimately poor approximations.
 @cindex Tile
 @cindex Gradient
 @cindex Mesh grid
-Some of the programs in Gnuastro will need to divide the pixels in an
-image into individual tiles or a mesh grid to be able to deal with
-gradients. In this section we will explain the concept in detail and
-how the user can check the grid. In the case of SubtractSky, if the
-image is completely uniform then one sky value will suffice for the
-whole image. See @ref{Sky value} for the definition of the sky value.
-Unfortunately though, as discussed in @ref{SubtractSky}, in most
-images taken with ground or space-based telescopes the sky value is
-not uniform. So we have to break the image up into small tiles on a
-mesh grid, assume the sky value is constant over them and find the sky
-value on those tiles.
-
-Meshes are considered to be a square with a side of
address@hidden pixels. The best mesh size is directly related to
-the gradient on the image. In practice we assume that the gradient is
-not significant over each mesh. So if there is a strong gradient (for
-example in long wavelength ground based images) or the image is of a
-crowded area where there isn't too much blank area, you have to choose
-a smaller mesh size. A larger mesh will give more pixels and so the
-scatter in the results will be less.
address@hidden Tile grid
address@hidden Tessellation
+Some of the programs in Gnuastro will need to divide the pixels in an image
+into a grid of individual, non-overlapping tiles. In the arts and
+mathematics, this process is formally known as
address@hidden://en.wikipedia.org/wiki/Tessellation, tessellation}. Defining a
+tessellation over a dataset (image) is necessary in many contexts. Their
+most common application is to account for local variations over the dataset
+(for example background sky gradients, which are more common as we go to
+longer/redder wavelengths/filters).
+
+In this section we will explain the concept in detail and how you can check
+the grid. In the case of SubtractSky, if the image is completely uniform
+then one sky value will suffice for the whole image. See @ref{Sky value}
+for the definition of the sky value.  Unfortunately though, as discussed in
address@hidden, in most images taken with ground or space-based
+telescopes the sky value is not uniform. So we have to break the image up
+into small tiles on a grid, assume the Sky value is constant over them and
+find the sky value on those tiles.
+
+The size of the regular tiles (in units of data-elements, or pixels in an
+image) can be defined with the @option{--tilesize} option. It takes
+multiple numbers (separated by a comma) which will be the length along the
+respective dimension (in Fortran/FITS order). For example
address@hidden,40} can be used for an image (a 2D dataset). The
+regular tile size along the first FITS axis (horizontal when viewed in SAO
+ds9) will be 30 pixels and along the second it will be 40 pixels. Ideally,
address@hidden should be selected such that all tiles in the image
+have exactly the same size (in other words, that the dataset length in each
+dimension is divisible by the tile size in that dimension).
+
+However, this is not always possible: the dataset can be any size and every
+pixel in it is valuable. In such cases, Gnuastro will look at the
+significance of the remainder length, if it is not significant (for example
+one or two pixels), then it will just increase the size of the first tile
+in the respective dimension and allow the rest of the tiles to have the
+required size. When the remainder is significant (for example one pixel
+less than the size along that dimension), the remainder will be added to
+one regular tile's size and the large tile will be cut in half and put in
+the two ends of the grid/tessellation. In this way, all the tiles in the
+central regions of the dataset will have the regular tile sizes and the
+tiles on the edge will be slightly smaller. You can set the significance
+fraction with the @option{--remainderfrac} option.
+
+The best tile size is directly related to the spatial properties of the
+gradient on the image. In practice we assume that the gradient is not
+present over each tile. So if there is a strong gradient (for example in
+long wavelength ground based images) or the image is of a crowded area
+where there isn't too much blank area, you have to choose a smaller tile
+size. A larger mesh will give more pixels and so the scatter in the results
+will be less.
 
 @cindex CCD
 @cindex Amplifier
@@ -10820,69 +10851,45 @@ scatter in the results will be less.
 @cindex Subaru Telescope
 @cindex Hyper Suprime-Cam
 @cindex Hubble Space Telescope
-For raw image processing, a simple mesh grid is not sufficient. Raw
-images are the unprocessed outputs of the camera detectors. Large
+For raw image processing, a single tesselation/grid is not sufficient. Raw
+images are the unprocessed outputs of the camera detectors. Modern
 detectors usually have multiple readout channels each with its own
 amplifier. For example the Hubble Space Telescope Advanced Camera for
-Surveys (ACS) has four amplifiers over its full detector area dividing
-the square field of view to four smaller squares. Ground based image
-detectors are not exempt, for example each CCD of Subaru Telescope's
-Hyper Suprime-Cam camera (which has 104 CCDs) has four amplifiers, but
-they have the same height of the CCD and divide the width by four
-parts.
+Surveys (ACS) has four amplifiers over its full detector area dividing the
+square field of view to four smaller squares. Ground based image detectors
+are not exempt, for example each CCD of Subaru Telescope's Hyper
+Suprime-Cam camera (which has 104 CCDs) has four amplifiers, but they have
+the same height of the CCD and divide the width by four parts.
 
 @cindex Channel
-The bias current on each amplifier is different, and normally bias
-subtraction is not accurately done. So even after subtracting the measured
-bias current, you can usually still identify the boundaries of different
+The bias current on each amplifier is different, and initial bias
+subtraction is not perfect. So even after subtracting the measured bias
+current, you can usually still identify the boundaries of different
 amplifiers by eye. See Figure 11(a) in Akhlaghi and Ichikawa (2015) for an
 example. This results in the final reduced data to have non-uniform
 amplifier-shaped regions with higher or lower background flux values. Such
 systematic biases will then propagate to all subsequent measurements we do
 on the data (for example photometry and subsequent stellar mass and star
-formation rate measurements in the case of galaxies). Therefore an accurate
-sky subtraction routine should also be able to account for such biases.
-
-To get an accurate result, the mesh boundaries should be located
-exactly on the amplifier boundaries. Otherwise, some meshes will
-contain pixels that have been read from two or four different
-amplifiers. These meshes are going to give very biased results and the
-amplifier boundary will still be present after sky subtraction.  So we
-define `channel's. A channel is an independent mesh grid that covers
-one amplifier to ensure that the meshes do not pass the amplifier
-boundary. They can also be used in subsequent steps as the area used
-to identify nearby neighbors to interpolate and smooth the final grid,
-see @ref{Grid interpolation and smoothing}. The number of channels
-along each axis can be specified by the user at run time through the
-command-line @option{--nch1} and @option{--nch2} options or in the
-configuration files, see @ref{Configuration files}. The area of each
-channel will then be tiled by meshes of the given size and subsequent
-processing will be done on those meshes. If the image is processed or
-the detector only has one amplifier, you can set the number of
-channels in both axes to 1.
-
-Unlike the channel size, that has to be an exact multiple of the image
-size, the mesh size can be any number. If it is not an exact multiple of
-the image side, the last mesh (rightest, for the first FITS dimension, and
-highest for the second when viewed in SAO ds9) will have a different size
-than the rest. If the remainder of the image size divided by mesh size is
-larger than a certain fraction (value to @option{--lastmeshfrac}) of the
-mesh size along each axis, a new (smaller) mesh will be put there instead
-of a larger mesh. This is done to avoid the last mesh becoming too large
-compared to the other meshes in the grid. Generally, it is best practice to
-choose the mesh size such that the last mesh is only a few (negligible)
-pixels wider or thinner than the rest.
-
-The final mesh grid can be seen on the image with the
address@hidden option that is available to all programs which
-use the mesh grid for localized operations. When this option is
-called, a multi-extension FITS file with a @file{_mesh.fits} suffix
-will be created along with the outputs, see @ref{Automatic
-output}. The first extension will be the input image. For each mesh
-grid the image produces, there will be a subsequent extension. Each
-pixel in the grid extensions is labeled to the mesh that it is part
-of. You can flip through the extensions to check the mesh sizes and
-locations compared to the input image.
+formation rate measurements in the case of galaxies).
+
+Therefore an accurate analysis requires a lower-level tesselation: with
+larger tiles, each large tile covering one amplifier channel. For clarity
+we'll call these larger tiles ``channels''. The number of channels along
+each dimension is defined through the @option{--numchannels}. Each channel
+will then be covered by its own individual smaller tessellation (with tile
+sizes determined by the @option{--tilesize} option). This will allow two
+adjacent pixels from different channels to never be considered together in
+the higher-level analysis. If the image is processed or the detector only
+has one amplifier, you can set the number of channels in both axes to 1.
+
+The final tesselation can be seen on the image with the
address@hidden option that is available to all programs which use
+tessellation for localized operations. When this option is called, a FITS
+file with a @file{_mesh.fits} suffix will be created along with the
+outputs, see @ref{Automatic output}. Each pixel in this image has the ID of
+the tile that covers it. If the number of channels in any dimension are
+larger than unity, you will notice that the tile IDs are defined such that
+the first channels is covered first, then the second and so on.
 
 @menu
 * Quantifying signal in a mesh::  Good meshs for gradients.
@@ -11124,7 +11131,7 @@ the `Invoking ProgramName' section within the list of 
options.
 @table @option
 
 @item -s INT
address@hidden --meshsize=INT
address@hidden --tilesize=INT
 The size of each mesh, see @ref{Tiling an image}. If the width of all
 channels are not an exact multiple of the specified size, then the last
 mesh on each axis will have a different size to cover the full channel.
diff --git a/lib/convolve.c b/lib/convolve.c
index 47cc530..81bf27e 100644
--- a/lib/convolve.c
+++ b/lib/convolve.c
@@ -56,9 +56,9 @@ convolve_tile_is_on_edge(size_t *h, size_t *start_end_coord, 
size_t *k,
      length along that dimension, then the tile is on the edge.
 
      If the end of the tile in this dimension added with the half-kernel
-     length along that dimension is larger than the host's size, then the
-     tile is on the edge. */
-  do if( (*s++ < *k/2) || (*e++ + *k/2 > *h++) ) return 1; while(++k<kf);
+     length along that dimension is equal to or larger than the host's
+     size, then the tile is on the edge. */
+  do if( (*s++ < *k/2) || (*e++ + *k/2 >= *h++) ) return 1; while(++k<kf);
 
   /* If control reaches here, then this is a central tile. */
   return 0;
@@ -97,6 +97,143 @@ 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. */
+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)
+{
+  size_t overlap_inc, ndim=kernel->ndim;
+  int full_overlap=1, dim_full_overlap, is_start, is_end;
+
+  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
+    {
+      /* 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
+         thing we do for the tiles that don't overlap for them. */
+      if(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;
+          if( is_start || is_end )
+            {
+              /* Overlapping with the outside.
+
+                 With the start: assume that in this dimension, the pixel
+                 is at position 2, while the kernel is 11 pixels wide (or 5
+                 pixels in half-width). As seen below, the kernel should
+                 start from pixel `5-2=3' in this dimension and the overlap
+                 size should decrease by the same amount.
+
+                    image:            0 1 2 3 4 5 6 7 8 9 ...
+                    pixel:                p
+                    kernel:     0 1 2 3 4 5 6 7 8 9 10
+
+                 With the end: Similar to above, but assume the pixel is
+                 two pixels away from the edge of a 100-pixel image. We are
+                 no longer worried about the overlap or kernel starting
+                 point, it is the width that we need to decrease it by:
+
+                   97 + 5 - 100 + 1 : The `1' is because we want the pixel
+                                      immediately after the end.
+
+                    image:        ... 92 93 94 95 96 97 98 99 | 100 101 102
+                    pixel:                            p
+                    kernel:                 0 1 2 3 4 5 6 7 8 | 9   10    */
+              *ks++ = is_start ? *k/2 - *p : 0;
+              *os++ = is_start ? 0         : *p - *k/2;
+
+              /* We will start with the full kernel width, then decrease it
+                 if the pixel is too close to the start or end along this
+                 dimension. Note that the host array/image might actually
+                 be smaller than kernel, so both cases might occur. */
+              *od = *k;
+              if(is_start) *od -= *k/2 - *p;
+              if(is_end)   *od -= *p + *k/2 - *h + 1;
+
+              /* Increment and finalize. */
+              ++h;
+              ++k;
+              ++od;
+              full_overlap=0;
+              dim_full_overlap=0;
+            }
+        }
+
+      /* There is full overlap for this pixel or tile over this
+         dimension.  */
+      if(dim_full_overlap)
+        {
+          /* Set the values. */
+          *ks++ = 0;
+          *od++ = *k;
+          *os++ = *p - *k/2;
+
+          /* Increment. */
+          ++h;
+          ++k;
+        }
+    }
+  while(++p<pf);
+  /*
+  if(tile_ind==2053)
+    {
+      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);
+    }
+  */
+  /* Set the increment to start working on the kernel. */
+  *k_start_inc = ( full_overlap
+                   ? 0
+                   : gal_multidim_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_multidim_add_coords(overlap_start, host_start, overlap_start, ndim);
+
+  /* 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_multidim_coord_to_index(ndim, block->dsize, overlap_start);
+  overlap->array=gal_data_ptr_increment(block->array, overlap_inc,
+                                        block->type);
+  return full_overlap;
+}
+
+
+
+
 
 /* Convolve over one tile that is not touching the edge. */
 static void
@@ -104,17 +241,20 @@ convolve_spatial_tile(struct spatial_params *cprm, size_t 
tile_ind,
                       size_t *parse_coords, gal_data_t *overlap)
 {
   gal_data_t *tile=&cprm->tiles[tile_ind];
+  gal_data_t *block=gal_tile_block(tile), *kernel=cprm->kernel;
 
-  int on_edge;
   double sum, ksum;
+  int on_edge, full_overlap;
   size_t i, ndim=tile->ndim, csize=tile->dsize[ndim-1];
-  gal_data_t *block=gal_tile_block(tile), *kernel=cprm->kernel;
+  gal_data_t *host=cprm->convoverch ? block : tile->block;
 
-  size_t i_inc, i_ninc, i_st_en[2], o_inc, o_ninc, o_st_en[2];
-  size_t *p, *pf, *k, *o, *e, o_st_inc, k_start, start_fastdim;
+  /* 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 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;
+  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
@@ -122,10 +262,9 @@ convolve_spatial_tile(struct spatial_params *cprm, size_t 
tile_ind,
      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;
-  size_t *end           = parse_coords + ndim;
-  size_t *overlap_start = parse_coords + ndim * 2;
-  size_t *o_c           = parse_coords + ndim * 3;
+  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;
 
 
   /* Set the starting and ending coordinates of this tile (recall that the
@@ -136,16 +275,21 @@ convolve_spatial_tile(struct spatial_params *cprm, size_t 
tile_ind,
   gal_tile_start_end_coord(tile, parse_coords, cprm->convoverch);
   start_fastdim = 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( ( cprm->convoverch
-                                      ? block->dsize
-                                      : tile->block->dsize ),
-                                    parse_coords, kernel->dsize, ndim);
-
-
-  /* Go over the tile, first find its limits, then loop over the fastest
-     dimension. */
+  on_edge=convolve_tile_is_on_edge(host->dsize, parse_coords, kernel->dsize,
+                                   ndim);
+  /*
+  if(tile_ind==2053)
+    {
+      printf("\ntile %zu...\n", tile_ind);
+      printf("\tpix: %zu, %zu\n", pix[0], pix[1]);
+      exit(0);
+    }
+  */
+  /* Go over the tile. */
   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] )
@@ -167,38 +311,32 @@ convolve_spatial_tile(struct spatial_params *cprm, size_t 
tile_ind,
             out[ in_v - in ]=NAN;
           else
             {
-              /* Coordinate to start convolution for this pixel. */
-              pf=(p=pix)+ndim; e=end;
-              o=overlap_start; k=kernel->dsize;
-              do
-                {
-                  if(on_edge)
-                    {
-                      return;
-                    }
-                  else { *o++ = *p - *k++/2; k_start=0; }
-                }
-              while(++p<pf);
-              o_st_inc=gal_multidim_coord_to_index(ndim, block->dsize,
-                                                   overlap_start);
-
-              /* Set the overlap parameters, then set the region. */
-              overlap->dsize=kernel->dsize;
-              overlap->array=gal_data_ptr_increment(block->array, o_st_inc,
-                                                    block->type);
+              /* 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. */
-              kv = ( k_start
-                     ? gal_data_ptr_increment(kernel->array, k_start,
-                                              kernel->type)
-                     : kernel->array );
-              ksum=sum=0.0f; o_inc=0; o_ninc=1;
+              ksum=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. */
+                  /* 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( !isnan(*iv) )
@@ -208,9 +346,13 @@ convolve_spatial_tile(struct spatial_params *cprm, size_t 
tile_ind,
                   while(++iv<ivf);
 
                   /* Update the incrementation to the next contiguous
-                     region of memory over this tile. */
+                     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. */
@@ -225,6 +367,10 @@ convolve_spatial_tile(struct spatial_params *cprm, size_t 
tile_ind,
          contiguous patch. */
       i_inc += gal_tile_block_increment(block, tile->dsize, i_ninc++, pix);
     }
+  /*
+  if(tile_ind==2053)
+    printf("... done.\n");
+  */
 }
 
 
@@ -241,7 +387,7 @@ convolve_spatial_on_thread(void *inparam)
   size_t i;
   gal_data_t overlap, *block=gal_tile_block(cprm->tiles);
   size_t *parse_coords=gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T,
-                                             4*block->ndim);
+                                             6*block->ndim);
 
   /* Initialize the necessary overlap parameters. */
   overlap.block=block;
@@ -256,6 +402,8 @@ convolve_spatial_on_thread(void *inparam)
 
   /* Clean up, wait until all other threads finish, then return. In a
      single thread situation, `tprm->b==NULL'. */
+  free(parse_coords);
+  free(overlap.dsize);
   if(tprm->b) pthread_barrier_wait(tprm->b);
   return NULL;
 }
@@ -286,7 +434,6 @@ gal_convolve_spatial(gal_data_t *tiles, gal_data_t *kernel,
   out=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT32, block->ndim, block->dsize,
                      block->wcs, 0, block->minmapsize, "CONVOLVED",
                      block->unit, NULL);
-  { float *f=out->array, *ff=f+out->size; do *f=NAN; while(++f<ff); }
 
   /* Set the pointers in the parameters structure. */
   params.out=out;
diff --git a/lib/tile.c b/lib/tile.c
index b2cad46..1fceecf 100644
--- a/lib/tile.c
+++ b/lib/tile.c
@@ -89,21 +89,36 @@ gal_tile_start_coord(gal_data_t *tile, size_t *start_coord)
 void
 gal_tile_start_end_coord(gal_data_t *tile, size_t *start_end, int rel_block)
 {
-  size_t start_ind, ndim=tile->ndim;
+  size_t *s, *sf, *h;
   gal_data_t *block=gal_tile_block(tile);
   gal_data_t *host=rel_block ? block : tile->block;
+  size_t *hcoord, start_ind, ndim=tile->ndim, *end=start_end+ndim;
 
   /* Get the starting index. Note that for the type we need the allocated
      block dataset and can't rely on the tiles. */
-  start_ind=gal_data_ptr_dist(host->array, tile->array, block->type);
+  start_ind=gal_data_ptr_dist(block->array, tile->array, block->type);
 
-  /* Get the coordinates of the starting point. */
-  gal_multidim_index_to_coord(start_ind, tile->ndim, host->dsize,
-                              start_end);
+  /* Get the coordinates of the starting point relative to the allocated
+     block. */
+  gal_multidim_index_to_coord(start_ind, ndim, block->dsize, start_end);
+
+  /* When the host is different from the block, the tile's starting
+     position needs to be corrected. */
+  if(host!=block)
+    {
+      /* Get the host's starting coordinates. */
+      start_ind=gal_data_ptr_dist(block->array, host->array, block->type);
+
+      /* Temporarily put the host's coordinates in the place held for the
+         ending coordinates. */
+      hcoord=end;
+      gal_multidim_index_to_coord(start_ind, ndim, block->dsize, hcoord);
+      sf=(s=start_end)+ndim; h=hcoord; do *s++ -= *h++; while(s<sf);
+    }
 
   /* Add the dimensions of the tile to the starting coordinate. Note that
      the ending coordinates are stored immediately after the start.*/
-  gal_multidim_add_coords(start_end, tile->dsize, start_end+ndim, ndim);
+  gal_multidim_add_coords(start_end, tile->dsize, end, ndim);
 }
 
 
@@ -264,13 +279,13 @@ gal_tile_block_increment(gal_data_t *block, size_t *tsize,
     /* 1D: the increment is just the tile size. */
     case 1:
       increment=t[0];
-      coord[0]+=increment;
+      if(coord) coord[0]+=increment;
       break;
 
     /* 2D: the increment is the block's number of fastest axis pixels. */
     case 2:
       increment=b[1];
-      ++coord[0];
+      if(coord) ++coord[0];
       break;
 
     /* Higher dimensions. */
@@ -278,13 +293,13 @@ gal_tile_block_increment(gal_data_t *block, size_t *tsize,
       if(num_increment % t[n-2])
         {
           increment=b[n-1];
-          ++coord[n-2];
+          if(coord) ++coord[n-2];
         }
       else
         {
           increment=(b[n-2] * b[n-1]) - ( (t[n-2]-1) * b[n-1] );
           ++coord[n-3];
-          coord[n-2]=coord[n-1]=0;
+          if(coord) coord[n-2]=coord[n-1]=0;
         }
       break;
     }
@@ -422,7 +437,8 @@ gal_tile_all_sanity_check(char *filename, char *hdu, 
gal_data_t *input,
      the dataset's dimensions). */
   if(i!=input->ndim)
     error(EXIT_FAILURE, 0, "%s (hdu: %s): has %zu dimensions, but only %zu "
-          "value(s) given for the tile size", filename, hdu, input->ndim, i);
+          "value(s) given for the tile size (`--tilesize' option).",
+          filename, hdu, input->ndim, i);
 
 
   /* Check the channel's dimensions. */