gnuastro-commits
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[gnuastro-commits] master e58ed29 108/125: One value per tile in multi-c


From: Mohammad Akhlaghi
Subject: [gnuastro-commits] master e58ed29 108/125: One value per tile in multi-channel tessellation
Date: Sun, 23 Apr 2017 22:36:49 -0400 (EDT)

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

    One value per tile in multi-channel tessellation
    
    To significantly facilitate the processing, tiles in each channel are
    contiguous. So the tile map would only correspond to the actual input
    dataset when there was only one channel or one dimension. The
    `gal_tile_full_values_for_disp' has been added to do the job: you give it
    an array (of any type) that is assumed to have one value for each tile. You
    also give it the tessellation structure. This function will then sort the
    values in the necessary way so that even when channels are present, the
    tile values will be displayed properly.
    
    This can later be used for interpolation (when interpolation is meant to be
    done within channels or over the whole image for example). Also, it makes
    it very easy to make a FITS file of the tile values for inspection (thanks
    to the new `gal_tile_full_values_write'). To make this second job more
    easily configurable the new `--oneelempertile' option has been added so if
    the user only wants one pixel/element per tile (and not a huge image with
    mostly constant values), they can do so easily with this option.
    
    Other work that has been done in this commit:
    
     - `gal_blank_initialize' is a new function to initialize an array with
       blank values corresponding to its type.
    
     - Until now, `GAL_TILE_PO_OISET' would need the C types as arguments, but
       now, it calculates the types of the input and output to find the proper
       type internally. This greatly simplified the code, since the caller
       doesn't have to worry about the type any more.
    
     - `gal_tile_block_write_const_value' is a new function to write the tile
       values (one value per tile) into an allocated block of memory (with the
       same size as the tile's allocated block of memory). It can write values
       of any type that don't necessarily cover the whole image.
    
     - `gal_tile_block_check_tiles' now uses `gal_tile_block_write_const_value'
       to write the tile IDs into the allocated block.
---
 bin/statistics/statistics.c |   2 +-
 doc/gnuastro.texi           |  39 +++++-
 lib/blank.c                 |  12 ++
 lib/commonopts.h            |  13 ++
 lib/gnuastro/blank.h        |   3 +
 lib/gnuastro/tile.h         | 218 ++++++++++++++++++++++++--------
 lib/options.h               |   1 +
 lib/statistics.c            | 155 ++---------------------
 lib/tile.c                  | 295 ++++++++++++++++++++++++++++++++++----------
 9 files changed, 471 insertions(+), 267 deletions(-)

diff --git a/bin/statistics/statistics.c b/bin/statistics/statistics.c
index f5743c6..f1fa283 100644
--- a/bin/statistics/statistics.c
+++ b/bin/statistics/statistics.c
@@ -379,7 +379,7 @@ statistics_on_tile(struct statisticsparams *p)
         }
 
       /* Save the values. */
-      gal_fits_img_write(values, cp->output, NULL, PROGRAM_STRING);
+      gal_tile_full_values_write(values, tl, cp->output, PROGRAM_STRING);
       gal_data_free(values);
 
       /* Clean up. */
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index fa627a6..0a8c8de 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -4414,11 +4414,12 @@ A FITS binary table (see @ref{Recognized table 
formats}).
 @node Processing options, Operating mode options, Input output options, Common 
options
 @subsubsection Processing options
 
-Some processing steps are common to several programs, so to help in using
-those steps in the different programs, they are defined as common options
-to all programs. Note that this class of common options is thus necessarily
-less common between all the programs than those described in @ref{Input
-output options}, or @ref{Operating mode options} options.
+Some processing steps are common to several programs, so they are defined
+as common options to all programs. Note that this class of common options
+is thus necessarily less common between all the programs than those
+described in @ref{Input output options}, or @ref{Operating mode options}
+options. Also, if they are irrelevant for a program, these options will not
+display in the @option{--help} output of the program.
 
 @table @option
 
@@ -4462,6 +4463,34 @@ Make a FITS file with the same dimensions as the input 
but each pixel is
 replaced with the ID of the tile that it is associated with. Note that the
 tile IDs start from 0. See @ref{Tessellation} for more on Tiling an image
 in Gnuastro.
+
address@hidden --oneelempertile
+When showing the tile values (for example with @option{--checktiles}, or
+when the program's output is tessellated) only use one element for each
+tile. This can be useful when only the relative values given to each tile
+compared to the rest are important or need to be checked. Since the tiles
+usually have a large number of pixels within them the output will be much
+smaller, and so easier to read, write, store, or send.
+
+Note that when the full input size in any dimension is not exactly
+divisible by the given @option{--tilesize} in that dimension, the edge
+tile(s) will have different sizes (in units of the input's size), see
address@hidden But with this option, all displayed values are
+going to have the (same) size of one data-element. Hence, in such cases,
+the image proportions are going to be slightly different with this
+option.
+
+If your input image is not exaclty divisible by the tile size and you want
+one value per tile for some higher-level processing, all is not lost
+though. You can see how many pixels were within each tile (for example to
+weight the values or discard some for later processing) with Gnuastro's
+Statistics (see @ref{Statistics}) as shown below. The output FITS file is
+going to have two extensions, one with the median calculated on each tile
+and one with the number of elements that each tile covers.
+
address@hidden
+$ aststatistics --median --number --ontile input.fits
address@hidden example
 @end table
 
 @node Operating mode options,  , Processing options, Common options
diff --git a/lib/blank.c b/lib/blank.c
index 5281787..776027b 100644
--- a/lib/blank.c
+++ b/lib/blank.c
@@ -78,6 +78,18 @@ gal_blank_write(void *ptr, uint8_t type)
 
 
 
+/* Initialize (set all the values in the array) with the blank value of the
+   given type. */
+void
+gal_blank_initialize(gal_data_t *input)
+{
+  GAL_TILE_PARSE_OPERATE({*i=b;}, input, NULL, 0, 0);
+}
+
+
+
+
+
 /* Allocate some space for the given type and put the blank value into
    it. */
 void *
diff --git a/lib/commonopts.h b/lib/commonopts.h
index 873ed2c..48eb245 100644
--- a/lib/commonopts.h
+++ b/lib/commonopts.h
@@ -159,6 +159,19 @@ struct argp_option gal_commonopts_options[] =
       GAL_OPTIONS_NOT_MANDATORY,
       GAL_OPTIONS_NOT_SET
     },
+    {
+      "oneelempertile",
+      GAL_OPTIONS_KEY_ONEELEMPERTILE,
+      0,
+      0,
+      "Only display 1 element/tile, not full input res.",
+      GAL_OPTIONS_GROUP_TESSELLATION,
+      &cp->tl.oneelempertile,
+      GAL_OPTIONS_NO_ARG_TYPE,
+      GAL_OPTIONS_RANGE_0_OR_1,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
+    },
 
 
 
diff --git a/lib/gnuastro/blank.h b/lib/gnuastro/blank.h
index e35110e..44c424a 100644
--- a/lib/gnuastro/blank.h
+++ b/lib/gnuastro/blank.h
@@ -70,6 +70,9 @@ __BEGIN_C_DECLS  /* From C++ preparations */
 void
 gal_blank_write(void *pointer, uint8_t type);
 
+void
+gal_blank_initialize(gal_data_t *input);
+
 void *
 gal_blank_alloc_write(uint8_t type);
 
diff --git a/lib/gnuastro/tile.h b/lib/gnuastro/tile.h
index f60fe15..eb4782a 100644
--- a/lib/gnuastro/tile.h
+++ b/lib/gnuastro/tile.h
@@ -74,6 +74,10 @@ gal_tile_block_increment(gal_data_t *block, size_t *tsize,
                          size_t num_increment, size_t *coord);
 
 gal_data_t *
+gal_tile_block_write_const_value(gal_data_t *tilevalues, gal_data_t *tilesll,
+                                 int initialize);
+
+gal_data_t *
 gal_tile_block_check_tiles(gal_data_t *tiles);
 
 
@@ -92,8 +96,10 @@ struct gal_tile_two_layer_params
   float          remainderfrac; /* Frac. of remainers in each dim to cut. */
   uint8_t           workoverch; /* Convolve over channel borders.         */
   uint8_t           checktiles; /* Tile IDs in an img, the size of input. */
+  uint8_t       oneelempertile; /* Only use one element for each tile.    */
 
   /* Internal parameters. */
+  size_t                  ndim; /* The number of dimensions.              */
   size_t              tottiles; /* Total number of tiles in all dim.      */
   size_t          tottilesinch; /* Number of tiles in one channel.        */
   size_t           totchannels; /* Total number of channels in all dim.   */
@@ -122,6 +128,17 @@ gal_tile_full_two_layers(gal_data_t *input,
 void
 gal_tile_full_free_contents(struct gal_tile_two_layer_params *tl);
 
+gal_data_t *
+gal_tile_full_values_for_disp(gal_data_t *tilevalues,
+                              struct gal_tile_two_layer_params *tl);
+
+void
+gal_tile_full_values_write(gal_data_t *tilevalues,
+                           struct gal_tile_two_layer_params *tl,
+                           char *filename, char *program_string);
+
+
+
 
 
 /*********************************************************************/
@@ -146,61 +163,14 @@ 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) {        \
-    size_t increment=0, num_increment=1;                                \
-    OT *ost, *o = OUT ? OUT->array : NULL;                              \
-    gal_data_t *iblock=gal_tile_block(IN);                              \
+#define GAL_TILE_PO_OISET(IT, OT, OP, IN, OUT, PARSE_OUT, CHECK_BLANK) { \
+    gal_data_t *out_w=OUT; /* Since `OUT' may be NULL. */               \
+    OT *ost, *o = out_w ? out_w->array : NULL;                          \
     IT b, *st, *i=IN->array, *f=i+IN->size;                             \
-    gal_data_t *oblock = OUT ? gal_tile_block(OUT) : NULL;              \
-    int hasblank=gal_blank_present(IN), parse_out=(OUT && PARSE_OUT);   \
-    size_t s_e_ind[2]={0,iblock->size-1}; /* -1: this is INCLUSIVE */   \
                                                                         \
     /* 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(hasblank) gal_blank_write(&b, iblock->type);                     \
                                                                         \
     /* If this is a tile, not a full block. */                          \
     if(IN!=iblock)                                                      \
@@ -213,6 +183,7 @@ gal_tile_function_on_threads(gal_data_t *tiles, void 
*(*meshfunc)(void *),
     /* Go over contiguous patches of memory. */                         \
     while( s_e_ind[0] + increment <= s_e_ind[1] )                       \
       {                                                                 \
+                                                                        \
         /* If we are on a tile, reset `i' and `f'. 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'.     */  \
@@ -238,8 +209,9 @@ gal_tile_function_on_threads(gal_data_t *tiles, void 
*(*meshfunc)(void *),
           }                                                             \
         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     */ \
+        /* `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   */ \
@@ -256,6 +228,148 @@ gal_tile_function_on_threads(gal_data_t *tiles, void 
*(*meshfunc)(void *),
   }
 
 
+#define GAL_TILE_PO_OSET(OT, OP, IN, OUT, PARSE_OUT, CHECK_BLANK) {     \
+  switch(iblock->type)                                                  \
+    {                                                                   \
+    case GAL_DATA_TYPE_UINT8:                                           \
+      GAL_TILE_PO_OISET(uint8_t,  OT,OP,IN,OUT,PARSE_OUT,CHECK_BLANK);  \
+      break;                                                            \
+    case GAL_DATA_TYPE_INT8:                                            \
+      GAL_TILE_PO_OISET(int8_t,   OT,OP,IN,OUT,PARSE_OUT,CHECK_BLANK);  \
+      break;                                                            \
+    case GAL_DATA_TYPE_UINT16:                                          \
+      GAL_TILE_PO_OISET(uint16_t, OT,OP,IN,OUT,PARSE_OUT,CHECK_BLANK);  \
+      break;                                                            \
+    case GAL_DATA_TYPE_INT16:                                           \
+      GAL_TILE_PO_OISET(int16_t,  OT,OP,IN,OUT,PARSE_OUT,CHECK_BLANK);  \
+      break;                                                            \
+    case GAL_DATA_TYPE_UINT32:                                          \
+      GAL_TILE_PO_OISET(uint32_t, OT,OP,IN,OUT,PARSE_OUT,CHECK_BLANK);  \
+      break;                                                            \
+    case GAL_DATA_TYPE_INT32:                                           \
+      GAL_TILE_PO_OISET(int32_t,  OT,OP,IN,OUT,PARSE_OUT,CHECK_BLANK);  \
+      break;                                                            \
+    case GAL_DATA_TYPE_UINT64:                                          \
+      GAL_TILE_PO_OISET(uint64_t, OT,OP,IN,OUT,PARSE_OUT,CHECK_BLANK);  \
+      break;                                                            \
+    case GAL_DATA_TYPE_INT64:                                           \
+      GAL_TILE_PO_OISET(int64_t,  OT,OP,IN,OUT,PARSE_OUT,CHECK_BLANK);  \
+      break;                                                            \
+    case GAL_DATA_TYPE_FLOAT32:                                         \
+      GAL_TILE_PO_OISET(float,    OT,OP,IN,OUT,PARSE_OUT,CHECK_BLANK);  \
+      break;                                                            \
+    case GAL_DATA_TYPE_FLOAT64:                                         \
+      GAL_TILE_PO_OISET(double,   OT,OP,IN,OUT,PARSE_OUT,CHECK_BLANK);  \
+      break;                                                            \
+    default:                                                            \
+      error(EXIT_FAILURE, 0, "type code %d not recognized in "          \
+            "`GAL_TILE_PO_OSET'", iblock->type);                        \
+    }                                                                   \
+  }
+
+
+/* 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 of the operation into the `gal_data_t'
+   that it points to. Note that OUT 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:
+
+       `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.
+
+       `CHECK_BLANK': If it is non-zero, then the input will be checked for
+             blank values.
+
+   For `OP' you have access to some other variables besides `i' and `o':
+
+      `i': Pointer to this element of input to parse.
+
+      `o': Pointer to corresponding element in output.
+
+      `b': blank value in input't type.
+
+   See `lib/statistics.c' for some example applications of this function.
+*/
+#define GAL_TILE_PARSE_OPERATE(OP, IN, OUT, PARSE_OUT, CHECK_BLANK) {   \
+    int parse_out=(OUT && PARSE_OUT);                                   \
+    size_t increment=0, num_increment=1;                                \
+    gal_data_t *iblock = gal_tile_block(IN);                            \
+    gal_data_t *oblock = OUT ? gal_tile_block(OUT) : NULL;              \
+    int hasblank = CHECK_BLANK ? gal_blank_present(IN) : 0;             \
+    size_t s_e_ind[2]={0,iblock->size-1}; /* -1: this is INCLUSIVE */   \
+                                                                        \
+    /* 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);                                              \
+                                                                        \
+    /* First set the OUTPUT type. */                                    \
+    if(OUT)                                                             \
+      switch(oblock->type)                                              \
+        {                                                               \
+        case GAL_DATA_TYPE_UINT8:                                       \
+          GAL_TILE_PO_OSET(uint8_t,  OP,IN,OUT,PARSE_OUT,CHECK_BLANK);  \
+          break;                                                        \
+        case GAL_DATA_TYPE_INT8:                                        \
+          GAL_TILE_PO_OSET(int8_t,   OP,IN,OUT,PARSE_OUT,CHECK_BLANK);  \
+          break;                                                        \
+        case GAL_DATA_TYPE_UINT16:                                      \
+          GAL_TILE_PO_OSET(uint16_t, OP,IN,OUT,PARSE_OUT,CHECK_BLANK);  \
+          break;                                                        \
+        case GAL_DATA_TYPE_INT16:                                       \
+          GAL_TILE_PO_OSET(int16_t,  OP,IN,OUT,PARSE_OUT,CHECK_BLANK);  \
+          break;                                                        \
+        case GAL_DATA_TYPE_UINT32:                                      \
+          GAL_TILE_PO_OSET(uint32_t, OP,IN,OUT,PARSE_OUT,CHECK_BLANK);  \
+          break;                                                        \
+        case GAL_DATA_TYPE_INT32:                                       \
+          GAL_TILE_PO_OSET(int32_t,  OP,IN,OUT,PARSE_OUT,CHECK_BLANK);  \
+          break;                                                        \
+        case GAL_DATA_TYPE_UINT64:                                      \
+          GAL_TILE_PO_OSET(uint64_t, OP,IN,OUT,PARSE_OUT,CHECK_BLANK);  \
+          break;                                                        \
+        case GAL_DATA_TYPE_INT64:                                       \
+          GAL_TILE_PO_OSET(int64_t,  OP,IN,OUT,PARSE_OUT,CHECK_BLANK);  \
+          break;                                                        \
+        case GAL_DATA_TYPE_FLOAT32:                                     \
+          GAL_TILE_PO_OSET(float,    OP,IN,OUT,PARSE_OUT,CHECK_BLANK);  \
+          break;                                                        \
+        case GAL_DATA_TYPE_FLOAT64:                                     \
+          GAL_TILE_PO_OSET(double,   OP,IN,OUT,PARSE_OUT,CHECK_BLANK);  \
+          break;                                                        \
+        default:                                                        \
+          error(EXIT_FAILURE, 0, "type code %d not recognized in "      \
+                "`GAL_TILE_PARSE_OPERATE'", oblock->type);              \
+        }                                                               \
+    else                                                                \
+      /* When `OUT==NULL', its type is irrelevant, we'll just use */    \
+      /*`int' as a place holder. */                                     \
+      GAL_TILE_PO_OSET(int,          OP,IN,OUT,PARSE_OUT,CHECK_BLANK);  \
+  }
 
 
 
diff --git a/lib/options.h b/lib/options.h
index fe76904..293bf6c 100644
--- a/lib/options.h
+++ b/lib/options.h
@@ -105,6 +105,7 @@ enum options_common_keys
   GAL_OPTIONS_KEY_ONLYVERSION,
   GAL_OPTIONS_KEY_WORKOVERCH,
   GAL_OPTIONS_KEY_CHECKTILES,
+  GAL_OPTIONS_KEY_ONEELEMPERTILE,
 };
 
 
diff --git a/lib/statistics.c b/lib/statistics.c
index db4a3af..71b6843 100644
--- a/lib/statistics.c
+++ b/lib/statistics.c
@@ -54,37 +54,17 @@ 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) {                                                 \
-    GAL_TILE_PARSE_OPERATE(IT, uint64_t, {++(*o);}, input, out, 0)      \
-  }
-
 gal_data_t *
 gal_statistics_number(gal_data_t *input)
 {
   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);
 
   /* 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);
-      }
+    { GAL_TILE_PARSE_OPERATE({++(*o);}, input, out, 0, 1) }
   else
     *((uint64_t *)(out->array)) = input->size;
 
@@ -99,39 +79,18 @@ gal_statistics_number(gal_data_t *input)
    the dataset. Note that a NaN (blank in floating point) will fail on any
    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) {                                                 \
-    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;
-  int type=gal_tile_block(input)->type;
-  gal_data_t *out=gal_data_alloc(NULL, type, 1, &dsize, NULL, 1, -1,
-                                 NULL, NULL, NULL);
+  gal_data_t *out=gal_data_alloc(NULL, gal_tile_block(input)->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);
 
   /* 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;
-    case GAL_DATA_TYPE_UINT16:    STATS_MIN( uint16_t );    break;
-    case GAL_DATA_TYPE_INT16:     STATS_MIN( int16_t  );    break;
-    case GAL_DATA_TYPE_UINT32:    STATS_MIN( uint32_t );    break;
-    case GAL_DATA_TYPE_INT32:     STATS_MIN( int32_t  );    break;
-    case GAL_DATA_TYPE_UINT64:    STATS_MIN( uint64_t );    break;
-    case GAL_DATA_TYPE_INT64:     STATS_MIN( int64_t  );    break;
-    case GAL_DATA_TYPE_FLOAT32:   STATS_MIN( float    );    break;
-    case GAL_DATA_TYPE_FLOAT64:   STATS_MIN( double   );    break;
-    default:
-      error(EXIT_FAILURE, 0, "type code %d not recognized in "
-            "`gal_statistics_minimum'", type);
-    }
+  GAL_TILE_PARSE_OPERATE({*o = *i < *o ? *i : *o; ++n;}, input, out, 0, 0);
 
   /* If there were no usable elements, set the output to blank, then
      return. */
@@ -145,39 +104,18 @@ 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) {                                                 \
-    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;
-  int type=gal_tile_block(input)->type;
-  gal_data_t *out=gal_data_alloc(NULL, type, 1, &dsize,
-                                 NULL, 1, -1, NULL, NULL, NULL);
+  gal_data_t *out=gal_data_alloc(NULL, gal_tile_block(input)->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);
 
   /* 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;
-    case GAL_DATA_TYPE_UINT16:    STATS_MAX( uint16_t );    break;
-    case GAL_DATA_TYPE_INT16:     STATS_MAX( int16_t  );    break;
-    case GAL_DATA_TYPE_UINT32:    STATS_MAX( uint32_t );    break;
-    case GAL_DATA_TYPE_INT32:     STATS_MAX( int32_t  );    break;
-    case GAL_DATA_TYPE_UINT64:    STATS_MAX( uint64_t );    break;
-    case GAL_DATA_TYPE_INT64:     STATS_MAX( int64_t  );    break;
-    case GAL_DATA_TYPE_FLOAT32:   STATS_MAX( float    );    break;
-    case GAL_DATA_TYPE_FLOAT64:   STATS_MAX( double   );    break;
-    default:
-      error(EXIT_FAILURE, 0, "type code %d not recognized in "
-            "`gal_statistics_maximum'", type);
-    }
+  GAL_TILE_PARSE_OPERATE({*o = *i > *o ? *i : *o; ++n;}, input, out, 0, 0);
 
   /* If there were no usable elements, set the output to blank, then
      return. */
@@ -191,35 +129,16 @@ 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) {                                             \
-    GAL_TILE_PARSE_OPERATE(IT, double, {++n; *o += *i;}, input, out, 0) \
-  }
 gal_data_t *
 gal_statistics_sum(gal_data_t *input)
 {
   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);
 
   /* 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;
-    case GAL_DATA_TYPE_UINT16:    STATS_SUM_NUM( uint16_t );    break;
-    case GAL_DATA_TYPE_INT16:     STATS_SUM_NUM( int16_t  );    break;
-    case GAL_DATA_TYPE_UINT32:    STATS_SUM_NUM( uint32_t );    break;
-    case GAL_DATA_TYPE_INT32:     STATS_SUM_NUM( int32_t  );    break;
-    case GAL_DATA_TYPE_UINT64:    STATS_SUM_NUM( uint64_t );    break;
-    case GAL_DATA_TYPE_INT64:     STATS_SUM_NUM( int64_t  );    break;
-    case GAL_DATA_TYPE_FLOAT32:   STATS_SUM_NUM( float    );    break;
-    case GAL_DATA_TYPE_FLOAT64:   STATS_SUM_NUM( double   );    break;
-    default:
-      error(EXIT_FAILURE, 0, "type code %d not recognized in "
-            "`gal_statistics_sum'", type);
-    }
+  GAL_TILE_PARSE_OPERATE({++n; *o += *i;}, input, out, 0, 1);
 
   /* If there were no usable elements, set the output to blank, then
      return. */
@@ -237,28 +156,12 @@ gal_data_t *
 gal_statistics_mean(gal_data_t *input)
 {
   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);
 
   /* 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;
-    case GAL_DATA_TYPE_UINT16:    STATS_SUM_NUM( uint16_t );    break;
-    case GAL_DATA_TYPE_INT16:     STATS_SUM_NUM( int16_t  );    break;
-    case GAL_DATA_TYPE_UINT32:    STATS_SUM_NUM( uint32_t );    break;
-    case GAL_DATA_TYPE_INT32:     STATS_SUM_NUM( int32_t  );    break;
-    case GAL_DATA_TYPE_UINT64:    STATS_SUM_NUM( uint64_t );    break;
-    case GAL_DATA_TYPE_INT64:     STATS_SUM_NUM( int64_t  );    break;
-    case GAL_DATA_TYPE_FLOAT32:   STATS_SUM_NUM( float    );    break;
-    case GAL_DATA_TYPE_FLOAT64:   STATS_SUM_NUM( double   );    break;
-    default:
-      error(EXIT_FAILURE, 0, "type code %d not recognized in "
-            "`gal_statistics_mean'", type);
-    }
+  GAL_TILE_PARSE_OPERATE({++n; *o += *i;}, input, out, 0, 1);
 
   /* 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
@@ -274,36 +177,16 @@ 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) {                                                \
-    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);
 
   /* 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;
-    case GAL_DATA_TYPE_UINT16:    STATS_NSS2( uint16_t );    break;
-    case GAL_DATA_TYPE_INT16:     STATS_NSS2( int16_t  );    break;
-    case GAL_DATA_TYPE_UINT32:    STATS_NSS2( uint32_t );    break;
-    case GAL_DATA_TYPE_INT32:     STATS_NSS2( int32_t  );    break;
-    case GAL_DATA_TYPE_UINT64:    STATS_NSS2( uint64_t );    break;
-    case GAL_DATA_TYPE_INT64:     STATS_NSS2( int64_t  );    break;
-    case GAL_DATA_TYPE_FLOAT32:   STATS_NSS2( float    );    break;
-    case GAL_DATA_TYPE_FLOAT64:   STATS_NSS2( double   );    break;
-    default:
-      error(EXIT_FAILURE, 0, "type code %d not recognized in "
-            "`gal_statistics_maximum'", type);
-    }
+  GAL_TILE_PARSE_OPERATE({++n; s += *i; s2 += *i * *i;}, input, out, 0, 1);
 
   /* Write the standard deviation into the output. */
   *((double *)(out->array)) = n ? sqrt( (s2-s*s/n)/n ) : GAL_BLANK_FLOAT64;
@@ -322,27 +205,11 @@ 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);
 
   /* 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;
-    case GAL_DATA_TYPE_UINT16:    STATS_NSS2( uint16_t );    break;
-    case GAL_DATA_TYPE_INT16:     STATS_NSS2( int16_t  );    break;
-    case GAL_DATA_TYPE_UINT32:    STATS_NSS2( uint32_t );    break;
-    case GAL_DATA_TYPE_INT32:     STATS_NSS2( int32_t  );    break;
-    case GAL_DATA_TYPE_UINT64:    STATS_NSS2( uint64_t );    break;
-    case GAL_DATA_TYPE_INT64:     STATS_NSS2( int64_t  );    break;
-    case GAL_DATA_TYPE_FLOAT32:   STATS_NSS2( float    );    break;
-    case GAL_DATA_TYPE_FLOAT64:   STATS_NSS2( double   );    break;
-    default:
-      error(EXIT_FAILURE, 0, "type code %d not recognized in "
-            "`gal_statistics_maximum'", type);
-    }
+  GAL_TILE_PARSE_OPERATE({++n; s += *i; s2 += *i * *i;}, input, out, 0, 1);
 
   /* Write in the output values and return. */
   ((double *)(out->array))[0] = n ? s/n                  : GAL_BLANK_FLOAT64;
diff --git a/lib/tile.c b/lib/tile.c
index 789e19f..0b28f8a 100644
--- a/lib/tile.c
+++ b/lib/tile.c
@@ -311,6 +311,71 @@ gal_tile_block_increment(gal_data_t *block, size_t *tsize,
 
 
 
+/* Write a constant value for each tile into each pixel covered by the
+   input tiles in an array the size of the block and return it.
+
+   Arguments
+   ---------
+
+     `tilevalues': This must be an array that has the same number of
+        elements as that in `tilesll' and in the same order that `tilesll'
+        elements are parsed (from first to last). As a result the
+        dimensionality of this array is irrelevant. Note that unlike
+        `tiles', `tilevalues' must be an array.
+
+     `tiles': This will be parsed as a linked list (using the `next'
+        element). Internally, it might be stored as an array, but this
+        function doesn't care! The position of the tile over its block will
+        be determined according to the `block' element and the pointer of
+        its `array' as fully described in `gnuastro/data.h'. This function
+        will not pop/free the list, it will only parse it from start to
+        end.
+
+     `initialize': Initialize the allocated space with blank values before
+        writing in the constant values. This can be useful when the tiles
+        don't cover the full allocated block. */
+gal_data_t *
+gal_tile_block_write_const_value(gal_data_t *tilevalues, gal_data_t *tilesll,
+                                 int initialize)
+{
+  void *in;
+  int type=tilevalues->type;
+  size_t tile_ind, nt=0, nv=tilevalues->size;
+  gal_data_t *tofill, *tile, *block=gal_tile_block(tilesll);
+
+  /* A small sanity check. */
+  for(tile=tilesll; tile!=NULL; tile=tile->next) ++nt;
+  if(nt!=nv)
+    error(EXIT_FAILURE, 0, "the number of elements in `tilevalues' (%zu) "
+          "and `tilesll' (%zu) must be the same in "
+          "`gal_tile_block_write_const_value'", nv, nt);
+
+  /* Allocate the output array. */
+  tofill=gal_data_alloc(NULL, type, block->ndim, block->dsize, block->wcs,
+                        0, block->minmapsize, tilevalues->name,
+                        tilevalues->unit, tilevalues->comment);
+
+  /* If requested, initialize `tofill'. */
+  if(initialize) gal_blank_initialize(tofill);
+
+  /* Go over the tiles and write the values in. Note Recall that `tofill'
+     has the same type as `tilevalues'. So we are using memcopy. */
+  tile_ind=0;
+  for(tile=tilesll; tile!=NULL; tile=tile->next)
+    {
+      /* Set the pointer to use as input. */
+      in=gal_data_ptr_increment(tilevalues->array, tile_ind++, type);;
+      GAL_TILE_PARSE_OPERATE({memcpy(o, in, gal_data_sizeof(type));},
+                             tile, tofill, 1, 0);
+    }
+
+  return tofill;
+}
+
+
+
+
+
 /* Make a copy of the memory block in integer type and fill it with the ID
    of each tile, the non-filled areas have blank values. Finally, save the
    final array into a FITS file, specified with `filename'. This is done
@@ -318,70 +383,25 @@ gal_tile_block_increment(gal_data_t *block, size_t *tsize,
    32-bit type because this is the standard FITS standard type for
    integers. */
 gal_data_t *
-gal_tile_block_check_tiles(gal_data_t *tiles)
+gal_tile_block_check_tiles(gal_data_t *tilesll)
 {
-  size_t num_increment;
-  gal_data_t *tofill, *tile;
-  gal_data_t *block=gal_tile_block(tiles);
-  int32_t *p, *pf, tile_index=0, *start=NULL;
-  size_t ndim=tiles->ndim, increment, start_end_inc[2];
-  size_t *coord=gal_data_calloc_array(GAL_DATA_TYPE_SIZE_T, ndim);
-
-  /***************************************************************/
-  /*************            For a check           ****************
-  float c=0;
-  block->wcs=NULL;
-  ndim=block->ndim=tiles->ndim=3;
-  block->dsize=gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, ndim);
-  block->dsize[0]=5; block->dsize[1]=5; block->dsize[2]=5;
-
-  tiles->dsize=gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, ndim);
-  tiles->dsize[0]=2; tiles->dsize[1]=3; tiles->dsize[2]=3;
-  tiles->array=gal_data_ptr_increment(block->array, 36, block->type);
-  tiles->next=NULL;
-  **************************************************************/
-
-  /* Allocate the output array. */
-  tofill=gal_data_alloc(NULL, GAL_DATA_TYPE_INT32, ndim, block->dsize,
-                        block->wcs, 0, block->minmapsize, "TILE_CHECK",
-                        "counts", "indexs of all tiles");
+  int32_t *arr;
+  size_t i, dsize=gal_data_num_in_ll(tilesll);
+  gal_data_t *ids, *out, *block=gal_tile_block(tilesll);
 
-  /* Initialize the allocated space with blank characters for this type. */
-  pf=(p=tofill->array)+tofill->size; do *p++=GAL_BLANK_INT32; while(p<pf);
+  /* Allocate the array to keep the IDs of each tile. */
+  ids=gal_data_alloc(NULL, GAL_DATA_TYPE_INT32, 1, &dsize,
+                     NULL, 0, block->minmapsize, NULL, NULL, NULL);
 
-  /* Fill in the labels of each tile. */
-  for(tile=tiles; tile!=NULL; tile=tile->next)
-    {
-      /* Set the starting and ending indexs of this tile over the allocated
-         block. */
-      start=gal_tile_start_end_ind_inclusive(tile, tofill, start_end_inc);
-
-      /* Go over the full area of this tile. The loop will stop as soon as
-         the incrementation will go over the last index of the tile. Note
-         that num_increment has to start from 1 because having a remainder
-         of zero is meaningful in the calculation of the increment. */
-      increment=0;
-      num_increment=1;
-      while( start_end_inc[0] + increment <= start_end_inc[1] )
-        {
-          /* Parse the elements in the fastest-dimension (the contiguous
-             patch of memory associated with this tile). */
-          pf = ( p = start + increment ) + tile->dsize[ndim-1];
-          do *p++=tile_index; while(p<pf);
-
-          /* Increase the increment from the start of the tile for the next
-             contiguous patch. */
-          increment += gal_tile_block_increment(block, tile->dsize,
-                                                num_increment++, coord);
-        }
+  /* Put the IDs into the array. */
+  arr=ids->array; for(i=0;i<dsize;++i) arr[i]=i;
 
-      /* Increment the index for the next tile. */
-      ++tile_index;
-    }
+  /* Make the output. */
+  out=gal_tile_block_write_const_value(ids, tilesll, 1);
 
-  /* Clean up. */
-  free(coord);
-  return tofill;
+  /* Clean up and return. */
+  gal_data_free(ids);
+  return out;
 }
 
 
@@ -757,7 +777,7 @@ gal_tile_full_two_layers(gal_data_t *input,
                          struct gal_tile_two_layer_params *tl)
 {
   gal_data_t *t;
-  size_t i, *junk;
+  size_t i, *junk, ndim=tl->ndim=input->ndim;
 
   /* Initialize. Note that `numchannels might have already been
      allocated. */
@@ -767,7 +787,7 @@ gal_tile_full_two_layers(gal_data_t *input,
   /* Initialize necessary values and do the channels tessellation. */
   tl->numchannels = gal_tile_full(input, tl->channelsize, tl->remainderfrac,
                                 &tl->channels, 1);
-  tl->totchannels = gal_multidim_total_size(input->ndim, tl->numchannels);
+  tl->totchannels = gal_multidim_total_size(ndim, tl->numchannels);
 
 
   /* Tile each channel. While tiling the first channel, we are also going
@@ -776,7 +796,7 @@ gal_tile_full_two_layers(gal_data_t *input,
   tl->numtilesinch = gal_tile_full(tl->channels, tl->tilesize,
                                    tl->remainderfrac, &tl->tiles,
                                    tl->totchannels);
-  tl->tottilesinch = gal_multidim_total_size(input->ndim, tl->numtilesinch);
+  tl->tottilesinch = gal_multidim_total_size(ndim, tl->numtilesinch);
   for(i=1; i<tl->totchannels; ++i)
     {
       /* Set the first tile in this channel. Then use it it fill the `next'
@@ -795,10 +815,155 @@ gal_tile_full_two_layers(gal_data_t *input,
   /* Multiply the number of tiles along each dimension OF ONE CHANNEL by
      the number of channels in each dimension to get the dimensionality of
      the full tile structure. */
-  tl->numtiles = gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, input->ndim);
-  for(i=0;i<input->ndim;++i)
+  tl->numtiles = gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, ndim);
+  for(i=0;i<ndim;++i)
     tl->numtiles[i] = tl->numtilesinch[i] * tl->numchannels[i];
-  tl->tottiles = gal_multidim_total_size(input->ndim, tl->numtiles);
+  tl->tottiles = gal_multidim_total_size(ndim, tl->numtiles);
+}
+
+
+
+
+
+/* If necessary (see below), make a new Re-ordered array that has one
+   pixel/element for each tile, such that it is ready for displaying. When
+   it isn't necessary, this function will just return the input data
+   structure. So you should free the output only if it is different from
+   the input.
+
+   When there is only one channel OR one dimension, you can just save the
+   array and everything will be fine. However, when there is more than one
+   channel AND more than one dimension this won't work and the values will
+   not be in the places corresponding to the original dataset. This happens
+   because the tiles (and thus the value you assign to each tile) in each
+   channel are taken to be a contiguous patch of memory. But when you look
+   at the whole tessellation, you want the tiles to be ordered based on
+   their position in the final dataset/image.
+
+   The example below may help clarify: assume you have a 6x6 tessellation
+   with two channels in the horizontal and one in the vertical. On the left
+   you can see how the tiles (and their values) are allocated in memory. On
+   the right is how they will be displayed if you just save it as a FITS
+   file with no modification (this will not correspond to your input
+   dataset).
+
+              Allocation map                  When displayed
+              --------------                  --------------
+              15 16 17 33 34 35               30 31 32 33 34 35
+              12 13 14 30 31 32               24 25 26 27 28 29
+              09 10 11 27 28 29               18 19 20 21 22 23
+              06 07 08 24 25 26               12 13 14 15 16 17
+              03 04 05 21 22 23               06 07 08 09 10 11
+              00 01 02 18 19 20               00 01 02 03 04 05       */
+gal_data_t *
+gal_tile_full_values_for_disp(gal_data_t *tilevalues,
+                              struct gal_tile_two_layer_params *tl)
+{
+  void *start;
+  gal_data_t *out, *fake;
+  size_t o_inc, num_increment, start_ind;
+  size_t i, ch_ind, contig, i_inc=0, ndim=tl->ndim;
+  size_t *chstart=gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, ndim);
+  size_t *s_e_inc=gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, ndim);
+
+  /* Make sure that the input and tessellation correspond to each
+     other. This is important, because the user might mistakenly give the
+     input dataset as input, not the dataset that has one value per tile.*/
+  if(tilevalues->ndim!=tl->ndim)
+    error(EXIT_FAILURE, 0, "the input (`tilevalues') and tessellation "
+          "(`tl') in `gal_tile_full_values_for_disp' have %zu and %zu "
+          "dimensions respectively! They must have the same number of "
+          "dimensions", tilevalues->ndim, tl->ndim);
+  for(i=0; i<ndim; ++i)
+    if(tilevalues->dsize[i]!=tl->numtiles[i])
+      error(EXIT_FAILURE, 0, "the total number of tiles (in all channels) "
+            "of dimension %zu does not correspond between the input "
+            "(`tilevalues': %zu) and tessellation (`tl': %zu)", ndim-i,
+            tilevalues->dsize[i], tl->numtiles[i]);
+
+  /* If there is only one dimension or one channel, then just return the
+     input `tilevalues'. */
+  if(tilevalues->ndim==1 || tl->totchannels==1)
+    return tilevalues;
+
+  /* Allocate the output. */
+  out=gal_data_alloc(NULL, tilevalues->type, tilevalues->ndim,
+                     tilevalues->dsize, tilevalues->wcs, 0,
+                     tilevalues->minmapsize, tilevalues->name,
+                     tilevalues->unit, tilevalues->comment);
+
+  /* Allocate a fake tile (with a size equal to the number of tiles in one
+     channel) for easy parsing with standard techniques. */
+  fake=gal_data_alloc(NULL, GAL_DATA_TYPE_UINT8, tilevalues->ndim,
+                      tl->numtilesinch, NULL, 0, tilevalues->minmapsize,
+                      NULL, NULL, NULL);
+
+  /* Initialize the fake array: we will set it's `block' to the output.*/
+  fake->block=out;
+  free(fake->array);
+  fake->type=GAL_DATA_TYPE_INVALID;
+
+  /* Put the values into the proper place. */
+  contig=tl->numtilesinch[ndim-1];
+  for(ch_ind=0; ch_ind<tl->totchannels; ++ch_ind)
+    {
+      /* Set the starting pointer of this channel for the `fake' tile that
+         represents it. */
+      gal_multidim_index_to_coord(ch_ind, ndim, tl->numchannels, chstart);
+      for(i=0;i<ndim;++i) chstart[i]*=tl->numtilesinch[i];
+      start_ind=gal_multidim_coord_to_index(ndim, tl->numtiles, chstart);
+      fake->array=gal_data_ptr_increment(out->array, start_ind, out->type);
+
+      /* Find the starting and ending tile indexs (inclusive in the output)
+         of this channel. */
+      start=gal_tile_start_end_ind_inclusive(fake, out, s_e_inc);
+
+      /* Start parsing the channel and putting it in the output. */
+      o_inc=0;
+      num_increment=1;
+      while(s_e_inc[0] + o_inc <= s_e_inc[1])
+        {
+          /* Copy the contiguous region from the `tilevalues' array to the
+             output array. Note that since they have the same type, we can
+             just use `memcpy' and don't actually have to parse it.*/
+          memcpy(gal_data_ptr_increment(start,             o_inc, out->type),
+                 gal_data_ptr_increment(tilevalues->array, i_inc, out->type),
+                 gal_data_sizeof(out->type)*contig);
+
+          /* Set the incrementation for the next contiguous patch. */
+          o_inc += gal_tile_block_increment(out, fake->dsize,
+                                            num_increment++, NULL);
+          i_inc += contig;
+        }
+    }
+
+  /* Clean up and return. */
+  fake->array=NULL;
+  gal_data_free(fake);
+  free(chstart);
+  free(s_e_inc);
+  return out;
+}
+
+
+
+
+
+void
+gal_tile_full_values_write(gal_data_t *tilevalues,
+                           struct gal_tile_two_layer_params *tl,
+                           char *filename, char *program_string)
+{
+  gal_data_t *disp;
+
+  /* Make the dataset to be displayed. */
+  disp = ( tl->oneelempertile
+           ? gal_tile_full_values_for_disp(tilevalues, tl)
+           : gal_tile_block_write_const_value(tilevalues, tl->tiles, 0) );
+
+  /* Write the array as a file and then clean up. */
+  gal_fits_img_write(disp, filename, NULL, program_string);
+  gal_data_free(disp);
 }
 
 



reply via email to

[Prev in Thread] Current Thread [Next in Thread]