[gnuastro-commits] master 63043e5 105/125: Statistical functions can accept tiled input

[Mohammad Akhlaghi]

branch: master
commit 63043e5cfa4e9900ef946fbffac2001973dd0b90
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>

    Statistical functions can accept tiled input
    
    The main statistical functions (currenly on the single measurement class in
    the Statistics program) are now equipt to deal with inputs as tiles. This
    will be correced in the next commit. This was made possible thanks to the
    new `GAL_TILE_PARSE_OPERATE' macro that will parse a tile and do a given
    operation. However, when channels are defined, the array keeping a single
    value for a tile will not be ordered properly.
    
    In the process, the following modifications were also made:
    
     - Statistics now takes a `mirrordist' option for mode-related
       calculations. Until now, this was hardcoded.
    
     - `gal_blank_present' can now work with tiles also, infact this was the
       first function to be modified for this job.
    
     - `gal_tile_to_contiguous' will copy the contents of a tile into a
       contiguous patch of memory.
    
     - `gal_statistics_no_blank_sorted' now also puts the tile's values into a
       contiguous piece of memory before removing blanks and sorting! In
       principle, these two operations cannot be defined in a non-contiguous
       patch of memory.
---
 bin/statistics/args.h             |  13 ++
 bin/statistics/aststatistics.conf |   3 +-
 bin/statistics/main.h             |   1 +
 bin/statistics/statistics.c       | 111 +++++++++++++++--
 bin/statistics/ui.c               |  23 +++-
 bin/statistics/ui.h               |   1 +
 lib/blank.c                       |  43 +++++--
 lib/gnuastro/statistics.h         |   2 +-
 lib/gnuastro/tile.h               | 121 +++++++++++++++++-
 lib/statistics.c                  | 250 +++++++++++++++++++++-----------------
 lib/tile.c                        |  57 +++++++++
 11 files changed, 486 insertions(+), 139 deletions(-)

diff --git a/bin/statistics/args.h b/bin/statistics/args.h
index ebb3d37..3d3f099 100644
--- a/bin/statistics/args.h
+++ b/bin/statistics/args.h
@@ -478,6 +478,19 @@ struct argp_option program_options[] =
       GAL_OPTIONS_NOT_MANDATORY,
       GAL_OPTIONS_NOT_SET
     },
+    {
+      "mirrordist",
+      ARGS_OPTION_KEY_MIRRORDIST,
+      "FLT",
+      0,
+      "Maximum dist. (in multip. of error) to find mode.",
+      ARGS_GROUP_HIST_CFP,
+      &p->mirrordist,
+      GAL_DATA_TYPE_FLOAT32,
+      GAL_OPTIONS_RANGE_ANY,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
+    },
 
 
     {0}
diff --git a/bin/statistics/aststatistics.conf 
b/bin/statistics/aststatistics.conf
index 46903fe..25810e3 100644
--- a/bin/statistics/aststatistics.conf
+++ b/bin/statistics/aststatistics.conf
@@ -22,4 +22,5 @@
 # Histogram and CFP settings
  numasciibins        70
  asciiheight         10
- numbins            100
\ No newline at end of file
+ numbins            100
+ mirrordist         1.5
\ No newline at end of file
diff --git a/bin/statistics/main.h b/bin/statistics/main.h
index 8c74ce5..a409f8c 100644
--- a/bin/statistics/main.h
+++ b/bin/statistics/main.h
@@ -80,6 +80,7 @@ struct statisticsparams
   uint8_t        normalize;  /* set the sum of all bins to 1.            */
   float        onebinstart;  /* Shift bins to start at this value.       */
   uint8_t        maxbinone;  /* Set the maximum bin to 1.                */
+  float         mirrordist;  /* Maximum distance after mirror for mode.  */
 
   /* Internal */
   uint8_t      inputformat;  /* Format of input dataset.                 */
diff --git a/bin/statistics/statistics.c b/bin/statistics/statistics.c
index f7d692d..fa04ad4 100644
--- a/bin/statistics/statistics.c
+++ b/bin/statistics/statistics.c
@@ -71,7 +71,7 @@ statistics_read_check_args(struct statisticsparams *p)
   double d;
   if(p->tp_args==NULL)
     error(EXIT_FAILURE, 0, "a bug! Not enough arguments for the "
-          "requested options to print in one row. Please contact "
+          "requested single measurement options. Please contact "
           "us at %s so we can address the problem",
           PACKAGE_BUGREPORT);
   gal_linkedlist_pop_from_dll(&p->tp_args, &d);
@@ -89,7 +89,6 @@ statistics_print_one_row(struct statisticsparams *p)
   char *toprint;
   size_t dsize=1;
   double arg, *d;
-  float mirrdist=1.5;
   struct gal_linkedlist_ill *tmp;
   gal_data_t *tmpv, *out=NULL, *num=NULL, *min=NULL, *max=NULL;
   gal_data_t *sum=NULL, *med=NULL, *meanstd=NULL, *modearr=NULL;
@@ -123,7 +122,7 @@ statistics_print_one_row(struct statisticsparams *p)
       case ARGS_OPTION_KEY_MODESYM:
       case ARGS_OPTION_KEY_MODESYMVALUE:
         modearr = ( modearr ? modearr
-                    : gal_statistics_mode(p->sorted, mirrdist, 0) );
+                    : gal_statistics_mode(p->sorted, p->mirrordist, 0) );
         d=modearr->array;
         if(d[2]<GAL_STATISTICS_MODE_GOOD_SYM) d[0]=d[1]=NAN;
         break;
@@ -172,7 +171,7 @@ statistics_print_one_row(struct statisticsparams *p)
         /* Not previously calculated. */
         case ARGS_OPTION_KEY_QUANTILE:
           arg = statistics_read_check_args(p);
-          out = gal_statistsics_quantile(p->sorted, arg, 0);
+          out = gal_statistics_quantile(p->sorted, arg, 0);
           break;
 
         case ARGS_OPTION_KEY_QUANTFUNC:
@@ -234,11 +233,18 @@ statistics_print_one_row(struct statisticsparams *p)
 void
 statistics_on_tile(struct statisticsparams *p)
 {
+  double arg;
   uint8_t type;
-  size_t i, numtiles=0, ntiles, *tsize;
+  gal_data_t *tmp, *tmpv, *ttmp;
   struct gal_linkedlist_ill *operation;
+  size_t ntiles, *tsize, tind, dsize=1, mind;
   gal_data_t *channels, *tiles, *check, *tile, *values;
 
+  /* Set the output name. */
+  if(p->cp.output==NULL)
+    p->cp.output=gal_checkset_automatic_output(&p->cp, p->inputname,
+                                               "_ontile.fits");
+
   /* Make the tile structure. */
   tsize=gal_tile_all_position_two_layers(p->input, p->cp.channelsize,
                                          p->cp.tilesize, p->cp.remainderfrac,
@@ -270,18 +276,18 @@ statistics_on_tile(struct statisticsparams *p)
 
         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_SUM:
         case ARGS_OPTION_KEY_MEAN:
         case ARGS_OPTION_KEY_STD:
+        case ARGS_OPTION_KEY_QUANTILE:
         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:
@@ -293,7 +299,93 @@ statistics_on_tile(struct statisticsparams *p)
       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);
+      /* Read the argument for those operations that need it. This is done
+         here, because below, the functions are repeated on each tile. */
+      switch(operation->v)
+        {
+        case ARGS_OPTION_KEY_QUANTILE:
+          arg = statistics_read_check_args(p);
+          break;
+        case ARGS_OPTION_KEY_QUANTFUNC:
+          arg = statistics_read_check_args(p);
+          tmpv = gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT64, 1, &dsize,
+                                NULL, 1, -1, NULL, NULL, NULL);
+          *((double *)(tmpv->array)) = arg;
+          tmpv = gal_data_copy_to_new_type_free(tmpv, p->input->type);
+        }
+
+      /* Do the operation on each tile. */
+      tind=0;
+      for(tile=tiles; tile!=NULL; tile=tile->next)
+        {
+          /* Do the proper operation. */
+          switch(operation->v)
+            {
+            case ARGS_OPTION_KEY_NUMBER:
+              tmp=gal_statistics_number(tile);                      break;
+
+            case ARGS_OPTION_KEY_MINIMUM:
+              tmp=gal_statistics_minimum(tile);                     break;
+
+            case ARGS_OPTION_KEY_MAXIMUM:
+              tmp=gal_statistics_maximum(tile);                     break;
+
+            case ARGS_OPTION_KEY_MEDIAN:
+              tmp=gal_statistics_median(tile, 1);                   break;
+
+            case ARGS_OPTION_KEY_QUANTFUNC:
+              tmp=gal_statistics_quantile_function(tile, tmpv, 1);  break;
+
+            case ARGS_OPTION_KEY_SUM:
+              tmp=gal_statistics_sum(tile);                         break;
+
+            case ARGS_OPTION_KEY_MEAN:
+              tmp=gal_statistics_mean(tile);                        break;
+
+            case ARGS_OPTION_KEY_STD:
+              tmp=gal_statistics_std(tile);                         break;
+
+            case ARGS_OPTION_KEY_QUANTILE:
+              tmp=gal_statistics_quantile(tile, arg, 1);            break;
+
+            case ARGS_OPTION_KEY_MODE:
+            case ARGS_OPTION_KEY_MODESYM:
+            case ARGS_OPTION_KEY_MODEQUANT:
+            case ARGS_OPTION_KEY_MODESYMVALUE:
+              switch(operation->v)
+                {
+                case ARGS_OPTION_KEY_MODE:         mind=0;  break;
+                case ARGS_OPTION_KEY_MODESYM:      mind=2;  break;
+                case ARGS_OPTION_KEY_MODEQUANT:    mind=1;  break;
+                case ARGS_OPTION_KEY_MODESYMVALUE: mind=3;  break;
+                }
+              tmp=gal_statistics_mode(tile, p->mirrordist, 1);
+              ttmp=statistics_pull_out_element(tmp, mind);
+              gal_data_free(tmp);
+              tmp=ttmp;
+              break;
+
+            default:
+              error(EXIT_FAILURE, 0, "a bug! please contact us at %s to "
+                    "fix the problem. THe operation code %d is not "
+                    "recognized in `statistics_on_tile'", PACKAGE_BUGREPORT,
+                    operation->v);
+            }
+
+          /* Put the output value into the `values' array and clean up. */
+          tmp=gal_data_copy_to_new_type_free(tmp, type);
+          memcpy(gal_data_ptr_increment(values->array, tind++,
+                                        values->type),
+                 tmp->array, gal_data_sizeof(type));
+          gal_data_free(tmp);
+        }
+
+      /* Save the values. */
+      gal_fits_img_write(values, p->cp.output, NULL, PROGRAM_STRING);
+      gal_data_free(values);
+
+      /* Clean up. */
+      if(operation->v==ARGS_OPTION_KEY_QUANTFUNC) gal_data_free(tmpv);
     }
 
 
@@ -301,9 +393,6 @@ statistics_on_tile(struct statisticsparams *p)
   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);
 }
 
 
diff --git a/bin/statistics/ui.c b/bin/statistics/ui.c
index 2cc0de5..9498ab0 100644
--- a/bin/statistics/ui.c
+++ b/bin/statistics/ui.c
@@ -133,6 +133,7 @@ ui_initialize_options(struct statisticsparams *p,
   p->quantmin            = NAN;
   p->quantmax            = NAN;
   p->mirror              = NAN;
+  p->mirrordist          = NAN;
   p->sigclipparam        = NAN;
   p->sigclipmultip       = NAN;
 
@@ -437,6 +438,8 @@ ui_parse_numbers(struct argp_option *option, char *arg,
 static void
 ui_read_check_only_options(struct statisticsparams *p)
 {
+  struct gal_linkedlist_ill *tmp;
+
   /* Check if the format of the output table is valid, given the type of
      the output. */
   gal_table_check_fits_format(p->cp.output, p->cp.tableformat);
@@ -458,6 +461,22 @@ ui_read_check_only_options(struct statisticsparams *p)
           "but in tessellation mode, the input dataset is broken up into "
           "individual tiles");
 
+  /* If any of the mode measurements are requested, then `mirrordist' is
+     mandatory. */
+  for(tmp=p->singlevalue; tmp!=NULL; tmp=tmp->next)
+    switch(tmp->v)
+      {
+      case ARGS_OPTION_KEY_MODE:
+      case ARGS_OPTION_KEY_MODESYM:
+      case ARGS_OPTION_KEY_MODEQUANT:
+      case ARGS_OPTION_KEY_MODESYMVALUE:
+        if( isnan(p->mirrordist) )
+          error(EXIT_FAILURE, 0, "`--mirrordist' is required for the "
+                "mode-related single measurements (`--mode', `--modequant', "
+                "`--modesym', and `--modesymvalue')");
+      }
+
+
   /* 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)
@@ -598,12 +617,12 @@ ui_out_of_range_to_blank(struct statisticsparams *p)
       if( isnan(p->quantmax) ) p->quantmax = 1 - p->quantmin;
 
       /* Set the greater-equal value. */
-      tmp=gal_statistsics_quantile(ref, p->quantmin, 1);
+      tmp=gal_statistics_quantile(ref, p->quantmin, 1);
       tmp=gal_data_copy_to_new_type_free(tmp, GAL_DATA_TYPE_FLOAT32);
       p->greaterequal=*((float *)(tmp->array));
 
       /* Set the lower-than value. */
-      tmp=gal_statistsics_quantile(ref, p->quantmax, 1);
+      tmp=gal_statistics_quantile(ref, p->quantmax, 1);
       tmp=gal_data_copy_to_new_type_free(tmp, GAL_DATA_TYPE_FLOAT32);
       p->lessthan=*((float *)(tmp->array));
     }
diff --git a/bin/statistics/ui.h b/bin/statistics/ui.h
index c385917..aa55c42 100644
--- a/bin/statistics/ui.h
+++ b/bin/statistics/ui.h
@@ -70,6 +70,7 @@ enum option_keys_enum
   ARGS_OPTION_KEY_LOWERBIN,
   ARGS_OPTION_KEY_ONEBINSTART,
   ARGS_OPTION_KEY_MAXBINONE,
+  ARGS_OPTION_KEY_MIRRORDIST,
 };
 
 
diff --git a/lib/blank.c b/lib/blank.c
index 2faed9b..27e7d21 100644
--- a/lib/blank.c
+++ b/lib/blank.c
@@ -29,6 +29,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <stdlib.h>
 
 #include <gnuastro/data.h>
+#include <gnuastro/tile.h>
 #include <gnuastro/blank.h>
 
 #include <checkset.h>
@@ -99,23 +100,41 @@ gal_blank_alloc_write(uint8_t type)
 
 
 /* Return 1 if the dataset has a blank value and zero if it doesn't. */
-#define HAS_BLANK(IT) {                                         \
-    IT b, *a=data->array, *af=a+data->size;                     \
-    gal_blank_write(&b, data->type);                            \
-    if(b==b) do if(*a==b)   return 1; while(++a<af);            \
-    else     do if(*a!=*a)  return 1; while(++a<af);            \
+#define HAS_BLANK(IT) {                                                 \
+    IT b, *a=input->array, *af=a+input->size, *start;                   \
+    gal_blank_write(&b, block->type);                                   \
+                                                                        \
+    /* If this is a tile, not a full block. */                          \
+    if(input!=block)                                                    \
+      start=gal_tile_start_end_ind_inclusive(input, block, start_end_inc); \
+                                                                        \
+    /* Go over all the elements. */                                     \
+    while( start_end_inc[0] + increment <= start_end_inc[1] )           \
+      {                                                                 \
+        if(input!=block)                                                \
+          af = ( a = start + increment ) + input->dsize[input->ndim-1]; \
+        if(b==b) do if(*a==b)  return 1; while(++a<af);                 \
+        else     do if(*a!=*a) return 1; while(++a<af);                 \
+        if(input!=block)                                                \
+          increment += gal_tile_block_increment(block, input->dsize,    \
+                                                num_increment++, NULL); \
+        else break;                                                     \
+      }                                                                 \
   }
 int
-gal_blank_present(gal_data_t *data)
+gal_blank_present(gal_data_t *input)
 {
-  char **str=data->array, **strf=str+data->size;
+  char **str, **strf;
+  size_t increment=0, num_increment=1;
+  gal_data_t *block=gal_tile_block(input);
+  size_t start_end_inc[2]={0,block->size};
 
   /* If there is nothing in the array (its size is zero), then return 0 (no
      blank is present. */
-  if(data->size==0) return 0;
+  if(input->size==0) return 0;
 
   /* Go over the pixels and check: */
-  switch(data->type)
+  switch(block->type)
     {
     /* Numeric types */
     case GAL_DATA_TYPE_UINT8:     HAS_BLANK( uint8_t  );    break;
@@ -131,6 +150,10 @@ gal_blank_present(gal_data_t *data)
 
     /* String. */
     case GAL_DATA_TYPE_STRING:
+      if(input!=block)
+        error(EXIT_FAILURE, 0, "tile mode is currently not supported for "
+              "strings");
+      strf = (str=input->array) + input->size;
       do if(!strcmp(*str++,GAL_BLANK_STRING)) return 1; while(str<strf);
       break;
 
@@ -147,7 +170,7 @@ gal_blank_present(gal_data_t *data)
 
     default:
       error(EXIT_FAILURE, 0, "a bug! type value (%d) not recognized "
-            "in `gal_blank_present'", data->type);
+            "in `gal_blank_present'", block->type);
     }
 
   /* If there was a blank value, then the function would have returned with
diff --git a/lib/gnuastro/statistics.h b/lib/gnuastro/statistics.h
index de6a9d0..cecab44 100644
--- a/lib/gnuastro/statistics.h
+++ b/lib/gnuastro/statistics.h
@@ -105,7 +105,7 @@ gal_data_t *
 gal_statistics_median(gal_data_t *input, int inplace);
 
 gal_data_t *
-gal_statistsics_quantile(gal_data_t *input, double quantile, int inplace);
+gal_statistics_quantile(gal_data_t *input, double quantile, int inplace);
 
 size_t
 gal_statistics_quantile_index(size_t size, double quant);
diff --git a/lib/gnuastro/tile.h b/lib/gnuastro/tile.h
index 30f3474..62a5f02 100644
--- a/lib/gnuastro/tile.h
+++ b/lib/gnuastro/tile.h
@@ -46,8 +46,6 @@ __BEGIN_C_DECLS  /* From C++ preparations */
 
 
 
-
-
 /***********************************************************************/
 /**************              Single tile              ******************/
 /***********************************************************************/
@@ -61,6 +59,9 @@ void *
 gal_tile_start_end_ind_inclusive(gal_data_t *tile, gal_data_t *work,
                                  size_t *start_end_inc);
 
+gal_data_t *
+gal_tile_to_contiguous(gal_data_t *input);
+
 
 
 /***********************************************************************/
@@ -116,6 +117,122 @@ gal_tile_function_on_threads(gal_data_t *tiles, void 
*(*meshfunc)(void *),
 
 
 
+
+/***********************************************************************/
+/**************           Function-like macros        ******************/
+/***********************************************************************/
+/* Parse over a region of memory (can be an n-dimensional tile or a fully
+   allocated block of memory) and do a certain operation. If `OUT' is not
+   NULL, it can also save the output into the `gal_data_t' that it points
+   to. Note that it must either have only one element (for the whole input)
+   or have exactly the same number of elements as the input (one value for
+   one pixel/element of the input). The input arguments are:
+
+       `IT': input type (C type) for example `int32_t' or `float'.
+
+       `OT': output type (C type) for example `int32_t' or `float'.
+
+       `OP': Operator: this can be any number of C experssions. This macro
+             is going to define an `IT *i' variable which will increment
+             over each element of the input array/tile. It will also define
+             an `OT *o' which you can use to access the output. If
+             `PARSE_OUT' is non-zero, then `o' will also be incremented to
+             the same index element but in the output array. You can use
+             these along with any other variable you define before this
+             macro to process the input and store it in the output. Note
+             that `i' and `o' will be incremented once in here, so don't
+             increment them in `OP'.
+
+       `IN': Input `gal_data_t', can be a tile or an allocated block of
+             memory.
+
+       `OUT': Output `gal_data_t'. It can be NULL. In that case, `o' will
+             be NULL and should not be used.
+
+       `PARSE_OUT': Parse the output along with the input. When this is
+             non-zero, then the `o' pointer (described in `OP') will be
+             incremented to cover the same region of the input in the
+             output.
+
+   See `lib/statistics.c' for some example applications of this function.
+*/
+#define GAL_TILE_PARSE_OPERATE(IT, OT, OP, IN, OUT, PARSE_OUT) {        \
+    OT *ost, *o = OUT ? OUT->array : NULL;                              \
+    gal_data_t *iblock=gal_tile_block(IN);                              \
+    IT b, *st, *i=IN->array, *f=i+IN->size;                             \
+    gal_data_t *oblock = OUT ? gal_tile_block(OUT) : NULL;              \
+    size_t increment=0, num_increment=1, s_e_ind[2]={0,iblock->size};   \
+    int hasblank=gal_blank_present(IN), parse_out=(OUT && PARSE_OUT);   \
+                                                                        \
+    /* Write the blank value for the input type into `b'). Note that */ \
+    /* a tile doesn't necessarily have to have a type. */               \
+    gal_blank_write(&b, iblock->type);                                  \
+                                                                        \
+    /* A small sanity check. */                                         \
+    if( parse_out && gal_data_dsize_is_different(iblock, oblock) )      \
+      error(EXIT_FAILURE, 0, "when `PARSE_OUT' is non-zero, the "       \
+            "allocated block size of the input and output of "          \
+            "`GAL_TILE_PARSE_OPERATE' must be equal, but they are "     \
+            "not: %zu and %zu elements respectively)", iblock->size,    \
+            oblock->size);                                              \
+                                                                        \
+    /* If this is a tile, not a full block. */                          \
+    if(IN!=iblock)                                                      \
+      {                                                                 \
+        st=gal_tile_start_end_ind_inclusive(IN, iblock, s_e_ind);       \
+        if(parse_out)                                                   \
+          ost = (OT *)(oblock->array) + ( st - (IT *)(iblock->array) ); \
+      }                                                                 \
+                                                                        \
+    /* Go over contiguous patches of memory. */                         \
+    while( s_e_ind[0] + increment <= s_e_ind[1] )                       \
+      {                                                                 \
+        /* If we are on a tile, reset `a' and  `af'. Also, if there */  \
+        /* is more than one element in `OUT', then set that. Note   */  \
+        /* that it is upto the caller to increment `o' in `OP'.     */  \
+        if(IN!=iblock)                                                  \
+          {                                                             \
+            f = ( i = st + increment ) + IN->dsize[IN->ndim-1];         \
+            if(parse_out) o = ost + increment;                          \
+          }                                                             \
+                                                                        \
+        /* Do the operation depending the nature of the blank value. */ \
+        /* Recall that for integer types, the blank value must be    */ \
+        /* checked with `=='. But for floats, the blank value can be */ \
+        /* a NaN. Recall that a NAN will fail any comparison         */ \
+        /* including `=='. So when the blank value is not equal to   */ \
+        /* itself, then it is floating point and is a NAN. In that   */ \
+        /* case, the only way to check if a data element is blank or */ \
+        /* not is to see if the value of each element is equal to    */ \
+        /* itself or not. */                                            \
+        if(hasblank)                                                    \
+          {                                                             \
+            if(b==b) do if(*i!=b)  {OP; if(parse_out) ++o;} while(++i<f); \
+            else     do if(*i==*i) {OP; if(parse_out) ++o;} while(++i<f); \
+          }                                                             \
+        else         do            {OP; if(parse_out) ++o;} while(++i<f); \
+                                                                        \
+        /* Set the incrementation. On a fully allocated iblock (when */ \
+        /* `IN==iblock'), we have already gone through the whole     */ \
+        /* array, so we'll set the incrementation to the size of the */ \
+        /* while block which will stop the `while' loop above. On a  */ \
+        /* tile, we need to increment to the next contiguous patch   */ \
+        /* of memory to continue parsing this tile. */                  \
+        increment += ( IN==iblock ? iblock->size                        \
+                       : gal_tile_block_increment(iblock, IN->dsize,    \
+                                                  num_increment++,      \
+                                                  NULL) );              \
+      }                                                                 \
+                                                                        \
+    /* This is done in case the caller doesn't need `o' to avoid */     \
+    /* compiler warnings. */                                            \
+    o = o ? o+0 : NULL;                                                 \
+  }
+
+
+
+
+
 __END_C_DECLS    /* From C++ preparations */
 
 #endif
diff --git a/lib/statistics.c b/lib/statistics.c
index b8a72e5..5669cc9 100644
--- a/lib/statistics.c
+++ b/lib/statistics.c
@@ -31,6 +31,8 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <stdlib.h>
 
 #include <gnuastro/data.h>
+#include <gnuastro/tile.h>
+#include <gnuastro/fits.h>
 #include <gnuastro/blank.h>
 #include <gnuastro/qsort.h>
 #include <gnuastro/arithmetic.h>
@@ -53,41 +55,39 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 /* Return the number of non-blank elements in an array as a single element,
    uint64 type data structure. */
 #define STATS_NUM(IT) {                                                 \
-    IT b, *a=input->array, *af=a+input->size;                           \
-    gal_blank_write(&b, input->type);                                   \
-    if( gal_blank_present(input) )                                      \
-      {                                                                 \
-        if(b==b)   /* Blank can be found with the equal sign. */        \
-          do if(*a!=b)       ++(*num);            while(++a<af);        \
-        else       /* Blank is NaN (fails on equal).  */                \
-          do if(*a==*a)     ++(*num);            while(++a<af);         \
-      }                                                                 \
-    else *num=input->size;                                              \
+    GAL_TILE_PARSE_OPERATE(IT, uint64_t, {++(*o);}, input, out, 0)      \
   }
+
 gal_data_t *
 gal_statistics_number(gal_data_t *input)
 {
-  uint64_t *num;
   size_t dsize=1;
+  int type=gal_tile_block(input)->type;
   gal_data_t *out=gal_data_alloc(NULL, GAL_DATA_TYPE_UINT64, 1, &dsize,
                                  NULL, 1, -1, NULL, NULL, NULL);
-  num=out->array;
-  switch(input->type)
-    {
-    case GAL_DATA_TYPE_UINT8:     STATS_NUM( uint8_t  );    break;
-    case GAL_DATA_TYPE_INT8:      STATS_NUM( int8_t   );    break;
-    case GAL_DATA_TYPE_UINT16:    STATS_NUM( uint16_t );    break;
-    case GAL_DATA_TYPE_INT16:     STATS_NUM( int16_t  );    break;
-    case GAL_DATA_TYPE_UINT32:    STATS_NUM( uint32_t );    break;
-    case GAL_DATA_TYPE_INT32:     STATS_NUM( int32_t  );    break;
-    case GAL_DATA_TYPE_UINT64:    STATS_NUM( uint64_t );    break;
-    case GAL_DATA_TYPE_INT64:     STATS_NUM( int64_t  );    break;
-    case GAL_DATA_TYPE_FLOAT32:   STATS_NUM( float    );    break;
-    case GAL_DATA_TYPE_FLOAT64:   STATS_NUM( double   );    break;
-    default:
-      error(EXIT_FAILURE, 0, "type code %d not recognized in "
-            "`gal_statistics_number'", input->type);
-    }
+
+  /* If there is no blank values in the input, then the total number is
+     just the size. */
+  if(gal_blank_present(input))
+    switch(type)
+      {
+      case GAL_DATA_TYPE_UINT8:     STATS_NUM( uint8_t  );    break;
+      case GAL_DATA_TYPE_INT8:      STATS_NUM( int8_t   );    break;
+      case GAL_DATA_TYPE_UINT16:    STATS_NUM( uint16_t );    break;
+      case GAL_DATA_TYPE_INT16:     STATS_NUM( int16_t  );    break;
+      case GAL_DATA_TYPE_UINT32:    STATS_NUM( uint32_t );    break;
+      case GAL_DATA_TYPE_INT32:     STATS_NUM( int32_t  );    break;
+      case GAL_DATA_TYPE_UINT64:    STATS_NUM( uint64_t );    break;
+      case GAL_DATA_TYPE_INT64:     STATS_NUM( int64_t  );    break;
+      case GAL_DATA_TYPE_FLOAT32:   STATS_NUM( float    );    break;
+      case GAL_DATA_TYPE_FLOAT64:   STATS_NUM( double   );    break;
+      default:
+        error(EXIT_FAILURE, 0, "type code %d not recognized in "
+              "`gal_statistics_number'", type);
+      }
+  else
+    *((uint64_t *)(out->array)) = input->size;
+
   return out;
 }
 
@@ -100,24 +100,23 @@ gal_statistics_number(gal_data_t *input)
    comparison. So when finding the minimum or maximum, when the blank value
    is NaN, we can safely assume there is no blank value at all. */
 #define STATS_MIN(IT) {                                                 \
-    IT b, *o=out->array, *a=input->array, *af=a+input->size;            \
-    gal_blank_write(&b, input->type);                                   \
-    if( b==b && gal_blank_present(input) )                              \
-      do if(*a!=b) { *o = *a < *o ? *a : *o; ++n; } while(++a<af);      \
-    else                                                                \
-      do           { *o = *a < *o ? *a : *o; ++n; } while(++a<af);      \
-                                                                        \
-    /* If there were no usable elements, set the output to blank. */    \
-    if(n==0) *o=b;                                                      \
+    GAL_TILE_PARSE_OPERATE(IT, IT, {*o = *i < *o ? *i : *o; ++n;},      \
+                           input, out, 0)                               \
   }
+
 gal_data_t *
 gal_statistics_minimum(gal_data_t *input)
 {
   size_t dsize=1, n=0;
-  gal_data_t *out=gal_data_alloc(NULL, input->type, 1, &dsize,
-                                 NULL, 1, -1, NULL, NULL, NULL);
+  int type=gal_tile_block(input)->type;
+  gal_data_t *out=gal_data_alloc(NULL, type, 1, &dsize, NULL, 1, -1,
+                                 NULL, NULL, NULL);
+
+  /* Initialize the output with the maximum possible value. */
   gal_data_type_max(out->type, out->array);
-  switch(input->type)
+
+  /* Parse the full input. */
+  switch(type)
     {
     case GAL_DATA_TYPE_UINT8:     STATS_MIN( uint8_t  );    break;
     case GAL_DATA_TYPE_INT8:      STATS_MIN( int8_t   );    break;
@@ -131,8 +130,12 @@ gal_statistics_minimum(gal_data_t *input)
     case GAL_DATA_TYPE_FLOAT64:   STATS_MIN( double   );    break;
     default:
       error(EXIT_FAILURE, 0, "type code %d not recognized in "
-            "`gal_statistics_minimum'", input->type);
+            "`gal_statistics_minimum'", type);
     }
+
+  /* If there were no usable elements, set the output to blank, then
+     return. */
+  if(n==0) gal_blank_write(out->array, out->type);
   return out;
 }
 
@@ -143,24 +146,23 @@ gal_statistics_minimum(gal_data_t *input)
 /* Return the maximum (non-blank) value of a dataset in the same type as
    the dataset. See explanations of `gal_statistics_minimum'. */
 #define STATS_MAX(IT) {                                                 \
-    IT b, *o=out->array, *a=input->array, *af=a+input->size;            \
-    gal_blank_write(&b, input->type);                                   \
-    if( b==b && gal_blank_present(input) )                              \
-      do if(*a!=b) { *o = *a > *o ? *a : *o; ++n; } while(++a<af);      \
-    else                                                                \
-      do           { *o = *a > *o ? *a : *o; ++n; } while(++a<af);      \
-                                                                        \
-    /* If there were no usable elements, set the output to blank. */    \
-    if(n==0) *o=b;                                                      \
+    GAL_TILE_PARSE_OPERATE(IT, IT, {*o = *i > *o ? *i : *o; ++n;},      \
+                           input, out, 0)                               \
   }
+
 gal_data_t *
 gal_statistics_maximum(gal_data_t *input)
 {
   size_t dsize=1, n=0;
-  gal_data_t *out=gal_data_alloc(NULL, input->type, 1, &dsize,
+  int type=gal_tile_block(input)->type;
+  gal_data_t *out=gal_data_alloc(NULL, type, 1, &dsize,
                                  NULL, 1, -1, NULL, NULL, NULL);
+
+  /* Initialize the output with the minimum possible value. */
   gal_data_type_min(out->type, out->array);
-  switch(input->type)
+
+  /* Parse the full input. */
+  switch(type)
     {
     case GAL_DATA_TYPE_UINT8:     STATS_MAX( uint8_t  );    break;
     case GAL_DATA_TYPE_INT8:      STATS_MAX( int8_t   );    break;
@@ -174,8 +176,12 @@ gal_statistics_maximum(gal_data_t *input)
     case GAL_DATA_TYPE_FLOAT64:   STATS_MAX( double   );    break;
     default:
       error(EXIT_FAILURE, 0, "type code %d not recognized in "
-            "`gal_statistics_maximum'", input->type);
+            "`gal_statistics_maximum'", type);
     }
+
+  /* If there were no usable elements, set the output to blank, then
+     return. */
+  if(n==0) gal_blank_write(out->array, out->type);
   return out;
 }
 
@@ -186,29 +192,19 @@ gal_statistics_maximum(gal_data_t *input)
 /* Return the sum of the input dataset as a single element dataset of type
    float64. */
 #define STATS_SUM_NUM(IT) {                                             \
-    IT b, *a=input->array, *af=a+input->size;                           \
-    gal_blank_write(&b, input->type);                                   \
-    if( gal_blank_present(input) )                                      \
-      {                                                                 \
-        if(b==b)       /* Blank value can be checked with an equals. */ \
-          do if(*a!=b)       { ++n; sum += *a; }  while(++a<af);        \
-        else           /* Blank value will fail an equal comparison. */ \
-          do if(*a==*a)      { ++n; sum += *a; }  while(++a<af);        \
-      }                                                                 \
-    else                                                                \
-      {                                                                 \
-        do                          sum += *a;    while(++a<af);        \
-        n = input->size;                                                \
-      }                                                                 \
+    GAL_TILE_PARSE_OPERATE(IT, double, {++n; *o += *i;}, input, out, 0) \
   }
 gal_data_t *
 gal_statistics_sum(gal_data_t *input)
 {
-  double sum=0.0f;
   size_t dsize=1, n=0;
+  int type=gal_tile_block(input)->type;
   gal_data_t *out=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT64, 1, &dsize,
                                  NULL, 1, -1, NULL, NULL, NULL);
-  switch(input->type)
+
+  /* Parse the dataset. Note that in `gal_data_alloc' we set the `clear'
+     flag to 1, so it will be 0.0f. */
+  switch(type)
     {
     case GAL_DATA_TYPE_UINT8:     STATS_SUM_NUM( uint8_t  );    break;
     case GAL_DATA_TYPE_INT8:      STATS_SUM_NUM( int8_t   );    break;
@@ -222,9 +218,12 @@ gal_statistics_sum(gal_data_t *input)
     case GAL_DATA_TYPE_FLOAT64:   STATS_SUM_NUM( double   );    break;
     default:
       error(EXIT_FAILURE, 0, "type code %d not recognized in "
-            "`gal_statistics_maximum'", input->type);
+            "`gal_statistics_sum'", type);
     }
-  *((double *)(out->array)) = n ? sum : GAL_BLANK_FLOAT64;
+
+  /* If there were no usable elements, set the output to blank, then
+     return. */
+  if(n==0) gal_blank_write(out->array, out->type);
   return out;
 }
 
@@ -237,11 +236,14 @@ gal_statistics_sum(gal_data_t *input)
 gal_data_t *
 gal_statistics_mean(gal_data_t *input)
 {
-  double sum=0.0f;
   size_t dsize=1, n=0;
+  int type=gal_tile_block(input)->type;
   gal_data_t *out=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT64, 1, &dsize,
                                  NULL, 1, -1, NULL, NULL, NULL);
-  switch(input->type)
+
+  /* Parse the dataset. Note that in `gal_data_alloc' we set the `clear'
+     flag to 1, so it will be 0.0f. */
+  switch(type)
     {
     case GAL_DATA_TYPE_UINT8:     STATS_SUM_NUM( uint8_t  );    break;
     case GAL_DATA_TYPE_INT8:      STATS_SUM_NUM( int8_t   );    break;
@@ -255,9 +257,14 @@ gal_statistics_mean(gal_data_t *input)
     case GAL_DATA_TYPE_FLOAT64:   STATS_SUM_NUM( double   );    break;
     default:
       error(EXIT_FAILURE, 0, "type code %d not recognized in "
-            "`gal_statistics_maximum'", input->type);
+            "`gal_statistics_mean'", type);
     }
-  *((double *)(out->array)) = n ? sum/n : GAL_BLANK_FLOAT64;
+
+  /* Above, we calculated the sum and number, so if there were any elements
+     in the dataset (`n!=0'), divide the sum by the number, otherwise, put
+     a blank value in the output. */
+  if(n) *((double *)(out->array)) /= n;
+  else gal_blank_write(out->array, out->type);
   return out;
 }
 
@@ -268,29 +275,20 @@ gal_statistics_mean(gal_data_t *input)
 /* Return the standard deviation of the input dataset as a single element
    dataset of type float64. */
 #define STATS_NSS2(IT) {                                                \
-    IT b, *a=input->array, *af=a+input->size;                           \
-    gal_blank_write(&b, input->type);                                   \
-    if( gal_blank_present(input) )                                      \
-      {                                                                 \
-        if(b==b)       /* Blank value can be checked with an equals. */ \
-          do if(*a!=b) { ++n; s += *a; s2 += *a * *a; }  while(++a<af); \
-        else           /* Blank value will fail an equal comparison. */ \
-          do if(*a==*a){ ++n; s += *a; s2 += *a * *a; }  while(++a<af); \
-      }                                                                 \
-    else                                                                \
-      {                                                                 \
-        do             {      s += *a; s2 += *a * *a; }  while(++a<af); \
-        n = input->size;                                                \
-      }                                                                 \
+    GAL_TILE_PARSE_OPERATE(IT, double, {++n; s += *i; s2 += *i * *i;},  \
+                           input, out, 0)                               \
   }
 gal_data_t *
 gal_statistics_std(gal_data_t *input)
 {
   size_t dsize=1, n=0;
   double s=0.0f, s2=0.0f;
+  int type=gal_tile_block(input)->type;
   gal_data_t *out=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT64, 1, &dsize,
                                  NULL, 1, -1, NULL, NULL, NULL);
-  switch(input->type)
+
+  /* Parse the input. */
+  switch(type)
     {
     case GAL_DATA_TYPE_UINT8:     STATS_NSS2( uint8_t  );    break;
     case GAL_DATA_TYPE_INT8:      STATS_NSS2( int8_t   );    break;
@@ -304,8 +302,10 @@ gal_statistics_std(gal_data_t *input)
     case GAL_DATA_TYPE_FLOAT64:   STATS_NSS2( double   );    break;
     default:
       error(EXIT_FAILURE, 0, "type code %d not recognized in "
-            "`gal_statistics_maximum'", input->type);
+            "`gal_statistics_maximum'", type);
     }
+
+  /* Write the standard deviation into the output. */
   *((double *)(out->array)) = n ? sqrt( (s2-s*s/n)/n ) : GAL_BLANK_FLOAT64;
   return out;
 }
@@ -322,9 +322,12 @@ gal_statistics_mean_std(gal_data_t *input)
 {
   size_t dsize=2, n=0;
   double s=0.0f, s2=0.0f;
+  int type=gal_tile_block(input)->type;
   gal_data_t *out=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT64, 1, &dsize,
                                  NULL, 1, -1, NULL, NULL, NULL);
-  switch(input->type)
+
+  /* Parse the input. */
+  switch(type)
     {
     case GAL_DATA_TYPE_UINT8:     STATS_NSS2( uint8_t  );    break;
     case GAL_DATA_TYPE_INT8:      STATS_NSS2( int8_t   );    break;
@@ -338,8 +341,10 @@ gal_statistics_mean_std(gal_data_t *input)
     case GAL_DATA_TYPE_FLOAT64:   STATS_NSS2( double   );    break;
     default:
       error(EXIT_FAILURE, 0, "type code %d not recognized in "
-            "`gal_statistics_maximum'", input->type);
+            "`gal_statistics_maximum'", type);
     }
+
+  /* Write in the output values and return. */
   ((double *)(out->array))[0] = n ? s/n                  : GAL_BLANK_FLOAT64;
   ((double *)(out->array))[1] = n ? sqrt( (s2-s*s/n)/n ) : GAL_BLANK_FLOAT64;
   return out;
@@ -392,8 +397,8 @@ gal_statistics_median(gal_data_t *input, int inplace)
 {
   size_t dsize=1;
   gal_data_t *nbs=gal_statistics_no_blank_sorted(input, inplace);
-  gal_data_t *out=gal_data_alloc(NULL, input->type, 1, &dsize,
-                                 NULL, 1, -1, NULL, NULL, NULL);
+  gal_data_t *out=gal_data_alloc(NULL, nbs->type, 1, &dsize, NULL, 1, -1,
+                                 NULL, NULL, NULL);
   /* Write the median. */
   statistics_median_in_sorted_no_blank(nbs, out->array);
 
@@ -438,11 +443,11 @@ gal_statistics_quantile_index(size_t size, double quant)
 /* Return a single element dataset of the same type as input keeping the
    value that has the given quantile. */
 gal_data_t *
-gal_statistsics_quantile(gal_data_t *input, double quantile, int inplace)
+gal_statistics_quantile(gal_data_t *input, double quantile, int inplace)
 {
   size_t dsize=1, index;
   gal_data_t *nbs=gal_statistics_no_blank_sorted(input, inplace);
-  gal_data_t *out=gal_data_alloc(NULL, input->type, 1, &dsize,
+  gal_data_t *out=gal_data_alloc(NULL, nbs->type, 1, &dsize,
                                  NULL, 1, -1, NULL, NULL, NULL);
 
   /* Find the index of the quantile. */
@@ -489,12 +494,12 @@ gal_statistics_quantile_function_index(gal_data_t *input, 
gal_data_t *value,
   gal_data_t *nbs=gal_statistics_no_blank_sorted(input, inplace);
 
   /* A sanity check. */
-  if(input->type!=value->type)
+  if(nbs->type!=value->type)
     error(EXIT_FAILURE, 0, "the types of the input dataset and value to "
           "`gal_statistics_quantile_function' have to be the same");
 
   /* Find the result: */
-  switch(input->type)
+  switch(nbs->type)
     {
     case GAL_DATA_TYPE_UINT8:     STATS_QFUNC( uint8_t  );     break;
     case GAL_DATA_TYPE_INT8:      STATS_QFUNC( int8_t   );     break;
@@ -508,7 +513,7 @@ gal_statistics_quantile_function_index(gal_data_t *input, 
gal_data_t *value,
     case GAL_DATA_TYPE_FLOAT64:   STATS_QFUNC( double   );     break;
     default:
       error(EXIT_FAILURE, 0, "type code %d not recognized in "
-            "`gal_statistics_quantile_function'", input->type);
+            "`gal_statistics_quantile_function'", nbs->type);
     }
 
   /* Clean up and return. */
@@ -951,10 +956,11 @@ gal_statistics_mode(gal_data_t *input, float mirrordist, 
int inplace)
   size_t modeindex;
   size_t dsize=4, mdsize=1;
   struct statistics_mode_params p;
-  gal_data_t *mode=gal_data_alloc(NULL, input->type, 1, &mdsize,
-                                 NULL, 1, -1, NULL, NULL, NULL);
-  gal_data_t *b_val=gal_data_alloc(NULL, input->type, 1, &mdsize,
-                                 NULL, 1, -1, NULL, NULL, NULL);
+  int type=gal_tile_block(input)->type;
+  gal_data_t *mode=gal_data_alloc(NULL, type, 1, &mdsize, NULL, 1, -1,
+                                  NULL, NULL, NULL);
+  gal_data_t *b_val=gal_data_alloc(NULL, type, 1, &mdsize, NULL, 1, -1,
+                                   NULL, NULL, NULL);
   gal_data_t *out=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT64, 1, &dsize,
                                  NULL, 1, -1, NULL, NULL, NULL);
 
@@ -1270,29 +1276,49 @@ gal_statistics_sort_decreasing(gal_data_t *data)
    the `inplace' value is set to 1, then the input array will be modified,
    otherwise, a new array will be allocated with the desired properties. So
    if it is already sorted and has blank values, the `inplace' variable is
-   irrelevant.*/
+   irrelevant.
+
+   This function can also work on tiles, in that case, `inplace' is
+   useless, because a tile doesn't own its dataset and the dataset is not
+   contiguous.
+*/
 gal_data_t *
 gal_statistics_no_blank_sorted(gal_data_t *input, int inplace)
 {
-  gal_data_t *noblank, *sorted;
+  gal_data_t *contig, *noblank, *sorted;
+
+  /* If this is a tile, then first we have to copy it into a contiguous
+     piece of memory. After this step, we will only be dealing with
+     `contig' (for a contiguous patch of memory). */
+  if(input->block)
+    {
+      /* Copy the input into a contiguous patch of memory. */
+      contig=gal_tile_to_contiguous(input);
+
+      /* When the data was a tile, we have already copied the array into a
+         separate allocated space. So to avoid any further copying, we will
+         just set the `inplace' variable to 1. */
+      inplace=1;
+    }
+  else contig=input;
 
   /* Make sure there is no blanks in the array that will be used. After
      this step, we won't be dealing with `input' any more, but with
      `noblank'.*/
-  if( gal_blank_present(input) )
+  if( gal_blank_present(contig) )
     {
       if(inplace)                     /* If we can free the input, then */
         {                             /* just remove the blank elements */
-          gal_blank_remove(input);    /* in place. */
-          noblank=input;
+          gal_blank_remove(contig);    /* in place. */
+          noblank=contig;
         }
       else
         {
-          noblank=gal_data_copy(input);   /* We aren't allowed to touch */
-          gal_blank_remove(input);        /* the input, so make a copy. */
+          noblank=gal_data_copy(contig);   /* We aren't allowed to touch */
+          gal_blank_remove(contig);        /* the input, so make a copy. */
         }
     }
-  else noblank=input;
+  else noblank=contig;
 
   /* Make sure the array is sorted. After this step, we won't be dealing
      with `noblank' any more but with `sorted'. */
diff --git a/lib/tile.c b/lib/tile.c
index bab614f..a40ef0b 100644
--- a/lib/tile.c
+++ b/lib/tile.c
@@ -189,6 +189,59 @@ gal_tile_start_end_ind_inclusive(gal_data_t *tile, 
gal_data_t *work,
 
 
 
+/* Return a contiguous patch of memory with the same contents as the
+   tile. If the input is not actually a tile, this function will return the
+   actual input untouched. */
+#define TO_CONTIGUOUS(IT) {                                             \
+    IT *i, *o=out->array, *f, *st;                                      \
+    st=gal_tile_start_end_ind_inclusive(input, block, s_e_ind);         \
+    while( s_e_ind[0] + increment <= s_e_ind[1] )                       \
+      {                                                                 \
+        f = ( i = st + increment ) + input->dsize[input->ndim-1];       \
+        do *o++=*i++; while(i<f);                                       \
+        increment += gal_tile_block_increment(block, input->dsize,      \
+                                              num_increment++, NULL);   \
+      }                                                                 \
+  }
+gal_data_t *
+gal_tile_to_contiguous(gal_data_t *input)
+{
+  gal_data_t *out, *block=gal_tile_block(input);
+  size_t s_e_ind[2], increment=0, num_increment=1;
+
+  /* Check if this is actually a tile. */
+  if(input->block)
+    {
+      /* Allocate the contiguous block of memory. */
+      out=gal_data_alloc(NULL, block->type, input->ndim, input->dsize,
+                         NULL, 0, input->minmapsize, NULL, input->unit,
+                         NULL);
+
+      /* Copy the tile's contents to the contiguous patch of memory. */
+      switch(block->type)
+        {
+        case GAL_DATA_TYPE_UINT8:     TO_CONTIGUOUS( uint8_t  );     break;
+        case GAL_DATA_TYPE_INT8:      TO_CONTIGUOUS( int8_t   );     break;
+        case GAL_DATA_TYPE_UINT16:    TO_CONTIGUOUS( uint16_t );     break;
+        case GAL_DATA_TYPE_INT16:     TO_CONTIGUOUS( int16_t  );     break;
+        case GAL_DATA_TYPE_UINT32:    TO_CONTIGUOUS( uint32_t );     break;
+        case GAL_DATA_TYPE_INT32:     TO_CONTIGUOUS( int32_t  );     break;
+        case GAL_DATA_TYPE_UINT64:    TO_CONTIGUOUS( uint64_t );     break;
+        case GAL_DATA_TYPE_INT64:     TO_CONTIGUOUS( int64_t  );     break;
+        case GAL_DATA_TYPE_FLOAT32:   TO_CONTIGUOUS( float    );     break;
+        case GAL_DATA_TYPE_FLOAT64:   TO_CONTIGUOUS( double   );     break;
+        }
+    }
+  else out=input;
+
+  /* Return. */
+  return out;
+}
+
+
+
+
+
 
 
 
@@ -691,6 +744,7 @@ gal_tile_all_position(gal_data_t *input, size_t *regular,
       tiles[i].array=gal_data_ptr_increment(block->array, tind, block->type);
 
       /* Set the sizes of the tile. */
+      tiles[i].size=1;
       tiles[i].ndim=input->ndim;
       tiles[i].dsize=gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T,input->ndim);
       for(d=0;d<input->ndim;++d)
@@ -707,6 +761,9 @@ gal_tile_all_position(gal_data_t *input, size_t *regular,
             }
           else
             tiles[i].dsize[d]=regular[d];
+
+          /* Set the size value. */
+          tiles[i].size *= tiles[i].dsize[d];
         }
 
       /* Set the block structure for this tile to the `input', and set the