[gnuastro-commits] master 4839a11 104/125: Preparations for tile interpolation

[Mohammad Akhlaghi]

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

    Preparations for tile interpolation
    
    The main change in this commit was in the Statistics program to prepare it
    for reading and interpolating single measured values on a tile. Everything
    is place for the job to be started, but many small issues came up that had
    to be resolved so they are being commited here. They are listed here:
    
     - In convolve (when we want to make a kernel), we need to divide both
       input images by their sum. However, since the sum is returned in
       `float64' type, the output would be stored in `float64'. However,
       afterwards, we have assumed that the raw input is `float32'. This made
       bad outputs. With this commit, the sum value is converted to `float32',
       so the division gives a `float32' type output.
    
     - While finding the cause of the above problem in Convolve, some function
       names were also re-written with underscores (`_') to be more clear and
       when allocation was necessary, the old direct calls to `malloc' were
       replaced with `gal_data_malloc_array'. Also, after convolution in the
       frequency domain, we would store the real (not complex) array in
       `pimg'. But this would cause confusion, because until then `pimg' is
       complex. So now we have a new `rpad' array to keep the real-but-padded
       array for explicity.
    
     - In the Statistics programs, all the names with `row' (like the
       `ui_add_to_print_in_row' function) have been changed to `single_value'
       (so that function is now `ui_add_to_single_value'). Since the same class
       of operations can be done on tiles now.
    
     - In Statistics, to print the tile numbers, we need to keep the format of
       the input. This is mainly a trouble in 1D arrays that can be both an
       image and a table column.
    
     - `statistics_on_tile' has been defined and now does all the preparations
       for the actual operation on each tile (from the next commit).
    
     - In the tiles library, `gal_tile_block_check_tiles' now returns the final
       data structure, it doesn't write it to a file. This is more clear and
       allows the caller to decide what to do with the labeled image.
    
     - The returned values of `gal_tile_all_position' and
       `gal_tile_all_position_two_layers' are now the array of how many tiles
       there are in each dimension. Previously they just returned the total
       number of tiles that was only useful for allocation but would loose the
       dimensionality. When we want to store/show the values on each tile (with
       one tile taking one pixel), we need the dimensionality for easy writing
       to a FITS file.
---
 bin/convolve/convolve.c     |  94 +++++++++++++++++++-------------------
 bin/convolve/main.h         |   2 +-
 bin/convolve/ui.c           |   9 ++--
 bin/fits/fits.c             |   4 +-
 bin/statistics/args.h       |  95 ++++++++++++++++++++++-----------------
 bin/statistics/main.h       |  17 ++++++-
 bin/statistics/statistics.c | 107 ++++++++++++++++++++++++++++++++++++++++++--
 bin/statistics/ui.c         |  93 +++++++++++++++++++++++++++++---------
 bin/statistics/ui.h         |   5 ++-
 doc/gnuastro.texi           |  29 ++++++------
 lib/commonopts.h            |   6 +--
 lib/gnuastro/tile.h         |   9 ++--
 lib/options.h               |   1 +
 lib/tile.c                  |  74 +++++++++++++++---------------
 14 files changed, 362 insertions(+), 183 deletions(-)

diff --git a/bin/convolve/convolve.c b/bin/convolve/convolve.c
index 83857de..8195342 100644
--- a/bin/convolve/convolve.c
+++ b/bin/convolve/convolve.c
@@ -58,12 +58,10 @@ complextoreal(double *c, size_t size, int action, double 
**output)
 {
   double *out, *o, *of;
 
-  errno=0;
-  *output=out=malloc(size*sizeof *out);
-  if(out==NULL)
-    error(EXIT_FAILURE, errno, "%zu bytes for out in complextoreal",
-          size*sizeof *out);
+  /* Allocate the space for the real array. */
+  *output=out=gal_data_malloc_array(GAL_DATA_TYPE_FLOAT64, size);
 
+  /* Fill the real array with the derived value from the complex array. */
   of=(o=out)+size;
   switch(action)
     {
@@ -187,7 +185,7 @@ complexarraydivide(double *a, double *b, size_t size, 
double minsharpspec)
 /*************      Padding and initializing      *****************/
 /******************************************************************/
 void
-makepaddedcomplex(struct convolveparams *p)
+frequency_make_padded_complex(struct convolveparams *p)
 {
   size_t i, ps0, ps1;
   double *o, *op, *pimg, *pker;
@@ -210,11 +208,7 @@ makepaddedcomplex(struct convolveparams *p)
 
 
   /* Allocate the space for the padded input image and fill it. */
-  errno=0;
-  pimg=p->pimg=malloc(2*ps0*ps1*sizeof *pimg);
-  if(pimg==NULL)
-    error(EXIT_FAILURE, errno, "%zu bytes for pimg",
-          2*ps0*ps1*sizeof *pimg);
+  pimg=p->pimg=gal_data_malloc_array(GAL_DATA_TYPE_FLOAT64, 2*ps0*ps1);
   for(i=0;i<ps0;++i)
     {
       op=(o=pimg+i*2*ps1)+2*ps1; /* pimg is complex.            */
@@ -228,11 +222,7 @@ makepaddedcomplex(struct convolveparams *p)
 
 
   /* Allocate the space for the padded Kernel and fill it. */
-  errno=0;
-  pker=p->pker=malloc(2*ps0*ps1*sizeof *pker);
-  if(pker==NULL)
-    error(EXIT_FAILURE, errno, "%zu bytes for pker",
-          2*ps0*ps1*sizeof *pker);
+  pker=p->pker=gal_data_malloc_array(GAL_DATA_TYPE_FLOAT64, 2*ps0*ps1);
   for(i=0;i<ps0;++i)
     {
       op=(o=pker+i*2*ps1)+2*ps1; /* pker is complex.            */
@@ -259,22 +249,22 @@ void
 removepaddingcorrectroundoff(struct convolveparams *p)
 {
   size_t ps1=p->ps1;
-  size_t i, hi0, hi1;
   size_t *isize=p->input->dsize;
   float *o, *input=p->input->array;
-  double *d, *df, *start, *pimg=p->pimg;
-
-  /* Set all the necessary parameters to crop the desired region. hi0
-     and hi1 are the coordinates of the first pixel in the output
-     image. In the case of deconvolution, if the maximum radius is
-     larger smaller than the input image, we will also only be using
-     region that contains non-zero rows and columns.*/
+  double *d, *df, *start, *rpad=p->rpad;
+  size_t i, hi0, hi1, mkwidth=2*p->makekernel-1;
+
+  /* Set all the necessary parameters to crop the desired region. hi0 and
+     hi1 are the coordinates of the first pixel in the output image. In the
+     case of deconvolution, if the maximum radius is larger than the input
+     image, we will also only be using region that contains non-zero rows
+     and columns.*/
   if(p->makekernel)
     {
-      hi0      = 2*p->makekernel-1 < isize[0] ? p->ps0/2-p->makekernel : 0;
-      hi1      = 2*p->makekernel-1 < isize[1] ? p->ps1/2-p->makekernel : 0;
-      isize[0] = 2*p->makekernel-1 < isize[0] ? 2*p->makekernel-1 : isize[0];
-      isize[1] = 2*p->makekernel-1 < isize[1] ? 2*p->makekernel-1 : isize[1];
+      hi0      = mkwidth < isize[0] ? p->ps0/2-p->makekernel : 0;
+      hi1      = mkwidth < isize[1] ? p->ps1/2-p->makekernel : 0;
+      isize[0] = mkwidth < isize[0] ? 2*p->makekernel-1 : isize[0];
+      isize[1] = mkwidth < isize[1] ? 2*p->makekernel-1 : isize[1];
     }
   else
     {
@@ -284,10 +274,11 @@ removepaddingcorrectroundoff(struct convolveparams *p)
 
   /* To start with, `start' points to the first pixel in the final
      image: */
-  start=&pimg[hi0*ps1+hi1];
+  start=&rpad[hi0*ps1+hi1];
   for(i=0;i<isize[0];++i)
     {
       o = &input[ i * isize[1] ];
+
       df = ( d = start + i * ps1 ) + isize[1];
       do
         *o++ = ( *d<-CONVFLOATINGPOINTERR || *d>CONVFLOATINGPOINTERR )
@@ -630,7 +621,7 @@ twodimensionfft(struct convolveparams *p, struct 
fftonthreadparams *fp,
 
 
 void
-frequencyconvolve(struct convolveparams *p)
+convolve_frequency(struct convolveparams *p)
 {
   double *tmp;
   size_t dsize[2];
@@ -641,7 +632,7 @@ frequencyconvolve(struct convolveparams *p)
 
   /* Make the padded arrays. */
   if(!p->cp.quiet) gettimeofday(&t1, NULL);
-  makepaddedcomplex(p);
+  frequency_make_padded_complex(p);
   if(!p->cp.quiet)
     gal_timing_report(&t1, "Input and Kernel images padded.", 1);
   if(p->checkfreqsteps)
@@ -715,16 +706,16 @@ frequencyconvolve(struct convolveparams *p)
   if(!p->cp.quiet) gettimeofday(&t1, NULL);
   twodimensionfft(p, fp, -1);
   if(p->makekernel)
-    correctdeconvolve(p, &tmp);
+    correctdeconvolve(p, &p->rpad);
   else
-    complextoreal(p->pimg, p->ps0*p->ps1, COMPLEX_TO_REAL_REAL, &tmp);
+    complextoreal(p->pimg, p->ps0*p->ps1, COMPLEX_TO_REAL_REAL, &p->rpad);
   if(!p->cp.quiet)
     gal_timing_report(&t1, "Converted back to the spatial domain.", 1);
   if(p->checkfreqsteps)
     {
-      data->array=tmp; data->name="padded output";
+      data->array=p->rpad; data->name="padded output";
       gal_fits_img_write(data, p->freqstepsname, NULL, PROGRAM_STRING);
-      data->name=NULL;
+      data->name=NULL; data->array=NULL;
     }
 
   /* Free the padded arrays (they are no longer needed) and put the
@@ -732,7 +723,6 @@ frequencyconvolve(struct convolveparams *p)
   gal_data_free(data);
   free(p->pimg);
   free(p->pker);
-  p->pimg=tmp;
 
   /* Crop out the center, numbers smaller than 10^{-17} are errors,
      remove them. */
@@ -769,29 +759,37 @@ frequencyconvolve(struct convolveparams *p)
 void
 convolve(struct convolveparams *p)
 {
-  gal_data_t *out;
-  gal_data_t *tiles, *channels=NULL; /* `channels' has to be initialized. */
+  size_t *tsize;
+  gal_data_t *out, *check;
+  gal_data_t *tiles, *channels;
+  struct gal_options_common_params *cp=&p->cp;
+
 
   /* Do the convolution. */
   if(p->domain==CONVOLVE_DOMAIN_SPATIAL)
     {
       /* Prepare the mesh structure. */
-      gal_tile_all_position_two_layers(p->input, p->cp.channelsize,
-                                       p->cp.tilesize, p->cp.remainderfrac,
-                                       &channels, &tiles);
+      tsize=gal_tile_all_position_two_layers(p->input, cp->channelsize,
+                                             cp->tilesize, cp->remainderfrac,
+                                             &channels, &tiles);
 
       /* Save the tile IDs if they are requested. */
-      if(p->tilesname)
-        gal_tile_block_check_tiles(tiles, p->tilesname, PROGRAM_NAME);
+      if(cp->tilecheckname)
+        {
+          check=gal_tile_block_check_tiles(tiles);
+          gal_fits_img_write(check, cp->tilecheckname, NULL, PROGRAM_NAME);
+          gal_data_free(check);
+        }
 
       /* Do the spatial convolution. One of the main reason someone would
          want to do spatial domain convolution with this Convolve program
          is edge correction. So by default we assume it and will only
          ignore it if the user asks.*/
-      out=gal_convolve_spatial(tiles, p->kernel, p->cp.numthreads,
-                               !p->noedgecorrection, p->cp.workoverch);
+      out=gal_convolve_spatial(tiles, p->kernel, cp->numthreads,
+                               !p->noedgecorrection, cp->workoverch);
 
       /* Clean up */
+      free(tsize);
       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);
@@ -800,9 +798,9 @@ convolve(struct convolveparams *p)
       p->input=out;
     }
   else
-    frequencyconvolve(p);
+    convolve_frequency(p);
 
   /* Save the output (which is in p->input) array. */
-  gal_fits_img_write_to_type(p->input, p->cp.output, NULL, PROGRAM_STRING,
-                             p->cp.type);
+  gal_fits_img_write_to_type(p->input, cp->output, NULL, PROGRAM_STRING,
+                             cp->type);
 }
diff --git a/bin/convolve/main.h b/bin/convolve/main.h
index bcfd67d..03edab5 100644
--- a/bin/convolve/main.h
+++ b/bin/convolve/main.h
@@ -89,10 +89,10 @@ struct convolveparams
   gal_data_t         *kernel;  /* Input Kernel array.                     */
   double               *pimg;  /* Padded image array.                     */
   double               *pker;  /* Padded kernel array.                    */
+  double               *rpad;  /* Real final image before removing pad'd. */
   size_t                 ps0;  /* Padded size along first C axis.         */
   size_t                 ps1;  /* Padded size along second C axis.        */
   char        *freqstepsname;  /* Name of file to check frequency steps.  */
-  char            *tilesname;  /* Name of file to check tiles.            */
   time_t             rawtime;  /* Starting time of the program.           */
 };
 
diff --git a/bin/convolve/ui.c b/bin/convolve/ui.c
index 483ee4a..2963ac6 100644
--- a/bin/convolve/ui.c
+++ b/bin/convolve/ui.c
@@ -333,9 +333,10 @@ ui_preparations(struct convolveparams *p)
     }
   if(cp->checktiles)
     {
-      p->tilesname=gal_checkset_automatic_output(&p->cp, p->filename,
-                                                 "_tiled.fits");
-      gal_checkset_check_remove_file(p->tilesname, 0, p->cp.dontdelete);
+      p->cp.tilecheckname=gal_checkset_automatic_output(&p->cp, p->filename,
+                                                        "_tiled.fits");
+      gal_checkset_check_remove_file(p->cp.tilecheckname, 0,
+                                     p->cp.dontdelete);
     }
 
 
@@ -391,9 +392,11 @@ ui_preparations(struct convolveparams *p)
          1 and their division (in the frequency domain) would be
          meaningful. */
       sum=gal_statistics_sum(p->input);
+      sum=gal_data_copy_to_new_type_free(sum, GAL_DATA_TYPE_FLOAT32);
       p->input = gal_arithmetic(GAL_ARITHMETIC_OP_DIVIDE,
                                 GAL_ARITHMETIC_FLAGS_ALL, p->input, sum);
       sum=gal_statistics_sum(p->kernel);
+      sum=gal_data_copy_to_new_type_free(sum, GAL_DATA_TYPE_FLOAT32);
       p->kernel = gal_arithmetic(GAL_ARITHMETIC_OP_DIVIDE,
                                 GAL_ARITHMETIC_FLAGS_ALL, p->kernel, sum);
     }
diff --git a/bin/fits/fits.c b/bin/fits/fits.c
index 6690246..9f141c9 100644
--- a/bin/fits/fits.c
+++ b/bin/fits/fits.c
@@ -224,8 +224,8 @@ fits_print_extension_info(struct fitsparams *p)
     {
       printf("%s\nRun on %s-----\n", PROGRAM_STRING, ctime(&p->rawtime));
       printf("HDU (extension) information: `%s'.\n", p->filename);
-      printf(" Column 1: Index (counting from 0).\n");
-      printf(" Column 2: Name (`EXTNAME' in FITS standard).\n");
+      printf(" Column 1: Index (counting from 0, usable with `--hdu').\n");
+      printf(" Column 2: Name (`EXTNAME' in FITS standard, usable with 
`--hdu').\n");
       printf(" Column 3: Image data type or `table' format (ASCII or "
              "binary).\n");
       printf(" Column 4: Size of data in HDU.\n");
diff --git a/bin/statistics/args.h b/bin/statistics/args.h
index 44091db..ebb3d37 100644
--- a/bin/statistics/args.h
+++ b/bin/statistics/args.h
@@ -104,8 +104,8 @@ struct argp_option program_options[] =
 
     {
       0, 0, 0, 0,
-      "Single value measurements (to print in one row)",
-      ARGS_GROUP_IN_ONE_ROW
+      "Single value measurements",
+      ARGS_GROUP_SINGLE_VALUE
     },
     {
       "number",
@@ -113,13 +113,13 @@ struct argp_option program_options[] =
       0,
       0,
       "Number (non-blank).",
-      ARGS_GROUP_IN_ONE_ROW,
-      &p->toprint,
+      ARGS_GROUP_SINGLE_VALUE,
+      &p->singlevalue,
       GAL_OPTIONS_NO_ARG_TYPE,
       GAL_OPTIONS_RANGE_ANY,
       GAL_OPTIONS_NOT_MANDATORY,
       GAL_OPTIONS_NOT_SET,
-      ui_add_to_print_in_row
+      ui_add_to_single_value
     },
     {
       "minimum",
@@ -127,13 +127,13 @@ struct argp_option program_options[] =
       0,
       0,
       "Minimum.",
-      ARGS_GROUP_IN_ONE_ROW,
-      &p->toprint,
+      ARGS_GROUP_SINGLE_VALUE,
+      &p->singlevalue,
       GAL_OPTIONS_NO_ARG_TYPE,
       GAL_OPTIONS_RANGE_0_OR_1,
       GAL_OPTIONS_NOT_MANDATORY,
       GAL_OPTIONS_NOT_SET,
-      ui_add_to_print_in_row
+      ui_add_to_single_value
     },
     {
       "maximum",
@@ -141,13 +141,13 @@ struct argp_option program_options[] =
       0,
       0,
       "Maximum.",
-      ARGS_GROUP_IN_ONE_ROW,
-      &p->toprint,
+      ARGS_GROUP_SINGLE_VALUE,
+      &p->singlevalue,
       GAL_OPTIONS_NO_ARG_TYPE,
       GAL_OPTIONS_RANGE_0_OR_1,
       GAL_OPTIONS_NOT_MANDATORY,
       GAL_OPTIONS_NOT_SET,
-      ui_add_to_print_in_row
+      ui_add_to_single_value
     },
     {
       "sum",
@@ -155,13 +155,13 @@ struct argp_option program_options[] =
       0,
       0,
       "Sum.",
-      ARGS_GROUP_IN_ONE_ROW,
-      &p->toprint,
+      ARGS_GROUP_SINGLE_VALUE,
+      &p->singlevalue,
       GAL_OPTIONS_NO_ARG_TYPE,
       GAL_OPTIONS_RANGE_0_OR_1,
       GAL_OPTIONS_NOT_MANDATORY,
       GAL_OPTIONS_NOT_SET,
-      ui_add_to_print_in_row
+      ui_add_to_single_value
     },
     {
       "mean",
@@ -169,13 +169,13 @@ struct argp_option program_options[] =
       0,
       0,
       "Mean.",
-      ARGS_GROUP_IN_ONE_ROW,
-      &p->toprint,
+      ARGS_GROUP_SINGLE_VALUE,
+      &p->singlevalue,
       GAL_OPTIONS_NO_ARG_TYPE,
       GAL_OPTIONS_RANGE_0_OR_1,
       GAL_OPTIONS_NOT_MANDATORY,
       GAL_OPTIONS_NOT_SET,
-      ui_add_to_print_in_row
+      ui_add_to_single_value
     },
     {
       "std",
@@ -183,13 +183,13 @@ struct argp_option program_options[] =
       0,
       0,
       "Standad deviation.",
-      ARGS_GROUP_IN_ONE_ROW,
-      &p->toprint,
+      ARGS_GROUP_SINGLE_VALUE,
+      &p->singlevalue,
       GAL_OPTIONS_NO_ARG_TYPE,
       GAL_OPTIONS_RANGE_0_OR_1,
       GAL_OPTIONS_NOT_MANDATORY,
       GAL_OPTIONS_NOT_SET,
-      ui_add_to_print_in_row
+      ui_add_to_single_value
     },
     {
       "median",
@@ -197,13 +197,13 @@ struct argp_option program_options[] =
       0,
       0,
       "Median.",
-      ARGS_GROUP_IN_ONE_ROW,
-      &p->toprint,
+      ARGS_GROUP_SINGLE_VALUE,
+      &p->singlevalue,
       GAL_OPTIONS_NO_ARG_TYPE,
       GAL_OPTIONS_RANGE_0_OR_1,
       GAL_OPTIONS_NOT_MANDATORY,
       GAL_OPTIONS_NOT_SET,
-      ui_add_to_print_in_row
+      ui_add_to_single_value
     },
     {
       "quantile",
@@ -211,13 +211,13 @@ struct argp_option program_options[] =
       "FLT[,...]",
       0,
       "Quantile (multiple values acceptable).",
-      ARGS_GROUP_IN_ONE_ROW,
-      &p->toprint,
+      ARGS_GROUP_SINGLE_VALUE,
+      &p->singlevalue,
       GAL_DATA_TYPE_FLOAT32,
       GAL_OPTIONS_RANGE_GE_0_LE_1,
       GAL_OPTIONS_NOT_MANDATORY,
       GAL_OPTIONS_NOT_SET,
-      ui_add_to_print_in_row
+      ui_add_to_single_value
     },
     {
       "quantfunc",
@@ -225,13 +225,13 @@ struct argp_option program_options[] =
       "FLT[,...]",
       0,
       "Quantile function (multiple values acceptable).",
-      ARGS_GROUP_IN_ONE_ROW,
-      &p->toprint,
+      ARGS_GROUP_SINGLE_VALUE,
+      &p->singlevalue,
       GAL_DATA_TYPE_FLOAT32,
       GAL_OPTIONS_RANGE_GE_0_LE_1,
       GAL_OPTIONS_NOT_MANDATORY,
       GAL_OPTIONS_NOT_SET,
-      ui_add_to_print_in_row
+      ui_add_to_single_value
     },
     {
       "mode",
@@ -239,13 +239,13 @@ struct argp_option program_options[] =
       0,
       0,
       "Mode (Appendix C of arXiv:1505.011664).",
-      ARGS_GROUP_IN_ONE_ROW,
-      &p->toprint,
+      ARGS_GROUP_SINGLE_VALUE,
+      &p->singlevalue,
       GAL_OPTIONS_NO_ARG_TYPE,
       GAL_OPTIONS_RANGE_0_OR_1,
       GAL_OPTIONS_NOT_MANDATORY,
       GAL_OPTIONS_NOT_SET,
-      ui_add_to_print_in_row
+      ui_add_to_single_value
     },
     {
       "modequant",
@@ -253,13 +253,13 @@ struct argp_option program_options[] =
       0,
       0,
       "Mode quantile (see --mode)",
-      ARGS_GROUP_IN_ONE_ROW,
-      &p->toprint,
+      ARGS_GROUP_SINGLE_VALUE,
+      &p->singlevalue,
       GAL_OPTIONS_NO_ARG_TYPE,
       GAL_OPTIONS_RANGE_0_OR_1,
       GAL_OPTIONS_NOT_MANDATORY,
       GAL_OPTIONS_NOT_SET,
-      ui_add_to_print_in_row
+      ui_add_to_single_value
     },
     {
       "modesym",
@@ -267,13 +267,13 @@ struct argp_option program_options[] =
       0,
       0,
       "Mode symmetricity (see --mode).",
-      ARGS_GROUP_IN_ONE_ROW,
-      &p->toprint,
+      ARGS_GROUP_SINGLE_VALUE,
+      &p->singlevalue,
       GAL_OPTIONS_NO_ARG_TYPE,
       GAL_OPTIONS_RANGE_0_OR_1,
       GAL_OPTIONS_NOT_MANDATORY,
       GAL_OPTIONS_NOT_SET,
-      ui_add_to_print_in_row
+      ui_add_to_single_value
     },
     {
       "modesymvalue",
@@ -281,13 +281,26 @@ struct argp_option program_options[] =
       0,
       0,
       "Value at mode symmetricity (see --mode).",
-      ARGS_GROUP_IN_ONE_ROW,
-      &p->toprint,
+      ARGS_GROUP_SINGLE_VALUE,
+      &p->singlevalue,
       GAL_OPTIONS_NO_ARG_TYPE,
       GAL_OPTIONS_RANGE_0_OR_1,
       GAL_OPTIONS_NOT_MANDATORY,
       GAL_OPTIONS_NOT_SET,
-      ui_add_to_print_in_row
+      ui_add_to_single_value
+    },
+    {
+      "ontile",
+      ARGS_OPTION_KEY_ONTILE,
+      0,
+      0,
+      "Do the measurements on separate tiles, not all.",
+      ARGS_GROUP_SINGLE_VALUE,
+      &p->ontile,
+      GAL_OPTIONS_NO_ARG_TYPE,
+      GAL_OPTIONS_RANGE_0_OR_1,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
     },
 
 
diff --git a/bin/statistics/main.h b/bin/statistics/main.h
index 59f6db1..8c74ce5 100644
--- a/bin/statistics/main.h
+++ b/bin/statistics/main.h
@@ -37,6 +37,17 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 
 
+/* Input formats. */
+enum statistics_input_format
+  {
+    INPUT_FORMAT_INVALID,
+
+    INPUT_FORMAT_TABLE,
+    INPUT_FORMAT_IMAGE,
+  };
+
+
+
 
 
 /* Main program parameters structure */
@@ -44,7 +55,7 @@ struct statisticsparams
 {
   /* From command-line */
   struct gal_options_common_params cp; /* Common parameters.             */
-  struct gal_linkedlist_ill  *toprint; /* Options to print in one row.   */
+  struct gal_linkedlist_ill  *singlevalue; /* Single value calculations. */
   struct gal_linkedlist_dll  *tp_args; /* Arguments for printing.        */
   char          *inputname;  /* Input filename.                          */
   char             *column;  /* Column name or number if input is table. */
@@ -53,6 +64,7 @@ struct statisticsparams
   float           lessthan;  /* Only use values <  this value.           */
   float           quantmin;  /* Quantile min or range: from Q to 1-Q.    */
   float           quantmax;  /* Quantile maximum.                        */
+  uint8_t           ontile;  /* Do single value calculations on tiles.   */
 
   uint8_t        asciihist;  /* Print an ASCII histogram.                */
   uint8_t         asciicfp;  /* Print an ASCII cumulative frequency plot.*/
@@ -70,11 +82,12 @@ struct statisticsparams
   uint8_t        maxbinone;  /* Set the maximum bin to 1.                */
 
   /* Internal */
+  uint8_t      inputformat;  /* Format of input dataset.                 */
   int          numoutfiles;  /* Number of output files made in this run. */
   uint8_t        needssort;  /* If sorting is needed.                    */
   gal_data_t        *input;  /* Input data structure.                    */
-  gal_data_t    *reference;  /* Reference data structure.                */
   gal_data_t       *sorted;  /* Sorted input data structure.             */
+  gal_data_t    *reference;  /* Reference data structure.                */
   int               isfits;  /* Input is a FITS file.                    */
   int             hdu_type;  /* Type of HDU (image or table).            */
   time_t           rawtime;  /* Starting time of the program.            */
diff --git a/bin/statistics/statistics.c b/bin/statistics/statistics.c
index f6deda4..f7d692d 100644
--- a/bin/statistics/statistics.c
+++ b/bin/statistics/statistics.c
@@ -31,6 +31,8 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <stdlib.h>
 
 #include <gnuastro/fits.h>
+#include <gnuastro/tile.h>
+#include <gnuastro/multidim.h>
 #include <gnuastro/arithmetic.h>
 #include <gnuastro/statistics.h>
 
@@ -96,7 +98,7 @@ statistics_print_one_row(struct statisticsparams *p)
   /* The user can ask for any of the operators more than once, also some
      operators might return more than one usable value (like mode). So we
      will calculate the desired values once, and then print them. */
-  for(tmp=p->toprint; tmp!=NULL; tmp=tmp->next)
+  for(tmp=p->singlevalue; tmp!=NULL; tmp=tmp->next)
     switch(tmp->v)
       {
       /* Calculate respective values. Checking with `if(num==NULL)' gives
@@ -140,7 +142,7 @@ statistics_print_one_row(struct statisticsparams *p)
 
 
   /* Print every requested number. */
-  for(tmp=p->toprint; tmp!=NULL; tmp=tmp->next)
+  for(tmp=p->singlevalue; tmp!=NULL; tmp=tmp->next)
     {
       /* By default don't free anything. */
       mustfree=0;
@@ -227,6 +229,102 @@ statistics_print_one_row(struct statisticsparams *p)
 
 
 /*******************************************************************/
+/**************         Single value on tile         ***************/
+/*******************************************************************/
+void
+statistics_on_tile(struct statisticsparams *p)
+{
+  uint8_t type;
+  size_t i, numtiles=0, ntiles, *tsize;
+  struct gal_linkedlist_ill *operation;
+  gal_data_t *channels, *tiles, *check, *tile, *values;
+
+  /* Make the tile structure. */
+  tsize=gal_tile_all_position_two_layers(p->input, p->cp.channelsize,
+                                         p->cp.tilesize, p->cp.remainderfrac,
+                                         &channels, &tiles);
+
+  /* Number of tiles. */
+  ntiles=gal_multidim_total_size(p->input->ndim, tsize);
+
+  /* Save the tile IDs if they are requested. */
+  if(p->cp.tilecheckname)
+    {
+      check=gal_tile_block_check_tiles(tiles);
+      if(p->inputformat==INPUT_FORMAT_IMAGE)
+        gal_fits_img_write(check, p->cp.tilecheckname, NULL, PROGRAM_NAME);
+      else
+        gal_table_write(check, NULL, p->cp.tableformat, p->cp.tilecheckname,
+                        p->cp.dontdelete);
+      gal_data_free(check);
+    }
+
+  /* Do the operation on each tile. */
+  for(operation=p->singlevalue; operation!=NULL; operation=operation->next)
+    {
+      /* Set the type of the output array. */
+      switch(operation->v)
+        {
+        case ARGS_OPTION_KEY_NUMBER:
+          type=GAL_DATA_TYPE_INT32; break;
+
+        case ARGS_OPTION_KEY_MINIMUM:
+        case ARGS_OPTION_KEY_MAXIMUM:
+        case ARGS_OPTION_KEY_SUM:
+        case ARGS_OPTION_KEY_MEDIAN:
+        case ARGS_OPTION_KEY_MODE:
+        case ARGS_OPTION_KEY_QUANTFUNC:
+          type=p->input->type; break;
+
+        case ARGS_OPTION_KEY_MEAN:
+        case ARGS_OPTION_KEY_STD:
+        case ARGS_OPTION_KEY_MODEQUANT:
+        case ARGS_OPTION_KEY_MODESYM:
+        case ARGS_OPTION_KEY_MODESYMVALUE:
+        case ARGS_OPTION_KEY_QUANTILE:
+          type=GAL_DATA_TYPE_FLOAT64; break;
+
+        default:
+          error(EXIT_FAILURE, 0, "a bug! %d is not a recognized operation "
+                "code in `statistics_on_tile'", operation->v);
+        }
+
+      /* Allocate the space necessary to keep the value for each tile. */
+      values=gal_data_alloc(NULL, type, p->input->ndim, tsize, NULL, 0,
+                            p->input->minmapsize, NULL, NULL, NULL);
+
+      for(tile=tiles; tile!=NULL; tile=tile->next);
+    }
+
+
+  /* Clean up. */
+  free(tsize);
+  gal_data_array_free(tiles, ntiles, 0);
+  gal_data_array_free(channels, gal_data_num_in_ll(channels), 0);
+
+  printf("\n... In statistics_on_tile ...\n");
+  exit(0);
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/*******************************************************************/
 /**************             ASCII plots              ***************/
 /*******************************************************************/
 static void
@@ -728,10 +826,11 @@ statistics(struct statisticsparams *p)
   int print_basic_info=1;
 
   /* Print the one-row numbers if the user asked for them. */
-  if(p->toprint)
+  if(p->singlevalue)
     {
       print_basic_info=0;
-      statistics_print_one_row(p);
+      if(p->ontile) statistics_on_tile(p);
+      else          statistics_print_one_row(p);
     }
 
   /* Print the ASCII plots if requested. */
diff --git a/bin/statistics/ui.c b/bin/statistics/ui.c
index 1abe260..2cc0de5 100644
--- a/bin/statistics/ui.c
+++ b/bin/statistics/ui.c
@@ -28,6 +28,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <stdio.h>
 
 #include <gnuastro/fits.h>
+#include <gnuastro/tile.h>
 #include <gnuastro/qsort.h>
 #include <gnuastro/blank.h>
 #include <gnuastro/table.h>
@@ -85,7 +86,7 @@ doc[] = GAL_STRINGS_TOP_HELP_INFO PROGRAM_NAME" will do 
statistical "
 /* Option groups particular to this program. */
 enum program_args_groups
 {
-  ARGS_GROUP_IN_ONE_ROW = GAL_OPTIONS_GROUP_AFTER_COMMON,
+  ARGS_GROUP_SINGLE_VALUE = GAL_OPTIONS_GROUP_AFTER_COMMON,
   ARGS_GROUP_PARTICULAR_STAT,
   ARGS_GROUP_HIST_CFP,
 };
@@ -203,8 +204,8 @@ parse_opt(int key, char *arg, struct argp_state *state)
 
 
 static void *
-ui_add_to_print_in_row(struct argp_option *option, char *arg,
-                       char *filename, size_t lineno, void *params)
+ui_add_to_single_value(struct argp_option *option, char *arg,
+                      char *filename, size_t lineno, void *params)
 {
   size_t i;
   double *d;
@@ -240,7 +241,7 @@ ui_add_to_print_in_row(struct argp_option *option, char 
*arg,
         }
 
       /* Add this option to the print list. */
-      gal_linkedlist_add_to_ill(&p->toprint, option->key);
+      gal_linkedlist_add_to_ill(&p->singlevalue, option->key);
     }
   else
     {
@@ -264,7 +265,7 @@ ui_add_to_print_in_row(struct argp_option *option, char 
*arg,
                               "and 1, you had asked for %g (read from `%s')",
                               d[i], arg);
               gal_linkedlist_add_to_dll(&p->tp_args, d[i]);
-              gal_linkedlist_add_to_ill(&p->toprint, option->key);
+              gal_linkedlist_add_to_ill(&p->singlevalue, option->key);
             }
           break;
 
@@ -441,6 +442,22 @@ ui_read_check_only_options(struct statisticsparams *p)
   gal_table_check_fits_format(p->cp.output, p->cp.tableformat);
 
 
+  /* If in tile-mode, we must have at least one single valued option. */
+  if(p->ontile && p->singlevalue==NULL)
+    error(EXIT_FAILURE, 0, "at least one of the single-value measurements "
+          "(for example `--median') must be requested with the `--ontile' "
+          "option: there is no value to put in each tile");
+
+  /* The tile mode cannot be called with any other modes. */
+  if(p->ontile && (p->asciihist || p->asciicfp || p->histogram
+                   || p->cumulative || !isnan(p->sigclipmultip)
+                   || !isnan(p->mirror) ) )
+    error(EXIT_FAILURE, 0, "`--ontile' cannot be called with any of the "
+          "`particular' calculation options, for example `--histogram'. "
+          "This is because these options deal with the whole dataset, "
+          "but in tessellation mode, the input dataset is broken up into "
+          "individual tiles");
+
   /* If less than and greater than are both given, make sure that the value
      to greater than is smaller than the value to less-than. */
   if( !isnan(p->lessthan) && !isnan(p->greaterequal)
@@ -476,7 +493,7 @@ ui_read_check_only_options(struct statisticsparams *p)
 
   /* Reverse the list of statistics to print in one row and also the
      arguments, so it has the same order the user wanted. */
-  gal_linkedlist_reverse_ill(&p->toprint);
+  gal_linkedlist_reverse_ill(&p->singlevalue);
   gal_linkedlist_reverse_dll(&p->tp_args);
 }
 
@@ -656,7 +673,7 @@ ui_make_sorted_if_necessary(struct statisticsparams *p)
   struct gal_linkedlist_ill *tmp;
 
   /* Check in the one-row outputs. */
-  for(tmp=p->toprint; tmp!=NULL; tmp=tmp->next)
+  for(tmp=p->singlevalue; tmp!=NULL; tmp=tmp->next)
     switch(tmp->v)
       {
       case ARGS_OPTION_KEY_MODE:
@@ -763,29 +780,61 @@ ui_read_columns(struct statisticsparams *p)
 void
 ui_preparations(struct statisticsparams *p)
 {
+  struct gal_options_common_params *cp=&p->cp;
+
   /* Read the input. */
   if(p->isfits && p->hdu_type==IMAGE_HDU)
-    p->input=gal_fits_img_read(p->inputname, p->cp.hdu, p->cp.minmapsize);
-  else ui_read_columns(p);
+    {
+      p->inputformat=INPUT_FORMAT_IMAGE;
+      p->input=gal_fits_img_read(p->inputname, cp->hdu, cp->minmapsize);
+    }
+  else
+    {
+      ui_read_columns(p);
+      p->inputformat=INPUT_FORMAT_TABLE;
+    }
+
+  /* Tile and channel sanity checks. */
+  if(p->ontile)
+    {
+      /* Check the tiles. */
+      cp->channelsize=gal_tile_all_sanity_check(p->inputname, p->cp.hdu,
+                                                p->input, cp->tilesize,
+                                                cp->numchannels);
+
+      /* Set the tile file name if the user wants to check the tiles. */
+      if(cp->checktiles)
+        {
+          cp->tilecheckname=gal_checkset_automatic_output(&p->cp,
+                                                          p->inputname,
+                                                          "_tiled.fits");
+          gal_checkset_check_remove_file(cp->tilecheckname, 0,
+                                         cp->dontdelete);
+        }
+    }
 
   /* Set the out-of-range values in the input to blank. */
   ui_out_of_range_to_blank(p);
 
-  /* Only keep the numbers we want. */
-  gal_blank_remove(p->input);
+  /* If we are not to work on tiles: */
+  if(!p->ontile)
+    {
+      /* Only keep the numbers we want. */
+      gal_blank_remove(p->input);
 
-  /* Make sure there is any data remaining: */
-  if(p->input->size==0)
-    error(EXIT_FAILURE, 0, "%s: no data, maybe the `--greaterequal' or "
-          "`--lessthan' options need to be adjusted",
-          gal_fits_name_save_as_string(p->inputname, p->cp.hdu) );
+      /* Make sure there is any data remaining: */
+      if(p->input->size==0)
+        error(EXIT_FAILURE, 0, "%s: no data, maybe the `--greaterequal' or "
+              "`--lessthan' options need to be adjusted",
+              gal_fits_name_save_as_string(p->inputname, cp->hdu) );
 
-  /* Make the sorted array if necessary. */
-  ui_make_sorted_if_necessary(p);
+      /* Make the sorted array if necessary. */
+      ui_make_sorted_if_necessary(p);
 
-  /* Set the number of output files. */
-  if( !isnan(p->mirror) )             ++p->numoutfiles;
-  if( p->histogram || p->cumulative ) ++p->numoutfiles;
+      /* Set the number of output files. */
+      if( !isnan(p->mirror) )             ++p->numoutfiles;
+      if( p->histogram || p->cumulative ) ++p->numoutfiles;
+    }
 }
 
 
@@ -894,6 +943,6 @@ ui_free_report(struct statisticsparams *p)
     gal_data_free(p->sorted);
   gal_data_free(p->input);
   gal_data_free(p->reference);
-  gal_linkedlist_free_ill(p->toprint);
   gal_linkedlist_free_dll(p->tp_args);
+  gal_linkedlist_free_ill(p->singlevalue);
 }
diff --git a/bin/statistics/ui.h b/bin/statistics/ui.h
index 7838859..c385917 100644
--- a/bin/statistics/ui.h
+++ b/bin/statistics/ui.h
@@ -29,7 +29,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 /* Available letters for short options:
 
-   a b d e f i j k p v w x y z
+   a b e f i j k p v w x y z
    B G J L R W X Y
 */
 enum option_keys_enum
@@ -41,7 +41,7 @@ enum option_keys_enum
   ARGS_OPTION_KEY_LESSTHAN     = 'l',
   ARGS_OPTION_KEY_QRANGE       = 'Q',
   ARGS_OPTION_KEY_MEAN         = 'm',
-  ARGS_OPTION_KEY_STD          = 't',
+  ARGS_OPTION_KEY_STD          = 'd',
   ARGS_OPTION_KEY_MEDIAN       = 'E',
   ARGS_OPTION_KEY_MODE         = 'O',
   ARGS_OPTION_KEY_QUANTILE     = 'u',
@@ -50,6 +50,7 @@ enum option_keys_enum
   ARGS_OPTION_KEY_CUMULATIVE   = 'C',
   ARGS_OPTION_KEY_SIGMACLIP    = 's',
   ARGS_OPTION_KEY_NORMALIZE    = 'n',
+  ARGS_OPTION_KEY_ONTILE       = 't',
 
   /* Only with long version (start with a value 1000, the rest will be set
      automatically). */
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 2622685..fa627a6 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -6396,9 +6396,9 @@ HDU, or modify the keywords, add new kewords or delete 
some others.
 @node Invoking astfits,  , Fits, Fits
 @subsection Invoking Fits
 
-Fits can print or manipulate the FITS file HDUs (extensions), or print an
-manipulate the meta-data keywords in a given HDU. The executable name is
address@hidden with the following general template
+Fits can print or manipulate the FITS file HDUs (extensions), meta-data
+keywords in a given HDU. The executable name is @file{astfits} with the
+following general template
 
 @example
 $ astfits [OPTION...] ASTRdata
@@ -6410,8 +6410,8 @@ One line examples:
 
 @example
 $ astfits image.fits
-$ astfits image.fits --printall
-$ astfits image.fits --printall | grep NAXIS
+$ astfits image.fits --printallkeys
+$ astfits image.fits --printallkeys | grep NAXIS
 $ astfits input.fits --copy=hdu-name --output=out.fits
 $ astfits --update=OLDKEY,153.034,"Old keyword comment"
 $ astfits --delete=COMMENT --comment="Anything you like ;-)."
@@ -6419,17 +6419,18 @@ $ astfits --write=MYKEY1,20.00,"An example keyword" 
--write=MYKEY2,fd
 @end example
 
 @cindex HDU
-When no paritcular HDU or keyword options are given, Fits will print a list
-of information about the extension(s) in the file. This information
-includes the HDU number, name, along with the type, and size the data it
-contains (dimension and size along each dimension). You can use this to get
-a general idea of the contents of the FITS file and what HDU to use for
-further processing, either with the Fits program or any other Gnuastro
-program.
+When no action is requested (and only a file name is given), Fits will
+print a list of information about the extension(s) in the file. This
+information includes the HDU number, HDU name (@code{EXTNAME} keyword),
+type of data (see @ref{Numeric data types}, and the number of data elements
+it contains (size along each dimension for images and table rows and
+columns). You can use this to get a general idea of the contents of the
+FITS file and what HDU to use for further processing, either with the Fits
+program or any other Gnuastro program.
 
 Here is one example of information about a FITS file with four exensions:
-the first extension has no data, it is just meta-data (commonly used to
-keep meta-data about the whole file, or grouping of extensions, see
+the first extension has no data, it is a purely meta-data HDU (commonly
+used to keep meta-data about the whole file, or grouping of extensions, see
 @ref{Fits}). The second extension is an image with name @code{IMAGE} and
 single precision floating point type (@code{float32}, see @ref{Numeric data
 types}), it has 4287 pixels along its first (horizontal) axis and 4286
diff --git a/lib/commonopts.h b/lib/commonopts.h
index 84581c6..59ea715 100644
--- a/lib/commonopts.h
+++ b/lib/commonopts.h
@@ -60,7 +60,7 @@ struct argp_option gal_commonopts_options[] =
       GAL_OPTIONS_KEY_SEARCHIN,
       "STR",
       0,
-      "Col. selection field: `name', `unit', `comment'.",
+      "Select column(s) in: `name', `unit', `comment'.",
       GAL_OPTIONS_GROUP_INPUT,
       &cp->searchin,
       GAL_DATA_TYPE_STRING,
@@ -97,7 +97,7 @@ struct argp_option gal_commonopts_options[] =
       GAL_OPTIONS_KEY_TILESIZE,
       "INT[,INT]",
       0,
-      "Size of regular tiles on each dim. (FITS order).",
+      "Regular tile size on each dim. (FITS order).",
       GAL_OPTIONS_GROUP_TESSELLATION,
       &cp->tilesize,
       GAL_DATA_TYPE_SIZE_T,
@@ -138,7 +138,7 @@ struct argp_option gal_commonopts_options[] =
       GAL_OPTIONS_KEY_WORKOVERCH,
       0,
       0,
-      "Convolve over channel borders.",
+      "Work (not tile) over channel edges.",
       GAL_OPTIONS_GROUP_TESSELLATION,
       &cp->workoverch,
       GAL_OPTIONS_NO_ARG_TYPE,
diff --git a/lib/gnuastro/tile.h b/lib/gnuastro/tile.h
index 45d6b63..30f3474 100644
--- a/lib/gnuastro/tile.h
+++ b/lib/gnuastro/tile.h
@@ -73,9 +73,8 @@ size_t
 gal_tile_block_increment(gal_data_t *block, size_t *tsize,
                          size_t num_increment, size_t *coord);
 
-void
-gal_tile_block_check_tiles(gal_data_t *tiles, char *filename,
-                           char *program_name);
+gal_data_t *
+gal_tile_block_check_tiles(gal_data_t *tiles);
 
 
 
@@ -87,11 +86,11 @@ size_t *
 gal_tile_all_sanity_check(char *filename, char *hdu, gal_data_t *input,
                           size_t *tile, size_t *numchannels);
 
-size_t
+size_t *
 gal_tile_all_position(gal_data_t *input, size_t *regular,
                       float remainderfrac, gal_data_t **out, size_t multiple);
 
-void
+size_t *
 gal_tile_all_position_two_layers(gal_data_t *input, size_t *channel_size,
                                  size_t *tile_size, float remainderfrac,
                                  gal_data_t **channels, gal_data_t **tiles);
diff --git a/lib/options.h b/lib/options.h
index 448a495..86ddb11 100644
--- a/lib/options.h
+++ b/lib/options.h
@@ -192,6 +192,7 @@ struct gal_options_common_params
   /* For internal (to option processing) purposes. */
   uint8_t                 keep; /* Output file can exist.                 */
   size_t          *channelsize; /* Size of channels along each dimension. */
+  char          *tilecheckname; /* Name of file to check tiles.           */
   void         *program_struct; /* Host program's main variable struct.   */
   char           *program_name; /* Official name to be used in text.      */
   char           *program_exec; /* Program's executable name.             */
diff --git a/lib/tile.c b/lib/tile.c
index 1fceecf..bab614f 100644
--- a/lib/tile.c
+++ b/lib/tile.c
@@ -318,9 +318,8 @@ gal_tile_block_increment(gal_data_t *block, size_t *tsize,
    mainly for inspecting the positioning of tiles. We are using a signed
    32-bit type because this is the standard FITS standard type for
    integers. */
-void
-gal_tile_block_check_tiles(gal_data_t *tiles, char *filename,
-                           char *program_name)
+gal_data_t *
+gal_tile_block_check_tiles(gal_data_t *tiles)
 {
   size_t num_increment;
   gal_data_t *tofill, *tile;
@@ -381,12 +380,9 @@ gal_tile_block_check_tiles(gal_data_t *tiles, char 
*filename,
       ++tile_index;
     }
 
-  /* Save the indexs into a file. */
-  gal_fits_img_write(tofill, filename, NULL, program_name);
-
   /* Clean up. */
   free(coord);
-  gal_data_free(tofill);
+  return tofill;
 }
 
 
@@ -531,7 +527,7 @@ gal_tile_all_sanity_check(char *filename, char *hdu, 
gal_data_t *input,
 static void
 gal_tile_all_regular_first(gal_data_t *parent, size_t *regular,
                            float significance, size_t *first, size_t *last,
-                           size_t *number)
+                           size_t *tsize)
 {
   size_t i, remainder, *dsize=parent->dsize;;
 
@@ -547,20 +543,20 @@ gal_tile_all_regular_first(gal_data_t *parent, size_t 
*regular,
           if( remainder > significance * regular[i] )
             {
               first[i]  = ( remainder + regular[i] )/2;
-              number[i] = dsize[i]/regular[i] + 1 ;
-              last[i]   = dsize[i] - ( first[i] + regular[i]*(number[i]-2) );
+              tsize[i] = dsize[i]/regular[i] + 1 ;
+              last[i]   = dsize[i] - ( first[i] + regular[i]*(tsize[i]-2) );
             }
           else
             {
               first[i]  = remainder + regular[i];
-              number[i] = dsize[i]/regular[i];
+              tsize[i]  = dsize[i]/regular[i];
               last[i]   = regular[i];
             }
         }
       else
         {
           first[i]  = last[i] = regular[i];
-          number[i] = dsize[i]/regular[i];
+          tsize[i] = dsize[i]/regular[i];
         }
     }
 }
@@ -586,6 +582,10 @@ gal_tile_all_regular_first(gal_data_t *parent, size_t 
*regular,
         dimensions. So it must have the same number of elements as the
         dimensions of `input'.
 
+     `remainderfrac' is the significant fraction of the remainder space if
+        the width of the input isn't an exact multiple of the tile size
+        along a dimension, see `gal_tile_all_regular_first'.
+
      `out' is the pointer to the array of data structures that is to keep
         the tile parameters. If `*out==NULL', then the necessary space will
         be allocated. If it is not NULL, then all the tile information will
@@ -602,9 +602,9 @@ gal_tile_all_regular_first(gal_data_t *parent, size_t 
*regular,
    Output
    ------
 
-     The returned output is an array of `gal_data_t' with `numtiles'
-     elements, note that you have to pass the pointer to `numtiles' as one
-     of the arguments.
+     The returned output is an array of numbers (the same size as the input
+        data structure's dimensions) keeping the number of tiles along each
+        dimension.
 
 
    Implementation
@@ -620,7 +620,7 @@ gal_tile_all_regular_first(gal_data_t *parent, size_t 
*regular,
         2. parent-size (or `psize') which is the size of the parent in each
            dimension (we don't want to go out of the paren't range).
 */
-size_t
+size_t *
 gal_tile_all_position(gal_data_t *input, size_t *regular,
                       float remainderfrac, gal_data_t **out, size_t multiple)
 {
@@ -725,12 +725,11 @@ gal_tile_all_position(gal_data_t *input, size_t *regular,
 
   /* Clean up and return. */
   free(last);
-  free(tsize);
   free(first);
   free(coord);
   free(tcoord);
   if(start) free(start);
-  return numtiles;
+  return tsize;
 }
 
 
@@ -749,27 +748,20 @@ gal_tile_all_position(gal_data_t *input, size_t *regular,
       Tiles: A combined tesselation of each channel with smaller
            tiles. These tiles can be used to calculate things like
            gradients over each channel and thus over the whole image.  */
-void
+size_t *
 gal_tile_all_position_two_layers(gal_data_t *input, size_t *channel_size,
                                  size_t *tile_size, float remainderfrac,
                                  gal_data_t **channels, gal_data_t **tiles)
 {
   gal_data_t *ch, *t;
-  size_t i, nch, ntiles_in_ch;
+  size_t i, nch=1, ntiles_in_ch;
+  size_t *chsize, *tsize, *ttsize;
 
-  /* First allocate the channels. Note that the channels tesselation might
-     have already been set. */
-  if(*channels)
-    {
-      nch=1;
-      for(i=0;i<input->ndim;++i)
-        nch *= input->dsize[i]/channel_size[i];
-    }
-  else
-    /* Note that the actual allocated input array will be the direct
-       `block' of each channel. */
-    nch=gal_tile_all_position(input, channel_size, remainderfrac,
-                              channels, 1);
+  /* First allocate the channels tessellation. */
+  *channels=NULL;
+  chsize=gal_tile_all_position(input, channel_size, remainderfrac,
+                               channels, 1);
+  for(i=0;i<input->ndim;++i) nch *= chsize[i];
 
 
   /* Now, tile each channel. While tiling the first channel, we are also
@@ -777,8 +769,9 @@ gal_tile_all_position_two_layers(gal_data_t *input, size_t 
*channel_size,
      pointers. */
   *tiles=NULL;
   ch=*channels;
-  ntiles_in_ch = gal_tile_all_position(ch, tile_size, remainderfrac,
-                                       tiles, nch);
+  tsize = gal_tile_all_position(ch, tile_size, remainderfrac,
+                                tiles, nch);
+  ntiles_in_ch=gal_multidim_total_size(input->ndim, tsize);
   for(i=1;i<nch;++i)
     {
       /* Set the first tile in this channel. Then use it it fill the `next'
@@ -788,8 +781,17 @@ gal_tile_all_position_two_layers(gal_data_t *input, size_t 
*channel_size,
       (*tiles)[ i * ntiles_in_ch - 1 ].next = t;
 
       /* Fill in the information for all the tiles in this channel. */
-      gal_tile_all_position(&ch[i], tile_size, remainderfrac, &t, 1);
+      ttsize=gal_tile_all_position(&ch[i], tile_size, remainderfrac, &t, 1);
+      free(ttsize);
     }
+
+  /* Multiply the number of tiles along each dimension OF ONE CHANNEL by
+     the number of channels in each dimension for the output. */
+  for(i=0;i<input->ndim;++i)
+    tsize[i] *= chsize[i];
+
+  /* Return the total number of tiles along each dimension. */
+  return tsize;
 }