gnuastro-commits
[Top][All Lists]
Advanced

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

[gnuastro-commits] master db39cee 088/125: Statistics returns basic resu


From: Mohammad Akhlaghi
Subject: [gnuastro-commits] master db39cee 088/125: Statistics returns basic results, no array lib, new blank lib
Date: Sun, 23 Apr 2017 22:36:44 -0400 (EDT)

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

    Statistics returns basic results, no array lib, new blank lib
    
    The work on this commit started with merging the Statistics program with
    the new data structure. However, it soon grew to the libraries. The full
    process is described below.
    
    The Statistics program can now read a table column as well as an image for
    input. It can also remove certain ranges of data from its internal analysis
    and as before also ignores blank pixels. Implementation of Statistics is
    not yet complete and its section will be updated in the book after it is
    complete (hopefully in the next commit). This commit was necessary because
    too many changes were made in the libraries.
    
    In order to do things clearly in statistics, many new operators were added
    to the Arithmetic library (and thus the Arithmetic program): `numvalue',
    `sumvalue', `meanvalue', `stdvalue', `medianvalue', `num', `mean' (instead
    of the old `average'), and `std'.
    
    Working with blank pixels until now was done in the main `data.h' library,
    causing a bloating of operations in that library. To be clear and easy to
    understand, the `data.h' is now only concerned with the most low-level
    dataset issues (like types, allocation, copying, or reading and writing
    single elements from strings). So we now have a new `blank.h' library for
    functions related to blank pixels.
    
    The old `array.h' library was a remnant of the past and with the new
    libraries is now obsolete. So it was removed (its files and also its
    sections in the book). Many of the old obsolete statistics functions have
    also been removed.
---
 bin/arithmetic/arithmetic.c       |  18 +-
 bin/convolve/ui.c                 |   5 +-
 bin/crop/onecrop.c                |   5 +-
 bin/crop/ui.c                     |  12 +-
 bin/mkprof/oneprofile.c           |  18 +-
 bin/mkprof/ui.c                   |  10 +-
 bin/statistics/args.h             | 697 ++++++++++------------------
 bin/statistics/aststatistics.conf |  17 +-
 bin/statistics/authors-cite.h     |   2 +-
 bin/statistics/main.c             |   8 +-
 bin/statistics/main.h             | 129 ++----
 bin/statistics/statistics.c       |  53 ++-
 bin/statistics/statistics.h       |   2 +-
 bin/statistics/ui.c               | 953 ++++++++++++++++++--------------------
 bin/statistics/ui.h               |  11 +-
 bin/warp/ui.c                     | 107 +----
 doc/gnuastro.texi                 | 490 +++-----------------
 lib/Makefile.am                   |  14 +-
 lib/arithmetic.c                  | 377 ++++++++++++---
 lib/array.c                       | 530 ---------------------
 lib/blank.c                       | 402 ++++++++++++++++
 lib/data.c                        | 430 +----------------
 lib/fits.c                        |  35 +-
 lib/gnuastro/arithmetic.h         |  30 +-
 lib/gnuastro/array.h              | 157 -------
 lib/gnuastro/blank.h              |  91 ++++
 lib/gnuastro/data.h               |  52 +--
 lib/gnuastro/qsort.h              |  46 +-
 lib/gnuastro/table.h              |  35 +-
 lib/mesh.c                        |   4 +-
 lib/options.c                     | 125 ++++-
 lib/options.h                     |   7 +-
 lib/qsort.c                       |  85 ++--
 lib/statistics.c                  | 783 +------------------------------
 lib/table.c                       | 118 +++--
 lib/txt.c                         |  33 +-
 tests/Makefile.am                 |   3 +-
 tests/lib/arraymanip.c            |  75 ---
 tmpfs-config-make                 |  10 +-
 39 files changed, 2103 insertions(+), 3876 deletions(-)

diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index 7ff824d..5c1e123 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -188,14 +188,28 @@ reversepolish(struct imgarithparams *p)
             { op=GAL_ARITHMETIC_OP_MINVAL;        nop=1;  }
           else if (!strcmp(token->v, "maxvalue"))
             { op=GAL_ARITHMETIC_OP_MAXVAL;        nop=1;  }
+          else if (!strcmp(token->v, "numvalue"))
+            { op=GAL_ARITHMETIC_OP_NUMVAL;        nop=1;  }
+          else if (!strcmp(token->v, "sumvalue"))
+            { op=GAL_ARITHMETIC_OP_SUMVAL;        nop=1;  }
+          else if (!strcmp(token->v, "meanvalue"))
+            { op=GAL_ARITHMETIC_OP_MEANVAL;       nop=1;  }
+          else if (!strcmp(token->v, "stdvalue"))
+            { op=GAL_ARITHMETIC_OP_STDVAL;        nop=1;  }
+          else if (!strcmp(token->v, "medianvalue"))
+            { op=GAL_ARITHMETIC_OP_MEDIANVAL;     nop=1;  }
           else if (!strcmp(token->v, "min"))
             { op=GAL_ARITHMETIC_OP_MIN;           nop=-1; }
           else if (!strcmp(token->v, "max"))
             { op=GAL_ARITHMETIC_OP_MAX;           nop=-1; }
+          else if (!strcmp(token->v, "num"))
+            { op=GAL_ARITHMETIC_OP_NUM;           nop=-1; }
           else if (!strcmp(token->v, "sum"))
             { op=GAL_ARITHMETIC_OP_SUM;           nop=-1; }
-          else if (!strcmp(token->v, "average"))
-            { op=GAL_ARITHMETIC_OP_AVERAGE;       nop=-1; }
+          else if (!strcmp(token->v, "mean"))
+            { op=GAL_ARITHMETIC_OP_MEAN;          nop=-1; }
+          else if (!strcmp(token->v, "std"))
+            { op=GAL_ARITHMETIC_OP_STD;           nop=-1; }
           else if (!strcmp(token->v, "median"))
             { op=GAL_ARITHMETIC_OP_MEDIAN;        nop=-1; }
 
diff --git a/bin/convolve/ui.c b/bin/convolve/ui.c
index 3f8dc34..12b4edd 100644
--- a/bin/convolve/ui.c
+++ b/bin/convolve/ui.c
@@ -30,6 +30,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 #include <gnuastro/wcs.h>
 #include <gnuastro/fits.h>
+#include <gnuastro/blank.h>
 #include <gnuastro/table.h>
 #include <gnuastro/linkedlist.h>
 
@@ -297,7 +298,7 @@ ui_read_kernel(struct convolveparams *p)
 
   /* Convert all the NaN pixels to zero if the kernel contains blank
      pixels. */
-  if(gal_data_has_blank(data))
+  if(gal_blank_present(data))
     {
       ff=(f=data->array)+data->size;
       do *f = isnan(*f) ? 0.0f : *f; while(++f<ff);
@@ -347,7 +348,7 @@ ui_preparations(struct convolveparams *p)
 
 
   /* See if there are any blank values. */
-  if(gal_data_has_blank(data) && p->domain==CONVOLVE_DOMAIN_FREQUENCY)
+  if(gal_blank_present(data) && p->domain==CONVOLVE_DOMAIN_FREQUENCY)
     fprintf(stderr, "\n----------------------------------------\n"
             "######## %s WARNING ########\n"
             "There are blank pixels in `%s' (hdu: `%s') and you have asked "
diff --git a/bin/crop/onecrop.c b/bin/crop/onecrop.c
index e8332bf..3ea6fae 100644
--- a/bin/crop/onecrop.c
+++ b/bin/crop/onecrop.c
@@ -33,6 +33,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 #include <gnuastro/box.h>
 #include <gnuastro/fits.h>
+#include <gnuastro/blank.h>
 #include <gnuastro/polygon.h>
 
 #include <timing.h>
@@ -307,7 +308,7 @@ imgpolygonflpixel(double *ipolygon, size_t nvertices, long 
*fpixel,
 
 
 #define POLYGON_MASK(CTYPE) {                                           \
-    CTYPE *ba=array, *bb=gal_data_alloc_blank(type);                    \
+    CTYPE *ba=array, *bb=gal_blank_alloc_write(type);                   \
     for(i=0;i<size;++i)                                                 \
       {                                                                 \
         point[0]=i%s1+1; point[1]=i/s1+1;                               \
@@ -805,7 +806,7 @@ iscenterfilled(struct onecropparams *crp)
   long naxes[2], fpixel[2], lpixel[2], inc[2]={1,1};
 
   /* If checkcenter is zero, then don't check. */
-  if(checkcenter==0) return GAL_DATA_BLANK_UINT8;
+  if(checkcenter==0) return GAL_BLANK_UINT8;
 
   /* Get the final size of the output image. */
   gal_fits_img_info(ofp, &type, &ndim, &dsize);
diff --git a/bin/crop/ui.c b/bin/crop/ui.c
index 682746e..da5038f 100644
--- a/bin/crop/ui.c
+++ b/bin/crop/ui.c
@@ -30,6 +30,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 #include <gnuastro/wcs.h>
 #include <gnuastro/fits.h>
+#include <gnuastro/blank.h>
 #include <gnuastro/table.h>
 #include <gnuastro/linkedlist.h>
 
@@ -584,7 +585,10 @@ ui_read_cols(struct cropparams *p)
 
       /* Print an error if there were too many columns. */
       if(toomanycols)
-        gal_table_too_many_columns(p->catname);
+        gal_table_error_col_selection(p->catname, p->cathdu, "too many "
+                                      "columns were selected by the given "
+                                      "values to the options ending in "
+                                      "`col'.");
 
       /* Sanity check and clean up.  Note that it might happen that the
          input structure is already freed. In that case, `corrtype' will be
@@ -592,7 +596,7 @@ ui_read_cols(struct cropparams *p)
       if(corrtype)
         {
           /* Make sure there are no blank values in this column. */
-          if( gal_data_has_blank(corrtype) )
+          if( gal_blank_present(corrtype) )
             error(EXIT_FAILURE, 0, "%s column has blank values. "
                   "Input columns cannot contain blank values", colname);
 
@@ -624,7 +628,7 @@ ui_make_log(struct cropparams *p)
 
   /* If central pixels are filled. */
   asprintf(&comment, "Are the central pixels filled? (1: yes, 0: no, "
-           "%u: not checked)", GAL_DATA_BLANK_UINT8);
+           "%u: not checked)", GAL_BLANK_UINT8);
   gal_data_add_to_ll(&p->log, NULL, GAL_DATA_TYPE_UINT8, 1, &p->numout,
                      NULL, 1, p->cp.minmapsize, "CENTER_FILLED", "bool",
                      comment);
@@ -730,7 +734,7 @@ ui_preparations(struct cropparams *p)
       if(firsttype==0)
         {
           firsttype=p->type;
-          p->bitnul = gal_data_alloc_blank(p->type);
+          p->bitnul = gal_blank_alloc_write(p->type);
         }
       else
         {
diff --git a/bin/mkprof/oneprofile.c b/bin/mkprof/oneprofile.c
index 33c3871..9824412 100644
--- a/bin/mkprof/oneprofile.c
+++ b/bin/mkprof/oneprofile.c
@@ -561,9 +561,9 @@ oneprofile_make(struct mkonthread *mkp)
 {
   struct mkprofparams *p=mkp->p;
 
-  float sum;
+  float *f, *ff;
   long os=p->oversample;
-  double pixfrac, intpart;
+  double sum, pixfrac, intpart;
   size_t size, id=mkp->ibq->id;
 
 
@@ -607,7 +607,7 @@ oneprofile_make(struct mkonthread *mkp)
   if(mkp->correction)
     {
       /* First get the sum of all the pixels in the profile. */
-      sum=gal_statistics_float_sum(mkp->ibq->img, size);
+      sum=0.0f; ff=(f=mkp->ibq->img)+size; do sum+=*f++; while(f<ff);
 
       /* Correct the fraction of brightness that was calculated
          accurately (not using the pixel center). */
@@ -615,10 +615,14 @@ oneprofile_make(struct mkonthread *mkp)
 
       /* Correct all the profile pixels. */
       if(p->magatpeak)
-        gal_array_fmultip_const(mkp->ibq->img, size,
-                                mkp->brightness/mkp->peakflux);
+        {
+          ff=(f=mkp->ibq->img)+size;
+          do *f++ *= mkp->brightness/mkp->peakflux; while(f<ff);
+        }
       else
-        gal_array_fmultip_const(mkp->ibq->img, size,
-                                mkp->brightness/sum);
+        {
+          ff=(f=mkp->ibq->img)+size;
+          do *f++ *= mkp->brightness/sum; while(f<ff);
+        }
     }
 }
diff --git a/bin/mkprof/ui.c b/bin/mkprof/ui.c
index 083a26b..f5fb997 100644
--- a/bin/mkprof/ui.c
+++ b/bin/mkprof/ui.c
@@ -31,6 +31,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/wcs.h>
 #include <gnuastro/box.h>
 #include <gnuastro/fits.h>
+#include <gnuastro/blank.h>
 #include <gnuastro/table.h>
 #include <gnuastro/linkedlist.h>
 
@@ -402,7 +403,7 @@ ui_read_profile_function(struct mkprofparams *p, char 
**strarr)
       else if ( !strcmp("circum", strarr[i]) )
         p->f[i]=PROFILE_CIRCUMFERENCE;
 
-      else if ( !strcmp(GAL_DATA_BLANK_STRING, strarr[i]) )
+      else if ( !strcmp(GAL_BLANK_STRING, strarr[i]) )
         error(EXIT_FAILURE, 0, "profile function column has blank values. "
               "Input columns cannot contain blank values");
       else
@@ -543,7 +544,10 @@ ui_read_cols(struct mkprofparams *p)
         /* If the index isn't recognized, then it is larger, showing that
            there was more than one match for the given criteria */
         default:
-          gal_table_too_many_columns(p->catname);
+          gal_table_error_col_selection(p->catname, p->cp.hdu, "too many "
+                                        "columns were selected by the given "
+                                        "values to the options ending in "
+                                        "`col'.");
         }
 
       /* Sanity check and clean up.  Note that it might happen that the
@@ -552,7 +556,7 @@ ui_read_cols(struct mkprofparams *p)
       if(corrtype)
         {
           /* Make sure there are no blank values in this column. */
-          if( checkblank && gal_data_has_blank(corrtype) )
+          if( checkblank && gal_blank_present(corrtype) )
             error(EXIT_FAILURE, 0, "%s column has blank values. "
                   "Input columns cannot contain blank values", colname);
 
diff --git a/bin/statistics/args.h b/bin/statistics/args.h
index f339dcd..7ec8bfd 100644
--- a/bin/statistics/args.h
+++ b/bin/statistics/args.h
@@ -1,5 +1,5 @@
 /*********************************************************************
-Statistics - Get general statistics about the image.
+Statistics - Statistical analysis on input dataset.
 Statistics is part of GNU Astronomy Utilities (Gnuastro) package.
 
 Original author:
@@ -23,102 +23,65 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #ifndef ARGS_H
 #define ARGS_H
 
-#include <argp.h>
 
-#include <gnuastro/linkedlist.h>
 
-#include <commonargs.h>
-#include <fixedstringmacros.h>
 
 
 
-
-
-
-
-
-
-
-/**************************************************************/
-/**************        argp.h definitions       ***************/
-/**************************************************************/
-
-
-
-
-/* Definition parameters for the argp: */
-const char *argp_program_version=SPACK_STRING"\n"GAL_STRINGS_COPYRIGHT
-  "\n\nWritten by Mohammad Akhlaghi";
-const char *argp_program_bug_address=PACKAGE_BUGREPORT;
-static char args_doc[] = "ASTRdata";
-
-
-
-
-
-const char doc[] =
-  /* Before the list of options: */
-  GAL_STRINGS_TOP_HELP_INFO
-  SPACK_NAME" will print the basic statistics of the input image pixel "
-  "flux distribution. All blank pixels or pixels specified by a mask "
-  "image will be ignored.\n"
-  GAL_STRINGS_MORE_HELP_INFO
-  /* After the list of options: */
-  "\v"
-  PACKAGE_NAME" home page: "PACKAGE_URL;
-
-
-
-
-
-/* Available letters for short options:
-
-   c e f j k m s t v w y z
-   C E F G I J L O R T W X Y Z
-
-   Number keys used: <=511
-
-   Options with keys (second structure element) larger than 500 do not
-   have a short version.
- */
-static struct argp_option options[] =
+/* Array of acceptable options. */
+struct argp_option program_options[] =
   {
     {
-      0, 0, 0, 0,
-      "Input:",
-      1
-    },
-    {
-      "mask",
-      'M',
+      "column",
+      ARGS_OPTION_KEY_COLUMN,
       "STR",
       0,
-      "Mask image file name.",
-      1
+      "Column name or number if input is a table.",
+      GAL_OPTIONS_GROUP_INPUT,
+      &p->column,
+      GAL_DATA_TYPE_STRING,
+      GAL_OPTIONS_RANGE_ANY,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
     },
     {
-      "mhdu",
-      'H',
-      "STR",
+      "greaterequal",
+      ARGS_OPTION_KEY_GREATEREQUAL,
+      "FLT",
       0,
-      "Mask image header name.",
-      1
+      "Only use values greater-equal than this.",
+      GAL_OPTIONS_GROUP_INPUT,
+      &p->greaterequal,
+      GAL_DATA_TYPE_FLOAT32,
+      GAL_OPTIONS_RANGE_ANY,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
     },
     {
-      "ignoremin",
-      'r',
-      0,
+      "lessthan",
+      ARGS_OPTION_KEY_LESSTHAN,
+      "FLT",
       0,
-      "Ignore data with values equal to minimum.",
-      1
+      "Only use values less than this.",
+      GAL_OPTIONS_GROUP_INPUT,
+      &p->lessthan,
+      GAL_DATA_TYPE_FLOAT32,
+      GAL_OPTIONS_RANGE_ANY,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
     },
     {
-      "mirrordist",
-      'd',
+      "quantrange",
+      ARGS_OPTION_KEY_QUANTRANGE,
       "FLT",
       0,
-      "Distance beyond mirror point. Multiple of std.",
-      1
+      "Quantile (Q) range: from Q to 1-Q.",
+      GAL_OPTIONS_GROUP_INPUT,
+      &p->quantilerange,
+      GAL_DATA_TYPE_FLOAT32,
+      GAL_OPTIONS_RANGE_ANY,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
     },
 
 
@@ -127,248 +90,270 @@ static struct argp_option options[] =
 
     {
       0, 0, 0, 0,
-      "Output:",
-      2
-    },
-    {
-      "lowerbin",
-      'l',
-      0,
-      0,
-      "Interval lower limit for column 1.",
-      2
-    },
-    {
-      "onebinvalue",
-      'B',
-      "FLT",
-      0,
-      "Shift bins so one bin starts on this value.",
-      2
+      "Values to print in one row",
+      ARGS_GROUP_IN_ONE_ROW
     },
     {
-      "noasciihist",
-      'A',
+      "number",
+      ARGS_OPTION_KEY_NUMBER,
       0,
       0,
-      "Do not show an ASCII histogram of the data.",
-      2
+      "Print the number of data-points.",
+      ARGS_GROUP_IN_ONE_ROW,
+      &p->toprint,
+      GAL_DATA_TYPE_INVALID,
+      GAL_OPTIONS_RANGE_ANY,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET,
+      ui_add_to_print_in_row
     },
     {
-      "checkmode",
-      509,
+      "minimum",
+      ARGS_OPTION_KEY_MINIMUM,
       0,
       0,
-      "Mode mirror plot. `_modehist.txt', `_modecfp.txt'",
-      2
+      "Print the minimum.",
+      ARGS_GROUP_IN_ONE_ROW,
+      &p->toprint,
+      GAL_OPTIONS_NO_ARG_TYPE,
+      GAL_OPTIONS_RANGE_0_OR_1,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET,
+      ui_add_to_print_in_row
     },
     {
-      "mirrorquant",
-      510,
-      "FLT",
+      "maximum",
+      ARGS_OPTION_KEY_MAXIMUM,
       0,
-      "Mirror quantile. `_mirhist.txt', `_mircfp.txt'.",
-      2
-    },
-    {
-      "histrangeformirror",
-      511,
       0,
-      0,
-      "Use input histogram range for mirror plots.",
-      2
-    },
-    {
-      "mirrorplotdist",
-      503,
-      "FLT",
-      0,
-      "Distance beyond mode to display.",
-      2
-    },
-
-
-
-    {
-      0, 0, 0, 0,
-      "Histogram (suffix: `_hist.txt'):",
-      3
+      "Print the maximum.",
+      ARGS_GROUP_IN_ONE_ROW,
+      &p->toprint,
+      GAL_OPTIONS_NO_ARG_TYPE,
+      GAL_OPTIONS_RANGE_0_OR_1,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET,
+      ui_add_to_print_in_row
     },
     {
-      "nohist",
-      500,
+      "sum",
+      ARGS_OPTION_KEY_SUM,
       0,
       0,
-      "Do not calculate histogram.",
-      3
+      "Print the sum of all elements.",
+      ARGS_GROUP_IN_ONE_ROW,
+      &p->toprint,
+      GAL_OPTIONS_NO_ARG_TYPE,
+      GAL_OPTIONS_RANGE_0_OR_1,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET,
+      ui_add_to_print_in_row
     },
     {
-      "normhist",
-      501,
+      "mean",
+      ARGS_OPTION_KEY_MEAN,
       0,
       0,
-      "Normalize the  histogram (sum of all bins 1).",
-      3
+      "Print the mean.",
+      ARGS_GROUP_IN_ONE_ROW,
+      &p->toprint,
+      GAL_OPTIONS_NO_ARG_TYPE,
+      GAL_OPTIONS_RANGE_0_OR_1,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET,
+      ui_add_to_print_in_row
     },
     {
-      "maxhistone",
-      502,
+      "std",
+      ARGS_OPTION_KEY_STD,
       0,
       0,
-      "Scale such that the maximum bin has value of one.",
-      3
+      "Print the standad deviation.",
+      ARGS_GROUP_IN_ONE_ROW,
+      &p->toprint,
+      GAL_OPTIONS_NO_ARG_TYPE,
+      GAL_OPTIONS_RANGE_0_OR_1,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET,
+      ui_add_to_print_in_row
     },
     {
-      "histnumbins",
-      'n',
-      "INT",
+      "median",
+      ARGS_OPTION_KEY_MEDIAN,
       0,
-      "Number of bins in the histogram.",
-      3
-    },
-    {
-      "histmin",
-      'i',
-      "FLT",
       0,
-      "The minimum value for the histogram.",
-      3
+      "Print the median.",
+      ARGS_GROUP_IN_ONE_ROW,
+      &p->toprint,
+      GAL_OPTIONS_NO_ARG_TYPE,
+      GAL_OPTIONS_RANGE_0_OR_1,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET,
+      ui_add_to_print_in_row
     },
     {
-      "histmax",
-      'x',
-      "FLT",
+      "mode",
+      ARGS_OPTION_KEY_MODE,
       0,
-      "The maximum value for the histogram.",
-      3
-    },
-    {
-      "histquant",
-      'Q',
-      "FLT",
       0,
-      "Quantile (Q) range. Histogram from Q to 1-Q.",
-      3
+      "Print the mode (for large datasets).",
+      ARGS_GROUP_IN_ONE_ROW,
+      &p->toprint,
+      GAL_OPTIONS_NO_ARG_TYPE,
+      GAL_OPTIONS_RANGE_0_OR_1,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET,
+      ui_add_to_print_in_row
     },
 
 
 
 
+
     {
       0, 0, 0, 0,
-      "Cumulative Frequency Plot (suffix: `_cfp.txt'):",
-      4
+      "Particular calculation",
+      ARGS_GROUP_PARTICULAR_STAT
     },
     {
-      "nocfp",
-      504,
+      "asciihist",
+      ARGS_OPTION_KEY_ASCIIHIST,
       0,
       0,
-      "No Cumulative Frequency Plot.",
-      4
+      "Print an ASCII histogram",
+      ARGS_GROUP_PARTICULAR_STAT,
+      &p->asciihist,
+      GAL_OPTIONS_NO_ARG_TYPE,
+      GAL_OPTIONS_RANGE_0_OR_1,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
     },
     {
-      "normcfp",
-      505,
+      "histogram",
+      ARGS_OPTION_KEY_HISTOGRAM,
       0,
       0,
-      "Normalize the CFP (sum of all bins 1).",
-      4
+      "Save the histogram in output.",
+      ARGS_GROUP_PARTICULAR_STAT,
+      &p->histogram,
+      GAL_OPTIONS_NO_ARG_TYPE,
+      GAL_OPTIONS_RANGE_0_OR_1,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
     },
     {
-      "maxcfpeqmaxhist",
-      506,
+      "cumulative",
+      ARGS_OPTION_KEY_CUMULATIVE,
       0,
       0,
-      "Set maximum of CFP to maximum of histogram.",
-      4
+      "Save the cumulative frequency plot.",
+      ARGS_GROUP_PARTICULAR_STAT,
+      &p->cumulative,
+      GAL_OPTIONS_NO_ARG_TYPE,
+      GAL_OPTIONS_RANGE_0_OR_1,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
     },
     {
-      "cfpsimhist",
-      507,
-      0,
+      "sigmaclip",
+      ARGS_OPTION_KEY_SIGMACLIP,
+      "FLT,FLT",
       0,
-      "Set CFP range and bins similar to histogram.",
-      4
+      "Multiple of sigma and tolerance or number.",
+      ARGS_GROUP_PARTICULAR_STAT,
+      &p->sigclipstr,
+      GAL_DATA_TYPE_STRING,
+      GAL_OPTIONS_RANGE_ANY,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
     },
     {
-      "cfpnum",
-      'p',
-      "INT",
-      0,
-      "Number of data points to find CFP.",
-      4
-    },
-    {
-      "cfpmin",
-      'a',
-      "FLT",
-      0,
-      "Minimum value to use in the CFP.",
-      4
-    },
-    {
-      "cfpmax",
-      'b',
+      "mirror",
+      ARGS_OPTION_KEY_MIRROR,
       "FLT",
       0,
-      "Maximum value to use in the CFP.",
-      4
-    },
-    {
-      "cfpquant",
-      'U',
-      "FLT",
-      0,
-      "Quantile of range: from U to 1-U.",
-      4
+      "Quantile of mirror distribution to save.",
+      ARGS_GROUP_PARTICULAR_STAT,
+      &p->mirrorquant,
+      GAL_DATA_TYPE_FLOAT32,
+      GAL_OPTIONS_RANGE_GT_0_LT_1,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
     },
 
 
 
+
+
     {
       0, 0, 0, 0,
-      "Sigma clipping:",
-      5
+      "Histogram or Cumulative frequency plot settings",
+      ARGS_GROUP_HIST_CFP
     },
     {
-      "nosigclip",
-      508,
-      0,
+      "numbins",
+      ARGS_OPTION_KEY_NUMBINS,
+      "INT",
       0,
-      "Do not preform sigma clipping.",
-      5
+      "Number of bins in histogram or CFP.",
+      ARGS_GROUP_HIST_CFP,
+      &p->numbins,
+      GAL_DATA_TYPE_SIZE_T,
+      GAL_OPTIONS_RANGE_GT_0,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
     },
     {
-      "sigclipmultip",
-      'u',
-      "FLT",
+      "lowerbin",
+      ARGS_OPTION_KEY_LOWERBIN,
       0,
-      "Multiple of standard deviation in sigma-clipping.",
-      5
+      0,
+      "Save interval lower limit, not center",
+      ARGS_GROUP_HIST_CFP,
+      &p->lowerbin,
+      GAL_OPTIONS_NO_ARG_TYPE,
+      GAL_OPTIONS_RANGE_GT_0,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
     },
     {
-      "sigcliptolerance",
-      't',
-      "FLT",
+      "normalize",
+      ARGS_OPTION_KEY_NORMALIZE,
+      0,
       0,
-      "Difference in STD tolerance to halt iteration.",
-      5
+      "Set sum of all bins to 1.",
+      ARGS_GROUP_HIST_CFP,
+      &p->normalize,
+      GAL_OPTIONS_NO_ARG_TYPE,
+      GAL_OPTIONS_RANGE_GT_0,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
     },
     {
-      "sigclipnum",
-      'g',
-      "INT",
+      "onebinstart",
+      ARGS_OPTION_KEY_ONEBINSTART,
+      "FLT",
       0,
-      "Number of times to do sigma clipping.",
-      5
+      "Shift bins so one bin starts on this value.",
+      ARGS_GROUP_HIST_CFP,
+      &p->onebinstart,
+      GAL_DATA_TYPE_FLOAT32,
+      GAL_OPTIONS_RANGE_ANY,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
     },
-
-
     {
-      0, 0, 0, 0,
-      "Operating modes:",
-      -1
+      "maxbinone",
+      ARGS_OPTION_KEY_MAXBINONE,
+      0,
+      0,
+      "Scale such that the maximum bin has value of one.",
+      ARGS_GROUP_HIST_CFP,
+      &p->maxbinone,
+      GAL_OPTIONS_NO_ARG_TYPE,
+      GAL_OPTIONS_RANGE_GT_0,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
     },
 
 
@@ -379,219 +364,21 @@ static struct argp_option options[] =
 
 
 
-/* Parse a single option: */
-static error_t
-parse_opt(int key, char *arg, struct argp_state *state)
-{
-  /* Save the arguments structure: */
-  struct statisticsparams *p = state->input;
-
-  /* Set the pointer to the common parameters for all programs
-     here: */
-  state->child_inputs[0]=&p->cp;
-
-  /* In case the user incorrectly uses the equal sign (for example
-     with a short format or with space in the long format, then `arg`
-     start with (if the short version was called) or be (if the long
-     version was called with a space) the equal sign. So, here we
-     check if the first character of arg is the equal sign, then the
-     user is warned and the program is stopped: */
-  if(arg && arg[0]=='=')
-    argp_error(state, "incorrect use of the equal sign (`=`). For short "
-               "options, `=` should not be used and for long options, "
-               "there should be no space between the option, equal sign "
-               "and value");
-
-  switch(key)
-    {
+/* Define the child argp structure. */
+struct argp
+gal_options_common_child = {gal_commonopts_options,
+                            gal_options_common_argp_parse,
+                            NULL, NULL, NULL, NULL, NULL};
 
-    /* Input: */
-    case 'M':
-      gal_checkset_allocate_copy_set(arg, &p->up.maskname, &p->up.masknameset);
-      break;
-    case 'H':
-      gal_checkset_allocate_copy_set(arg, &p->up.mhdu, &p->up.mhduset);
-      break;
-    case 'r':
-      p->ignoremin=1;
-      break;
-    case 'd':
-      gal_checkset_float_l_0(arg, &p->mirrordist, "mirrordist", key, SPACK,
-                             NULL, 0);
-      p->up.mirrordistset=1;
-      break;
-
-    /* Output: */
-    case 'l':
-      p->lowerbin=1;
-      break;
-    case 'B':
-      gal_checkset_any_float(arg, &p->onebinvalue, "onebinvalue", key, SPACK,
-                             NULL, 0);
-      p->up.onebinvalueset=1;
-      break;
-    case 'A':
-      p->asciihist=0;
-      break;
-    case 509:
-      p->mhistname="a";
-      break;
-    case 510:
-      gal_checkset_float_l_0_s_1(arg, &p->mirror, "mirrorquant", key, SPACK,
-                                 NULL, 0);
-      break;
-    case 511:
-      p->histrangeformirror=1;
-      break;
-    case 503:
-      gal_checkset_float_l_0(arg, &p->mirrorplotdist, "mirrorplotdist", key,
-                             SPACK, NULL, 0);
-      p->up.mirrorplotdistset=1;
-      break;
-
-    /* Histogram */
-    case 500:
-      p->histname=NULL;
-      break;
-    case 501:
-      p->normhist=1;
-      break;
-    case 502:
-      p->maxhistone=1;
-      break;
-    case 'n':
-      gal_checkset_sizet_l_zero(arg, &p->histnumbins, "histnumbins", key, 
SPACK,
-                                NULL, 0);
-      p->up.histnumbinsset=1;
-      break;
-    case 'i':
-      gal_checkset_any_float(arg, &p->histmin, "histmin", key, SPACK, NULL, 0);
-      p->up.histminset=1;
-      break;
-    case 'x':
-      gal_checkset_any_float(arg, &p->histmax, "histmax", key, SPACK, NULL, 0);
-      p->up.histmaxset=1;
-      break;
-    case 'Q':
-      gal_checkset_float_l_0_s_1(arg, &p->histquant, "histquant", key, SPACK,
-                                 NULL, 0);
-      p->up.histquantset=1;
-      break;
-
-    /* Cumulative frequency plot: */
-    case 504:
-      p->cfpname=NULL;
-      break;
-    case 505:
-      p->normcfp=1;
-      break;
-    case 506:
-      p->maxcfpeqmaxhist=1;
-      break;
-    case 507:
-      p->cfpsimhist=1;
-      break;
-    case 'p':
-      gal_checkset_sizet_l_zero(arg, &p->cfpnum, "cfpnum", key, SPACK, NULL, 
0);
-      p->up.cfpnumset=1;
-      break;
-    case 'a':
-      gal_checkset_any_float(arg, &p->cfpmin, "cfpmin", key, SPACK, NULL, 0);
-      p->up.cfpminset=1;
-      break;
-    case 'b':
-      gal_checkset_any_float(arg, &p->cfpmax, "cfpmax", key, SPACK, NULL, 0);
-      p->up.cfpmaxset=1;
-      break;
-    case 'U':
-      gal_checkset_float_l_0_s_1(arg, &p->cfpquant, "cfpquant", key, SPACK,
-                                 NULL, 0);
-      p->up.cfpquantset=1;
-      break;
-
-
-    /* Sigma clipping: */
-    case 508:
-      p->sigclip=0;
-      break;
-    case 'u':
-      gal_checkset_float_l_0(arg, &p->sigclipmultip, "sigclipmultip", key,
-                             SPACK, NULL, 0);
-      p->up.sigclipmultipset=1;
-      break;
-    case 't':
-      gal_checkset_float_l_0(arg, &p->sigcliptolerance, "sigcliptolerance", 
key,
-                             SPACK, NULL, 0);
-      p->up.sigcliptoleranceset=1;
-      break;
-    case 'g':
-      gal_checkset_sizet_l_zero(arg, &p->sigclipnum, "sigclipnum", key, SPACK,
-                                NULL, 0);
-      p->up.sigclipnumset=1;
-      break;
-
-    /* Operating modes: */
-
-
-
-    /* Read the non-option arguments: */
-    case ARGP_KEY_ARG:
-
-      /* See what type of input value it is and put it in. */
-      if( gal_fits_name_is_fits(arg) )
-        {
-          if(p->up.inputname)
-            argp_error(state, "only one input image should be given");
-          else
-            p->up.inputname=arg;
-        }
-      else
-        argp_error(state, "%s is not a valid file type", arg);
-      break;
-
-
-
-
-
-    /* The command line options and arguments are finished. */
-    case ARGP_KEY_END:
-      if(p->cp.setdirconf==0 && p->cp.setusrconf==0
-         && p->cp.printparams==0)
-        {
-          if(state->arg_num==0)
-            argp_error(state, "no argument given");
-          if(p->up.inputname==NULL)
-            argp_error(state, "no input FITS image(s) provided");
-        }
-      break;
-
-
-
-
-
-    default:
-      return ARGP_ERR_UNKNOWN;
-    }
-  return 0;
-}
-
-
-
-
-
-/* Specify the children parsers: */
-struct argp_child children[]=
-  {
-    {&commonargp, 0, NULL, 0},
-    {0, 0, 0, 0}
-  };
-
-
-
-
-
-/* Basic structure defining the whole argument reading process. */
-static struct argp thisargp = {options, parse_opt, args_doc,
-                               doc, children, NULL, NULL};
+/* Use the child argp structure in list of children (only one for now). */
+struct argp_child
+children[]=
+{
+  {&gal_options_common_child, 0, NULL, 0},
+  {0, 0, 0, 0}
+};
 
+/* Set all the necessary argp parameters. */
+struct argp
+thisargp = {program_options, parse_opt, args_doc, doc, children, NULL, NULL};
 #endif
diff --git a/bin/statistics/aststatistics.conf 
b/bin/statistics/aststatistics.conf
index 45e9b0f..dd22c54 100644
--- a/bin/statistics/aststatistics.conf
+++ b/bin/statistics/aststatistics.conf
@@ -19,17 +19,8 @@
 
 # Input image:
  hdu                 0
- mirrordist          1.5
 
-# Output:
- mirrorplotdist      2
- onebinvalue         nan
-
-# Histogram:
- histnumbins         200
-
-# Cumulative frequency plot:
- cfpnum              200
- sigclipmultip       4.000
- sigcliptolerance    0.200
- sigclipnum          5
\ No newline at end of file
+# Common options
+ minmapsize          1000000000
+ tableformat         fits-binary
+ searchin            name
\ No newline at end of file
diff --git a/bin/statistics/authors-cite.h b/bin/statistics/authors-cite.h
index e90da98..185a930 100644
--- a/bin/statistics/authors-cite.h
+++ b/bin/statistics/authors-cite.h
@@ -1,5 +1,5 @@
 /*********************************************************************
-Statistics - Get general statistics about the image.
+Statistics - Statistical analysis on input dataset.
 Statistics is part of GNU Astronomy Utilities (Gnuastro) package.
 
 Original author:
diff --git a/bin/statistics/main.c b/bin/statistics/main.c
index 61cffec..93ce912 100644
--- a/bin/statistics/main.c
+++ b/bin/statistics/main.c
@@ -1,5 +1,5 @@
 /*********************************************************************
-Statistics - Get general statistics about the image.
+Statistics - Statistical analysis on input dataset.
 Statistics is part of GNU Astronomy Utilities (Gnuastro) package.
 
 Original author:
@@ -36,20 +36,20 @@ int
 main (int argc, char *argv[])
 {
   struct timeval t1;
-  struct statisticsparams p={{0}, {0}, 0};
+  struct statisticsparams p={{0}, 0};
 
   /* Set the starting time. */
   time(&p.rawtime);
   gettimeofday(&t1, NULL);
 
   /* Read the input parameters. */
-  setparams(argc, argv, &p);
+  ui_read_check_inputs_setup(argc, argv, &p);
 
   /* Run MakeProfiles */
   statistics(&p);
 
   /* Free all non-freed allocations. */
-  freeandreport(&p, &t1);
+  ui_free_report(&p);
 
   /* Return successfully.*/
   return 0;
diff --git a/bin/statistics/main.h b/bin/statistics/main.h
index 6c7fba1..2f9a796 100644
--- a/bin/statistics/main.h
+++ b/bin/statistics/main.h
@@ -1,5 +1,5 @@
 /*********************************************************************
-Statistics - Get general statistics about the image.
+Statistics - Statistical analysis on input dataset.
 Statistics is part of GNU Astronomy Utilities (Gnuastro) package.
 
 Original author:
@@ -23,111 +23,50 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #ifndef MAIN_H
 #define MAIN_H
 
-#include <commonparams.h>
+/* Include necessary headers */
+#include <gnuastro/data.h>
 
-/* Progarm name macros: */
-#define SPACK           "aststatistics"    /* Subpackage executable name. */
-#define SPACK_NAME      "Statistics"       /* Subpackage full name.       */
-#define SPACK_STRING    SPACK_NAME" ("PACKAGE_NAME") "PACKAGE_VERSION
+#include <options.h>
 
+/* Progarm names.  */
+#define PROGRAM_NAME   "Statistics"    /* Program full name.       */
+#define PROGRAM_EXEC   "aststatistics" /* Program executable name. */
+#define PROGRAM_STRING PROGRAM_NAME" (" PACKAGE_NAME ") " PACKAGE_VERSION
 
 
-#define PRINTFLT  "%.6f\n"
-#define PRINTINT  "%.0f\n"
-#define STRVAL    "   -- %-45s%s\n"
-#define FNAMEVAL  "   -- %-45s%f\n"
-#define SNAMEVAL  "   -- %-45s%zu\n"
-#define ASCIIHISTNUMBINS    60
-#define ASCIIHISTHEIGHT     10
-#define HISTSTRING     "Histogram"
-#define CFPSTRING      "Cumulative Frequency Plot"
-
-
-
-struct uiparams
-{
-  char             *inputname;  /* Name of input file.               */
-  char              *maskname;  /* Name of mask image.               */
-  char                  *mhdu;  /* Mask image HDU.                   */
-  int             masknameset;
-  int       masknameallocated;
-  int                 mhduset;
-  int           mirrordistset;
-  int          onebinvalueset;
-  int          histnumbinsset;
-  int              histminset;
-  int              histmaxset;
-  int            histquantset;
-  int               cfpnumset;
-  int               cfpminset;
-  int               cfpmaxset;
-  int             cfpquantset;
-  int       mirrorplotdistset;
-  int        sigclipmultipset;
-  int     sigcliptoleranceset;
-  int           sigclipnumset;
-
-};
 
 
 
 
 
+/* Main program parameters structure */
 struct statisticsparams
 {
-  /* Other structures: */
-  struct uiparams         up; /* User interface parameters.             */
-  struct gal_commonparams cp; /* Common parameters.                     */
-
-  /* Input: */
-  float               *img;  /* Input image array.                      */
-  float            *sorted;  /* Sorted input data.                      */
-  size_t              size;  /* Number of non-blank data elements.      */
-  int            ignoremin;  /* Ignore all data with minimum value.     */
-  float         mirrordist;  /* Distance to go out for mode check.      */
-
-  /* Output: */
-  int            asciihist;  /* ==1: print an ASCII histogram.          */
-  float        onebinvalue;  /* Shift bins so one bin starts with this. */
-  int             lowerbin;  /* Interval lower limit as column 1.       */
-  float             mirror;  /* Mirror quantile.                        */
-  char         *mirrorhist;  /* Name of mirror histogram.               */
-  char          *mirrorcfp;  /* Name of mirror CFP.                     */
-  char          *mhistname;  /* Name of mode mirror histogram.          */
-  float     mirrorplotdist;  /* Dist. after mode to display on plot.    */
-  char           *mcfpname;  /* Name of mode mirror CFP.                */
-  int   histrangeformirror;  /* ==1: Use histogram range for mirror.    */
-
-  /* Histogram: */
-  char           *histname;  /* Histogram file name.                    */
-  int             normhist;  /* ==1: Normalize the histogram.           */
-  int           maxhistone;  /* Scale such that max bin is one.         */
-  size_t       histnumbins;  /* Number of bins in the histogram.        */
-  float            histmin;  /* Minimum value to use in histogram.      */
-  float            histmax;  /* Maximum value to use in histogram.      */
-  float          histquant;  /* Quantile range of histogram.            */
-
-  /* Cumulative frequency plot: */
-  char            *cfpname;  /* Cumultiave frequency plot file name.    */
-  int              normcfp;  /* ==1: Normalize the CFP.                 */
-  int      maxcfpeqmaxhist;  /* ==1: CFP max equal to historgram max.   */
-  int           cfpsimhist;  /* ==1: CFP range equal to hist range.     */
-  size_t            cfpnum;  /* The number of points to sample CFP.     */
-  float             cfpmin;  /* Minimum value for CFP.                  */
-  float             cfpmax;  /* Maximum value for CFP.                  */
-  float           cfpquant;  /* Quantile range of CFP.                  */
-
-  /* Sigma clipping: */
-  int              sigclip;  /* ==1: Do sigma clipping.                 */
-  float      sigclipmultip;  /* Multiple of STD in sigma clipping.      */
-  float   sigcliptolerance;  /* Tolerance level in sigma clipping.      */
-  size_t        sigclipnum;  /* Number of times to sigma clip.          */
-
-
-  /* Operating mode: */
-
-  /* Internal: */
-  time_t           rawtime;  /* Starting time of the program.           */
+  /* From command-line */
+  struct gal_options_common_params cp; /* Common parameters.             */
+  struct gal_linkedlist_ill  *toprint; /* Values to print in one row.    */
+  char          *inputname;  /* Input filename.                          */
+  char             *column;  /* Column name or number if input is table. */
+  float       greaterequal;  /* Only use values >= this value.           */
+  float           lessthan;  /* Only use values <  this value.           */
+  float      quantilerange;  /* Quantile (Q) range: from Q to 1-Q.       */
+
+  uint8_t        asciihist;  /* Print an ASCII histogram.                */
+  uint8_t        histogram;  /* Save histogram in output.                */
+  uint8_t       cumulative;  /* Save cumulative distibution in output.   */
+  char         *sigclipstr;  /* Multiple of sigma, and tolerance or num. */
+  float        mirrorquant;  /* Quantile of mirror distribution to save. */
+  size_t           numbins;  /* Number of bins in histogram or CFP.      */
+  uint8_t         lowerbin;  /* Save interval lower limit, not center.   */
+  uint8_t        normalize;  /* set the sum of all bins to 1.            */
+  float        onebinstart;  /* Shift bins to start at this value.       */
+  uint8_t        maxbinone;  /* Set the maximum bin to 1.                */
+
+  /* Internal */
+  gal_data_t        *input;  /* Input data structure.                    */
+  int               isfits;  /* Input is a FITS file.                    */
+  int             hdu_type;  /* Type of HDU (image or table).            */
+  time_t           rawtime;  /* Starting time of the program.            */
 };
 
 #endif
diff --git a/bin/statistics/statistics.c b/bin/statistics/statistics.c
index d05b9f1..83dad99 100644
--- a/bin/statistics/statistics.c
+++ b/bin/statistics/statistics.c
@@ -1,5 +1,5 @@
 /*********************************************************************
-Statistics - Get general statistics about the image.
+Statistics - Statistical analysis on input dataset.
 Statistics is part of GNU Astronomy Utilities (Gnuastro) package.
 
 Original author:
@@ -30,6 +30,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <string.h>
 #include <stdlib.h>
 
+#include <gnuastro/arithmetic.h>
 #include <gnuastro/statistics.h>
 
 #include <timing.h>
@@ -44,6 +45,11 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 void
 reportsimplestats(struct statisticsparams *p)
 {
+  printf("\n... in reportsimplestats ...\n");
+  exit(1);
+
+#if 0
+
   double sum;
   size_t modeindex;
   float modequant, symvalue;
@@ -84,6 +90,7 @@ reportsimplestats(struct statisticsparams *p)
                                      p->mcfpname, (p->histrangeformirror
                                                    ? 0.0f
                                                    : p->mirrorplotdist) );
+#endif
 }
 
 
@@ -91,31 +98,32 @@ reportsimplestats(struct statisticsparams *p)
 
 
 void
-printasciihist(struct statisticsparams *p)
+ascii_hist(struct statisticsparams *p)
 {
   size_t j;
-  float *bins;
-  float quant=-1.0f; /* histmin and histmax were already set before. */
+  size_t numbins=60, histheight=10;
+  float *bins, *sorted=p->input->array;
+  float quant=0.0f; /* histmin and histmax were already set before. */
   int i, binonzero=0, normhist=0, maxhistone=1;
 
   /* Find the histogram for the ASCII plot: */
-  gal_statistics_set_bins(p->sorted, p->size, ASCIIHISTNUMBINS, p->histmin,
-                          p->histmax, binonzero, quant, &bins);
-  gal_statistics_histogram(p->sorted, p->size, bins, ASCIIHISTNUMBINS,
+  gal_statistics_set_bins(sorted, p->input->size, numbins, 0, 0,
+                          binonzero, quant, &bins);
+  gal_statistics_histogram(sorted, p->input->size, bins, numbins,
                            normhist, maxhistone);
 
   /* It's maximum value is set to one. Multiply that by the desired
      height in pixels. */
-  for(j=0;j<ASCIIHISTNUMBINS;++j)
-    bins[j*2+1]*=ASCIIHISTHEIGHT;
+  for(j=0;j<numbins;++j)
+    bins[j*2+1]*=histheight;
 
   /* Plot the ASCII histogram: */
-  printf("   -- ASCII histogram in the range: %f  --  %f:\n",
-         p->histmin, p->histmax);
-  for(i=ASCIIHISTHEIGHT;i>=0;--i)
+  printf("\nRange of histogram: %g to %g\n\n", sorted[0],
+         sorted[p->input->size-1]);
+  for(i=histheight;i>=0;--i)
     {
       printf("    |");
-      for(j=0;j<ASCIIHISTNUMBINS;++j)
+      for(j=0;j<numbins;++j)
         {
           if(bins[j*2+1]>=((float)i-0.5f)
              && bins[j*2+1]>0.0f) printf("*");
@@ -124,7 +132,7 @@ printasciihist(struct statisticsparams *p)
       printf("\n");
     }
   printf("    |");
-  for(j=0;j<ASCIIHISTNUMBINS;++j) printf("-");
+  for(j=0;j<numbins;++j) printf("-");
   printf("\n\n");
 
   /* Clean up.*/
@@ -138,6 +146,10 @@ void
 printhistcfp(struct statisticsparams *p, float *bins, size_t numbins,
              char *filename, char *outputtype)
 {
+  printf("\n... in printhistcfp ...\n");
+  exit(1);
+
+#if 0
   float d;
   size_t i;
   FILE *out;
@@ -201,6 +213,7 @@ printhistcfp(struct statisticsparams *p, float *bins, 
size_t numbins,
     for(i=0;i<numbins;++i)
       fprintf(out, "%-20.6f"PRINTINT, bins[i*2]+d, bins[i*2+1]);
   fclose(out);
+#endif
 }
 
 
@@ -208,8 +221,13 @@ printhistcfp(struct statisticsparams *p, float *bins, 
size_t numbins,
 
 
 void
-imgstat(struct statisticsparams *p)
+statistics(struct statisticsparams *p)
 {
+
+  /* Print the ASCII histogram. */
+  if(p->asciihist) ascii_hist(p);
+
+#if 0
   int r;
   size_t i;
   float quant=-1.0f;                   /* The quantile was already   */
@@ -270,8 +288,8 @@ imgstat(struct statisticsparams *p)
   /* Make the mirror distribution if asked for: */
   if(isnan(p->mirror)==0)
     gal_statistics_mode_mirror_plots(p->sorted, p->size,
-                                     
gal_statistics_index_from_quantile(p->size,
-                                                                        
p->mirror),
+                 gal_statistics_index_from_quantile(p->size,
+                                                    p->mirror),
                                      p->histmin, p->histmax, p->histnumbins,
                                      p->mirrorhist, p->mirrorcfp,
                                      (p->histrangeformirror ? 0.0f
@@ -298,4 +316,5 @@ imgstat(struct statisticsparams *p)
 
   /* Free the allocated arrays: */
   if(p->histname || p->cfpname) free(bins);
+#endif
 }
diff --git a/bin/statistics/statistics.h b/bin/statistics/statistics.h
index 3bddbb0..bc7b59a 100644
--- a/bin/statistics/statistics.h
+++ b/bin/statistics/statistics.h
@@ -1,5 +1,5 @@
 /*********************************************************************
-Statistics - Get general statistics about the image.
+Statistics - Statistical analysis on input dataset.
 Statistics is part of GNU Astronomy Utilities (Gnuastro) package.
 
 Original author:
diff --git a/bin/statistics/ui.c b/bin/statistics/ui.c
index 70f83c0..735bf78 100644
--- a/bin/statistics/ui.c
+++ b/bin/statistics/ui.c
@@ -1,5 +1,5 @@
 /*********************************************************************
-Statistics - Get general statistics about the image.
+Statistics - Statistical analysis on input dataset.
 Statistics is part of GNU Astronomy Utilities (Gnuastro) package.
 
 Original author:
@@ -22,333 +22,238 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 **********************************************************************/
 #include <config.h>
 
-#include <math.h>
-#include <stdio.h>
+#include <argp.h>
 #include <errno.h>
 #include <error.h>
-#include <stdlib.h>
-#include <string.h>
-#include <fitsio.h>
+#include <stdio.h>
 
 #include <gnuastro/fits.h>
 #include <gnuastro/qsort.h>
-#include <gnuastro/array.h>
-#include <gnuastro/txtarray.h>
-#include <gnuastro/statistics.h>
+#include <gnuastro/table.h>
+#include <gnuastro/table.h>
+#include <gnuastro/arithmetic.h>
+#include <gnuastro/linkedlist.h>
 
-#include <nproc.h>               /* From Gnulib.                     */
-#include <timing.h>              /* Includes time.h and sys/time.h   */
+#include <timing.h>
+#include <options.h>
 #include <checkset.h>
-#include <commonargs.h>
-#include <configfiles.h>
 #include <fixedstringmacros.h>
 
 #include "main.h"
 
-#include "ui.h"                  /* Needs main.h                   */
-#include "args.h"                /* Needs main.h, includes argp.h. */
+#include "ui.h"
+#include "authors-cite.h"
 
 
-/* Set the file names of the places where the default parameters are
-   put. */
-#define CONFIG_FILE SPACK CONF_POSTFIX
-#define SYSCONFIG_FILE SYSCONFIG_DIR "/" CONFIG_FILE
-#define USERCONFIG_FILEEND USERCONFIG_DIR CONFIG_FILE
-#define CURDIRCONFIG_FILE CURDIRCONFIG_DIR CONFIG_FILE
 
 
 
+/**************************************************************/
+/*********      Argp necessary global entities     ************/
+/**************************************************************/
+/* Definition parameters for the Argp: */
+const char *
+argp_program_version = PROGRAM_STRING "\n"
+                       GAL_STRINGS_COPYRIGHT
+                       "\n\nWritten/developed by "PROGRAM_AUTHORS;
 
+const char *
+argp_program_bug_address = PACKAGE_BUGREPORT;
 
+static char
+args_doc[] = "ASTRdata";
 
+const char
+doc[] = GAL_STRINGS_TOP_HELP_INFO PROGRAM_NAME" will print the basic "
+  "statistics of the input image pixel flux distribution. All blank pixels "
+  "or pixels specified by a mask image will be ignored.\n"
+  GAL_STRINGS_MORE_HELP_INFO
+  /* After the list of options: */
+  "\v"
+  PACKAGE_NAME" home page: "PACKAGE_URL;
 
 
 
 
-/**************************************************************/
-/**************       Options and parameters    ***************/
-/**************************************************************/
-void
-readconfig(char *filename, struct statisticsparams *p)
+
+/* Option groups particular to this program. */
+enum program_args_groups
 {
-  FILE *fp;
-  size_t lineno=0, len=200;
-  char *line, *name, *value;
-  struct uiparams *up=&p->up;
-  struct gal_commonparams *cp=&p->cp;
-  char key='a';        /* Not used, just a place holder. */
-
-  /* When the file doesn't exist or can't be opened, it is ignored. It
-     might be intentional, so there is no error. If a parameter is
-     missing, it will be reported after all defaults are read. */
-  fp=fopen(filename, "r");
-  if (fp==NULL) return;
-
-
-  /* Allocate some space for `line` with `len` elements so it can
-     easily be freed later on. The value of `len` is arbitarary at
-     this point, during the run, getline will change it along with the
-     pointer to line. */
-  errno=0;
-  line=malloc(len*sizeof *line);
-  if(line==NULL)
-    error(EXIT_FAILURE, errno, "ui.c: %zu bytes in readdefaults",
-          len * sizeof *line);
+  ARGS_GROUP_IN_ONE_ROW = GAL_OPTIONS_GROUP_AFTER_COMMON,
+  ARGS_GROUP_PARTICULAR_STAT,
+  ARGS_GROUP_HIST_CFP,
+};
 
-  /* Read the tokens in the file:  */
-  while(getline(&line, &len, fp) != -1)
-    {
-      /* Prepare the "name" and "value" strings, also set lineno. */
-      GAL_CONFIGFILES_START_READING_LINE;
 
 
 
+/* Available letters for short options:
+
+   a b c d e f i j k p r u v w x y z
+   B E F G J L R W X Y Z
+*/
+enum option_keys_enum
+{
+  /* With short-option version. */
+  ARGS_OPTION_KEY_COLUMN       = 'c',
+  ARGS_OPTION_KEY_GREATEREQUAL = 'g',
+  ARGS_OPTION_KEY_LESSTHAN     = 'l',
+  ARGS_OPTION_KEY_QUANTRANGE   = 'Q',
+  ARGS_OPTION_KEY_MEAN         = 'm',
+  ARGS_OPTION_KEY_STD          = 't',
+  ARGS_OPTION_KEY_MEDIAN       = 'M',
+  ARGS_OPTION_KEY_MODE         = 'O',
+  ARGS_OPTION_KEY_ASCIIHIST    = 'A',
+  ARGS_OPTION_KEY_HISTOGRAM    = 'H',
+  ARGS_OPTION_KEY_CUMULATIVE   = 'C',
+  ARGS_OPTION_KEY_SIGMACLIP    = 's',
+  ARGS_OPTION_KEY_NORMALIZE    = 'n',
+
+  /* Only with long version (start with a value 1000, the rest will be set
+     automatically). */
+  ARGS_OPTION_KEY_NUMBER       = 1000,
+  ARGS_OPTION_KEY_MINIMUM,
+  ARGS_OPTION_KEY_MAXIMUM,
+  ARGS_OPTION_KEY_SUM,
+  ARGS_OPTION_KEY_MIRROR,
+  ARGS_OPTION_KEY_NUMBINS,
+  ARGS_OPTION_KEY_LOWERBIN,
+  ARGS_OPTION_KEY_ONEBINSTART,
+  ARGS_OPTION_KEY_MAXBINONE,
+};
+
+
 
-      /* Inputs: */
-      if(strcmp(name, "hdu")==0)
-        gal_checkset_allocate_copy_set(value, &cp->hdu, &cp->hduset);
 
-      else if (strcmp(name, "mask")==0)
-        gal_checkset_allocate_copy_set(value, &up->maskname,
-                                       &up->masknameset);
 
-      else if (strcmp(name, "mhdu")==0)
-        gal_checkset_allocate_copy_set(value, &up->mhdu, &up->mhduset);
-      else if(strcmp(name, "mirrordist")==0)
-        {
-          if(up->mirrordistset) continue;
-          gal_checkset_float_l_0(value, &p->mirrordist, name, key, SPACK,
-                                 filename, lineno);
-          up->mirrordistset=1;
-        }
 
 
 
-      /* Outputs */
-      else if(strcmp(name, "output")==0)
-        gal_checkset_allocate_copy_set(value, &cp->output, &cp->outputset);
 
-      else if(strcmp(name, "mirrorplotdist")==0)
-        {
-          if(up->mirrorplotdistset) continue;
-          gal_checkset_float_l_0(value, &p->mirrorplotdist, name, key,
-                                 SPACK, filename, lineno);
-          up->mirrorplotdistset=1;
-        }
-      else if(strcmp(name, "onebinvalue")==0)
-        {
-          if(up->onebinvalueset) continue;
-          gal_checkset_any_float(value, &p->onebinvalue, name, key, SPACK,
-                                 filename, lineno);
-          up->onebinvalueset=1;
-        }
 
 
-      /* Histogram: */
-      else if(strcmp(name, "histnumbins")==0)
-        {
-          if(up->histnumbinsset) continue;
-          gal_checkset_sizet_l_zero(value, &p->histnumbins, name, key,
-                                    SPACK, filename, lineno);
-          up->histnumbinsset=1;
-        }
-      else if(strcmp(name, "histmin")==0)
-        {
-          if(up->histminset) continue;
-          gal_checkset_any_float(value, &p->histmin, name, key, SPACK,
-                                 filename, lineno);
-          up->histminset=1;
-        }
-      else if(strcmp(name, "histmax")==0)
-        {
-          if(up->histmaxset) continue;
-          gal_checkset_any_float(value, &p->histmax, name, key, SPACK,
-                                 filename, lineno);
-          up->histmaxset=1;
-        }
-      else if(strcmp(name, "histquant")==0)
-        {
-          if(up->histquantset) continue;
-          gal_checkset_float_l_0_s_1(value, &p->histquant, name, key,
-                                     SPACK, filename, lineno);
-          up->histquantset=1;
-        }
 
 
-      /* Cumulative Frequency Plot: */
-      else if(strcmp(name, "cfpnum")==0)
-        {
-          if(up->cfpnumset) continue;
-          gal_checkset_sizet_l_zero(value, &p->cfpnum, name, key, SPACK,
-                                    filename, lineno);
-          up->cfpnumset=1;
-        }
-      else if(strcmp(name, "cfpmin")==0)
-        {
-          if(up->cfpminset) continue;
-          gal_checkset_any_float(value, &p->cfpmin, name, key, SPACK,
-                                 filename, lineno);
-          up->cfpminset=1;
-        }
-      else if(strcmp(name, "cfpmax")==0)
-        {
-          if(up->cfpmaxset) continue;
-          gal_checkset_any_float(value, &p->cfpmax, name, key, SPACK,
-                                 filename, lineno);
-          up->cfpmaxset=1;
-        }
-      else if(strcmp(name, "cfpquant")==0)
-        {
-          if(up->cfpquantset) continue;
-          gal_checkset_float_l_0_s_1(value, &p->cfpquant, name, key,
-                                     SPACK, filename, lineno);
-          up->cfpquantset=1;
-        }
 
-      /* Sigma clipping: */
-      else if(strcmp(name, "sigclipmultip")==0)
-        {
-          if(up->sigclipmultipset) continue;
-          gal_checkset_float_l_0(value, &p->sigclipmultip, name, key,
-                                 SPACK, filename, lineno);
-          up->sigclipmultipset=1;
-        }
-      else if(strcmp(name, "sigcliptolerance")==0)
-        {
-          if(up->sigcliptoleranceset) continue;
-          gal_checkset_float_l_0(value, &p->sigcliptolerance, name, key,
-                                 SPACK, filename, lineno);
-          up->sigcliptoleranceset=1;
-        }
-      else if(strcmp(name, "sigclipnum")==0)
-        {
-          if(up->sigclipnumset) continue;
-          gal_checkset_sizet_l_zero(value, &p->sigclipnum, name, key,
-                                    SPACK, filename, lineno);
-          up->sigclipnumset=1;
-        }
 
 
-      /* Operating modes: */
-      /* Read options common to all programs */
-      GAL_CONFIGFILES_READ_COMMONOPTIONS_FROM_CONF
 
 
-      else
-        error_at_line(EXIT_FAILURE, 0, filename, lineno,
-                      "`%s` not recognized.\n", name);
-    }
 
-  free(line);
-  fclose(fp);
+/**************************************************************/
+/*********    Initialize & Parse command-line    **************/
+/**************************************************************/
+static void
+ui_initialize_options(struct statisticsparams *p,
+                      struct argp_option *program_options,
+                      struct argp_option *gal_commonopts_options)
+{
+  size_t i;
+  struct gal_options_common_params *cp=&p->cp;
+
+  /* Set the necessary common parameters structure. */
+  cp->program_struct     = p;
+  cp->poptions           = program_options;
+  cp->program_name       = PROGRAM_NAME;
+  cp->program_exec       = PROGRAM_EXEC;
+  cp->program_bibtex     = PROGRAM_BIBTEX;
+  cp->program_authors    = PROGRAM_AUTHORS;
+  cp->coptions           = gal_commonopts_options;
+
+  /* Program-specific initializers */
+  p->lessthan            = NAN;
+  p->onebinstart         = NAN;
+  p->greaterequal        = NAN;
+  p->quantilerange       = NAN;
+
+  /* Set the mandatory common options. */
+  for(i=0; !gal_options_is_last(&cp->coptions[i]); ++i)
+    switch(cp->coptions[i].key)
+      {
+      case GAL_OPTIONS_KEY_LOG:
+      case GAL_OPTIONS_KEY_TYPE:
+        cp->coptions[i].flags=OPTION_HIDDEN;
+        break;
+
+      case GAL_OPTIONS_KEY_SEARCHIN:
+      case GAL_OPTIONS_KEY_MINMAPSIZE:
+      case GAL_OPTIONS_KEY_TABLEFORMAT:
+        cp->coptions[i].mandatory=GAL_OPTIONS_MANDATORY;
+        break;
+      }
 }
 
 
 
 
 
-void
-printvalues(FILE *fp, struct statisticsparams *p)
+/* Parse a single option: */
+error_t
+parse_opt(int key, char *arg, struct argp_state *state)
 {
-  struct uiparams *up=&p->up;
-  struct gal_commonparams *cp=&p->cp;
-
-  /* Print all the options that are set. Separate each group with a
-     commented line explaining the options in that group. */
-  fprintf(fp, "\n# Input image:\n");
-  if(cp->hduset)
-    GAL_CHECKSET_PRINT_STRING_MAYBE_WITH_SPACE("hdu", cp->hdu);
-  if(up->masknameset)
-    GAL_CHECKSET_PRINT_STRING_MAYBE_WITH_SPACE("mask", up->maskname);
-  if(up->mhdu)
-    GAL_CHECKSET_PRINT_STRING_MAYBE_WITH_SPACE("mhdu", up->mhdu);
-  if(up->mirrordistset)
-    fprintf(fp, CONF_SHOWFMT"%.2f\n", "mirrordist", p->mirrordist);
-
-  /* Output: */
-  fprintf(fp, "\n# Output:\n");
-  if(cp->outputset)
-    fprintf(fp, CONF_SHOWFMT"%s\n", "output", cp->output);
-  if(up->mirrorplotdistset)
-    fprintf(fp, CONF_SHOWFMT"%.2f\n", "mirrorplotdist", p->mirrorplotdist);
-  if(up->onebinvalueset)
-    fprintf(fp, CONF_SHOWFMT"%.5f\n", "onebinvalue", p->onebinvalue);
-
-  /* Histogram: */
-  fprintf(fp, "\n# Histogram:\n");
-  if(up->histnumbinsset)
-    fprintf(fp, CONF_SHOWFMT"%zu\n", "histnumbins", p->histnumbins);
-  if(up->histminset)
-    fprintf(fp, CONF_SHOWFMT"%.5f\n", "histmin", p->histmin);
-  if(up->histmaxset)
-    fprintf(fp, CONF_SHOWFMT"%.5f\n", "histmax", p->histmax);
-  if(up->histquantset)
-    fprintf(fp, CONF_SHOWFMT"%.3f\n", "histquant", p->histquant);
-
-  /* Cumulative frequency plot: */
-  fprintf(fp, "\n# Cumulative frequency plot:\n");
-  if(up->cfpnumset)
-    fprintf(fp, CONF_SHOWFMT"%zu\n", "cfpnum", p->cfpnum);
-  if(up->cfpminset)
-    fprintf(fp, CONF_SHOWFMT"%.5f\n", "cfpmin", p->cfpmin);
-  if(up->cfpmaxset)
-    fprintf(fp, CONF_SHOWFMT"%.5f\n", "cfpmax", p->cfpmax);
-  if(up->cfpquantset)
-    fprintf(fp, CONF_SHOWFMT"%.3f\n", "cfpquant", p->cfpquant);
-
-  /* Sigma clipping: */
-  if(up->sigclipmultipset)
-    fprintf(fp, CONF_SHOWFMT"%.3f\n", "sigclipmultip", p->sigclipmultip);
-  if(up->sigcliptoleranceset)
-    fprintf(fp, CONF_SHOWFMT"%.3f\n", "sigcliptolerance",
-            p->sigcliptolerance);
-  if(up->sigclipnumset)
-    fprintf(fp, CONF_SHOWFMT"%zu\n", "sigclipnum", p->sigclipnum);
-
-
-  /* For the operating mode, first put the macro to print the common
-     options, then the (possible options particular to this
-     program). */
-  fprintf(fp, "\n# Operating mode:\n");
-  GAL_CONFIGFILES_PRINT_COMMONOPTIONS;
-}
+  struct statisticsparams *p = state->input;
+
+  /* Pass `gal_options_common_params' into the child parser.  */
+  state->child_inputs[0] = &p->cp;
+
+  /* In case the user incorrectly uses the equal sign (for example
+     with a short format or with space in the long format, then `arg`
+     start with (if the short version was called) or be (if the long
+     version was called with a space) the equal sign. So, here we
+     check if the first character of arg is the equal sign, then the
+     user is warned and the program is stopped: */
+  if(arg && arg[0]=='=')
+    argp_error(state, "incorrect use of the equal sign (`=`). For short "
+               "options, `=` should not be used and for long options, "
+               "there should be no space between the option, equal sign "
+               "and value");
+
+  /* Set the key to this option. */
+  switch(key)
+    {
+
+    /* Read the non-option tokens (arguments): */
+    case ARGP_KEY_ARG:
+      if(p->inputname)
+        argp_error(state, "only one argument (input file) should be given");
+      else
+        p->inputname=arg;
+      break;
+
 
+    /* This is an option, set its value. */
+    default:
+      return gal_options_set_from_key(key, arg, p->cp.poptions, &p->cp);
+    }
 
+  return 0;
+}
 
 
 
 
-/* Note that numthreads will be used automatically based on the
-   configure time. Note that those options which are not mandatory
-   must not be listed here. */
-void
-checkifset(struct statisticsparams *p)
+
+static void *
+ui_add_to_print_in_row(struct argp_option *option, char *arg,
+                       char *filename, size_t lineno, void *params)
 {
-  struct uiparams *up=&p->up;
-  struct gal_commonparams *cp=&p->cp;
-
-  int intro=0;
-  if(cp->hduset==0)
-    GAL_CONFIGFILES_REPORT_NOTSET("hdu");
-  if(up->mirrordistset==0)
-    GAL_CONFIGFILES_REPORT_NOTSET("mirrordist");
-  if(up->mirrorplotdistset==0)
-    GAL_CONFIGFILES_REPORT_NOTSET("mirrorplotdist");
-  if(up->onebinvalueset==0)
-    GAL_CONFIGFILES_REPORT_NOTSET("onebinvalue");
-  if(up->histnumbinsset==0)
-    GAL_CONFIGFILES_REPORT_NOTSET("histnumbins");
-  if(up->cfpnumset==0)
-    GAL_CONFIGFILES_REPORT_NOTSET("cfpnum");
-  if(up->sigclipmultipset==0)
-    GAL_CONFIGFILES_REPORT_NOTSET("sigclipmultip");
-  if(up->sigcliptoleranceset==0)
-    GAL_CONFIGFILES_REPORT_NOTSET("sigcliptolerance");
-  if(up->sigclipnumset==0)
-    GAL_CONFIGFILES_REPORT_NOTSET("sigclipnum");
-
-
-  GAL_CONFIGFILES_END_OF_NOTSET_REPORT;
+  struct statisticsparams *p=(struct statisticsparams *)params;
+
+  /* If this option is given in a configuration file, then `arg' will not
+     be NULL and we don't want to do anything if it is `0'. */
+  if( arg && arg[1]!='\0' && *arg!='0' && *arg!='1' )
+    error_at_line(EXIT_FAILURE, 0, filename, lineno, "the `--%s' "
+                  "option takes no arguments. In a configuration file "
+                  "it can only have the values `1' or `0', indicating "
+                  "if it should be used or not", option->name);
+
+  /* Only proceed if the (possibly given) argument is 1. */
+  if(arg && *arg=='0') return NULL;
+
+  /* Add this option to the print list. */
+  gal_linkedlist_add_to_ill(&p->toprint, option->key);
+
+  return NULL;
 }
 
 
@@ -373,80 +278,85 @@ checkifset(struct statisticsparams *p)
 /**************************************************************/
 /***************       Sanity Check         *******************/
 /**************************************************************/
-void
-sanitycheck(struct statisticsparams *p)
+/* Read and check ONLY the options. When arguments are involved, do the
+   check in `ui_check_options_and_arguments'. */
+static void
+ui_read_check_only_options(struct statisticsparams *p)
 {
-  char *basename;
+  /* Check if the format of the output table is valid, given the type of
+     the output. */
+  gal_table_check_fits_format(p->cp.output, p->cp.tableformat);
+
+  /* If less than and greater than are both given, make sure that the value
+     to greater than is smaller than the value to less-than. */
+  if( !isnan(p->lessthan) && !isnan(p->greaterequal)
+      && p->lessthan < p->greaterequal )
+    error(EXIT_FAILURE, 0, "the value to `--lessthan' (%g) must be larger "
+          "than the value to `--greaterequal' (%g)", p->lessthan,
+          p->greaterequal);
+
+  /* Less-than and greater-equal cannot be called together with
+     quantrange. */
+  if( ( !isnan(p->lessthan) || !isnan(p->greaterequal) )
+      && !isnan(p->quantilerange) )
+    error(EXIT_FAILURE, 0, "`--lessthan' and/or `--greaterequal' cannot "
+          "be called together with `--quantrange'");
+}
 
-  /* Make sure the input file exists. */
-  gal_checkset_check_file(p->up.inputname);
 
-  /* Set the p->up.maskname accordingly: */
-  gal_fits_file_or_ext_name(p->up.inputname, p->cp.hdu, p->up.masknameset,
-                                 &p->up.maskname, p->up.mhdu,
-                                 p->up.mhduset, "mask");
 
-  /* Set the names of the output files: */
-  if(p->cp.outputset) basename=p->cp.output;
-  else                basename=p->up.inputname;
-  if(p->histname)
-    {
-      p->histname=NULL;         /* It wasn't allocated. */
-      gal_checkset_automatic_output(basename, "_hist.txt",
-                                    p->cp.removedirinfo,
-                                    p->cp.dontdelete, &p->histname);
-    }
-  if(p->cfpname)
-    {
-      p->cfpname=NULL;         /* It wasn't allocated. */
-      gal_checkset_automatic_output(basename, "_cfp.txt",
-                                    p->cp.removedirinfo,
-                                    p->cp.dontdelete, &p->cfpname);
-    }
-  if(p->mhistname)             /* The mode mirror distribution will need */
-    {                          /* both a histogram and cfp.              */
-      p->mcfpname=p->mhistname=NULL;
-      gal_checkset_automatic_output(basename, "_modehist.txt",
-                                    p->cp.removedirinfo, p->cp.dontdelete,
-                                    &p->mhistname);
-      gal_checkset_automatic_output(basename, "_modecfp.txt",
-                                    p->cp.removedirinfo, p->cp.dontdelete,
-                                    &p->mcfpname);
-    }
-  if(isnan(p->mirror)==0)
-    {
-      p->mirrorhist=p->mirrorcfp=NULL;
-      gal_checkset_automatic_output(basename, "_mirrorhist.txt",
-                                    p->cp.removedirinfo, p->cp.dontdelete,
-                                    &p->mirrorhist);
-      gal_checkset_automatic_output(basename, "_mirrorcfp.txt",
-                                    p->cp.removedirinfo, p->cp.dontdelete,
-                                    &p->mirrorcfp);
-    }
 
 
-  /* If the cumulative frequency plot parameters are to depend on the
-     histogram, then make sure that the histogram will be created.*/
-  if(p->cfpname && p->histname==NULL)
+static void
+ui_check_options_and_arguments(struct statisticsparams *p)
+{
+  char *name=NULL;
+
+  if(p->inputname)
     {
-      if(p->cfpsimhist)
-        error(EXIT_FAILURE, 0, "without a histogram, `--cfpsimhist` is "
-              "meaningless");
-      if (p->maxcfpeqmaxhist)
-        error(EXIT_FAILURE, 0, "without a histogram, `--maxcfpeqmaxhist` "
-              "is meaningless");
-    }
+      /* If input is FITS. */
+      if( (p->isfits=gal_fits_name_is_fits(p->inputname)) )
+        {
+          /* Check if a HDU is given. */
+          if( p->cp.hdu==NULL )
+            error(EXIT_FAILURE, 0, "no HDU specified. When the input is a "
+                  "FITS file, a HDU must also be specified, you can use "
+                  "the `--hdu' (`-h') option and give it the HDU number "
+                  "(starting from zero), extension name, or anything "
+                  "acceptable by CFITSIO");
+
+          /* If its a table, make sure a column is also specified. */
+          p->hdu_type=gal_fits_hdu_type(p->inputname, p->cp.hdu);
+          if(p->hdu_type==IMAGE_HDU)
+            {
+              if(p->column)
+                error(EXIT_FAILURE, 0, "%s (hdu: %s): is a FITS image "
+                      "extension. The `--column' option is only applicable "
+                      "to tables.", p->inputname, p->cp.hdu);
+            }
+          else if(p->column==NULL)
+            asprintf(&name, "%s (hdu: %s)", p->inputname, p->cp.hdu);
+        }
 
-  /* Check if `--maxcfpeqmaxhist` and `--normcfp` are not called
-     together: */
-  if(p->normcfp && p->maxcfpeqmaxhist)
-    error(EXIT_FAILURE, 0, "`--normcfp` and `--maxcfpeqmaxhist` "
-          "cannot be called together\n");
+      /* If its not FITS, it must be a table. */
+      else
+        {
+          if(p->column==NULL) name=p->inputname;
+        }
 
-  /* Check if `normhist` and `maxhistone` are not called together: */
-  if(p->normhist && p->maxhistone)
-    error(EXIT_FAILURE, 0, "`--normhist` and `--histnumbins` cannot be "
-          "called together\n");
+      /* If a column was necessary, but not given, print an error. */
+      if(name)
+        error(EXIT_FAILURE, 0, "%s is a table but no column is "
+              "specified. Please use the `--column' (`-c') option to "
+              "specify a column.\n\nYou can either give it the column number "
+              "(couting from 1), or a match/search in its meta-data (e.g., "
+              "column names). For more information, please run the "
+              "following command (press the `SPACE' key to go down and "
+              "`q' to return to the command-line):\n\n"
+              "    $ info gnuastro \"Selecting table columns\"\n", name);
+    }
+  else
+    error(EXIT_FAILURE, 0, "no input file is specified");
 }
 
 
@@ -467,142 +377,190 @@ sanitycheck(struct statisticsparams *p)
 
 
 
+
 /**************************************************************/
 /***************       Preparations         *******************/
 /**************************************************************/
+static void
+ui_print_one_number(gal_data_t *data, int operator)
+{
+  char *toprint;
+  gal_data_t *out;
+  int flags=GAL_ARITHMETIC_NUMOK;
+  out=gal_arithmetic(operator, flags, data);
+  toprint=gal_data_write_to_string(out->array, out->type, 0);
+  printf("%s ", toprint);
+  gal_data_free(out);
+  free(toprint);
+}
+
+
+
+
+
 void
-preparearrays(struct statisticsparams *p)
+ui_print_one_row(struct statisticsparams *p)
 {
-  float min;
-  size_t s0, s1;
-  int bitpix, anyblank;
-  struct uiparams *up=&p->up;
-
-  /* Read the input and mask arrays: */
-  gal_fits_file_to_float(up->inputname, up->maskname, p->cp.hdu, up->mhdu,
-                              &p->img, &bitpix, &anyblank, &s0, &s1);
-  p->size=s0*s1;
-
-  /* If the minimum value is to be used as a mask then do it: */
-  if(p->ignoremin)
-    {
-      gal_statistics_float_min(p->img, p->size, &min);
-      gal_array_freplace_value(p->img, p->size, min, NAN);
-    }
+  struct gal_linkedlist_ill *tmp;
+
+  /* Reverse the list to have the same order the user wanted. */
+  gal_linkedlist_reverse_ill(&p->toprint);
+
+  /* Print the numbers. */
+  for(tmp=p->toprint; tmp!=NULL; tmp=tmp->next)
+    switch(tmp->v)
+      {
+      case ARGS_OPTION_KEY_NUMBER:
+        ui_print_one_number(p->input, GAL_ARITHMETIC_OP_NUMVAL); break;
+      case ARGS_OPTION_KEY_MINIMUM:
+        ui_print_one_number(p->input, GAL_ARITHMETIC_OP_MINVAL); break;
+      case ARGS_OPTION_KEY_MAXIMUM:
+        ui_print_one_number(p->input, GAL_ARITHMETIC_OP_MAXVAL); break;
+      case ARGS_OPTION_KEY_SUM:
+        ui_print_one_number(p->input, GAL_ARITHMETIC_OP_SUMVAL); break;
+      case ARGS_OPTION_KEY_MEAN:
+        ui_print_one_number(p->input, GAL_ARITHMETIC_OP_MEANVAL); break;
+      case ARGS_OPTION_KEY_STD:
+        ui_print_one_number(p->input, GAL_ARITHMETIC_OP_STDVAL); break;
+      case ARGS_OPTION_KEY_MEDIAN:
+        ui_print_one_number(p->input, GAL_ARITHMETIC_OP_MEDIANVAL); break;
+      case ARGS_OPTION_KEY_MODE:
+        error(EXIT_FAILURE, 0, "mode isn't implemented yet!");
+        break;
+      default:
+        error(EXIT_FAILURE, 0, "A bug! Operation code %d not recognized in "
+              "`ui_print_row'. Please contact us at %s so we can address "
+              "the problem", tmp->v, PACKAGE_BUGREPORT);
+      }
+
+  /* Print a new line. */
+  printf("\n");
+}
+
+
 
-  /* Move all the non-nan elements to the start of the array: */
-  gal_array_no_nans(p->img, &p->size);
 
-  /* Make a sorted array for most of the jobs: */
-  gal_array_float_copy(p->img, p->size, &p->sorted);
-  qsort(p->sorted, p->size, sizeof *p->sorted, gal_qsort_float_increasing);
 
-  /* Check the given range: */
-  if(p->histname || p->asciihist || p->mhistname || p->mirrorhist)
+static void
+ui_out_of_range_to_blank(struct statisticsparams *p)
+{
+  size_t one=1;
+  unsigned char flags=GAL_ARITHMETIC_NUMOK;
+  gal_data_t *tmp, *cond_g=NULL, *cond_l=NULL, *cond;
+  unsigned char flagsor = ( GAL_ARITHMETIC_FREE
+                            | GAL_ARITHMETIC_INPLACE
+                            | GAL_ARITHMETIC_NUMOK );
+
+  /* Set the condition. Note that the `greaterequal' name is for the data
+     we want. So we will set the condition based on those that are
+     less-than  */
+  if(!isnan(p->greaterequal))
     {
-      if(up->histquantset)
-        {
-          if(p->histquant>=0.5)
-            error(EXIT_FAILURE, 0, "the value to `--histquant' (-Q) must "
-                  "Be smaller than 0.5, because it sets the lower limit of "
-                  "the value range. The higher limit will be 1-Q");
-          p->histmin=
-            p->sorted[gal_statistics_index_from_quantile(p->size,
-                                                         p->histquant)];
-          p->histmax=
-            p->sorted[gal_statistics_index_from_quantile(p->size,
-                                                         1 - p->histquant)];
-        }
-      else
-        {
-          switch(up->histminset+up->histmaxset)
-            {
-            case 0:
-              p->histmin=p->sorted[0];
-              p->histmax=p->sorted[p->size-1];
-              break;
-            case 1:
-              error(EXIT_FAILURE, 0, "the options `--histmin' (-i) and "
-                    "`--histmax' (-x) should both be specified. You have "
-                    "only given the %s"GAL_STRINGS_HOW_TO_CHECK_VALUES,
-                    up->histminset==1 ? "former" : "latter");
-              break;
-            case 2:
-              if(p->histmin>=p->histmax)
-                error(EXIT_FAILURE, 0, "the value to `--histmin' (-i) (%f) "
-                      "is larger or equal to that of `--histmax' (-x) (%f)"
-                      GAL_STRINGS_HOW_TO_CHECK_VALUES, p->histmin, p->histmax);
-              if(p->histmin>p->sorted[p->size-1] || p->histmax<p->sorted[0])
-                error(EXIT_FAILURE, 0, "the range of data is %.5f to %.5f. "
-                      "However, you have set `--histmin' (-i) and "
-                      "`--histmax' (-x) to %.5f and %.5f respectively. "
-                      "They do not overlap"GAL_STRINGS_HOW_TO_CHECK_VALUES,
-                      p->sorted[0], p->sorted[p->size-1], p->histmin,
-                      p->histmax);
-              break;
-            default:
-              error(EXIT_FAILURE, 0, "a bug! Please contact us at "
-                    PACKAGE_BUGREPORT" So we can solve the problem. the "
-                    "value of up->histminset+up->histmaxset is not 0, 1 or "
-                    "2");
-            }
-        }
+      tmp=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT32, 1, &one, NULL, 0, -1,
+                        NULL, NULL, NULL);
+      *((float *)(tmp->array)) = p->greaterequal;
+      cond_g=gal_arithmetic(GAL_ARITHMETIC_OP_LT, flags, p->input, tmp);
+      gal_data_free(tmp);
     }
-  else                          /* For the ascii histogram. */
+
+  /* Same reasoning as above for `p->greaterthan'. */
+  if(!isnan(p->lessthan))
     {
-      p->histmin=p->sorted[0];
-      p->histmax=p->sorted[p->size-1];
+      tmp=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT32, 1, &one, NULL, 0, -1,
+                        NULL, NULL, NULL);
+      *((float *)(tmp->array)) = p->lessthan;
+      cond_l=gal_arithmetic(GAL_ARITHMETIC_OP_GE, flags, p->input, tmp);
+      gal_data_free(tmp);
     }
 
-  if(p->cfpname && p->cfpsimhist==0)
+  /* Now, set the final condition. If both values were specified, then use
+     the GAL_ARITHMETIC_OP_OR to merge them into one. */
+  switch( !isnan(p->greaterequal) + !isnan(p->lessthan) )
     {
-      if(up->cfpquantset)
-        {
-          if(p->cfpquant>=0.5)
-            error(EXIT_FAILURE, 0, "the value to `--cfpquant' (-U) must "
-                  "Be smaller than 0.5, because it sets the lower limit of "
-                  "the value range. The higher limit will be 1-U");
-          p->cfpmin=p->sorted[gal_statistics_index_from_quantile(p->size,
-                                                                 p->cfpquant)];
-          p->cfpmax=
-            p->sorted[gal_statistics_index_from_quantile(p->size,
-                                                         1 - p->cfpquant)];
-        }
-      else
+    case 0: return;             /* No condition was specified, return.  */
+    case 1:                     /* Only one condition was specified.    */
+      cond = isnan(p->greaterequal) ? cond_l : cond_g;
+      break;
+    case 2:
+      cond = gal_arithmetic(GAL_ARITHMETIC_OP_OR, flagsor, cond_l, cond_g);
+      break;
+    }
+
+  /* Allocate a blank value to mask all input pixels. */
+  tmp=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT32, 1, &one, NULL,
+                     0, -1, NULL, NULL, NULL);
+  *((float *)(tmp->array)) = NAN;
+
+  /* Set all the pixels that satisfy the condition to blank. */
+  gal_arithmetic(GAL_ARITHMETIC_OP_WHERE, flagsor, p->input, cond, tmp);
+}
+
+
+
+
+
+void
+ui_preparations(struct statisticsparams *p)
+{
+  gal_data_t *tmp;
+  char *errorstring;
+  float *f, *ff, *fo;
+  size_t num, numcolmatches=0;
+  struct gal_linkedlist_stll *column=NULL;
+
+  /* Read the input: no matter if its an image or a table column. */
+  if(p->isfits && p->hdu_type==IMAGE_HDU)
+    p->input=gal_fits_img_read(p->inputname, p->cp.hdu, p->cp.minmapsize);
+  else
+    {
+      /* Read the input column. */
+      gal_linkedlist_add_to_stll(&column, p->column, 0);
+      p->input=gal_table_read(p->inputname, p->cp.hdu, column, p->cp.searchin,
+                              p->cp.ignorecase, p->cp.minmapsize);
+      if(p->input->next)
         {
-          switch(up->cfpminset+up->cfpmaxset)
-            {
-            case 0:
-              p->cfpmin=p->sorted[0];
-              p->cfpmax=p->sorted[p->size-1];
-              break;
-            case 1:
-              error(EXIT_FAILURE, 0, "the options `--cfpmin' (-a) and "
-                    "`--cfpmax' (-b) should both be specified. You have "
-                    "only given the %s"GAL_STRINGS_HOW_TO_CHECK_VALUES,
-                    up->cfpminset==1 ? "former" : "latter");
-              break;
-            case 2:
-              if(p->cfpmin>p->cfpmax)
-                error(EXIT_FAILURE, 0, "the value to `--cfpmin' (-a) (%.f) "
-                      "is larger than that of `--cfpmax' (-b) (%f)"
-                      GAL_STRINGS_HOW_TO_CHECK_VALUES, p->cfpmin, p->cfpmax);
-              if(p->cfpmin>p->sorted[p->size-1] || p->cfpmax<p->sorted[0])
-                error(EXIT_FAILURE, 0, "the range of data is %.5f to %.5f. "
-                      "However, you have set `--cfpmin' (-a) and "
-                      "`--cfpmax' (-b) to %.5f and %.5f respectively. "
-                      "They do not overlap"GAL_STRINGS_HOW_TO_CHECK_VALUES,
-                      p->sorted[0], p->sorted[p->size-1], p->cfpmin,
-                      p->cfpmax);
-              break;
-            default:
-              error(EXIT_FAILURE, 0, "a bug! Please contact us at "
-                    PACKAGE_BUGREPORT" So we can solve the problem. the "
-                    "value of up->cfpminset+up->cfpmaxset is not 0, 1 or "
-                    "2");
-            }
+          for(tmp=p->input;tmp!=NULL;tmp=tmp->next) ++numcolmatches;
+          asprintf(&errorstring, "%zu columns were selected with `%s' "
+                   "(value to `--column' option). In this context, "
+                   "Statistics can only work on one data-set (column in a "
+                   "table).", numcolmatches, p->column);
+          gal_table_error_col_selection(p->inputname, p->cp.hdu, errorstring);
         }
+
+      /* Clean up. */
+      gal_linkedlist_free_stll(column, 0);
     }
+
+  /* Set the out-of-range values in the input to blank. */
+  ui_out_of_range_to_blank(p);
+
+  /* Print the one-row numbers if the user asked for them. */
+  if(p->toprint) ui_print_one_row(p);
+
+  /* For the main body of the program, we will assume the data is in
+     float32. */
+  p->input=gal_data_copy_to_new_type_free(p->input, GAL_DATA_TYPE_FLOAT32);
+
+  /* Put all the blank elements in the end for easier sorting */
+  ff=(fo=f=p->input->array)+p->input->size;
+  num=0; do if(!isnan(*f)) {*fo++=*f; ++num;} while(++f<ff);
+  f=((float *)(p->input->array))+num; do *f++=NAN; while(f<ff);
+
+  /* Sort the non-blank elements. */
+  qsort(p->input->array, num, sizeof(float), gal_qsort_float32_increasing);
+
+  /* For a check:
+  {
+    size_t i;
+    float *f=p->input->array;
+    for(i=0;i<p->input->size;++i) printf("%f\n", f[i]);
+  }
+  */
+
+  /* Remove the dimensionality of the array and correct its size to only be
+     the non-blank elements. */
+  p->input->ndim=1;
+  p->input->dsize[0]=p->input->size=num;
 }
 
 
@@ -626,53 +584,57 @@ preparearrays(struct statisticsparams *p)
 /**************************************************************/
 /************         Set the parameters          *************/
 /**************************************************************/
+
 void
-setparams(int argc, char *argv[], struct statisticsparams *p)
+ui_read_check_inputs_setup(int argc, char *argv[], struct statisticsparams *p)
 {
-  struct gal_commonparams *cp=&p->cp;
-
-  /* Set the non-zero initial values, the structure was initialized to
-     have a zero value for all elements. */
-  cp->spack         = SPACK;
-  cp->verb          = 1;
-  cp->numthreads    = num_processors(NPROC_CURRENT);
-  cp->removedirinfo = 1;
-
-  p->asciihist      = 1;
-  p->sigclip        = 1;
-  p->mirror         = NAN;
-  p->onebinvalue    = NAN;
-  p->histname=p->cfpname="a";   /* Will be set later, just a sign that */
-                                /* they should be output.              */
-  /* Read the arguments. */
+  struct gal_options_common_params *cp=&p->cp;
+
+
+  /* Include the parameters necessary for argp from this program (`args.h')
+     and for the common options to all Gnuastro (`commonopts.h'). We want
+     to directly put the pointers to the fields in `p' and `cp', so we are
+     simply including the header here to not have to use long macros in
+     those headers which make them hard to read and modify. This also helps
+     in having a clean environment: everything in those headers is only
+     available within the scope of this function. */
+#include <commonopts.h>
+#include "args.h"
+
+
+  /* Initialize the options and necessary information.  */
+  ui_initialize_options(p, program_options, gal_commonopts_options);
+
+
+  /* Read the command-line options and arguments. */
   errno=0;
   if(argp_parse(&thisargp, argc, argv, 0, 0, p))
     error(EXIT_FAILURE, errno, "parsing arguments");
 
-  /* Add the user default values and save them if asked. */
-  GAL_CONFIGFILES_CHECK_SET_CONFIG;
 
-  /* Check if all the required parameters are set. */
-  checkifset(p);
+  /* Read the configuration files and set the common values. */
+  gal_options_read_config_set(&p->cp);
 
-  /* Print the values for each parameter. */
-  if(cp->printparams)
-    GAL_CONFIGFILES_REPORT_PARAMETERS_SET;
 
-  /* Do a sanity check. */
-  sanitycheck(p);
+  /* Read the options into the program's structure, and check them and
+     their relations prior to printing. */
+  ui_read_check_only_options(p);
 
-  /* Make the array of input images. */
-  preparearrays(p);
 
-  /* Everything is ready, notify the user of the program starting. */
-  if(cp->verb)
-    {
-      printf(SPACK_NAME" started on %s", ctime(&p->rawtime));
-      printf("  - Input read: %s (hdu: %s)\n", p->up.inputname, p->cp.hdu);
-      if(p->up.maskname)
-        printf("  - Mask read: %s (hdu: %s)\n", p->up.maskname, p->up.mhdu);
-    }
+  /* Print the option values if asked. Note that this needs to be done
+     after the option checks so un-sane values are not printed in the
+     output state. */
+  gal_options_print_state(&p->cp);
+
+
+  /* Check that the options and arguments fit well with each other. Note
+     that arguments don't go in a configuration file. So this test should
+     be done after (possibly) printing the option values. */
+  ui_check_options_and_arguments(p);
+
+
+  /* Read/allocate all the necessary starting arrays. */
+  ui_preparations(p);
 }
 
 
@@ -698,20 +660,9 @@ setparams(int argc, char *argv[], struct statisticsparams 
*p)
 /************      Free allocated, report         *************/
 /**************************************************************/
 void
-freeandreport(struct statisticsparams *p, struct timeval *t1)
+ui_free_report(struct statisticsparams *p)
 {
   /* Free the allocated arrays: */
-  free(p->img);
-  free(p->sorted);
   free(p->cp.hdu);
-  free(p->cfpname);
-  free(p->histname);
-  free(p->mcfpname);
-  free(p->mhistname);
   free(p->cp.output);
-  if(p->up.masknameallocated) free(p->up.maskname);
-
-  /* Print the final message. */
-  if(p->cp.verb)
-    gal_timing_report(t1, SPACK_NAME" finished in: ", 0);
 }
diff --git a/bin/statistics/ui.h b/bin/statistics/ui.h
index dcfa17d..88bd088 100644
--- a/bin/statistics/ui.h
+++ b/bin/statistics/ui.h
@@ -1,5 +1,5 @@
 /*********************************************************************
-Statistics - Get general statistics about the image.
+Statistics - Statistical analysis on input dataset.
 Statistics is part of GNU Astronomy Utilities (Gnuastro) package.
 
 Original author:
@@ -20,13 +20,14 @@ 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 IMCROPUI_H
-#define IMCROPUI_H
+#ifndef UI_H
+#define UI_H
 
 void
-setparams(int argc, char *argv[], struct statisticsparams *p);
+ui_read_check_inputs_setup(int argc, char *argv[],
+                           struct statisticsparams *p);
 
 void
-freeandreport(struct statisticsparams *p, struct timeval *t1);
+ui_free_report(struct statisticsparams *p);
 
 #endif
diff --git a/bin/warp/ui.c b/bin/warp/ui.c
index 10120a2..771ff5c 100644
--- a/bin/warp/ui.c
+++ b/bin/warp/ui.c
@@ -239,17 +239,11 @@ static void *
 ui_add_to_modular_warps_ll(struct argp_option *option, char *arg,
                            char *filename, size_t lineno, void *params)
 {
-  size_t i, num=0;
+  size_t i;
+  double tmp;
   gal_data_t *new;
-  char *tailptr, *c=arg;
-  double numerator=NAN, denominator=NAN, tmp;
-  struct gal_linkedlist_dll *list=NULL, *tdll;
   struct warpparams *p=(struct warpparams *)params;
 
-  /* The nature of the arrays/numbers read here is very small, so since
-     `p->cp.minmapsize' might not have been read yet, we will set it to -1
-     (largest size_t number), so the values are kept in memory. */
-  size_t minmapsize=-1;
 
   /* Parse the (possible) arguments. */
   if(option->key==ARGS_OPTION_KEY_ALIGN)
@@ -268,98 +262,11 @@ ui_add_to_modular_warps_ll(struct argp_option *option, 
char *arg,
       if(arg && *arg=='0') return NULL;
 
       /* Allocate the data structure. */
-      new=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT64, 0, &num, NULL, 0,
-                         minmapsize, NULL, NULL, NULL);
+      new=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT64, 0, NULL, NULL, 0,
+                         -1, NULL, NULL, NULL);
     }
-  else
-    {
-      /* Go through the input character by character. */
-      while(*c!='\0')
-        {
-        switch(*c)
-          {
-          /* Ignore space or tab. */
-          case ' ':
-          case '\t':
-            ++c;
-            break;
-
-
-          /* Comma marks the transition to the next number. */
-          case ',':
-            if(isnan(numerator))
-              error_at_line(EXIT_FAILURE, 0, filename, lineno, "a number "
-                            "must be given before `,'. You have given: `%s'",
-                            arg);
-            gal_linkedlist_add_to_dll(&list, isnan(denominator)
-                                      ? numerator : numerator/denominator);
-            numerator=denominator=NAN;
-            ++num;
-            ++c;
-            break;
-
-
-          /* Divide two numbers. */
-          case '/':
-            if( isnan(numerator) || !isnan(denominator) )
-              error_at_line(EXIT_FAILURE, 0, filename, lineno, "`/' must "
-                            "only be between two numbers and used for "
-                            "division. But you have given `%s'. This "
-                            "was a value to the `%s' option", arg,
-                            option->name);
-            ++c;
-            break;
-
-
-          /* Read the number. */
-          default:
-
-            /* Parse the string. */
-            tmp=strtod(c, &tailptr);
-            if(tailptr==c)
-              error_at_line(EXIT_FAILURE, 0, filename, lineno, "the first "
-                            "part of `%s' couldn't be read as a number. This "
-                            "was part of `%s' (value to the `%s' option)", c,
-                            arg, option->name);
-
-            /* See if the number should be put in the numerator or
-               denominator. */
-            if(isnan(numerator)) numerator=tmp;
-            else
-              {
-                if(isnan(denominator)) denominator=tmp;
-                else error_at_line(EXIT_FAILURE, 0, filename, lineno, "more "
-                                "than two numbers in each element.");
-              }
-
-            /* Set `c' to tailptr. */
-            c=tailptr;
-          }
-        }
+  else new=gal_options_parse_list_of_numbers(arg, filename, lineno);
 
-      /* If the last number, wasn't finished by a `,', add the read value
-         to the list */
-      if( !isnan(numerator) )
-        {
-          ++num;
-          gal_linkedlist_add_to_dll(&list, isnan(denominator)
-                                    ? numerator : numerator/denominator);
-        }
-
-      /* Allocate the new data structure and fill it up. */
-      i=num;
-      new=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT64, 1, &num, NULL, 0,
-                         minmapsize, NULL, NULL, NULL);
-      for(tdll=list;tdll!=NULL;tdll=tdll->next)
-        ((double *)new->array)[--i]=tdll->v;
-      gal_linkedlist_free_dll(list);
-    }
-
-  /* For a check.
-  printf("%s (%s): %zu number(s)\n", option->name, arg, num);
-  for(i=0;i<num;++i)
-    printf("\t%.15f\n", ((double *)new->array)[i]);
-  */
 
   /* If this was a matrix, then put it in the matrix element of the main
      data structure. Otherwise, add the list of given values to the modular
@@ -370,10 +277,10 @@ ui_add_to_modular_warps_ll(struct argp_option *option, 
char *arg,
       if(p->matrix)
         error_at_line(EXIT_FAILURE, 0, filename, lineno, "only one matrix "
                       "may be given, you can use multiple modular warpings");
-      if(num!=4 && num!=9)
+      if(new->size!=4 && new->size!=9)
         error_at_line(EXIT_FAILURE, 0, filename, lineno, "only a 4 or 9 "
                       "element `matrix' is currently acceptable. `%s' has "
-                      "%zu elements", arg, num);
+                      "%zu elements", arg, new->size);
 
       /* Keep the matrix in the main structure. */
       p->matrix=new;
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 9a92e15..f495497 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -522,7 +522,6 @@ Review of library fundamentals
 Gnuastro library
 
 * Configuration information::   General information about library config.
-* Array manipulation::          Functions for manipulating arrays.
 * Bounding box::                Finding the bounding box.
 * FITS files::                  Working with FITS data.
 * Git wrappers::                Wrappers for functions in libgit2.
@@ -533,7 +532,6 @@ Gnuastro library
 * Spatial convolution::         Doing spatial convolution on an image
 * Statistical operations::      Functions for basic statistics.
 * Multithreaded programming::   Facilitating use of pthreads.
-* Text table into a C array::   Read an ASCII text table into a C array.
 * World Coordinate System::     Dealing with the world coordinate system.
 
 FITS files (@file{fits.h})
@@ -7878,15 +7876,34 @@ to @mymath{\log(4)}. The output type is determined from 
the input, see the
 explanation under @command{sqrt} for more.
 
 @item minvalue
-Minimum value in the top operand on the stack, so address@hidden
-minvalue}'' will push the the minimum pixel value in this image onto the
-stack. Therefore this operator is mainly intended for data (for example
-images), if the top operand is a number, this operator just returns it
-without any change. So note that when this operator acts on a single image,
-the output will no longer be an image, but a number.
+Minimum (non-blank) value in the top operand on the stack, so
address@hidden minvalue}'' will push the the minimum pixel value in this
+image onto the stack. Therefore this operator is mainly intended for data
+(for example images), if the top operand is a number, this operator just
+returns it without any change. So note that when this operator acts on a
+single image, the output will no longer be an image, but a number.
 
 @item maxvalue
-Maximum value of first operand, similar to @command{minvalue}.
+Maximum (non-blank) value of first operand, similar to @command{minvalue}.
+
address@hidden numvalue
+Number of non-blank elements in first operand, similar to
address@hidden
+
address@hidden sumvalue
+Sum of non-blank elements in first operand, similar to @command{minvalue}.
+
address@hidden meanvalue
+Mean value of non-blank elements in first operand, similar to
address@hidden
+
address@hidden stdvalue
+Standard deviation of non-blank elements in first operand, similar to
address@hidden
+
address@hidden medianvalue
+Median of non-blank elements in first operand, similar to
address@hidden
 
 @cindex NaN
 @item min
@@ -7925,13 +7942,21 @@ conversion will be used.
 Similar to @command{min}, but the pixels of the output will contain
 the maximum of the respective pixels in all operands in the stack.
 
address@hidden num
+Similar to @command{min}, but the pixels of the output will contain the
+number of the respective non-blank pixels in all input operands.
+
 @item sum
 Similar to @command{min}, but the pixels of the output will contain the sum
 of the respective pixels in all input operands.
 
address@hidden average
-Similar to @command{min}, but the pixels of the output will contain
-the average of the respective pixels in all operands in the stack.
address@hidden mean
+Similar to @command{min}, but the pixels of the output will contain the
+mean (average) of the respective pixels in all operands in the stack.
+
address@hidden std
+Similar to @command{min}, but the pixels of the output will contain the
+standard deviation of the respective pixels in all operands in the stack.
 
 @item median
 Similar to @command{min}, but the pixels of the output will contain
@@ -8005,20 +8030,22 @@ format is demonstrated in this simplified example:
 $ astarithmetic modify.fits cond.fits if-true.fits where
 @end example
 
-The third and second popped operands (@file{modify.fits} and
address@hidden above, see @ref{Reverse polish notation}) have to be
-images or all three operands have to be numbers. In the former case (when
-the first two are images), the third popped operand can be an image or a
-number. The value of any pixel in @file{out.fits} that corresponds to a
-non-zero pixel of @file{cond.fits} will be changed to the value of the same
-pixel in @file{if-true.fits} (or a fixed number). Note than the
address@hidden has to have @code{unsigned char} type.
-
-Commonly you won't be dealing with an actual FITS file of a conditional
-image. You will probably define the condition in the same run based on some
-other reference image and use the conditional and logical operators above
-to make a true/false (or one/zero) image for you internally. For example
-the case below:
+The value of any pixel in @file{out.fits} that corresponds to a non-zero
+pixel of @file{cond.fits} will be changed to the value of the same pixel in
address@hidden (or a fixed number). The third and second popped
+operands (@file{modify.fits} and @file{cond.fits} above, see @ref{Reverse
+polish notation}) have to have the same
+dimensions/size. @file{if-true.fits} can be either a number, or have the
+same dimension/size as the other two.
+
+That @file{cond.fits} has to have @code{uint8} (or @code{unsigned char} in
+standard C) type (see @ref{Numeric data types} and the type conversion
+operators below). However, ommonly this is no problem because in most cases
+you won't be dealing with an actual FITS file of a conditional image. You
+will probably define the condition in the same run based on some other
+reference image and use the conditional and logical operators above to make
+a true/false (or one/zero) image for you internally. For example the case
+below:
 
 @example
 $ astarithmetic in.fits reference.fits 100 gt new.fits where
@@ -8030,6 +8057,16 @@ corresponding pixel in @file{new.fits}. Finally, the 
input operands are
 read and used independently, so you can use the same file more than once as
 any of the operands.
 
+When the first popped operand (@file{if-true.fits}) is a single number, it
+can be a NaN value (or any blank value, depending on its type) like the
+example below (see @ref{Blank pixels}). In this example, all the pixels in
address@hidden that have a value greater than 100, will become blank
+in the natural data type of @file{in.fits}.
+
address@hidden
+$ astarithmetic in.fits reference.fits 100 gt nan where
address@hidden example
+
 @item bitand
 Bitwise AND operator: only bits with values of 1 in both popped operands
 will get the value of 1, the rest will be set to 0. For example (assuming
@@ -15862,7 +15899,6 @@ problems. It will stabilize with the removal of this 
notice. Check the
 
 @menu
 * Configuration information::   General information about library config.
-* Array manipulation::          Functions for manipulating arrays.
 * Bounding box::                Finding the bounding box.
 * FITS files::                  Working with FITS data.
 * Git wrappers::                Wrappers for functions in libgit2.
@@ -15873,11 +15909,10 @@ problems. It will stabilize with the removal of this 
notice. Check the
 * Spatial convolution::         Doing spatial convolution on an image
 * Statistical operations::      Functions for basic statistics.
 * Multithreaded programming::   Facilitating use of pthreads.
-* Text table into a C array::   Read an ASCII text table into a C array.
 * World Coordinate System::     Dealing with the world coordinate system.
 @end menu
 
address@hidden Configuration information, Array manipulation, Gnuastro library, 
Gnuastro library
address@hidden Configuration information, Bounding box, Gnuastro library, 
Gnuastro library
 @subsection Configuration information (@file{config.h})
 
 The @file{gnuastro/config.h} header contains information about the full
@@ -15925,186 +15960,10 @@ value of @code{0}. see @ref{Implementation of 
pthread_barrier} for more.
 @end deffn
 
 
address@hidden Array manipulation, Bounding box, Configuration information, 
Gnuastro library
address@hidden Array manipulation (@file{array.h})
-
-When working on arrays, certain operations are commonly necessary. It is
-very easy to write a loop for such operations, however such loops can
-unnecessarily make the code harder to write, read, and debug. So, it is
-much more cleaner if these operations are defined elsewhere. The
-declarations of these functions are available in @file{gnuastro/array.h}.
-
-In C, an array is a contiguous region of memory. So for most operations
-here, the dimentionality of the data is irrelevant, we just need the total
-size of the array. Please note that these functions are still very
-rudimentary, primarily intended for Gnuastro's internal operations. In
-future releases, these functions will be updated to be more suitable for a
-generic library.
-
address@hidden void gal_array_uchar_init_on_region (unsigned char @code{*in}, 
const unsigned char @code{v}, size_t @code{start}, size_t @code{s0}, size_t 
@code{s1}, size_t @code{is1})
-Initialize a region of the array @code{in} to the value @code{v}. The
-region is defined by the 1D index of the first pixel (@code{start}) and its
-C array (opposite of FITS array) width (@code{s0}) and height
-(@code{s1}). The full C array width is defined by @code{is1}.
address@hidden deftypefun
-
address@hidden void gal_array_long_init (long @code{*in}, size_t @code{size}, 
const long @code{v})
-Initialize the long type array @code{in} with @code{size} elements to the
-value @code{v}.
address@hidden deftypefun
-
address@hidden void gal_array_long_init_on_region (long @code{*in}, const long 
@code{v}, size_t @code{start}, size_t @code{s0}, size_t @code{s1}, size_t 
@code{is1})
-Similar to @code{gal_array_uchar_init_on_region} but for long type arrays.
address@hidden deftypefun
-
address@hidden void gal_array_uchar_copy (unsigned char @code{*in}, size_t 
@code{size}, unsigned char @code{**out})
-Allocate space for and copy the input array @code{in} with @code{size}
-elements into an output array (@code{out}).
address@hidden deftypefun
-
address@hidden void gal_array_float_copy (float @code{*in}, size_t @code{size}, 
float @code{**out})
-Similar to @code{gal_array_uchar_copy} but for the single precision
-floating point type.
address@hidden deftypefun
-
address@hidden void gal_array_float_copy_noalloc (float @code{*in}, size_t 
@code{size}, float @code{*out})
-Similar to @code{gal_array_float_copy} but the @code{out} array must
-already be allocated.
address@hidden deftypefun
-
address@hidden void gal_array_fset_const (float @code{*in}, size_t @code{size}, 
float @code{a})
-Set the input array @code{in} with @code{size} elements to the constant
-value @code{a}.
address@hidden deftypefun
-
address@hidden void gal_array_freplace_value (float @code{*in}, size_t 
@code{size}, float @code{from}, float @code{to})
address@hidden NaN
-Replace all elements with value @code{from} to value @code{to} in the
address@hidden array with @code{size} elements. @code{from} can be a NaN value,
-in this case, C's @code{isnan} function is used, not @code{==}.
address@hidden deftypefun
-
address@hidden void gal_array_freplace_nonnans (float @code{*in}, size_t 
@code{size}, float @code{to})
-Replace any non-NaN value in the array @code{in} with @code{size} elements
-to the value @code{to}.
address@hidden deftypefun
-
address@hidden void gal_array_no_nans (float @code{*in}, size_t @code{*size})
-Move all the non-NaN elements in the array @code{in} to the start of the
-array so that the non-NaN elements are contiguous. This is useful for cases
-where you want to sort the data. Note that before this function,
address@hidden must point to the initial size of the array. After this
-function returns, size will point to the new size of the array, with only
-non-Nan elements.
address@hidden deftypefun
-
address@hidden void gal_array_no_nans_double (double @code{*in}, size_t 
@code{*size})
-Similar to @code{gal_array_no_nans} but for double types.
address@hidden deftypefun
-
address@hidden void gal_array_uchar_replace (unsigned char @code{*in}, size_t 
@code{size}, unsigned char @code{from}, unsigned char @code{to})
-Replace all elements with value @code{from} to value @code{to} in the
address@hidden array with @code{size} elements.
address@hidden deftypefun
-
address@hidden void gal_array_fmultip_const (float @code{*in}, size_t 
@code{size}, float @code{a})
-Multiply all elements in the float array @code{in} with size @code{size} by
-the constant @code{a}.
address@hidden deftypefun
-
address@hidden void gal_array_fsum_const (float @code{*in}, size_t @code{size}, 
float @code{a})
-Add all elements in the float array @code{in} with size @code{size} by the
-constant @code{a}.
address@hidden deftypefun
-
address@hidden {float *} gal_array_fsum_arrays (float @code{*in1}, float 
@code{*in2}, size_t @code{size})
-Add the two arrays @code{in1} and @code{in2} and keep the result in a newly
-allocated array. The returned value is a pointer to that array.
address@hidden deftypefun
-
address@hidden void gal_array_dmultip_const (double @code{*in}, size_t 
@code{size}, double @code{a})
-Multiply the double type array @code{in} with @code{size} elements with the
-constant @code{a}.
address@hidden deftypefun
-
address@hidden void gal_array_dmultip_arrays (double @code{*in1}, double 
@code{*in2}, size_t @code{size})
-Multiply the two double type arrays @code{in1} with @code{in2} (both with
address@hidden elements) with each other and store the output in the first
-array.
address@hidden deftypefun
-
address@hidden void gal_array_ddivide_const (double @code{*in}, size_t 
@code{size}, double @code{a})
-Divide the double type array @code{in} with @code{size} elements by the
-constant @code{a}.
address@hidden deftypefun
-
address@hidden void gal_array_dconst_divide (double @code{*in}, size_t 
@code{size}, double @code{a})
-Divide the constant @code{a} by the array @code{in} with @code{size} and
-store the result for each pixel in the array.
address@hidden deftypefun
-
address@hidden void gal_array_ddivide_arrays (double @code{*in1}, double 
@code{*in2}, size_t @code{size})
-Divide each element of the two double array @code{in1} by the elements in
address@hidden (both with @code{size} elements) and store the output in the
-first array.
address@hidden deftypefun
 
address@hidden void gal_array_dsum_const (double @code{*in}, size_t 
@code{size}, double @code{a})
-Add the double type array @code{in} with @code{size} elements with the
-constant @code{a}.
address@hidden deftypefun
 
address@hidden void gal_array_dsum_arrays (double @code{*in1}, double 
@code{*in2}, size_t @code{size})
-Add the two double type arrays @code{in1} with @code{in2} (both with
address@hidden elements) with each other and store the output in the first
-array.
address@hidden deftypefun
-
address@hidden void gal_array_dsubtract_const (double @code{*in}, size_t 
@code{size}, double @code{a})
-Subtract the double type array @code{in} with @code{size} elements from the
-constant @code{a}.
address@hidden deftypefun
 
address@hidden void gal_array_dconst_subtract (double @code{*in}, size_t 
@code{size}, double @code{a})
-Subtract the constant @code{a} from the array @code{in} with @code{size}
-and store the result for each pixel in the array.
address@hidden deftypefun
-
address@hidden void gal_array_dsubtract_arrays (double @code{*in1}, double 
@code{*in2}, size_t @code{size})
-Subtract each element of th the double type array @code{in2} from
address@hidden and put the result in the first array.
address@hidden deftypefun
-
address@hidden void gal_array_dpower_const (double @code{*in}, size_t 
@code{size}, double @code{a})
-Each element in the intput array is taken to the power of of @code{a} and
-stored in place.
address@hidden deftypefun
-
address@hidden void gal_array_dconst_power (double @code{*in}, size_t 
@code{size}, double @code{a})
-The constant @code{a} is taken to the power of each element in the intput
-array (@code{in}) and stored in place.
address@hidden deftypefun
-
address@hidden void gal_array_dpower_arrays (double @code{*in1}, double 
@code{*in2}, size_t @code{size})
-Each element in the first input array is taken to the power of the same
-element in the second array.
address@hidden deftypefun
-
address@hidden void gal_array_dlog_array (double @code{*in1}, size_t 
@code{size})
-The natural logarithm of the input array is taken and stored in place.
address@hidden deftypefun
-
address@hidden void gal_array_dlog10_array (double @code{*in1}, size_t 
@code{size})
-The base-10 logarithm of the input array is taken and stored in place.
address@hidden deftypefun
-
address@hidden void gal_array_dabs_array (double @code{*in1}, size_t 
@code{size})
-Replace each element of the input array with its absolute value.
address@hidden deftypefun
-
-
-
address@hidden Bounding box, FITS files, Array manipulation, Gnuastro library
address@hidden Bounding box, FITS files, Configuration information, Gnuastro 
library
 @subsection Bounding box (@file{box.h})
 
 Functions related to reporting a the bounding box of certain inputs are
@@ -17459,161 +17318,6 @@ The maximum number of times to try for 
@mymath{\sigma}-clipping (see
 @ref{Sigma clipping}).
 @end deffn
 
address@hidden void gal_statistics_long_non_blank_min (long @code{*in}, size_t 
@code{size}, long @code{*min}, long  @code{blank})
-Find the the minimum (non-blank) value in the @code{in} array (with
address@hidden elements). The long type doesn't have a NaN value like the
-float types, see @ref{Blank pixels}. So as blank pixels, a value in the
-range of acceptable values (@code{blank} must be given so it is explicitly
-ignored. You can use @code{GAL_FITS_LONG_BLANK} in @file{gnuastro/fits.h}.
address@hidden deftypefun
-
address@hidden void gal_statistics_long_non_blank_max (long @code{*in}, size_t 
@code{size}, long @code{*max}, long @code{blank})
-Similar to @code{gal_statistics_long_non_blank_min}, but find the maximum.
address@hidden deftypefun
-
address@hidden void gal_statistics_float_min (float @code{*in}, size_t 
@code{size}, float @code{*min})
-Find the minimum value in the single precision floating point @code{in}
-array (with @code{size} elements) and store the value in @code{min}. Note
-that for all floating point types, any blank (NaN) value will be
-automatically ignored in this function.
address@hidden deftypefun
-
address@hidden void gal_statistics_float_max (float @code{*in}, size_t 
@code{size}, float @code{*max})
-Similar to @code{gal_statistics_float_min} but find the maximum.
address@hidden deftypefun
-
address@hidden void gal_statistics_double_min (double @code{*in}, size_t 
@code{size}, double @code{*min})
-Similar to @code{gal_statistics_float_min} but for double precision
-floating point types.
address@hidden deftypefun
-
address@hidden double gal_statistics_double_min_return (double @code{*in}, 
size_t @code{size})
-Similar to @code{gal_statistics_double_min} but minimum will be returned,
-not put in a pointer.
address@hidden deftypefun
-
address@hidden void gal_statistics_double_max (double @code{*in}, size_t 
@code{size}, double @code{*max})
-Similar to @code{gal_statistics_double_min}, but find maximum value.
address@hidden deftypefun
-
address@hidden double gal_statistics_double_max_return (double @code{*in}, 
size_t @code{size})
-Similar to @code{gal_statistics_double_max_return}, but return maximum.
address@hidden deftypefun
-
address@hidden void gal_statistics_float_max_masked (float @code{*in}, unsigned 
char @code{*mask}, size_t @code{size}, float @code{*max})
-Only find maximum where @code{mask} has a value of @code{0}.
address@hidden deftypefun
-
address@hidden void gal_statistics_float_second_max (float @code{*in}, size_t 
@code{size}, float @code{*secondmax})
-Find the second largest value in the array.
address@hidden deftypefun
-
address@hidden void gal_statistics_float_second_min (float @code{*in}, size_t 
@code{size}, float @code{*secondmin})
-Find the second smallest value in the array.
address@hidden deftypefun
-
address@hidden void gal_statistics_f_min_max (float @code{*in}, size_t 
@code{size}, float @code{*min}, float @code{*max})
-Find the minimum and maximum simultaneously for single precision floating
-point types.
address@hidden deftypefun
-
address@hidden void gal_statistics_d_min_max (double @code{*in}, size_t 
@code{size}, double @code{*min}, double @code{*max})
-Similar to @code{gal_statistics_f_min_max}, but for double precision
-floating point types.
address@hidden deftypefun
-
address@hidden void gal_statistics_f_max_with_index (float @code{*in}, size_t 
@code{size}, float @code{*max}, size_t @code{*index})
-Find the maximum value in the array along with its index for a single
-precision floating point array.
address@hidden deftypefun
-
address@hidden void gal_statistics_d_max_with_index (double @code{*in}, size_t 
@code{size}, double @code{*max}, size_t @code{*index})
-Similar to @code{gal_statistics_f_max_with_index}, but for a double
-precision floating point array.
address@hidden deftypefun
-
address@hidden void gal_statistics_d_min_with_index (double @code{*in}, size_t 
@code{size}, double @code{*min}, size_t @code{*index})
address@hidden deftypefun
-
address@hidden void gal_statistics_f_min_with_index (float @code{*in}, size_t 
@code{size}, float @code{*min}, size_t @code{*index})
address@hidden deftypefun
-
address@hidden float gal_statistics_float_sum (float @code{*in}, size_t 
@code{size})
-Return the sum of the single precision floating point arrays. Note that NaN
-values will be ignored.
address@hidden deftypefun
-
address@hidden float gal_statistics_float_sum_num (float @code{*in}, size_t 
@code{*size})
-Return the sum of elements in the floating point array along with the
-number of used (non-Nan) elements.
address@hidden deftypefun
-
address@hidden float gal_statistics_float_sum_squared (float @code{*in}, size_t 
@code{size})
-Return the sum of squares of non-NaN elements in the array.
address@hidden deftypefun
-
address@hidden float gal_statistics_float_sum_mask (float @code{*in}, unsigned 
char @code{*mask}, size_t @code{size}, size_t @code{*nsize})
-Given a mask (@code{mask}), return the sum of elements with a @code{mask}
-value of @code{0} (masked pixels are considered to have a non-zero value)
-and put the number of used elements in @code{nsize}.
address@hidden deftypefun
-
address@hidden float gal_statistics_float_sum_mask_l (float @code{*in}, long 
@code{*mask}, size_t @code{size}, size_t @code{*nsize})
-Similar to @code{gal_statistics_float_sum_mask}, but when the mask is a
-long type.
address@hidden deftypefun
-
address@hidden float gal_statistics_float_sum_squared_mask (float @code{*in}, 
unsigned char @code{*mask}, size_t @code{size}, size_t @code{*nsize})
-Return the sum of the squared elements in @code{in} that are not masked
-(corresponding element in @code{mask} is @code{0}.
address@hidden deftypefun
-
address@hidden float gal_statistics_float_sum_squared_mask_l (float @code{*in}, 
long @code{*mask}, size_t @code{size}, size_t @code{*nsize})
-Similar to @code{gal_statistics_float_sum_squared_mask} but when the mask
-is a long type.
address@hidden deftypefun
-
address@hidden float gal_statistics_float_average (float @code{*in}, size_t 
@code{size})
-Return the average of the single precision floating point array @code{in}
-with @code{size} elements.
address@hidden deftypefun
-
address@hidden double gal_statistics_double_average (double @code{*in}, size_t 
@code{size})
-Return the average of the double precision floating point array @code{in}
-with @code{size} elements.
address@hidden deftypefun
-
address@hidden void gal_statistics_f_ave (float @code{*in}, size_t @code{size}, 
float @code{*ave}, unsigned char @code{*mask})
-Return the average of the single precision floating point array @code{in}
-with @code{size} elements, do not use those elements where @code{mask} is
-non-zero.
address@hidden deftypefun
-
address@hidden void gal_statistics_f_ave_l (float @code{*in}, size_t 
@code{size}, float @code{*ave}, long @code{*mask})
-Similar to @code{gal_statistics_f_ave}, but when the mask is of type long.
address@hidden deftypefun
-
address@hidden void gal_statistics_f_ave_std (float @code{*in}, size_t 
@code{size}, float @code{*ave}, float @code{*std}, unsigned char @code{*mask})
-Find the average and standard deviation of the @code{in} array for pixels
-where @code{mask} has value of @code{0}.
address@hidden deftypefun
-
address@hidden void gal_statistics_f_ave_std_l (float @code{*in}, size_t 
@code{size}, float @code{*ave}, float @code{*std}, long @code{*mask})
-Similar to @code{gal_statistics_f_ave_std}, but when the mask has a type of
-long.
address@hidden deftypefun
-
address@hidden float gal_statistics_median (float @code{*array}, size_t 
@code{insize})
-Return the median value of a single precision floating point array
address@hidden (which has @code{insize} elements).
address@hidden deftypefun
-
address@hidden double gal_statistics_median_double_in_place (double 
@code{*array}, size_t @code{insize})
-Return the median value of a double precision floating point array
address@hidden (which has @code{insize} elements). After this function, the
-elements of @code{array} are also sorted.
address@hidden deftypefun
-
 @deftypefun void gal_statistics_set_bins (float @code{*sorted}, size_t 
@code{size}, size_t @code{numbins}, float @code{min}, float @code{max}, float 
@code{onebinvalue}, float @code{quant}, float @code{**obins})
 Given the input @code{sorted} array with @code{size} elements, a given
 number of @code{numbins}, an optional @code{min} and @code{max}, find the
@@ -17742,7 +17446,7 @@ ApJS 220, 1. arXiv:1505.01664).
 
 
 
address@hidden Multithreaded programming, Text table into a C array, 
Statistical operations, Gnuastro library
address@hidden Multithreaded programming, World Coordinate System, Statistical 
operations, Gnuastro library
 @subsection Multithreaded programming (@file{threads.h})
 
 @cindex Multithreaded programming
@@ -17918,68 +17622,8 @@ program in @file{tests/lib/multithread.c} for a 
demonstration.
 
 
 
address@hidden Text table into a C array, World Coordinate System, 
Multithreaded programming, Gnuastro library
address@hidden Text table into a C array (@file{txtarray.h})
-
-One of the most common formats to store and transfer normally sized data
-are ASCII files (which can be seen on a text editor). During your
-processing experiments (with the data, or what programmers call
-``developing'') you will also often want to see how your data is being
-processed and one of the easiest ways to check out your arrays is to write
-them out as an ASCII file and check them. The functions here are mainly
-focused on tables that contain numbers.
-
address@hidden Macro GAL_TXTARRAY_LOG
-When an element of the input table cannot be read as a number, the column
-and row number of the element are stored along with the element into one
-row of a log file. This macro contains the name of the log-file.
address@hidden deffn
-
-
address@hidden void gal_txtarray_txt_to_array (char @code{*filename}, double 
@code{**array}, size_t @code{*s0}, size_t @code{*s1})
-Store the ASCII text table in @file{filename} into the C array @code{array}
-which has @code{s0} rows and @code{s1} columns. Any blank lines or lines
-starting with @code{#} (comments) will be ignored. The column delimiters
-(separating one column from the previous and next) are address@hidden }' 
(space),
address@hidden,}', `TAB'.
address@hidden deftypefun
-
address@hidden void gal_txtarray_printf_format (int @code{numcols}, char 
@code{**fmt}, int @code{*int_cols}, int @code{*accu_cols}, int @code{*space}, 
int @code{*prec}, char @code{forg})
-Return the @code{printf} formatting string (@code{fmt}) for @code{numcols}
-columns. Three types of output columns are assumed: integers, normal
-accuracy floating points and extra-precision floating points. Normally most
-columns will not need too much accuracy when printing (for example only 3
-decimals might be enough). However usually a very limited number of columns
-will need more precision (like RA and Dec). So to store space (each
-character is a byte) and also make the table more readable, we allow these
-two types of floating points.
-
address@hidden and @code{accu_cols} are arrays containing the column
-numbers (starting from zero) that contain integers and extra-accuracy
-columns. The ending of the column indexes in @code{int_cols} and
address@hidden is defined by a negative number. @code{space} is a three
-element array which will keep the number of characters that must be used
-for the integer, normal-accuracy and extra-accuracy columns. @code{prec} is
-a two element array which contains the number of decimals to print for the
-normal and extra accuracy columns. @code{forg} (read as `f-or-g') is the
address@hidden type for the floating point numbers: either @code{f} (as in
address@hidden, which will only print decimals) and @code{g} (as in @code{%g},
-for either printing decimals or exponentials). See
address@hidden for more.
address@hidden deftypefun
-
address@hidden void gal_txtarray_array_to_txt (double @code{*array}, size_t 
@code{s0}, size_t @code{s1}, char @code{*comments}, int @code{*int_cols}, int 
@code{*accu_cols}, int @code{*space}, int @code{*prec}, char @code{forg}, const 
char @code{*filename})
-Print the C array @code{array} (with @code{s0} rows and @code{s1} columns)
-as an ASCII text table in @code{filename}. The string pointed to by
address@hidden will be printed at the top of the table. See the
-explanation of @code{gal_txtarray_printf_format} for the other arguments.
address@hidden deftypefun
-
-
-
-
 
address@hidden World Coordinate System,  , Text table into a C array, Gnuastro 
library
address@hidden World Coordinate System,  , Multithreaded programming, Gnuastro 
library
 @subsection World Coordinate System (@file{wcs.h})
 
 The FITS world coordinate system is a mechanism to give physical values to
diff --git a/lib/Makefile.am b/lib/Makefile.am
index b4dce9f..cb587d6 100644
--- a/lib/Makefile.am
+++ b/lib/Makefile.am
@@ -27,8 +27,8 @@
 ## during the tests since the bootstrapped libraries are not installed.
 ##
 ## The SYSCONFIG_DIR is necessary in `options.c'
-AM_CPPFLAGS = -I\$(top_srcdir)/bootstrapped/lib      \
-             -DSYSCONFIG_DIR=\"$(sysconfdir)\"
+AM_CPPFLAGS = -I\$(top_srcdir)/bootstrapped/lib  \
+              -DSYSCONFIG_DIR=\"$(sysconfdir)\"
 
 
 
@@ -47,10 +47,10 @@ libgnuastro_la_LIBADD = 
$(top_builddir)/bootstrapped/lib/libgnu.la
 
 
 # Specify the library .c files
-libgnuastro_la_SOURCES = arithmetic.c arithmetic-binary.c               \
-  arithmetic-onlyint.c array.c box.c checkset.c data.c fits.c git.c     \
-  linkedlist.c mesh.c mode.c options.c polygon.c qsort.c                \
-  spatialconvolve.c statistics.c table.c threads.c timing.c txt.c wcs.c
+libgnuastro_la_SOURCES = arithmetic.c arithmetic-binary.c                  \
+  arithmetic-onlyint.c blank.c box.c checkset.c data.c fits.c git.c        \
+  linkedlist.c mesh.c mode.c options.c polygon.c qsort.c spatialconvolve.c \
+  statistics.c table.c threads.c timing.c txt.c wcs.c
 
 
 
@@ -63,7 +63,7 @@ libgnuastro_la_SOURCES = arithmetic.c arithmetic-binary.c     
          \
 # installed.
 headersdir=$(top_srcdir)/lib/gnuastro
 pkginclude_HEADERS = gnuastro/config.h $(headersdir)/arithmetic.h       \
-  $(headersdir)/array.h $(headersdir)/box.h $(headersdir)/data.h        \
+  $(headersdir)/blank.h $(headersdir)/box.h $(headersdir)/data.h        \
   $(headersdir)/fits.h $(headersdir)/git.h $(headersdir)/linkedlist.h   \
   $(headersdir)/mesh.h $(headersdir)/polygon.h $(headersdir)/qsort.h    \
   $(headersdir)/spatialconvolve.h $(headersdir)/statistics.h            \
diff --git a/lib/arithmetic.c b/lib/arithmetic.c
index 531db06..f6f1c19 100644
--- a/lib/arithmetic.c
+++ b/lib/arithmetic.c
@@ -24,9 +24,11 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 #include <errno.h>
 #include <error.h>
+#include <stdio.h>
 #include <stdlib.h>
 #include <stdarg.h>
 
+#include <gnuastro/blank.h>
 #include <gnuastro/qsort.h>
 #include <gnuastro/arithmetic.h>
 
@@ -261,22 +263,129 @@ arithmetic_check_float_input(gal_data_t *in, int 
operator, char *numstr)
    such cases, even if there are blank values, the smaller and larger
    condition checked will fail, therefore the final result will be what we
    want (to ignore the blank values). */
-#define UNIFUNC_MINVALUE {                                              \
-    int blankeq = (*b==*b && gal_data_has_blank(in)) ? 1 : 0;           \
+#define UNIFUNC_MINVALUE(IT) {                                          \
+    IT *oa=o->array;                                                    \
+    int blankeq = (*b==*b && gal_blank_present(in)) ? 1 : 0;            \
     if(blankeq)                                                         \
       do if(*ia!=*b) *oa = *ia < *oa ? *ia : *oa; while(++ia<iaf);      \
     else                                                                \
       do *oa = *ia < *oa ? *ia : *oa; while(++ia<iaf);                  \
   }
 
-#define UNIFUNC_MAXVALUE {                                              \
-    int blankeq = (*b==*b && gal_data_has_blank(in)) ? 1 : 0;           \
+#define UNIFUNC_MAXVALUE(IT) {                                          \
+    IT *oa=o->array;                                                    \
+    int blankeq = (*b==*b && gal_blank_present(in)) ? 1 : 0;            \
     if(blankeq)                                                         \
       do if(*ia!=*b) *oa = *ia > *oa ? *ia : *oa; while(++ia<iaf);      \
     else                                                                \
       do *oa = *ia > *oa ? *ia : *oa; while(++ia<iaf);                  \
   }
 
+#define UNIFUNC_NUMVALUE {                                              \
+    uint64_t *oa=o->array, num=0;                                       \
+    if( gal_blank_present(in) )                                         \
+      {                                                                 \
+        if(*b==*b) /* Is integer type. */                               \
+          do if(*ia!=*b)       ++num;               while(++ia<iaf);    \
+        else       /* Is float type with NaN blank.   */                \
+          do if(*ia==*ia)      ++num;               while(++ia<iaf);    \
+      }                                                                 \
+    else num=in->size;                                                  \
+    *oa=num;                                                            \
+  }
+
+#define UNIFUNC_SUMVALUE {                                              \
+    double sum=0.0f;                                                    \
+    float *oa=o->array;                                                 \
+    if( gal_blank_present(in) )                                         \
+      {                                                                 \
+        if(*b==*b) /* Is integer type. */                               \
+          do if(*ia!=*b)              sum += *ia;   while(++ia<iaf);    \
+        else       /* Is float type with NaN blank.   */                \
+          do if(*ia==*ia)             sum += *ia;   while(++ia<iaf);    \
+      }                                                                 \
+    else                                                                \
+      do                              sum += *ia;   while(++ia<iaf);    \
+    *oa=sum;                                                            \
+  }
+
+#define UNIFUNC_MEANVALUE {                                             \
+    int64_t num=0;                                                      \
+    double sum=0.0f;                                                    \
+    float *oa=o->array;                                                 \
+    if( gal_blank_present(in) )                                         \
+      {                                                                 \
+        if(*b==*b) /* Is integer type. */                               \
+          do if(*ia!=*b)     { ++num; sum += *ia; } while(++ia<iaf);    \
+        else       /* Is float type with NaN blank.   */                \
+          do if(*ia==*ia)    { ++num; sum += *ia; } while(++ia<iaf);    \
+      }                                                                 \
+    else                                                                \
+      {                                                                 \
+        do                            sum += *ia;   while(++ia<iaf);    \
+        num=in->size;                                                   \
+      }                                                                 \
+    *oa=sum/num;                                                        \
+  }
+
+#define UNIFUNC_STDVALUE {                                              \
+    int64_t n=0;                                                        \
+    float *oa=o->array;                                                 \
+    double s=0.0f, s2=0.0f;                                             \
+    if( gal_blank_present(in) )                                         \
+      {                                                                 \
+        if(*b==*b) /* Is integer type. */                               \
+          do if(*ia!=*b)                                                \
+               { ++n; s += *ia; s2 += *ia * *ia; } while(++ia<iaf);     \
+        else       /* Is float type with NaN blank.   */                \
+          do if(*ia==*ia)                                               \
+               { ++n; s += *ia; s2 += *ia * *ia; } while(++ia<iaf);     \
+      }                                                                 \
+    else                                                                \
+      {                                                                 \
+        do     {      s += *ia; s2 += *ia * *ia; }  while(++ia<iaf);    \
+        n=in->size;                                                     \
+      }                                                                 \
+    *oa=sqrt( (s2-s*s/n)/n );                                           \
+  }
+
+#define UNIFUNC_MEDIANVALUE(IT, QSORT_F) {                              \
+    size_t n=0;                                                         \
+    IT *a, *noblank, *oa=o->array;                                      \
+                                                                        \
+    /* Prepare space for a clean (having no blanks) dataset. If the */  \
+    /* input array is to be freed later, then just use its own space */ \
+    if( flags & GAL_ARITHMETIC_FREE ) a=noblank=in->array;              \
+    else                                                                \
+      {                                                                 \
+        errno=0;                                                        \
+        a=noblank=malloc(in->size*sizeof *noblank);                     \
+        if(noblank==NULL)                                               \
+          error(EXIT_FAILURE, 0, "%zu bytes in UNIFUNC_MEDIANVALUE",    \
+                in->size*sizeof *noblank);                              \
+      }                                                                 \
+                                                                        \
+    /* Fill it in with the elements. */                                 \
+    if(gal_blank_present(in))                                           \
+      {                                                                 \
+        if(*b==*b) do if (*ia!=*b)  { *a++=*ia; ++n;} while(++ia<iaf);  \
+        else       do if (*ia==*ia) { *a++=*ia; ++n;} while(++ia<iaf);  \
+      }                                                                 \
+    else                                                                \
+      {                                                                 \
+        n=in->size;                                                     \
+        if( (flags & GAL_ARITHMETIC_FREE)==0 )                          \
+          do *a++=*ia; while(++ia<iaf);                                 \
+      }                                                                 \
+                                                                        \
+    /* Sort the array and put the median value in the output. */        \
+    qsort(noblank, n, sizeof *noblank, QSORT_F);                        \
+    *oa = n%2 ? noblank[n/2] : (noblank[n/2] + noblank[n/2-1])/2 ;      \
+                                                                        \
+    /* Clean up. */                                                     \
+    if( (flags & GAL_ARITHMETIC_FREE)==0 ) free(noblank);               \
+  }
+
 
 
 
@@ -286,16 +395,31 @@ arithmetic_check_float_input(gal_data_t *in, int 
operator, char *numstr)
     do *oa++ = OP(*ia++); while(ia<iaf);                                \
   }
 
-#define UNIFUNC_RUN_FUNCTION_ON_ARRAY(IT){                              \
-    IT *b = gal_data_alloc_blank(in->type);                             \
-    IT *ia=in->array, *oa=o->array, *iaf=ia + in->size;                 \
+#define UNIFUNC_RUN_FUNCTION_ON_ARRAY(IT, QSORT_F){                     \
+    IT *ia=in->array, *iaf=ia + in->size;                               \
+    IT *b = gal_blank_alloc_write(in->type);                            \
     switch(operator)                                                    \
       {                                                                 \
       case GAL_ARITHMETIC_OP_MINVAL:                                    \
-        UNIFUNC_MINVALUE;                                               \
+        UNIFUNC_MINVALUE(IT);                                           \
         break;                                                          \
       case GAL_ARITHMETIC_OP_MAXVAL:                                    \
-        UNIFUNC_MAXVALUE;                                               \
+        UNIFUNC_MAXVALUE(IT);                                           \
+        break;                                                          \
+      case GAL_ARITHMETIC_OP_NUMVAL:                                    \
+        UNIFUNC_NUMVALUE;                                               \
+        break;                                                          \
+      case GAL_ARITHMETIC_OP_SUMVAL:                                    \
+        UNIFUNC_SUMVALUE;                                               \
+        break;                                                          \
+      case GAL_ARITHMETIC_OP_MEANVAL:                                   \
+        UNIFUNC_MEANVALUE;                                              \
+        break;                                                          \
+      case GAL_ARITHMETIC_OP_STDVAL:                                    \
+        UNIFUNC_STDVALUE;                                               \
+        break;                                                          \
+      case GAL_ARITHMETIC_OP_MEDIANVAL:                                 \
+        UNIFUNC_MEDIANVALUE(IT, QSORT_F);                               \
         break;                                                          \
       default:                                                          \
         error(EXIT_FAILURE, 0, "the operator code %d is not "           \
@@ -350,42 +474,42 @@ arithmetic_check_float_input(gal_data_t *in, int 
operator, char *numstr)
 
 
 
-#define UNIARY_FUNCTION_ON_ARRAY                                        \
-  switch(in->type)                                                      \
-    {                                                                   \
-    case GAL_DATA_TYPE_UINT8:                                           \
-      UNIFUNC_RUN_FUNCTION_ON_ARRAY(uint8_t)                            \
-      break;                                                            \
-    case GAL_DATA_TYPE_INT8:                                            \
-      UNIFUNC_RUN_FUNCTION_ON_ARRAY(int16_t)                            \
-      break;                                                            \
-    case GAL_DATA_TYPE_UINT16:                                          \
-      UNIFUNC_RUN_FUNCTION_ON_ARRAY(uint16_t)                           \
-      break;                                                            \
-    case GAL_DATA_TYPE_INT16:                                           \
-      UNIFUNC_RUN_FUNCTION_ON_ARRAY(int16_t)                            \
-        break;                                                          \
-    case GAL_DATA_TYPE_UINT32:                                          \
-      UNIFUNC_RUN_FUNCTION_ON_ARRAY(uint32_t)                           \
-        break;                                                          \
-    case GAL_DATA_TYPE_INT32:                                           \
-      UNIFUNC_RUN_FUNCTION_ON_ARRAY(int32_t)                            \
-        break;                                                          \
-    case GAL_DATA_TYPE_UINT64:                                          \
-      UNIFUNC_RUN_FUNCTION_ON_ARRAY(uint64_t)                           \
-        break;                                                          \
-    case GAL_DATA_TYPE_INT64:                                           \
-      UNIFUNC_RUN_FUNCTION_ON_ARRAY(int64_t)                            \
-        break;                                                          \
-    case GAL_DATA_TYPE_FLOAT32:                                         \
-      UNIFUNC_RUN_FUNCTION_ON_ARRAY(float)                              \
-      break;                                                            \
-    case GAL_DATA_TYPE_FLOAT64:                                         \
-      UNIFUNC_RUN_FUNCTION_ON_ARRAY(double)                             \
-        break;                                                          \
-    default:                                                            \
-      error(EXIT_FAILURE, 0, "type code %d not recognized in "          \
-            "`UNIFUNC_PER_ELEMENT'", in->type);                         \
+#define UNIARY_FUNCTION_ON_ARRAY                                          \
+  switch(in->type)                                                        \
+    {                                                                     \
+    case GAL_DATA_TYPE_UINT8:                                             \
+      UNIFUNC_RUN_FUNCTION_ON_ARRAY(uint8_t, gal_qsort_uint8_increasing)  \
+      break;                                                              \
+    case GAL_DATA_TYPE_INT8:                                              \
+      UNIFUNC_RUN_FUNCTION_ON_ARRAY(int8_t, gal_qsort_int8_increasing)    \
+      break;                                                              \
+    case GAL_DATA_TYPE_UINT16:                                            \
+      UNIFUNC_RUN_FUNCTION_ON_ARRAY(uint16_t, gal_qsort_uint16_increasing)\
+      break;                                                              \
+    case GAL_DATA_TYPE_INT16:                                             \
+      UNIFUNC_RUN_FUNCTION_ON_ARRAY(int16_t, gal_qsort_int16_increasing)  \
+        break;                                                            \
+    case GAL_DATA_TYPE_UINT32:                                            \
+      UNIFUNC_RUN_FUNCTION_ON_ARRAY(uint32_t, gal_qsort_uint32_increasing)\
+        break;                                                            \
+    case GAL_DATA_TYPE_INT32:                                             \
+      UNIFUNC_RUN_FUNCTION_ON_ARRAY(int32_t, gal_qsort_int32_increasing)  \
+        break;                                                            \
+    case GAL_DATA_TYPE_UINT64:                                            \
+      UNIFUNC_RUN_FUNCTION_ON_ARRAY(uint64_t, gal_qsort_uint64_increasing)\
+        break;                                                            \
+    case GAL_DATA_TYPE_INT64:                                             \
+      UNIFUNC_RUN_FUNCTION_ON_ARRAY(int64_t, gal_qsort_int64_increasing)  \
+        break;                                                            \
+    case GAL_DATA_TYPE_FLOAT32:                                           \
+      UNIFUNC_RUN_FUNCTION_ON_ARRAY(float, gal_qsort_float32_increasing)  \
+      break;                                                              \
+    case GAL_DATA_TYPE_FLOAT64:                                           \
+      UNIFUNC_RUN_FUNCTION_ON_ARRAY(double, gal_qsort_float64_increasing) \
+        break;                                                            \
+    default:                                                              \
+      error(EXIT_FAILURE, 0, "type code %d not recognized in "            \
+            "`UNIFUNC_PER_ELEMENT'", in->type);                           \
     }
 
 
@@ -405,14 +529,24 @@ arithmetic_unary_function(int operator, unsigned char 
flags, gal_data_t *in)
 
     /* Operators with only one value as output. */
     case GAL_ARITHMETIC_OP_MINVAL:
-      o = gal_data_alloc(NULL, in->type, 1, &dsize, NULL, 0, -1,
-                         NULL, NULL, NULL);
-      gal_data_type_max(o->type, o->array);
-      break;
     case GAL_ARITHMETIC_OP_MAXVAL:
+    case GAL_ARITHMETIC_OP_MEDIANVAL:
       o = gal_data_alloc(NULL, in->type, 1, &dsize, NULL, 0, -1,
                          NULL, NULL, NULL);
-      gal_data_type_min(o->type, o->array);
+      if(operator==GAL_ARITHMETIC_OP_MINVAL)       /* For the min and max  */
+        gal_data_type_max(o->type, o->array);      /* operators we need to */
+      else if (operator==GAL_ARITHMETIC_OP_MAXVAL) /* start with the max   */
+        gal_data_type_min(o->type, o->array);      /* and min values       */
+      break;                                       /* respectively.        */
+    case GAL_ARITHMETIC_OP_NUMVAL:
+      o = gal_data_alloc(NULL, GAL_DATA_TYPE_UINT64, 1, &dsize,
+                         NULL, 0, -1, NULL, NULL, NULL);
+      break;
+    case GAL_ARITHMETIC_OP_SUMVAL:
+    case GAL_ARITHMETIC_OP_MEANVAL:
+    case GAL_ARITHMETIC_OP_STDVAL:
+      o = gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT32, 1, &dsize,
+                         NULL, 0, -1, NULL, NULL, NULL);
       break;
 
     /* The other operators  */
@@ -441,6 +575,11 @@ arithmetic_unary_function(int operator, unsigned char 
flags, gal_data_t *in)
 
     case GAL_ARITHMETIC_OP_MINVAL:
     case GAL_ARITHMETIC_OP_MAXVAL:
+    case GAL_ARITHMETIC_OP_NUMVAL:
+    case GAL_ARITHMETIC_OP_SUMVAL:
+    case GAL_ARITHMETIC_OP_MEANVAL:
+    case GAL_ARITHMETIC_OP_STDVAL:
+    case GAL_ARITHMETIC_OP_MEDIANVAL:
       UNIARY_FUNCTION_ON_ARRAY;
       break;
 
@@ -659,13 +798,26 @@ arithmetic_binary_function_flt(int operator, unsigned 
char flags,
 /***********************************************************************/
 /***************                  Where                   **************/
 /***********************************************************************/
+
+/* When the `iftrue' dataset only has one element and the element is blank,
+   then it will be replaced with the blank value of the type of the output
+   data. */
 #define DO_WHERE_OPERATION(ITT, OT) {                                \
     ITT *it=iftrue->array;                                           \
-    OT *o=out->array, *of=o+out->size;                               \
+    OT *b, *o=out->array, *of=o+out->size;                           \
     if(iftrue->size==1)                                              \
-      do { *o = *c++ ? *it : *o;       } while(++o<of);              \
+      {                                                              \
+        if( gal_blank_present(iftrue) )                              \
+          {                                                          \
+            b=gal_blank_alloc_write(out->type);                      \
+            do { *o = *c++ ? *b : *o;        } while(++o<of);        \
+            free(b);                                                 \
+          }                                                          \
+        else                                                         \
+          do   { *o = *c++ ? *it : *o;       } while(++o<of);        \
+      }                                                              \
     else                                                             \
-      do { *o = *c++ ? *it : *o; ++it; } while(++o<of);              \
+      do       { *o = *c++ ? *it : *o; ++it; } while(++o<of);        \
 }
 
 
@@ -844,6 +996,32 @@ arithmetic_where(unsigned char flags, gal_data_t *out, 
gal_data_t *cond,
 
 
 
+#define MULTIOPERAND_NUM {                                              \
+    int n, use;                                                         \
+    do    /* Loop over each pixel */                                    \
+      {                                                                 \
+        n=0;                                                            \
+        use=1;                                                          \
+        for(i=0;i<dnum;++i)  /* Loop over each array. */                \
+          {                                                             \
+            /* Only integers and non-NaN floats: v==v is 1. */          \
+            if(hasblank[i])                                             \
+              use = ( *b==*b                                            \
+                      ? ( *a[i]!=*b    ? 1 : 0 )          /* Integer */ \
+                      : ( *a[i]==*a[i] ? 1 : 0 ) );       /* Float   */ \
+                                                                        \
+            /* a[i] must be incremented to next pixel in any case. */   \
+            if(use) ++n; else ++a[i];                                   \
+          }                                                             \
+        *o++ = n;                                                       \
+      }                                                                 \
+    while(o<of);                                                        \
+  }
+
+
+
+
+
 #define MULTIOPERAND_SUM {                                              \
     double sum;                                                         \
     int n, use;                                                         \
@@ -872,7 +1050,7 @@ arithmetic_where(unsigned char flags, gal_data_t *out, 
gal_data_t *cond,
 
 
 
-#define MULTIOPERAND_AVERAGE {                                          \
+#define MULTIOPERAND_MEAN {                                             \
     double sum;                                                         \
     int n, use;                                                         \
     do    /* Loop over each pixel */                                    \
@@ -900,6 +1078,40 @@ arithmetic_where(unsigned char flags, gal_data_t *out, 
gal_data_t *cond,
 
 
 
+#define MULTIOPERAND_STD {                                              \
+    int n, use;                                                         \
+    double sum, sum2;                                                   \
+    do    /* Loop over each pixel */                                    \
+      {                                                                 \
+        n=0;                                                            \
+        use=1;                                                          \
+        sum=sum2=0.0f;                                                  \
+        for(i=0;i<dnum;++i)  /* Loop over each array. */                \
+          {                                                             \
+            /* Only integers and non-NaN floats: v==v is 1. */          \
+            if(hasblank[i])                                             \
+              use = ( *b==*b                                            \
+                      ? ( *a[i]!=*b    ? 1 : 0 )          /* Integer */ \
+                      : ( *a[i]==*a[i] ? 1 : 0 ) );       /* Float   */ \
+                                                                        \
+            /* a[i] must be incremented to next pixel in any case. */   \
+            if(use)                                                     \
+              {                                                         \
+                sum2  = *a[i] * *a[i];                                  \
+                sum  += *a[i]++;                                        \
+                ++n;                                                    \
+              }                                                         \
+            else ++a[i];                                                \
+          }                                                             \
+        *o++ = n ? sqrt( (sum2-sum*sum/n)/n ) : *b;                     \
+      }                                                                 \
+    while(o<of);                                                        \
+  }
+
+
+
+
+
 #define MULTIOPERAND_MEDIAN(TYPE, QSORT_F) {                            \
     TYPE *pixs;                                                         \
     int n, use;                                                         \
@@ -944,7 +1156,7 @@ arithmetic_where(unsigned char flags, gal_data_t *out, 
gal_data_t *cond,
 
 #define MULTIOPERAND_TYPE_SET(TYPE, QSORT_F) {                          \
     TYPE *o=out->array, *of=out->array+out->size;                       \
-    TYPE **a, *b=gal_data_alloc_blank(list->type);                      \
+    TYPE **a, *b=gal_blank_alloc_write(list->type);                     \
                                                                         \
     /* Allocate space to keep the pointers to the arrays of each. */    \
     /* Input data structure. The operators will increment these */      \
@@ -970,12 +1182,20 @@ arithmetic_where(unsigned char flags, gal_data_t *out, 
gal_data_t *cond,
         MULTIOPERAND_MAX(TYPE);                                         \
         break;                                                          \
                                                                         \
+      case GAL_ARITHMETIC_OP_NUM:                                       \
+        MULTIOPERAND_NUM;                                               \
+        break;                                                          \
+                                                                        \
       case GAL_ARITHMETIC_OP_SUM:                                       \
         MULTIOPERAND_SUM;                                               \
         break;                                                          \
                                                                         \
-      case GAL_ARITHMETIC_OP_AVERAGE:                                   \
-        MULTIOPERAND_AVERAGE;                                           \
+      case GAL_ARITHMETIC_OP_MEAN:                                      \
+        MULTIOPERAND_MEAN;                                              \
+        break;                                                          \
+                                                                        \
+      case GAL_ARITHMETIC_OP_STD:                                       \
+        MULTIOPERAND_STD;                                               \
         break;                                                          \
                                                                         \
       case GAL_ARITHMETIC_OP_MEDIAN:                                    \
@@ -1051,41 +1271,41 @@ arithmetic_multioperand(int operator, unsigned char 
flags, gal_data_t *list)
 
   /* Fill in the hasblank array. */
   for(tmp=list;tmp!=NULL;tmp=tmp->next)
-    hasblank[i++]=gal_data_has_blank(tmp);
+    hasblank[i++]=gal_blank_present(tmp);
 
 
   /* Start the operation. */
   switch(list->type)
     {
     case GAL_DATA_TYPE_UINT8:
-      MULTIOPERAND_TYPE_SET(uint8_t,   gal_qsort_uchar_increasing);
+      MULTIOPERAND_TYPE_SET(uint8_t,   gal_qsort_uint8_increasing);
       break;
     case GAL_DATA_TYPE_INT8:
-      MULTIOPERAND_TYPE_SET(int8_t,    gal_qsort_char_increasing);
+      MULTIOPERAND_TYPE_SET(int8_t,    gal_qsort_int8_increasing);
       break;
     case GAL_DATA_TYPE_UINT16:
-      MULTIOPERAND_TYPE_SET(uint16_t,  gal_qsort_ushort_increasing);
+      MULTIOPERAND_TYPE_SET(uint16_t,  gal_qsort_uint16_increasing);
       break;
     case GAL_DATA_TYPE_INT16:
-      MULTIOPERAND_TYPE_SET(int16_t,   gal_qsort_short_increasing);
+      MULTIOPERAND_TYPE_SET(int16_t,   gal_qsort_int16_increasing);
       break;
     case GAL_DATA_TYPE_UINT32:
-      MULTIOPERAND_TYPE_SET(uint32_t,  gal_qsort_uint_increasing);
+      MULTIOPERAND_TYPE_SET(uint32_t,  gal_qsort_uint32_increasing);
       break;
     case GAL_DATA_TYPE_INT32:
-      MULTIOPERAND_TYPE_SET(int32_t,   gal_qsort_int_increasing);
+      MULTIOPERAND_TYPE_SET(int32_t,   gal_qsort_int32_increasing);
       break;
     case GAL_DATA_TYPE_UINT64:
-      MULTIOPERAND_TYPE_SET(uint64_t,  gal_qsort_ulong_increasing);
+      MULTIOPERAND_TYPE_SET(uint64_t,  gal_qsort_uint64_increasing);
       break;
     case GAL_DATA_TYPE_INT64:
-      MULTIOPERAND_TYPE_SET(int64_t,   gal_qsort_long_increasing);
+      MULTIOPERAND_TYPE_SET(int64_t,   gal_qsort_int64_increasing);
       break;
     case GAL_DATA_TYPE_FLOAT32:
-      MULTIOPERAND_TYPE_SET(float,     gal_qsort_float_increasing);
+      MULTIOPERAND_TYPE_SET(float,     gal_qsort_float32_increasing);
       break;
     case GAL_DATA_TYPE_FLOAT64:
-      MULTIOPERAND_TYPE_SET(double,    gal_qsort_double_increasing);
+      MULTIOPERAND_TYPE_SET(double,    gal_qsort_float64_increasing);
       break;
     default:
       error(EXIT_FAILURE, 0, "type code %d not recognized in "
@@ -1328,9 +1548,17 @@ gal_arithmetic_operator_string(int operator)
 
     case GAL_ARITHMETIC_OP_MINVAL:       return "minvalue";
     case GAL_ARITHMETIC_OP_MAXVAL:       return "maxvalue";
+    case GAL_ARITHMETIC_OP_NUMVAL:       return "numvalue";
+    case GAL_ARITHMETIC_OP_SUMVAL:       return "sumvalue";
+    case GAL_ARITHMETIC_OP_MEANVAL:      return "meanvalue";
+    case GAL_ARITHMETIC_OP_STDVAL:       return "stdvalue";
+    case GAL_ARITHMETIC_OP_MEDIANVAL:    return "medianvalue";
     case GAL_ARITHMETIC_OP_MIN:          return "min";
     case GAL_ARITHMETIC_OP_MAX:          return "max";
-    case GAL_ARITHMETIC_OP_AVERAGE:      return "average";
+    case GAL_ARITHMETIC_OP_NUM:          return "num";
+    case GAL_ARITHMETIC_OP_SUM:          return "sum";
+    case GAL_ARITHMETIC_OP_MEAN:         return "mean";
+    case GAL_ARITHMETIC_OP_STD:          return "std";
     case GAL_ARITHMETIC_OP_MEDIAN:       return "median";
 
     case GAL_ARITHMETIC_OP_TO_UINT8:     return "uchar";
@@ -1452,7 +1680,7 @@ gal_arithmetic(int operator, unsigned char flags, ...)
 
     case GAL_ARITHMETIC_OP_ISBLANK:
       d1 = va_arg(va, gal_data_t *);
-      out = gal_data_flag_blank(d1);
+      out = gal_blank_flag(d1);
       if(flags & GAL_ARITHMETIC_FREE) gal_data_free(d1);
       break;
 
@@ -1471,6 +1699,11 @@ gal_arithmetic(int operator, unsigned char flags, ...)
     case GAL_ARITHMETIC_OP_LOG10:
     case GAL_ARITHMETIC_OP_MINVAL:
     case GAL_ARITHMETIC_OP_MAXVAL:
+    case GAL_ARITHMETIC_OP_NUMVAL:
+    case GAL_ARITHMETIC_OP_SUMVAL:
+    case GAL_ARITHMETIC_OP_MEANVAL:
+    case GAL_ARITHMETIC_OP_STDVAL:
+    case GAL_ARITHMETIC_OP_MEDIANVAL:
       d1 = va_arg(va, gal_data_t *);
       out=arithmetic_unary_function(operator, flags, d1);
       break;
@@ -1482,8 +1715,10 @@ gal_arithmetic(int operator, unsigned char flags, ...)
 
     case GAL_ARITHMETIC_OP_MIN:
     case GAL_ARITHMETIC_OP_MAX:
+    case GAL_ARITHMETIC_OP_NUM:
     case GAL_ARITHMETIC_OP_SUM:
-    case GAL_ARITHMETIC_OP_AVERAGE:
+    case GAL_ARITHMETIC_OP_MEAN:
+    case GAL_ARITHMETIC_OP_STD:
     case GAL_ARITHMETIC_OP_MEDIAN:
       d1 = va_arg(va, gal_data_t *);
       out=arithmetic_multioperand(operator, flags, d1);
diff --git a/lib/array.c b/lib/array.c
deleted file mode 100644
index 301286a..0000000
--- a/lib/array.c
+++ /dev/null
@@ -1,530 +0,0 @@
-/*********************************************************************
-Functions to manipulate arrays.
-This is part of GNU Astronomy Utilities (Gnuastro) package.
-
-Original author:
-     Mohammad Akhlaghi <address@hidden>
-Contributing author(s):
-Copyright (C) 2015, 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 <math.h>
-#include <errno.h>
-#include <error.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <pthread.h>
-
-#include <gnuastro/array.h>
-
-
-
-
-/*********************************************************************
- **********************       Initialize        **********************
- *********************************************************************/
-void
-gal_array_uchar_init_on_region(unsigned char *in, const unsigned char v,
-                               size_t start, size_t s0, size_t s1,
-                               size_t is1)
-{
-  size_t r;
-  unsigned char *p, *fp;
-
-  for(r=0;r<s0;++r)
-    {
-      fp = (p=in+start) + s1;
-      do
-        *p=v;
-      while(++p<fp);
-      start+=is1;
-    }
-}
-
-
-
-
-void
-gal_array_long_init(long *in, size_t size, const long v)
-{
-  long *end=in+size;
-  do *in++=v; while(in<end);
-}
-
-
-
-
-
-void
-gal_array_long_init_on_region(long *in, const long v, size_t start,
-                              size_t s0, size_t s1, size_t is1)
-{
-  size_t r;
-  long *p, *fp;
-
-  for(r=0;r<s0;++r)
-    {
-      fp = (p=in+start) + s1;
-      do
-        *p=v;
-      while(++p<fp);
-      start+=is1;
-    }
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-/*********************************************************************
- **********************       Copy array        **********************
- *********************************************************************/
-void
-gal_array_uchar_copy(unsigned char *in, size_t size, unsigned char **out)
-{
-  unsigned char *fp=in+size, *o;
-
-  errno=0;
-  o=*out=malloc(size*sizeof *out);
-  if(*out==NULL)
-    error(EXIT_FAILURE, errno, "%zu bytes for copying", size);
-  do *o++=*in; while(++in<fp);
-}
-
-
-
-
-
-void
-gal_array_float_copy(float *in, size_t size, float **out)
-{
-  float *fp=in+size, *o;
-
-  errno=0;
-  o=*out=malloc(size*sizeof *out);
-  if(*out==NULL)
-    error(EXIT_FAILURE, errno, "%zu bytes for copying", size);
-  do *o++=*in; while(++in<fp);
-}
-
-
-
-
-
-void
-gal_array_float_copy_noalloc(float *in, size_t size, float *out)
-{
-  float *fp=in+size;
-  do *out++=*in; while(++in<fp);
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-/*********************************************************************
- **********************         Values          **********************
- *********************************************************************/
-void
-gal_array_fset_const(float *in, size_t size, float a)
-{
-  float *fpt;
-  fpt=in+size;
-  do
-    *in = a;
-  while(++in<fpt);
-}
-
-
-
-
-
-void
-gal_array_freplace_value(float *in, size_t size, float from, float to)
-{
-  float *fpt;
-  fpt=in+size;
-  if(isnan(from))
-    {
-      do
-        *in = isnan(*in) ? to : *in;
-      while(++in<fpt);
-    }
-  else
-    {
-      do
-        *in = *in==from ? to : *in;
-      while(++in<fpt);
-    }
-}
-
-
-
-
-/* Only replace non-NaN values in a float array. */
-void
-gal_array_freplace_nonnans(float *in, size_t size, float to)
-{
-  float *fpt=in+size;
-  do
-    *in = isnan(*in) ? *in : to;
-  while(++in<fpt);
-}
-
-
-
-
-
-/* Move all the non-NaN elements in the array to the start of the
-   array so that the non-NaN alements are contiguous. This is useful
-   for cases where you want to sort the data.*/
-void
-gal_array_no_nans(float *in, size_t *size)
-{
-  size_t outsize=0;
-  float *f=in, *fp=in+*size;
-  do
-    if(!isnan(*in))
-      {
-        *f++=*in;
-        ++outsize;
-      }
-  while(++in<fp);
-  *size=outsize;
-}
-
-
-
-
-
-void
-gal_array_no_nans_double(double *in, size_t *size)
-{
-  size_t outsize=0;
-  double *f=in, *fp=in+*size;
-  do
-    if(!isnan(*in))
-      {
-        *f++=*in;
-        ++outsize;
-      }
-  while(++in<fp);
-  *size=outsize;
-}
-
-
-
-
-
-void
-gal_array_uchar_replace(unsigned char *in, size_t size,
-                        unsigned char from, unsigned char to)
-{
-  unsigned char *fpt=in+size;
-  do
-    *in = *in==from ? to : *in;
-  while(++in<fpt);
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-/*********************************************************************
- **********************   Multiply or Sum with  **********************
- *********************************************************************/
-void
-gal_array_fmultip_const(float *in, size_t size, float a)
-{
-  float *fpt;
-  fpt=in+size;
-  do
-    *in *= a;
-  while(++in<fpt);
-}
-
-
-
-
-
-void
-gal_array_fsum_const(float *in, size_t size, float a)
-{
-  float *fpt;
-  fpt=in+size;
-  do
-    *in += a;
-  while(++in<fpt);
-}
-
-
-
-
-
-float *
-gal_array_fsum_arrays(float *in1, float *in2, size_t size)
-{
-  float *out, *o, *op;
-
-  errno=0;
-  o=out=malloc(size*sizeof *out);
-  if(out==NULL)
-    error(EXIT_FAILURE, errno, "%zu bytes for out in gal_array_fsum_arrays"
-          " (array.c)", size*sizeof *out);
-
-  op=o+size;
-  do *o = *in1++ + *in2++; while(++o<op);
-
-  return out;
-}
-
-
-
-
-
-void
-gal_array_dmultip_const(double *in, size_t size, double a)
-{
-  double *dpt;
-  dpt=in+size;
-  do *in *= a; while(++in<dpt);
-}
-
-
-
-
-void
-gal_array_dmultip_arrays(double *in1, double *in2, size_t size)
-{
-  double *dpt;
-  dpt=in1+size;
-  do *in1 *= *in2++; while(++in1<dpt);
-}
-
-
-
-
-
-void
-gal_array_ddivide_const(double *in, size_t size, double a)
-{
-  double *dpt;
-  dpt=in+size;
-  do *in /= a; while(++in<dpt);
-}
-
-
-
-
-
-void
-gal_array_dconst_divide(double *in, size_t size, double a)
-{
-  double *dpt;
-  dpt=in+size;
-  do *in = a / *in; while(++in<dpt);
-}
-
-
-
-
-void
-gal_array_ddivide_arrays(double *in1, double *in2, size_t size)
-{
-  double *dpt;
-  dpt=in1+size;
-  do *in1 /= *in2++; while(++in1<dpt);
-}
-
-
-
-
-void
-gal_array_dsum_const(double *in, size_t size, double a)
-{
-  double *dpt;
-  dpt=in+size;
-  do *in += a; while(++in<dpt);
-}
-
-
-
-
-
-void
-gal_array_dsum_arrays(double *in1, double *in2, size_t size)
-{
-  double *dpt;
-  dpt=in1+size;
-  do *in1 += *in2++; while(++in1<dpt);
-}
-
-
-
-void
-gal_array_dsubtract_const(double *in, size_t size, double a)
-{
-  double *dpt;
-  dpt=in+size;
-  do *in -= a; while(++in<dpt);
-}
-
-
-
-
-
-void
-gal_array_dconst_subtract(double *in, size_t size, double a)
-{
-  double *dpt;
-  dpt=in+size;
-  do *in = a - *in; while(++in<dpt);
-}
-
-
-
-
-
-void
-gal_array_dsubtract_arrays(double *in1, double *in2, size_t size)
-{
-  double *dpt;
-  dpt=in1+size;
-  do *in1 -= *in2++; while(++in1<dpt);
-}
-
-
-
-void
-gal_array_dpower_const(double *in, size_t size, double a)
-{
-  double *dpt;
-  dpt=in+size;
-
-  /* Since simply multiplying or taking the square root are more
-     efficient than calling the power function. */
-  if(a==2.0f)      do *in *= *in; while(++in<dpt);
-  else if(a==0.5f) do *in = sqrt(*in); while(++in<dpt);
-  else             do *in = pow(*in, a); while(++in<dpt);
-}
-
-
-
-
-
-void
-gal_array_dconst_power(double *in, size_t size, double a)
-{
-  double *dpt;
-  dpt=in+size;
-  do *in = pow(a, *in); while(++in<dpt);
-}
-
-
-
-
-
-void
-gal_array_dpower_arrays(double *in1, double *in2, size_t size)
-{
-  double *dpt;
-  dpt=in1+size;
-  do *in1 = pow(*in1, *in2++); while(++in1<dpt);
-}
-
-
-
-
-
-void
-gal_array_dlog_array(double *in, size_t size)
-{
-  double *dpt;
-  dpt=in+size;
-
-  do *in = log(*in); while(++in<dpt);
-}
-
-
-
-
-
-void
-gal_array_dlog10_array(double *in, size_t size)
-{
-  double *dpt;
-  dpt=in+size;
-
-  do *in = log10(*in); while(++in<dpt);
-}
-
-
-
-
-
-void
-gal_array_dabs_array(double *in, size_t size)
-{
-  double *dpt;
-  dpt=in+size;
-
-  do *in = fabs(*in); while(++in<dpt);
-}
diff --git a/lib/blank.c b/lib/blank.c
new file mode 100644
index 0000000..4474a09
--- /dev/null
+++ b/lib/blank.c
@@ -0,0 +1,402 @@
+/*********************************************************************
+blank -- Deal with blank values in a datasets
+This is part of GNU Astronomy Utilities (Gnuastro) package.
+
+Original author:
+     Mohammad Akhlaghi <address@hidden>
+Contributing author(s):
+Copyright (C) 2017, 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 <errno.h>
+#include <error.h>
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+
+#include <gnuastro/data.h>
+#include <gnuastro/blank.h>
+
+#include <checkset.h>
+
+
+
+
+/* Write the blank value of the type into an already allocate space.
+
+   Note that for strings, pointer should actually be `char **'. */
+void
+gal_blank_write(void *pointer, uint8_t type)
+{
+  switch(type)
+    {
+    case GAL_DATA_TYPE_STRING:
+      gal_checkset_allocate_copy(GAL_BLANK_STRING, pointer);
+      break;
+
+    case GAL_DATA_TYPE_UINT8:
+      *(uint8_t *)pointer       = GAL_BLANK_UINT8;
+      break;
+
+    case GAL_DATA_TYPE_INT8:
+      *(int8_t *)pointer        = GAL_BLANK_INT8;
+      break;
+
+    case GAL_DATA_TYPE_UINT16:
+      *(uint16_t *)pointer      = GAL_BLANK_UINT16;
+      break;
+
+    case GAL_DATA_TYPE_INT16:
+      *(int16_t *)pointer       = GAL_BLANK_INT16;
+      break;
+
+    case GAL_DATA_TYPE_UINT32:
+      *(uint32_t *)pointer      = GAL_BLANK_UINT32;
+      break;
+
+    case GAL_DATA_TYPE_INT32:
+      *(int32_t *)pointer       = GAL_BLANK_INT32;
+      break;
+
+    case GAL_DATA_TYPE_UINT64:
+      *(uint64_t *)pointer      = GAL_BLANK_UINT64;
+      break;
+
+    case GAL_DATA_TYPE_INT64:
+      *(int64_t *)pointer       = GAL_BLANK_INT64;
+      break;
+
+    case GAL_DATA_TYPE_FLOAT32:
+      *(float *)pointer         = GAL_BLANK_FLOAT32;
+      break;
+
+    case GAL_DATA_TYPE_FLOAT64:
+      *(double *)pointer        = GAL_BLANK_FLOAT64;
+      break;
+
+    default:
+      error(EXIT_FAILURE, 0, "type code %d not recognized in "
+            "`gal_blank_write'", type);
+    }
+}
+
+
+
+
+
+/* Allocate some space for the given type and put the blank value into
+   it. */
+void *
+gal_blank_alloc_write(uint8_t type)
+{
+  void *out;
+
+  /* Allocate the space to keep the blank value. */
+  out=gal_data_malloc_array(type, 1);
+
+  /* Put the blank value in the allcated space. */
+  gal_blank_write(out, type);
+
+  /* Return the allocated space. */
+  return out;
+}
+
+
+
+
+
+/* Print the blank value as a string. */
+char *
+gal_blank_as_string(uint8_t type, int width)
+{
+  char *blank;
+  switch(type)
+    {
+    case GAL_DATA_TYPE_BIT:
+      error(EXIT_FAILURE, 0, "bit types are not implemented in "
+            "`gal_data_blank_as_string' yet.");
+      break;
+
+    case GAL_DATA_TYPE_STRING:
+      if(width)
+        asprintf(&blank, "%*s", width, GAL_BLANK_STRING);
+      else
+        asprintf(&blank, "%s", GAL_BLANK_STRING);
+      break;
+
+    case GAL_DATA_TYPE_UINT8:
+      if(width)
+        asprintf(&blank, "%*u", width, (uint8_t)GAL_BLANK_UINT8);
+      else
+        asprintf(&blank, "%u",         (uint8_t)GAL_BLANK_UINT8);
+      break;
+
+    case GAL_DATA_TYPE_INT8:
+      if(width)
+        asprintf(&blank, "%*d", width, (int8_t)GAL_BLANK_INT8);
+      else
+        asprintf(&blank, "%d",         (int8_t)GAL_BLANK_INT8);
+      break;
+
+    case GAL_DATA_TYPE_UINT16:
+      if(width)
+        asprintf(&blank, "%*u", width, (uint16_t)GAL_BLANK_UINT16);
+      else
+        asprintf(&blank, "%u",         (uint16_t)GAL_BLANK_UINT16);
+      break;
+
+    case GAL_DATA_TYPE_INT16:
+      if(width)
+        asprintf(&blank, "%*d", width, (int16_t)GAL_BLANK_INT16);
+      else
+        asprintf(&blank, "%d",         (int16_t)GAL_BLANK_INT16);
+      break;
+
+    case GAL_DATA_TYPE_UINT32:
+      if(width)
+        asprintf(&blank, "%*u", width, (uint32_t)GAL_BLANK_UINT32);
+      else
+        asprintf(&blank, "%u",         (uint32_t)GAL_BLANK_UINT32);
+      break;
+
+    case GAL_DATA_TYPE_INT32:
+      if(width)
+        asprintf(&blank, "%*d", width, (int32_t)GAL_BLANK_INT32);
+      else
+        asprintf(&blank, "%d",         (int32_t)GAL_BLANK_INT32);
+      break;
+
+    case GAL_DATA_TYPE_UINT64:
+      if(width)
+        asprintf(&blank, "%*lu", width, (uint64_t)GAL_BLANK_UINT64);
+      else
+        asprintf(&blank, "%lu",         (uint64_t)GAL_BLANK_UINT64);
+      break;
+
+    case GAL_DATA_TYPE_INT64:
+      if(width)
+        asprintf(&blank, "%*ld", width, (int64_t)GAL_BLANK_INT64);
+      else
+        asprintf(&blank, "%ld",         (int64_t)GAL_BLANK_INT64);
+      break;
+
+    case GAL_DATA_TYPE_FLOAT32:
+      if(width)
+        asprintf(&blank, "%*f", width, (float)GAL_BLANK_FLOAT32);
+      else
+        asprintf(&blank, "%f",         (float)GAL_BLANK_FLOAT32);
+      break;
+
+    case GAL_DATA_TYPE_FLOAT64:
+      if(width)
+        asprintf(&blank, "%*f", width, (double)GAL_BLANK_FLOAT64);
+      else
+        asprintf(&blank, "%f",         (double)GAL_BLANK_FLOAT64);
+      break;
+
+    default:
+      error(EXIT_FAILURE, 0, "type code %d not recognized in "
+            "`gal_blank_as_string'", type);
+    }
+  return blank;
+}
+
+
+
+
+
+
+/* For integers a simple equality is enough. */
+#define HAS_BLANK_INT(CTYPE, BLANK) {                                   \
+    CTYPE *a=data->array, *af=a+data->size;                             \
+    do if(*a++ == BLANK) return 1; while(a<af);                         \
+  }
+
+/* Note that a NaN value is not equal to another NaN value, so we can't use
+   the easy check for cases were the blank value is NaN. Also note that
+   `isnan' is actually a macro, so it works for both float and double
+   types.*/
+#define HAS_BLANK_FLT(CTYPE, BLANK, MULTIP) {                           \
+    CTYPE *a=data->array, *af=a+(MULTIP*data->size);                    \
+    if(isnan(BLANK)) do if(isnan(*a++)) return 1; while(a<af);          \
+    else             do if(*a++==BLANK) return 1; while(a<af);          \
+  }
+
+/* Return 1 if the dataset has a blank value and zero if it doesn't. */
+int
+gal_blank_present(gal_data_t *data)
+{
+  char **str=data->array, **strf=str+data->size;
+
+  /* Go over the pixels and check: */
+  switch(data->type)
+    {
+    case GAL_DATA_TYPE_BIT:
+      error(EXIT_FAILURE, 0, "Currently Gnuastro doesn't support bit "
+            "datatype, please get in touch with us to implement it.");
+
+    case GAL_DATA_TYPE_UINT8:
+      HAS_BLANK_INT(uint8_t, GAL_BLANK_UINT8);     break;
+
+    case GAL_DATA_TYPE_INT8:
+      HAS_BLANK_INT(int8_t, GAL_BLANK_INT8);       break;
+
+    case GAL_DATA_TYPE_UINT16:
+      HAS_BLANK_INT(uint16_t, GAL_BLANK_UINT16);   break;
+
+    case GAL_DATA_TYPE_INT16:
+      HAS_BLANK_INT(int16_t, GAL_BLANK_INT16);     break;
+
+    case GAL_DATA_TYPE_UINT32:
+      HAS_BLANK_INT(uint32_t, GAL_BLANK_UINT32);   break;
+
+    case GAL_DATA_TYPE_INT32:
+      HAS_BLANK_INT(int32_t, GAL_BLANK_INT32);     break;
+
+    case GAL_DATA_TYPE_UINT64:
+      HAS_BLANK_INT(uint64_t, GAL_BLANK_UINT64);   break;
+
+    case GAL_DATA_TYPE_INT64:
+      HAS_BLANK_INT(int64_t, GAL_BLANK_INT64);     break;
+
+    case GAL_DATA_TYPE_FLOAT32:
+      HAS_BLANK_FLT(float, GAL_BLANK_FLOAT32, 1);  break;
+
+    case GAL_DATA_TYPE_FLOAT64:
+      HAS_BLANK_FLT(double, GAL_BLANK_FLOAT64, 1); break;
+
+    case GAL_DATA_TYPE_COMPLEX32:
+      HAS_BLANK_FLT(float, GAL_BLANK_FLOAT32, 2);  break;
+
+    case GAL_DATA_TYPE_COMPLEX64:
+      HAS_BLANK_FLT(double, GAL_BLANK_FLOAT64, 2); break;
+
+    case GAL_DATA_TYPE_STRING:
+      do if(!strcmp(*str++,GAL_BLANK_STRING)) return 1; while(str<strf);
+      break;
+
+    default:
+      error(EXIT_FAILURE, 0, "a bug! type value (%d) not recognized "
+            "in `gal_data_has_blank'", data->type);
+    }
+
+  /* If there was a blank value, then the function would have returned with
+     a value of 1. So if it reaches here, then we can be sure that there
+     was no blank values, hence, return 0. */
+  return 0;
+}
+
+
+
+
+
+/* For integers a simple equality is enough. */
+#define FLAG_BLANK_INT(CTYPE, BLANK) {                                  \
+    CTYPE *a=data->array; do *o = (*a==BLANK); while(++o<of);           \
+  }
+
+/* Note that a NaN value is not equal to another NaN value, so we can't use
+   the easy check for cases were the blank value is NaN. Also note that
+   `isnan' is actually a macro, so it works for both float and double
+   types.*/
+#define FLAG_BLANK_FLT(CTYPE, BLANK) {                                  \
+    CTYPE *a=data->array;                                               \
+    if(isnan(BLANK)) do *o = isnan(*a++);   while(++o<of);              \
+    else             do *o = (*a++==BLANK); while(++o<of);              \
+  }
+
+#define FLAG_BLANK_COMPLEX(CTYPE, BLANK) {                              \
+    CTYPE *a=data->array;                                               \
+    if(isnan(BLANK))                                                    \
+      do { *o=(isnan(*a) || isnan(*(a+1))); a+=2; } while(++o<of);      \
+    else                                                                \
+      do { *o=(*a==BLANK || *(a+1)==BLANK); a+=2; } while(++o<of);      \
+  }
+
+/* Output a data-set of the the same size as the input, but with an uint8_t
+   type that has a value of 1 for data that are blank and 0 for those that
+   aren't. */
+gal_data_t *
+gal_blank_flag(gal_data_t *data)
+{
+  uint8_t *o, *of;
+  gal_data_t *out;
+  char **str=data->array, **strf=str+data->size;
+
+  /* Allocate the output array. */
+  out=gal_data_alloc(NULL, GAL_DATA_TYPE_UINT8, data->ndim, data->dsize,
+                     data->wcs, 0, data->minmapsize, data->name, data->unit,
+                     data->comment);
+
+  /* Set the pointers for easy looping. */
+  of=(o=out->array)+data->size;
+
+  /* Go over the pixels and set the output values. */
+  switch(data->type)
+    {
+    case GAL_DATA_TYPE_BIT:
+      error(EXIT_FAILURE, 0, "Currently Gnuastro doesn't support bit "
+            "datatype, please get in touch with us to implement it.");
+
+    case GAL_DATA_TYPE_UINT8:
+      FLAG_BLANK_INT(uint8_t, GAL_BLANK_UINT8);      break;
+
+    case GAL_DATA_TYPE_INT8:
+      FLAG_BLANK_INT(int8_t, GAL_BLANK_INT8);        break;
+
+    case GAL_DATA_TYPE_UINT16:
+      FLAG_BLANK_INT(uint16_t, GAL_BLANK_UINT16);    break;
+
+    case GAL_DATA_TYPE_INT16:
+      FLAG_BLANK_INT(int16_t, GAL_BLANK_INT16);      break;
+
+    case GAL_DATA_TYPE_UINT32:
+      FLAG_BLANK_INT(uint32_t, GAL_BLANK_UINT32);    break;
+
+    case GAL_DATA_TYPE_INT32:
+      FLAG_BLANK_INT(int32_t, GAL_BLANK_INT32);      break;
+
+    case GAL_DATA_TYPE_UINT64:
+      FLAG_BLANK_INT(uint64_t, GAL_BLANK_UINT64);    break;
+
+    case GAL_DATA_TYPE_INT64:
+      FLAG_BLANK_INT(int64_t, GAL_BLANK_INT64);      break;
+
+    case GAL_DATA_TYPE_FLOAT32:
+      FLAG_BLANK_FLT(float, GAL_BLANK_FLOAT32);      break;
+
+    case GAL_DATA_TYPE_FLOAT64:
+      FLAG_BLANK_FLT(double, GAL_BLANK_FLOAT64);     break;
+
+    case GAL_DATA_TYPE_COMPLEX32:
+      FLAG_BLANK_COMPLEX(float, GAL_BLANK_FLOAT32);  break;
+
+    case GAL_DATA_TYPE_COMPLEX64:
+      FLAG_BLANK_COMPLEX(double, GAL_BLANK_FLOAT64); break;
+
+    case GAL_DATA_TYPE_STRING:
+      do *o++ = !strcmp(*str,GAL_BLANK_STRING); while(++str<strf);
+      break;
+
+    default:
+      error(EXIT_FAILURE, 0, "type value (%d) not recognized "
+            "in `gal_blank_flag'", data->type);
+    }
+
+  /* Return */
+  return out;
+}
diff --git a/lib/data.c b/lib/data.c
index 9087004..7fd3c9f 100644
--- a/lib/data.c
+++ b/lib/data.c
@@ -34,6 +34,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <sys/mman.h>
 
 #include <gnuastro/data.h>
+#include <gnuastro/blank.h>
 #include <gnuastro/table.h>
 #include <gnuastro/linkedlist.h>
 
@@ -1051,391 +1052,6 @@ gal_data_free_ll(gal_data_t *list)
 
 
 
-/*************************************************************
- **************          Blank data            ***************
- *************************************************************/
-/* Write the blank value of the type into an already allocate space.
-
-   Note that for strings, pointer should actually be `char **'. */
-void
-gal_data_set_blank(void *pointer, uint8_t type)
-{
-  switch(type)
-    {
-    case GAL_DATA_TYPE_STRING:
-      gal_checkset_allocate_copy(GAL_DATA_BLANK_STRING, pointer);
-      break;
-
-    case GAL_DATA_TYPE_UINT8:
-      *(uint8_t *)pointer       = GAL_DATA_BLANK_UINT8;
-      break;
-
-    case GAL_DATA_TYPE_INT8:
-      *(int8_t *)pointer        = GAL_DATA_BLANK_INT8;
-      break;
-
-    case GAL_DATA_TYPE_UINT16:
-      *(uint16_t *)pointer      = GAL_DATA_BLANK_UINT16;
-      break;
-
-    case GAL_DATA_TYPE_INT16:
-      *(int16_t *)pointer       = GAL_DATA_BLANK_INT16;
-      break;
-
-    case GAL_DATA_TYPE_UINT32:
-      *(uint32_t *)pointer      = GAL_DATA_BLANK_UINT32;
-      break;
-
-    case GAL_DATA_TYPE_INT32:
-      *(int32_t *)pointer       = GAL_DATA_BLANK_INT32;
-      break;
-
-    case GAL_DATA_TYPE_UINT64:
-      *(uint64_t *)pointer      = GAL_DATA_BLANK_UINT64;
-      break;
-
-    case GAL_DATA_TYPE_INT64:
-      *(int64_t *)pointer       = GAL_DATA_BLANK_INT64;
-      break;
-
-    case GAL_DATA_TYPE_FLOAT32:
-      *(float *)pointer         = GAL_DATA_BLANK_FLOAT32;
-      break;
-
-    case GAL_DATA_TYPE_FLOAT64:
-      *(double *)pointer        = GAL_DATA_BLANK_FLOAT64;
-      break;
-
-    default:
-      error(EXIT_FAILURE, 0, "type code %d not recognized in "
-            "`gal_data_set_blank'", type);
-    }
-}
-
-
-
-
-
-/* Allocate some space for the given type and put the blank value into
-   it. */
-void *
-gal_data_alloc_blank(uint8_t type)
-{
-  void *out;
-
-  /* Allocate the space to keep the blank value. */
-  out=gal_data_malloc_array(type, 1);
-
-  /* Put the blank value in the allcated space. */
-  gal_data_set_blank(out, type);
-
-  /* Return the allocated space. */
-  return out;
-}
-
-
-
-
-
-/* Print the blank value as a string. */
-char *
-gal_data_blank_as_string(uint8_t type, int width)
-{
-  char *blank;
-  switch(type)
-    {
-    case GAL_DATA_TYPE_BIT:
-      error(EXIT_FAILURE, 0, "bit types are not implemented in "
-            "`gal_data_blank_as_string' yet.");
-      break;
-
-    case GAL_DATA_TYPE_STRING:
-      if(width)
-        asprintf(&blank, "%*s", width, GAL_DATA_BLANK_STRING);
-      else
-        asprintf(&blank, "%s", GAL_DATA_BLANK_STRING);
-      break;
-
-    case GAL_DATA_TYPE_UINT8:
-      if(width)
-        asprintf(&blank, "%*u", width, (uint8_t)GAL_DATA_BLANK_UINT8);
-      else
-        asprintf(&blank, "%u",         (uint8_t)GAL_DATA_BLANK_UINT8);
-      break;
-
-    case GAL_DATA_TYPE_INT8:
-      if(width)
-        asprintf(&blank, "%*d", width, (int8_t)GAL_DATA_BLANK_INT8);
-      else
-        asprintf(&blank, "%d",         (int8_t)GAL_DATA_BLANK_INT8);
-      break;
-
-    case GAL_DATA_TYPE_UINT16:
-      if(width)
-        asprintf(&blank, "%*u", width, (uint16_t)GAL_DATA_BLANK_UINT16);
-      else
-        asprintf(&blank, "%u",         (uint16_t)GAL_DATA_BLANK_UINT16);
-      break;
-
-    case GAL_DATA_TYPE_INT16:
-      if(width)
-        asprintf(&blank, "%*d", width, (int16_t)GAL_DATA_BLANK_INT16);
-      else
-        asprintf(&blank, "%d",         (int16_t)GAL_DATA_BLANK_INT16);
-      break;
-
-    case GAL_DATA_TYPE_UINT32:
-      if(width)
-        asprintf(&blank, "%*u", width, (uint32_t)GAL_DATA_BLANK_UINT32);
-      else
-        asprintf(&blank, "%u",         (uint32_t)GAL_DATA_BLANK_UINT32);
-      break;
-
-    case GAL_DATA_TYPE_INT32:
-      if(width)
-        asprintf(&blank, "%*d", width, (int32_t)GAL_DATA_BLANK_INT32);
-      else
-        asprintf(&blank, "%d",         (int32_t)GAL_DATA_BLANK_INT32);
-      break;
-
-    case GAL_DATA_TYPE_UINT64:
-      if(width)
-        asprintf(&blank, "%*lu", width, (uint64_t)GAL_DATA_BLANK_UINT64);
-      else
-        asprintf(&blank, "%lu",         (uint64_t)GAL_DATA_BLANK_UINT64);
-      break;
-
-    case GAL_DATA_TYPE_INT64:
-      if(width)
-        asprintf(&blank, "%*ld", width, (int64_t)GAL_DATA_BLANK_INT64);
-      else
-        asprintf(&blank, "%ld",         (int64_t)GAL_DATA_BLANK_INT64);
-      break;
-
-    case GAL_DATA_TYPE_FLOAT32:
-      if(width)
-        asprintf(&blank, "%*f", width, (float)GAL_DATA_BLANK_FLOAT32);
-      else
-        asprintf(&blank, "%f",         (float)GAL_DATA_BLANK_FLOAT32);
-      break;
-
-    case GAL_DATA_TYPE_FLOAT64:
-      if(width)
-        asprintf(&blank, "%*f", width, (double)GAL_DATA_BLANK_FLOAT64);
-      else
-        asprintf(&blank, "%f",         (double)GAL_DATA_BLANK_FLOAT64);
-      break;
-
-    default:
-      error(EXIT_FAILURE, 0, "type code %d not recognized in "
-            "`gal_data_blank_as_string'", type);
-    }
-  return blank;
-}
-
-
-
-
-
-
-/* For integers a simple equality is enough. */
-#define HAS_BLANK_INT(CTYPE, BLANK) {                                   \
-    CTYPE *a=data->array, *af=a+data->size;                             \
-    do if(*a++ == BLANK) return 1; while(a<af);                         \
-  }
-
-/* Note that a NaN value is not equal to another NaN value, so we can't use
-   the easy check for cases were the blank value is NaN. Also note that
-   `isnan' is actually a macro, so it works for both float and double
-   types.*/
-#define HAS_BLANK_FLT(CTYPE, BLANK, MULTIP) {                           \
-    CTYPE *a=data->array, *af=a+(MULTIP*data->size);                    \
-    if(isnan(BLANK)) do if(isnan(*a++)) return 1; while(a<af);          \
-    else             do if(*a++==BLANK) return 1; while(a<af);          \
-  }
-
-/* Return 1 if the dataset has a blank value and zero if it doesn't. */
-int
-gal_data_has_blank(gal_data_t *data)
-{
-  char **str=data->array, **strf=str+data->size;
-
-  /* Go over the pixels and check: */
-  switch(data->type)
-    {
-    case GAL_DATA_TYPE_BIT:
-      error(EXIT_FAILURE, 0, "Currently Gnuastro doesn't support bit "
-            "datatype, please get in touch with us to implement it.");
-
-    case GAL_DATA_TYPE_UINT8:
-      HAS_BLANK_INT(uint8_t, GAL_DATA_BLANK_UINT8);     break;
-
-    case GAL_DATA_TYPE_INT8:
-      HAS_BLANK_INT(int8_t, GAL_DATA_BLANK_INT8);       break;
-
-    case GAL_DATA_TYPE_UINT16:
-      HAS_BLANK_INT(uint16_t, GAL_DATA_BLANK_UINT16);   break;
-
-    case GAL_DATA_TYPE_INT16:
-      HAS_BLANK_INT(int16_t, GAL_DATA_BLANK_INT16);     break;
-
-    case GAL_DATA_TYPE_UINT32:
-      HAS_BLANK_INT(uint32_t, GAL_DATA_BLANK_UINT32);   break;
-
-    case GAL_DATA_TYPE_INT32:
-      HAS_BLANK_INT(int32_t, GAL_DATA_BLANK_INT32);     break;
-
-    case GAL_DATA_TYPE_UINT64:
-      HAS_BLANK_INT(uint64_t, GAL_DATA_BLANK_UINT64);   break;
-
-    case GAL_DATA_TYPE_INT64:
-      HAS_BLANK_INT(int64_t, GAL_DATA_BLANK_INT64);     break;
-
-    case GAL_DATA_TYPE_FLOAT32:
-      HAS_BLANK_FLT(float, GAL_DATA_BLANK_FLOAT32, 1);  break;
-
-    case GAL_DATA_TYPE_FLOAT64:
-      HAS_BLANK_FLT(double, GAL_DATA_BLANK_FLOAT64, 1); break;
-
-    case GAL_DATA_TYPE_COMPLEX32:
-      HAS_BLANK_FLT(float, GAL_DATA_BLANK_FLOAT32, 2);  break;
-
-    case GAL_DATA_TYPE_COMPLEX64:
-      HAS_BLANK_FLT(double, GAL_DATA_BLANK_FLOAT64, 2); break;
-
-    case GAL_DATA_TYPE_STRING:
-      do if(!strcmp(*str++,GAL_DATA_BLANK_STRING)) return 1; while(str<strf);
-      break;
-
-    default:
-      error(EXIT_FAILURE, 0, "a bug! type value (%d) not recognized "
-            "in `gal_data_blank_to_value'", data->type);
-    }
-
-  /* If there was a blank value, then the function would have returned with
-     a value of 1. So if it reaches here, then we can be sure that there
-     was no blank values, hence, return 0. */
-  return 0;
-}
-
-
-
-
-
-/* For integers a simple equality is enough. */
-#define FLAG_BLANK_INT(CTYPE, BLANK) {                                  \
-    CTYPE *a=data->array; do *o = (*a==BLANK); while(++o<of);           \
-  }
-
-/* Note that a NaN value is not equal to another NaN value, so we can't use
-   the easy check for cases were the blank value is NaN. Also note that
-   `isnan' is actually a macro, so it works for both float and double
-   types.*/
-#define FLAG_BLANK_FLT(CTYPE, BLANK) {                                  \
-    CTYPE *a=data->array;                                               \
-    if(isnan(BLANK)) do *o = isnan(*a++);   while(++o<of);              \
-    else             do *o = (*a++==BLANK); while(++o<of);              \
-  }
-
-#define FLAG_BLANK_COMPLEX(CTYPE, BLANK) {                              \
-    CTYPE *a=data->array;                                               \
-    if(isnan(BLANK))                                                    \
-      do { *o=(isnan(*a) || isnan(*(a+1))); a+=2; } while(++o<of);      \
-    else                                                                \
-      do { *o=(*a==BLANK || *(a+1)==BLANK); a+=2; } while(++o<of);      \
-  }
-
-/* Output a data-set of the the same size as the input, but with an uint8_t
-   type that has a value of 1 for data that are blank and 0 for those that
-   aren't. */
-gal_data_t *
-gal_data_flag_blank(gal_data_t *data)
-{
-  uint8_t *o, *of;
-  gal_data_t *out;
-  char **str=data->array, **strf=str+data->size;
-
-  /* Allocate the output array. */
-  out=gal_data_alloc(NULL, GAL_DATA_TYPE_UINT8, data->ndim, data->dsize,
-                     data->wcs, 0, data->minmapsize, data->name, data->unit,
-                     data->comment);
-
-  /* Set the pointers for easy looping. */
-  of=(o=out->array)+data->size;
-
-  /* Go over the pixels and set the output values. */
-  switch(data->type)
-    {
-    case GAL_DATA_TYPE_BIT:
-      error(EXIT_FAILURE, 0, "Currently Gnuastro doesn't support bit "
-            "datatype, please get in touch with us to implement it.");
-
-    case GAL_DATA_TYPE_UINT8:
-      FLAG_BLANK_INT(uint8_t, GAL_DATA_BLANK_UINT8);      break;
-
-    case GAL_DATA_TYPE_INT8:
-      FLAG_BLANK_INT(int8_t, GAL_DATA_BLANK_INT8);        break;
-
-    case GAL_DATA_TYPE_UINT16:
-      FLAG_BLANK_INT(uint16_t, GAL_DATA_BLANK_UINT16);    break;
-
-    case GAL_DATA_TYPE_INT16:
-      FLAG_BLANK_INT(int16_t, GAL_DATA_BLANK_INT16);      break;
-
-    case GAL_DATA_TYPE_UINT32:
-      FLAG_BLANK_INT(uint32_t, GAL_DATA_BLANK_UINT32);    break;
-
-    case GAL_DATA_TYPE_INT32:
-      FLAG_BLANK_INT(int32_t, GAL_DATA_BLANK_INT32);      break;
-
-    case GAL_DATA_TYPE_UINT64:
-      FLAG_BLANK_INT(uint64_t, GAL_DATA_BLANK_UINT64);    break;
-
-    case GAL_DATA_TYPE_INT64:
-      FLAG_BLANK_INT(int64_t, GAL_DATA_BLANK_INT64);      break;
-
-    case GAL_DATA_TYPE_FLOAT32:
-      FLAG_BLANK_FLT(float, GAL_DATA_BLANK_FLOAT32);      break;
-
-    case GAL_DATA_TYPE_FLOAT64:
-      FLAG_BLANK_FLT(double, GAL_DATA_BLANK_FLOAT64);     break;
-
-    case GAL_DATA_TYPE_COMPLEX32:
-      FLAG_BLANK_COMPLEX(float, GAL_DATA_BLANK_FLOAT32);  break;
-
-    case GAL_DATA_TYPE_COMPLEX64:
-      FLAG_BLANK_COMPLEX(double, GAL_DATA_BLANK_FLOAT64); break;
-
-    case GAL_DATA_TYPE_STRING:
-      do *o++ = !strcmp(*str,GAL_DATA_BLANK_STRING); while(++str<strf);
-      break;
-
-    default:
-      error(EXIT_FAILURE, 0, "type value (%d) not recognized "
-            "in `gal_data_flag_blank'", data->type);
-    }
-
-  /* Return */
-  return out;
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
 
 
 
@@ -1447,8 +1063,8 @@ gal_data_flag_blank(gal_data_t *data)
 static void
 data_copy_to_string_not_parsed(char *string, void *to, uint8_t type)
 {
-  if( strcmp(string, GAL_DATA_BLANK_STRING) )
-    gal_data_set_blank(to, type);
+  if( strcmp(string, GAL_BLANK_STRING) )
+    gal_blank_write(to, type);
   else
     error(EXIT_FAILURE, 0, "`%s' couldn't be parsed as `%s' type",
           string, gal_data_type_as_string(type, 1));
@@ -1536,7 +1152,7 @@ data_copy_from_string(gal_data_t *from, gal_data_t *to)
       {                                                                 \
         if(a[i]!=BLANK) asprintf(&strarr[i], FMT, a[i]);                \
         else                                                            \
-          gal_checkset_allocate_copy(GAL_DATA_BLANK_STRING, &strarr[i]); \
+          gal_checkset_allocate_copy(GAL_BLANK_STRING, &strarr[i]); \
       }                                                                 \
   }
 
@@ -1547,7 +1163,7 @@ data_copy_from_string(gal_data_t *from, gal_data_t *to)
         if(isnan(BLANK)) isblank = isnan(a[i]) ? 1 : 0;                 \
         else             isblank = a[i]==BLANK ? 1 : 0;                 \
         if(isblank==0) asprintf(&strarr[i], "%f", a[i]);                \
-        else gal_checkset_allocate_copy(GAL_DATA_BLANK_STRING, &strarr[i]); \
+        else gal_checkset_allocate_copy(GAL_BLANK_STRING, &strarr[i]); \
       }                                                                 \
   }
 
@@ -1569,34 +1185,34 @@ data_copy_to_string(gal_data_t *from, gal_data_t *to)
   switch(from->type)
     {
     case GAL_DATA_TYPE_UINT8:
-      COPY_TO_STR_INT(uint8_t,  GAL_DATA_BLANK_UINT8, "%u");    break;
+      COPY_TO_STR_INT(uint8_t,  GAL_BLANK_UINT8, "%u");    break;
 
     case GAL_DATA_TYPE_INT8:
-      COPY_TO_STR_INT(int8_t,   GAL_DATA_BLANK_INT8, "%d");     break;
+      COPY_TO_STR_INT(int8_t,   GAL_BLANK_INT8, "%d");     break;
 
     case GAL_DATA_TYPE_UINT16:
-      COPY_TO_STR_INT(uint16_t, GAL_DATA_BLANK_UINT16, "%u");   break;
+      COPY_TO_STR_INT(uint16_t, GAL_BLANK_UINT16, "%u");   break;
 
     case GAL_DATA_TYPE_INT16:
-      COPY_TO_STR_INT(int16_t,  GAL_DATA_BLANK_INT16, "%d");    break;
+      COPY_TO_STR_INT(int16_t,  GAL_BLANK_INT16, "%d");    break;
 
     case GAL_DATA_TYPE_UINT32:
-      COPY_TO_STR_INT(uint32_t, GAL_DATA_BLANK_UINT32, "%u");   break;
+      COPY_TO_STR_INT(uint32_t, GAL_BLANK_UINT32, "%u");   break;
 
     case GAL_DATA_TYPE_INT32:
-      COPY_TO_STR_INT(int32_t,  GAL_DATA_BLANK_INT32, "%d");    break;
+      COPY_TO_STR_INT(int32_t,  GAL_BLANK_INT32, "%d");    break;
 
     case GAL_DATA_TYPE_UINT64:
-      COPY_TO_STR_INT(uint64_t, GAL_DATA_BLANK_UINT64, "%lu");  break;
+      COPY_TO_STR_INT(uint64_t, GAL_BLANK_UINT64, "%lu");  break;
 
     case GAL_DATA_TYPE_INT64:
-      COPY_TO_STR_INT(int64_t,  GAL_DATA_BLANK_INT64, "%ld");   break;
+      COPY_TO_STR_INT(int64_t,  GAL_BLANK_INT64, "%ld");   break;
 
     case GAL_DATA_TYPE_FLOAT32:
-      COPY_TO_STR_FLT(float, GAL_DATA_BLANK_FLOAT32);           break;
+      COPY_TO_STR_FLT(float, GAL_BLANK_FLOAT32);           break;
 
     case GAL_DATA_TYPE_FLOAT64:
-      COPY_TO_STR_FLT(double, GAL_DATA_BLANK_FLOAT32);          break;
+      COPY_TO_STR_FLT(double, GAL_BLANK_FLOAT32);          break;
 
     case GAL_DATA_TYPE_STRING:
       for(i=0;i<from->size;++i)
@@ -1629,11 +1245,11 @@ data_copy_to_string(gal_data_t *from, gal_data_t *to)
                                                                         \
     /* Check if there are blank values in the input array and that */   \
     /* the types of the two structures are different. */                \
-    if( in->type!=newtype && gal_data_has_blank(in) )                   \
+    if( in->type!=newtype && gal_blank_present(in) )                    \
       {                                                                 \
         /* Set the blank values */                                      \
-        gal_data_set_blank(&iblank, in->type);                          \
-        gal_data_set_blank(&oblank, newtype);                           \
+        gal_blank_write(&iblank, in->type);                             \
+        gal_blank_write(&oblank, newtype);                              \
                                                                         \
         /* Copy the input to the output. */                             \
         do { *oa = *ia==iblank ? oblank : *ia; ia++; } while(++oa<of);  \
@@ -1656,11 +1272,11 @@ data_copy_to_string(gal_data_t *from, gal_data_t *to)
                                                                         \
     /* Check if there are blank values in the input array and that */   \
     /* the types of the two structures are different. */                \
-    if( in->type!=newtype && gal_data_has_blank(in) )                   \
+    if( in->type!=newtype && gal_blank_present(in) )                    \
       {                                                                 \
         /* Set the blank values */                                      \
-        gal_data_set_blank(&iblank, in->type);                          \
-        gal_data_set_blank(&oblank, newtype);                           \
+        gal_blank_write(&iblank, in->type);                             \
+        gal_blank_write(&oblank, newtype);                              \
                                                                         \
         /* When the blank value isn't NaN, then we should use the */    \
         /* equal operator to check for blank values. */                 \
@@ -1930,8 +1546,8 @@ gal_data_write_to_string(void *ptr, uint8_t type, int 
quote_if_str_has_space)
     case GAL_DATA_TYPE_INT32:   WRITE_TO_STRING(int32_t,   "%d");  break;
     case GAL_DATA_TYPE_UINT64:  WRITE_TO_STRING(uint64_t, "%lu");  break;
     case GAL_DATA_TYPE_INT64:   WRITE_TO_STRING(int64_t,  "%ld");  break;
-    case GAL_DATA_TYPE_FLOAT32: WRITE_TO_STRING(float,   "%.6f");  break;
-    case GAL_DATA_TYPE_FLOAT64: WRITE_TO_STRING(double, "%.10f");  break;
+    case GAL_DATA_TYPE_FLOAT32: WRITE_TO_STRING(float,   "%.6g");  break;
+    case GAL_DATA_TYPE_FLOAT64: WRITE_TO_STRING(double, "%.10g");  break;
 
     default:
       error(EXIT_FAILURE, 0, "type code %d not recognized in "
diff --git a/lib/fits.c b/lib/fits.c
index b3cb9ca..4423254 100644
--- a/lib/fits.c
+++ b/lib/fits.c
@@ -36,6 +36,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/git.h>
 #include <gnuastro/wcs.h>
 #include <gnuastro/fits.h>
+#include <gnuastro/blank.h>
 
 #include "checkset.h"
 #include "fixedstringmacros.h"
@@ -1182,10 +1183,12 @@ gal_fits_img_read(char *filename, char *hdu, size_t 
minmapsize)
      no images). */
   if(ndim==0)
     error(EXIT_FAILURE, 0, "%s (hdu: %s) has 0 dimensions! The most common "
-          "cause for this is a wrongly specified HDU: in some FITS images, "
-          "the first HDU doesn't have any data. So probably reading the "
-          "second HDU (with `--hdu=1' or `-h1') will solve the problem. Note "
-          "that currently HDU counting starts from 0." , filename, hdu);
+          "cause for this is a wrongly specified HDU. In some FITS images, "
+          "the first HDU doesn't have any data, the data is in subsequent "
+          "extensions. So probably reading the second HDU (with `--hdu=1' "
+          "or `-h1') will solve the problem (following CFITSIO's "
+          "convention, currently HDU counting starts from 0)." , filename,
+          hdu);
 
 
   /* Set the fpixel array (first pixel in all dimensions): */
@@ -1210,7 +1213,7 @@ gal_fits_img_read(char *filename, char *hdu, size_t 
minmapsize)
   /* Allocate the space for the array and for the blank values. */
   img=gal_data_alloc(NULL, type, ndim, dsize, NULL, 0, minmapsize,
                      name, unit, NULL);
-  blank=gal_data_alloc_blank(type);
+  blank=gal_blank_alloc_write(type);
   gal_data_free_ll(keysll);
   free(dsize);
 
@@ -1361,7 +1364,7 @@ gal_fits_img_write_to_ptr(gal_data_t *data, char 
*filename)
 
   /* If we have blank pixels, we need to define a BLANK keyword when we are
      dealing with integer types. */
-  if(gal_data_has_blank(data))
+  if(gal_blank_present(data))
     switch(data->type)
       {
       case GAL_DATA_TYPE_FLOAT32:
@@ -1371,7 +1374,7 @@ gal_fits_img_write_to_ptr(gal_data_t *data, char 
*filename)
         break;
 
       default:
-        blank=gal_data_alloc_blank(data->type);
+        blank=gal_blank_alloc_write(data->type);
         if(fits_write_key(fptr, datatype, "BLANK", blank,
                           "Pixels with no data.", &status) )
           gal_fits_io_error(status, "adding the BLANK keyword");
@@ -1816,10 +1819,6 @@ gal_fits_tab_info(char *filename, char *hdu, size_t 
*numcols,
               /* Write the type into the data structure. */
               allcols[index].type=gal_fits_datatype_to_type(datatype, 1);
 
-              printf("%zu: %s (%d)\n", index,
-                     gal_data_type_as_string(allcols[index].type, 1),
-                     datatype);
-
               /* If we are dealing with a string type, we need to know the
                  number of bytes in both cases for printing later. */
               if( allcols[index].type==GAL_DATA_TYPE_STRING )
@@ -1992,7 +1991,7 @@ gal_fits_tab_read(char *filename, char *hdu, size_t 
numrows,
 
       /* Allocate a blank value for the given type and read/store the
          column using CFITSIO. Afterwards, free the blank value. */
-      blank=gal_data_alloc_blank(out->type);
+      blank=gal_blank_alloc_write(out->type);
       fits_read_col(fptr, gal_fits_type_to_datatype(out->type), ind->v+1,
                     1, 1, out->size, blank, out->array, &anynul, &status);
       gal_fits_io_error(status, NULL);
@@ -2060,8 +2059,8 @@ fits_table_prepare_arrays(gal_data_t *cols, size_t 
numcols, int tabletype,
                expected width or not. Its initial width is set the output
                of the function above, but if the value is larger,
                `asprintf' (which is used) will make it wider. */
-            blank = ( gal_data_has_blank(col)
-                      ? gal_data_blank_as_string(col->type, col->disp_width)
+            blank = ( gal_blank_present(col)
+                      ? gal_blank_as_string(col->type, col->disp_width)
                       : NULL );
 
             /* Adjust the width. */
@@ -2147,7 +2146,7 @@ fits_write_tnull_tcomm(fitsfile *fptr, gal_data_t *col, 
int tabletype,
 
       /* Print the keyword and value. */
       asprintf(&keyname, "TNULL%zu", colnum);
-      blank=gal_data_blank_as_string(col->type, col->disp_width);
+      blank=gal_blank_as_string(col->type, col->disp_width);
 
       /* When in exponential form (`tform' starting with `E'), CFITSIO
          writes a NaN value as `NAN', but when in floating point form
@@ -2173,7 +2172,7 @@ fits_write_tnull_tcomm(fitsfile *fptr, gal_data_t *col, 
int tabletype,
           && col->type!=GAL_DATA_TYPE_FLOAT64
           && col->type!=GAL_DATA_TYPE_STRING )
         {
-          blank=gal_data_alloc_blank(col->type);
+          blank=gal_blank_alloc_write(col->type);
           asprintf(&keyname, "TNULL%zu", colnum);
           fits_write_key(fptr, gal_fits_type_to_datatype(col->type),
                          keyname, blank, "blank value for this column",
@@ -2265,8 +2264,8 @@ gal_fits_tab_write(gal_data_t *cols, char *comments, int 
tabletype,
 
       /* Set the blank pointer if its necessary, note that strings don't
          need a blank pointer in a FITS ASCII table.*/
-      blank = ( gal_data_has_blank(col)
-                ? gal_data_alloc_blank(col->type) : NULL );
+      blank = ( gal_blank_present(col)
+                ? gal_blank_alloc_write(col->type) : NULL );
       if(tabletype==GAL_TABLE_FORMAT_AFITS && col->type==GAL_DATA_TYPE_STRING)
         { if(blank) free(blank); blank=NULL; }
 
diff --git a/lib/gnuastro/arithmetic.h b/lib/gnuastro/arithmetic.h
index ec2ab40..9243f18 100644
--- a/lib/gnuastro/arithmetic.h
+++ b/lib/gnuastro/arithmetic.h
@@ -106,22 +106,30 @@ enum gal_arithmetic_operators
 
   GAL_ARITHMETIC_OP_MINVAL,       /* Minimum value of array.               */
   GAL_ARITHMETIC_OP_MAXVAL,       /* Maximum value of array.               */
+  GAL_ARITHMETIC_OP_NUMVAL,       /* Number of (non-blank) elements.       */
+  GAL_ARITHMETIC_OP_SUMVAL,       /* Sum of (non-blank) elements.          */
+  GAL_ARITHMETIC_OP_MEANVAL,      /* Mean value of array.                  */
+  GAL_ARITHMETIC_OP_STDVAL,       /* Standard deviation value of array.    */
+  GAL_ARITHMETIC_OP_MEDIANVAL,    /* Median value of array.                */
+
   GAL_ARITHMETIC_OP_MIN,          /* Minimum per pixel of multiple arrays. */
   GAL_ARITHMETIC_OP_MAX,          /* Maximum per pixel of multiple arrays. */
+  GAL_ARITHMETIC_OP_NUM,          /* Non-blank number of pixels in arrays. */
   GAL_ARITHMETIC_OP_SUM,          /* Sum per pixel of multiple arrays.     */
-  GAL_ARITHMETIC_OP_AVERAGE,      /* Average per pixel of multiple arrays. */
+  GAL_ARITHMETIC_OP_MEAN,         /* Mean per pixel of multiple arrays.    */
+  GAL_ARITHMETIC_OP_STD,          /* STD per pixel of multiple arrays.     */
   GAL_ARITHMETIC_OP_MEDIAN,       /* Median per pixel of multiple arrays.  */
 
-  GAL_ARITHMETIC_OP_TO_UINT8,     /* Convert to unsigned char.             */
-  GAL_ARITHMETIC_OP_TO_INT8,      /* Convert to char.                      */
-  GAL_ARITHMETIC_OP_TO_UINT16,    /* Convert to unsigned short.            */
-  GAL_ARITHMETIC_OP_TO_INT16,     /* Convert to short.                     */
-  GAL_ARITHMETIC_OP_TO_UINT32,    /* Convert to unsigned int.              */
-  GAL_ARITHMETIC_OP_TO_INT32,     /* Convert to int.                       */
-  GAL_ARITHMETIC_OP_TO_UINT64,    /* Convert to unsigned long.             */
-  GAL_ARITHMETIC_OP_TO_INT64,     /* Convert to long.                      */
-  GAL_ARITHMETIC_OP_TO_FLOAT32,   /* Convert to float.                     */
-  GAL_ARITHMETIC_OP_TO_FLOAT64,   /* Convert to double.                    */
+  GAL_ARITHMETIC_OP_TO_UINT8,     /* Convert to uint8_t.                   */
+  GAL_ARITHMETIC_OP_TO_INT8,      /* Convert to int8_t.                    */
+  GAL_ARITHMETIC_OP_TO_UINT16,    /* Convert to uint16_t.                  */
+  GAL_ARITHMETIC_OP_TO_INT16,     /* Convert to int16_t.                   */
+  GAL_ARITHMETIC_OP_TO_UINT32,    /* Convert to uint32_t.                  */
+  GAL_ARITHMETIC_OP_TO_INT32,     /* Convert to int32_t.                   */
+  GAL_ARITHMETIC_OP_TO_UINT64,    /* Convert to uint64_t.                  */
+  GAL_ARITHMETIC_OP_TO_INT64,     /* Convert to int64_t.                   */
+  GAL_ARITHMETIC_OP_TO_FLOAT32,   /* Convert to float32.                   */
+  GAL_ARITHMETIC_OP_TO_FLOAT64,   /* Convert to float64.                   */
 };
 
 
diff --git a/lib/gnuastro/array.h b/lib/gnuastro/array.h
deleted file mode 100644
index 5975078..0000000
--- a/lib/gnuastro/array.h
+++ /dev/null
@@ -1,157 +0,0 @@
-/*********************************************************************
-Functions to manipulate arrays.
-This is part of GNU Astronomy Utilities (Gnuastro) package.
-
-Original author:
-     Mohammad Akhlaghi <address@hidden>
-Contributing author(s):
-Copyright (C) 2015, 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_ARRAY_H__
-#define __GAL_ARRAY_H__
-
-/* Include other headers if necessary here. Note that other header files
-   must be included before the C++ preparations below */
-
-
-
-/* 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_array_uchar_init_on_region(unsigned char *in, const unsigned char v,
-                                    size_t start, size_t s0, size_t s1,
-                                    size_t is1);
-
-void
-gal_array_long_init(long *in, size_t size, const long v);
-
-void
-gal_array_long_init_on_region(long *in, const long v, size_t start,
-                              size_t s0, size_t s1, size_t is1);
-
-void
-gal_array_uchar_copy(unsigned char *in, size_t size, unsigned char **out);
-
-void
-gal_array_float_copy(float *in, size_t size, float **out);
-
-void
-gal_array_float_copy_noalloc(float *in, size_t size, float *out);
-
-void
-gal_array_fset_const(float *in, size_t size, float a);
-
-void
-gal_array_freplace_value(float *in, size_t size, float from, float to);
-
-void
-gal_array_freplace_nonnans(float *in, size_t size, float to);
-
-void
-gal_array_no_nans(float *in, size_t *size);
-
-void
-gal_array_no_nans_double(double *in, size_t *size);
-
-void
-gal_array_uchar_replace(unsigned char *in, size_t size,
-                        unsigned char from, unsigned char to);
-
-void
-gal_array_fmultip_const(float *in, size_t size, float a);
-
-void
-gal_array_fsum_const(float *in, size_t size, float a);
-
-float *
-gal_array_fsum_arrays(float *in1, float *in2, size_t size);
-
-
-void
-gal_array_dmultip_const(double *in, size_t size, double a);
-
-void
-gal_array_dmultip_arrays(double *in1, double *in2, size_t size);
-
-
-void
-gal_array_ddivide_const(double *in, size_t size, double a);
-
-void
-gal_array_dconst_divide(double *in, size_t size, double a);
-
-void
-gal_array_ddivide_arrays(double *in1, double *in2, size_t size);
-
-
-void
-gal_array_dsum_const(double *in, size_t size, double a);
-
-void
-gal_array_dsum_arrays(double *in1, double *in2, size_t size);
-
-
-void
-gal_array_dsubtract_const(double *in, size_t size, double a);
-
-void
-gal_array_dconst_subtract(double *in, size_t size, double a);
-
-void
-gal_array_dsubtract_arrays(double *in1, double *in2, size_t size);
-
-
-void
-gal_array_dpower_const(double *in, size_t size, double a);
-
-void
-gal_array_dconst_power(double *in, size_t size, double a);
-
-void
-gal_array_dpower_arrays(double *in1, double *in2, size_t size);
-
-
-void
-gal_array_dlog_array(double *in1, size_t size);
-
-void
-gal_array_dlog10_array(double *in1, size_t size);
-
-void
-gal_array_dabs_array(double *in1, size_t size);
-
-
-
-__END_C_DECLS    /* From C++ preparations */
-
-#endif           /* __GAL_ARRAY_H__ */
diff --git a/lib/gnuastro/blank.h b/lib/gnuastro/blank.h
new file mode 100644
index 0000000..2ee0716
--- /dev/null
+++ b/lib/gnuastro/blank.h
@@ -0,0 +1,91 @@
+/*********************************************************************
+blank -- Deal with blank values in a datasets
+This is part of GNU Astronomy Utilities (Gnuastro) package.
+
+Original author:
+     Mohammad Akhlaghi <address@hidden>
+Contributing author(s):
+Copyright (C) 2017, 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_BLANK_H__
+#define __GAL_BLANK_H__
+
+/* Include other headers if necessary here. Note that other header files
+   must be included before the C++ preparations below */
+#include <gnuastro/data.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 */
+
+
+/* Macros: */
+
+/* Blank values: Note that for the unsigned types or small types (like
+   char), the maximum value is considered as a blank value, since the
+   minimum value of an unsigned type is zero and zero is often meaningful
+   in contexts were unsigned values are used. */
+#define GAL_BLANK_UINT8      UINT8_MAX
+#define GAL_BLANK_INT8       INT8_MIN
+#define GAL_BLANK_UINT16     UINT16_MAX
+#define GAL_BLANK_INT16      INT16_MIN
+#define GAL_BLANK_UINT32     UINT32_MAX
+#define GAL_BLANK_INT32      INT32_MIN
+#define GAL_BLANK_UINT64     UINT64_MAX
+#define GAL_BLANK_INT64      INT64_MIN
+#define GAL_BLANK_FLOAT32    NAN
+#define GAL_BLANK_FLOAT64    NAN
+#define GAL_BLANK_STRING     "n/a"
+
+
+
+
+/* Functions. */
+void
+gal_blank_write(void *pointer, uint8_t type);
+
+void *
+gal_blank_alloc_write(uint8_t type);
+
+char *
+gal_blank_as_string(uint8_t type, int width);
+
+int
+gal_blank_present(gal_data_t *data);
+
+gal_data_t *
+gal_blank_flag(gal_data_t *data);
+
+
+
+
+
+__END_C_DECLS    /* From C++ preparations */
+
+#endif           /* __GAL_DATA_H__ */
diff --git a/lib/gnuastro/data.h b/lib/gnuastro/data.h
index c9f16ec..7c7fae5 100644
--- a/lib/gnuastro/data.h
+++ b/lib/gnuastro/data.h
@@ -67,28 +67,6 @@ __BEGIN_C_DECLS  /* From C++ preparations */
 
 
 
-/* Macros: */
-
-/* Blank values: Note that for the unsigned types or small types (like
-   char), the maximum value is considered as a blank value, since the
-   minimum value of an unsigned type is zero and zero is often meaningful
-   in contexts were unsigned values are used. */
-#define GAL_DATA_BLANK_UINT8      UINT8_MAX
-#define GAL_DATA_BLANK_INT8       INT8_MIN
-#define GAL_DATA_BLANK_UINT16     UINT16_MAX
-#define GAL_DATA_BLANK_INT16      INT16_MIN
-#define GAL_DATA_BLANK_UINT32     UINT32_MAX
-#define GAL_DATA_BLANK_INT32      INT32_MIN
-#define GAL_DATA_BLANK_UINT64     UINT64_MAX
-#define GAL_DATA_BLANK_INT64      INT64_MIN
-#define GAL_DATA_BLANK_FLOAT32    NAN
-#define GAL_DATA_BLANK_FLOAT64    NAN
-#define GAL_DATA_BLANK_STRING     "n/a"
-
-
-
-
-
 /* Macros to identify the type of data. */
 enum gal_data_types
 {
@@ -266,35 +244,6 @@ gal_data_free_ll(gal_data_t *list);
 
 
 
-
-/*************************************************************
- **************          Blank data            ***************
- *************************************************************/
-void
-gal_data_set_blank(void *pointer, uint8_t type);
-
-void *
-gal_data_alloc_blank(uint8_t type);
-
-char *
-gal_data_blank_as_string(uint8_t type, int width);
-
-void
-gal_data_apply_mask(gal_data_t *in, gal_data_t *mask);
-
-void
-gal_data_blank_to_value(gal_data_t *data, void *value);
-
-int
-gal_data_has_blank(gal_data_t *data);
-
-gal_data_t *
-gal_data_flag_blank(gal_data_t *data);
-
-
-
-
-
 /*************************************************************
  **************            Copying             ***************
  *************************************************************/
@@ -315,6 +264,7 @@ gal_data_to_same_type(gal_data_t *f, gal_data_t *s, 
gal_data_t **of,
                       gal_data_t **os, uint8_t type, int freeinputs);
 
 
+
 /*************************************************************
  **************              Write             ***************
  *************************************************************/
diff --git a/lib/gnuastro/qsort.h b/lib/gnuastro/qsort.h
index cd59910..0cb532d 100644
--- a/lib/gnuastro/qsort.h
+++ b/lib/gnuastro/qsort.h
@@ -60,70 +60,64 @@ gal_qsort_index_float_decreasing(const void * a, const void 
* b);
 
 
 int
-gal_qsort_uchar_decreasing(const void * a, const void * b);
+gal_qsort_uint8_decreasing(const void *a, const void *b);
 
 int
-gal_qsort_uchar_increasing(const void * a, const void * b);
+gal_qsort_uint8_increasing(const void *a, const void *b);
 
 int
-gal_qsort_char_decreasing(const void * a, const void * b);
+gal_qsort_int8_decreasing(const void *a, const void *b);
 
 int
-gal_qsort_char_increasing(const void * a, const void * b);
+gal_qsort_int8_increasing(const void *a, const void *b);
 
 int
-gal_qsort_ushort_decreasing(const void * a, const void * b);
+gal_qsort_uint16_decreasing(const void *a, const void *b);
 
 int
-gal_qsort_ushort_increasing(const void * a, const void * b);
+gal_qsort_uint16_increasing(const void *a, const void *b);
 
 int
-gal_qsort_short_decreasing(const void * a, const void * b);
+gal_qsort_int16_decreasing(const void *a, const void *b);
 
 int
-gal_qsort_short_increasing(const void * a, const void * b);
+gal_qsort_int16_increasing(const void *a, const void *b);
 
 int
-gal_qsort_uint_decreasing(const void * a, const void * b);
+gal_qsort_uint32_decreasing(const void *a, const void *b);
 
 int
-gal_qsort_uint_increasing(const void * a, const void * b);
+gal_qsort_uint32_increasing(const void *a, const void *b);
 
 int
-gal_qsort_int_decreasing(const void * a, const void * b);
+gal_qsort_int32_decreasing(const void *a, const void *b);
 
 int
-gal_qsort_int_increasing(const void * a, const void * b);
+gal_qsort_int32_increasing(const void *a, const void *b);
 
 int
-gal_qsort_ulong_decreasing(const void * a, const void * b);
+gal_qsort_uint64_decreasing(const void *a, const void *b);
 
 int
-gal_qsort_ulong_increasing(const void * a, const void * b);
+gal_qsort_uint64_increasing(const void *a, const void *b);
 
 int
-gal_qsort_long_decreasing(const void * a, const void * b);
+gal_qsort_int64_decreasing(const void *a, const void *b);
 
 int
-gal_qsort_long_increasing(const void * a, const void * b);
+gal_qsort_int64_increasing(const void *a, const void *b);
 
 int
-gal_qsort_longlong_decreasing(const void * a, const void * b);
+gal_qsort_float32_decreasing(const void *a, const void *b);
 
 int
-gal_qsort_longlong_increasing(const void * a, const void * b);
+gal_qsort_float32_increasing(const void *a, const void *b);
 
 int
-gal_qsort_float_decreasing(const void * a, const void * b);
+gal_qsort_float64_decreasing(const void *a, const void *b);
 
 int
-gal_qsort_float_increasing(const void * a, const void * b);
-
-int
-gal_qsort_double_decreasing(const void * a, const void * b);
-
-int
-gal_qsort_double_increasing(const void * a, const void * b);
+gal_qsort_float64_increasing(const void *a, const void *b);
 
 
 
diff --git a/lib/gnuastro/table.h b/lib/gnuastro/table.h
index 711d912..4471144 100644
--- a/lib/gnuastro/table.h
+++ b/lib/gnuastro/table.h
@@ -121,7 +121,18 @@ enum gal_table_diplay_formats
 
 
 
-/* Functions */
+/************************************************************************/
+/***************              Error messages              ***************/
+/************************************************************************/
+void
+gal_table_error_col_selection(char *filename, char *hdu, char *errorstring);
+
+
+
+
+/************************************************************************/
+/***************                 Formats                  ***************/
+/************************************************************************/
 uint8_t
 gal_table_string_to_format(char *string);
 
@@ -137,9 +148,11 @@ gal_table_searchin_as_string(uint8_t searchin);
 void
 gal_table_check_fits_format(char *filename, int tableformat);
 
-void
-gal_table_too_many_columns(char *filename);
 
+
+/************************************************************************/
+/***************          Printing information            ***************/
+/************************************************************************/
 void
 gal_table_print_info(gal_data_t *allcols, size_t numcols, size_t numrows);
 
@@ -150,14 +163,30 @@ gal_table_col_print_info(gal_data_t *col, int tableformat,
 void
 gal_table_read_blank(gal_data_t *col, char *blank);
 
+
+
+
+/************************************************************************/
+/***************         Information about a table        ***************/
+/************************************************************************/
 gal_data_t *
 gal_table_info(char *filename, char *hdu, size_t *numcols,
                size_t *numrows, int *tableformat);
 
+
+
+/************************************************************************/
+/***************               Read a table               ***************/
+/************************************************************************/
 gal_data_t *
 gal_table_read(char *filename, char *hdu, struct gal_linkedlist_stll *cols,
                int searchin, int ignorecase, int minmapsize);
 
+
+
+/************************************************************************/
+/***************              Write a table               ***************/
+/************************************************************************/
 void
 gal_table_write(gal_data_t *cols, char *comments, int tableformat,
                 char *filename, int dontdelete);
diff --git a/lib/mesh.c b/lib/mesh.c
index 4bed78f..0b81c7c 100644
--- a/lib/mesh.c
+++ b/lib/mesh.c
@@ -1209,7 +1209,7 @@ meshinterponthread(void *inparams)
 
       /* Find the median of the nearest neighbors and put it in: */
       qsort(nearest1, numnearest, sizeof *nearest1,
-            gal_qsort_float_increasing);
+            gal_qsort_float32_increasing);
       outgarray1[thisind] = ( numnearest%2 ?
                               nearest1[numnearest/2] : /* Odd.  */
                               (nearest1[numnearest/2]  /* Even. */
@@ -1217,7 +1217,7 @@ meshinterponthread(void *inparams)
       if(ngarrays==2)
         {
           qsort(nearest2, numnearest, sizeof *nearest2,
-                gal_qsort_float_increasing);
+                gal_qsort_float32_increasing);
           outgarray2[thisind] = ( numnearest%2 ?
                                   nearest2[numnearest/2] : /* Odd.  */
                                   (nearest2[numnearest/2]  /* Even. */
diff --git a/lib/options.c b/lib/options.c
index 7a69271..3f0885c 100644
--- a/lib/options.c
+++ b/lib/options.c
@@ -174,6 +174,114 @@ options_get_home()
 
 
 
+/* The input to this function is a string of any number of numbers
+   separated by a comma (`,') and possibly containing fractions, for
+   example: `1,2/3, 4.95'. The output `gal_data_t' contains the array of
+   given values in `double' type. You can read the number from its `size'
+   element. */
+gal_data_t *
+gal_options_parse_list_of_numbers(char *string, char *filename, size_t lineno)
+{
+  size_t i, num=0;
+  gal_data_t *out;
+  char *c=string, *tailptr;
+  double numerator=NAN, denominator=NAN, tmp;
+  struct gal_linkedlist_dll *list=NULL, *tdll;
+
+
+  /* The nature of the arrays/numbers read here is very small, so since
+     `p->cp.minmapsize' might not have been read yet, we will set it to -1
+     (largest size_t number), so the values are kept in memory. */
+  size_t minmapsize=-1;
+
+
+  /* Go through the input character by character. */
+  while(*c!='\0')
+    {
+      switch(*c)
+        {
+
+        /* Ignore space or tab. */
+        case ' ':
+        case '\t':
+          ++c;
+          break;
+
+        /* Comma marks the transition to the next number. */
+        case ',':
+          if(isnan(numerator))
+            error_at_line(EXIT_FAILURE, 0, filename, lineno, "a number "
+                          "must be given before `,'. You have given: `%s'",
+                          string);
+          gal_linkedlist_add_to_dll(&list, isnan(denominator)
+                                    ? numerator : numerator/denominator);
+          numerator=denominator=NAN;
+          ++num;
+          ++c;
+          break;
+
+        /* Divide two numbers. */
+        case '/':
+          if( isnan(numerator) || !isnan(denominator) )
+            error_at_line(EXIT_FAILURE, 0, filename, lineno, "`/' must "
+                          "only be between two numbers and used for "
+                          "division. But you have given `%s'", string);
+          ++c;
+          break;
+
+        /* Read the number. */
+        default:
+
+          /* Parse the string. */
+          tmp=strtod(c, &tailptr);
+          if(tailptr==c)
+            error_at_line(EXIT_FAILURE, 0, filename, lineno, "the first "
+                          "part of `%s' couldn't be read as a number. This "
+                          "was part of `%s'", c, string);
+
+          /* See if the number should be put in the numerator or
+             denominator. */
+          if(isnan(numerator)) numerator=tmp;
+          else
+            {
+              if(isnan(denominator)) denominator=tmp;
+              else error_at_line(EXIT_FAILURE, 0, filename, lineno, "more "
+                                 "than two numbers in each element.");
+            }
+
+          /* Set `c' to tailptr. */
+          c=tailptr;
+        }
+    }
+
+
+  /* If the last number wasn't finished by a `,', add the read value to the
+     list */
+  if( !isnan(numerator) )
+    {
+      ++num;
+      gal_linkedlist_add_to_dll(&list, isnan(denominator)
+                                ? numerator : numerator/denominator);
+    }
+
+
+  /* Allocate the output data structure and fill it up. */
+  i=num;
+  out=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT64, 1, &num, NULL, 0,
+                     minmapsize, NULL, NULL, NULL);
+  for(tdll=list;tdll!=NULL;tdll=tdll->next)
+    ((double *)out->array)[--i]=tdll->v;
+
+
+  /* Clean up and return. */
+  gal_linkedlist_free_dll(list);
+  return out;
+}
+
+
+
+
+
 
 
 
@@ -492,7 +600,7 @@ options_sanity_check(struct argp_option *option, char *arg,
 
 
     case GAL_OPTIONS_RANGE_GE_0_LE_1:
-      message="between zero and one";
+      message="between zero and one (inclusive)";
       ref1=gal_data_alloc(NULL, GAL_DATA_TYPE_UINT8, 1, &dsize, NULL,
                           0, -1, NULL, NULL, NULL);
       ref2=gal_data_alloc(NULL, GAL_DATA_TYPE_UINT8, 1, &dsize, NULL,
@@ -506,6 +614,21 @@ options_sanity_check(struct argp_option *option, char *arg,
       break;
 
 
+    case GAL_OPTIONS_RANGE_GT_0_LT_1:
+      message="between zero and one (not inclusive)";
+      ref1=gal_data_alloc(NULL, GAL_DATA_TYPE_UINT8, 1, &dsize, NULL,
+                          0, -1, NULL, NULL, NULL);
+      ref2=gal_data_alloc(NULL, GAL_DATA_TYPE_UINT8, 1, &dsize, NULL,
+                          0, -1, NULL, NULL, NULL);
+      *(unsigned char *)(ref1->array)=0;
+      *(unsigned char *)(ref2->array)=1;
+
+      operator1=GAL_ARITHMETIC_OP_GT;
+      operator2=GAL_ARITHMETIC_OP_LT;
+      multicheckop=GAL_ARITHMETIC_OP_AND;
+      break;
+
+
     case GAL_OPTIONS_RANGE_GT_0_ODD:
       message="greater than zero and odd";
       ref1=gal_data_alloc(NULL, GAL_DATA_TYPE_UINT8, 1, &dsize, NULL,
diff --git a/lib/options.h b/lib/options.h
index ff8c131..ebbe7bc 100644
--- a/lib/options.h
+++ b/lib/options.h
@@ -70,7 +70,7 @@ enum options_standard_groups
 /* Key values for the common options, the free alphabetical keys are listed
    below. You can use any of the free letters for an option in a
    program. Note that `-V', which is used by GNU and implemented by Argp,
-   is also removed.
+   is also removed from this list.
 
    a b c d e f g i j k l m n p r s t u v w x y z
    A B C E F G H J L M O Q R W X Y Z
@@ -114,6 +114,7 @@ enum gal_options_range_values
   GAL_OPTIONS_RANGE_GE_0,
   GAL_OPTIONS_RANGE_0_OR_1,
   GAL_OPTIONS_RANGE_GE_0_LE_1,
+  GAL_OPTIONS_RANGE_GT_0_LT_1,
 
   GAL_OPTIONS_RANGE_GT_0_ODD,
   GAL_OPTIONS_RANGE_0_OR_ODD,
@@ -208,6 +209,10 @@ gal_options_add_to_not_given(struct 
gal_options_common_params *cp,
 void
 gal_options_abort_if_mandatory_missing(struct gal_options_common_params *cp);
 
+gal_data_t *
+gal_options_parse_list_of_numbers(char *string, char *filename,
+                                  size_t lineno);
+
 
 
 
diff --git a/lib/qsort.c b/lib/qsort.c
index 889cd94..d02a892 100644
--- a/lib/qsort.c
+++ b/lib/qsort.c
@@ -23,6 +23,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <config.h>
 
 #include <stdlib.h>
+#include <stdint.h>
 
 #include <fitsio.h>
 
@@ -49,116 +50,104 @@ gal_qsort_index_float_decreasing(const void * a, const 
void * b)
 
 
 int
-gal_qsort_uchar_decreasing(const void * a, const void * b)
+gal_qsort_uint8_decreasing(const void *a, const void *b)
 {
-  return ( *(unsigned char*)b - *(unsigned char*)a );
+  return ( *(uint8_t *)b - *(uint8_t *)a );
 }
 
 int
-gal_qsort_uchar_increasing(const void * a, const void * b)
+gal_qsort_uint8_increasing(const void *a, const void *b)
 {
-  return ( *(unsigned char*)a - *(unsigned char*)b );
+  return ( *(uint8_t *)a - *(uint8_t *)b );
 }
 
 int
-gal_qsort_char_decreasing(const void * a, const void * b)
+gal_qsort_int8_decreasing(const void *a, const void *b)
 {
-  return ( *(char*)b - *(char*)a );
+  return ( *(int8_t *)b - *(int8_t *)a );
 }
 
 int
-gal_qsort_char_increasing(const void * a, const void * b)
+gal_qsort_int8_increasing(const void *a, const void *b)
 {
-  return ( *(char*)a - *(char*)b );
+  return ( *(int8_t *)a - *(int8_t *)b );
 }
 
 int
-gal_qsort_ushort_decreasing(const void * a, const void * b)
+gal_qsort_uint16_decreasing(const void *a, const void *b)
 {
-  return ( *(unsigned short*)b - *(unsigned short*)a );
+  return ( *(uint16_t *)b - *(uint16_t *)a );
 }
 
 int
-gal_qsort_ushort_increasing(const void * a, const void * b)
+gal_qsort_uint16_increasing(const void *a, const void *b)
 {
-  return ( *(unsigned short*)a - *(unsigned short*)b );
+  return ( *(uint16_t *)a - *(uint16_t *)b );
 }
 
 int
-gal_qsort_short_decreasing(const void * a, const void * b)
+gal_qsort_in16_decreasing(const void *a, const void *b)
 {
-  return ( *(short*)b - *(short*)a );
+  return ( *(int16_t *)b - *(int16_t *)a );
 }
 
 int
-gal_qsort_short_increasing(const void * a, const void * b)
+gal_qsort_int16_increasing(const void *a, const void *b)
 {
-  return ( *(short*)a - *(short*)b );
+  return ( *(int16_t *)a - *(int16_t *)b );
 }
 
 int
-gal_qsort_uint_decreasing(const void * a, const void * b)
+gal_qsort_uint32_decreasing(const void *a, const void *b)
 {
-  return ( *(unsigned int*)b - *(unsigned int*)a );
+  return ( *(uint32_t *)b - *(uint32_t *)a );
 }
 
 int
-gal_qsort_uint_increasing(const void * a, const void * b)
+gal_qsort_uint32_increasing(const void *a, const void *b)
 {
-  return ( *(unsigned int*)a - *(unsigned int*)b );
+  return ( *(uint32_t *)a - *(uint32_t *)b );
 }
 
 int
-gal_qsort_int_decreasing(const void * a, const void * b)
+gal_qsort_int32_decreasing(const void *a, const void *b)
 {
-  return ( *(int*)b - *(int*)a );
+  return ( *(int32_t *)b - *(int32_t *)a );
 }
 
 int
-gal_qsort_int_increasing(const void * a, const void * b)
+gal_qsort_int32_increasing(const void *a, const void *b)
 {
-  return ( *(int*)a - *(int*)b );
+  return ( *(int32_t *)a - *(int32_t *)b );
 }
 
 int
-gal_qsort_ulong_decreasing(const void * a, const void * b)
+gal_qsort_uint64_decreasing(const void *a, const void *b)
 {
-  return ( *(unsigned long*)b - *(unsigned long*)a );
+  return ( *(uint64_t *)b - *(uint64_t *)a );
 }
 
 int
-gal_qsort_ulong_increasing(const void * a, const void * b)
+gal_qsort_uint64_increasing(const void *a, const void *b)
 {
-  return ( *(unsigned long*)a - *(unsigned long*)b );
+  return ( *(uint64_t *)a - *(uint64_t *)b );
 }
 
 
 int
-gal_qsort_long_decreasing(const void * a, const void * b)
+gal_qsort_int64_decreasing(const void *a, const void *b)
 {
-  return ( *(long*)b - *(long*)a );
+  return ( *(int64_t *)b - *(int64_t *)a );
 }
 
 int
-gal_qsort_long_increasing(const void * a, const void * b)
+gal_qsort_int64_increasing(const void *a, const void *b)
 {
-  return ( *(long*)a - *(long*)b );
+  return ( *(int64_t *)a - *(int64_t *)b );
 }
 
 int
-gal_qsort_longlong_decreasing(const void * a, const void * b)
-{
-  return ( *(LONGLONG*)b - *(LONGLONG*)a );
-}
-
-int
-gal_qsort_longlong_increasing(const void * a, const void * b)
-{
-  return ( *(LONGLONG*)a - *(LONGLONG*)b );
-}
-
-int
-gal_qsort_float_decreasing(const void * a, const void * b)
+gal_qsort_float32_decreasing(const void *a, const void *b)
 {
   float ta=*(float*)a;
   float tb=*(float*)b;
@@ -166,7 +155,7 @@ gal_qsort_float_decreasing(const void * a, const void * b)
 }
 
 int
-gal_qsort_float_increasing(const void * a, const void * b)
+gal_qsort_float32_increasing(const void *a, const void *b)
 {
   float ta=*(float*)a;
   float tb=*(float*)b;
@@ -174,7 +163,7 @@ gal_qsort_float_increasing(const void * a, const void * b)
 }
 
 int
-gal_qsort_double_decreasing(const void * a, const void * b)
+gal_qsort_float64_decreasing(const void *a, const void *b)
 {
   double ta=*(double*)a;
   double tb=*(double*)b;
@@ -182,7 +171,7 @@ gal_qsort_double_decreasing(const void * a, const void * b)
 }
 
 int
-gal_qsort_double_increasing(const void * a, const void * b)
+gal_qsort_float64_increasing(const void *a, const void *b)
 {
   double ta=*(double*)a;
   double tb=*(double*)b;
diff --git a/lib/statistics.c b/lib/statistics.c
index 598346c..4d26261 100644
--- a/lib/statistics.c
+++ b/lib/statistics.c
@@ -30,8 +30,8 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <stdint.h>
 #include <stdlib.h>
 
+#include <gnuastro/data.h>
 #include <gnuastro/qsort.h>
-#include <gnuastro/array.h>
 #include <gnuastro/statistics.h>
 
 #include "mode.h"
@@ -41,770 +41,17 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 
 
-
-/****************************************************************
- *****************    Mininum and Maximum    ********************
- ****************************************************************/
-/* Find the the minimum (non-blank) value in an array of type long. Note
-   that the long type doesn't have a NaN value like the floats above. So as
-   blank pixels, a value in the range of acceptable values is given. so we
-   have to explicitly ignore those values.*/
-void
-gal_statistics_long_non_blank_min(long *in, size_t size, long *min,
-                                  long blank)
-{
-  long tmin=INT32_MAX, *lpt;
-  lpt=in+size;
-  do
-    if(*in!=blank && *in<tmin)
-      tmin=*in;
-  while(++in<lpt);
-  *min=tmin;
-}
-
-
-
-
-
-/* Find the the minimum (non-blank) value in an array of type long. See the
-   explanation above `gal_statistics_long_min'. */
-void
-gal_statistics_long_non_blank_max(long *in, size_t size, long *max,
-                                  long blank)
-{
-  long tmax=INT32_MIN, *lpt;
-  lpt=in+size;
-  do
-    if(*in!=blank && *in>tmax)
-      tmax=*in;
-  while(++in<lpt);
-  *max=tmax;
-}
-
-
-
-
-
-void
-gal_statistics_float_min(float *in, size_t size, float *min)
-{
-  float tmin=FLT_MAX, *fpt;
-  fpt=in+size;
-  do   /* Works for NAN, since NAN is not smaller than any number. */
-    if(*in<tmin) tmin=*in;
-  while(++in<fpt);
-  *min=tmin;
-}
-
-
-
-
-
-void
-gal_statistics_float_max(float *in, size_t size, float *max)
-{
-  float tmax=-FLT_MAX, *fpt;
-  fpt=in+size;
-  do /* Works for NAN, since NAN is not larger than any number. */
-    if(*in>tmax) tmax=*in;
-  while(++in<fpt);
-  *max=tmax;
-}
-
-
-
-
-
-void
-gal_statistics_double_min(double *in, size_t size, double *min)
-{
-  double tmin=FLT_MAX, *fpt;
-  fpt=in+size;
-  do   /* Works for NAN, since NAN is not smaller than any number. */
-    if(*in<tmin) tmin=*in;
-  while(++in<fpt);
-  *min=tmin;
-}
-
-
-
-
-
-double
-gal_statistics_double_min_return(double *in, size_t size)
-{
-  double tmin=FLT_MAX, *fpt;
-  fpt=in+size;
-  do   /* Works for NAN, since NAN is not smaller than any number. */
-    if(*in<tmin) tmin=*in;
-  while(++in<fpt);
-  return tmin;
-}
-
-
-
-
-
-void
-gal_statistics_double_max(double *in, size_t size, double *max)
-{
-  double tmax=-FLT_MAX, *fpt;
-  fpt=in+size;
-  do   /* Works for NAN, since NAN is not larger than any number. */
-    if(*in>tmax) tmax=*in;
-  while(++in<fpt);
-  *max=tmax;
-}
-
-
-
-
-
-double
-gal_statistics_double_max_return(double *in, size_t size)
-{
-  double tmax=-FLT_MAX, *fpt;
-  fpt=in+size;
-  do   /* Works for NAN, since NAN is not larger than any number. */
-    if(*in>tmax) tmax=*in;
-  while(++in<fpt);
-  return tmax;
-}
-
-
-
-
-void
-gal_statistics_float_max_masked(float *in, unsigned char *mask, size_t size,
-                                float *max)
-{
-  float tmax=-FLT_MAX, *fpt;
-  fpt=in+size;
-  do
-    if(*mask++==0 && *in>tmax)
-      tmax=*in;
-  while(++in<fpt);
-  *max=tmax;
-}
-
-
-
-
-
-void
-gal_statistics_float_second_max(float *in, size_t size, float *secondmax)
-{
-  float smax=-FLT_MAX, max=-FLT_MAX, *fpt;
-  fpt=in+size;
-  do
-    { /* Works for NAN, since NAN is not larger than any number. */
-      if(*in>max)
-        {
-          smax=max;
-          max=*in;
-        }
-      else if(*in>smax) smax=*in;
-    }
-  while(++in<fpt);
-  *secondmax=smax;
-}
-
-
-
-
-
-void
-gal_statistics_float_second_min(float *in, size_t size, float *secondmin)
-{
-  float smin=FLT_MAX, min=FLT_MAX, *fpt;
-  fpt=in+size;
-  do
-    { /* Works for NAN, since NAN is not smaller than any number. */
-      if(*in<min)
-        {
-          smin=min;
-          min=*in;
-        }
-      else if(*in<smin) smin=*in;
-    }
-  while(++in<fpt);
-  *secondmin=smin;
-}
-
-
-
-
-
-void
-gal_statistics_f_min_max(float *in, size_t size, float *min, float *max)
-{
-  float tmin=FLT_MAX, tmax=-FLT_MAX, *f, *fpt;
-
-  fpt=(f=in)+size;
-  do
-    {   /* Works for NAN, because NaN values are not greater or
-           smaller than any number. */
-      if (*f>tmax) tmax=*f;
-      if (*f<tmin) tmin=*f;
-    }
-  while(++f<fpt);
-
-  /* If the whole data was a NaN, then tmin and tmax did not change
-     from their initial values. */
-  if(tmin==FLT_MAX || tmax==-FLT_MAX)
-    *min=*max=NAN;
-  else
-    {
-      *max=tmax;
-      *min=tmin;
-    }
-}
-
-
-
-
-
-void
-gal_statistics_d_min_max(double *in, size_t size, double *min, double *max)
-{
-  double tmin=FLT_MAX, tmax=-FLT_MAX, *d, *dpt;
-
-  dpt=(d=in)+size;
-  do
-    { /* Works for NAN */
-      if (*d>tmax) tmax=*d;
-      if (*d<tmin) tmin=*d;
-    }
-  while(++d<dpt);
-
-  /* If all the data was a NaN */
-  if(tmin==FLT_MAX || tmax==-FLT_MAX)
-    *min=*max=NAN;
-  else
-    {
-      *max=tmax;
-      *min=tmin;
-    }
-}
-
-
-
-
-void
-gal_statistics_d_max_with_index(double *in, size_t size, double *max,
-                                size_t *index)
-{
-  size_t tindex=0;
-  double *fpt, *pt=in, tmax=-FLT_MAX;
-
-  fpt=pt+size;
-  do  /*  Works for NAN, see comments above. */
-    if(*pt>tmax)
-      {
-        tmax=*pt;
-        tindex=pt-in;
-      }
-  while(++pt<fpt);
-  *index=tindex;
-  *max=tmax;
-}
-
-
-
-
-
-void
-gal_statistics_f_max_with_index(float *in, size_t size,
-                                float *max, size_t *index)
-{
-  size_t tindex=0;
-  float *pt=in, *fpt, tmax=-FLT_MAX;
-
-  fpt=pt+size;
-  do  /* Works for NAN, see comments above.x */
-    if(*pt>tmax)
-      {
-        tmax=*pt;
-        tindex=pt-in;
-      }
-  while(++pt<fpt);
-  *index=tindex;
-  *max=tmax;
-}
-
-
-
-
-
-void
-gal_statistics_d_min_with_index(double *in, size_t size,
-                                double *min, size_t *index)
-{
-  size_t tindex=0;
-  double *pt=in, *fpt, tmin=FLT_MAX;
-
-  fpt=pt+size;
-  do  /* Works for NAN, see comments above. */
-    if(*pt<tmin)
-      {
-        tmin=*pt;
-        tindex=pt-in;
-      }
-  while(++pt<fpt);
-  *index=tindex;
-  *min=tmin;
-}
-
-
-
-
-
-void
-gal_statistics_f_min_with_index(float *in, size_t size,
-                                float *min, size_t *index)
-{
-  size_t tindex=0;
-  float *pt=in, *fpt, tmin=FLT_MAX;
-
-  fpt=pt+size;
-  do /* Works for NAN, see comments above. */
-    if(*pt<tmin)
-      {
-        tmin=*pt;
-        tindex=pt-in;
-      }
-  while(++pt<fpt);
-  *index=tindex;
-  *min=tmin;
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-/****************************************************************
- *****************            Sum            ********************
- ****************************************************************/
-float
-gal_statistics_float_sum(float *in, size_t size)
-{
-  float *fpt;
-  double sum=0;
-  fpt=in+size;
-  do
-    if(!isnan(*in))
-      sum+=*in;
-  while(++in<fpt);
-  return sum;
-}
-
-
-
-
-
-float
-gal_statistics_float_sum_num(float *in, size_t *size)
-{
-  float *fpt;
-  double sum=0;
-  fpt=in+*size;
-  *size=0;
-  do
-    if(!isnan(*in))
-      {
-        sum+=*in;
-        ++(*size);
-      }
-  while(++in<fpt);
-  return sum;
-}
-
-
-
-
-
-static double
-doublesumnum(double *in, size_t *size)
-{
-  double *fpt;
-  double sum=0;
-
-  /* If size is initially zero, then return 0. */
-  if(*size==0) return 0.0f;
-
-  /* Go through all the array: */
-  fpt=in+*size;
-  *size=0;
-  do
-    if(!isnan(*in))
-      {
-        sum+=*in;
-        ++(*size);
-      }
-  while(++in<fpt);
-
-  /* If the size was not initially zero, but after going through the
-     whole array, it is still zero, then the whole array had NaN
-     values. */
-  return *size ? sum : NAN;
-}
-
-
-
-
-float
-gal_statistics_float_sum_squared(float *in, size_t size)
-{
-  float *fpt;
-  double sum=0;
-  fpt=in+size;
-  do
-    if(!isnan(*in))
-      sum+=*in * *in;
-  while(++in<fpt);
-  return sum;
-}
-
-
-
-
-
-/* Sum over all elements of the array that are not covered by a
-   mask. Any non-zero masked pixel is considered to be a masked
-   pixel. */
-float
-gal_statistics_float_sum_mask(float *in, unsigned char *mask,
-                              size_t size, size_t *nsize)
-{
-  double sum=0;
-  size_t counter=0;
-  unsigned char *pt=mask, *fpt;
-
-  fpt=pt+size;
-  do
-    if(*pt==0)
-      {
-        sum+=in[pt-mask];
-        ++counter;
-      }
-  while(++pt<fpt);
-
-  *nsize=counter;
-  return sum;
-}
-
-
-
-
-
-float
-gal_statistics_float_sum_mask_l(float *in, long *mask,
-                                size_t size, size_t *nsize)
-{
-  double sum=0;
-  size_t counter=0;
-  long *pt=mask, *fpt;
-
-  fpt=pt+size;
-  do
-    if(*pt==0)
-      {
-        sum+=in[pt-mask];
-        ++counter;
-      }
-  while(++pt<fpt);
-
-  *nsize=counter;
-  return sum;
-}
-
-
-
-
-
-float
-gal_statistics_float_sum_squared_mask(float *in, unsigned char *mask,
-                                      size_t size, size_t *nsize)
-{
-  double sum=0;
-  size_t counter=0;
-  unsigned char *pt=mask, *fpt;
-
-  fpt=pt+size;
-  do
-    if(*pt==0)
-      {
-        sum+=in[pt-mask] * in[pt-mask];
-        ++counter;
-      }
-  while(++pt<fpt);
-
-  *nsize=counter;
-  return sum;
-}
-
-
-
-
-
-
-float
-gal_statistics_float_sum_squared_mask_l(float *in, long *mask,
-                                        size_t size, size_t *nsize)
-{
-  double sum=0;
-  size_t counter=0;
-  long *pt=mask, *fpt;
-
-  fpt=pt+size;
-  do
-    if(*pt==0)
-      {
-        sum+=in[pt-mask] * in[pt-mask];
-        ++counter;
-      }
-  while(++pt<fpt);
-
-  *nsize=counter;
-  return sum;
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
 /****************************************************************
- *****************      Average and          ********************
- ****************    Standard deviation      ********************
+ ********                      Sort                       *******
  ****************************************************************/
-float
-gal_statistics_float_average(float *in, size_t size)
-{
-  float sum;
-  sum=gal_statistics_float_sum_num(in, &size);
-  return sum/size;
-}
-
-
-
-
-
-double
-gal_statistics_double_average(double *in, size_t size)
-{
-  double sum;
-  sum=doublesumnum(in, &size);
-  return sum/size;
-}
-
-
-
-
-
 void
-gal_statistics_f_ave_l(float *in, size_t size, float *ave, long *mask)
+gal_statistics_sort()
 {
-  float sum;
-  size_t nsize;
-  if(mask==NULL)
-    sum=gal_statistics_float_sum(in, size);
-  else
-    {
-      sum=gal_statistics_float_sum_mask_l(in, mask, size, &nsize);
-      size=nsize;
-    }
-  *ave=sum/size;
-}
 
-
-
-
-
-/* Find the average and standard deviation of an array, assuming that
-   there is a mask array. Any mask array pixel that is not zero will
-   not be included in the average and standard deviations.  Here the
-   mask is assumed to be unsigned char.  */
-void
-gal_statistics_f_ave_std(float *in, size_t size, float *ave,
-                         float *std, unsigned char *mask)
-{
-  size_t nsize1, nsize2;
-  float sum, sum2;
-  if(mask)
-    {
-      sum=gal_statistics_float_sum_mask(in, mask, size, &nsize1);
-      sum2=gal_statistics_float_sum_squared_mask(in, mask, size, &nsize2);
-      if(nsize1!=nsize2)
-        error(EXIT_FAILURE, 0, "a bug in gal_statistics_f_ave_std "
-              "(lib/statistics.h).  Somehow the number of masked pixels is "
-              "measured differently.  Please contact us so we can find the "
-              "cause");
-      size=nsize1;
-    }
-  else
-    {
-      sum=gal_statistics_float_sum(in, size);
-      sum2=gal_statistics_float_sum_squared(in, size);
-    }
-  *ave=sum/size;
-  *std=sqrt( (sum2-sum*sum/size)/size );
 }
 
 
 
-
-
-/* Similar to gal_statistics_f_ave_std, but when the mask is assumed to be a
-   long array.  */
-void
-gal_statistics_f_ave_std_l(float *in, size_t size, float *ave,
-                           float *std, long *mask)
-{
-  size_t nsize1, nsize2;
-  float sum, sum2;
-  if(mask==NULL)
-    {
-      sum=gal_statistics_float_sum(in, size);
-      sum2=gal_statistics_float_sum_squared(in, size);
-    }
-  else
-    {
-      sum=gal_statistics_float_sum_mask_l(in, mask, size, &nsize1);
-      sum2=gal_statistics_float_sum_squared_mask_l(in, mask, size, &nsize2);
-      if(nsize1!=nsize2)
-        error(EXIT_FAILURE, 0, "a bug in favestl (lib/statistics.h). "
-              "Somehow the number of masked pixels is measured "
-              "differently. Please contact us so we can find the cause");
-      size=nsize1;
-    }
-  *ave=sum/size;
-  *std=sqrt( (sum2-sum*sum/size)/size );
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-/****************************************************************
- *****************           Median            ******************
- ****************************************************************/
-float
-gal_statistics_median(float *array, size_t insize)
-{
-  float *copy, median;
-  size_t size=insize, medind;
-
-  /* Make a copy of the input, shift all its non-NaN elements to the
-     start of the array, then sort them and find the median. */
-  gal_array_float_copy(array, insize, &copy);
-  gal_array_no_nans(copy, &size);
-  medind=gal_statistics_index_from_quantile(size, 0.5);
-  qsort(copy, size, sizeof*copy, gal_qsort_float_increasing);
-  median=copy[medind];
-
-  /* Clean up */
-  free(copy);
-  return median;
-}
-
-
-
-
-
-double
-gal_statistics_median_double_in_place(double *array, size_t insize)
-{
-  size_t size=insize, medind;
-
-  /* Shift all its non-NaN elements to the start of the array, then
-     sort them and find the median. */
-  gal_array_no_nans_double(array, &size);
-
-  /* If all the elements are NaN (size==0), then return NaN,
-     otherwise, find the median. */
-  if(size)
-    {
-      qsort(array, size, sizeof*array, gal_qsort_double_increasing);
-      medind=gal_statistics_index_from_quantile(size, 0.5);
-
-      if(size%2)
-        return array[medind];
-      else
-        return (array[medind]+array[medind-1])/2;
-    }
-  else
-    return NAN;
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
 /****************************************************************
  ********     Histogram and Cumulative Frequency Plot     *******
  ****************************************************************/
@@ -1134,9 +381,12 @@ gal_statistics_index_from_quantile(size_t size, float 
quant)
 int
 gal_statistics_sigma_clip_converge(float *array, int o1_n0, size_t num_elem,
                                    float sigma_multiple, float accuracy,
-                                   float *outave, float *outmed, float *outstd,
-                                   int print)
+                                   float *outave, float *outmed,
+                                   float *outstd, int print)
 {
+  printf("\n ... in gal_statistics_sigma_clip_converge ... \n");
+  exit(1);
+#if 0
   size_t counter=0;
   float *start, *oldstart, *dpt;
   float ave=*outave, med=*outmed, std=*outstd;
@@ -1146,7 +396,7 @@ gal_statistics_sigma_clip_converge(float *array, int 
o1_n0, size_t num_elem,
     {
       gal_array_float_copy(array, num_elem, &orderedarray);
       qsort(orderedarray, num_elem, sizeof*orderedarray,
-            gal_qsort_float_increasing);
+            gal_qsort_float32_increasing);
     }
   else orderedarray=array;
 
@@ -1196,6 +446,7 @@ gal_statistics_sigma_clip_converge(float *array, int 
o1_n0, size_t num_elem,
       oldstd=std;
       ++counter;
     }
+#endif
   return 0;
 }
 
@@ -1212,6 +463,9 @@ gal_statistics_sigma_clip_certain_num(float *array, int 
o1_n0, size_t num_elem,
                                       float *outave, float *outmed,
                                       float *outstd, int print)
 {
+  printf("\n ... in gal_statistics_sigma_clip_certain_num ... \n");
+  exit(1);
+#if 0
   size_t counter=0;
   float ave=*outave, med=*outmed, std=*outstd;
   float *start, *oldstart, *dpt, *orderedarray;
@@ -1220,7 +474,7 @@ gal_statistics_sigma_clip_certain_num(float *array, int 
o1_n0, size_t num_elem,
     {
       gal_array_float_copy(array, num_elem, &orderedarray);
       qsort(orderedarray, num_elem, sizeof*orderedarray,
-            gal_qsort_float_increasing);
+            gal_qsort_float32_increasing);
     }
   else orderedarray=array;
 
@@ -1265,6 +519,7 @@ gal_statistics_sigma_clip_certain_num(float *array, int 
o1_n0, size_t num_elem,
   *outave=ave;
   *outmed=med;
   *outstd=std;
+#endif
   return 1;
 }
 
@@ -1294,6 +549,9 @@ gal_statistics_sigma_clip_certain_num(float *array, int 
o1_n0, size_t num_elem,
 void
 gal_statistics_remove_outliers_flat_cdf(float *sorted, size_t *outsize)
 {
+  printf("\n ... in gal_statistics_remove_outliers_flat_cdf ... \n");
+  exit(1);
+#if 0
   int firstfound=0;
   size_t size=*outsize, i, maxind;
   float *slopes, minslope, maxslope;
@@ -1343,6 +601,7 @@ gal_statistics_remove_outliers_flat_cdf(float *sorted, 
size_t *outsize)
   */
 
   free(slopes);
+#endif
 }
 
 
@@ -1373,6 +632,9 @@ gal_statistics_mode_mirror_plots(float *sorted, size_t 
size, size_t mirrorindex,
                                  char *histsname, char *cfpsname,
                                  float mirrorplotdist)
 {
+  printf("\n... in gal_statistics_mode_mirror_plots ...\n");
+  exit(1);
+#if 0
   FILE *fp;
   size_t i, msize;
   float *out, maxhist=-FLT_MAX, maxcfp, d;
@@ -1490,6 +752,7 @@ gal_statistics_mode_mirror_plots(float *sorted, size_t 
size, size_t mirrorindex,
   free(bins);
   free(mirror);
   free(actual);
+#endif
 }
 
 
diff --git a/lib/table.c b/lib/table.c
index 6da9e3c..597880c 100644
--- a/lib/table.c
+++ b/lib/table.c
@@ -30,6 +30,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <string.h>
 
 #include <gnuastro/txt.h>
+#include <gnuastro/blank.h>
 #include <gnuastro/table.h>
 
 #include <checkset.h>
@@ -39,9 +40,56 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 
 
+/************************************************************************/
+/***************              Error messages              ***************/
+/************************************************************************/
+void
+gal_table_error_col_selection(char *filename, char *hdu, char *errorstring)
+{
+  char *c, *name, *command;
+
+  /* Set the proper pointers. */
+  if(gal_fits_name_is_fits(filename))
+    {
+      asprintf(&name, "%s (hdu: %s)", filename, hdu);
+      c=hdu; while(*c!='\0') if(isspace(*c++)) break;
+      asprintf(&command, *c=='\0' ? "%s --hdu=%s" : "%s --hdu=\"%s\"",
+               filename, hdu);
+    }
+  else command=name=filename;
+
+  /* Abort with with the proper error. */
+  error(EXIT_FAILURE, 0, "%s: %s\n\nFor more information on selecting "
+        "columns in Gnuastro, please run the following command (press "
+        "`SPACE' to go down and `q' to return to the command-line):\n\n"
+        "    $ info gnuastro \"Selecting table columns\"\n\n"
+        "To define a better column selection criteria, you can see "
+        "the list of column meta-data in this table, with the following "
+        "command:\n\n"
+        "    $ asttable %s --info\n", name, errorstring, command);
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
 
 /************************************************************************/
-/***************               Table types                ***************/
+/***************                 Formats                  ***************/
 /************************************************************************/
 /* Return the type of desired table based on a standard string. */
 uint8_t
@@ -152,30 +200,6 @@ gal_table_check_fits_format(char *filename, int 
tableformat)
 
 
 
-/* Programs that read columns might match more than one column for a given
-   property, for example it might happen (due to bad table design) that
-   multiple columns have the same name, or the user uses regular
-   expressions badly. In such cases, the following error message will be
-   useful and commonly repeated by the programs, so we are putting it here
-   for easiness and a clean code. */
-void
-gal_table_too_many_columns(char *filename)
-{
-  error(EXIT_FAILURE, 0, "there was more than one match for at least one of "
-        "the input columns from `%s'. You can check the column information "
-        "with the following command, and correct the options ending with "
-        "`col' appropriately.\n\n"
-        "   $ asttable --information %s\n\n"
-        "The current values to all options can be checked by calling the "
-        "`--printparams' (or `-P') option. To learn more about how the "
-        "columns are selected, please try the following command:\n\n"
-        "   $ info gnuastro \"Selecting table columns\" ", filename,
-        filename);
-}
-
-
-
-
 
 
 
@@ -235,10 +259,10 @@ gal_table_print_info(gal_data_t *allcols, size_t numcols, 
size_t numrows)
       unit    = allcols[i].unit;       /* readability. The compiler   */
       comment = allcols[i].comment;    /* optimizer will remove them. */
       printf("%-*zu%-*s%-*s%-*s%s\n", Nw, i+1,
-             nw, name ? name : GAL_DATA_BLANK_STRING ,
-             uw, unit ? unit : GAL_DATA_BLANK_STRING ,
+             nw, name ? name : GAL_BLANK_STRING ,
+             uw, unit ? unit : GAL_BLANK_STRING ,
              tw, gal_data_type_as_string(allcols[i].type, 1),
-             comment ? comment : GAL_DATA_BLANK_STRING);
+             comment ? comment : GAL_BLANK_STRING);
     }
 
   /* Print the number of rows. */
@@ -522,7 +546,7 @@ gal_table_info(char *filename, char *hdu, size_t *numcols, 
size_t *numrows,
 
 /* Function to print regular expression error. This is taken from the GNU C
    library manual, with small modifications to fit out style, */
-void
+static void
 regexerrorexit(int errcode, regex_t *compiled, char *input)
 {
   char *regexerrbuf;
@@ -580,24 +604,25 @@ make_list_of_indexs(struct gal_linkedlist_stll *cols, 
gal_data_t *allcols,
   long tlong;
   int regreturn;
   regex_t *regex;
-  size_t i, nummatch;
+  size_t i, nummatch, len;
   struct gal_linkedlist_stll *tmp;
-  char *str, *colts, *strcheck, *tailptr;
   struct gal_linkedlist_sll *indexll=NULL;
+  char *str, *strcheck, *tailptr, *errorstring;
 
   for(tmp=cols; tmp!=NULL; tmp=tmp->next)
     {
-      /* Counter for number of columns matched. */
+      /* Counter for number of columns matched, and length of name. */
       nummatch=0;
+      len=strlen(tmp->v);
 
       /* REGULAR EXPRESSION: When the first and last characters are `/'. */
-      if( tmp->v[0]=='/' && tmp->v[strlen(tmp->v)-1]=='/' )
+      if( tmp->v[0]=='/' && tmp->v[len-1]=='/' )
         {
           /* Remove the slashes, note that we don't want to change `tmp->v'
              (because it should be freed later). So first we set the last
              character to `\0', then define a new string from the first
              element. */
-          tmp->v[strlen(tmp->v)-1]='\0';
+          tmp->v[len-1]='\0';
           str = tmp->v + 1;
 
           /* Allocate the regex_t structure: */
@@ -642,6 +667,11 @@ make_list_of_indexs(struct gal_linkedlist_stll *cols, 
gal_data_t *allcols,
 
           /* Free the regex_t structure: */
           regfree(regex);
+
+          /* Put the `/' back into the input string. This is done because
+             after this function, the calling program might want to inform
+             the user of their exact input string. */
+          tmp->v[len-1]='/';
         }
 
 
@@ -714,24 +744,12 @@ make_list_of_indexs(struct gal_linkedlist_stll *cols, 
gal_data_t *allcols,
 
       /* If there was no match, then report an error. This can only happen
          for string matches, not column numbers, for numbers, the checks
-         are done before the reading.*/
+         are done (and program is aborted) before this step. */
       if(nummatch==0)
         {
-          colts = ( searchin==GAL_TABLE_SEARCH_NAME ? "names"
-                    : ( searchin==GAL_TABLE_SEARCH_UNIT ? "units"
-                        : "comments") );
-          error(EXIT_FAILURE, 0, "`%s' didn't match with any of the column "
-                "%s in `%s'%s%s%s. You can check the available column "
-                "information by running `asttable %s%s%s --information'. "
-                "%s", tmp->v, colts, filename,
-                gal_fits_name_is_fits(filename) ? " (hdu: " : "",
-                gal_fits_name_is_fits(filename) ? hdu : "",
-                gal_fits_name_is_fits(filename) ? ")" : "",
-                filename,
-                gal_fits_name_is_fits(filename) ? " --hdu" : "",
-                gal_fits_name_is_fits(filename) ? hdu : "",
-                ignorecase ? "" : "For a case-insensitive match, "
-                "run again with the `--ignorecase' option. " );
+          asprintf(&errorstring, "`%s' didn't match any of the column %ss.",
+                   tmp->v, gal_table_searchin_as_string(searchin));
+          gal_table_error_col_selection(filename, hdu, errorstring);
         }
     }
 
diff --git a/lib/txt.c b/lib/txt.c
index 8c93cb4..684972d 100644
--- a/lib/txt.c
+++ b/lib/txt.c
@@ -31,6 +31,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <stdlib.h>
 
 #include <gnuastro/txt.h>
+#include <gnuastro/blank.h>
 #include <gnuastro/table.h>
 
 #include <checkset.h>
@@ -665,56 +666,56 @@ txt_read_token(gal_data_t *data, gal_data_t *info, char 
*token,
       if( (strb=info->array) && !strcmp( *strb, str[i] ) )
         {
           free(str[i]);
-          gal_checkset_allocate_copy(GAL_DATA_BLANK_STRING, &str[i]);
+          gal_checkset_allocate_copy(GAL_BLANK_STRING, &str[i]);
         }
       break;
 
     case GAL_DATA_TYPE_UINT8:
       uc[i]=strtol(token, &tailptr, 0);
       if( (ucb=info->array) && *ucb==uc[i] )
-        uc[i]=GAL_DATA_BLANK_UINT8;
+        uc[i]=GAL_BLANK_UINT8;
       break;
 
     case GAL_DATA_TYPE_INT8:
       c[i]=strtol(token, &tailptr, 0);
       if( (cb=info->array) && *cb==c[i] )
-        c[i]=GAL_DATA_BLANK_INT8;
+        c[i]=GAL_BLANK_INT8;
       break;
 
     case GAL_DATA_TYPE_UINT16:
       us[i]=strtol(token, &tailptr, 0);
       if( (usb=info->array) && *usb==us[i] )
-        us[i]=GAL_DATA_BLANK_UINT16;
+        us[i]=GAL_BLANK_UINT16;
       break;
 
     case GAL_DATA_TYPE_INT16:
       s[i]=strtol(token, &tailptr, 0);
       if( (sb=info->array) && *sb==s[i] )
-        s[i]=GAL_DATA_BLANK_INT16;
+        s[i]=GAL_BLANK_INT16;
       break;
 
     case GAL_DATA_TYPE_UINT32:
       ui[i]=strtol(token, &tailptr, 0);
       if( (uib=info->array) && *uib==ui[i] )
-        ui[i]=GAL_DATA_BLANK_UINT32;
+        ui[i]=GAL_BLANK_UINT32;
       break;
 
     case GAL_DATA_TYPE_INT32:
       ii[i]=strtol(token, &tailptr, 0);
       if( (ib=info->array) && *ib==ii[i] )
-        ii[i]=GAL_DATA_BLANK_INT32;
+        ii[i]=GAL_BLANK_INT32;
       break;
 
     case GAL_DATA_TYPE_UINT64:
       ul[i]=strtoul(token, &tailptr, 0);
       if( (ulb=info->array) && *ulb==ul[i] )
-        ul[i]=GAL_DATA_BLANK_UINT64;
+        ul[i]=GAL_BLANK_UINT64;
       break;
 
     case GAL_DATA_TYPE_INT64:
       l[i]=strtol(token, &tailptr, 0);
       if( (lb=info->array) && *lb==l[i] )
-        l[i]=GAL_DATA_BLANK_INT64;
+        l[i]=GAL_BLANK_INT64;
       break;
 
       /* For the blank value of floating point types, we need to make
@@ -725,14 +726,14 @@ txt_read_token(gal_data_t *data, gal_data_t *info, char 
*token,
       f[i]=strtod(token, &tailptr);
       if( (fb=info->array)
           && ( (isnan(*fb) && isnan(f[i])) || *fb==f[i] ) )
-        f[i]=GAL_DATA_BLANK_FLOAT64;
+        f[i]=GAL_BLANK_FLOAT64;
       break;
 
     case GAL_DATA_TYPE_FLOAT64:
       d[i]=strtod(token, &tailptr);
       if( (db=info->array)
           && ( (isnan(*db) && isnan(d[i])) || *db==d[i] ) )
-        d[i]=GAL_DATA_BLANK_FLOAT64;
+        d[i]=GAL_BLANK_FLOAT64;
       break;
 
     default:
@@ -1039,16 +1040,16 @@ make_fmts_for_printf(gal_data_t *datall, int 
leftadjust, size_t *len)
 
       /* If we have a blank value, get the blank value as a string and
          adjust the width */
-      if(gal_data_has_blank(data)==0)
+      if(gal_blank_present(data)==0)
         fmts[i*FMTS_COLS+2]=NULL;
       else
         {
           /* Set the blank value. */
           if(data->type==GAL_DATA_TYPE_STRING)
-            gal_checkset_allocate_copy(GAL_DATA_BLANK_STRING,
+            gal_checkset_allocate_copy(GAL_BLANK_STRING,
                                        &fmts[i*FMTS_COLS+2]);
           else
-            fmts[i*FMTS_COLS+2]=gal_data_blank_as_string(data->type, 0);
+            fmts[i*FMTS_COLS+2]=gal_blank_as_string(data->type, 0);
         }
 
 
@@ -1146,8 +1147,8 @@ txt_print_value(FILE *fp, void *array, int type, size_t 
ind, char *fmt)
 
       /* Special consideration for strings. */
     case GAL_DATA_TYPE_STRING:
-      if( !strcmp( ((char **)array)[ind], GAL_DATA_BLANK_STRING ) )
-        fprintf(fp, fmt, GAL_DATA_BLANK_STRING);
+      if( !strcmp( ((char **)array)[ind], GAL_BLANK_STRING ) )
+        fprintf(fp, fmt, GAL_BLANK_STRING);
       else
         fprintf(fp, fmt, ((char **)array)[ind]);
       break;
diff --git a/tests/Makefile.am b/tests/Makefile.am
index f8980fa..4f000b4 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -56,8 +56,7 @@ export hasghostscript=$(MAYBE_HASGHOSTSCRIPT);
 
 # Compilations that are to be done with `make check'.
 LDADD = -lgnuastro
-check_PROGRAMS = versioncpp arraymanip multithread
-arraymanip_SOURCES = lib/arraymanip.c
+check_PROGRAMS = versioncpp multithread
 versioncpp_SOURCES = lib/versioncpp.cpp
 multithread_SOURCES = lib/multithread.c
 
diff --git a/tests/lib/arraymanip.c b/tests/lib/arraymanip.c
deleted file mode 100644
index 114e3f4..0000000
--- a/tests/lib/arraymanip.c
+++ /dev/null
@@ -1,75 +0,0 @@
-/*********************************************************************
-A test program to test array functions in Gnuastro.
-
-Original author:
-     Mohammad Akhlaghi <address@hidden>
-Contributing author(s):
-Copyright (C) 2016, 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 <stdio.h>
-#include <stdlib.h>
-
-#include "gnuastro/array.h"
-#include "gnuastro/config.h"
-
-int
-main(void)
-{
-  double *array;
-  size_t i, size=6;
-
-  /* Allocate the array, report an error if it wasn't allocated. */
-  array=malloc(size * sizeof *array);
-  if(array==NULL)
-    {
-      fprintf(stderr, "%zu bytes for d.\n", size);
-      exit(EXIT_FAILURE);
-    }
-
-  /* Print the version of Gnuastro being used: */
-  printf("Test of Gnuastro %s\n", GAL_CONFIG_VERSION);
-
-  /* Fill in the test array and report its contents at the same time. */
-  printf("Input array: ");
-  for(i=0;i<size;++i)
-    {
-      array[i]=i;
-      printf("%.3f, ", array[i]);
-    }
-
-  /* delete the last two characters, add a . and newline */
-  printf("\b\b.\n");
-
-  /* Now do some simple operations, and report each */
-  printf("Multipyling by 2.0...\n");
-  gal_array_dmultip_const(array, size, 2.0f);
-
-  printf("Dividing by 3.0...\n");
-  gal_array_ddivide_const(array, size, 3.0f);
-
-  printf("Taking natural logarithm...\n");
-  gal_array_dlog_array(array, size);
-
-  /* Report the final array: */
-  printf("Output array: ");
-  for(i=0;i<size;++i)
-    printf("%.3f, ", array[i]);
-  printf("\b\b.\n");
-
-  /* Cleanup and return */
-  free(array);
-  return EXIT_SUCCESS;
-}
diff --git a/tmpfs-config-make b/tmpfs-config-make
index 09579f4..8dde40a 100755
--- a/tmpfs-config-make
+++ b/tmpfs-config-make
@@ -130,11 +130,11 @@ cd $build_dir
 #
 # ####################################
 if [ ! -f Makefile ]; then
-    $srcdir/configure --srcdir=$srcdir --disable-shared CFLAGS="-g -O0"      \
-                      --enable-arithmetic --enable-convertt --enable-convolve\
-                      --enable-cosmiccal  --enable-crop     --enable-header  \
-                      --enable-mknoise    --enable-mkprof   --enable-table   \
-                      --enable-warp
+    $srcdir/configure --srcdir=$srcdir --disable-shared CFLAGS="-g -O0"        
\
+                      --enable-arithmetic --enable-convertt   
--enable-convolve\
+                      --enable-cosmiccal  --enable-crop       --enable-header  
\
+                      --enable-mknoise    --enable-statistics --enable-mkprof  
\
+                      --enable-table      --enable-warp
 fi
 
 



reply via email to

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