[gnuastro-commits] master 31b3b76 099/125: Spatial convolution now using new tessellation

[Mohammad Akhlaghi]

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

    Spatial convolution now using new tessellation
    
    Spatial convolution (on multiple threads) is one of the most basic and
    common applications of the tessellation management that was designed in the
    previous commits. With this commit convolution on central tiles of the
    image has been applied. Only the edge tiles remain. The new system has
    significantly simplified the code (and probably the efficiency) of spatial
    convolution in Gnuastro.
    
    Also, reading of the input image has now been changed to float32
    type. Previously it was float64, while the kernel was initially read as
    float64. For frequency domain convolution this makes no difference, because
    the input is put into a complex float64 type array before the processing is
    done. For spatial domain convolution it also makes no difference because
    the sum is stored as float64.
    
    Finally, as before with the `mesh' library, a new
    `gal_tile_function_on_threads' function has been defined to make it easy to
    spin off threads working on separate tiles.
---
 bin/convolve/convolve.c |  17 ++-
 bin/convolve/ui.c       |   2 +-
 lib/convolve.c          | 258 +++++++++++++++++++++++++++++++++++++++++-
 lib/gnuastro/convolve.h |   5 +-
 lib/gnuastro/tile.h     |  41 +++++--
 lib/tile.c              | 291 ++++++++++++++++++++++++++++++++++++++----------
 6 files changed, 535 insertions(+), 79 deletions(-)

diff --git a/bin/convolve/convolve.c b/bin/convolve/convolve.c
index 3c7d1c0..ff523d1 100644
--- a/bin/convolve/convolve.c
+++ b/bin/convolve/convolve.c
@@ -190,10 +190,10 @@ void
 makepaddedcomplex(struct convolveparams *p)
 {
   size_t i, ps0, ps1;
-  float *f, *ff, *kernel=p->kernel->array;
+  double *o, *op, *pimg, *pker;
   size_t is0=p->input->dsize[0],  is1=p->input->dsize[1];
   size_t ks0=p->kernel->dsize[0], ks1=p->kernel->dsize[1];
-  double *o, *op, *d, *df, *pimg, *pker, *input=p->input->array;
+  float *f, *ff, *input=p->input->array, *kernel=p->kernel->array;
 
 
   /* Find the sizes of the padded image, note that since the kernel
@@ -220,8 +220,8 @@ makepaddedcomplex(struct convolveparams *p)
       op=(o=pimg+i*2*ps1)+2*ps1; /* pimg is complex.            */
       if(i<is0)
         {
-          df=(d=input+i*is1)+is1;
-          do {*o++=*d; *o++=0.0f;} while(++d<df);
+          ff=(f=input+i*is1)+is1;
+          do {*o++=*f; *o++=0.0f;} while(++f<ff);
         }
       do *o++=0.0f; while(o<op);
     }
@@ -769,6 +769,7 @@ frequencyconvolve(struct convolveparams *p)
 void
 convolve(struct convolveparams *p)
 {
+  gal_data_t *out;
   gal_data_t *tiles, *channels=NULL; /* `channels' has to be initialized. */
 
   /* Do the convolution. */
@@ -782,12 +783,18 @@ convolve(struct convolveparams *p)
       if(p->tilesname)
         gal_tile_block_check_tiles(tiles, p->tilesname, PROGRAM_NAME);
 
+
       /* Do the spatial convolution. */
-      gal_convolve_spatial(tiles, p->kernel);
+      out=gal_convolve_spatial(tiles, p->kernel, p->cp.numthreads,
+                               p->convoverch);
 
       /* Clean up */
+      gal_data_free(p->input);
       gal_data_array_free(tiles, gal_data_num_in_ll(tiles), 0);
       gal_data_array_free(channels, gal_data_num_in_ll(channels), 0);
+
+      /* Put the convolved dataset into the `input' to save as output. */
+      p->input=out;
     }
   else
     frequencyconvolve(p);
diff --git a/bin/convolve/ui.c b/bin/convolve/ui.c
index 8ccf534..513c3ee 100644
--- a/bin/convolve/ui.c
+++ b/bin/convolve/ui.c
@@ -339,7 +339,7 @@ ui_preparations(struct convolveparams *p)
 
   /* Read the input image as a float64 array and its WCS info. */
   p->input=gal_fits_img_read_to_type(p->filename, p->cp.hdu,
-                                     GAL_DATA_TYPE_FLOAT64, p->cp.minmapsize);
+                                     GAL_DATA_TYPE_FLOAT32, p->cp.minmapsize);
   gal_wcs_read(p->filename, p->cp.hdu, 0, 0, &p->input->nwcs, &p->input->wcs);
 
 
diff --git a/lib/convolve.c b/lib/convolve.c
index 802c50d..47cc530 100644
--- a/lib/convolve.c
+++ b/lib/convolve.c
@@ -25,16 +25,238 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <stdio.h>
 #include <errno.h>
 #include <error.h>
+#include <string.h>
 #include <stdlib.h>
 
 #include <gnuastro/tile.h>
+#include <gnuastro/threads.h>
+#include <gnuastro/multidim.h>
+#include <gnuastro/convolve.h>
+
+
+
+
+
+
+
+
+/*********************************************************************/
+/********************          Utilities          ********************/
+/*********************************************************************/
+/* See if the tile is on the edge of the hosted region or not. It doesn't
+   matter if the host is the allocated block of memory or a region in it (a
+   channel). */
+static int
+convolve_tile_is_on_edge(size_t *h, size_t *start_end_coord, size_t *k,
+                         size_t ndim)
+{
+  size_t *s=start_end_coord, *e=start_end_coord+ndim, *kf=k+ndim;
+
+  /* If the starting point of the tile is smaller than the half-kernel
+     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);
+
+  /* If control reaches here, then this is a central tile. */
+  return 0;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/*********************************************************************/
+/********************     Spatial convolution     ********************/
+/*********************************************************************/
+struct spatial_params
+{
+  gal_data_t    *out;
+  gal_data_t  *tiles;
+  gal_data_t *kernel;
+  int     convoverch;
+};
+
+
+
+
+
+/* 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)
+{
+  gal_data_t *tile=&cprm->tiles[tile_ind];
+
+  int on_edge;
+  double sum, ksum;
+  size_t i, ndim=tile->ndim, csize=tile->dsize[ndim-1];
+  gal_data_t *block=gal_tile_block(tile), *kernel=cprm->kernel;
+
+  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;
+
+  /* These variables depend on the type of the input. */
+  float *kv, *iv, *ivf, *i_start, *o_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;
+  size_t *end           = parse_coords + ndim;
+  size_t *overlap_start = parse_coords + ndim * 2;
+  size_t *o_c           = parse_coords + ndim * 3;
+
+
+  /* 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'. */
+  gal_tile_start_end_coord(tile, parse_coords, cprm->convoverch);
+  start_fastdim = pix[ndim-1];
+
+
+  /* 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. */
+  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] )
+    {
+      /* Initialize the value along the fastest dimension (it is not
+         incremented during `gal_tile_block_increment'). */
+      pix[ndim-1]=start_fastdim;
+
+      /* Go over each pixel to convolve. */
+      for(i=0;i<csize;++i)
+        {
+          /* Pointer to the pixel under consideration. */
+          in_v = i_start + i_inc + i;
+
+          /* 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
+             on this pixel. */
+          if( isnan(*in_v) )
+            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);
+              o_start=gal_tile_start_end_ind_inclusive(overlap, block,
+                                                       o_st_en);
+
+              /* 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;
+              while( o_st_en[0] + o_inc <= o_st_en[1] )
+                {
+                  /* Go over the contiguous region. */
+                  ivf = ( iv = o_start + o_inc ) + overlap->dsize[ndim-1];
+                  do
+                    {
+                      if( !isnan(*iv) )
+                        { ksum += *kv; sum += *iv * *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++, o_c);
+                }
+
+              /* Set the output value. */
+              out[ in_v - in ] = ksum==0.0f ? NAN : sum/ksum;
+            }
+
+          /* Increment the last coordinate. */
+          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);
+    }
+}
+
+
 
 
 
 /* Do spatial convolution on each mesh. */
 static void *
-gal_convolve_spatial_on_thread(void *inparm)
+convolve_spatial_on_thread(void *inparam)
 {
+  struct gal_tile_thread_param *tprm=(struct gal_tile_thread_param *)inparam;
+  struct spatial_params *cprm=(struct spatial_params *)(tprm->params);
+
+  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);
+
+  /* Initialize the necessary overlap parameters. */
+  overlap.block=block;
+  overlap.ndim=block->ndim;
+  overlap.dsize=gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, block->ndim);
+
+
+  /* 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);
+
+
+  /* Clean up, wait until all other threads finish, then return. In a
+     single thread situation, `tprm->b==NULL'. */
+  if(tprm->b) pthread_barrier_wait(tprm->b);
   return NULL;
 }
 
@@ -44,8 +266,38 @@ gal_convolve_spatial_on_thread(void *inparm)
 
 /* Convolve a dataset with a given kernel in the spatial domain. Since
    spatial convolution can be very series of tiles arranged as an array. */
-void
-gal_convolve_spatial(gal_data_t *tiles, gal_data_t *kernel)
+gal_data_t *
+gal_convolve_spatial(gal_data_t *tiles, gal_data_t *kernel,
+                     size_t numthreads, int convoverch)
 {
+  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 "
+          "input should be the same in `gal_convolve_spatial'");
+  if( block->type!=GAL_DATA_TYPE_FLOAT32
+      || kernel->type!=GAL_DATA_TYPE_FLOAT32 )
+    error(EXIT_FAILURE, 0, "`gal_convolve_spatial' currently only works on "
+          "`float32' type input and kernel");
+
+  /* Allocate the output convolved dataset. */
+  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;
+  params.tiles=tiles;
+  params.kernel=kernel;
+  params.convoverch=convoverch;
+
+  /* Do the spatial convolution on threads. */
+  gal_tile_function_on_threads(tiles, convolve_spatial_on_thread,
+                               numthreads, &params);
 
+  /* Clean up and return the output array. */
+  return out;
 }
diff --git a/lib/gnuastro/convolve.h b/lib/gnuastro/convolve.h
index 9ed035f..5751168 100644
--- a/lib/gnuastro/convolve.h
+++ b/lib/gnuastro/convolve.h
@@ -44,8 +44,9 @@ __BEGIN_C_DECLS  /* From C++ preparations */
 
 
 
-void
-gal_convolve_spatial(gal_data_t *tiles, gal_data_t *kernel);
+gal_data_t *
+gal_convolve_spatial(gal_data_t *tiles, gal_data_t *kernel,
+                     size_t numthreads, int convoverch);
 
 
 
diff --git a/lib/gnuastro/tile.h b/lib/gnuastro/tile.h
index 6891d4d..45d6b63 100644
--- a/lib/gnuastro/tile.h
+++ b/lib/gnuastro/tile.h
@@ -48,23 +48,30 @@ __BEGIN_C_DECLS  /* From C++ preparations */
 
 
 
-
 /***********************************************************************/
-/**************             About block           **********************/
+/**************              Single tile              ******************/
 /***********************************************************************/
-gal_data_t *
-gal_tile_block(gal_data_t *input);
+void
+gal_tile_start_coord(gal_data_t *tile, size_t *start_coord);
 
 void
-gal_tile_block_start_coord(gal_data_t *tile, size_t *start_coord);
+gal_tile_start_end_coord(gal_data_t *tile, size_t *start_end, int rel_block);
 
 void *
-gal_tile_block_start_end(gal_data_t *tile, gal_data_t *work,
-                         size_t *start_end);
+gal_tile_start_end_ind_inclusive(gal_data_t *tile, gal_data_t *work,
+                                 size_t *start_end_inc);
+
+
+
+/***********************************************************************/
+/**************           Allocated block         **********************/
+/***********************************************************************/
+gal_data_t *
+gal_tile_block(gal_data_t *input);
 
 size_t
 gal_tile_block_increment(gal_data_t *block, size_t *tsize,
-                         size_t num_increment);
+                         size_t num_increment, size_t *coord);
 
 void
 gal_tile_block_check_tiles(gal_data_t *tiles, char *filename,
@@ -92,6 +99,24 @@ gal_tile_all_position_two_layers(gal_data_t *input, size_t 
*channel_size,
 
 
 
+/*********************************************************************/
+/********************         On threads          ********************/
+/*********************************************************************/
+struct gal_tile_thread_param
+{
+  size_t            id; /* Id of this thread.                            */
+  void         *params; /* Input structure for higher-level settings.    */
+  size_t       *indexs; /* Indexes of actions to be done in this thread. */
+  pthread_barrier_t *b; /* Pointer the barrier for all threads.          */
+};
+
+void
+gal_tile_function_on_threads(gal_data_t *tiles, void *(*meshfunc)(void *),
+                             size_t numthreads, void *caller_params);
+
+
+
+
 __END_C_DECLS    /* From C++ preparations */
 
 #endif
diff --git a/lib/tile.c b/lib/tile.c
index f1f742d..b2cad46 100644
--- a/lib/tile.c
+++ b/lib/tile.c
@@ -31,6 +31,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/fits.h>
 #include <gnuastro/tile.h>
 #include <gnuastro/blank.h>
+#include <gnuastro/threads.h>
 #include <gnuastro/convolve.h>
 #include <gnuastro/multidim.h>
 
@@ -39,28 +40,18 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 
 
-/***********************************************************************/
-/**************        Allocated block of memory      ******************/
-/***********************************************************************/
-/* When you are working on an array, it important to know the size of the
-   allocated space in each dimension. This simple function will just follow
-   the block pointer and return the `dsize' element of lowest-level
-   structure. */
-gal_data_t *
-gal_tile_block(gal_data_t *input)
-{
-  while(input->block!=NULL) input=input->block;
-  return input;
-}
 
 
 
 
 
+/***********************************************************************/
+/**************              Single tile              ******************/
+/***********************************************************************/
 /* Calculate the starting coordinates of a tile in the allocated block of
    memory. */
 void
-gal_tile_block_start_coord(gal_data_t *tile, size_t *start_coord)
+gal_tile_start_coord(gal_data_t *tile, size_t *start_coord)
 {
   size_t ind, ndim=tile->ndim;
   gal_data_t *block=gal_tile_block(tile);
@@ -80,13 +71,55 @@ gal_tile_block_start_coord(gal_data_t *tile, size_t 
*start_coord)
 
 
 
-/* Given a tile  */
+
+/* Put the starting and ending (end point is not inclusive) coordinates of
+   a tile into the `start_end' array. It is assumed that a space of
+   `2*tile->ndim' has been already allocated (static or dynamic) before
+   this function is called.
+
+   `rel_block' (or relative-to-block) is only relevant when the tile has an
+   intermediate tile between it and the allocated space (like a channel,
+   see `gal_tile_all_position_two_layers'). If it doesn't (`tile->block'
+   points the allocated dataset), then the value to `rel_block' is
+   irrelevant.
+
+   However, when `tile->block' is its self a larger block and `rel_block'
+   is set to 0, then the starting and ending positions will be based on the
+   position within `tile->block', not the allocated space. */
+void
+gal_tile_start_end_coord(gal_data_t *tile, size_t *start_end, int rel_block)
+{
+  size_t start_ind, ndim=tile->ndim;
+  gal_data_t *block=gal_tile_block(tile);
+  gal_data_t *host=rel_block ? block : tile->block;
+
+  /* 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);
+
+  /* Get the coordinates of the starting point. */
+  gal_multidim_index_to_coord(start_ind, tile->ndim, host->dsize,
+                              start_end);
+
+  /* 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);
+}
+
+
+
+
+
+/* Put the indexs of the first/start and last/end pixels (inclusive) in a
+   tile into the `start_end' array (that has two elements). It will then
+   return the pointer to the start of the tile in the `work' data
+   structure. */
 void *
-gal_tile_block_start_end(gal_data_t *tile, gal_data_t *work,
-                         size_t *start_end)
+gal_tile_start_end_ind_inclusive(gal_data_t *tile, gal_data_t *work,
+                                 size_t *start_end_inc)
 {
-  size_t ndim=tile->ndim, *s, *sf;
   gal_data_t *block=gal_tile_block(tile);
+  size_t ndim=tile->ndim, *s, *e, *l, *sf;
   size_t *start_coord = gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, ndim);
   size_t *end_coord   = gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, ndim);
 
@@ -94,17 +127,13 @@ gal_tile_block_start_end(gal_data_t *tile, gal_data_t 
*work,
   /* The starting index can be found from the distance of the `tile->array'
      pointer and `block->array' pointer. IMPORTANT: with the type of the
      block array.  */
-  start_end[0]=gal_data_ptr_dist(block->array, tile->array, block->type);
+  start_end_inc[0]=gal_data_ptr_dist(block->array, tile->array, block->type);
 
 
   /* To find the end index, we need to know the coordinates of the starting
      point in the allocated block.  */
-  gal_multidim_index_to_coord(start_end[0], ndim, block->dsize, start_coord);
-
-
-  /* Adding the coordinates of the starting point with the tile's
-     dimensions, we will get the ending point's coordinates. */
-  gal_multidim_add_coords(start_coord, tile->dsize, end_coord, ndim);
+  gal_multidim_index_to_coord(start_end_inc[0], ndim, block->dsize,
+                              start_coord);
 
 
   /* `end_coord' is one unit ahead of the last element in the tile in every
@@ -112,16 +141,13 @@ gal_tile_block_start_end(gal_data_t *tile, gal_data_t 
*work,
      value, so we get the coordinates of the last pixel in the tile
      (inclusive). We will finally, increment that value by one to get to
      the pixel immediately outside of the tile.*/
-  sf=(s=end_coord)+ndim; do *s-=1; while(++s<sf);
+  e=end_coord;
+  l=tile->dsize;
+  sf=(s=start_coord)+ndim; do *e++ = *s + *l++ - 1; while(++s<sf);
 
 
   /* Convert the (inclusive) ending point's coordinates into an index. */
-  start_end[1]=gal_multidim_coord_to_index(ndim, block->dsize, end_coord);
-
-
-  /* Increment the (inclusive) ending point's index by one to be
-     immediately outside the tile. */
-  ++start_end[1];
+  start_end_inc[1]=gal_multidim_coord_to_index(ndim, block->dsize, end_coord);
 
 
   /* For a check:
@@ -131,8 +157,8 @@ gal_tile_block_start_end(gal_data_t *tile, gal_data_t *work,
          start_coord[2]);
   printf("end_coord: %zu, %zu, %zu\n", end_coord[0], end_coord[1],
          end_coord[2]);
-  printf("start_index: %zu\n", start_end[0]);
-  printf("end_index: %zu\n", start_end[1]);
+  printf("start_index: %zu\n", start_end_inc[0]);
+  printf("end_index: %zu\n", start_end_inc[1]);
   exit(1);
   */
 
@@ -141,7 +167,40 @@ gal_tile_block_start_end(gal_data_t *tile, gal_data_t 
*work,
      from. */
   free(end_coord);
   free(start_coord);
-  return gal_data_ptr_increment(work->array, start_end[0], work->type);
+  return gal_data_ptr_increment(work->array, start_end_inc[0], work->type);
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/***********************************************************************/
+/**************        Allocated block of memory      ******************/
+/***********************************************************************/
+/* When you are working on an array, it important to know the size of the
+   allocated space in each dimension. This simple function will just follow
+   the block pointer and return the `dsize' element of lowest-level
+   structure. */
+gal_data_t *
+gal_tile_block(gal_data_t *input)
+{
+  while(input->block!=NULL) input=input->block;
+  return input;
 }
 
 
@@ -167,27 +226,28 @@ gal_tile_block_start_end(gal_data_t *tile, gal_data_t 
*work,
    the tile has been parsed), simply incrementing by `b[n-2] * b[n-1]' will
    take us to the last row of
 
-        num_increment              increment
-        -------------              ---------
-             0               b[n-1]: fastest dimension of the block.
-             1               Similar to previous
-             .                         .
-             .                         .
-           t[n-2]            (b[n-2] * b[n-1]) - ( (t[n-2]-1) * b[n-1] )
-           t[n-2] + 1        b[n-1]
-             .                         .
-             .                         .
-          2 * t[n-2]         b[n-2] * b[n-1]
-           t[n-2]+1          b[n-1]
-             .                         .
-             .                         .
-      t[n-3] * t[n-2]        b[n-3] * b[n-2] * b[n-1]
+  num_increment      coord         increment
+  -------------      -----         ---------
+         0          (...0,0,0)     b[n-1]: fastest dimension of the block.
+         1          (...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] )
+      t[n-2] + 1    (...1,1,0)      b[n-1]
+         .              .               .
+         .              .               .
+      2 * t[n-2]    (...2,0,0)     b[n-2] * b[n-1]
+      t[n-2]+1      (...2,1,0)     b[n-1]
+         .              .                .
+         .              .                .
+   t[n-3] * t[n-2]  (..1,0,0,0)    b[n-3] * b[n-2] * b[n-1]
 
  */
 size_t
 gal_tile_block_increment(gal_data_t *block, size_t *tsize,
-                         size_t num_increment)
+                         size_t num_increment, size_t *coord)
 {
+  size_t increment;
   size_t n=block->ndim;
   size_t *b=block->dsize, *t=tsize;
 
@@ -203,20 +263,34 @@ gal_tile_block_increment(gal_data_t *block, size_t *tsize,
 
     /* 1D: the increment is just the tile size. */
     case 1:
-      return t[0];
+      increment=t[0];
+      coord[0]+=increment;
       break;
 
     /* 2D: the increment is the block's number of fastest axis pixels. */
     case 2:
-      return b[1];
+      increment=b[1];
+      ++coord[0];
       break;
 
     /* Higher dimensions. */
     default:
-      if(num_increment % t[n-2]) return b[n-1];
-      else return (b[n-2] * b[n-1]) - ( (t[n-2]-1) * b[n-1] );
+      if(num_increment % t[n-2])
+        {
+          increment=b[n-1];
+          ++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;
+        }
       break;
     }
+
+  /* Return the final increment value. */
+  return increment;
 }
 
 
@@ -237,7 +311,8 @@ gal_tile_block_check_tiles(gal_data_t *tiles, char 
*filename,
   gal_data_t *tofill, *tile;
   gal_data_t *block=gal_tile_block(tiles);
   int32_t *p, *pf, tile_index=0, *start=NULL;
-  size_t ndim=tiles->ndim, increment, start_end[2];
+  size_t ndim=tiles->ndim, increment, start_end_inc[2];
+  size_t *coord=gal_data_calloc_array(GAL_DATA_TYPE_SIZE_T, ndim);
 
   /***************************************************************/
   /*************            For a check           ****************
@@ -266,7 +341,7 @@ gal_tile_block_check_tiles(gal_data_t *tiles, char 
*filename,
     {
       /* Set the starting and ending indexs of this tile over the allocated
          block. */
-      start=gal_tile_block_start_end(tile, tofill, start_end);
+      start=gal_tile_start_end_ind_inclusive(tile, tofill, start_end_inc);
 
       /* Go over the full area of this tile. The loop will stop as soon as
          the incrementation will go over the last index of the tile. Note
@@ -274,7 +349,7 @@ gal_tile_block_check_tiles(gal_data_t *tiles, char 
*filename,
          of zero is meaningful in the calculation of the increment. */
       increment=0;
       num_increment=1;
-      while( start_end[0] + increment < start_end[1] )
+      while( start_end_inc[0] + increment <= start_end_inc[1] )
         {
           /* Parse the elements in the fastest-dimension (the contiguous
              patch of memory associated with this tile). */
@@ -284,7 +359,7 @@ gal_tile_block_check_tiles(gal_data_t *tiles, char 
*filename,
           /* Increase the increment from the start of the tile for the next
              contiguous patch. */
           increment += gal_tile_block_increment(block, tile->dsize,
-                                                num_increment++);
+                                                num_increment++, coord);
         }
 
       /* Increment the index for the next tile. */
@@ -295,6 +370,7 @@ gal_tile_block_check_tiles(gal_data_t *tiles, char 
*filename,
   gal_fits_img_write(tofill, filename, NULL, program_name);
 
   /* Clean up. */
+  free(coord);
   gal_data_free(tofill);
 }
 
@@ -559,7 +635,7 @@ gal_tile_all_position(gal_data_t *input, size_t *regular,
   if(input->block)
     {
       start=gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, input->ndim);
-      gal_tile_block_start_coord(input, start);
+      gal_tile_start_coord(input, start);
     }
 
 
@@ -699,3 +775,98 @@ gal_tile_all_position_two_layers(gal_data_t *input, size_t 
*channel_size,
       gal_tile_all_position(&ch[i], tile_size, remainderfrac, &t, 1);
     }
 }
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/*********************************************************************/
+/********************         On threads          ********************/
+/*********************************************************************/
+/* Run a given function on the given tiles. */
+void
+gal_tile_function_on_threads(gal_data_t *tiles, void *(*function)(void *),
+                             size_t numthreads, void *caller_params)
+{
+  int err;
+  pthread_t t;          /* All thread ids saved in this, not used. */
+  pthread_attr_t attr;
+  pthread_barrier_t b;
+  struct gal_tile_thread_param *prm;
+  size_t i, *indexs, thrdcols, numbarriers;
+  size_t numtiles=gal_data_num_in_ll(tiles);
+
+  /* Allocate the array of parameters structure structures. */
+  prm=malloc(numthreads*sizeof *prm);
+  if(prm==NULL)
+    {
+      fprintf(stderr, "%zu bytes could not be allocated for prm.",
+              numthreads*sizeof *prm);
+      exit(EXIT_FAILURE);
+    }
+
+  /* Distribute the actions into the threads: */
+  gal_threads_dist_in_threads(numtiles, numthreads, &indexs, &thrdcols);
+
+  /* Do the job: when only one thread is necessary, there is no need to
+     spin off one thread, just call the function directly (spinning off
+     threads is expensive). This is for the generic thread spinner
+     function, not this simple function where `numthreads' is a
+     constant. */
+  if(numthreads==1)
+    {
+      prm[0].id=0;
+      prm[0].b=NULL;
+      prm[0].indexs=indexs;
+      prm[0].params=caller_params;
+      function(&prm[0]);
+    }
+  else
+    {
+      /* Initialize the attributes. Note that this running thread
+         (that spinns off the nt threads) is also a thread, so the
+         number the barriers should be one more than the number of
+         threads spinned off. */
+      numbarriers = (numtiles<numthreads ? numtiles : numthreads) + 1;
+      gal_threads_attr_barrier_init(&attr, &b, numbarriers);
+
+      /* Spin off the threads: */
+      for(i=0;i<numthreads;++i)
+        if(indexs[i*thrdcols]!=GAL_THREADS_NON_THRD_INDEX)
+          {
+            prm[i].id=i;
+            prm[i].b=&b;
+            prm[i].params=caller_params;
+            prm[i].indexs=&indexs[i*thrdcols];
+            err=pthread_create(&t, &attr, function, &prm[i]);
+            if(err)
+              {
+                fprintf(stderr, "can't create thread %zu", i);
+                exit(EXIT_FAILURE);
+              }
+          }
+
+      /* Wait for all threads to finish and free the spaces. */
+      pthread_barrier_wait(&b);
+      pthread_attr_destroy(&attr);
+      pthread_barrier_destroy(&b);
+    }
+
+  /* Clean up. */
+  free(indexs);
+}