gnuastro-commits
[Top][All Lists]
Advanced

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

[gnuastro-commits] master 486f915 067/113: Imported recent work in maste


From: Mohammad Akhlaghi
Subject: [gnuastro-commits] master 486f915 067/113: Imported recent work in master, conflicts fixed
Date: Fri, 16 Apr 2021 10:33:48 -0400 (EDT)

branch: master
commit 486f9151348e34e5aae49698e4ede0cd280b0ddb
Merge: f537368 8dab7b8
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    Imported recent work in master, conflicts fixed
    
    Several conflicts came up regarding the new allocations functions and
    `bin/noisechisel/ui.c', all were fixed.
---
 NEWS                         |  24 +-
 THANKS                       |   1 +
 bin/arithmetic/arithmetic.c  |   7 +-
 bin/arithmetic/operands.c    |   2 +-
 bin/buildprog/buildprog.c    |   3 +-
 bin/convolve/convolve.c      |  13 +-
 bin/crop/onecrop.c           |   5 +-
 bin/crop/ui.c                |   7 +-
 bin/crop/wcsmode.c           |   5 +-
 bin/fits/fits.c              |   4 +-
 bin/match/match.c            |  10 +-
 bin/mkcatalog/columns.c      |  10 +-
 bin/mkcatalog/mkcatalog.c    |  17 +-
 bin/mkcatalog/parse.c        |  26 +-
 bin/mkcatalog/ui.c           |  31 ++-
 bin/mkcatalog/upperlimit.c   |  21 +-
 bin/mkprof/mkprof.c          |  11 +-
 bin/mkprof/oneprofile.c      |   7 +-
 bin/mkprof/ui.c              |  56 ++--
 bin/noisechisel/detection.c  |  36 +--
 bin/noisechisel/sky.c        |  52 ++--
 bin/noisechisel/threshold.c  |  29 +-
 bin/noisechisel/ui.c         |  49 +++-
 bin/segment/clumps.c         |  44 +--
 bin/segment/segment.c        |  11 +-
 bin/segment/ui.c             |   8 +-
 bin/statistics/sky.c         |  13 +-
 bin/statistics/statistics.c  |   5 +-
 bin/warp/ui.c                |  11 +-
 doc/announce-acknowledge.txt |   1 +
 doc/gnuastro.texi            | 643 ++++++++++++++++++++++++++-----------------
 lib/Makefile.am              |  13 +-
 lib/arithmetic.c             |  36 +--
 lib/binary.c                 |  17 +-
 lib/blank.c                  |   3 +-
 lib/convolve.c               |  32 +--
 lib/data.c                   | 212 +-------------
 lib/dimension.c              |  25 +-
 lib/fits.c                   |  30 +-
 lib/gnuastro/data.h          |  19 +-
 lib/gnuastro/dimension.h     |   3 +
 lib/gnuastro/pointer.h       |  70 +++++
 lib/gnuastro/qsort.h         |  79 ++++--
 lib/gnuastro/tile.h          |   5 +-
 lib/interpolate.c            |  45 +--
 lib/label.c                  |  23 +-
 lib/list.c                   |   9 +-
 lib/match.c                  |   9 +-
 lib/options.c                |   4 +-
 lib/permutation.c            |   7 +-
 lib/pointer.c                | 178 ++++++++++++
 lib/qsort.c                  | 140 +++++++---
 lib/statistics.c             |  47 ++--
 lib/tableintern.c            |   5 +-
 lib/tiff.c                   |   5 +-
 lib/tile.c                   |  92 ++++---
 lib/txt.c                    |   2 +-
 lib/type.c                   |   8 +-
 lib/wcs.c                    |  35 ++-
 59 files changed, 1347 insertions(+), 968 deletions(-)

diff --git a/NEWS b/NEWS
index ba39e32..a870f12 100644
--- a/NEWS
+++ b/NEWS
@@ -90,7 +90,11 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
     gal_pdf_name_is_pdf: Returns 1 if given filename is PDF.
     gal_pdf_suffix_is_pdf: Returns 1 if given suffix is PDF.
     gal_pdf_write: Writes a dataset into an PDF file.
-    gal_qsort_index_float_increasing: Sort indexs in increasing value order.
+    gal_pointer_allocate_mmap: Allocate space in a file, not in RAM.
+    gal_qsort_index_single_d: Sort indexs of single array in decreasing order.
+    gal_qsort_index_single_i: Sort indexs of single array in decreasing order.
+    gal_qsort_index_multi_d: Sort indexs in multiple arrays (on different 
threads).
+    gal_qsort_index_multi_i: Sort indexs in multiple arrays (on different 
threads).
     gal_tiff_name_is_tiff: check if name contains a TIFF suffix.
     gal_tiff_suffix_is_tiff: check if suffix is a TIFF suffix.
     gal_tiff_dir_string_read: convert a string to a TIFF directory number.
@@ -143,7 +147,6 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
     - With no output name, the output has a `_detected.fits' suffix.
 
   MakeCatalog:
-
     - The `WCLUMPS' keyword in the objects labeled image is no longer used
          to see if a clumps catalog should also be made. To build a clumps
          catalog, you can now use the `--clumpscat' option.
@@ -158,26 +161,28 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
          still acceptable also).
 
   Libraries:
+    gal_binary_holes_fill: new name for `gal_binary_fill_holes'.
+    gal_dimension_is_different: new name for `gal_data_dsize_is_different'.
     gal_fits_img_read: now only reads the data not the WCS, therefore it no
         longer needs the last two arguments. A subsequent call to
         `gal_wcs_read' can be used to read the WCS information in the file.
-
+    gal_pointer_increment: new name for `gal_data_ptr_increment'.
+    gal_pointer_num_between: new name for `gal_data_ptr_dist'.
+    gal_pointer_allocate: replaces `gal_data_malloc_array' and
+        `gal_data_calloc_array', through an argument you can ask for the
+        allocated memory to be cleared or not.
+    gal_qsort_TYPE_i: new name for gal_qsort_TYPE_increasing.
+    gal_qsort_TYPE_d: new name for gal_qsort_TYPE_decreasing.
     gal_statistics_is_sorted: can now also update the bit-flags regarding
         the sorted nature of the input (to optimize future calls to the
         function).
-
     gal_statistics_quantile_function: returns `inf' or `-inf' if the given
         value is smaller than the minimum or larger than the maximum of the
         input dataset's range. Until now, it would return blank in such
         cases.
-
     gal_statistics_number: the output dataset now has a `size_t' type. Until
         now it was `uint64_t'.
 
-  Library:
-    - gal_binary_holes_fill: new name for `gal_binary_fill_holes'.
-    - gal_data_num_between: new name for `gal_data_ptr_dist'.
-
 ** Bug fixes
 
   bug #52979: Many unused result warnings for asprintf in some compilers.
@@ -193,6 +198,7 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
   bug #53407: Instrumental noise in MakeNoise should be squared.
   bug #53424: Sigma-clipping seg-faults when there are no more elements.
   bug #53580: Warp crash when no WCS present.
+  bug #53825: NoiseChisel not accounting for fully zero-valued tiles.
 
 
 
diff --git a/THANKS b/THANKS
index 5103035..22fafa8 100644
--- a/THANKS
+++ b/THANKS
@@ -34,6 +34,7 @@ support in Gnuastro. The list is ordered alphabetically (by 
family name).
     Brandon Invergo                      brandon@gnu.org
     Aurélien Jarno                       aurelien.jarno@univ-lyon1.fr
     Lee Kelvin                           l.s.kelvin@ljmu.ac.uk
+    Brandon Kelly                        b.k.kelly@2017.ljmu.ac.uk
     Mohammad-Reza Khellat                moha.khe@gmail.com
     Floriane Leclercq                    floriane.leclercq@univ-lyon1.fr
     Alan Lefor                           alefor@astr.tohoku.ac.jp
diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index ac1a183..44ab2e5 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -34,6 +34,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/array.h>
 #include <gnuastro/binary.h>
 #include <gnuastro/threads.h>
+#include <gnuastro/pointer.h>
 #include <gnuastro/dimension.h>
 #include <gnuastro/statistics.h>
 #include <gnuastro/arithmetic.h>
@@ -219,8 +220,7 @@ arithmetic_filter(void *in_prm)
 
       /* Set the tile's starting pointer. */
       index=gal_dimension_coord_to_index(ndim, dsize, start);
-      tile->array=gal_data_ptr_increment(input->array, index,
-                                         input->type);
+      tile->array=gal_pointer_increment(input->array, index, input->type);
 
       /* Do the necessary calculation. */
       switch(afp->operator)
@@ -278,8 +278,7 @@ arithmetic_filter(void *in_prm)
 
 
       /* Copy the result into the output array. */
-      memcpy(gal_data_ptr_increment(afp->out->array, ind,
-                                    afp->out->type),
+      memcpy(gal_pointer_increment(afp->out->array, ind, afp->out->type),
              result->array, gal_type_sizeof(afp->out->type));
 
       /* Clean up for this pixel. */
diff --git a/bin/arithmetic/operands.c b/bin/arithmetic/operands.c
index 48d3632..8329b5e 100644
--- a/bin/arithmetic/operands.c
+++ b/bin/arithmetic/operands.c
@@ -138,7 +138,7 @@ operands_pop(struct arithmeticparams *p, char *operator)
          checks. */
       if(p->refdata.ndim)
         {
-          if(gal_data_dsize_is_different(&p->refdata, data))
+          if( gal_dimension_is_different(&p->refdata, data) )
             error(EXIT_FAILURE, 0, "%s (hdu=%s): has a different size "
                   "compared to previous images. All the images must be "
                   "the same size in order for Arithmetic to work",
diff --git a/bin/buildprog/buildprog.c b/bin/buildprog/buildprog.c
index 79beb3f..92f956d 100644
--- a/bin/buildprog/buildprog.c
+++ b/bin/buildprog/buildprog.c
@@ -29,6 +29,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <string.h>
 
 #include <gnuastro/list.h>
+#include <gnuastro/pointer.h>
 
 #include <main.h>
 
@@ -50,7 +51,7 @@ buildprog_as_one_string(char *opt, gal_list_str_t *list)
         len += 1 + (opt ? strlen(opt) : 0) + strlen(tmp->v);
 
       /* Allocate space for the string. */
-      out=gal_data_malloc_array(GAL_TYPE_UINT8, len+1, __func__, "out");
+      out=gal_pointer_allocate(GAL_TYPE_UINT8, len+1, 0, __func__, "out");
 
       /* Write all the strings into the allocated space. */
       len=0;
diff --git a/bin/convolve/convolve.c b/bin/convolve/convolve.c
index 94354a3..ae23711 100644
--- a/bin/convolve/convolve.c
+++ b/bin/convolve/convolve.c
@@ -32,6 +32,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/wcs.h>
 #include <gnuastro/tile.h>
 #include <gnuastro/fits.h>
+#include <gnuastro/pointer.h>
 #include <gnuastro/threads.h>
 #include <gnuastro/convolve.h>
 
@@ -58,8 +59,8 @@ complextoreal(double *c, size_t size, int action, double 
**output)
   double *out, *o, *of;
 
   /* Allocate the space for the real array. */
-  *output=out=gal_data_malloc_array(GAL_TYPE_FLOAT64, size, __func__,
-                                    "output");
+  *output=out=gal_pointer_allocate(GAL_TYPE_FLOAT64, size, 0, __func__,
+                                   "output");
 
   /* Fill the real array with the derived value from the complex array. */
   of=(o=out)+size;
@@ -207,8 +208,8 @@ frequency_make_padded_complex(struct convolveparams *p)
 
 
   /* Allocate the space for the padded input image and fill it. */
-  pimg=p->pimg=gal_data_malloc_array(GAL_TYPE_FLOAT64, 2*ps0*ps1, __func__,
-                                     "pimg");
+  pimg=p->pimg=gal_pointer_allocate(GAL_TYPE_FLOAT64, 2*ps0*ps1, 0,
+                                    __func__, "pimg");
   for(i=0;i<ps0;++i)
     {
       op=(o=pimg+i*2*ps1)+2*ps1; /* pimg is complex.            */
@@ -222,8 +223,8 @@ frequency_make_padded_complex(struct convolveparams *p)
 
 
   /* Allocate the space for the padded Kernel and fill it. */
-  pker=p->pker=gal_data_malloc_array(GAL_TYPE_FLOAT64, 2*ps0*ps1, __func__,
-                                     "pker");
+  pker=p->pker=gal_pointer_allocate(GAL_TYPE_FLOAT64, 2*ps0*ps1, 0,
+                                    __func__, "pker");
   for(i=0;i<ps0;++i)
     {
       op=(o=pker+i*2*ps1)+2*ps1; /* pker is complex.            */
diff --git a/bin/crop/onecrop.c b/bin/crop/onecrop.c
index cabe7b9..2cd6eba 100644
--- a/bin/crop/onecrop.c
+++ b/bin/crop/onecrop.c
@@ -35,6 +35,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/fits.h>
 #include <gnuastro/blank.h>
 #include <gnuastro/polygon.h>
+#include <gnuastro/pointer.h>
 
 #include <gnuastro-internal/timing.h>
 #include <gnuastro-internal/checkset.h>
@@ -746,7 +747,7 @@ onecrop(struct onecropparams *crp)
          the desired pixels into it. */
       status=0;
       for(i=0;i<ndim;++i) cropsize *= ( lpixel_i[i] - fpixel_i[i] + 1 );
-      array=gal_data_malloc_array(p->type, cropsize, __func__, "array");
+      array=gal_pointer_allocate(p->type, cropsize, 0, __func__, "array");
       if(fits_read_subset(ifp, gal_fits_type_to_datatype(p->type),
                           fpixel_i, lpixel_i, inc, p->bitnul, array,
                           &anynul, &status))
@@ -900,7 +901,7 @@ onecrop_center_filled(struct onecropparams *crp)
   */
 
   /* Allocate the array and read in the pixels. */
-  array=gal_data_malloc_array(type, size, __func__, "array");
+  array=gal_pointer_allocate(type, size, 0, __func__, "array");
   if( fits_read_subset(ofp, gal_fits_type_to_datatype(type), fpixel, lpixel,
                        inc, p->bitnul, array, &anynul, &status) )
      gal_fits_io_error(status, NULL);
diff --git a/bin/crop/ui.c b/bin/crop/ui.c
index f03eed3..13f2428 100644
--- a/bin/crop/ui.c
+++ b/bin/crop/ui.c
@@ -33,6 +33,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/fits.h>
 #include <gnuastro/blank.h>
 #include <gnuastro/table.h>
+#include <gnuastro/pointer.h>
 #include <gnuastro/dimension.h>
 
 #include <gnuastro-internal/timing.h>
@@ -703,9 +704,9 @@ ui_prepare_center(struct cropparams *p)
       carray=p->center->array;
       for(i=0;i<ndim;++i)
         {
-          p->centercoords[i]=gal_data_malloc_array(GAL_TYPE_FLOAT64,
-                                                   1, __func__,
-                                                   "p->centercoords[i]");
+          p->centercoords[i]=gal_pointer_allocate(GAL_TYPE_FLOAT64,
+                                                  1, 0, __func__,
+                                                  "p->centercoords[i]");
           p->centercoords[i][0]=carray[i];
         }
     }
diff --git a/bin/crop/wcsmode.c b/bin/crop/wcsmode.c
index c4aa8f3..049cbc9 100644
--- a/bin/crop/wcsmode.c
+++ b/bin/crop/wcsmode.c
@@ -31,6 +31,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <stdlib.h>
 
 #include <gnuastro/wcs.h>
+#include <gnuastro/pointer.h>
 
 #include "main.h"
 
@@ -463,8 +464,8 @@ fillcrpipolygon(struct onecropparams *crp)
 
   /* Allocate the image polygon array, and put the image polygon vertice
      values into it. */
-  crp->ipolygon=gal_data_malloc_array(GAL_TYPE_FLOAT64, ndim*p->nvertices,
-                                      __func__, "crp->ipolygon");
+  crp->ipolygon=gal_pointer_allocate(GAL_TYPE_FLOAT64, ndim*p->nvertices, 0,
+                                     __func__, "crp->ipolygon");
   for(i=0;i<p->nvertices;++i)
     {
       d=0;
diff --git a/bin/fits/fits.c b/bin/fits/fits.c
index 895e7c0..0381d72 100644
--- a/bin/fits/fits.c
+++ b/bin/fits/fits.c
@@ -29,6 +29,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/list.h>
 #include <gnuastro/fits.h>
 #include <gnuastro/blank.h>
+#include <gnuastro/pointer.h>
 
 #include <gnuastro-internal/timing.h>
 #include <gnuastro-internal/checkset.h>
@@ -146,7 +147,8 @@ fits_print_extension_info(struct fitsparams *p)
         case BINARY_TBL:
           ndim=2;
           tstr = hdutype==ASCII_TBL ? "table_ascii" : "table_binary";
-          dsize=gal_data_malloc_array(GAL_TYPE_SIZE_T, 2, __func__, "dsize");
+          dsize=gal_pointer_allocate(GAL_TYPE_SIZE_T, 2, 0, __func__,
+                                     "dsize");
           gal_fits_tab_size(fptr, dsize+1, dsize);
           break;
 
diff --git a/bin/match/match.c b/bin/match/match.c
index 3841afd..bfbbf21 100644
--- a/bin/match/match.c
+++ b/bin/match/match.c
@@ -30,6 +30,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 #include <gnuastro/match.h>
 #include <gnuastro/table.h>
+#include <gnuastro/pointer.h>
 #include <gnuastro/permutation.h>
 
 #include <gnuastro-internal/checkset.h>
@@ -60,9 +61,9 @@ match_catalog_read_write_all(struct matchparams *p, size_t 
*permutation,
   /* When the output contains columns from both inputs, we need to keep the
      number of columns matched against each column identifier. */
   if(p->outcols)
-    *numcolmatch=gal_data_malloc_array(GAL_TYPE_SIZE_T,
-                                       gal_list_str_number(cols), __func__,
-                                       "numcolmatch");
+    *numcolmatch=gal_pointer_allocate(GAL_TYPE_SIZE_T,
+                                      gal_list_str_number(cols), 0,
+                                      __func__, "numcolmatch");
 
   /* Read the full table. */
   cat=gal_table_read(filename, hdu, cols, p->cp.searchin, p->cp.ignorecase,
@@ -87,8 +88,7 @@ match_catalog_read_write_all(struct matchparams *p, size_t 
*permutation,
 
           /* Reset the data structure's array element to start where the
              non-matched elements start. */
-          tmp->array=gal_data_ptr_increment(tmp->array, nummatched,
-                                            tmp->type);
+          tmp->array=gal_pointer_increment(tmp->array, nummatched, tmp->type);
 
           /* Correct the size of the tile. */
           tmp->size = tmp->dsize[0] = tmp->size - nummatched;
diff --git a/bin/mkcatalog/columns.c b/bin/mkcatalog/columns.c
index b09d64e..7fcf0b2 100644
--- a/bin/mkcatalog/columns.c
+++ b/bin/mkcatalog/columns.c
@@ -30,6 +30,8 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <string.h>
 #include <pthread.h>
 
+#include <gnuastro/pointer.h>
+
 #include <gnuastro-internal/checkset.h>
 
 #include "main.h"
@@ -335,10 +337,10 @@ columns_define_alloc(struct mkcatalogparams *p)
      smaller domain of raw measurements. So to avoid having to calculate
      something multiple times, each parameter will flag the intermediate
      parameters it requires in these arrays. */
-  oiflag = p->oiflag = gal_data_calloc_array(GAL_TYPE_UINT8, OCOL_NUMCOLS,
-                                             __func__, "oiflag");
-  ciflag = p->ciflag = gal_data_calloc_array(GAL_TYPE_UINT8, CCOL_NUMCOLS,
-                                             __func__, "ciflag");
+  oiflag = p->oiflag = gal_pointer_allocate(GAL_TYPE_UINT8, OCOL_NUMCOLS, 1,
+                                            __func__, "oiflag");
+  ciflag = p->ciflag = gal_pointer_allocate(GAL_TYPE_UINT8, CCOL_NUMCOLS, 1,
+                                            __func__, "ciflag");
 
   /* Allocate the columns. */
   for(colcode=p->columnids; colcode!=NULL; colcode=colcode->next)
diff --git a/bin/mkcatalog/mkcatalog.c b/bin/mkcatalog/mkcatalog.c
index 683b0e6..3d39c67 100644
--- a/bin/mkcatalog/mkcatalog.c
+++ b/bin/mkcatalog/mkcatalog.c
@@ -36,6 +36,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/data.h>
 #include <gnuastro/fits.h>
 #include <gnuastro/threads.h>
+#include <gnuastro/pointer.h>
 #include <gnuastro/dimension.h>
 #include <gnuastro/statistics.h>
 #include <gnuastro/permutation.h>
@@ -105,8 +106,8 @@ mkcatalog_single_object(void *in_prm)
   pp.p               = p;
   pp.clumpstartindex = 0;
   pp.rng             = p->rng ? gsl_rng_clone(p->rng) : NULL;
-  pp.oi              = gal_data_malloc_array(GAL_TYPE_FLOAT64, OCOL_NUMCOLS,
-                                             __func__, "pp.oi");
+  pp.oi              = gal_pointer_allocate(GAL_TYPE_FLOAT64, OCOL_NUMCOLS,
+                                            0, __func__, "pp.oi");
 
   /* If we have second order measurements, allocate the array keeping the
      temporary shift values for each object of this thread. Note that the
@@ -118,7 +119,7 @@ mkcatalog_single_object(void *in_prm)
                  || oif[ OCOL_VXX ]
                  || oif[ OCOL_VYY ]
                  || oif[ OCOL_VXY ] )
-               ? gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim, __func__,
+               ? gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__,
                                        "pp.shift")
                : NULL );
 
@@ -150,9 +151,9 @@ mkcatalog_single_object(void *in_prm)
       if(p->clumps)
         {
           /* Allocate space for the properties of each clump. */
-          pp.ci = gal_data_calloc_array(GAL_TYPE_FLOAT64,
-                                        pp.clumpsinobj * CCOL_NUMCOLS,
-                                        __func__, "pp.ci");
+          pp.ci = gal_pointer_allocate(GAL_TYPE_FLOAT64,
+                                       pp.clumpsinobj * CCOL_NUMCOLS, 1,
+                                       __func__, "pp.ci");
 
           /* Get the starting row of this object's clumps in the final
              catalog. This index is also necessary for the unique random
@@ -560,9 +561,9 @@ sort_clumps_by_objid(struct mkcatalogparams *p)
 
 
   /* Allocate the necessary arrays. */
-  rowstart=gal_data_malloc_array(GAL_TYPE_SIZE_T, p->numobjects, __func__,
+  rowstart=gal_pointer_allocate(GAL_TYPE_SIZE_T, p->numobjects, 0, __func__,
                                  "rowstart");
-  permute=gal_data_malloc_array(GAL_TYPE_SIZE_T, p->numclumps, __func__,
+  permute=gal_pointer_allocate(GAL_TYPE_SIZE_T, p->numclumps, 0, __func__,
                                 "permute");
 
 
diff --git a/bin/mkcatalog/parse.c b/bin/mkcatalog/parse.c
index c5a22e2..c3bdd0d 100644
--- a/bin/mkcatalog/parse.c
+++ b/bin/mkcatalog/parse.c
@@ -31,6 +31,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <stdlib.h>
 
 #include <gnuastro/data.h>
+#include <gnuastro/pointer.h>
 #include <gnuastro/dimension.h>
 #include <gnuastro/statistics.h>
 
@@ -125,8 +126,8 @@ parse_objects(struct mkcatalog_passparams *pp)
 
   /* Coordinate shift. */
   size_t *sc = ( pp->shift
-                 ? gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim, __func__,
-                                         "sc")
+                 ? gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__,
+                                        "sc")
                  : NULL );
 
   /* If any coordinate columns are requested. */
@@ -146,7 +147,7 @@ parse_objects(struct mkcatalog_passparams *pp)
                     the coordinate to find which tile a pixel belongs
                     to. */
                  || tid==GAL_BLANK_SIZE_T )
-               ? gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim, __func__, "c")
+               ? gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__, "c")
                : NULL );
 
 
@@ -340,8 +341,8 @@ parse_clumps(struct mkcatalog_passparams *pp)
 
   /* Coordinate shift. */
   size_t *sc = ( pp->shift
-                 ? gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim, __func__,
-                                         "sc")
+                 ? gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__,
+                                        "sc")
                  : NULL );
 
   /* If any coordinate columns are requested. */
@@ -353,15 +354,16 @@ parse_clumps(struct mkcatalog_passparams *pp)
                   || cif[ CCOL_VZ ]
                   || sc
                   || tid==GAL_BLANK_SIZE_T )
-                ? gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim, __func__, "c")
+                ? gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__,
+                                       "c")
                 : NULL );
 
   /* Preparations for neighbor parsing. */
   int32_t *ngblabs=( ( cif[    CCOL_RIV_NUM     ]
                        || cif[ CCOL_RIV_SUM     ]
                        || cif[ CCOL_RIV_SUM_VAR ] )
-                     ? gal_data_malloc_array(GAL_TYPE_INT32, nngb, __func__,
-                                             "ngblabs")
+                     ? gal_pointer_allocate(GAL_TYPE_INT32, nngb, 0,
+                                             __func__, "ngblabs")
                      : NULL );
   size_t *dinc = ngblabs ? gal_dimension_increment(ndim, dsize) : NULL;
 
@@ -590,8 +592,8 @@ parse_median(struct mkcatalog_passparams *pp)
 
 
       /* Allocate the array necessary to keep the values of each clump. */
-      ccounter=gal_data_calloc_array(GAL_TYPE_SIZE_T, pp->clumpsinobj,
-                                     __func__, "ccounter");
+      ccounter=gal_pointer_allocate(GAL_TYPE_SIZE_T, pp->clumpsinobj, 1,
+                                    __func__, "ccounter");
       for(i=0;i<pp->clumpsinobj;++i)
         {
           tsize=pp->ci[ i * CCOL_NUMCOLS + CCOL_NUM ];
@@ -619,13 +621,13 @@ parse_median(struct mkcatalog_passparams *pp)
           if( *O==pp->object && !( p->hasblank && isnan(*V) ) )
             {
               /* Copy the value for the whole object. */
-              memcpy( gal_data_ptr_increment(objmed->array, counter++,
+              memcpy( gal_pointer_increment(objmed->array, counter++,
                                              p->values->type), V,
                       gal_type_sizeof(p->values->type) );
 
               /* We are also on a clump. */
               if(p->clumps && *C>0)
-                memcpy( gal_data_ptr_increment(clumpsmed[*C-1]->array,
+                memcpy( gal_pointer_increment(clumpsmed[*C-1]->array,
                                                ccounter[*C-1]++,
                                                p->values->type), V,
                         gal_type_sizeof(p->values->type) );
diff --git a/bin/mkcatalog/ui.c b/bin/mkcatalog/ui.c
index 65c1142..0422049 100644
--- a/bin/mkcatalog/ui.c
+++ b/bin/mkcatalog/ui.c
@@ -34,6 +34,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/blank.h>
 #include <gnuastro/array.h>
 #include <gnuastro/threads.h>
+#include <gnuastro/pointer.h>
 #include <gnuastro/dimension.h>
 #include <gnuastro/arithmetic.h>
 #include <gnuastro/statistics.h>
@@ -606,7 +607,7 @@ ui_read_labels(struct mkcatalogparams *p)
                                         p->cp.minmapsize);
 
       /* Check its size. */
-      if( gal_data_dsize_is_different(p->objects, p->clumps) )
+      if( gal_dimension_is_different(p->objects, p->clumps) )
         error(EXIT_FAILURE, 0, "`%s' (hdu: %s) and `%s' (hdu: %s) have a"
               "different dimension/size", p->usedclumpsfile, p->clumpshdu,
               p->objectsfile, p->cp.hdu);
@@ -796,7 +797,7 @@ ui_preparation_check_size_read_tiles(struct mkcatalogparams 
*p,
   struct gal_tile_two_layer_params *tl=&p->cp.tl;
 
   /* See if we should treat this dataset as tile values or not. */
-  if( gal_data_dsize_is_different(p->objects, in) )
+  if( gal_dimension_is_different(p->objects, in) )
     {
       /* The `tl' structure is initialized here. But this function may be
          called multiple times. So, first check if the `tl' structure has
@@ -845,7 +846,7 @@ ui_subtract_sky(struct mkcatalogparams *p)
   struct gal_tile_two_layer_params *tl=&p->cp.tl;
 
   /* It is the same size as the input or a single value. */
-  if( gal_data_dsize_is_different(p->values, p->sky)==0 || p->sky->size==1)
+  if( gal_dimension_is_different(p->values, p->sky)==0 || p->sky->size==1)
     {
       s=p->sky->array;
       ff = (f=p->values->array) + p->values->size;
@@ -910,7 +911,7 @@ ui_preparations_read_inputs(struct mkcatalogparams *p)
                                               p->cp.minmapsize);
 
       /* Make sure it has the correct size. */
-      if( gal_data_dsize_is_different(p->objects, p->values) )
+      if( gal_dimension_is_different(p->objects, p->values) )
         error(EXIT_FAILURE, 0, "`%s' (hdu: %s) and `%s' (hdu: %s) have a"
               "different dimension/size", p->usedvaluesfile, p->valueshdu,
               p->objectsfile, p->cp.hdu);
@@ -1019,7 +1020,7 @@ ui_preparations_read_inputs(struct mkcatalogparams *p)
                                             p->cp.minmapsize);
 
           /* Check its size. */
-          if( gal_data_dsize_is_different(p->objects, p->upmask) )
+          if( gal_dimension_is_different(p->objects, p->upmask) )
             error(EXIT_FAILURE, 0, "`%s' (hdu: %s) and `%s' (hdu: %s) have a"
                   "different dimension/size", p->upmaskfile, p->upmaskhdu,
                   p->objectsfile, p->cp.hdu);
@@ -1238,10 +1239,10 @@ ui_one_tile_per_object(struct mkcatalogparams *p)
 
   int32_t *l, *lf, *start;
   size_t i, d, *min, *max, width=2*ndim;
-  size_t *minmax=gal_data_malloc_array(GAL_TYPE_SIZE_T,
-                                       width*p->numobjects, __func__,
-                                       "minmax");
-  size_t *coord=gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim, __func__,
+  size_t *minmax=gal_pointer_allocate(GAL_TYPE_SIZE_T,
+                                      width*p->numobjects, 0, __func__,
+                                      "minmax");
+  size_t *coord=gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__,
                                       "coord");
 
 
@@ -1411,12 +1412,12 @@ ui_preparations(struct mkcatalogparams *p)
      together. */
   if(p->clumps && !p->noclumpsort && p->cp.numthreads>1)
     {
-      p->hostobjid_c=gal_data_malloc_array(GAL_TYPE_SIZE_T,
-                                           p->clumpcols->size, __func__,
-                                           "p->hostobjid_c");
-      p->numclumps_c=gal_data_malloc_array(GAL_TYPE_SIZE_T,
-                                           p->objectcols->size, __func__,
-                                           "p->numclumps_c");
+      p->hostobjid_c=gal_pointer_allocate(GAL_TYPE_SIZE_T,
+                                          p->clumpcols->size, 0, __func__,
+                                          "p->hostobjid_c");
+      p->numclumps_c=gal_pointer_allocate(GAL_TYPE_SIZE_T,
+                                          p->objectcols->size, 0, __func__,
+                                          "p->numclumps_c");
     }
 }
 
diff --git a/bin/mkcatalog/upperlimit.c b/bin/mkcatalog/upperlimit.c
index e571641..eb29a04 100644
--- a/bin/mkcatalog/upperlimit.c
+++ b/bin/mkcatalog/upperlimit.c
@@ -31,6 +31,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 #include <gnuastro/tile.h>
 #include <gnuastro/threads.h>
+#include <gnuastro/pointer.h>
 #include <gnuastro/dimension.h>
 #include <gnuastro/statistics.h>
 
@@ -54,11 +55,11 @@ upperlimit_make_clump_tiles(struct mkcatalog_passparams *pp)
   size_t increment=0, num_increment=1;
   size_t i, d, *min, *max, width=2*ndim;
   int32_t *O, *OO, *C, *start=objects->array;
-  size_t *coord=gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim, __func__,
-                                      "coord");
-  size_t *minmax=gal_data_malloc_array(GAL_TYPE_SIZE_T,
-                                       width*pp->clumpsinobj, __func__,
-                                       "minmax");
+  size_t *coord=gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__,
+                                     "coord");
+  size_t *minmax=gal_pointer_allocate(GAL_TYPE_SIZE_T,
+                                      width*pp->clumpsinobj, 0, __func__,
+                                      "minmax");
 
   /* Initialize the minimum and maximum position for each tile/clump. So,
      we'll initialize the minimum coordinates to the maximum possible
@@ -160,8 +161,8 @@ upperlimit_random_range(struct mkcatalog_passparams *pp, 
gal_data_t *tile,
   /* Set the minimum and maximum acceptable value for the range.  */
   if(p->uprange)
     {
-      tstart=gal_data_num_between(tile->block->array, tile->array,
-                                  p->objects->type);
+      tstart=gal_pointer_num_between(tile->block->array, tile->array,
+                                     p->objects->type);
       gal_dimension_index_to_coord(tstart, ndim, dsize, coord);
     }
 
@@ -576,8 +577,8 @@ upperlimit_one_tile(struct mkcatalog_passparams *pp, 
gal_data_t *tile,
   int32_t *O, *OO, *oO, *st_o, *st_oo, *st_oc, *oC=NULL;
   size_t maxcount = p->upnum * MKCATALOG_UPPERLIMIT_STOP_MULTIP;
   struct gal_list_sizet_t *check_x=NULL, *check_y=NULL, *check_z=NULL;
-  size_t *rcoord=gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim, __func__,
-                                       "rcoord");
+  size_t *rcoord=gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__,
+                                      "rcoord");
 
   /* See if a check table must be created for this distribution. */
   if( p->checkupperlimit[0]==pp->object )
@@ -619,7 +620,7 @@ upperlimit_one_tile(struct mkcatalog_passparams *pp, 
gal_data_t *tile,
         rcoord[d] = upperlimit_random_position(pp, tile, d, min, max);
 
       /* Set the tile's new starting pointer. */
-      tile->array = gal_data_ptr_increment(p->objects->array,
+      tile->array = gal_pointer_increment(p->objects->array,
                           gal_dimension_coord_to_index(ndim, dsize, rcoord),
                                            p->objects->type);
 
diff --git a/bin/mkprof/mkprof.c b/bin/mkprof/mkprof.c
index 740e12d..c5ad3a7 100644
--- a/bin/mkprof/mkprof.c
+++ b/bin/mkprof/mkprof.c
@@ -33,6 +33,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/git.h>
 #include <gnuastro/fits.h>
 #include <gnuastro/threads.h>
+#include <gnuastro/pointer.h>
 #include <gnuastro/dimension.h>
 #include <gnuastro/statistics.h>
 
@@ -152,7 +153,8 @@ saveindividual(struct mkonthread *mkp)
     {
       /* Allocate space for the corrected crpix and fill it in. Both
          `crpix' and `fpixel_i' are in FITS order. */
-      crpix=gal_data_malloc_array(GAL_TYPE_FLOAT64, ndim, __func__, "crpix");
+      crpix=gal_pointer_allocate(GAL_TYPE_FLOAT64, ndim, 0, __func__,
+                                 "crpix");
       for(i=0;i<ndim;++i)
         crpix[i] = ((double *)(p->crpix->array))[i] - os*(mkp->fpixel_i[i]-1);
 
@@ -259,8 +261,7 @@ mkprof_build_single(struct mkonthread *mkp, long *fpixel_i, 
long *lpixel_i,
           /* If a crop is needed, set the starting pointer. */
           ind=gal_dimension_coord_to_index(ndim, ibq->image->dsize,
                                            start_indiv);
-          ptr=gal_data_ptr_increment(ibq->image->array, ind,
-                                     ibq->image->type);
+          ptr=gal_pointer_increment(ibq->image->array, ind, ibq->image->type);
         }
       else ptr=ibq->image->array;
       ibq->overlap_i=gal_data_alloc(ptr, ibq->image->type, ndim, dsize, NULL,
@@ -270,7 +271,7 @@ mkprof_build_single(struct mkonthread *mkp, long *fpixel_i, 
long *lpixel_i,
 
       /* Define the merged overlap tile. */
       ind=gal_dimension_coord_to_index(ndim, p->out->dsize, start_mrg);
-      ptr=gal_data_ptr_increment(p->out->array, ind, p->out->type);
+      ptr=gal_pointer_increment(p->out->array, ind, p->out->type);
       ibq->overlap_m=gal_data_alloc(ptr, p->out->type, ndim, dsize, NULL,
                                     0, -1, NULL, NULL, NULL);
       ibq->overlap_m->block=p->out;
@@ -665,7 +666,7 @@ mkprof(struct mkprofparams *p)
      ignore it. */
   if(p->out)
     {
-      onaxes=gal_data_malloc_array(GAL_TYPE_LONG, ndim, __func__, "onaxes");
+      onaxes=gal_pointer_allocate(GAL_TYPE_LONG, ndim, 0, __func__, "onaxes");
       for(fi=0; fi < ndim; ++fi)
         {
           i=ndim-fi-1;
diff --git a/bin/mkprof/oneprofile.c b/bin/mkprof/oneprofile.c
index 79b3d29..06935ab 100644
--- a/bin/mkprof/oneprofile.c
+++ b/bin/mkprof/oneprofile.c
@@ -34,6 +34,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gsl/gsl_integration.h> /* gsl_integration_qng  */
 
 #include <gnuastro/fits.h>
+#include <gnuastro/pointer.h>
 #include <gnuastro/dimension.h>
 #include <gnuastro/statistics.h>
 
@@ -381,9 +382,9 @@ oneprofile_pix_by_pix(struct mkonthread *mkp)
 
   /* Allocate the `byt' array. It is used as a flag to make sure that we
      don't re-calculate the profile value on a pixel more than once. */
-  byt = gal_data_calloc_array(GAL_TYPE_UINT8,
-                              gal_dimension_total_size(ndim, dsize),
-                              __func__, "byt");
+  byt = gal_pointer_allocate(GAL_TYPE_UINT8,
+                             gal_dimension_total_size(ndim, dsize), 1,
+                             __func__, "byt");
 
   /* Start the queue: */
   byt[p]=1;
diff --git a/bin/mkprof/ui.c b/bin/mkprof/ui.c
index 6e3ffb9..5fc7d22 100644
--- a/bin/mkprof/ui.c
+++ b/bin/mkprof/ui.c
@@ -34,6 +34,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/array.h>
 #include <gnuastro/blank.h>
 #include <gnuastro/table.h>
+#include <gnuastro/pointer.h>
 
 #include <gnuastro-internal/timing.h>
 #include <gnuastro-internal/options.h>
@@ -689,8 +690,8 @@ ui_read_cols_2d(struct mkprofparams *p)
         case 3:
           if(tmp->type==GAL_TYPE_STRING)
             {
-              p->f=gal_data_malloc_array(GAL_TYPE_UINT8, p->num,
-                                         __func__, "p->f");
+              p->f=gal_pointer_allocate(GAL_TYPE_UINT8, p->num, 0,
+                                        __func__, "p->f");
               strarr=tmp->array;
               for(i=0;i<p->num;++i)
                 p->f[i]=ui_profile_name_read(strarr[i], i+1);
@@ -900,8 +901,8 @@ ui_read_cols_3d(struct mkprofparams *p)
         case 4:
           if(tmp->type==GAL_TYPE_STRING)
             {
-              p->f=gal_data_malloc_array(GAL_TYPE_UINT8, p->num,
-                                         __func__, "p->f");
+              p->f=gal_pointer_allocate(GAL_TYPE_UINT8, p->num, 0,
+                                        __func__, "p->f");
               strarr=tmp->array;
               for(i=0;i<p->num;++i)
                 p->f[i]=ui_profile_name_read(strarr[i], i+1);
@@ -1041,7 +1042,6 @@ ui_read_cols_3d(struct mkprofparams *p)
 
 
 
-
 /* It is possible to define the internal catalog through a catalog or the
    `--kernel' option. This function will do the job. */
 static void
@@ -1058,21 +1058,21 @@ ui_prepare_columns(struct mkprofparams *p)
       p->num=1;
 
       /* Allocate the necessary columns. */
-      p->x  = gal_data_calloc_array(GAL_TYPE_FLOAT64, 1, __func__, "p->x");
-      p->y  = gal_data_calloc_array(GAL_TYPE_FLOAT64, 1, __func__, "p->y");
-      p->f  = gal_data_calloc_array(GAL_TYPE_UINT8,   1, __func__, "p->f");
-      p->r  = gal_data_calloc_array(GAL_TYPE_FLOAT32, 1, __func__, "p->r");
-      p->n  = gal_data_calloc_array(GAL_TYPE_FLOAT32, 1, __func__, "p->n");
-      p->p1 = gal_data_calloc_array(GAL_TYPE_FLOAT32, 1, __func__, "p->p1");
-      p->q1 = gal_data_calloc_array(GAL_TYPE_FLOAT32, 1, __func__, "p->q1");
-      p->m  = gal_data_calloc_array(GAL_TYPE_FLOAT32, 1, __func__, "p->m");
-      p->t  = gal_data_calloc_array(GAL_TYPE_FLOAT32, 1, __func__, "p->t");
+      p->x  = gal_pointer_allocate(GAL_TYPE_FLOAT64, 1, 1, __func__, "p->x");
+      p->y  = gal_pointer_allocate(GAL_TYPE_FLOAT64, 1, 1, __func__, "p->y");
+      p->f  = gal_pointer_allocate(GAL_TYPE_UINT8,   1, 1, __func__, "p->f");
+      p->r  = gal_pointer_allocate(GAL_TYPE_FLOAT32, 1, 1, __func__, "p->r");
+      p->n  = gal_pointer_allocate(GAL_TYPE_FLOAT32, 1, 1, __func__, "p->n");
+      p->p1 = gal_pointer_allocate(GAL_TYPE_FLOAT32, 1, 1, __func__, "p->p1");
+      p->q1 = gal_pointer_allocate(GAL_TYPE_FLOAT32, 1, 1, __func__, "p->q1");
+      p->m  = gal_pointer_allocate(GAL_TYPE_FLOAT32, 1, 1, __func__, "p->m");
+      p->t  = gal_pointer_allocate(GAL_TYPE_FLOAT32, 1, 1, __func__, "p->t");
       if(p->ndim==3)
         {
-          p->z =gal_data_calloc_array(GAL_TYPE_FLOAT64, 1, __func__, "p->z");
-          p->p2=gal_data_calloc_array(GAL_TYPE_FLOAT32, 1, __func__, "p->p2");
-          p->p3=gal_data_calloc_array(GAL_TYPE_FLOAT32, 1, __func__, "p->p3");
-          p->q2=gal_data_calloc_array(GAL_TYPE_FLOAT32, 1, __func__, "p->q2");
+          p->z =gal_pointer_allocate(GAL_TYPE_FLOAT64, 1, 1, __func__, "p->z");
+          p->p2=gal_pointer_allocate(GAL_TYPE_FLOAT32, 1, 1, __func__, 
"p->p2");
+          p->p3=gal_pointer_allocate(GAL_TYPE_FLOAT32, 1, 1, __func__, 
"p->p3");
+          p->q2=gal_pointer_allocate(GAL_TYPE_FLOAT32, 1, 1, __func__, 
"p->q2");
         }
 
       /* For profiles that need a different number of input values. Note
@@ -1348,8 +1348,8 @@ ui_prepare_canvas(struct mkprofparams *p)
           if( p->dsize ) free(p->dsize);
 
           /* Write the size of the background image into `dsize'. */
-          p->dsize=gal_data_malloc_array(GAL_TYPE_SIZE_T, p->ndim, __func__,
-                                         "p->dsize");
+          p->dsize=gal_pointer_allocate(GAL_TYPE_SIZE_T, p->ndim, 0,
+                                        __func__, "p->dsize");
           for(i=0;i<p->ndim;++i) p->dsize[i] = p->out->dsize[i];
 
           /* Set all pixels to zero if the user wanted a clear canvas. */
@@ -1369,8 +1369,8 @@ ui_prepare_canvas(struct mkprofparams *p)
          there is no shifts. */
       p->oversample=1;
       if(p->shift) free(p->shift);
-      p->shift=gal_data_calloc_array(GAL_TYPE_SIZE_T, p->ndim, __func__,
-                                     "p->shift (1)");
+      p->shift=gal_pointer_allocate(GAL_TYPE_SIZE_T, p->ndim, 1, __func__,
+                                    "p->shift (1)");
     }
   else
     {
@@ -1441,8 +1441,8 @@ ui_prepare_canvas(struct mkprofparams *p)
                  shifts (from zero). So, we'll just free it and reset
                  it. */
               if(p->shift) free(p->shift);
-              p->shift=gal_data_calloc_array(GAL_TYPE_SIZE_T, p->ndim,
-                                             __func__, "p->shift (2)");
+              p->shift=gal_pointer_allocate(GAL_TYPE_SIZE_T, p->ndim, 1,
+                                            __func__, "p->shift (2)");
               if(setshift)
                 {
                   p->shift[0]  = (width[0]/2)*p->oversample;
@@ -1454,8 +1454,8 @@ ui_prepare_canvas(struct mkprofparams *p)
 
       /* If shift has not been set until now, set it. */
       if(p->shift==NULL)
-        p->shift=gal_data_calloc_array(GAL_TYPE_SIZE_T, p->ndim, __func__,
-                                       "p->shift (3)");
+        p->shift=gal_pointer_allocate(GAL_TYPE_SIZE_T, p->ndim, 1,
+                                      __func__, "p->shift (3)");
 
       /* Prepare the sizes of the final merged image (if it is to be
          made). Note that even if we don't want a merged image, we still
@@ -1656,8 +1656,8 @@ ui_preparations(struct mkprofparams *p)
       p->ndim=p->kernel->flag;
 
       /* Set the shift array. */
-      p->shift=gal_data_calloc_array(GAL_TYPE_SIZE_T, p->ndim,
-                                     __func__, "p->shift");
+      p->shift=gal_pointer_allocate(GAL_TYPE_SIZE_T, p->ndim, 1,
+                                    __func__, "p->shift");
     }
   else
     ui_prepare_canvas(p);
diff --git a/bin/noisechisel/detection.c b/bin/noisechisel/detection.c
index e0ebf9e..891e41a 100644
--- a/bin/noisechisel/detection.c
+++ b/bin/noisechisel/detection.c
@@ -32,6 +32,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/label.h>
 #include <gnuastro/binary.h>
 #include <gnuastro/threads.h>
+#include <gnuastro/pointer.h>
 #include <gnuastro/dimension.h>
 #include <gnuastro/statistics.h>
 
@@ -346,9 +347,9 @@ detection_pseudo_find(struct noisechiselparams *p, 
gal_data_t *workbin,
      thread. `maxltcontig' is the maximum contiguous patch of memory needed
      to store all tiles. Finally, since we are working on a `uint8_t' type,
      the size of each element is only 1 byte. */
-  fho_prm.copyspace=gal_data_malloc_array(GAL_TYPE_UINT8,
-                                          p->cp.numthreads*p->maxltcontig,
-                                          __func__, "fho_prm.copyspace");
+  fho_prm.copyspace=gal_pointer_allocate(GAL_TYPE_UINT8,
+                                         p->cp.numthreads*p->maxltcontig, 0,
+                                         __func__, "fho_prm.copyspace");
 
 
   /* Fill the holes and open on each large tile. When no check image is
@@ -500,7 +501,7 @@ detection_sn(struct noisechiselparams *p, gal_data_t 
*worklab, size_t num,
   size_t i, j, *area, counter=0, *dsize=p->input->dsize;
   float *img=p->input->array, *f=p->input->array, *ff=f+p->input->size;
   int32_t *plab = worklab->array, *dlab = s0d1D2 ? NULL : p->olabel->array;
-  size_t *coord=gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim, __func__,
+  size_t *coord=gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__,
                                       "coord");
 
 
@@ -520,14 +521,14 @@ detection_sn(struct noisechiselparams *p, gal_data_t 
*worklab, size_t num,
   /* Allocate all the necessary arrays, note that since we want to put each
      object's information into the same index, the number of allocated
      spaces has to be `tablen=num+1'. */
-  area       = gal_data_calloc_array(GAL_TYPE_SIZE_T,  tablen, __func__,
-                                     "area");
-  brightness = gal_data_calloc_array(GAL_TYPE_FLOAT64, tablen, __func__,
+  area       = gal_pointer_allocate(GAL_TYPE_SIZE_T,  tablen, 1, __func__,
+                                    "area");
+  brightness = gal_pointer_allocate(GAL_TYPE_FLOAT64, tablen, 1, __func__,
                                      "brightness");
-  pos        = gal_data_calloc_array(GAL_TYPE_FLOAT64, pcols*tablen,
-                                     __func__, "pos");
+  pos        = gal_pointer_allocate(GAL_TYPE_FLOAT64, pcols*tablen, 1,
+                                    __func__, "pos");
   flag       = ( s0d1D2==0
-                 ? gal_data_calloc_array(GAL_TYPE_UINT8, tablen, __func__,
+                 ? gal_pointer_allocate(GAL_TYPE_UINT8, tablen, 1, __func__,
                                          "flag")
                  : NULL );
   sn         = gal_data_alloc(NULL, GAL_TYPE_FLOAT32, 1, &tablen, NULL, 1,
@@ -658,7 +659,6 @@ detection_sn(struct noisechiselparams *p, gal_data_t 
*worklab, size_t num,
           }
     }
 
-
   /* If we are in Sky mode, the sizes have to be corrected */
   if(s0d1D2==0)
     {
@@ -697,8 +697,8 @@ detection_pseudo_remove_low_sn(struct noisechiselparams *p,
   float *snarr=sn->array;
   uint8_t *b=workbin->array;
   int32_t *l=worklab->array, *lf=l+worklab->size;
-  uint8_t *keep=gal_data_calloc_array(GAL_TYPE_UINT8, sn->size, __func__,
-                                      "keep");
+  uint8_t *keep=gal_pointer_allocate(GAL_TYPE_UINT8, sn->size, 1, __func__,
+                                     "keep");
 
   /* Specify the new labels for those that must be kept/changed. Note that
      when an object didn't have an S/N, its S/N was given a value of NaN
@@ -816,8 +816,8 @@ detection_final_remove_small_sn(struct noisechiselparams *p,
   gal_data_t *sn, *snind;
   int32_t *l, *lf, curlab=1;
   gal_list_str_t *comments=NULL;
-  int32_t *newlabs=gal_data_calloc_array(GAL_TYPE_INT32, num+1, __func__,
-                                         "newlabs");
+  int32_t *newlabs=gal_pointer_allocate(GAL_TYPE_INT32, num+1, 1, __func__,
+                                        "newlabs");
 
   /* Get the Signal to noise ratio of all detections. */
   sn=detection_sn(p, p->olabel, num, 2, "DILATED");
@@ -906,9 +906,9 @@ detection_remove_false_initial(struct noisechiselparams *p,
   uint8_t *b=workbin->array;
   float *e_th, *arr=p->conv->array;
   int32_t *l=p->olabel->array, *lf=l+p->olabel->size, curlab=1;
-  int32_t *newlabels=gal_data_calloc_array(GAL_TYPE_UINT32,
-                                           p->numinitialdets+1, __func__,
-                                           "newlabels");
+  int32_t *newlabels=gal_pointer_allocate(GAL_TYPE_UINT32,
+                                          p->numinitialdets+1, 1, __func__,
+                                          "newlabels");
 
   /* Find the new labels for all the existing labels. Recall that
      `newlabels' was initialized to zero, so any label that is not given a
diff --git a/bin/noisechisel/sky.c b/bin/noisechisel/sky.c
index ca0f35b..2d003c9 100644
--- a/bin/noisechisel/sky.c
+++ b/bin/noisechisel/sky.c
@@ -31,6 +31,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/tile.h>
 #include <gnuastro/fits.h>
 #include <gnuastro/blank.h>
+#include <gnuastro/pointer.h>
 #include <gnuastro/threads.h>
 #include <gnuastro/statistics.h>
 
@@ -56,7 +57,7 @@ sky_mean_std_undetected(void *in_prm)
   struct noisechiselparams *p=(struct noisechiselparams *)tprm->params;
 
   double *darr, s, s2;
-  int type=p->sky->type;
+  int setblank, type=p->sky->type;
   size_t i, tind, numsky, dsize=2;
   gal_data_t *tile, *meanstd_d, *meanstd, *bintile;
 
@@ -93,6 +94,7 @@ sky_mean_std_undetected(void *in_prm)
 
       /* Only continue, if the fraction of Sky values are less than the
          requested fraction. */
+      setblank=0;
       if( (float)(numsky)/(float)(tile->size) > p->minskyfrac)
         {
           /* Calculate the mean and STD over this tile. */
@@ -107,26 +109,41 @@ sky_mean_std_undetected(void *in_prm)
           darr[0]=s/numsky;
           darr[1]=sqrt( (s2-s*s/numsky)/numsky );
 
-          /* Convert the mean and std into the same type as the sky and std
-             arrays. */
-          meanstd=gal_data_copy_to_new_type(meanstd_d, type);
-
-          /* Copy the mean and STD to their respective places in the tile
-             arrays. */
-          memcpy(gal_data_ptr_increment(p->sky->array, tind, type),
-                 meanstd->array, gal_type_sizeof(type));
-          memcpy(gal_data_ptr_increment(p->std->array, tind, type),
-                 gal_data_ptr_increment(meanstd->array, 1, type),
-                 gal_type_sizeof(type));
-
-          /* Clean up. */
-          gal_data_free(meanstd);
+          /* When there are zero-valued pixels on the edges of the dataset
+             (that have not been set to NaN/blank), given special
+             conditions, the whole zero-valued region can get a binary
+             value of 1 and so the Sky and its standard deviation can
+             become zero. So, we need ignore such tiles. */
+          if(darr[1]==0.0f)
+            setblank=1;
+          else
+            {
+              /* Convert the mean and std into the same type as the sky and
+                 std arrays. */
+              meanstd=gal_data_copy_to_new_type(meanstd_d, type);
+
+              /* Copy the mean and STD to their respective places in the
+                 tile arrays. */
+              memcpy(gal_pointer_increment(p->sky->array, tind, type),
+                     meanstd->array, gal_type_sizeof(type));
+              memcpy(gal_pointer_increment(p->std->array, tind, type),
+                     gal_pointer_increment(meanstd->array, 1, type),
+                     gal_type_sizeof(type));
+
+              /* Clean up. */
+              gal_data_free(meanstd);
+            }
         }
       else
+        setblank=1;
+
+      /* If the tile is marked as being blank, write blank values into
+         it. */
+      if(setblank==1)
         {
-          gal_blank_write(gal_data_ptr_increment(p->sky->array, tind, type),
+          gal_blank_write(gal_pointer_increment(p->sky->array, tind, type),
                           type);
-          gal_blank_write(gal_data_ptr_increment(p->std->array, tind, type),
+          gal_blank_write(gal_pointer_increment(p->std->array, tind, type),
                           type);
         }
     }
@@ -203,6 +220,7 @@ sky_and_std(struct noisechiselparams *p, char *checkname)
      factor here. */
   p->cpscorr = p->minstd>1 ? 1.0f : p->minstd;
 
+
   /* Interpolate and smooth the derived values. */
   threshold_interp_smooth(p, &p->sky, &p->std, NULL, checkname);
 
diff --git a/bin/noisechisel/threshold.c b/bin/noisechisel/threshold.c
index 592a745..c489f54 100644
--- a/bin/noisechisel/threshold.c
+++ b/bin/noisechisel/threshold.c
@@ -31,6 +31,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/fits.h>
 #include <gnuastro/blank.h>
 #include <gnuastro/threads.h>
+#include <gnuastro/pointer.h>
 #include <gnuastro/statistics.h>
 #include <gnuastro/interpolate.h>
 
@@ -400,8 +401,8 @@ qthresh_on_tile(void *in_prm)
 
   /* Put the temporary usage space for this thread into a data set for easy
      processing. */
-  usage=gal_data_alloc(gal_data_ptr_increment(qprm->usage,
-                                              tprm->id*p->maxtcontig, type),
+  usage=gal_data_alloc(gal_pointer_increment(qprm->usage,
+                                             tprm->id*p->maxtcontig, type),
                        type, ndim, p->maxtsize, NULL, 0, p->cp.minmapsize,
                        NULL, NULL, NULL);
 
@@ -462,13 +463,13 @@ qthresh_on_tile(void *in_prm)
           /* Get the erosion quantile for this tile and save it. Note that
              the type of `qvalue' is the same as the input dataset. */
           qvalue=gal_statistics_quantile(usage, p->qthresh, 1);
-          memcpy(gal_data_ptr_increment(qprm->erode_th->array, tind, type),
+          memcpy(gal_pointer_increment(qprm->erode_th->array, tind, type),
                  qvalue->array, twidth);
           gal_data_free(qvalue);
 
           /* Same for the no-erode quantile. */
           qvalue=gal_statistics_quantile(usage, p->noerodequant, 1);
-          memcpy(gal_data_ptr_increment(qprm->noerode_th->array, tind, type),
+          memcpy(gal_pointer_increment(qprm->noerode_th->array, tind, type),
                  qvalue->array, twidth);
           gal_data_free(qvalue);
 
@@ -476,7 +477,7 @@ qthresh_on_tile(void *in_prm)
           if(p->detgrowquant!=1.0f)
             {
               qvalue=gal_statistics_quantile(usage, p->detgrowquant, 1);
-              memcpy(gal_data_ptr_increment(qprm->expand_th->array, tind,
+              memcpy(gal_pointer_increment(qprm->expand_th->array, tind,
                                             type),
                      qvalue->array, twidth);
               gal_data_free(qvalue);
@@ -484,12 +485,12 @@ qthresh_on_tile(void *in_prm)
         }
       else
         {
-          gal_blank_write(gal_data_ptr_increment(qprm->erode_th->array,
+          gal_blank_write(gal_pointer_increment(qprm->erode_th->array,
                                                  tind, type), type);
-          gal_blank_write(gal_data_ptr_increment(qprm->noerode_th->array,
+          gal_blank_write(gal_pointer_increment(qprm->noerode_th->array,
                                                  tind, type), type);
           if(p->detgrowquant!=1.0f)
-            gal_blank_write(gal_data_ptr_increment(qprm->expand_th->array,
+            gal_blank_write(gal_pointer_increment(qprm->expand_th->array,
                                                    tind, type), type);
         }
 
@@ -538,11 +539,11 @@ threshold_qthresh_clean_work(struct noisechiselparams *p, 
gal_data_t *first,
       if(third) oa3=third->array;
 
       /* Increment the array pointers. */
-      first->array=gal_data_ptr_increment(first->array, start, first->type);
-      second->array=gal_data_ptr_increment(second->array, start,
+      first->array=gal_pointer_increment(first->array, start, first->type);
+      second->array=gal_pointer_increment(second->array, start,
                                            second->type);
       if(third)
-        third->array=gal_data_ptr_increment(third->array, start, third->type);
+        third->array=gal_pointer_increment(third->array, start, third->type);
 
       /* Correct their sizes. */
       first->size=number;
@@ -697,9 +698,9 @@ threshold_quantile_find_apply(struct noisechiselparams *p)
 
 
   /* Allocate temporary space for processing in each tile. */
-  qprm.usage=gal_data_malloc_array(p->input->type,
-                                   cp->numthreads * p->maxtcontig,
-                                   __func__, "qprm.usage");
+  qprm.usage=gal_pointer_allocate(p->input->type,
+                                  cp->numthreads * p->maxtcontig, 0,
+                                  __func__, "qprm.usage");
 
 
   /* Find the threshold on each tile, free the temporary processing space
diff --git a/bin/noisechisel/ui.c b/bin/noisechisel/ui.c
index 4bb3d12..70204ba 100644
--- a/bin/noisechisel/ui.c
+++ b/bin/noisechisel/ui.c
@@ -33,6 +33,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/array.h>
 #include <gnuastro/blank.h>
 #include <gnuastro/threads.h>
+#include <gnuastro/pointer.h>
 #include <gnuastro/dimension.h>
 
 #include <gnuastro-internal/timing.h>
@@ -67,8 +68,8 @@ args_doc[] = "ASTRdata";
 const char
 doc[] = GAL_STRINGS_TOP_HELP_INFO PROGRAM_NAME" Detects and segments signal "
   "that is deeply burried in noise. It employs a noise-based detection and "
-  "segmentation method enabling it to be very resilient to the rich diversity "
-  "of shapes in astronomical targets.\n"
+  "segmentation method enabling it to be very resilient to the rich "
+  "diversity of shapes in astronomical targets.\n"
   GAL_STRINGS_MORE_HELP_INFO
   /* After the list of options: */
   "\v"
@@ -471,8 +472,8 @@ ui_prepare_tiles(struct noisechiselparams *p)
   /* Check the tile parameters for the small tile sizes and make the tile
      structure. We will also need the dimensions of the tile with the
      maximum required memory. */
-  p->maxtsize=gal_data_malloc_array(GAL_TYPE_SIZE_T, p->input->ndim,
-                                    __func__, "p->maxtsize");
+  p->maxtsize=gal_pointer_allocate(GAL_TYPE_SIZE_T, p->input->ndim, 0,
+                                   __func__, "p->maxtsize");
   gal_tile_full_sanity_check(p->inputname, p->cp.hdu, p->input, tl);
   gal_tile_full_two_layers(p->input, tl);
   gal_tile_full_permutation(tl);
@@ -491,8 +492,8 @@ ui_prepare_tiles(struct noisechiselparams *p)
   ltl->workoverch     = tl->workoverch;
   ltl->checktiles     = tl->checktiles;
   ltl->oneelempertile = tl->oneelempertile;
-  p->maxltsize=gal_data_malloc_array(GAL_TYPE_SIZE_T, p->input->ndim,
-                                     __func__, "p->maxltsize");
+  p->maxltsize=gal_pointer_allocate(GAL_TYPE_SIZE_T, p->input->ndim, 0,
+                                    __func__, "p->maxltsize");
   gal_tile_full_sanity_check(p->inputname, p->cp.hdu, p->input, ltl);
   gal_tile_full_two_layers(p->input, ltl);
   gal_tile_full_permutation(ltl);
@@ -544,6 +545,7 @@ ui_prepare_tiles(struct noisechiselparams *p)
 static void
 ui_preparations_read_input(struct noisechiselparams *p)
 {
+  float *f;
   size_t value;
   char *option_name=NULL, *good_values=NULL;
 
@@ -556,7 +558,6 @@ ui_preparations_read_input(struct noisechiselparams *p)
   if(p->input->name==NULL)
     gal_checkset_allocate_copy("INPUT", &p->input->name);
 
-
   /* Check dimensionality and neighbors. */
   switch(p->input->ndim)
     {
@@ -596,13 +597,36 @@ ui_preparations_read_input(struct noisechiselparams *p)
             p->input->ndim);
     }
 
-
   /* Abort with an error message if necessary. */
   if(option_name)
     error(EXIT_FAILURE, 0, "%zu not acceptable for `--%s' for the "
           "input %zu-dimensional dataset in `%s' (hdu %s). It must be %s "
           "(specifying the type of connectivity)", value, option_name,
           p->input->ndim, p->inputname, p->cp.hdu, good_values);
+
+  /* A small check to see if the edges of the dataset aren't zero valued:
+     they should be masked. */
+  f=p->input->array;
+  if( (f[0]==0.0 && f[1]==0.0)
+      || (f[ p->input->size-1 ]==0.0 && f[ p->input->size-2 ]==0.0) )
+    error(0, 0, "%s (hdu %s): [*** WARNING ***]: The first and/or last few "
+          "pixels have a value of 0.0. As described below, the result of "
+          "this run may thus not be reasonable/optimal.\n\n"
+          "Some data reduction pipelines put 0.0 where there isn't data "
+          "(most commonly on the edges). However, NoiseChisel's "
+          "noise-based detection paradigm starts from the lower values of "
+          "the dataset (not high S/N peaks): its initial threshold is "
+          "mostly below the Sky value (0.0 in processed images). Therefore "
+          "0.0 is meaningful for NoiseChisel and must not be used for a "
+          "blank value.\n\n"
+          "To ignore certain pixels, they must have a blank/NaN value. "
+          "To mask (set to blank/NaN) the 0.0 valued elements, you can use "
+          "Gnuastro's Arithmetic program with a command like this:\n\n"
+          "    $ astarithmetic %s %s 0.0 eq nan where -g%s\n\n"
+          "If the few 0.0 valued pixels on the edges are meaningful for "
+          "your analysis, please ignore this warming message.\n"
+          "--------------------------",
+          p->inputname, p->cp.hdu, p->inputname, p->inputname, p->cp.hdu);
 }
 
 
@@ -615,11 +639,9 @@ ui_preparations(struct noisechiselparams *p)
   /* Prepare the names of the outputs. */
   ui_set_output_names(p);
 
-
-  /* Read the input dataset and check the dimensions. */
+  /* Read the input datasets and do the basic checks.*/
   ui_preparations_read_input(p);
 
-
   /* If a convolved image was given, read it in. Otherwise, read the given
      kernel. */
   if(p->convolvedname)
@@ -630,7 +652,7 @@ ui_preparations(struct noisechiselparams *p)
                                               p->cp.minmapsize);
 
       /* Make sure the convolved image is the same size as the input. */
-      if( gal_data_dsize_is_different(p->input, p->conv) )
+      if( gal_dimension_is_different(p->input, p->conv) )
         error(EXIT_FAILURE, 0, "%s (hdu %s), given to `--convolved' and "
               "`--convolvehdu', is not the same size as NoiseChisel's "
               "input: %s (hdu: %s)", p->convolvedname, p->chdu,
@@ -639,15 +661,12 @@ ui_preparations(struct noisechiselparams *p)
   else
     ui_prepare_kernel(p);
 
-
   /* Check for blank values to help later processing.  */
   gal_blank_present(p->input, 1);
 
-
   /* Prepare the tessellation. */
   ui_prepare_tiles(p);
 
-
   /* Allocate space for the over-all necessary arrays. */
   p->binary=gal_data_alloc(NULL, GAL_TYPE_UINT8, p->input->ndim,
                            p->input->dsize, p->input->wcs, 0,
diff --git a/bin/segment/clumps.c b/bin/segment/clumps.c
index 5877f88..5055216 100644
--- a/bin/segment/clumps.c
+++ b/bin/segment/clumps.c
@@ -32,6 +32,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/blank.h>
 #include <gnuastro/label.h>
 #include <gnuastro/threads.h>
+#include <gnuastro/pointer.h>
 #include <gnuastro/dimension.h>
 #include <gnuastro/statistics.h>
 
@@ -255,7 +256,7 @@ clumps_get_raw_info(struct clumps_thread_params *cltprm)
   int32_t lab, nlab, *ngblabs, *clabel=p->clabel->array;
 
   /* Allocate the array to keep the neighbor labels of river pixels. */
-  ngblabs=gal_data_malloc_array(GAL_TYPE_INT32, nngb, __func__, "ngblabs");
+  ngblabs=gal_pointer_allocate(GAL_TYPE_INT32, nngb, 0, __func__, "ngblabs");
 
   /* Go over all the pixels in this region. */
   af=(a=cltprm->indexs->array)+cltprm->indexs->size;
@@ -398,22 +399,22 @@ clumps_make_sn_table(struct clumps_thread_params *cltprm)
   cltprm->sn        = &cltprm->clprm->sn[ cltprm->id ];
   cltprm->sn->ndim  = 1;                        /* Depends on `cltprm->sn' */
   cltprm->sn->type  = GAL_TYPE_FLOAT32;
-  cltprm->sn->dsize = gal_data_malloc_array(GAL_TYPE_SIZE_T, 1, __func__,
-                                            "cltprm->sn->dsize");
-  cltprm->sn->array = gal_data_malloc_array(cltprm->sn->type, tablen,
-                                            __func__, "cltprm->sn->array");
+  cltprm->sn->dsize = gal_pointer_allocate(GAL_TYPE_SIZE_T, 1, 0, __func__,
+                                           "cltprm->sn->dsize");
+  cltprm->sn->array = gal_pointer_allocate(cltprm->sn->type, tablen, 0,
+                                           __func__, "cltprm->sn->array");
   cltprm->sn->size  = cltprm->sn->dsize[0] = tablen;       /* After dsize. */
   if( cltprm->clprm->snind )
     {
       cltprm->snind        = &cltprm->clprm->snind [ cltprm->id ];
       cltprm->snind->ndim  = 1;              /* Depends on `cltprm->snind' */
       cltprm->snind->type  = GAL_TYPE_INT32;
-      cltprm->snind->dsize = gal_data_malloc_array(GAL_TYPE_SIZE_T, 1,
-                                                   __func__,
-                                                   "cltprm->snind->dsize");
+      cltprm->snind->dsize = gal_pointer_allocate(GAL_TYPE_SIZE_T, 1, 0,
+                                                  __func__,
+                                                  "cltprm->snind->dsize");
       cltprm->snind->size  = cltprm->snind->dsize[0]=tablen;/* After dsize */
-      cltprm->snind->array = gal_data_malloc_array(cltprm->snind->type,
-                                                   tablen, __func__,
+      cltprm->snind->array = gal_pointer_allocate(cltprm->snind->type,
+                                                   tablen, 0, __func__,
                                                    "cltprm->snind->array");
     }
   else cltprm->snind=NULL;
@@ -571,10 +572,10 @@ clumps_find_make_sn_table(void *in_prm)
   uint8_t *binary=p->binary->array;
   struct clumps_thread_params cltprm;
   size_t i, j, c, ind, tind, num, numsky, *indarr;
-  size_t *scoord=gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim, __func__,
-                                       "scoord");
-  size_t *icoord=gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim, __func__,
-                                       "icoord");
+  size_t *scoord=gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__,
+                                      "scoord");
+  size_t *icoord=gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__,
+                                      "icoord");
 
 
   /* Initialize the parameters for this thread. */
@@ -645,9 +646,8 @@ clumps_find_make_sn_table(void *in_prm)
              rivers and not include them in the list of indexs to set
              clumps. To do that, we need this tile's starting
              coordinates. */
-          gal_dimension_index_to_coord(gal_data_num_between(p->clabel->array,
-                                                            tile->array,
-                                                            p->clabel->type),
+          gal_dimension_index_to_coord(gal_pointer_num_between(
+                           p->clabel->array, tile->array, p->clabel->type),
                                        ndim, dsize, scoord);
 
 
@@ -692,8 +692,8 @@ clumps_find_make_sn_table(void *in_prm)
                     {
                       if(cltprm.id==282) *i+=2;
                   */
-                      indarr[c++]=gal_data_num_between(p->clabel->array, i,
-                                                       p->clabel->type);
+                      indarr[c++]=gal_pointer_num_between(p->clabel->array,
+                                                          i, p->clabel->type);
                   /*
                     }
                   else
@@ -1075,9 +1075,9 @@ clumps_det_keep_true_relabel(struct clumps_thread_params 
*cltprm)
   if(cltprm->sn)
     {
       /* Allocate the necessary arrays. */
-      newlabs=gal_data_malloc_array(GAL_TYPE_INT32,
-                                    cltprm->numinitclumps+1, __func__,
-                                    "newlabs");
+      newlabs=gal_pointer_allocate(GAL_TYPE_INT32,
+                                   cltprm->numinitclumps+1, 0, __func__,
+                                   "newlabs");
       dinc=gal_dimension_increment(ndim, dsize);
 
       /* Initialize the new labels with GAL_LABEL_INIT (so the diffuse area
diff --git a/bin/segment/segment.c b/bin/segment/segment.c
index 1cf0154..8dca3ec 100644
--- a/bin/segment/segment.c
+++ b/bin/segment/segment.c
@@ -34,6 +34,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/label.h>
 #include <gnuastro/binary.h>
 #include <gnuastro/threads.h>
+#include <gnuastro/pointer.h>
 #include <gnuastro/convolve.h>
 #include <gnuastro/dimension.h>
 #include <gnuastro/statistics.h>
@@ -229,7 +230,7 @@ segment_relab_to_objects(struct clumps_thread_params 
*cltprm)
   size_t *dinc=gal_dimension_increment(ndim, dsize);
   size_t *s, *sf, i, j, ii, rpnum, *nums=nums_d->array;
   double ave, rpsum, c=sqrt(1/p->cpscorr), *sums=sums_d->array;
-  int32_t *ngblabs=gal_data_malloc_array(GAL_TYPE_UINT32, nngb, __func__,
+  int32_t *ngblabs=gal_pointer_allocate(GAL_TYPE_UINT32, nngb, 0, __func__,
                                          "ngblabs");
 
 
@@ -421,10 +422,10 @@ segment_relab_clumps_in_objects(struct 
clumps_thread_params *cltprm)
   int32_t *clumptoobj=cltprm->clumptoobj->array;
   int32_t *clabel=cltprm->clprm->p->clabel->array;
   size_t i, *s=cltprm->indexs->array, *sf=s+cltprm->indexs->size;
-  size_t *nclumpsinobj=gal_data_calloc_array(GAL_TYPE_SIZE_T, numobjects+1,
-                                             __func__, "nclumpsinobj");
-  int32_t *newlabs=gal_data_calloc_array(GAL_TYPE_UINT32, numtrueclumps+1,
-                                         __func__, "newlabs");
+  size_t *nclumpsinobj=gal_pointer_allocate(GAL_TYPE_SIZE_T, numobjects+1,
+                                             1, __func__, "nclumpsinobj");
+  int32_t *newlabs=gal_pointer_allocate(GAL_TYPE_UINT32, numtrueclumps+1,
+                                         1, __func__, "newlabs");
 
   /* Fill both arrays. */
   for(i=1;i<numtrueclumps+1;++i)
diff --git a/bin/segment/ui.c b/bin/segment/ui.c
index 8231dc3..228190b 100644
--- a/bin/segment/ui.c
+++ b/bin/segment/ui.c
@@ -417,7 +417,7 @@ ui_prepare_inputs(struct segmentparams *p)
       p->conv->wcs=gal_wcs_copy(p->input->wcs);
 
       /* Make sure it is the same size as the input. */
-      if( gal_data_dsize_is_different(p->input, p->conv) )
+      if( gal_dimension_is_different(p->input, p->conv) )
         error(EXIT_FAILURE, 0, "%s (hdu %s), given to `--convolved' and "
               "`--chdu', is not the same size as the input (%s, hdu: %s)",
               p->convolvedname, p->chdu, p->inputname, p->cp.hdu);
@@ -432,7 +432,7 @@ ui_prepare_inputs(struct segmentparams *p)
       /* Read the dataset into memory. */
       p->olabel = gal_array_read_one_ch(p->useddetectionname, p->dhdu,
                                         p->cp.minmapsize);
-      if( gal_data_dsize_is_different(p->input, p->olabel) )
+      if( gal_dimension_is_different(p->input, p->olabel) )
         error(EXIT_FAILURE, 0, "`%s' (hdu: %s) and `%s' (hdu: %s) have a"
               "different dimension/size", p->useddetectionname, p->dhdu,
               p->inputname, p->cp.hdu);
@@ -614,7 +614,7 @@ static void
 ui_check_size(gal_data_t *base, gal_data_t *comp, size_t numtiles,
               char *bname, char *bhdu, char *cname, char *chdu)
 {
-  if( gal_data_dsize_is_different(base, comp) && numtiles!=comp->size)
+  if( gal_dimension_is_different(base, comp) && numtiles!=comp->size )
     error(EXIT_FAILURE, 0, "%s (hdu: %s): doesn't have the right size "
           "(%zu elements or pixels).\n\n"
           "It must either be the same size as `%s' (hdu: `%s'), or "
@@ -646,7 +646,7 @@ ui_subtract_sky(gal_data_t *in, gal_data_t *sky,
   float *s, *f, *ff, *skyarr=sky->array;
 
   /* It is the same size as the input or a single value. */
-  if( gal_data_dsize_is_different(in, sky)==0 || sky->size==1)
+  if( gal_dimension_is_different(in, sky)==0 || sky->size==1)
     {
       s=sky->array;
       ff=(f=in->array)+in->size;
diff --git a/bin/statistics/sky.c b/bin/statistics/sky.c
index f50359d..8b5bc62 100644
--- a/bin/statistics/sky.c
+++ b/bin/statistics/sky.c
@@ -31,6 +31,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/fits.h>
 #include <gnuastro/qsort.h>
 #include <gnuastro/blank.h>
+#include <gnuastro/pointer.h>
 #include <gnuastro/threads.h>
 #include <gnuastro/convolve.h>
 #include <gnuastro/statistics.h>
@@ -95,19 +96,19 @@ sky_on_thread(void *in_prm)
           /* Put the mean and its standard deviation into the respective
              place for this tile. */
           sigmaclip=gal_data_copy_to_new_type_free(sigmaclip, stype);
-          memcpy(gal_data_ptr_increment(p->sky_t->array, tind, stype),
-                 gal_data_ptr_increment(sigmaclip->array, 2, stype), twidth);
-          memcpy(gal_data_ptr_increment(p->std_t->array, tind, stype),
-                 gal_data_ptr_increment(sigmaclip->array, 3, stype), twidth);
+          memcpy(gal_pointer_increment(p->sky_t->array, tind, stype),
+                 gal_pointer_increment(sigmaclip->array, 2, stype), twidth);
+          memcpy(gal_pointer_increment(p->std_t->array, tind, stype),
+                 gal_pointer_increment(sigmaclip->array, 3, stype), twidth);
 
           /* Clean up. */
           gal_data_free(sigmaclip);
         }
       else
         {
-          gal_blank_write(gal_data_ptr_increment(p->sky_t->array, tind,
+          gal_blank_write(gal_pointer_increment(p->sky_t->array, tind,
                                                  stype), stype);
-          gal_blank_write(gal_data_ptr_increment(p->std_t->array, tind,
+          gal_blank_write(gal_pointer_increment(p->std_t->array, tind,
                                                  stype), stype);
         }
 
diff --git a/bin/statistics/statistics.c b/bin/statistics/statistics.c
index 1d10f8e..3ecdf48 100644
--- a/bin/statistics/statistics.c
+++ b/bin/statistics/statistics.c
@@ -33,6 +33,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/fits.h>
 #include <gnuastro/tile.h>
 #include <gnuastro/blank.h>
+#include <gnuastro/pointer.h>
 #include <gnuastro/arithmetic.h>
 #include <gnuastro/statistics.h>
 #include <gnuastro/interpolate.h>
@@ -61,7 +62,7 @@ statistics_pull_out_element(gal_data_t *input, size_t index)
   gal_data_t *out=gal_data_alloc(NULL, input->type, 1, &dsize,
                                  NULL, 1, -1, NULL, NULL, NULL);
   memcpy( out->array,
-          gal_data_ptr_increment(input->array, index, input->type),
+          gal_pointer_increment(input->array, index, input->type),
           gal_type_sizeof(input->type) );
   return out;
 }
@@ -384,7 +385,7 @@ statistics_on_tile(struct statisticsparams *p)
 
           /* 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),
+          memcpy(gal_pointer_increment(values->array, tind++, values->type),
                  tmp->array, gal_type_sizeof(type));
           gal_data_free(tmp);
         }
diff --git a/bin/warp/ui.c b/bin/warp/ui.c
index 632715a..1e052ad 100644
--- a/bin/warp/ui.c
+++ b/bin/warp/ui.c
@@ -33,6 +33,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/array.h>
 #include <gnuastro/table.h>
 #include <gnuastro/threads.h>
+#include <gnuastro/pointer.h>
 
 #include <gnuastro-internal/timing.h>
 #include <gnuastro-internal/options.h>
@@ -413,7 +414,7 @@ ui_matrix_prepare_raw(struct warpparams *p)
   if(p->matrix->size==4)
     {
       /* Allocate the final matrix. */
-      final=gal_data_malloc_array(GAL_TYPE_FLOAT64, 9, __func__, "final");
+      final=gal_pointer_allocate(GAL_TYPE_FLOAT64, 9, 0, __func__, "final");
 
       /* Fill in the final 3x3 matrix from the 2x2 matrix. */
       final[0]=in[0];    final[1]=in[1];   final[2]=0.0f;
@@ -429,8 +430,8 @@ ui_matrix_prepare_raw(struct warpparams *p)
   /* Correct the dimensional information, because the matrix was read as a
      single dimensional list of numbers. */
   free(p->matrix->dsize);
-  dsize=p->matrix->dsize=gal_data_malloc_array(GAL_TYPE_SIZE_T, 2, __func__,
-                                               "dsize");
+  dsize=p->matrix->dsize=gal_pointer_allocate(GAL_TYPE_SIZE_T, 2, 0,
+                                               __func__, "dsize");
   dsize[0]=dsize[1]=3;
   p->matrix->ndim=2;
 }
@@ -756,8 +757,8 @@ ui_matrix_finalize(struct warpparams *p)
      efficient processing. */
 
    /* Make the inverse matrix: */
-  inv=p->inverse=gal_data_malloc_array(GAL_TYPE_FLOAT64, 9, __func__,
-                                       "p->inverse");
+  inv=p->inverse=gal_pointer_allocate(GAL_TYPE_FLOAT64, 9, 0, __func__,
+                                      "p->inverse");
   inv[0] = d[4]*d[8] - d[5]*d[7];
   inv[1] = d[2]*d[7] - d[1]*d[8];
   inv[2] = d[1]*d[5] - d[2]*d[4];
diff --git a/doc/announce-acknowledge.txt b/doc/announce-acknowledge.txt
index c962921..c74b8fa 100644
--- a/doc/announce-acknowledge.txt
+++ b/doc/announce-acknowledge.txt
@@ -4,6 +4,7 @@ Leindert Boogaard
 Nima Dehdilani
 Antonio Diaz Diaz
 Lee Kelvin
+Brandon Kelly
 Guillaume Mahler
 Bertrand Pain
 Ole Streicher
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index b5da0f1..ae3b209 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -563,6 +563,7 @@ Gnuastro library
 * Configuration information::   General information about library config.
 * Multithreaded programming::   Tools for easy multi-threaded operations.
 * Library data types::          Definitions and functions for types.
+* Pointers::                    Wrappers for easy working with 
pointers.@strong{}
 * Library blank values::        Blank values and functions to deal with them.
 * Library data container::      General data container in Gnuastro.
 * Dimensions::                  Dealing with coordinates and dimensions.
@@ -595,7 +596,7 @@ Multithreaded programming (@file{threads.h})
 Data container (@file{data.h})
 
 * Generic data container::      Definition of Gnuastro's generic container.
-* Dataset size and allocation::  Functions for size and allocation.
+* Dataset allocation::          Allocate, initialize and free a dataset.
 * Arrays of datasets::          Functions to help with array of datasets.
 * Copying datasets::            Functions to copy a dataset to a new one.
 
@@ -4321,10 +4322,10 @@ added or changed, also some updates that are not yet 
officially released
 might be in it.
 
 We use Git for the version control of Gnuastro. For those who are not
-familiar with it, we recommend the Pro
-Git@footnote{@url{https://progit.org/}} book. The whole book is publicly
-available for online reading and downloading and does a wonderful job at
-explaining the concepts and best practices.
+familiar with it, we recommend the @url{https://git-scm.com/book/en, ProGit
+book}. The whole book is publicly available for online reading and
+downloading and does a wonderful job at explaining the concepts and best
+practices.
 
 Let's assume you want to keep Gnuastro in the @file{TOPGNUASTRO} directory
 (can be any directory, change the value below). The full version controlled
@@ -4350,13 +4351,14 @@ and @file{README}) or mainly written in small-caps (for 
example
 @file{configure.ac} and @file{Makefile.am}). The former are
 non-programming, standard writing for human readers containing high-level
 information about the whole package. The latter are instructions to
-customize the GNU build system for Gnuastro.
+customize the GNU build system for Gnuastro. For more on Gnuastro's source
+code structure, please see @ref{Developing}. We won't go any deeper here.
 
 The cloned Gnuastro source cannot immediately be configured, compiled, or
 installed since it only contains hand-written files, not automatically
 generated or imported files which do all the hard work of the build
 process. See @ref{Bootstrapping} for the process of generating and
-importing those files (it is very easy!). Once you have bootstrapped
+importing those files (its not too hard!). Once you have bootstrapped
 Gnuastro, you can run the standard procedures (in @ref{Quick start}). Very
 soon after you have cloned it, Gnuastro's main @file{master} branch will be
 updated on the main repository (since the developers are actively working
@@ -13753,15 +13755,17 @@ one for the Sky standard deviation).
 Once instrumental signatures are removed from the raw data (image) in the
 initial reduction process (see @ref{Data manipulation}). You are naturally
 eager to start answering the scientific questions that motivated the data
-collection. However, the raw dataset/image is just an array of
-values/pixels, that is all! These raw values cannot directly be used to
-answer your scientific questions (for example ``how many galaxies are there
-in the image?''). Therefore, the first high-level step will be to classify,
-or label, the dataset elements (pixels) into two classes: 1) noise, where
-random effects are the major contributor to the value, and 2) signal, where
-non-random factors (for example light from a distant galaxy) are
-present. This classification of a dataset is formally known as
-@emph{detection}.
+collection in the first place. However, the raw dataset/image is just an
+array of values/pixels, that is all! These raw values cannot directly be
+used to answer your scientific questions: for example ``how many galaxies
+are there in the image?''.
+
+The first high-level step in your analysis of your targets will thus be to
+classify, or label, the dataset elements (pixels) into two classes: 1)
+noise, where random effects are the major contributor to the value, and 2)
+signal, where non-random factors (for example light from a distant galaxy)
+are present. This classification of the elements in a dataset is formally
+known as @emph{detection}.
 
 In an observational/experimental dataset, signal is always burried in
 noise: only mock/simulated datasets are free of noise. Therefore detection,
@@ -13771,62 +13775,69 @@ on them. Detection is thus the most important step of 
any analysis and is
 not trivial. In particular, the most scientifically interesting
 astronomical targets are faint, can have a large variety of morphologies,
 along with a large distribution in brightness and size. Therefore when
-noise is significant, detection is the uniquely decisive step of your
-scientific result.
+noise is significant, proper detection of your targets is a uniquely
+decisive step in your final scientific analysis/result.
 
 @cindex Erosion
 NoiseChisel is Gnuastro's program for detection of targets that don't have
 a sharp border (almost all astronomical objects). When the targets have a
-sharp edges/border, a simple threshold is enough to separate them from
-noise and each other (if they are not touching). For detection of such
-targets, you can use Gnuastro's Arithmetic program in a command like below
-(assuming the threshold is @code{100}, see @ref{Arithmetic}):
+sharp edges/border (for example cells in biological imaging), a simple
+threshold is enough to separate them from noise and each other (if they are
+not touching). To detect such sharp-edged targets, you can use Gnuastro's
+Arithmetic program in a command like below (assuming the threshold is
+@code{100}, see @ref{Arithmetic}):
 
 @example
 $ astarithmetic in.fits 100 gt 2 connected-components
 @end example
 
-NoiseChisel uses a new noise-based paradigm for detection of very exteded
-and diffuse targets that are drowned deeply in the ocean of noise. It was
-initially introduced in @url{https://arxiv.org/abs/1505.01664, Akhlaghi and
-Ichikawa [2015]}. The name of NoiseChisel is derived from the first thing
-it does after thresholding the dataset: to erode it. In mathematical
-morphology, erosion on pixels can be pictured as carving-off boundary
-pixels. Hence, what NoiseChisel does is similar to what a wood chisel or
-stone chisel do. It is just not a hardware, but a software. In fact,
-looking at it as a chisel and your dataset as a solid cube of rock will
-greatly help in effectively understanding and optimally using it: with
-NoiseChisel you literally carve your targets out of the noise. Try running
-it with the @option{--checkdetection} option to see each step of the
-carving process on your input dataset.
+Since almost no astronomical target has such sharp edges, we need a more
+advanced detection methodology. NoiseChisel uses a new noise-based paradigm
+for detection of very exteded and diffuse targets that are drowned deeply
+in the ocean of noise. It was initially introduced in
+@url{https://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa [2015]}. The
+name of NoiseChisel is derived from the first thing it does after
+thresholding the dataset: to erode it. In mathematical morphology, erosion
+on pixels can be pictured as carving-off boundary pixels. Hence, what
+NoiseChisel does is similar to what a wood chisel or stone chisel do. It is
+just not a hardware, but a software. In fact, looking at it as a chisel and
+your dataset as a solid cube of rock will greatly help in effectively
+understanding and optimally using it: with NoiseChisel you literally carve
+your targets out of the noise. Try running it with the
+@option{--checkdetection} option to see each step of the carving process on
+your input dataset.
 
 @cindex Segmentation
 NoiseChisel's primary output is a binary detection map with the same size
 as the input but only with two values: 0 and 1. Pixels that don't harbor
 any detected signal (noise) are given a label (or value) of zero and those
-with a value of 1 have been identified as hosting signal. Segmentation is
-the process of classifing the signal into higher-level constructs, for
-example if you have two separate galaxies in one image, after segmentation,
-they will get separate labels. NoiseChisel is only focused on detection
-(separating signal from noise), to segment the signal (into separate
-galaxies for example), Gnuastro has a separate specialized program
-@ref{Segment}. NoiseChisel's output can be directly/readily fed into
-Segment.
+with a value of 1 have been identified as hosting signal.
+
+Segmentation is the process of classifing the signal into higher-level
+constructs. For example if you have two separate galaxies in one image, by
+default NoiseChisel will give a value of 1 to the pixels of both, but after
+segmentation, the pixels in each will get separate labels. NoiseChisel is
+only focused on detection (separating signal from noise), to @emph{segment}
+the signal (into separate galaxies for example), Gnuastro has a separate
+specialized program @ref{Segment}. NoiseChisel's output can be
+directly/readily fed into Segment.
 
 For more on NoiseChisel's output format and its benefits (especially in
 conjunction with @ref{Segment} and later @ref{MakeCatalog}), please see
 @url{https://arxiv.org/abs/1611.06387, Akhlaghi [2016]}. Just note that
 when that paper was published, Segment was not yet spinned-off, and
-NoiseChisel done both detection and segmentation. The output is designed to
-be generic enough to be easily used in any higher-level analysis. If your
-targets are not touching after running NoiseChisel and you aren't
-interested in their sub-structure, you don't need the Segment program at
-all. You can ask NoiseChisel to find the connected pixels in the output
-with the @option{--label} option. In this case, the output won't be a
-binary image any more, the signal will have counters/labels starting from 1
-for each connected group of pixels. You can then directly feed
-NoiseChisel's output into MakeCatalog for measurements over the detections
-and the production of a catalog (see @ref{MakeCatalog}).
+NoiseChisel done both detection and segmentation.
+
+NoiseChisel's output is designed to be generic enough to be easily used in
+any higher-level analysis. If your targets are not touching after running
+NoiseChisel and you aren't interested in their sub-structure, you don't
+need the Segment program at all. You can ask NoiseChisel to find the
+connected pixels in the output with the @option{--label} option. In this
+case, the output won't be a binary image any more, the signal will have
+counters/labels starting from 1 for each connected group of pixels. You can
+then directly feed NoiseChisel's output into MakeCatalog for measurements
+over the detections and the production of a catalog (see
+@ref{MakeCatalog}).
 
 Thanks to the published papers mentioned above, there is no need to provide
 a more complete introduction to NoiseChisel in this book. However,
@@ -13836,15 +13847,15 @@ evolved/changed. The changes since publication are 
documented in
 astnoisechisel}, the details of running NoiseChisel and its options are
 discussed.
 
-As discussed above, detection is the single most important step for your
+As discussed above, detection is one of the most important steps for your
 scientific result. It is therefore very important to obtain a good
 understanding of NoiseChisel (and afterwards @ref{Segment} and
 @ref{MakeCatalog}). We thus strongly recommend that after reading the
 papers above and the respective sections of Gnuastro's book, you play a
-little with the settings (in the order presented in @ref{Invoking
-astnoisechisel}) on a dataset you are familiar with and inspect all the
-check images (options starting with @option{--check}) to see the effect of
-each parameter.
+little with the settings (in the order presented in the paper and
+@ref{Invoking astnoisechisel}) on a dataset you are familiar with and
+inspect all the check images (options starting with @option{--check}) to
+see the effect of each parameter.
 
 @ref{General program usage tutorial} is also a good place to get a feeling
 of the modular principle behind Gnuastro's programs and how they are built
@@ -13856,10 +13867,10 @@ that tutorial for optimal usage of NoiseChisel in 
conjunction with all the
 other Gnuastro programs.
 
 In @ref{NoiseChisel changes after publication}, we'll review the changes in
-NoiseChisel since the publication of @url{https://arxiv.org/abs/1611.06387,
-Akhlaghi [2016]}. We will then review NoiseChisel's input, detection, and
-output options in @ref{NoiseChisel input}, @ref{Detection options}, and
-@ref{NoiseChisel output}.
+NoiseChisel since the publication of @url{https://arxiv.org/abs/1505.01664,
+Akhlaghi and Ichikawa [2015]}. We will then review NoiseChisel's input,
+detection, and output options in @ref{NoiseChisel input}, @ref{Detection
+options}, and @ref{NoiseChisel output}.
 
 @menu
 * NoiseChisel changes after publication::  Changes to the software after 
publication.
@@ -13871,10 +13882,10 @@ output options in @ref{NoiseChisel input}, 
@ref{Detection options}, and
 
 NoiseChisel was initially introduced in
 @url{https://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa [2015]}. It is
-thus strongly recommended to read it for a good understanding of what it
-does and how each parameter influences the output. To help in understanding
-how it works, that paper has a large number of figures showing every step
-on multiple mock and real examples.
+thus strongly recommended to read this paper for a good understanding of
+what it does and how each parameter influences the output. To help in
+understanding how it works, that paper has a large number of figures
+showing every step on multiple mock and real examples.
 
 However, the paper cannot be updated anymore, but NoiseChisel has evolved
 (and will continue to do so): better algorithms or steps have been found,
@@ -19878,57 +19889,71 @@ this function is actually defined here, @code{U} says 
that it is undefined
 here). Try opening the @file{.c} file to check some of these functions for
 your self. Run @command{info nm} for more information.
 
-@cindex Static linking
-@cindex Linking: Static
+@cindex Linking
 To recap, the @emph{compiler} created the separate object files mentioned
 above for each @file{.c} file. The @emph{linker} will then combine all the
 symbols of the various object files (and libraries) into one program or
 library. In the case of Arithmetic (a program) the contents of the object
 files in @file{bin/arithmetic/} are copied (and re-ordered) into one final
-executable file which we can run from the operating system. When the
-symbols (computer-readable function definitions in most cases) are copied
-into the output like this, we call the process @emph{static} linking. Let's
-have a closer look at static linking: we'll assume you have installed
+executable file which we can run from the operating system.
+
+@cindex Static linking
+@cindex Linking: Static
+@cindex Dynamic linking
+@cindex Linking: Dynamic
+There are two ways to @emph{link} all the necessary symbols: static and
+dynamic/shared. When the symbols (computer-readable function definitions in
+most cases) are copied into the output, it is called @emph{static}
+linking. When the symbols are kept in their original file and only a
+reference to them is kept in the executable, it is called @emph{dynamic},
+or @emph{shared} linking.
+
+Let's have a closer look at the executable to understand this better: we'll
+assume you have built Gnuastro without any customization and installed
 Gnuastro into the default @file{/usr/local/} directory (see
 @ref{Installation directory}). If you tried the @command{nm} command on one
 of Arithmetic's object files above, then with the command below you can
-confirm that all the functions that were defined in the object files (had a
-@code{T} in the second column) are also defined in the @file{astarithmetic}
-executable:
+confirm that all the functions that were defined in the object file above
+(had a @code{T} in the second column) are also defined in the
+@file{astarithmetic} executable:
 
 @example
 $ nm /usr/local/bin/astarithmetic
 @end example
 
 @noindent
-But you will notice that there are still many undefined symbols (have a
-@code{U} in the second column). One class of such functions are Gnuastro's
-own library functions that start with `@code{gal_}':
+These symbols/function have been statically linked (copied) in the final
+executable. But you will notice that there are still many undefined symbols
+in the executable (those with a @code{U} in the second column). One class
+of such functions are Gnuastro's own library functions that start with
+`@code{gal_}':
 
 @example
 $ nm /usr/local/bin/astarithmetic | grep gal_
 @end example
 
+@cindex Plugin
 @cindex GNU Libtool
 @cindex Shared library
 @cindex Library: shared
 @cindex Dynamic linking
 @cindex Linking: dynamic
-These undefined symbols (functions) will be linked to the executable
-every time you run arithmetic. Therefore they are known as dynamically
-@emph{linked} libraries @footnote{Do not confuse dynamically @emph{linked}
-libraries with dynamically @emph{loaded} libraries. The former (that is
-discussed here) are only loaded once at the program startup. However, the
-latter can be loaded anytime during the program's execution, they are also
-known as plugins.}. When the functions of a library need to be dynamically
-linked, the library is known as a shared library. As we saw above, static
-linking is done when the executable is being built. However, when a library
-is linked dynamically, its symbols are only checked with the available
-libraries at build time: they are not actually copied into the
-executable. Every time you run the program, the linker will be activated and
-will try to link the program to the installed library before it starts. If
-you want all the libraries to be statically linked to the executables, you
-have to tell Libtool (which Gnuastro uses for the linking) to disable
+These undefined symbols (functions) are present in another file and will be
+linked to the Arithmetic program every time you run it. Therefore they are
+known as dynamically @emph{linked} libraries @footnote{Do not confuse
+dynamically @emph{linked} libraries with dynamically @emph{loaded}
+libraries. The former (that is discussed here) are only loaded once at the
+program startup. However, the latter can be loaded anytime during the
+program's execution, they are also known as plugins.}. As we saw above,
+static linking is done when the executable is being built. However, when a
+program is dynamically linked to a library, at build-time, the library's
+symbols are only checked with the available libraries: they are not
+actually copied into the program's executable. Every time you run the
+program, the (dynamic) linker will be activated and will try to link the
+program to the installed library before the program starts.
+
+If you want all the libraries to be statically linked to the executables,
+you have to tell Libtool (which Gnuastro uses for the linking) to disable
 shared libraries at configure time@footnote{Libtool is very common and is
 commonly used. Therefore, you can use this option to configure on most
 programs using the GNU build system if you want static linking.}:
@@ -19938,12 +19963,12 @@ $ configure --disable-shared
 @end example
 
 @noindent
-Try configuring, statically building and installing Gnuastro with the
-command above. Then check the @code{gal_} symbols in the installed
-Arithmetic executable like before. You will see that they are actually
-copied this time (have a @code{T} in the second column). If the second
-column doesn't convince you, look at the executable file size with the
-following command:
+Try configuring Gnuastro with the command above, then build and install it
+(as described in @ref{Quick start}). Afterwards, check the @code{gal_}
+symbols in the installed Arithmetic executable like before. You will see
+that they are actually copied this time (have a @code{T} in the second
+column). If the second column doesn't convince you, look at the executable
+file size with the following command:
 
 @example
 $ ls -lh /usr/local/bin/astarithmetic
@@ -19952,19 +19977,20 @@ $ ls -lh /usr/local/bin/astarithmetic
 @noindent
 It should be around 4.2 Megabytes with this static linking. If you
 configure and build Gnuastro again with shared libraries enabled (which is
-the default), you will notice that it is roughly 100 Kilobytes! This huge
-difference would have been very significant in the old days, but with the
-roughly Terabyte storage drives commonly in use today, it is
+the default), you will notice that it is roughly 100 Kilobytes!
+
+This huge difference would have been very significant in the old days, but
+with the roughly Terabyte storage drives commonly in use today, it is
 negligible. Fortunately, output file size is not the only benefit of
 dynamic linking: since it links to the libraries at run-time (rather than
 build-time), you don't have to re-build a higher-level program or library
 when an update comes for one of the lower-level libraries it depends
 on. You just install the new low-level library and it will automatically be
-used next time in your higher-level tools. To be fair, this also creates a
-few complications@footnote{Both of these can be avoided by joining the
-mailing lists of the lower-level libraries and checking the changes in
-newer versions before installing them. Updates that result in such
-behaviors are generally heavily emphasized in the release notes.}:
+used/linked next time in the programs that use it. To be fair, this also
+creates a few complications@footnote{Both of these can be avoided by
+joining the mailing lists of the lower-level libraries and checking the
+changes in newer versions before installing them. Updates that result in
+such behaviors are generally heavily emphasized in the release notes.}:
 
 @itemize
 @item
@@ -19978,7 +20004,7 @@ you need to re-build your higher-level program or 
library.
 
 @cindex GNU C library
 To see a list of all the shared libraries that are needed for a program or
-a shared library to run, you can use the GNU C library's
+a shared library to run, you can use GNU C library's
 @command{ldd}@footnote{If your operating system is not using the GNU C
 library, you might need another tool.} program, for example:
 
@@ -19986,46 +20012,61 @@ library, you might need another tool.} program, for 
example:
 $ ldd /usr/local/bin/astarithmetic
 @end example
 
-Library file names start with a @file{lib} and end with a suffix depending
-on their type as described below. In between these two is the name of the
-library, for example @file{libgnuastro.a} (Gnuastro's static library) and
-@file{libgsl.so.0.0.0} (GSL's shared library).
+Library file names (in their installation directory) start with a
+@file{lib} and their ending (suffix) shows if they are static (@file{.a})
+or dynamic (@file{.so}), as described below. The name of the library is in
+the middle of these two, for example @file{libgsl.a} or
+@file{libgnuastro.a} (GSL and Gnuastro's static libraries), and
+@file{libgsl.so.23.0.0} or @file{libgnuastro.so.4.0.0} (GSL and Gnuastro's
+shared library, the numbers may be different).
 
 @itemize
 @item
-A static library is known as an archive file and has a @file{.a} suffix. A
-static library is not an executable file.
+A static library is known as an archive file and has the @file{.a}
+suffix. A static library is not an executable file.
 
 @item
 @cindex Shared library versioning
 @cindex Versioning: Shared library
-A shared library ends with a @file{.so.X.Y.Z} suffix and is executable. The
-three numbers in the prefix are the version of the shared library. Shared
-library versions are defined to allow multiple versions of a shared library
-simultaneously on a system and to help detect possible updates in the
-library and programs that depend on it by the linker. It is very important
-to mention that this version number is different from from the software
-version number (see @ref{Version numbering}), so do not confuse the
-two. See the ``Library interface versions'' chapter of GNU Libtool for
-more.
+A shared library ends with the @file{.so.X.Y.Z} suffix and is
+executable. The three numbers in the suffix, describe the version of the
+shared library. Shared library versions are defined to allow multiple
+versions of a shared library simultaneously on a system and to help detect
+possible updates in the library and programs that depend on it by the
+linker.
+
+It is very important to mention that this version number is different from
+the software version number (see @ref{Version numbering}), so do not
+confuse the two. See the ``Library interface versions'' chapter of GNU
+Libtool for more.
 
 For each shared library, we also have two symbolic links ending with
 @file{.so.X} and @file{.so}. They are automatically set by the installer,
 but you can change them (point them to another version of the library) when
-you have multiple versions on your system.
+you have multiple versions of a library on your system.
 
 @end itemize
 
 @cindex GNU Libtool
-For those libraries that use GNU Libtool (including Gnuastro and its
-dependencies), both static and dynamic libraries are built and installed in
-the @file{prefix/lib/} directory (see @ref{Installation directory}). In
-this way other programs can make which ever kind of link that they want.
+Libraries that are built with GNU Libtool (including Gnuastro and its
+dependencies), build both static and dynamic libraries by default and
+install them in @file{prefix/lib/} directory (for more on @file{prefix},
+see @ref{Installation directory}). In this way, programs depending on the
+libraries can link with them however they prefer. See the contents of
+@file{/usr/local/lib} with the command below to see both the static and
+shared libraries available there, along with their executable nature and
+the symbolic links:
+
+@example
+$ ls -l /usr/local/lib/
+@end example
 
 To link with a library, the linker needs to know where to find the
 library. @emph{At compilation time}, these locations can be passed to the
 linker with two separate options (see @ref{Summary and example on
-libraries} for an example) as described below.
+libraries} for an example) as described below. You can see these options
+and their usage in practice while building Gnuastro (after running
+@command{make}):
 
 @table @option
 @item -L DIR
@@ -20033,21 +20074,22 @@ Will tell the linker to look into @file{DIR} for the 
libraries. For example
 @file{-L/usr/local/lib}, or @file{-L/home/yourname/.local/lib}. You can
 make multiple calls to this option, so the linker looks into several
 directories at compilation time. Note that the space between @key{L} and
-the directory is optional and commonly not used.
+the directory is optional and commonly ignored (written as @option{-LDIR}).
 
 @item -lLIBRARY
-Specify the unique name of a library to be linked. As discussed above,
-library file names have fixed parts which must not be given to this
+Specify the unique library identifier/name (not containing directory or
+shared/dynamic nature) to be linked with the executable. As discussed
+above, library file names have fixed parts which must not be given to this
 option. So @option{-lgsl} will guide the linker to either look for
 @file{libgsl.a} or @file{libgsl.so} (depending on the type of linking it is
 suppose to do). You can link many libraries by repeated calls to this
 option.
 
-@strong{Very important: } The place of this option on the command line
-matters. This is often a source of confusion for beginners, so let's assume
-you have asked the linker to link with library A using this option. As soon
-as the linker confronts this option, it looks into the list of the
-undefined symbols it has found until that point and does a search in
+@strong{Very important: } The place of this option on the compiler's
+command matters. This is often a source of confusion for beginners, so
+let's assume you have asked the linker to link with library A using this
+option. As soon as the linker confronts this option, it looks into the list
+of the undefined symbols it has found until that point and does a search in
 library A for any of those symbols. If any pending undefined symbol is
 found in library A, it is used. After the search in undefined symbols is
 complete, the contents of library A are completely discarded from the
@@ -20055,22 +20097,30 @@ linker's memory. Therefore, if a later object file or 
library uses an
 unlinked symbol in library A, the linker will abort after it has finished
 its search in all the input libraries or object files.
 
-As an example, Gnuastro's @code{gal_array_dlog10_array} function depends on
-the @code{log10} function of the C Math library (specified with
-@option{-lm}). So the proper way to link something that uses this function
-is @option{-lgnuastro -lm}. If instead, you give: @option{-lm -lgnuastro}
-the linker will complain and abort.
+As an example, Gnuastro's @code{gal_fits_img_read} function depends on the
+@code{fits_read_pix} function of CFITSIO (specified with
+@option{-lcfitsio}, which in turn depends on the cURL library, called with
+@option{-lcurl}). So the proper way to link something that uses this
+function is @option{-lgnuastro -lcfitsio -lcurl}. If instead, you give:
+@option{-lcfitsio -lgnuastro} the linker will complain and abort. To avoid
+such linking complexities when using Gnuastro's library, we recommend using
+@ref{BuildProgram}.
 
 @end table
 
 If you have compiled and linked your program with a dynamic library, then
-the linker needs to know the location of the libraries @emph{every time}
-the final program is run. For this purpose, the linker looks into the
-@code{LD_LIBRARY_PATH} environment variable. Therefore, if you don't get
-any errors when compiling/linking, but are unable to run your program
-because of a failure to find a library, it is because the linker hasn't
-found the dynamic library @emph{at run time}. See @ref{Installation
-directory} on how to add any directory to @code{LD_LIBRARY_PATH}.
+the dynamic linker also needs to know the location of the libraries after
+building the program: @emph{every time} the program is run
+afterwards. Therefore, it may happen that you don't get any errors when
+compiling/linking a program, but are unable to run your program because of
+a failure to find a library. This happens because the dynamic linker hasn't
+found the dynamic library @emph{at run time}.
+
+To find the dynamic libraries at run-time, the linker looks into the paths,
+or directories, in the @code{LD_LIBRARY_PATH} environment variable. For a
+discussion on environment variables, especially search paths like
+@code{LD_LIBRARY_PATH}, and how you can add new directories to them, see
+@ref{Installation directory}.
 
 
 
@@ -20425,6 +20475,7 @@ documentation will correspond to your installed version.
 * Configuration information::   General information about library config.
 * Multithreaded programming::   Tools for easy multi-threaded operations.
 * Library data types::          Definitions and functions for types.
+* Pointers::                    Wrappers for easy working with 
pointers.@strong{}
 * Library blank values::        Blank values and functions to deal with them.
 * Library data container::      General data container in Gnuastro.
 * Dimensions::                  Dealing with coordinates and dimensions.
@@ -20733,7 +20784,7 @@ program in @file{tests/lib/multithread.c} for a 
demonstration.
 
 @end deftypefun
 
-@node Library data types, Library blank values, Multithreaded programming, 
Gnuastro library
+@node Library data types, Pointers, Multithreaded programming, Gnuastro library
 @subsection Library data types (@file{type.h})
 
 Data in astronomy can have many types, numeric (numbers) and strings
@@ -21005,8 +21056,6 @@ if( gal_type_from_string(&out, string, 
GAL_TYPE_FLOAT32) )
 @end example
 @end deftypefun
 
-
-
 @deftypefun {void *} gal_type_string_to_number (char @code{*string}, uint8_t 
@code{*type})
 Read @code{string} into smallest type that can host the number, the
 allocated space for the number will be returned and the type of the number
@@ -21019,7 +21068,84 @@ the number of significant digits and determine if the 
given string is
 single or double precision as described in that section.
 @end deftypefun
 
-@node Library blank values, Library data container, Library data types, 
Gnuastro library
+@node Pointers, Library blank values, Library data types, Gnuastro library
+@subsection Pointers (@file{pointer.h})
+
+@cindex Pointers
+Pointers play an important role in the C programming language. As the name
+suggests, they @emph{point} to a byte in memory (like an address in a
+city). The C programming language gives you complete freedom in how to use
+the byte (and the bytes that follow it). Pointers are thus a very powerful
+feature of C. However, as the saying goes: ``With great power comes great
+responsability'', so they must be approached with care. The functions in
+this header are not very complex, they are just wrappers over some basic
+pointer functionality regarding pointer arithmetic and allocation (in
+memory or HDD/SSD).
+
+@deftypefun {void *} gal_pointer_increment (void @code{*pointer}, size_t 
@code{increment}, uint8_t @code{type})
+Return a pointer to an element that is @code{increment} elements ahead of
+@code{pointer}, assuming each element has type of @code{type}. For the type
+codes, see @ref{Library data types}.
+
+When working with the @code{array} elements of @code{gal_data_t}, we are
+actually dealing with @code{void *} pointers. However, pointer arithmetic
+doesn't apply to @code{void *}, because the system doesn't know how many
+bytes there are in each element to increment the pointer respectively. This
+function will use the given @code{type} to calculate where the incremented
+element is located in memory.
+@end deftypefun
+
+@deftypefun size_t gal_pointer_num_between (void @code{*earlier}, void 
@code{*later}, uint8_t @code{type})
+Return the number of elements (in the given @code{type}) between
+@code{earlier} and @code{later}. For the type codes, see @ref{Library data
+types}).
+@end deftypefun
+
+@deftypefun {void *} gal_pointer_allocate (uint8_t @code{type}, size_t 
@code{size}, int @code{clear}, const char @code{*funcname}, const char 
@code{*varname})
+Allocate an array of type @code{type} with @code{size} elements in RAM (for
+the type codes, see @ref{Library data types}). If @code{clear!=0}, then the
+allocated space is set to zero (cleared). This is effectively just a
+wrapper around C's @code{malloc} or @code{calloc} functions but takes
+Gnuastro's integer type codes and will also abort with a clear error if
+there the allocation was not successful.
+
+@cindex C99
+When space cannot be allocated, this function will abort the program with a
+message containing the reason for the failure. @code{funcname} (name of the
+function calling this function) and @code{varname} (name of variable that
+needs this space) will be used in this error message if they are not
+@code{NULL}. In most modern compilers, you can use the generic
+@code{__func__} variable for @code{funcname}. In this way, you don't have
+to manually copy and paste the function name or worry about it changing
+later (@code{__func__} was standardized in C99).
+@end deftypefun
+
+@deftypefun {void *} gal_pointer_allocate_mmap (size_t @code{size}, uint8_t 
@code{type}, int @code{clear}, char @code{**mmapname})
+Allocate the necessary space to keep @code{size} elements of type
+@code{type} in HDD/SSD (a file, not in RAM). for the type codes, see
+@ref{Library data types}. If @code{clear!=0}, then the allocated space will
+also be cleared. The allocation is done using C's @code{mmap} function. The
+name of the file containing the allocated space is an allocated string that
+will be put in @code{*mmapname}.
+
+Note that the kernel doesn't allow an infinite number of memory mappings to
+files. So it is not recommended to use this function with every
+allocation. The best case scenario to use this function is for large arrays
+that are very large and can fill up the RAM. Keep the smaller arrays in
+RAM, which is faster and can have a (theoretically) unlimited number of
+allocations.
+
+When you are done with the dataset and don't need it anymore, don't use
+@code{free} (the dataset isn't in RAM). Just delete the file (and the
+allocated space for the filename) with the commands below:
+
+@example
+remove(mmapname);
+free(mmapname);
+@end example
+@end deftypefun
+
+@node Library blank values, Library data container, Pointers, Gnuastro library
 @subsection Library blank values (@file{blank.h})
 When the position of an element in a dataset is important (for example a
 pixel in an image), a place-holder is necessary for the element if we don't
@@ -21209,12 +21335,12 @@ and some comments. To deal with any generic dataset, 
Gnuastro defines the
 
 @menu
 * Generic data container::      Definition of Gnuastro's generic container.
-* Dataset size and allocation::  Functions for size and allocation.
+* Dataset allocation::          Allocate, initialize and free a dataset.
 * Arrays of datasets::          Functions to help with array of datasets.
 * Copying datasets::            Functions to copy a dataset to a new one.
 @end menu
 
-@node Generic data container, Dataset size and allocation, Library data 
container, Library data container
+@node Generic data container, Dataset allocation, Library data container, 
Library data container
 @subsubsection Generic data container (@code{gal_data_t})
 
 To be able to deal with any dataset (various dimensions, numeric data
@@ -21477,62 +21603,14 @@ tiles that are created from this pointer.
 @end table
 
 
-@node Dataset size and allocation, Arrays of datasets, Generic data container, 
Library data container
-@subsubsection Dataset size and allocation
+@node Dataset allocation, Arrays of datasets, Generic data container, Library 
data container
+@subsubsection Dataset allocation
 
 Gnuastro's main data container was defined in @ref{Generic data container}.
 The functions listed in this section describe the most basic operations on
-@code{gal_data_t}: those related to the size, pointers, allocation and
-freeing. These functions are declared in @file{gnuastro/data.h} which is
-also visible from the function names (see @ref{Gnuastro library}).
-
-@deftypefun int gal_data_dsize_is_different (gal_data_t @code{*first}, 
gal_data_t @code{*second})
-Return @code{1} (one) if the two datasets don't have the same size along
-all dimensions. This function will also return @code{1} when the number of
-dimensions of the two datasets are different.
-@end deftypefun
-
-@deftypefun {void *} gal_data_ptr_increment (void @code{*pointer}, size_t 
@code{increment}, uint8_t @code{type})
-Return a pointer to an element that is @code{increment} elements ahead of
-@code{pointer}, assuming each element has type of @code{type} (for the type
-codes, see @ref{Library data types}).
-
-When working with the @code{array} elements of @code{gal_data_t}, we are
-actually dealing with @code{void *} pointers. However, pointer arithmetic
-doesn't apply to @code{void *}, because the system doesn't know how many
-bytes there are in each element to increment the pointer respectively. This
-function will use the given @code{type} to calculate where the incremented
-element is located in memory.
-@end deftypefun
-
-@deftypefun size_t gal_data_num_between (void @code{*earlier}, void 
@code{*later}, uint8_t @code{type})
-Return the number of elements (in the given @code{type}) between
-@code{earlier} and @code{later}. For the type codes, see @ref{Library data
-types}).
-@end deftypefun
-
-@deftypefun {void *} gal_data_malloc_array (uint8_t @code{type}, size_t 
@code{size}, const char @code{*funcname}, const char @code{*varname})
-Allocate an array of type @code{type} with @code{size} elements in RAM (for
-the type codes, see @ref{Library data types}). This is effectively just a
-wrapper around C's @code{malloc} function but takes Gnuastro's integer type
-codes and will also abort with an error if there the allocation was not
-successful.
-
-@cindex C99
-When space cannot be allocated, this function will abort the program with a
-message containing the reason for the failure. @code{funcname} (name of the
-function calling this function) and @code{varname} (name of variable that
-needs this space) will be used in this error message if they are not
-@code{NULL}. In most modern compilers, you can use the generic
-@code{__func__} variable for @code{funcname}. In this way, you don't have
-to manually copy and paste the function name or worry about it changing
-later (@code{__func__} was standardized in C99).
-@end deftypefun
-
-@deftypefun {void *} gal_data_calloc_array (uint8_t @code{type}, size_t 
@code{size}, const char @code{*funcname}, const char @code{*varname})
-Similar to @code{gal_data_malloc_array}, but the space is cleared (set to
-0) after allocation.
-@end deftypefun
+@code{gal_data_t}: those related to allocation and freeing. These functions
+are declared in @file{gnuastro/data.h} which is also visible from the
+function names (see @ref{Gnuastro library}).
 
 @deftypefun void gal_data_initialize (gal_data_t @code{*data}, void 
@code{*array}, uint8_t @code{type}, size_t @code{ndim}, size_t @code{*dsize}, 
struct wcsprm @code{*wcs}, int @code{clear}, size_t @code{minmapsize}, char 
@code{*name}, char @code{*unit}, char @code{*comment})
 
@@ -21587,7 +21665,7 @@ Free all the non-@code{NULL} pointers in 
@code{gal_data_t}, then free the
 actual data structure.
 @end deftypefun
 
-@node Arrays of datasets, Copying datasets, Dataset size and allocation, 
Library data container
+@node Arrays of datasets, Copying datasets, Dataset allocation, Library data 
container
 @subsubsection Arrays of datasets
 
 Gnuastro's generic data container (@code{gal_data_t}) is a very versatile
@@ -21723,6 +21801,12 @@ Return the total number of elements for a dataset with 
@code{ndim}
 dimensions that has @code{dsize} elements along each dimension.
 @end deftypefun
 
+@deftypefun int gal_dimension_is_different (gal_data_t @code{*first}, 
gal_data_t @code{*second})
+Return @code{1} (one) if the two datasets don't have the same size along
+all dimensions. This function will also return @code{1} when the number of
+dimensions of the two datasets are different.
+@end deftypefun
+
 @deftypefun {size_t *} gal_dimension_increment (size_t @code{ndim}, size_t 
@code{*dsize})
 Return an allocated array that has the number of elements necessary to
 increment an index along every dimension. For example along the fastest
@@ -22594,7 +22678,7 @@ gal_list_data_add( &list, tmp );
 @end deftypefun
 
 @deftypefun void gal_list_data_add_alloc (gal_data_t @code{**list}, void 
@code{*array}, uint8_t @code{type}, size_t @code{ndim}, size_t @code{*dsize}, 
struct wcsprm @code{*wcs}, int @code{clear}, size_t @code{minmapsize}, char 
@code{*name}, char @code{*unit}, char @code{*comment})
-Allocate a new dataset (with @code{gal_data_alloc} in @ref{Dataset size and
+Allocate a new dataset (with @code{gal_data_alloc} in @ref{Dataset
 allocation}) and put it as the first element of @code{list}. Note that if
 this is the first node to be added to the list, @code{list} must be
 @code{NULL}.
@@ -24571,8 +24655,8 @@ the main parameters of this two-layer tessellation and 
help in benefiting
 from it.
 
 @deftp {Type (C @code{struct})} gal_tile_two_layer_params
-This is the general structure to keep all the necessary parameters for a
-two-layer tessellation.
+The general structure to keep all the necessary parameters for a two-layer
+tessellation.
 
 @example
 struct gal_tile_two_layer_params
@@ -24994,22 +25078,65 @@ both polygons have to be sorted in an anti-clock-wise 
manner.
 @node Qsort functions, Permutations, Polygons, Gnuastro library
 @subsection Qsort functions (@file{qsort.h})
 
-The C programming language comes with the @code{qsort} (Quick sort)
-function. @code{qsort} is a generic function which allows you to sort any
-kind of data structure (not just a single number). To define greater and
-smaller (for sorting), @code{qsort} needs another function, even for simple
-numerical types. To facilitate numerical sorting for Gnuastro's
-programs/libraries, Gnuastro defines a function for each type's increasing
-and decreasing function. You can pass these functions to @code{qsort} when
-your array has the respective type (see @ref{Numeric data types}).
+When sorting a dataset is necessary, the C programming language provides
+the @code{qsort} (Quick sort) function. @code{qsort} is a generic function
+which allows you to sort any kind of data structure (not just a single
+array of numbers). To define ``greater'' and ``smaller'' (for sorting),
+@code{qsort} needs another function, even for simple numerical types. The
+functions introduced in this section are to passed onto @code{qsort}.
+
+The first class of functions below (with @code{TYPE} in their names) can be
+used for sorting a simple numeric array. Just replace @code{TYPE} with the
+dataset's numeric datatype. The second set of functions can be used to sort
+indexs (leave the actual numbers untouched). To use the second set of
+functions, a global variable or structure are also necessary as described
+below.
 
-@deffn {Global variable} {gal_qsort_index_arr}
+@deffn {Global variable} {gal_qsort_index_single}
+@cindex Thread-safety
+@cindex Multi-threaded operation
 Pointer to a floating point array (@code{float *}) to use as a reference in
-@code{gal_qsort_index_float_decreasing}, see the explanation there for
-more.
+@code{gal_qsort_index_single_d} or @code{gal_qsort_index_single_i}, see the
+explanation of these functions for more. Note that if more than one array
+is to be sorted in a multi-threaded operation, the two functions will not
+work as expected. Therefore, if the operation is multi-threaded, but all
+the threads just sort the indexs based on a single array, this global
+variable can safely be used.
 @end deffn
 
-@deftypefun int gal_qsort_index_float_decreasing (const void @code{*a}, const 
void @code{*b})
+@deftp {Type (C @code{struct})} gal_qsort_index_multi
+Structure to get the sorted indexs of multiple datasets on multiple threads
+with @code{gal_qsort_index_multi_d} or @code{gal_qsort_index_multi_i}. Note
+that the @code{values} array will not be changed by these functions, it is
+only read. Therefore all the @code{values} elements in the (to be sorted)
+array of @code{gal_qsort_index_multi} must point to the same place.
+
+@example
+struct gal_qsort_index_multi
+@{
+  float *values;         /* Array of values (same in all).      */
+  size_t  index;         /* Index of each element to be sorted. */
+@};
+@end example
+@end deftp
+
+@deftypefun int gal_qsort_TYPE_d (const void @code{*a}, const void @code{*b})
+When passed to @code{qsort}, this function will sort a @code{TYPE} array in
+decreasing order (first element will be the largest). Please replace
+@code{TYPE} (in the function name) with one of the @ref{Numeric data
+types}, for example @code{gal_qsort_int32_d}, or
+@code{gal_qsort_float64_d}.
+@end deftypefun
+
+@deftypefun int gal_qsort_TYPE_i (const void @code{*a}, const void @code{*b})
+When passed to @code{qsort}, this function will sort a @code{TYPE} array in
+increasing order (first element will be the smallest). Please replace
+@code{TYPE} (in the function name) with one of the @ref{Numeric data
+types}, for example @code{gal_qsort_int32_i}, or
+@code{gal_qsort_float64_i}.
+@end deftypefun
+
+@deftypefun int gal_qsort_index_single_d (const void @code{*a}, const void 
@code{*b})
 When passed to @code{qsort}, this function will sort a @code{size_t} array
 based on decreasing values in the @code{gal_qsort_index_arr} single
 precision floating point array. The floating point array will not be
@@ -25027,7 +25154,7 @@ main (void)
   size_t s[4]=@{0, 1, 2, 3@};
   float f[4]=@{1.3,0.2,1.8,0.1@};
   gal_qsort_index_arr=f;
-  qsort(s, 4, sizeof(size_t), gal_qsort_index_float_decreasing);
+  qsort(s, 4, sizeof(size_t), gal_qsort_index_single_d);
   printf("%zu, %zu, %zu, %zu\n", s[0], s[1], s[2], s[3]);
   return EXIT_SUCCESS;
 @}
@@ -25037,27 +25164,28 @@ main (void)
 The output will be: @code{2, 0, 1, 3}.
 @end deftypefun
 
-@deftypefun int gal_qsort_index_float_increasing (const void @code{*a}, const 
void @code{*b})
-Similar to @code{gal_qsort_index_float_decreasing}, but will sort the
-indexes such that the values of @code{gal_qsort_index_arr} are in
-increasing order.
+@deftypefun int gal_qsort_index_single_i (const void @code{*a}, const void 
@code{*b})
+Similar to @code{gal_qsort_index_single_d}, but will sort the indexes such
+that the values of @code{gal_qsort_index_arr} are in increasing order.
 @end deftypefun
 
+@deftypefun int gal_qsort_index_multi_d (const void @code{*a}, const void 
@code{*b})
+When passed to @code{qsort} with an array of @code{gal_qsort_index_multi},
+this function will sort the array based on the values of the given
+indexs. The sorting will be ordered according to the @code{values} pointer
+of @code{gal_qsort_index_multi}. Note that @code{values} must point to the
+same place in all the structures of the @code{gal_qsort_index_multi} array.
 
-@deftypefun int gal_qsort_TYPE_increasing (const void @code{*a}, const void 
@code{*b})
-When passed to @code{qsort}, this function will sort an @code{TYPE} array
-in increasing order (first element will be the smallest). Please replace
-@code{TYPE} (in the function name) with one of the @ref{Numeric data
-types}, for example @code{gal_qsort_int32_increasing}, or
-@code{gal_qsort_float64_increasing}.
+This function is only useful when the the indexs of multiple arrays on
+multiple threads are to be sorted. If your program is single threaded, or
+all the indexs belong to a single array (sorting different sub-sets of
+indexs in a single array on multiple threads), it is recommended to use
+@code{gal_qsort_index_single_d}.
 @end deftypefun
 
-@deftypefun int gal_qsort_TYPE_decreasing (const void @code{*a}, const void 
@code{*b})
-When passed to @code{qsort}, this function will sort an @code{TYPE} array
-in decreasing order (first element will be the largest). Please replace
-@code{TYPE} (in the function name) with one of the @ref{Numeric data
-types}, for example @code{gal_qsort_int32_decreasing}, or
-@code{gal_qsort_float64_decreasing}.
+@deftypefun int gal_qsort_index_multi_i (const void @code{*a}, const void 
*@code{b})
+Similar to @code{gal_qsort_index_multi_d}, but the result will be sorted in
+increasing order (first element will have the smallest value).
 @end deftypefun
 
 
@@ -25824,13 +25952,26 @@ The @code{values} dataset must have a 32-bit floating 
point type
 read by this function. @code{indexs} must contain the indexs of the
 elements/pixels that will be over-segmented by this function and have a
 @code{GAL_TYPE_SIZE_T} type, see the description of
-@code{gal_label_indexs}, above. After this function, @code{indexs} will be
-sorted by the values of the respective pixel in @code{values}. The final
-labels will be written in the respective positions of @code{labels}, which
-must have a @code{GAL_TYPE_INT32} type and be the same size as
-@code{values}. All local minima (maxima), when @code{min0_max1} is @code{1}
-(@code{0}), are considered rivers (watersheds) and given a label of
-@code{GAL_LABEL_RIVER} (see above).
+@code{gal_label_indexs}, above. The final labels will be written in the
+respective positions of @code{labels}, which must have a
+@code{GAL_TYPE_INT32} type and be the same size as @code{values}.
+
+When @code{indexs} is already sorted, this function will ignore
+@code{min0_max1}. To judge if the dataset is sorted or not (by the values
+the indexs correspond to in @code{values}, not the actual indexs), this
+function will look into the bits of @code{indexs->flag}, for the respective
+bit flags, see @ref{Generic data container}. If @code{indexs} is not
+already sorted, this function will sort it according to the values of the
+respective pixel in @code{values}. The increasing/decreasing order will be
+determined by @code{min0_max1}. Note that if this function is called on
+multiple threads @emph{and} @code{values} points to a different array on
+each thread, this function will not return a reasonable result. In this
+case, please sort @code{indexs} prior to calling this function (see
+@code{gal_qsort_index_multi_d} in @ref{Qsort functions}).
+
+When @code{indexs} is decreasing (increasing), or @code{min0_max1} is
+@code{1} (@code{0}), local minima (maxima), are considered rivers
+(watersheds) and given a label of @code{GAL_LABEL_RIVER} (see above).
 
 Note that rivers/watersheds will also be formed on the edges of the labeled
 regions or when the labeled pixels touch a blank pixel. Therefore this
@@ -25839,7 +25980,7 @@ efficient, it is thus recommended to use 
@code{gal_blank_present} (with
 @code{updateflag=1}) prior to calling this function (see @ref{Library blank
 values}. Once the flag has been set, no other function (including this one)
 that needs special behavior for blank pixels will have to parse the dataset
-any more.
+to see if it has blank values any more.
 
 If you are sure your dataset doesn't have blank values (by the design of
 your software), to avoid an extra parsing of the dataset and improve
@@ -26279,8 +26420,8 @@ system and you can safely use it as a guide.
 #include <stdio.h>
 #include <stdlib.h>
 
-#include "gnuastro/fits.h"
-#include "gnuastro/threads.h"
+#include <gnuastro/fits.h>
+#include <gnuastro/threads.h>
 
 
 /* This structure can keep all information you want to pass onto the
diff --git a/lib/Makefile.am b/lib/Makefile.am
index c978b4b..4445f97 100644
--- a/lib/Makefile.am
+++ b/lib/Makefile.am
@@ -59,8 +59,9 @@ libgnuastro_la_SOURCES = arithmetic.c arithmetic-and.c 
arithmetic-bitand.c \
   arithmetic-or.c arithmetic-plus.c array.c binary.c blank.c box.c         \
   checkset.c convolve.c cosmology.c data.c eps.c fits.c git.c              \
   interpolate.c jpeg.c label.c list.c match.c options.c pdf.c              \
-  permutation.c polygon.c qsort.c dimension.c statistics.c table.c         \
-  tableintern.c threads.c tiff.c tile.c timing.c txt.c type.c wcs.c
+  permutation.c pointer.c polygon.c qsort.c dimension.c statistics.c       \
+  table.c tableintern.c threads.c tiff.c tile.c timing.c txt.c type.c      \
+  wcs.c
 
 
 
@@ -77,10 +78,10 @@ pkginclude_HEADERS = gnuastro/config.h 
$(headersdir)/arithmetic.h          \
   $(headersdir)/fits.h $(headersdir)/git.h $(headersdir)/interpolate.h     \
   $(headersdir)/jpeg.h $(headersdir)/label.h $(headersdir)/list.h          \
   $(headersdir)/match.h $(headersdir)/pdf.h $(headersdir)/permutation.h    \
-  $(headersdir)/polygon.h $(headersdir)/qsort.h $(headersdir)/statistics.h \
-  $(headersdir)/table.h $(headersdir)/threads.h $(headersdir)/tiff.h       \
-  $(headersdir)/tile.h $(headersdir)/txt.h $(headersdir)/type.h            \
-  $(headersdir)/wcs.h
+  $(headersdir)/pointer.h $(headersdir)/polygon.h $(headersdir)/qsort.h    \
+  $(headersdir)/statistics.h $(headersdir)/table.h $(headersdir)/threads.h \
+  $(headersdir)/tiff.h $(headersdir)/tile.h $(headersdir)/txt.h            \
+  $(headersdir)/type.h $(headersdir)/wcs.h
 
 
 
diff --git a/lib/arithmetic.c b/lib/arithmetic.c
index 42a2300..6f31f51 100644
--- a/lib/arithmetic.c
+++ b/lib/arithmetic.c
@@ -30,6 +30,8 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 #include <gnuastro/blank.h>
 #include <gnuastro/qsort.h>
+#include <gnuastro/pointer.h>
+#include <gnuastro/dimension.h>
 #include <gnuastro/statistics.h>
 #include <gnuastro/arithmetic.h>
 
@@ -553,7 +555,7 @@ arithmetic_where(int flags, gal_data_t *out, gal_data_t 
*cond,
 
   /* The dimension and sizes of the out and condition data sets must be the
      same. */
-  if(gal_data_dsize_is_different(out, cond))
+  if( gal_dimension_is_different(out, cond) )
     error(EXIT_FAILURE, 0, "%s: the output and condition data sets of the "
           "must be the same size", __func__);
 
@@ -785,7 +787,8 @@ arithmetic_where(int flags, gal_data_t *out, gal_data_t 
*cond,
 #define MULTIOPERAND_MEDIAN(TYPE, QSORT_F) {                            \
     int use;                                                            \
     size_t n, j=0;                                                      \
-    TYPE *pixs=gal_data_malloc_array(list->type, dnum, __func__, "pixs"); \
+    TYPE *pixs=gal_pointer_allocate(list->type, dnum, 0, __func__,      \
+                                    "pixs");                            \
                                                                         \
     /* Loop over each pixel */                                          \
     do                                                                  \
@@ -919,7 +922,7 @@ arithmetic_multioperand(int operator, int flags, gal_data_t 
*list)
               gal_arithmetic_operator_string(operator));
 
       /* Check the sizes. */
-      if( gal_data_dsize_is_different(list, tmp) )
+      if( gal_dimension_is_different(list, tmp) )
         error(EXIT_FAILURE, 0, "%s: the sizes of all operands to the %s "
               "operator must be same", __func__,
               gal_arithmetic_operator_string(operator));
@@ -936,7 +939,8 @@ arithmetic_multioperand(int operator, int flags, gal_data_t 
*list)
 
   /* hasblank is used to see if a blank value should be checked for each
      list element or not. */
-  hasblank=gal_data_malloc_array(GAL_TYPE_UINT8, dnum, __func__, "hasblank");
+  hasblank=gal_pointer_allocate(GAL_TYPE_UINT8, dnum, 0, __func__,
+                                "hasblank");
   for(tmp=list;tmp!=NULL;tmp=tmp->next)
     hasblank[i++]=gal_blank_present(tmp, 0);
 
@@ -945,34 +949,34 @@ arithmetic_multioperand(int operator, int flags, 
gal_data_t *list)
   switch(list->type)
     {
     case GAL_TYPE_UINT8:
-      MULTIOPERAND_TYPE_SET(uint8_t,   gal_qsort_uint8_increasing);
+      MULTIOPERAND_TYPE_SET(uint8_t,   gal_qsort_uint8_i);
       break;
     case GAL_TYPE_INT8:
-      MULTIOPERAND_TYPE_SET(int8_t,    gal_qsort_int8_increasing);
+      MULTIOPERAND_TYPE_SET(int8_t,    gal_qsort_int8_i);
       break;
     case GAL_TYPE_UINT16:
-      MULTIOPERAND_TYPE_SET(uint16_t,  gal_qsort_uint16_increasing);
+      MULTIOPERAND_TYPE_SET(uint16_t,  gal_qsort_uint16_i);
       break;
     case GAL_TYPE_INT16:
-      MULTIOPERAND_TYPE_SET(int16_t,   gal_qsort_int16_increasing);
+      MULTIOPERAND_TYPE_SET(int16_t,   gal_qsort_int16_i);
       break;
     case GAL_TYPE_UINT32:
-      MULTIOPERAND_TYPE_SET(uint32_t,  gal_qsort_uint32_increasing);
+      MULTIOPERAND_TYPE_SET(uint32_t,  gal_qsort_uint32_i);
       break;
     case GAL_TYPE_INT32:
-      MULTIOPERAND_TYPE_SET(int32_t,   gal_qsort_int32_increasing);
+      MULTIOPERAND_TYPE_SET(int32_t,   gal_qsort_int32_i);
       break;
     case GAL_TYPE_UINT64:
-      MULTIOPERAND_TYPE_SET(uint64_t,  gal_qsort_uint64_increasing);
+      MULTIOPERAND_TYPE_SET(uint64_t,  gal_qsort_uint64_i);
       break;
     case GAL_TYPE_INT64:
-      MULTIOPERAND_TYPE_SET(int64_t,   gal_qsort_int64_increasing);
+      MULTIOPERAND_TYPE_SET(int64_t,   gal_qsort_int64_i);
       break;
     case GAL_TYPE_FLOAT32:
-      MULTIOPERAND_TYPE_SET(float,     gal_qsort_float32_increasing);
+      MULTIOPERAND_TYPE_SET(float,     gal_qsort_float32_i);
       break;
     case GAL_TYPE_FLOAT64:
-      MULTIOPERAND_TYPE_SET(double,    gal_qsort_float64_increasing);
+      MULTIOPERAND_TYPE_SET(double,    gal_qsort_float64_i);
       break;
     default:
       error(EXIT_FAILURE, 0, "%s: type code %d not recognized",
@@ -1071,7 +1075,7 @@ arithmetic_binary(int operator, int flags, gal_data_t *l, 
gal_data_t *r)
 
   /* Simple sanity check on the input sizes */
   if( !( (flags & GAL_ARITHMETIC_NUMOK) && (l->size==1 || r->size==1))
-      && gal_data_dsize_is_different(l, r) )
+      && gal_dimension_is_different(l, r) )
     error(EXIT_FAILURE, 0, "%s: the non-number inputs to %s don't have the "
           "same dimension/size", __func__,
           gal_arithmetic_operator_string(operator));
@@ -1230,7 +1234,7 @@ arithmetic_binary_function_flt(int operator, int flags, 
gal_data_t *l,
 
   /* Simple sanity check on the input sizes */
   if( !( (flags & GAL_ARITHMETIC_NUMOK) && (l->size==1 || r->size==1))
-      && gal_data_dsize_is_different(l, r) )
+      && gal_dimension_is_different(l, r) )
     error(EXIT_FAILURE, 0, "%s: the input datasets don't have the same "
           "dimension/size", __func__);
 
diff --git a/lib/binary.c b/lib/binary.c
index 39369c9..3eed8c1 100644
--- a/lib/binary.c
+++ b/lib/binary.c
@@ -32,6 +32,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/tile.h>
 #include <gnuastro/blank.h>
 #include <gnuastro/binary.h>
+#include <gnuastro/pointer.h>
 #include <gnuastro/dimension.h>
 
 
@@ -410,7 +411,7 @@ gal_binary_connected_components(gal_data_t *binary, 
gal_data_t **out,
       lab=*out;
 
       /* Make sure the given dataset has the same size as the input. */
-      if( gal_data_dsize_is_different(binary, lab) )
+      if( gal_dimension_is_different(binary, lab) )
         error(EXIT_FAILURE, 0, "%s: the `binary' and `out' datasets must "
               "have the same size", __func__);
 
@@ -604,10 +605,10 @@ binary_make_padded_inverse(gal_data_t *input, gal_data_t 
**outtile)
   uint8_t *in;
   size_t i, startind;
   gal_data_t *inv, *tile;
-  size_t *startcoord=gal_data_malloc_array(GAL_TYPE_SIZE_T, input->ndim,
-                                           __func__, "startcoord");
-  size_t *dsize=gal_data_malloc_array(GAL_TYPE_SIZE_T, input->ndim, __func__,
-                                      "dsize");
+  size_t *startcoord=gal_pointer_allocate(GAL_TYPE_SIZE_T, input->ndim, 0,
+                                          __func__, "startcoord");
+  size_t *dsize=gal_pointer_allocate(GAL_TYPE_SIZE_T, input->ndim, 0,
+                                     __func__, "dsize");
 
 
   /* Set the size of the padded inverse image and the coordinates of the
@@ -634,7 +635,7 @@ binary_make_padded_inverse(gal_data_t *input, gal_data_t 
**outtile)
 
   /* Define a tile to fill the central regions of the inverse. */
   startind=gal_dimension_coord_to_index(input->ndim, inv->dsize, startcoord);
-  tile=gal_data_alloc(gal_data_ptr_increment(inv->array, startind, inv->type),
+  tile=gal_data_alloc(gal_pointer_increment(inv->array, startind, inv->type),
                       inv->type, input->ndim, input->dsize, NULL, 0, 0, NULL,
                       NULL, NULL);
   *outtile=tile;
@@ -794,8 +795,8 @@ gal_binary_holes_fill(gal_data_t *input, int connectivity, 
size_t maxsize)
   if(maxsize<-1)
     {
       /* Allocate space to keep the size of each hole: */
-      sizes=gal_data_calloc_array(GAL_TYPE_SIZE_T, numholes+1, __func__,
-                                  "sizes");
+      sizes=gal_pointer_allocate(GAL_TYPE_SIZE_T, numholes+1, 1, __func__,
+                                 "sizes");
       fi=(i=holelabs->array)+holelabs->size; do ++sizes[*i]; while(++i<fi);
 
       /* Set those labels with a larger size to 1 (treat it as
diff --git a/lib/blank.c b/lib/blank.c
index a3ebc19..b84a5ce 100644
--- a/lib/blank.c
+++ b/lib/blank.c
@@ -32,6 +32,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/data.h>
 #include <gnuastro/tile.h>
 #include <gnuastro/blank.h>
+#include <gnuastro/pointer.h>
 
 #include <gnuastro-internal/checkset.h>
 
@@ -87,7 +88,7 @@ gal_blank_alloc_write(uint8_t type)
   void *out;
 
   /* Allocate the space to keep the blank value. */
-  out=gal_data_malloc_array(type, 1, __func__, "out");
+  out=gal_pointer_allocate(type, 1, 0, __func__, "out");
 
   /* Put the blank value in the allcated space. */
   gal_blank_write(out, type);
diff --git a/lib/convolve.c b/lib/convolve.c
index 0df670c..1f558b4 100644
--- a/lib/convolve.c
+++ b/lib/convolve.c
@@ -31,6 +31,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/list.h>
 #include <gnuastro/tile.h>
 #include <gnuastro/threads.h>
+#include <gnuastro/pointer.h>
 #include <gnuastro/convolve.h>
 #include <gnuastro/dimension.h>
 
@@ -271,8 +272,8 @@ convolve_spatial_overlap(struct per_thread_spatial_prm 
*pprm, int tocorrect)
   /* Set the starting point of the dataset overlap tile. */
   increment=gal_dimension_coord_to_index(ndim, block->dsize,
                                          pprm->overlap_start);
-  pprm->i_overlap->array=gal_data_ptr_increment(block->array, increment,
-                                                block->type);
+  pprm->i_overlap->array=gal_pointer_increment(block->array, increment,
+                                               block->type);
 
 
   /* Set the starting point of the kernel overlap tile. */
@@ -281,8 +282,8 @@ convolve_spatial_overlap(struct per_thread_spatial_prm 
*pprm, int tocorrect)
                 : gal_dimension_coord_to_index(ndim,
                                                kernel->dsize,
                                                pprm->kernel_start) );
-  pprm->k_overlap->array=gal_data_ptr_increment(kernel->array, increment,
-                                                kernel->type);
+  pprm->k_overlap->array=gal_pointer_increment(kernel->array, increment,
+                                               kernel->type);
   return full_overlap;
 }
 
@@ -423,8 +424,8 @@ convolve_spatial_on_thread(void *inparam)
   size_t i;
   size_t ndim=block->ndim;
   struct per_thread_spatial_prm *pprm=&cprm->pprm[tprm->id];
-  size_t *dsize=gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim, __func__,
-                                      "dsize");
+  size_t *dsize=gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__,
+                                     "dsize");
 
 
   /* Set all dsize values to 1 (the values within `overlap->dsize' will be
@@ -434,15 +435,14 @@ convolve_spatial_on_thread(void *inparam)
 
   /* Initialize/Allocate necessary items for this thread. */
   pprm->cprm          = cprm;
-  pprm->pix           = gal_data_malloc_array(GAL_TYPE_SIZE_T, 2*ndim,
-                                              __func__, "pprm->pix");
-  pprm->host_start    = gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim,
-                                              __func__, "pprm->host_start");
-  pprm->kernel_start  = gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim,
-                                              __func__, "pprm->kernel_start");
-  pprm->overlap_start = gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim,
-                                               __func__,
-                                              "pprm->overlap_start");
+  pprm->pix           = gal_pointer_allocate(GAL_TYPE_SIZE_T, 2*ndim, 0,
+                                             __func__, "pprm->pix");
+  pprm->host_start    = gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0,
+                                             __func__, "pprm->host_start");
+  pprm->kernel_start  = gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0,
+                                             __func__, "pprm->kernel_start");
+  pprm->overlap_start = gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0,
+                                             __func__, "pprm->overlap_start");
   pprm->i_overlap     = gal_data_alloc(NULL, block->type, ndim, dsize,
                                        NULL, 0, -1, NULL, NULL, NULL);
   pprm->k_overlap     = gal_data_alloc(NULL, cprm->kernel->type, ndim, dsize,
@@ -582,7 +582,7 @@ gal_convolve_spatial_correct_ch_edge(gal_data_t *tiles, 
gal_data_t *kernel,
   gal_data_t *block=gal_tile_block(tiles);
 
   /* Some small sanity checks. */
-  if( gal_data_dsize_is_different(block, tocorrect) )
+  if( gal_dimension_is_different(block, tocorrect) )
     error(EXIT_FAILURE, 0, "%s: the `tocorrect' dataset has to have the "
           "same dimensions/size as the block of the `tiles' input", __func__);
   if( block->type != tocorrect->type )
diff --git a/lib/data.c b/lib/data.c
index 58f03d8..b1e40a5 100644
--- a/lib/data.c
+++ b/lib/data.c
@@ -27,18 +27,17 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <fcntl.h>
 #include <float.h>
 #include <ctype.h>
-#include <unistd.h>
 #include <stdlib.h>
 #include <string.h>
 #include <dirent.h>
 #include <inttypes.h>
-#include <sys/mman.h>
 
 #include <gnuastro/wcs.h>
 #include <gnuastro/data.h>
 #include <gnuastro/tile.h>
 #include <gnuastro/blank.h>
 #include <gnuastro/table.h>
+#include <gnuastro/pointer.h>
 
 #include <gnuastro-internal/checkset.h>
 
@@ -62,202 +61,8 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 
 /*********************************************************************/
-/*************          Size and allocation        *******************/
+/*************              Allocation             *******************/
 /*********************************************************************/
-int
-gal_data_dsize_is_different(gal_data_t *first, gal_data_t *second)
-{
-  size_t i;
-
-  /* First make sure that the dimensionality is the same. */
-  if(first->ndim!=second->ndim)
-    return 1;
-
-  /* Check if the sizes along all dimensions are the same: */
-  for(i=0;i<first->ndim;++i)
-    if( first->dsize[i] != second->dsize[i] )
-      return 1;
-
-  /* If it got to here, we know the dimensions have the same length. */
-  return 0;
-}
-
-
-
-
-
-/* Increment a give pointer depending on the given type.
-
-   When working with the `array' elements of `gal_data_t', we are actually
-   dealing with `void *' pointers. Pointer arithmetic doesn't apply to
-   `void *', because the system doesn't know how much space each element
-   has to increment the pointer respectively.
-
-   So, here, we will use the type information to find the increment. This
-   is mainly useful when dealing with the `block' pointer of a tile over a
-   larger image. This function reads the address as a `char *' type (note
-   that `char' is guaranteed to have a size of 1 (byte)). It then
-   increments the `char *' by `increment*sizeof(type)' */
-void *
-gal_data_ptr_increment(void *pointer, size_t increment, uint8_t type)
-{
-  char *p=(char *)pointer;
-  return p + increment * gal_type_sizeof(type);
-}
-
-
-
-
-
-/* Find the number of values between two void pointers with a given
-   type. See the explanations before `gal_data_ptr_increment'. */
-size_t
-gal_data_num_between(void *earlier, void *later, uint8_t type)
-{
-  char *e=(char *)earlier, *l=(char *)later;
-  return (l-e)/gal_type_sizeof(type);
-}
-
-
-
-
-
-/* Allocate an array based on the value of type. Note that the argument
-   `size' is the number of elements, necessary in the array, the number of
-   bytes each element needs will be determined internaly by this function
-   using the datatype argument, so you don't have to worry about it. */
-void *
-gal_data_malloc_array(uint8_t type, size_t size, const char *funcname,
-                      const char *varname)
-{
-  void *array;
-
-  errno=0;
-  array=malloc( size * gal_type_sizeof(type) );
-  if(array==NULL)
-    {
-      if(varname)
-        error(EXIT_FAILURE, errno, "%s: %zu bytes couldn't be allocated "
-              "for variable `%s'", funcname ? funcname : __func__,
-              size * gal_type_sizeof(type), varname);
-      else
-        error(EXIT_FAILURE, errno, "%s: %zu bytes couldn't be allocated",
-              funcname ? funcname : __func__, size * gal_type_sizeof(type));
-    }
-
-  return array;
-}
-
-
-
-
-
-void *
-gal_data_calloc_array(uint8_t type, size_t size, const char *funcname,
-                      const char *varname)
-{
-  void *array;
-
-  errno=0;
-  array=calloc( size, gal_type_sizeof(type) );
-  if(array==NULL)
-    {
-      if(varname)
-        error(EXIT_FAILURE, errno, "%s: %zu bytes couldn't be allocated "
-              "for variable `%s'", funcname ? funcname : __func__,
-              size * gal_type_sizeof(type), varname);
-      else
-        error(EXIT_FAILURE, errno, "%s: %zu bytes couldn't be allocated",
-              funcname ? funcname : __func__, size * gal_type_sizeof(type));
-    }
-
-  return array;
-}
-
-
-
-
-
-static void
-gal_data_mmap(gal_data_t *data, int clear, size_t minmapsize)
-{
-  int filedes;
-  uint8_t uc=0;
-  char *filename;
-  size_t bsize=data->size*gal_type_sizeof(data->type);
-
-
-  /* Check if the .gnuastro folder exists, write the file there. If it
-     doesn't exist, then make the .gnuastro directory.*/
-  gal_checkset_mkdir(".gnuastro");
-
-
-  /* Set the filename */
-  gal_checkset_allocate_copy("./.gnuastro/mmap_XXXXXX", &filename);
-
-
-  /* Create a zero-sized file and keep its descriptor.  */
-  errno=0;
-  /*filedes=open(filename, O_RDWR | O_CREAT | O_EXCL | O_TRUNC );*/
-  filedes=mkstemp(filename);
-  if(filedes==-1)
-    error(EXIT_FAILURE, errno, "%s: %s couldn't be created",
-          __func__, filename);
-
-
-  /* Make enough space to keep the array data. */
-  errno=0;
-  if( lseek(filedes, bsize, SEEK_SET) == -1 )
-    error(EXIT_FAILURE, errno, "%s: %s: unable to change file position by "
-          "%zu bytes", __func__, filename, bsize);
-
-
-  /* Write to the newly set file position so the space is allocated. To do
-     this, we are simply writing `uc' (a byte with value 0) into the space
-     we identified by `lseek' (above). This will ensure that this space is
-     set a side for this array and prepare us to use `mmap'. */
-  if( write(filedes, &uc, 1) == -1)
-    error(EXIT_FAILURE, errno, "%s: %s: unable to write one byte at the "
-          "%zu-th position", __func__, filename, bsize);
-
-
-  /* Map the memory. */
-  errno=0;
-  data->array=mmap(NULL, bsize, PROT_READ | PROT_WRITE, MAP_SHARED,
-                   filedes, 0);
-  if(data->array==MAP_FAILED)
-    {
-      if(minmapsize<10000u)
-        fprintf(stderr, "\nIf the processing involves many small mappings "
-                "(along with larger ones), the following error may be "
-                "corrected with a larger value to `minmapsize' (minimum "
-                "number of bytes to use mapping instead of RAM for each "
-                "patch of memory), for example 10000. In this way, mapping "
-                "will only be reserved for larger sizes. The current value is "
-                "%zu.\n\n", minmapsize);
-      error(EXIT_FAILURE, errno, "couldn't map %zu bytes into the file `%s'",
-            bsize, filename);
-    }
-
-
-  /* Close the file. */
-  if( close(filedes) == -1 )
-    error(EXIT_FAILURE, errno, "%s: %s couldn't be closed",
-          __func__, filename);
-
-
-  /* Keep the filename. */
-  data->mmapname=filename;
-
-
-  /* If it was supposed to be cleared, then clear the memory. */
-  if(clear) memset(data->array, 0, bsize);
-}
-
-
-
-
-
 /* Initialize the data structure.
 
    Some notes:
@@ -357,16 +162,13 @@ gal_data_initialize(gal_data_t *data, void *array, 
uint8_t type,
             {
               if( gal_type_sizeof(type)*data->size > minmapsize )
                 /* Allocate the space into disk (HDD/SSD). */
-                gal_data_mmap(data, clear, minmapsize);
+                data->array=gal_pointer_allocate_mmap(data->type, data->size,
+                                                      clear, &data->mmapname);
               else
                 /* Allocate the space in RAM. */
-                data->array = ( clear
-                                ? gal_data_calloc_array(data->type,
-                                                        data->size, __func__,
-                                                        "data->array")
-                                : gal_data_malloc_array(data->type,
-                                                        data->size, __func__,
-                                                        "data->array") );
+                data->array = gal_pointer_allocate(data->type, data->size,
+                                                   clear, __func__,
+                                                   "data->array");
             }
           else data->array=NULL; /* The given size was zero! */
         }
diff --git a/lib/dimension.c b/lib/dimension.c
index 00f568e..3d8d8e2 100644
--- a/lib/dimension.c
+++ b/lib/dimension.c
@@ -28,6 +28,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <error.h>
 #include <stdlib.h>
 
+#include <gnuastro/pointer.h>
 #include <gnuastro/dimension.h>
 
 
@@ -49,13 +50,35 @@ gal_dimension_total_size(size_t ndim, size_t *dsize)
 
 
 
+int
+gal_dimension_is_different(gal_data_t *first, gal_data_t *second)
+{
+  size_t i;
+
+  /* First make sure that the dimensionality is the same. */
+  if(first->ndim!=second->ndim)
+    return 1;
+
+  /* Check if the sizes along all dimensions are the same: */
+  for(i=0;i<first->ndim;++i)
+    if( first->dsize[i] != second->dsize[i] )
+      return 1;
+
+  /* If it got to here, we know the dimensions have the same length. */
+  return 0;
+}
+
+
+
+
+
 /* Calculate the values necessary to increment/decrement along each
    dimension of a dataset with size `dsize'. */
 size_t *
 gal_dimension_increment(size_t ndim, size_t *dsize)
 {
   int i;
-  size_t *out=gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim, __func__, "out");
+  size_t *out=gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__, "out");
 
   /* Along the fastest dimension, it is 1. */
   out[ndim-1]=1;
diff --git a/lib/fits.c b/lib/fits.c
index e584704..18ee668 100644
--- a/lib/fits.c
+++ b/lib/fits.c
@@ -39,6 +39,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/fits.h>
 #include <gnuastro/tile.h>
 #include <gnuastro/blank.h>
+#include <gnuastro/pointer.h>
 
 #include <gnuastro-internal/checkset.h>
 #include <gnuastro-internal/tableintern.h>
@@ -843,7 +844,7 @@ gal_fits_key_img_blank(uint8_t type)
      type. */
   if(tocopy)
     {
-      out = gal_data_malloc_array(type, 1, __func__, "out");
+      out = gal_pointer_allocate(type, 1, 0, __func__, "out");
       memcpy(out, tocopy, gal_type_sizeof(type));
     }
 
@@ -955,8 +956,8 @@ gal_fits_key_read_from_ptr(fitsfile *fptr, gal_data_t 
*keysll,
            set the size and ndim to 1. But first allocate dsize if it
            wasn't already allocated. */
         if(tmp->dsize==NULL)
-          tmp->dsize=gal_data_malloc_array(GAL_TYPE_SIZE_T, 1, __func__,
-                                           "tmp->dsize");
+          tmp->dsize=gal_pointer_allocate(GAL_TYPE_SIZE_T, 1, 0, __func__,
+                                          "tmp->dsize");
         tmp->ndim=tmp->size=tmp->dsize[0]=1;
 
         /* When the type is a string, `tmp->array' is an array of pointers
@@ -968,9 +969,9 @@ gal_fits_key_read_from_ptr(fitsfile *fptr, gal_data_t 
*keysll,
           case GAL_TYPE_STRING:
             tmp->array=strarray=( tmp->array
                                   ? tmp->array
-                                  : gal_data_malloc_array(tmp->type, 1,
-                                                          __func__,
-                                                          "tmp->array") );
+                                  : gal_pointer_allocate(tmp->type, 1, 0,
+                                                         __func__,
+                                                         "tmp->array") );
             errno=0;
             valueptr=strarray[0]=malloc(FLEN_VALUE * sizeof *strarray[0]);
             if(strarray[0]==NULL)
@@ -981,9 +982,9 @@ gal_fits_key_read_from_ptr(fitsfile *fptr, gal_data_t 
*keysll,
           default:
             tmp->array=valueptr=( tmp->array
                                   ? tmp->array
-                                  : gal_data_malloc_array(tmp->type, 1,
-                                                          __func__,
-                                                          "tmp->array") );
+                                  : gal_pointer_allocate(tmp->type, 1, 0,
+                                                         __func__,
+                                                         "tmp->array") );
           }
 
         /* Allocate space for the keyword comment if necessary.*/
@@ -1489,7 +1490,7 @@ gal_fits_img_info(fitsfile *fptr, int *type, size_t 
*ndim, size_t **dsize,
 
   /* Allocate the array to keep the dimension size and fill it in, note
      that its order is the opposite of naxes. */
-  *dsize=gal_data_malloc_array(GAL_TYPE_INT64, *ndim, __func__, "dsize");
+  *dsize=gal_pointer_allocate(GAL_TYPE_INT64, *ndim, 0, __func__, "dsize");
   for(i=0; i<*ndim; ++i)
     (*dsize)[i]=naxes[*ndim-1-i];
 
@@ -1542,7 +1543,7 @@ gal_fits_img_read(char *filename, char *hdu, size_t 
minmapsize)
      problem, the space will be written/read as 32-bit `long' any way,
      we'll just have a few empty bytes that will be freed anyway at the end
      of this function. */
-  fpixel=gal_data_malloc_array(GAL_TYPE_INT64, ndim, __func__, "fpixel");
+  fpixel=gal_pointer_allocate(GAL_TYPE_INT64, ndim, 0, __func__, "fpixel");
   for(i=0;i<ndim;++i) fpixel[i]=1;
 
 
@@ -1680,9 +1681,10 @@ gal_fits_img_write_to_ptr(gal_data_t *input, char 
*filename)
   hasblank=gal_blank_present(towrite, 0);
 
   /* Allocate the naxis area. */
-  naxes=gal_data_malloc_array( ( sizeof(long)==8
-                                 ? GAL_TYPE_INT64
-                                 : GAL_TYPE_INT32 ), ndim, __func__, "naxes");
+  naxes=gal_pointer_allocate( ( sizeof(long)==8
+                                ? GAL_TYPE_INT64
+                                : GAL_TYPE_INT32 ), ndim, 0, __func__,
+                              "naxes");
 
 
   /* Open the file for writing */
diff --git a/lib/gnuastro/data.h b/lib/gnuastro/data.h
index 2e0aeee..590917f 100644
--- a/lib/gnuastro/data.h
+++ b/lib/gnuastro/data.h
@@ -230,25 +230,8 @@ typedef struct gal_data_t
 
 
 /*********************************************************************/
-/*************         Size and allocation         *******************/
+/*************              allocation             *******************/
 /*********************************************************************/
-int
-gal_data_dsize_is_different(gal_data_t *first, gal_data_t *second);
-
-void *
-gal_data_ptr_increment(void *pointer, size_t increment, uint8_t type);
-
-size_t
-gal_data_num_between(void *earlier, void *later, uint8_t type);
-
-void *
-gal_data_malloc_array(uint8_t type, size_t size, const char *funcname,
-                      const char *varname);
-
-void *
-gal_data_calloc_array(uint8_t type, size_t size, const char *funcname,
-                      const char *varname);
-
 void
 gal_data_initialize(gal_data_t *data, void *array, uint8_t type, size_t ndim,
                     size_t *dsize, struct wcsprm *wcs, int clear,
diff --git a/lib/gnuastro/dimension.h b/lib/gnuastro/dimension.h
index 90c2c41..5c83a0d 100644
--- a/lib/gnuastro/dimension.h
+++ b/lib/gnuastro/dimension.h
@@ -53,6 +53,9 @@ __BEGIN_C_DECLS  /* From C++ preparations */
 size_t
 gal_dimension_total_size(size_t ndim, size_t *dsize);
 
+int
+gal_dimension_is_different(gal_data_t *first, gal_data_t *second);
+
 size_t *
 gal_dimension_increment(size_t ndim, size_t *dsize);
 
diff --git a/lib/gnuastro/pointer.h b/lib/gnuastro/pointer.h
new file mode 100644
index 0000000..e4292e6
--- /dev/null
+++ b/lib/gnuastro/pointer.h
@@ -0,0 +1,70 @@
+/*********************************************************************
+pointer -- facilitate working with pointers and allocation.
+This is part of GNU Astronomy Utilities (Gnuastro) package.
+
+Original author:
+     Mohammad Akhlaghi <mohammad@akhlaghi.org>
+Contributing author(s):
+Copyright (C) 2017-2018, Free Software Foundation, Inc.
+
+Gnuastro is free software: you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation, either version 3 of the License, or (at your
+option) any later version.
+
+Gnuastro is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
+**********************************************************************/
+#ifndef __GAL_POINTER_H__
+#define __GAL_POINTER_H__
+
+/* Include other headers if necessary here. Note that other header files
+   must be included before the C++ preparations below */
+#include <stdint.h>
+
+
+/* C++ Preparations */
+#undef __BEGIN_C_DECLS
+#undef __END_C_DECLS
+#ifdef __cplusplus
+# define __BEGIN_C_DECLS extern "C" {
+# define __END_C_DECLS }
+#else
+# define __BEGIN_C_DECLS                /* empty */
+# define __END_C_DECLS                  /* empty */
+#endif
+/* End of C++ preparations */
+
+
+/* Actual header contants (the above were for the Pre-processor). */
+__BEGIN_C_DECLS  /* From C++ preparations */
+
+
+
+
+
+void *
+gal_pointer_increment(void *pointer, size_t increment, uint8_t type);
+
+size_t
+gal_pointer_num_between(void *earlier, void *later, uint8_t type);
+
+void *
+gal_pointer_allocate(uint8_t type, size_t size, int clear,
+                     const char *funcname, const char *varname);
+
+void *
+gal_pointer_allocate_mmap(uint8_t type, size_t size, int clear,
+                          char **filename);
+
+
+
+
+
+__END_C_DECLS    /* From C++ preparations */
+#endif
diff --git a/lib/gnuastro/qsort.h b/lib/gnuastro/qsort.h
index 2b20876..c6e8b76 100644
--- a/lib/gnuastro/qsort.h
+++ b/lib/gnuastro/qsort.h
@@ -1,5 +1,5 @@
 /*********************************************************************
-forqsort -- Functions used by qsort to sort an array.
+qsort -- Functions used by qsort to sort an array.
 This is part of GNU Astronomy Utilities (Gnuastro) package.
 
 Original author:
@@ -48,77 +48,104 @@ __BEGIN_C_DECLS  /* From C++ preparations */
 
 
 
-/* Pointer used to sort the indexs of an array based on their flux
-   (value in this array). */
-extern float *gal_qsort_index_arr;
 
+
+/*****************************************************************/
+/**********        Sorting of actual array        ****************/
+/*****************************************************************/
 int
-gal_qsort_index_float_decreasing(const void * a, const void * b);
+gal_qsort_uint8_d(const void *a, const void *b);
 
 int
-gal_qsort_index_float_increasing(const void * a, const void * b);
+gal_qsort_uint8_i(const void *a, const void *b);
 
+int
+gal_qsort_int8_d(const void *a, const void *b);
 
+int
+gal_qsort_int8_i(const void *a, const void *b);
 
 int
-gal_qsort_uint8_decreasing(const void *a, const void *b);
+gal_qsort_uint16_d(const void *a, const void *b);
 
 int
-gal_qsort_uint8_increasing(const void *a, const void *b);
+gal_qsort_uint16_i(const void *a, const void *b);
 
 int
-gal_qsort_int8_decreasing(const void *a, const void *b);
+gal_qsort_int16_d(const void *a, const void *b);
 
 int
-gal_qsort_int8_increasing(const void *a, const void *b);
+gal_qsort_int16_i(const void *a, const void *b);
 
 int
-gal_qsort_uint16_decreasing(const void *a, const void *b);
+gal_qsort_uint32_d(const void *a, const void *b);
 
 int
-gal_qsort_uint16_increasing(const void *a, const void *b);
+gal_qsort_uint32_i(const void *a, const void *b);
 
 int
-gal_qsort_int16_decreasing(const void *a, const void *b);
+gal_qsort_int32_d(const void *a, const void *b);
 
 int
-gal_qsort_int16_increasing(const void *a, const void *b);
+gal_qsort_int32_i(const void *a, const void *b);
 
 int
-gal_qsort_uint32_decreasing(const void *a, const void *b);
+gal_qsort_uint64_d(const void *a, const void *b);
 
 int
-gal_qsort_uint32_increasing(const void *a, const void *b);
+gal_qsort_uint64_i(const void *a, const void *b);
 
 int
-gal_qsort_int32_decreasing(const void *a, const void *b);
+gal_qsort_int64_d(const void *a, const void *b);
 
 int
-gal_qsort_int32_increasing(const void *a, const void *b);
+gal_qsort_int64_i(const void *a, const void *b);
 
 int
-gal_qsort_uint64_decreasing(const void *a, const void *b);
+gal_qsort_float32_d(const void *a, const void *b);
 
 int
-gal_qsort_uint64_increasing(const void *a, const void *b);
+gal_qsort_float32_i(const void *a, const void *b);
 
 int
-gal_qsort_int64_decreasing(const void *a, const void *b);
+gal_qsort_float64_d(const void *a, const void *b);
 
 int
-gal_qsort_int64_increasing(const void *a, const void *b);
+gal_qsort_float64_i(const void *a, const void *b);
+
+
+
+
+
+/*****************************************************************/
+/***************          Sorting indexs        ******************/
+/*****************************************************************/
+/* Pointer used to sort the indexs of an array based on their flux (value
+   in this array). Note: when EACH THREAD USES A DIFFERENT ARRAY, this is
+   not thread-safe . */
+extern float *gal_qsort_index_single;
+
+
+/* When each thread is working on a different array, we'll need to keep the
+   pointer to the array in question for every index. */
+struct gal_qsort_index_multi
+{
+  float *values; /* Array of values (pointer, so original is not touched). */
+                 /* This should be identical in all elements.              */
+  size_t  index; /* Index of each element to be sorted.                    */
+};
 
 int
-gal_qsort_float32_decreasing(const void *a, const void *b);
+gal_qsort_index_single_d(const void *a, const void *b);
 
 int
-gal_qsort_float32_increasing(const void *a, const void *b);
+gal_qsort_index_single_i(const void *a, const void *b);
 
 int
-gal_qsort_float64_decreasing(const void *a, const void *b);
+gal_qsort_index_multi_d(const void *a, const void *b);
 
 int
-gal_qsort_float64_increasing(const void *a, const void *b);
+gal_qsort_index_multi_i(const void *a, const void *b);
 
 
 
diff --git a/lib/gnuastro/tile.h b/lib/gnuastro/tile.h
index 6397268..2699142 100644
--- a/lib/gnuastro/tile.h
+++ b/lib/gnuastro/tile.h
@@ -27,6 +27,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
    must be included before the C++ preparations below */
 #include <gnuastro/data.h>
 #include <gnuastro/fits.h>
+#include <gnuastro/dimension.h>
 
 /* C++ Preparations */
 #undef __BEGIN_C_DECLS
@@ -203,7 +204,7 @@ gal_tile_full_free_contents(struct 
gal_tile_two_layer_params *tl);
       {                                                                 \
         if( OTHER==tpo_oblock )    /* `OTHER' is a block. */            \
           {                                                             \
-            if(gal_data_dsize_is_different(tpo_iblock, tpo_oblock) )    \
+            if( gal_dimension_is_different(tpo_iblock, tpo_oblock) )    \
               {                                                         \
                 /* `error' function, is a GNU extension, see above. */  \
                 fprintf(stderr, "GAL_TILE_PO_OISET: when "              \
@@ -216,7 +217,7 @@ gal_tile_full_free_contents(struct 
gal_tile_two_layer_params *tl);
               }                                                         \
           }                                                             \
         else                                                            \
-          if(gal_data_dsize_is_different(IN, OTHER) )                   \
+          if( gal_dimension_is_different(IN, OTHER) )                   \
             {                                                           \
               /* The `error' function, is a GNU extension and this */   \
               /* is a header, not a library which the user has to  */   \
diff --git a/lib/interpolate.c b/lib/interpolate.c
index 8d66deb..4d9dceb 100644
--- a/lib/interpolate.c
+++ b/lib/interpolate.c
@@ -31,6 +31,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/list.h>
 #include <gnuastro/fits.h>
 #include <gnuastro/blank.h>
+#include <gnuastro/pointer.h>
 #include <gnuastro/threads.h>
 #include <gnuastro/dimension.h>
 #include <gnuastro/statistics.h>
@@ -92,12 +93,12 @@ interpolate_close_neighbors_on_thread(void *in_prm)
   size_t ngb_counter, dist, pind, *dinc;
   size_t i, index, fullind, chstart=0, ndim=input->ndim;
   gal_data_t *median, *tin, *tout, *tnear, *nearest=NULL;
-  size_t *icoord=gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim, __func__,
-                                       "icoord");
-  size_t *ncoord=gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim, __func__,
-                                       "ncoord");
   size_t size = (correct_index ? tl->tottilesinch : input->size);
   size_t *dsize = (correct_index ? tl->numtilesinch : input->dsize);
+  size_t *icoord=gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__,
+                                      "icoord");
+  size_t *ncoord=gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__,
+                                      "ncoord");
   uint8_t *fullflag=&prm->thread_flags[tprm->id*input->size], *flag=fullflag;
 
 
@@ -113,8 +114,8 @@ interpolate_close_neighbors_on_thread(void *in_prm)
   tin=input;
   for(tvll=prm->ngb_vals; tvll!=NULL; tvll=tvll->next)
     {
-      nv=gal_data_ptr_increment(tvll->v, tprm->id*prm->numneighbors,
-                                input->type);
+      nv=gal_pointer_increment(tvll->v, tprm->id*prm->numneighbors,
+                               input->type);
       gal_list_data_add_alloc(&nearest, nv, tin->type, 1, &prm->numneighbors,
                               NULL, 0, -1, NULL, NULL, NULL);
       tin=tin->next;
@@ -138,8 +139,8 @@ interpolate_close_neighbors_on_thread(void *in_prm)
           tin=input;
           for(tout=prm->out; tout!=NULL; tout=tout->next)
             {
-              memcpy(gal_data_ptr_increment(tout->array, fullind, tin->type),
-                     gal_data_ptr_increment(tin->array,  fullind, tin->type),
+              memcpy(gal_pointer_increment(tout->array, fullind, tin->type),
+                     gal_pointer_increment(tin->array,  fullind, tin->type),
                      gal_type_sizeof(tin->type));
               tin=tin->next;
             }
@@ -163,7 +164,7 @@ interpolate_close_neighbors_on_thread(void *in_prm)
           chstart = (fullind / tl->tottilesinch) * tl->tottilesinch;
 
           /* Set the channel's starting pointer for the flags. */
-          flag = gal_data_ptr_increment(fullflag, chstart, GAL_TYPE_UINT8);
+          flag = gal_pointer_increment(fullflag, chstart, GAL_TYPE_UINT8);
         }
       else
         {
@@ -200,10 +201,10 @@ interpolate_close_neighbors_on_thread(void *in_prm)
               tin=input;
               for(tnear=nearest; tnear!=NULL; tnear=tnear->next)
                 {
-                  memcpy(gal_data_ptr_increment(tnear->array, ngb_counter,
-                                                tin->type),
-                         gal_data_ptr_increment(tin->array, chstart+pind,
-                                                tin->type),
+                  memcpy(gal_pointer_increment(tnear->array, ngb_counter,
+                                               tin->type),
+                         gal_pointer_increment(tin->array, chstart+pind,
+                                               tin->type),
                          gal_type_sizeof(tin->type));
                   tin=tin->next;
                 }
@@ -257,7 +258,7 @@ interpolate_close_neighbors_on_thread(void *in_prm)
         {
           /* Find the median and copy it. */
           median=gal_statistics_median(tnear, 1);
-          memcpy(gal_data_ptr_increment(tout->array, fullind, tout->type),
+          memcpy(gal_pointer_increment(tout->array, fullind, tout->type),
                  median->array, gal_type_sizeof(tout->type));
 
           /* Clean up and go to next array. */
@@ -377,8 +378,8 @@ gal_interpolate_close_neighbors(gal_data_t *input,
                          input->wcs, 0, input->minmapsize, NULL,
                          input->unit, NULL);
   gal_list_void_add(&prm.ngb_vals,
-                    gal_data_malloc_array(input->type, ngbvnum, __func__,
-                                          "prm.ngb_vals"));
+                    gal_pointer_allocate(input->type, ngbvnum, 0, __func__,
+                                         "prm.ngb_vals"));
 
 
   /* If we are given a list of datasets, make the necessary
@@ -390,7 +391,7 @@ gal_interpolate_close_neighbors(gal_data_t *input,
     for(tin=input->next; tin!=NULL; tin=tin->next)
       {
         /* A small sanity check. */
-        if( gal_data_dsize_is_different(input, tin) )
+        if( gal_dimension_is_different(input, tin) )
           error(EXIT_FAILURE, 0, "%s: all datasets in the list must have "
                 "the same dimension and size", __func__);
 
@@ -401,8 +402,8 @@ gal_interpolate_close_neighbors(gal_data_t *input,
 
         /* Allocate the space for the neighbor values of this input. */
         gal_list_void_add(&prm.ngb_vals,
-                          gal_data_malloc_array(tin->type, ngbvnum, __func__,
-                                                "prm.ngb_vals"));
+                          gal_pointer_allocate(tin->type, ngbvnum, 0,
+                                               __func__, "prm.ngb_vals"));
       }
   gal_list_data_reverse(&prm.out);
   gal_list_void_reverse(&prm.ngb_vals);
@@ -410,9 +411,9 @@ gal_interpolate_close_neighbors(gal_data_t *input,
 
   /* Allocate space for all the flag values of all the threads here (memory
      in each thread is limited) and this is cleaner. */
-  prm.thread_flags=gal_data_malloc_array(GAL_TYPE_UINT8,
-                                         numthreads*input->size, __func__,
-                                         "prm.thread_flags");
+  prm.thread_flags=gal_pointer_allocate(GAL_TYPE_UINT8,
+                                        numthreads*input->size, 0, __func__,
+                                        "prm.thread_flags");
 
 
   /* Spin off the threads. */
diff --git a/lib/label.c b/lib/label.c
index 44f7b54..4fd8e8d 100644
--- a/lib/label.c
+++ b/lib/label.c
@@ -31,6 +31,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/list.h>
 #include <gnuastro/qsort.h>
 #include <gnuastro/label.h>
+#include <gnuastro/pointer.h>
 #include <gnuastro/dimension.h>
 #include <gnuastro/statistics.h>
 
@@ -107,7 +108,8 @@ gal_label_indexs(gal_data_t *labels, size_t numlabs, size_t 
minmapsize)
      to allocate). If blank values are present, an extra check is
      necessary, so to get faster results when there aren't any blank
      values, we'll also do a check. */
-  areas=gal_data_calloc_array(GAL_TYPE_SIZE_T, numlabs+1, __func__, "areas");
+  areas=gal_pointer_allocate(GAL_TYPE_SIZE_T, numlabs+1, 1, __func__,
+                             "areas");
   lf=(l=labels->array)+labels->size;
   do
     if(*l>0)  /* Only labeled regions: *l==0 (undetected), *l<0 (blank). */
@@ -192,7 +194,7 @@ gal_label_oversegment(gal_data_t *values, gal_data_t 
*indexs,
   label_check_type(values, GAL_TYPE_FLOAT32, "values", __func__);
   label_check_type(indexs, GAL_TYPE_SIZE_T,  "indexs", __func__);
   label_check_type(labels, GAL_TYPE_INT32,   "labels", __func__);
-  if( gal_data_dsize_is_different(values, labels) )
+  if( gal_dimension_is_different(values, labels) )
     error(EXIT_FAILURE, 0, "%s: the `values' and `labels' arguments must "
           "have the same size", __func__);
   if(indexs->ndim!=1)
@@ -234,13 +236,18 @@ gal_label_oversegment(gal_data_t *values, gal_data_t 
*indexs,
   if(indexs->size==0) return 0;
 
 
-  /* Sort the given indexs based on their flux (`gal_qsort_index_arr' is
+  /* If the indexs aren't already sorted (by the value they correspond to),
+     sort them given indexs based on their flux (`gal_qsort_index_arr' is
      defined as static in `gnuastro/qsort.h') */
-  gal_qsort_index_arr=values->array;
-  qsort(indexs->array, indexs->size, sizeof(size_t),
-        min0_max1
-        ? gal_qsort_index_float_decreasing
-        : gal_qsort_index_float_increasing );
+  if( !( (indexs->flag & GAL_DATA_FLAG_SORT_CH)
+        && ( indexs->flag
+             & (GAL_DATA_FLAG_SORTED_I
+                | GAL_DATA_FLAG_SORTED_D) ) ) )
+    {
+      gal_qsort_index_single=values->array;
+      qsort(indexs->array, indexs->size, sizeof(size_t),
+            min0_max1 ? gal_qsort_index_single_d : gal_qsort_index_single_i);
+    }
 
 
   /* Initialize the region we want to over-segment. */
diff --git a/lib/list.c b/lib/list.c
index 9af66fd..4197690 100644
--- a/lib/list.c
+++ b/lib/list.c
@@ -31,6 +31,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 #include <gnuastro/list.h>
 #include <gnuastro/blank.h>
+#include <gnuastro/pointer.h>
 
 #include <gnuastro-internal/checkset.h>
 
@@ -268,7 +269,7 @@ gal_list_i32_to_array(gal_list_i32_t *list, int reverse, 
size_t *num)
 
   if(*num)
     {
-      out=gal_data_malloc_array(GAL_TYPE_SIZE_T, *num, __func__, "out");
+      out=gal_pointer_allocate(GAL_TYPE_SIZE_T, *num, 0, __func__, "out");
 
       i = reverse ? *num-1: 0;
       if(reverse)
@@ -417,7 +418,7 @@ gal_list_sizet_to_array(gal_list_sizet_t *list, int 
reverse, size_t *num)
 
   if(*num)
     {
-      out=gal_data_malloc_array(GAL_TYPE_SIZE_T, *num, __func__, "out");
+      out=gal_pointer_allocate(GAL_TYPE_SIZE_T, *num, 0, __func__, "out");
 
       i = reverse ? *num-1: 0;
       if(reverse)
@@ -572,7 +573,7 @@ gal_list_f32_to_array(gal_list_f32_t *list, int reverse, 
size_t *num)
   if(*num)
     {
       /* Allocate the space: */
-      out=gal_data_malloc_array(GAL_TYPE_FLOAT32, *num, __func__, "out");
+      out=gal_pointer_allocate(GAL_TYPE_FLOAT32, *num, 0, __func__, "out");
 
       /* Fill in the array. */
       i = reverse ? *num-1: 0;
@@ -729,7 +730,7 @@ gal_list_f64_to_array(gal_list_f64_t *list, int reverse, 
size_t *num)
   if(*num)
     {
       /* Allocate the space: */
-      out=gal_data_malloc_array(GAL_TYPE_FLOAT64, *num, __func__, "out");
+      out=gal_pointer_allocate(GAL_TYPE_FLOAT64, *num, 0, __func__, "out");
 
       /* Fill in the array. */
       i = reverse ? *num-1: 0;
diff --git a/lib/match.c b/lib/match.c
index f156339..1f7b6dd 100644
--- a/lib/match.c
+++ b/lib/match.c
@@ -31,6 +31,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 #include <gnuastro/box.h>
 #include <gnuastro/list.h>
+#include <gnuastro/pointer.h>
 #include <gnuastro/permutation.h>
 
 
@@ -217,8 +218,8 @@ static size_t *
 match_coordinates_prepare_sort(gal_data_t *coords, size_t minmapsize)
 {
   gal_data_t *tmp;
-  size_t *permutation=gal_data_malloc_array(GAL_TYPE_SIZE_T, coords->size,
-                                            __func__, "permutation");
+  size_t *permutation=gal_pointer_allocate(GAL_TYPE_SIZE_T, coords->size, 0,
+                                           __func__, "permutation");
 
   /* Get the permutation necessary to sort all the columns (based on the
      first column). */
@@ -765,8 +766,8 @@ gal_match_coordinates_output(gal_data_t *A, gal_data_t *B, 
size_t *A_perm,
   /* Allocate the `Bmatched' array which is a flag for which rows of the
      second catalog were matched. The columns that had a match will get a
      value of one while we are parsing them below. */
-  Bmatched=gal_data_calloc_array(GAL_TYPE_UINT8, B->size, __func__,
-                                 "Bmatched");
+  Bmatched=gal_pointer_allocate(GAL_TYPE_UINT8, B->size, 1, __func__,
+                                "Bmatched");
 
 
   /* Initialize the indexs. We want the first `nummatched' indexs in both
diff --git a/lib/options.c b/lib/options.c
index c1df2fb..ef10f35 100644
--- a/lib/options.c
+++ b/lib/options.c
@@ -35,6 +35,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/table.h>
 #include <gnuastro/blank.h>
 #include <gnuastro/threads.h>
+#include <gnuastro/pointer.h>
 #include <gnuastro/arithmetic.h>
 
 #include <gnuastro-internal/timing.h>
@@ -814,7 +815,8 @@ gal_options_parse_sizes_reverse(struct argp_option *option, 
char *arg,
       /* Write the values into an allocated size_t array and finish it with
          a `-1' so the total number can be found later.*/
       num=values->size;
-      array=gal_data_malloc_array(GAL_TYPE_SIZE_T, num+1, __func__, "array");
+      array=gal_pointer_allocate(GAL_TYPE_SIZE_T, num+1, 0, __func__,
+                                 "array");
       for(i=0;i<num;++i) array[num-1-i]=v[i];
       array[num] = GAL_BLANK_SIZE_T;
 
diff --git a/lib/permutation.c b/lib/permutation.c
index 3885e65..363cd6a 100644
--- a/lib/permutation.c
+++ b/lib/permutation.c
@@ -26,6 +26,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <string.h>
 #include <stdlib.h>
 
+#include <gnuastro/pointer.h>
 #include <gnuastro/permutation.h>
 
 
@@ -100,7 +101,7 @@ gal_permutation_apply(gal_data_t *input, size_t 
*permutation)
     {
       /* Necessary initializations. */
       width=gal_type_sizeof(input->type);
-      tmp=gal_data_malloc_array(input->type, 1, __func__, "tmp");
+      tmp=gal_pointer_allocate(input->type, 1, 0, __func__, "tmp");
 
       /* Do the permutation. */
       for(i=0;i<input->size;++i)
@@ -150,8 +151,8 @@ gal_permutation_apply_inverse(gal_data_t *input, size_t 
*permutation)
     {
       /* Initializations */
       width=gal_type_sizeof(input->type);
-      tmp=gal_data_malloc_array(input->type, 1, __func__, "tmp");
-      ttmp=gal_data_malloc_array(input->type, 1, __func__, "ttmp");
+      tmp=gal_pointer_allocate(input->type, 1, 0, __func__, "tmp");
+      ttmp=gal_pointer_allocate(input->type, 1, 0, __func__, "ttmp");
 
       /* Re-order the values. */
       for(i=0;i<input->size;++i)
diff --git a/lib/pointer.c b/lib/pointer.c
new file mode 100644
index 0000000..9727874
--- /dev/null
+++ b/lib/pointer.c
@@ -0,0 +1,178 @@
+/*********************************************************************
+pointer -- facilitate working with pointers and allocation.
+This is part of GNU Astronomy Utilities (Gnuastro) package.
+
+Original author:
+     Mohammad Akhlaghi <mohammad@akhlaghi.org>
+Contributing author(s):
+Copyright (C) 2016-2018, Free Software Foundation, Inc.
+
+Gnuastro is free software: you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation, either version 3 of the License, or (at your
+option) any later version.
+
+Gnuastro is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
+**********************************************************************/
+#include <config.h>
+
+#include <stdio.h>
+#include <errno.h>
+#include <error.h>
+#include <unistd.h>
+#include <stdlib.h>
+#include <string.h>
+#include <sys/mman.h>
+
+#include <gnuastro/type.h>
+#include <gnuastro/pointer.h>
+
+#include <gnuastro-internal/checkset.h>
+
+
+/* Increment a give pointer depending on the given type.
+
+   When working with the `array' elements of `gal_data_t', we are actually
+   dealing with `void *' pointers. Pointer arithmetic doesn't apply to
+   `void *', because the system doesn't know how much space each element
+   has to increment the pointer respectively.
+
+   So, here, we will use the type information to find the increment. This
+   is mainly useful when dealing with the `block' pointer of a tile over a
+   larger image. This function reads the address as a `char *' type (note
+   that `char' is guaranteed to have a size of 1 (byte)). It then
+   increments the `char *' by `increment*sizeof(type)' */
+void *
+gal_pointer_increment(void *pointer, size_t increment, uint8_t type)
+{
+  char *p=(char *)pointer;
+  return p + increment * gal_type_sizeof(type);
+}
+
+
+
+
+
+/* Find the number of values between two void pointers with a given
+   type. See the explanations before `gal_data_ptr_increment'. */
+size_t
+gal_pointer_num_between(void *earlier, void *later, uint8_t type)
+{
+  char *e=(char *)earlier, *l=(char *)later;
+  return (l-e)/gal_type_sizeof(type);
+}
+
+
+
+
+
+/* Allocate an array based on the value of type. Note that the argument
+   `size' is the number of elements, necessary in the array, the number of
+   bytes each element needs will be determined internaly by this function
+   using the datatype argument, so you don't have to worry about it. */
+void *
+gal_pointer_allocate(uint8_t type, size_t size, int clear,
+                     const char *funcname, const char *varname)
+{
+  void *array;
+
+  errno=0;
+  array = ( clear
+            ? calloc( size,  gal_type_sizeof(type) )
+            : malloc( size * gal_type_sizeof(type) ) );
+  if(array==NULL)
+    {
+      if(varname)
+        error(EXIT_FAILURE, errno, "%s: %zu bytes couldn't be allocated "
+              "for variable `%s'", funcname ? funcname : __func__,
+              size * gal_type_sizeof(type), varname);
+      else
+        error(EXIT_FAILURE, errno, "%s: %zu bytes couldn't be allocated",
+              funcname ? funcname : __func__, size * gal_type_sizeof(type));
+    }
+
+  return array;
+}
+
+
+
+
+
+void *
+gal_pointer_allocate_mmap(uint8_t type, size_t size, int clear,
+                          char **filename)
+{
+  void *out;
+  int filedes;
+  uint8_t uc=0;
+  size_t bsize=size*gal_type_sizeof(type);
+
+
+  /* Check if the .gnuastro folder exists, write the file there. If it
+     doesn't exist, then make the .gnuastro directory.*/
+  gal_checkset_mkdir(".gnuastro");
+
+
+  /* Set the filename */
+  gal_checkset_allocate_copy("./.gnuastro/mmap_XXXXXX", filename);
+
+
+  /* Create a zero-sized file and keep its descriptor.  */
+  errno=0;
+  /*filedes=open(filename, O_RDWR | O_CREAT | O_EXCL | O_TRUNC );*/
+  filedes=mkstemp(*filename);
+  if(filedes==-1)
+    error(EXIT_FAILURE, errno, "%s: %s couldn't be created",
+          __func__, *filename);
+
+
+  /* Make the necessary space on the file. */
+  errno=0;
+  if( lseek(filedes, bsize, SEEK_SET) == -1 )
+    error(EXIT_FAILURE, errno, "%s: %s: unable to change file position by "
+          "%zu bytes", __func__, *filename, bsize);
+
+
+  /* Write to the newly set file position so the space is allocated. To do
+     this, we are simply writing `uc' (a byte with value 0) into the space
+     we identified by `lseek' (above). This will ensure that this space is
+     set a side for this array and prepare us to use `mmap'. */
+  if( write(filedes, &uc, 1) == -1)
+    error(EXIT_FAILURE, errno, "%s: %s: unable to write one byte at the "
+          "%zu-th position", __func__, *filename, bsize);
+
+
+  /* Map the memory. */
+  errno=0;
+  out=mmap(NULL, bsize, PROT_READ | PROT_WRITE, MAP_SHARED, filedes, 0);
+  if(out==MAP_FAILED)
+    {
+      fprintf(stderr, "\n%s: WARNING: the following error may be due to "
+              "many mmap allocations. Recall that the kernel only allows "
+              "finite number of mmap allocations. It is recommended to use "
+              "ordinary RAM allocation for smaller arrays and keep mmap'd "
+              "allocation only for the large volumes.\n\n", __func__);
+      error(EXIT_FAILURE, errno, "couldn't map %zu bytes into the file `%s'",
+            bsize, *filename);
+    }
+
+
+  /* Close the file. */
+  if( close(filedes) == -1 )
+    error(EXIT_FAILURE, errno, "%s: %s couldn't be closed",
+          __func__, *filename);
+
+
+  /* If it was supposed to be cleared, then clear the memory. */
+  if(clear) memset(out, 0, bsize);
+
+
+  /* Return the mmap'd pointer and save the file name. */
+  return out;
+}
diff --git a/lib/qsort.c b/lib/qsort.c
index 9800b5c..59e5a59 100644
--- a/lib/qsort.c
+++ b/lib/qsort.c
@@ -1,5 +1,5 @@
 /*********************************************************************
-forqsort -- Functions used by qsort to sort an array.
+qsort -- Functions used by qsort to sort an array.
 This is part of GNU Astronomy Utilities (Gnuastro) package.
 
 Original author:
@@ -29,133 +29,112 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 #include <gnuastro/qsort.h>
 
-/* Initialize the array for sorting indexs to NULL. */
-float *gal_qsort_index_arr;
-
-int
-gal_qsort_index_float_decreasing(const void * a, const void * b)
-{
-  float ta=gal_qsort_index_arr[ *(size_t *)a ];
-  float tb=gal_qsort_index_arr[ *(size_t *)b ];
-  return (tb > ta) - (tb < ta);
-}
-
-int
-gal_qsort_index_float_increasing(const void * a, const void * b)
-{
-  float ta=gal_qsort_index_arr[ *(size_t *)a ];
-  float tb=gal_qsort_index_arr[ *(size_t *)b ];
-  return (ta > tb) - (ta < tb);
-}
-
-
-
-
-
-
 
 
 
 
+/*****************************************************************/
+/**********        Sorting of actual array        ****************/
+/*****************************************************************/
 int
-gal_qsort_uint8_decreasing(const void *a, const void *b)
+gal_qsort_uint8_d(const void *a, const void *b)
 {
   return ( *(uint8_t *)b - *(uint8_t *)a );
 }
 
 int
-gal_qsort_uint8_increasing(const void *a, const void *b)
+gal_qsort_uint8_i(const void *a, const void *b)
 {
   return ( *(uint8_t *)a - *(uint8_t *)b );
 }
 
 int
-gal_qsort_int8_decreasing(const void *a, const void *b)
+gal_qsort_int8_d(const void *a, const void *b)
 {
   return ( *(int8_t *)b - *(int8_t *)a );
 }
 
 int
-gal_qsort_int8_increasing(const void *a, const void *b)
+gal_qsort_int8_i(const void *a, const void *b)
 {
   return ( *(int8_t *)a - *(int8_t *)b );
 }
 
 int
-gal_qsort_uint16_decreasing(const void *a, const void *b)
+gal_qsort_uint16_d(const void *a, const void *b)
 {
   return ( *(uint16_t *)b - *(uint16_t *)a );
 }
 
 int
-gal_qsort_uint16_increasing(const void *a, const void *b)
+gal_qsort_uint16_i(const void *a, const void *b)
 {
   return ( *(uint16_t *)a - *(uint16_t *)b );
 }
 
 int
-gal_qsort_int16_decreasing(const void *a, const void *b)
+gal_qsort_int16_d(const void *a, const void *b)
 {
   return ( *(int16_t *)b - *(int16_t *)a );
 }
 
 int
-gal_qsort_int16_increasing(const void *a, const void *b)
+gal_qsort_int16_i(const void *a, const void *b)
 {
   return ( *(int16_t *)a - *(int16_t *)b );
 }
 
 int
-gal_qsort_uint32_decreasing(const void *a, const void *b)
+gal_qsort_uint32_d(const void *a, const void *b)
 {
   return ( *(uint32_t *)b - *(uint32_t *)a );
 }
 
 int
-gal_qsort_uint32_increasing(const void *a, const void *b)
+gal_qsort_uint32_i(const void *a, const void *b)
 {
   return ( *(uint32_t *)a - *(uint32_t *)b );
 }
 
 int
-gal_qsort_int32_decreasing(const void *a, const void *b)
+gal_qsort_int32_d(const void *a, const void *b)
 {
   return ( *(int32_t *)b - *(int32_t *)a );
 }
 
 int
-gal_qsort_int32_increasing(const void *a, const void *b)
+gal_qsort_int32_i(const void *a, const void *b)
 {
   return ( *(int32_t *)a - *(int32_t *)b );
 }
 
 int
-gal_qsort_uint64_decreasing(const void *a, const void *b)
+gal_qsort_uint64_d(const void *a, const void *b)
 {
   return ( *(uint64_t *)b - *(uint64_t *)a );
 }
 
 int
-gal_qsort_uint64_increasing(const void *a, const void *b)
+gal_qsort_uint64_i(const void *a, const void *b)
 {
   return ( *(uint64_t *)a - *(uint64_t *)b );
 }
 
 
 int
-gal_qsort_int64_decreasing(const void *a, const void *b)
+gal_qsort_int64_d(const void *a, const void *b)
 {
   return ( *(int64_t *)b - *(int64_t *)a );
 }
 
 int
-gal_qsort_int64_increasing(const void *a, const void *b)
+gal_qsort_int64_i(const void *a, const void *b)
 {
   return ( *(int64_t *)a - *(int64_t *)b );
 }
 
 int
-gal_qsort_float32_decreasing(const void *a, const void *b)
+gal_qsort_float32_d(const void *a, const void *b)
 {
   float ta=*(float*)a;
   float tb=*(float*)b;
@@ -163,7 +142,7 @@ gal_qsort_float32_decreasing(const void *a, const void *b)
 }
 
 int
-gal_qsort_float32_increasing(const void *a, const void *b)
+gal_qsort_float32_i(const void *a, const void *b)
 {
   float ta=*(float*)a;
   float tb=*(float*)b;
@@ -171,7 +150,7 @@ gal_qsort_float32_increasing(const void *a, const void *b)
 }
 
 int
-gal_qsort_float64_decreasing(const void *a, const void *b)
+gal_qsort_float64_d(const void *a, const void *b)
 {
   double ta=*(double*)a;
   double tb=*(double*)b;
@@ -179,9 +158,80 @@ gal_qsort_float64_decreasing(const void *a, const void *b)
 }
 
 int
-gal_qsort_float64_increasing(const void *a, const void *b)
+gal_qsort_float64_i(const void *a, const void *b)
 {
   double ta=*(double*)a;
   double tb=*(double*)b;
   return (ta > tb) - (ta < tb);
 }
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/*****************************************************************/
+/***************          Sorting indexs        ******************/
+/*****************************************************************/
+/* Initialize the array for sorting indexs to NULL. */
+float *gal_qsort_index_single=NULL;
+
+int
+gal_qsort_index_single_d(const void *a, const void *b)
+{
+  float ta=gal_qsort_index_single[ *(size_t *)a ];
+  float tb=gal_qsort_index_single[ *(size_t *)b ];
+  return (tb > ta) - (tb < ta);
+}
+
+int
+gal_qsort_index_single_i(const void *a, const void *b)
+{
+  float ta=gal_qsort_index_single[ *(size_t *)a ];
+  float tb=gal_qsort_index_single[ *(size_t *)b ];
+  return (ta > tb) - (ta < tb);
+}
+
+int
+gal_qsort_index_multi_d(const void *a, const void *b)
+{
+  /* Define the structures. */
+  struct gal_qsort_index_multi *A = (struct gal_qsort_index_multi *)a;
+  struct gal_qsort_index_multi *B = (struct gal_qsort_index_multi *)b;
+
+  /* For easy reading. */
+  float ta=A->values[ A->index ];
+  float tb=B->values[ B->index ];
+
+  /* Return the result. */
+  return (tb > ta) - (tb < ta);
+}
+
+int
+gal_qsort_index_multi_i(const void *a, const void *b)
+{
+  /* Define the structures. */
+  struct gal_qsort_index_multi *A = (struct gal_qsort_index_multi *)a;
+  struct gal_qsort_index_multi *B = (struct gal_qsort_index_multi *)b;
+
+  /* For easy reading. */
+  float ta=A->values[ A->index ];
+  float tb=B->values[ B->index ];
+
+  /* Return the result. */
+  return (ta > tb) - (ta < tb);
+}
diff --git a/lib/statistics.c b/lib/statistics.c
index 34630bc..b476f7e 100644
--- a/lib/statistics.c
+++ b/lib/statistics.c
@@ -36,6 +36,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/fits.h>
 #include <gnuastro/blank.h>
 #include <gnuastro/qsort.h>
+#include <gnuastro/pointer.h>
 #include <gnuastro/arithmetic.h>
 #include <gnuastro/statistics.h>
 
@@ -358,13 +359,13 @@ gal_statistics_quantile(gal_data_t *input, double 
quantile, int inplace)
       /* Write the value at this index into the output. */
       if(index==GAL_BLANK_SIZE_T)
         {
-          blank=gal_data_malloc_array(nbs->type, 1, __func__, "blank");
+          blank=gal_pointer_allocate(nbs->type, 1, 0, __func__, "blank");
           memcpy(out->array, blank, gal_type_sizeof(nbs->type));
           free(blank);
         }
       else
         memcpy(out->array,
-               gal_data_ptr_increment(nbs->array, index, nbs->type),
+               gal_pointer_increment(nbs->array, index, nbs->type),
                gal_type_sizeof(nbs->type));
     }
   else
@@ -997,7 +998,7 @@ gal_statistics_mode(gal_data_t *input, float mirrordist, 
int inplace)
      same type as the input. */
   modeindex = mode_golden_section(&p);
   memcpy( tmptype->array,
-          gal_data_ptr_increment(p.data->array, modeindex, p.data->type),
+          gal_pointer_increment(p.data->array, modeindex, p.data->type),
           gal_type_sizeof(p.data->type) );
 
 
@@ -1263,25 +1264,25 @@ gal_statistics_sort_increasing(gal_data_t *input)
     switch(input->type)
       {
       case GAL_TYPE_UINT8:
-        STATISTICS_SORT(gal_qsort_uint8_increasing);    break;
+        STATISTICS_SORT(gal_qsort_uint8_i);    break;
       case GAL_TYPE_INT8:
-        STATISTICS_SORT(gal_qsort_int8_increasing);     break;
+        STATISTICS_SORT(gal_qsort_int8_i);     break;
       case GAL_TYPE_UINT16:
-        STATISTICS_SORT(gal_qsort_uint16_increasing);   break;
+        STATISTICS_SORT(gal_qsort_uint16_i);   break;
       case GAL_TYPE_INT16:
-        STATISTICS_SORT(gal_qsort_int16_increasing);    break;
+        STATISTICS_SORT(gal_qsort_int16_i);    break;
       case GAL_TYPE_UINT32:
-        STATISTICS_SORT(gal_qsort_uint32_increasing);   break;
+        STATISTICS_SORT(gal_qsort_uint32_i);   break;
       case GAL_TYPE_INT32:
-        STATISTICS_SORT(gal_qsort_int32_increasing);    break;
+        STATISTICS_SORT(gal_qsort_int32_i);    break;
       case GAL_TYPE_UINT64:
-        STATISTICS_SORT(gal_qsort_uint64_increasing);   break;
+        STATISTICS_SORT(gal_qsort_uint64_i);   break;
       case GAL_TYPE_INT64:
-        STATISTICS_SORT(gal_qsort_int64_increasing);    break;
+        STATISTICS_SORT(gal_qsort_int64_i);    break;
       case GAL_TYPE_FLOAT32:
-        STATISTICS_SORT(gal_qsort_float32_increasing);  break;
+        STATISTICS_SORT(gal_qsort_float32_i);  break;
       case GAL_TYPE_FLOAT64:
-        STATISTICS_SORT(gal_qsort_float64_increasing);  break;
+        STATISTICS_SORT(gal_qsort_float64_i);  break;
       default:
         error(EXIT_FAILURE, 0, "%s: type code %d not recognized",
               __func__, input->type);
@@ -1306,25 +1307,25 @@ gal_statistics_sort_decreasing(gal_data_t *input)
     switch(input->type)
       {
       case GAL_TYPE_UINT8:
-        STATISTICS_SORT(gal_qsort_uint8_decreasing);    break;
+        STATISTICS_SORT(gal_qsort_uint8_d);    break;
       case GAL_TYPE_INT8:
-        STATISTICS_SORT(gal_qsort_int8_decreasing);     break;
+        STATISTICS_SORT(gal_qsort_int8_d);     break;
       case GAL_TYPE_UINT16:
-        STATISTICS_SORT(gal_qsort_uint16_decreasing);   break;
+        STATISTICS_SORT(gal_qsort_uint16_d);   break;
       case GAL_TYPE_INT16:
-        STATISTICS_SORT(gal_qsort_int16_decreasing);    break;
+        STATISTICS_SORT(gal_qsort_int16_d);    break;
       case GAL_TYPE_UINT32:
-        STATISTICS_SORT(gal_qsort_uint32_decreasing);   break;
+        STATISTICS_SORT(gal_qsort_uint32_d);   break;
       case GAL_TYPE_INT32:
-        STATISTICS_SORT(gal_qsort_int32_decreasing);    break;
+        STATISTICS_SORT(gal_qsort_int32_d);    break;
       case GAL_TYPE_UINT64:
-        STATISTICS_SORT(gal_qsort_uint64_decreasing);   break;
+        STATISTICS_SORT(gal_qsort_uint64_d);   break;
       case GAL_TYPE_INT64:
-        STATISTICS_SORT(gal_qsort_int64_decreasing);    break;
+        STATISTICS_SORT(gal_qsort_int64_d);    break;
       case GAL_TYPE_FLOAT32:
-        STATISTICS_SORT(gal_qsort_float32_decreasing);  break;
+        STATISTICS_SORT(gal_qsort_float32_d);  break;
       case GAL_TYPE_FLOAT64:
-        STATISTICS_SORT(gal_qsort_float64_decreasing);  break;
+        STATISTICS_SORT(gal_qsort_float64_d);  break;
       default:
         error(EXIT_FAILURE, 0, "%s: type code %d not recognized",
               __func__, input->type);
diff --git a/lib/tableintern.c b/lib/tableintern.c
index b6e0bdb..a014bcb 100644
--- a/lib/tableintern.c
+++ b/lib/tableintern.c
@@ -33,6 +33,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/txt.h>
 #include <gnuastro/blank.h>
 #include <gnuastro/table.h>
+#include <gnuastro/pointer.h>
 
 #include <gnuastro-internal/timing.h>
 #include <gnuastro-internal/checkset.h>
@@ -426,8 +427,8 @@ gal_tableintern_read_blank(gal_data_t *col, char *blank)
      correctly. */
   if( !gal_type_from_string((void **)(&col->array), blank, col->type) )
     {
-      col->dsize=gal_data_malloc_array(GAL_TYPE_SIZE_T, 1, __func__,
-                                       "col->dsize");
+      col->dsize=gal_pointer_allocate(GAL_TYPE_SIZE_T, 1, 0, __func__,
+                                      "col->dsize");
       col->dsize[0]=col->ndim=col->size=1;
     }
 }
diff --git a/lib/tiff.c b/lib/tiff.c
index 05526d3..e9a455d 100644
--- a/lib/tiff.c
+++ b/lib/tiff.c
@@ -35,7 +35,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/data.h>
 #include <gnuastro/list.h>
 #include <gnuastro/tiff.h>
-
+#include <gnuastro/pointer.h>
 
 
 
@@ -433,7 +433,8 @@ tiff_reverse_rows(gal_data_t *out)
   gal_data_t *ch=out;
   size_t c, i, j, numch=gal_list_data_number(out);
   size_t width=out->dsize[1]*gal_type_sizeof(out->type);
-  void *tmp=gal_data_malloc_array(out->type, out->dsize[1], __func__, "tmp");
+  void *tmp=gal_pointer_allocate(out->type, out->dsize[1], 0, __func__,
+                                 "tmp");
 
   /* A small sanity check. */
   if(out->ndim==3)
diff --git a/lib/tile.c b/lib/tile.c
index 71c036b..a37ece9 100644
--- a/lib/tile.c
+++ b/lib/tile.c
@@ -32,6 +32,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/tile.h>
 #include <gnuastro/blank.h>
 #include <gnuastro/threads.h>
+#include <gnuastro/pointer.h>
 #include <gnuastro/convolve.h>
 #include <gnuastro/dimension.h>
 #include <gnuastro/interpolate.h>
@@ -65,7 +66,7 @@ gal_tile_start_coord(gal_data_t *tile, size_t *start_coord)
   else
     {
       /* Calculate the coordinates of the first pixel of the tile. */
-      ind = gal_data_num_between(block->array, tile->array, block->type);
+      ind = gal_pointer_num_between(block->array, tile->array, block->type);
       gal_dimension_index_to_coord(ind, ndim, block->dsize, start_coord);
     }
 }
@@ -97,7 +98,7 @@ gal_tile_start_end_coord(gal_data_t *tile, size_t *start_end, 
int rel_block)
 
   /* Get the starting index. Note that for the type we need the allocated
      block dataset and can't rely on the tiles. */
-  start_ind=gal_data_num_between(block->array, tile->array, block->type);
+  start_ind=gal_pointer_num_between(block->array, tile->array, block->type);
 
   /* Get the coordinates of the starting point relative to the allocated
      block. */
@@ -108,7 +109,8 @@ gal_tile_start_end_coord(gal_data_t *tile, size_t 
*start_end, int rel_block)
   if(host!=block)
     {
       /* Get the host's starting coordinates. */
-      start_ind=gal_data_num_between(block->array, host->array, block->type);
+      start_ind=gal_pointer_num_between(block->array, host->array,
+                                        block->type);
 
       /* Temporarily put the host's coordinates in the place held for the
          ending coordinates. */
@@ -136,17 +138,17 @@ gal_tile_start_end_ind_inclusive(gal_data_t *tile, 
gal_data_t *work,
 {
   gal_data_t *block=gal_tile_block(tile);
   size_t ndim=tile->ndim, *s, *e, *l, *sf;
-  size_t *start_coord = gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim,
-                                              __func__, "start_coord");
-  size_t *end_coord   = gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim,
-                                              __func__, "end_coord");
+  size_t *start_coord = gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0,
+                                             __func__, "start_coord");
+  size_t *end_coord   = gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0,
+                                             __func__, "end_coord");
 
 
   /* The starting index can be found from the distance of the `tile->array'
      pointer and `block->array' pointer. IMPORTANT: with the type of the
      block array.  */
-  start_end_inc[0]=gal_data_num_between(block->array, tile->array,
-                                        block->type);
+  start_end_inc[0]=gal_pointer_num_between(block->array, tile->array,
+                                           block->type);
 
 
   /* To find the end index, we need to know the coordinates of the starting
@@ -187,7 +189,7 @@ gal_tile_start_end_ind_inclusive(gal_data_t *tile, 
gal_data_t *work,
      from. */
   free(end_coord);
   free(start_coord);
-  return gal_data_ptr_increment(work->array, start_end_inc[0], work->type);
+  return gal_pointer_increment(work->array, start_end_inc[0], work->type);
 }
 
 
@@ -249,14 +251,14 @@ gal_tile_series_from_minmax(gal_data_t *block, size_t 
*minmax, size_t number)
       /* Set the size related constants. */
       size = 1;
       tiles[i].ndim  = ndim;
-      tiles[i].dsize = gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim, __func__,
-                                             "tiles[i].dsize");
+      tiles[i].dsize = gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0,
+                                            __func__, "tiles[i].dsize");
       for(d=0;d<ndim;++d) size *= tiles[i].dsize[d] = max[d] - min[d] + 1;
       tiles[i].size  = size;
 
       /* Tile's array pointer. */
       ind=gal_dimension_coord_to_index(ndim, block->dsize, min);
-      tiles[i].array = gal_data_ptr_increment(block->array, ind, block->type);
+      tiles[i].array = gal_pointer_increment(block->array, ind, block->type);
     }
 
   /* For a check (put all the objects in an extension of a test file).
@@ -467,7 +469,7 @@ gal_tile_block_write_const_value(gal_data_t *tilevalues, 
gal_data_t *tilesll,
       /* Set the pointer to use as input. The `if(o)' statement is set
          because GCC 7.1.1 complained about the possiblity of the first
          argument of `memcpy' being NULL. Recall that `o' is a pointer. */
-      in=gal_data_ptr_increment(tilevalues->array, tile_ind++, type);;
+      in=gal_pointer_increment(tilevalues->array, tile_ind++, type);
       GAL_TILE_PARSE_OPERATE( tile, tofill, 1, withblank, {
           if(o) memcpy(o, in, gal_type_sizeof(type));
         } );
@@ -515,11 +517,11 @@ void *
 gal_tile_block_relative_to_other(gal_data_t *tile, gal_data_t *other)
 {
   gal_data_t *block=gal_tile_block(tile);
-  return gal_data_ptr_increment(other->array,
-                                gal_data_num_between(block->array,
-                                                     tile->array,
-                                                     block->type),
-                                other->type);
+  return gal_pointer_increment(other->array,
+                               gal_pointer_num_between(block->array,
+                                                       tile->array,
+                                                       block->type),
+                               other->type);
 }
 
 
@@ -737,16 +739,16 @@ gal_tile_full(gal_data_t *input, size_t *regular,
 {
   size_t i, d, tind, numtiles, *start=NULL;
   gal_data_t *tiles, *block=gal_tile_block(input);
-  size_t *last   = gal_data_malloc_array(GAL_TYPE_SIZE_T, input->ndim,
-                                         __func__, "last");
-  size_t *first  = gal_data_malloc_array(GAL_TYPE_SIZE_T, input->ndim,
-                                         __func__, "first");
-  size_t *coord  = gal_data_malloc_array(GAL_TYPE_SIZE_T, input->ndim,
-                                         __func__, "coord");
-  size_t *tcoord = gal_data_malloc_array(GAL_TYPE_SIZE_T, input->ndim,
-                                         __func__, "tcoord");
-  size_t *tsize  = gal_data_malloc_array(GAL_TYPE_SIZE_T, input->ndim+1,
-                                         __func__, "tsize");
+  size_t *last   = gal_pointer_allocate(GAL_TYPE_SIZE_T, input->ndim, 0,
+                                      __func__, "last");
+  size_t *first  = gal_pointer_allocate(GAL_TYPE_SIZE_T, input->ndim, 0,
+                                      __func__, "first");
+  size_t *coord  = gal_pointer_allocate(GAL_TYPE_SIZE_T, input->ndim, 0,
+                                      __func__, "coord");
+  size_t *tcoord = gal_pointer_allocate(GAL_TYPE_SIZE_T, input->ndim, 0,
+                                      __func__, "tcoord");
+  size_t *tsize  = gal_pointer_allocate(GAL_TYPE_SIZE_T, input->ndim+1, 0,
+                                      __func__, "tsize");
 
 
   /* Set the first tile size and total number of tiles along each
@@ -765,8 +767,8 @@ gal_tile_full(gal_data_t *input, size_t *regular,
      the block's dimensions when calculating the position of this block. */
   if(input->block)
     {
-      start=gal_data_malloc_array(GAL_TYPE_SIZE_T, input->ndim, __func__,
-                                  "start");
+      start=gal_pointer_allocate(GAL_TYPE_SIZE_T, input->ndim, 0, __func__,
+                                 "start");
       gal_tile_start_coord(input, start);
     }
 
@@ -804,14 +806,14 @@ gal_tile_full(gal_data_t *input, size_t *regular,
       /* Now that we have the index of this tile's starting point compared
          to the allocated block, put it in to the tile's `array'
          pointer. */
-      tiles[i].array=gal_data_ptr_increment(block->array, tind, block->type);
+      tiles[i].array=gal_pointer_increment(block->array, tind, block->type);
 
       /* Set the sizes of the tile. */
       tiles[i].size=1; /* Just an initializer, will be changed. */
       tiles[i].ndim=input->ndim;
       tiles[i].minmapsize=input->minmapsize;
-      tiles[i].dsize=gal_data_malloc_array(GAL_TYPE_SIZE_T,input->ndim,
-                                           __func__, "tiles[i].dsize");
+      tiles[i].dsize=gal_pointer_allocate(GAL_TYPE_SIZE_T,input->ndim, 0,
+                                          __func__, "tiles[i].dsize");
       for(d=0;d<input->ndim;++d)
         {
           /* The size of the first and last tiles can be different from the
@@ -907,8 +909,8 @@ gal_tile_full_sanity_check(char *filename, char *hdu, 
gal_data_t *input,
 
 
   /* Allocate space for the channel sizes. */
-  tl->channelsize=gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim,
-                                        __func__, "tl->channelsize");
+  tl->channelsize=gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__,
+                                       "tl->channelsize");
 
 
   /* Check if the channels are exactly divisible by the input's size along
@@ -1006,8 +1008,8 @@ 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_TYPE_SIZE_T, ndim, __func__,
-                                       "tl->numtiles");
+  tl->numtiles = gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__,
+                                      "tl->numtiles");
   for(i=0;i<ndim;++i)
     tl->numtiles[i] = tl->numtilesinch[i] * tl->numchannels[i];
   tl->tottiles = gal_dimension_total_size(ndim, tl->numtiles);
@@ -1042,11 +1044,12 @@ gal_tile_full_permutation(struct 
gal_tile_two_layer_params *tl)
   if( ndim==1 || tl->totchannels==1) return;
 
   /* Allocate the space for the permutation and coordinates. */
-  ch_coord=gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim, __func__, "ch_coord");
-  tinch_coord=gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim, __func__,
-                                    "tinch_coord");
-  tl->permutation=gal_data_malloc_array(GAL_TYPE_SIZE_T, tl->tottiles,
-                                        __func__, "tl->permutation");
+  ch_coord=gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__,
+                                "ch_coord");
+  tinch_coord=gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__,
+                                 "tinch_coord");
+  tl->permutation=gal_pointer_allocate(GAL_TYPE_SIZE_T, tl->tottiles, 0,
+                                       __func__, "tl->permutation");
 
   /* Fill in the permutation, we use the fact that the tiles are filled
      from the first channel to the last. */
@@ -1146,7 +1149,8 @@ gal_tile_full_values_smooth(gal_data_t *tilevalues,
 
 
   /* Prepare the kernel size along every dimension. */
-  kdsize=gal_data_malloc_array(GAL_TYPE_SIZE_T, tl->ndim, __func__, "kdsize");
+  kdsize=gal_pointer_allocate(GAL_TYPE_SIZE_T, tl->ndim, 0, __func__,
+                              "kdsize");
   for(i=0;i<tl->ndim;++i) kdsize[i]=width;
 
 
diff --git a/lib/txt.c b/lib/txt.c
index 5ea795e..3005d53 100644
--- a/lib/txt.c
+++ b/lib/txt.c
@@ -1276,7 +1276,7 @@ gal_txt_write(gal_data_t *input, gal_list_str_t *comment, 
char *filename)
 
       /* Check if the dimensionality and size is the same for all the
          elements. */
-      if( input!=data && gal_data_dsize_is_different(input, data) )
+      if( input!=data && gal_dimension_is_different(input, data) )
         error(EXIT_FAILURE, 0, "%s: the input list of datasets must have the "
               "same sizes (dimentionality and length along each dimension)",
               __func__);
diff --git a/lib/type.c b/lib/type.c
index f59981d..5d89a72 100644
--- a/lib/type.c
+++ b/lib/type.c
@@ -34,6 +34,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/type.h>
 #include <gnuastro/data.h>
 #include <gnuastro/list.h>
+#include <gnuastro/pointer.h>
 #include <gnuastro-internal/checkset.h>
 
 
@@ -340,7 +341,8 @@ gal_type_bit_string(void *in, size_t size)
 {
   size_t i;
   char *byte=in;
-  char *str=gal_data_malloc_array(GAL_TYPE_UINT8, 8*size+1, __func__, "str");
+  char *str=gal_pointer_allocate(GAL_TYPE_UINT8, 8*size+1, 0, __func__,
+                                 "str");
 
   /* Print the bits into the allocated string. This was inspired from
 
@@ -449,7 +451,7 @@ gal_type_from_string(void **out, char *string, uint8_t type)
   if( *out==NULL && !gal_type_is_list(type) )
     {
       allocated=1;
-      *out=gal_data_malloc_array(type, 1, __func__, "out");
+      *out=gal_pointer_allocate(type, 1, 0, __func__, "out");
     }
   value=*out;
 
@@ -609,7 +611,7 @@ gal_type_string_to_number(char *string, uint8_t *type)
     }
 
   /* Allocate a one-element dataset, then copy the number into it. */
-  out=gal_data_malloc_array(*type, 1, __func__, "out");
+  out=gal_pointer_allocate(*type, 1, 0, __func__, "out");
   memcpy(out, ptr, gal_type_sizeof(*type));
   return out;
 }
diff --git a/lib/wcs.c b/lib/wcs.c
index 3d12fea..47a4420 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -36,6 +36,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/wcs.h>
 #include <gnuastro/tile.h>
 #include <gnuastro/fits.h>
+#include <gnuastro/pointer.h>
 #include <gnuastro/dimension.h>
 #include <gnuastro/permutation.h>
 
@@ -269,8 +270,8 @@ gal_wcs_on_tile(gal_data_t *tile)
 {
   size_t i, start_ind, ndim=tile->ndim;
   gal_data_t *block=gal_tile_block(tile);
-  size_t *coord=gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim, __func__,
-                                      "coord");
+  size_t *coord=gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__,
+                                     "coord");
 
   /* If the tile already has a WCS structure, don't do anything. */
   if(tile->wcs) return;
@@ -280,7 +281,8 @@ gal_wcs_on_tile(gal_data_t *tile)
       tile->wcs=gal_wcs_copy(block->wcs);
 
       /* Find the coordinates of the tile's starting index. */
-      start_ind=gal_data_num_between(block->array, tile->array, block->type);
+      start_ind=gal_pointer_num_between(block->array, tile->array,
+                                        block->type);
       gal_dimension_index_to_coord(start_ind, ndim, block->dsize, coord);
 
       /* Correct the copied WCS structure. Note that crpix is indexed in
@@ -445,9 +447,9 @@ gal_wcs_pixel_scale(struct wcsprm *wcs)
   int warning_printed;
   size_t i, j, maxj, n=wcs->naxis;
   double jvmax, *a, *out, maxrow, minrow;
-  double *v=gal_data_malloc_array(GAL_TYPE_FLOAT64, n*n, __func__, "v");
-  size_t *permutation=gal_data_malloc_array(GAL_TYPE_SIZE_T, n, __func__,
-                                            "permutation");
+  double *v=gal_pointer_allocate(GAL_TYPE_FLOAT64, n*n, 0, __func__, "v");
+  size_t *permutation=gal_pointer_allocate(GAL_TYPE_SIZE_T, n, 0, __func__,
+                                           "permutation");
   gal_data_t *pixscale=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &n, NULL,
                                       0, -1, NULL, NULL, NULL);
 
@@ -665,15 +667,18 @@ wcs_convert_sanity_check_alloc(gal_data_t *coords, struct 
wcsprm *wcs,
           ndim, wcs->naxis);
 
   /* Allocate all the necessary arrays. */
-  *phi    = gal_data_malloc_array( GAL_TYPE_FLOAT64, size, __func__, "phi");
-  *stat   = gal_data_calloc_array( GAL_TYPE_INT32,   size, __func__, "stat");
-  *theta  = gal_data_malloc_array( GAL_TYPE_FLOAT64, size, __func__, "theta");
-  *world  = gal_data_malloc_array( GAL_TYPE_FLOAT64, ndim*size, __func__,
-                                   "world");
-  *imgcrd = gal_data_malloc_array( GAL_TYPE_FLOAT64, ndim*size, __func__,
-                                   "imgcrd");
-  *pixcrd = gal_data_malloc_array( GAL_TYPE_FLOAT64, ndim*size, __func__,
-                                   "pixcrd");
+  *phi    = gal_pointer_allocate( GAL_TYPE_FLOAT64, size,      0, __func__,
+                                  "phi");
+  *stat   = gal_pointer_allocate( GAL_TYPE_INT32,   size,      1, __func__,
+                                  "stat");
+  *theta  = gal_pointer_allocate( GAL_TYPE_FLOAT64, size,      0, __func__,
+                                  "theta");
+  *world  = gal_pointer_allocate( GAL_TYPE_FLOAT64, ndim*size, 0, __func__,
+                                  "world");
+  *imgcrd = gal_pointer_allocate( GAL_TYPE_FLOAT64, ndim*size, 0, __func__,
+                                  "imgcrd");
+  *pixcrd = gal_pointer_allocate( GAL_TYPE_FLOAT64, ndim*size, 0, __func__,
+                                  "pixcrd");
 }
 
 



reply via email to

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