[gnuastro-commits] master 5831e9e 052/113: Recent work in master imported, no conflicts

[Mohammad Akhlaghi]

branch: master
commit 5831e9ecc734b3b361e40fab7ea6e8371b79de4d
Merge: 1b0046b 27a61e0
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    Recent work in master imported, no conflicts
    
    Recent features and bug-fixes in the main branch are now imported here with
    no conflicts.
---
 NEWS                        |   7 ++
 bin/match/ui.c              |  13 +++
 bin/mkcatalog/args.h        |  28 ++++++
 bin/mkcatalog/columns.c     |  66 +++++++++++++
 bin/mkcatalog/main.h        |   6 +-
 bin/mkcatalog/mkcatalog.c   | 126 ++++++++++++++++++++++-
 bin/mkcatalog/ui.h          |   2 +
 bin/statistics/statistics.c |   2 -
 bin/statistics/ui.c         |   7 +-
 doc/gnuastro.texi           | 236 ++++++++++++++++++++++++++------------------
 lib/statistics.c            |  10 +-
 11 files changed, 393 insertions(+), 110 deletions(-)

diff --git a/NEWS b/NEWS
index 19c3ca3..7e83d09 100644
--- a/NEWS
+++ b/NEWS
@@ -9,6 +9,9 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
   cumulative frequency plots to be set outside the minimum or maximum
   values of the dataset.
 
+  MakeCatalog: The new `--mean' and `--median' options will calculate the
+  mean and median pixel value within an object or clump respectively.
+
 ** Removed features
 
 ** Changed features
@@ -23,6 +26,10 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
 
   NULL value of onlyversion option causing a crash (bug #53147).
 
+  Match output directory ignored when making multiple files (bug #53226).
+
+  Statistics program bad results on integer columns with limits (bug #53230).
+
 
 
 
diff --git a/bin/match/ui.c b/bin/match/ui.c
index 8cb6942..c4f4123 100644
--- a/bin/match/ui.c
+++ b/bin/match/ui.c
@@ -691,6 +691,10 @@ ui_preparations_out_cols(struct matchparams *p)
 static void
 ui_preparations_out_name(struct matchparams *p)
 {
+  /* To temporarily keep the original value. */
+  uint8_t keepinputdir_orig;
+
+  /* Set the output file(s) name(s). */
   if(p->logasoutput)
     {
       /* Set the logname (as output). */
@@ -731,12 +735,21 @@ ui_preparations_out_name(struct matchparams *p)
                 }
               else
                 {
+                  /* Here, we are be using the output name as input to the
+                     automatic output generating function (usually it is
+                     the input name, not the output name). Therefore, the
+                     `keepinputdir' variable should be 1. So we will
+                     temporarily change it here, then set it back to what
+                     it was. */
+                  keepinputdir_orig=p->cp.keepinputdir;
+                  p->cp.keepinputdir=1;
                   p->out1name=gal_checkset_automatic_output(&p->cp,
                                                             p->cp.output,
                                                             "_matched_1.txt");
                   p->out2name=gal_checkset_automatic_output(&p->cp,
                                                             p->cp.output,
                                                             "_matched_2.txt");
+                  p->cp.keepinputdir=keepinputdir_orig;
                 }
             }
           else
diff --git a/bin/mkcatalog/args.h b/bin/mkcatalog/args.h
index 288aa51..5602b8a 100644
--- a/bin/mkcatalog/args.h
+++ b/bin/mkcatalog/args.h
@@ -815,6 +815,34 @@ struct argp_option program_options[] =
       ui_column_codes_ll
     },
     {
+      "mean",
+      UI_KEY_MEAN,
+      0,
+      0,
+      "Mean of values in object/clump.",
+      UI_GROUP_COLUMNS_BRIGHTNESS,
+      0,
+      GAL_TYPE_INVALID,
+      GAL_OPTIONS_RANGE_ANY,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET,
+      ui_column_codes_ll
+    },
+    {
+      "median",
+      UI_KEY_MEDIAN,
+      0,
+      0,
+      "Median of values in object/clump.",
+      UI_GROUP_COLUMNS_BRIGHTNESS,
+      0,
+      GAL_TYPE_INVALID,
+      GAL_OPTIONS_RANGE_ANY,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET,
+      ui_column_codes_ll
+    },
+    {
       "magnitude",
       UI_KEY_MAGNITUDE,
       0,
diff --git a/bin/mkcatalog/columns.c b/bin/mkcatalog/columns.c
index b8061cb..035acae 100644
--- a/bin/mkcatalog/columns.c
+++ b/bin/mkcatalog/columns.c
@@ -855,6 +855,42 @@ columns_define_alloc(struct mkcatalogparams *p)
           ciflag[ CCOL_SUM ] = 1;
           break;
 
+        case UI_KEY_MEAN:
+          name           = "MEAN";
+          unit           = p->input->unit ? p->input->unit : "pixelunit";
+          ocomment       = "Mean of sky subtracted values.";
+          ccomment       = "Mean of pixels subtracted by rivers.";
+          otype          = GAL_TYPE_FLOAT32;
+          ctype          = GAL_TYPE_FLOAT32;
+          disp_fmt       = GAL_TABLE_DISPLAY_FMT_GENERAL;
+          disp_width     = 10;
+          disp_precision = 4;
+          oiflag[ OCOL_NUM ] = 1;
+          oiflag[ OCOL_SUM ] = 1;
+          ciflag[ CCOL_NUM ] = 1;
+          ciflag[ CCOL_SUM ] = 1;
+          ciflag[ CCOL_RIV_NUM ] = 1;
+          ciflag[ CCOL_RIV_SUM ] = 1;
+          break;
+
+        case UI_KEY_MEDIAN:
+          name           = "MEDIAN";
+          unit           = p->input->unit ? p->input->unit : "pixelunit";
+          ocomment       = "Median of sky subtracted values.";
+          ccomment       = "Median of pixels subtracted by rivers.";
+          otype          = GAL_TYPE_FLOAT32;
+          ctype          = GAL_TYPE_FLOAT32;
+          disp_fmt       = GAL_TABLE_DISPLAY_FMT_GENERAL;
+          disp_width     = 10;
+          disp_precision = 4;
+          oiflag[ OCOL_MEDIAN  ] = 1;
+          oiflag[ OCOL_NUMALL  ] = 1;
+          ciflag[ CCOL_MEDIAN  ] = 1;
+          ciflag[ CCOL_NUMALL  ] = 1;
+          ciflag[ CCOL_RIV_NUM ] = 1;
+          ciflag[ CCOL_RIV_SUM ] = 1;
+          break;
+
         case UI_KEY_MAGNITUDE:
           name           = "MAGNITUDE";
           unit           = "log";
@@ -1585,6 +1621,18 @@ columns_fill(struct mkcatalog_passparams *pp)
                                       : NAN );
           break;
 
+        case UI_KEY_MEAN:
+          ((float *)colarr)[oind] = ( oi[ OCOL_NUM ]>0.0f
+                                      ? oi[ OCOL_SUM ] / oi[ OCOL_NUM ]
+                                      : NAN );
+          break;
+
+        case UI_KEY_MEDIAN:
+          ((float *)colarr)[oind] = ( oi[ OCOL_NUM ]>0.0f
+                                      ? oi[ OCOL_MEDIAN ]
+                                      : NAN );
+          break;
+
         case UI_KEY_MAGNITUDE:
           ((float *)colarr)[oind] = MKC_MAG(oi[ OCOL_SUM ]);
           break;
@@ -1760,6 +1808,24 @@ columns_fill(struct mkcatalog_passparams *pp)
                                         ? ci[ CCOL_SUM ] : NAN );
             break;
 
+          case UI_KEY_MEAN:
+            /* Similar to brightness. */
+            tmp = ( ci[ CCOL_RIV_NUM ]>0.0f
+                    ? ci[ CCOL_RIV_SUM ]/ci[ CCOL_RIV_NUM ]
+                    : 0 );
+
+            /* Subtract it from the clump's mean. */
+            ((float *)colarr)[cind] = ( ci[ CCOL_NUM ]>0.0f
+                                        ? (ci[CCOL_SUM]/ci[CCOL_NUM] - tmp)
+                                        : NAN );
+
+            break;
+
+          case UI_KEY_MEDIAN:
+            ((float *)colarr)[cind] = ( ci[ CCOL_NUM ]>0.0f
+                                        ? ci[ CCOL_MEDIAN ] : NAN );
+            break;
+
           case UI_KEY_MAGNITUDE: /* Similar: brightness for clumps */
             tmp = ( ci[ CCOL_RIV_NUM ]
                     ? ci[ CCOL_RIV_SUM ]/ci[ CCOL_RIV_NUM ]*ci[ CCOL_NUM ]
diff --git a/bin/mkcatalog/main.h b/bin/mkcatalog/main.h
index 6ae8b49..7b5f0c2 100644
--- a/bin/mkcatalog/main.h
+++ b/bin/mkcatalog/main.h
@@ -69,6 +69,7 @@ enum objectcols
     OCOL_NUMALL,         /* Number of all pixels with this label.     */
     OCOL_NUM,            /* Number of pixels above threshold.         */
     OCOL_SUM,            /* Sum of (value-sky) in object.             */
+    OCOL_MEDIAN,         /* Median of value in object.                */
     OCOL_VX,             /* Sum of (value-sky) * x.                   */
     OCOL_VY,             /* Sum of (value-sky) * y.                   */
     OCOL_VZ,             /* Sum of (value-sky) * z.                   */
@@ -91,7 +92,7 @@ enum objectcols
     OCOL_UPPERLIMIT_B,   /* Upper limit magnitude.                    */
     OCOL_C_NUMALL,       /* Value independent no. of pixels in clumps.*/
     OCOL_C_NUM,          /* Area of clumps in this object.            */
-    OCOL_C_SUM,          /* Brightness  in object clumps.             */
+    OCOL_C_SUM,          /* Brightness in object clumps.              */
     OCOL_C_VX,           /* Sum of (value-sky)*x on clumps.           */
     OCOL_C_VY,           /* Sum of (value-sky)*y on obj. clumps.      */
     OCOL_C_VZ,           /* Sum of (value-sky)*z on obj. clumps.      */
@@ -115,6 +116,7 @@ enum clumpcols
     CCOL_VYY,            /* Sum of flux*y*y of this clump.            */
     CCOL_VXY,            /* Sum of flux*x*y of this clump.            */
     CCOL_SUM,            /* River subtracted brightness.              */
+    CCOL_MEDIAN,         /* Median of values in clump.                */
     CCOL_RIV_SUM,        /* Sum of rivers around clump.               */
     CCOL_RIV_NUM,        /* Num river pixels around this clump.       */
     CCOL_SUMSKY,         /* Sum of sky value on this object.          */
@@ -163,7 +165,7 @@ struct mkcatalogparams
   char            *upmaskfile;  /* Name of upper limit mask file.       */
   char             *upmaskhdu;  /* HDU of upper limit mask file.        */
   size_t                upnum;  /* Number of upper-limit random samples.*/
-  size_t             *uprange;  /* Range of random positions about target. */
+  size_t             *uprange;  /* Range of random positions about target.*/
   uint8_t             envseed;  /* Use the environment for random seed. */
   double       upsigmaclip[2];  /* Sigma clip to measure upper limit.   */
   float              upnsigma;  /* Multiple of sigma to define up-lim.  */
diff --git a/bin/mkcatalog/mkcatalog.c b/bin/mkcatalog/mkcatalog.c
index 8e5e02c..2fafb30 100644
--- a/bin/mkcatalog/mkcatalog.c
+++ b/bin/mkcatalog/mkcatalog.c
@@ -37,6 +37,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/fits.h>
 #include <gnuastro/threads.h>
 #include <gnuastro/dimension.h>
+#include <gnuastro/statistics.h>
 
 #include <gnuastro-internal/timing.h>
 
@@ -417,6 +418,125 @@ mkcatalog_second_pass(struct mkcatalog_passparams *pp)
 
 
 
+static void
+mkcatalog_median_pass(struct mkcatalog_passparams *pp)
+{
+  struct mkcatalogparams *p=pp->p;
+  size_t ndim=p->input->ndim, *dsize=p->input->dsize;
+
+  double *ci;
+  gal_data_t *median;
+  int32_t *O, *C=NULL;
+  gal_data_t **clumpsmed=NULL;
+  float ss, *I, *II, *SK, *ST;
+  size_t i, increment=0, num_increment=1;
+  size_t counter=0, *ccounter=NULL, tsize=pp->oi[OCOL_NUM];
+  gal_data_t *objmed=gal_data_alloc(NULL, p->input->type, 1, &tsize, NULL, 0,
+                                    p->cp.minmapsize, NULL, NULL, NULL);
+
+  /* A small sanity check. */
+  if(p->input->type!=GAL_TYPE_FLOAT32)
+    error(EXIT_FAILURE, 0, "%s: the input has to be float32 type", __func__);
+
+
+  /* Allocate space for the clump medians. */
+  if(p->clumps)
+    {
+      errno=0;
+      clumpsmed=malloc(pp->clumpsinobj * sizeof *clumpsmed);
+      if(clumpsmed==NULL)
+        error(EXIT_FAILURE, errno, "%s: couldn't allocate `clumpsmed' for "
+              "%zu clumps", __func__, pp->clumpsinobj);
+
+
+      /* Allocate the array necessary to keep the values of each clump. */
+      ccounter=gal_data_calloc_array(GAL_TYPE_SIZE_T, pp->clumpsinobj,
+                                     __func__, "ccounter");
+      for(i=0;i<pp->clumpsinobj;++i)
+        {
+          tsize=pp->ci[ i * CCOL_NUMCOLS + CCOL_NUM ];
+          clumpsmed[i]=gal_data_alloc(NULL, p->input->type, 1, &tsize, NULL,
+                                      0, p->cp.minmapsize, NULL, NULL, NULL);
+        }
+    }
+
+
+  /* Parse each contiguous patch of memory covered by this object. */
+  while( pp->start_end_inc[0] + increment <= pp->start_end_inc[1] )
+    {
+      /* Set the contiguous range to parse, we will check the count
+         over the `I' pointer and just increment the rest. */
+      O  = pp->st_o   + increment;
+      SK = pp->st_sky + increment;
+      ST = pp->st_std + increment;
+      if(p->clumps) C = pp->st_c + increment;
+      II = ( I = pp->st_i + increment ) + pp->tile->dsize[ndim-1];
+
+      /* Parse the next contiguous region of this tile. */
+      do
+        {
+          /* If this pixel belongs to the requested object, then do the
+             processing. `hasblank' is constant, so when the input doesn't
+             have any blank values, the `isnan' will never be checked. */
+          if( *O==pp->object
+              && !( p->hasblank && isnan(*I) )
+              && !( (ss = *I - *SK) < p->threshold * *ST ))
+            {
+              /* Copy the value for the whole object. */
+              ss = *I - *SK;
+              memcpy( gal_data_ptr_increment(objmed->array, counter++,
+                                             p->input->type),
+                      &ss, gal_type_sizeof(p->input->type) );
+
+              /* We are also on a clump. */
+              if(p->clumps && *C>0)
+                memcpy( gal_data_ptr_increment(clumpsmed[*C-1]->array,
+                                               ccounter[*C-1]++,
+                                               p->input->type),
+                        &ss, gal_type_sizeof(p->input->type) );
+            }
+
+          /* Increment the other pointers. */
+          ++O; ++SK; ++ST; if(p->clumps) ++C;
+        }
+      while(++I<II);
+
+      /* Increment to the next contiguous region of this tile. */
+      increment += ( gal_tile_block_increment(p->input, dsize,
+                                              num_increment++, NULL) );
+    }
+
+
+  /* Calculate the final medians for objects. */
+  median=gal_data_copy_to_new_type_free(gal_statistics_median(objmed, 1),
+                                        GAL_TYPE_FLOAT64);
+  pp->oi[OCOL_MEDIAN]=*((double *)(median->array));
+  gal_data_free(objmed);
+  gal_data_free(median);
+
+
+  /* Calculate the median for clumps. */
+  if(p->clumps)
+    {
+      for(i=0;i<pp->clumpsinobj;++i)
+        {
+          ci=&pp->ci[ i * CCOL_NUMCOLS ];
+          median=gal_statistics_median(clumpsmed[i], 1);
+          median=gal_data_copy_to_new_type_free(median, GAL_TYPE_FLOAT64);
+          ci[ CCOL_MEDIAN ] = ( *((double *)(median->array))
+                                - (ci[ CCOL_RIV_SUM ]/ci[ CCOL_RIV_NUM ]) );
+          gal_data_free(clumpsmed[i]);
+          gal_data_free(median);
+        }
+      free(clumpsmed);
+      free(ccounter);
+    }
+}
+
+
+
+
+
 
 
 
@@ -504,7 +624,7 @@ mkcatalog_single_object(void *in_prm)
       mkcatalog_first_pass(&pp);
 
       /* Currently the second pass is only necessary when there is a clumps
-         image. */
+         image or the median is requested. */
       if(p->clumps)
         {
           /* Allocate space for the properties of each clump. */
@@ -521,6 +641,10 @@ mkcatalog_single_object(void *in_prm)
           mkcatalog_second_pass(&pp);
         }
 
+      /* If the median is requested, another pass is necessary. */
+      if( p->oiflag[ OCOL_MEDIAN ] )
+        mkcatalog_median_pass(&pp);
+
       /* Calculate the upper limit magnitude (if necessary). */
       if(p->upperlimit) upperlimit_calculate(&pp);
 
diff --git a/bin/mkcatalog/ui.h b/bin/mkcatalog/ui.h
index e5a32a5..7ecc919 100644
--- a/bin/mkcatalog/ui.h
+++ b/bin/mkcatalog/ui.h
@@ -123,6 +123,8 @@ enum option_keys_enum
   UI_KEY_CLUMPSGEOW3,
   UI_KEY_CLUMPSBRIGHTNESS,
   UI_KEY_NORIVERBRIGHTNESS,
+  UI_KEY_MEAN,
+  UI_KEY_MEDIAN,
   UI_KEY_CLUMPSMAGNITUDE,
   UI_KEY_UPPERLIMIT,
   UI_KEY_RIVERAVE,
diff --git a/bin/statistics/statistics.c b/bin/statistics/statistics.c
index 62bf0ec..1d10f8e 100644
--- a/bin/statistics/statistics.c
+++ b/bin/statistics/statistics.c
@@ -97,7 +97,6 @@ statistics_print_one_row(struct statisticsparams *p)
   gal_data_t *tmpv, *out=NULL, *num=NULL, *min=NULL, *max=NULL;
   gal_data_t *sum=NULL, *med=NULL, *meanstd=NULL, *modearr=NULL;
 
-
   /* The user can ask for any of the operators more than once, also some
      operators might return more than one usable value (like mode). So we
      will calculate the desired values once, and then print them. */
@@ -1016,5 +1015,4 @@ statistics(struct statisticsparams *p)
   /* If nothing was requested print the simple statistics. */
   if(print_basic_info)
     print_basics(p);
-
 }
diff --git a/bin/statistics/ui.c b/bin/statistics/ui.c
index 2e1d57f..bed2975 100644
--- a/bin/statistics/ui.c
+++ b/bin/statistics/ui.c
@@ -654,7 +654,7 @@ ui_out_of_range_to_blank(struct statisticsparams *p)
   /* Allocate a blank value to mask all pixels that don't satisfy the
      condition. */
   blank=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, 1, &one, NULL,
-                     0, -1, NULL, NULL, NULL);
+                       0, -1, NULL, NULL, NULL);
   *((float *)(blank->array)) = NAN;
 
 
@@ -662,6 +662,11 @@ ui_out_of_range_to_blank(struct statisticsparams *p)
      blank value will be used in the proper type of the input in the
      `where' operator.*/
   gal_arithmetic(GAL_ARITHMETIC_OP_WHERE, flagsor, p->input, cond, blank);
+
+
+  /* Reset the blank flags so they are checked again if necessary. */
+  p->input->flag &= ~GAL_DATA_FLAG_BLANK_CH;
+  p->input->flag &= ~GAL_DATA_FLAG_HASBLANK;
 }
 
 
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index c0b76b6..2abbdf5 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -830,21 +830,45 @@ was changed.
 demonstrates how four data sets with widely different shapes (when plotted)
 give nearly identical output from standard regression techniques. Anscombe
 uses this (now famous) quartet, which was introduced in the paper quoted
-above, to argue that ``Good statistical analysis is not a purely routine
-matter, and generally calls for more than one pass through the
-computer''. Anscombe's quartet can be generalized to say that users of a
-software cannot claim to understand how it works only based on the
-experience they have gained by frequently using it. This kind of subjective
-experience is prone to very serious mis-understandings about the data, what
-the software/statistical-method really does (especially as it gets more
+above, to argue that ``@emph{Good statistical analysis is not a purely
+routine matter, and generally calls for more than one pass through the
+computer}''. Echoing Anscombe's concern after 44 years, some of the highly
+recognized statisticians of our time (Leek, McShane, Gelman, Colquhoun,
+Nuijten and Goodman), wrote in Nature that:
+
+@quotation
+We need to appreciate that data analysis is not purely computational and
+algorithmic — it is a human behaviour....Researchers who hunt hard enough
+will turn up a result that fits statistical criteria — but their discovery
+will probably be a false positive.
+@author Five ways to fix statistics, Nature, 551, Nov 2017.
+@end quotation
+
+Users of statistical (scientific) methods (software) are therefore not
+passive (objective) agents in their result. Therefore, it is necessary to
+actually understand the method not just use it as a black box. The
+subjective experience gained by frequently using a method/software is not
+sufficient to claim an understanding of how the tool/method works and how
+relevant it is to the data and analysis. This kind of subjective experience
+is prone to very serious mis-understandings about the data, what the
+software/statistical-method really does (especially as it gets more
 complicated), and thus the scientific interpretation of the result. This
 attitude is further encouraged through non-free
-software@footnote{@url{https://www.gnu.org/philosophy/free-sw.html}}. This
-approach to scientific software only helps in producing dogmas and an
-``obscurantist faith in the expert's special skill, and in his personal
-knowledge and authority''@footnote{Karl Popper. The logic of scientific
-discovery. 1959. Larger quote is given at the start of the PDF (for print)
-version of this book.}.
+software@footnote{@url{https://www.gnu.org/philosophy/free-sw.html}},
+poorly written (or non-existant) scientific software manuals, and
+non-reproducible papers@footnote{Where the authors omit many of the
+analysis/processing ``details'' from the paper by arguing that they would
+make the paper too long/unreadable. However, software methods do allows us
+to supplement papers with all the details necessary to exactly reproduce
+the result. For example see @url{https://doi.org/10.5281/zenodo.1163746,
+zenodo.1163746} and @url{https://doi.org/10.5281/zenodo.1164774,
+zenodo.1164774} and this @url{
+http://akhlaghi.org/reproducible-science.html, general discussion}.}. This
+approach to scientific software and methods only helps in producing dogmas
+and an ``@emph{obscurantist faith in the expert's special skill, and in his
+personal knowledge and authority}''@footnote{Karl Popper. The logic of
+scientific discovery. 1959. Larger quote is given at the start of the PDF
+(for print) version of this book.}.
 
 @quotation
 @cindex Douglas Rushkoff
@@ -907,8 +931,8 @@ Imagine if Galileo did not have the technical knowledge to 
build a
 telescope. Astronomical objects could not be seen with the Dutch military
 design of the telescope. In the beginning of his ``The Sidereal Messenger''
 (1610) he cautions the readers on this issue and instructs them on how to
-build a suitable instrument: without a detailed description of ``how'' he
-made his observations, no one would believe him. The same is true today,
+build a suitable instrument: without a detailed description of @emph{how}
+he made his observations, no one would believe him. The same is true today,
 science cannot progress with a black box. Before he actually saw the moons
 of Jupiter, the mountains on the Moon or the crescent of Venus, he was
 “evasive” to Kepler@footnote{Galileo G. (Translated by Maurice
@@ -917,17 +941,19 @@ first edition, 2008.}. Science is not independent of its 
tools.
 
 @cindex Ken Thomson
 @cindex Stroustrup, Bjarne
-Bjarne Stroustrup (creator of the C++ language) says: ``Without
-understanding software, you are reduced to believing in magic''.  Ken
-Thomson (the designer or the Unix operating system) says ``I abhor a system
-designed for the `user' if that word is a coded pejorative meaning `stupid
-and unsophisticated'.'' Certainly no scientist (user of a scientific
-software) would want to be considered a believer in magic, or `stupid and
-unsophisticated'. However, this can happen when scientists get too distant
-from the raw data and are mainly indulging themselves in their own
-high-level (abstract) models (creations). For example, roughly five years
-before special relativity and about two decades before quantum mechanics
-fundamentally changed Physics, Kelvin is quoted as saying:
+Bjarne Stroustrup (creator of the C++ language) says: ``@emph{Without
+understanding software, you are reduced to believing in magic}''.  Ken
+Thomson (the designer or the Unix operating system) says ``@emph{I abhor a
+system designed for the `user' if that word is a coded pejorative meaning
+`stupid and unsophisticated'}.'' Certainly no scientist (user of a
+scientific software) would want to be considered a believer in magic, or
+stupid and unsophisticated. However, this often happen when scientists get
+too distant from the raw data and methods and are mainly indulging
+themselves in their own high-level (abstract) models (creations).
+
+Roughly five years before special relativity and about two decades before
+quantum mechanics fundamentally changed Physics, Kelvin is quoted as
+saying:
 
 @quotation
 @cindex Lord Kelvin
@@ -2391,22 +2417,24 @@ unambiguous. In short, you often don't need to type the 
full name of the
 program you want to run.
 
 Gnuastro contains a large number of programs and it is natural to forget
-the details of each program's options. Therefore, before starting the
-analysis, let's review how you can refresh your memory any time you want on
-the command-line while working (without taking your hands off the
-keyboard). This book comes with your installation, so it will always
-correspond to your installed version of Gnuastro. Please see @ref{Info} for
-more.
+the details of each program's options or inputs and outputs. Therefore,
+before starting the analysis, let's review how you can access this book to
+refresh your memory any time you want. For example when working on the
+command-line, without having to take your hands off the keyboard. When you
+install Gnuastro, this book is also installed on your system along with all
+the programs and libraries, so you don't need an internet connection to to
+access/read it. Also, by accessing this book as described below, you can be
+sure that it corresponds to your installed version of Gnuastro.
 
 @cindex GNU Info
 GNU Info@footnote{GNU Info is already available on almost all Unix-like
 operating systems.} is the program in charge of displaying the manual on
-the command-line. To see this whole book on your command-line, please run
-the following command and subsequent keys. Info has its own
-mini-environment, therefore we'll show the keys that must be pressed in the
-mini-environment after a @code{->} sign. You can also ignore anything after
-the @code{#} sign in the middle of the line, they are only for your
-information.
+the command-line (for more, see @ref{Info}). To see this whole book on your
+command-line, please run the following command and press subsequent
+keys. Info has its own mini-environment, therefore we'll show the keys that
+must be pressed in the mini-environment after a @code{->} sign. You can
+also ignore anything after the @code{#} sign in the middle of the line,
+they are only for your information.
 
 @example
 $ info gnuastro                # Open the top of the manual.
@@ -2722,22 +2750,22 @@ $ for z in $(seq 0.5 0.1 5); do                         
         \
   done
 @end example
 
-If you want a fast result and commonly need such processing for a larger
-number of redshifts, the command above can be slow. This is because the
-CosmicCalculator program has a lot of overhead: it has to parse the
-command-line and all configuration files (see below). These are both
-human-readable characters, not computer-readable bits. CosmicCalculator
-then has to check the sanity of its inputs and check which of its many
-options you have asked for. It has to do all of these for every
-redshift. To greatly speed up the processing, Gnuastro gives you direct
-access to the root work-horse of CosmicCalculator without all that
-overhead: @ref{Gnuastro library}.
+This is a fast and simple way for this repeated calculation when it is only
+necessary once. However, if you commonly need this calculation and possibly
+for a larger number of redshifts, the command above can be slow. This is
+because the CosmicCalculator program has a lot of overhead. To be generic
+and easy to operate, it has to parse the command-line and all configuration
+files (see below) which contain human-readable characters and need a lot of
+processing to be ready for processing by the computer. Afterwards,
+CosmicCalculator has to check the sanity of its inputs and check which of
+its many options you have asked for. It has to do all of these for every
+redshift in the loop above.
 
-Using Gnuastro's library, you can write your own tiny little program for
-this same calculation, without all that extra overhead of CosmicCalculator
-(or any of Gnuastro's programs, see @ref{Library}). For this particular
-job, you want Gnuastro's @ref{Cosmology library}. Here is one example: put
-the following small C program in a file called @file{myprogram.c}.
+To greatly speed up the processing, you can directly access the root
+work-horse of CosmicCalculator without all that overhead. Using Gnuastro's
+library, you can write your own tiny program particularly designed for this
+exact calculation (and nothing else!). To do that, copy and paste the
+following C program in a file called @file{myprogram.c}.
 
 @example
 #include <math.h>
@@ -2774,45 +2802,46 @@ main(void)
 @end example
 
 @noindent
-To build and run this C program, you can use @ref{BuildProgram}. It is
-designed to manage Gnuastro's dependencies, compile the program you give it
-and then run it. In short, it hides all the complexities of compiling,
-linking and running C programs based on Gnuastro's library. The library
-will later be usable higher level languages like Python or Julia, for now
-it is only usable by C and C++ programs.
+To greatly simplify the compilation, linking and running of simple C
+programs like this that use Gnuastro's library, Gnuastro has
+@ref{BuildProgram}. This program designed to manage Gnuastro's
+dependencies, compile and link the program and then run the new program. To
+build and run the program above, simply run the following command:
 
 @example
 $ astbuildprog myprogram.c
 @end example
 
-See how much faster this is compared to the shell loop we wrote above? You
-might have noticed that a new file called @file{myprogram} is also created
-in the directory. This is the compiled program that was created and run by
-the command above (its in binary machine code, not human-readable any
-more). You can run it again to get the same results with a command like
-this:
+Did you notice how much faster this was compared to the shell loop we wrote
+above? You might have noticed that a new file called @file{myprogram} is
+also created in the directory. This is the compiled program that was
+created and run by the command above (its in binary machine code format,
+not human-readable any more). You can run it again to get the same results
+with a command like this:
 
 @example
 $ ./myprogram
 @end example
 
 The efficiency of @file{myprogram} compared to CosmicCalculator is because
-the requested processing is faster/comparable to the overheads. For other
-programs that take large input datasets (images for example), the overhead
-is usually negligible compared to the processing. In such cases, the
-libraries are only useful if you want a different/new processing compared
-to the functionalities in Gnuastro's existing programs.
-
-Gnuastro has a large library which all the programs use for various steps
-their processing. For the full list of available functions classified by
-context, please see @ref{Gnuastro library}. Gnuastro's library and
-BuildProgram are created to make it easy for you to use these powerful
-features and get to your scientific results as efficiently (fast) and
-accurately as possible. Several other complete working examples (involving
-images and tables) of Gnuastro's libraries can be see in @ref{Library demo
-programs}. We'll stop with the libraries at this point in this tutorial and
-get back to Gnuastro's already built programs which ready to be used on the
-command-line.
+the requested processing is faster/comparable to the overheads necessary
+for each processing. For other programs that take large input datasets
+(images for example), the overhead is usually negligible compared to the
+processing. In such cases, the libraries are only useful if you want a
+different/new processing compared to the functionalities in Gnuastro's
+existing programs.
+
+Gnuastro has a large library which is heavily used by all the programs. In
+other words, the library is like the skeleton of Gnuastro. For the full
+list of available functions classified by context, please see @ref{Gnuastro
+library}. Gnuastro's library and BuildProgram are created to make it easy
+for you to use these powerful features as you like. This gives you a high
+level of creativity, while also providing efficiency and
+robustness. Several other complete working examples (involving images and
+tables) of Gnuastro's libraries can be see in @ref{Library demo
+programs}. Let's stop the discussion on libraries at this point in this
+tutorial and get back to Gnuastro's already built programs which were the
+main purpose of this tutorial.
 
 None of Gnuastro's programs keep a default value internally within their
 code. However, when you ran CosmicCalculator with the @option{-z2} option
@@ -15147,7 +15176,7 @@ variable to the @code{mkcatalogparams} structure.
 
 @item ui.h
 The @code{option_keys_enum} associates a unique value for each option to
-MakeProfiles. The options that have a short option version, the single
+MakeCatalog. The options that have a short option version, the single
 character short comment is used for the value. Those that don't have a
 short option version, get a large integer automatically. You should add a
 variable here to identify your desired column.
@@ -15167,15 +15196,14 @@ sanity checks and preparations for it here. 
Otherwise, you can ignore this
 file.
 
 @item columns.c
-This file will contain the main definition and high-level calculation of
-your new column through the @code{columns_define_alloc} and
-@code{columns_fill} functions. In the first, you specify the basic
-information about the column: its name, units, comments, type (see
-@ref{Numeric data types}) and how it should be printed if the output is a
-text file. You should also specify the raw/internal columns that are
-necessary for this column here as the many existing examples show. Through
-the types for objects and rows, you can specify if this column is only for
-clumps, objects or both.
+This file contains the main definition and high-level calculation of your
+new column through the @code{columns_define_alloc} and @code{columns_fill}
+functions. In the first, you specify the basic information about the
+column: its name, units, comments, type (see @ref{Numeric data types}) and
+how it should be printed if the output is a text file. You should also
+specify the raw/internal columns that are necessary for this column here as
+the many existing examples show. Through the types for objects and rows,
+you can specify if this column is only for clumps, objects or both.
 
 The second main function (@code{columns_fill}) writes the final value into
 the appropriate column for each object and clump. As you can see in the
@@ -15188,6 +15216,9 @@ work-horses: @code{mkcatalog_first_pass} and 
@code{mkcatalog_second_pass},
 their names are descriptive enough and their internals are also clear and
 heavily commented.
 
+@item doc/gnuastro.texi
+Update this manual and add a description for the new column.
+
 @end table
 
 
@@ -15729,6 +15760,16 @@ usable pixels (blank or below the threshold) are 
present over the clump or
 object, the stored value will be NaN, because zero (note that zero is
 meaningful).
 
+@item --mean
+The mean sky subtracted value of pixels within the object or clump. For
+clumps, the average river flux is subtracted from the sky subtracted
+mean.
+
+@item --median
+The median sky subtracted value of pixels within the object or clump. For
+clumps, the average river flux is subtracted from the sky subtracted
+median.
+
 @item --noriverbrightness
 [Clumps] The Sky (not river) subtracted clump brightness. By definition,
 for the clumps, the average brightness of the rivers surrounding it are
@@ -23011,15 +23052,14 @@ use to access an element of the @code{OTHER} dataset 
(if
 (similar to @code{itype} and @code{INPUT} discussed above). If
 @code{PARSE_OTHER} is non-zero, then @code{o} will also be incremented to
 the same index element but in the other array. You can use these along with
-any other variable you define before this macro to process the input and
-store it in the other.
+any other variable you define before this macro to process the input and/or
+the other.
 
 All variables within this function-like macro begin with @code{tpo_} except
-for the three variables listed above. So as long as you don't start the
-names of your variables with this prefix everything will be fine. Note that
-@code{i} (and possibly @code{o}) will be incremented once by this
-function-like macro, so don't increment them within @code{OP}. As a
-summary, the three variables you have access to within @code{OP} are:
+for the three variables listed below. Therefore, as long as you don't start
+the names of your variables with this prefix everything will be fine. Note
+that @code{i} (and possibly @code{o}) will be incremented once by this
+function-like macro, so don't increment them within @code{OP}.
 
 @table @code
 @item i
diff --git a/lib/statistics.c b/lib/statistics.c
index 4ff9b79..ce71667 100644
--- a/lib/statistics.c
+++ b/lib/statistics.c
@@ -80,9 +80,7 @@ gal_statistics_number(gal_data_t *input)
 
 
 /* Return the minimum (non-blank) value of a dataset in the same type as
-   the dataset. Note that a NaN (blank in floating point) will fail on any
-   comparison. So when finding the minimum or maximum, when the blank value
-   is NaN, we can safely assume there is no blank value at all. */
+   the dataset. */
 gal_data_t *
 gal_statistics_minimum(gal_data_t *input)
 {
@@ -94,7 +92,7 @@ gal_statistics_minimum(gal_data_t *input)
   gal_type_max(out->type, out->array);
 
   /* Parse the full input. */
-  GAL_TILE_PARSE_OPERATE(input, out, 0, 0, {*o = *i < *o ? *i : *o; ++n;});
+  GAL_TILE_PARSE_OPERATE(input, out, 0, 1, {*o = *i < *o ? *i : *o; ++n;});
 
   /* If there were no usable elements, set the output to blank, then
      return. */
@@ -107,7 +105,7 @@ gal_statistics_minimum(gal_data_t *input)
 
 
 /* Return the maximum (non-blank) value of a dataset in the same type as
-   the dataset. See explanations of `gal_statistics_minimum'. */
+   the dataset. */
 gal_data_t *
 gal_statistics_maximum(gal_data_t *input)
 {
@@ -119,7 +117,7 @@ gal_statistics_maximum(gal_data_t *input)
   gal_type_min(out->type, out->array);
 
   /* Parse the full input. */
-  GAL_TILE_PARSE_OPERATE(input, out, 0, 0, {*o = *i > *o ? *i : *o; ++n;});
+  GAL_TILE_PARSE_OPERATE(input, out, 0, 1, {*o = *i > *o ? *i : *o; ++n;});
 
   /* If there were no usable elements, set the output to blank, then
      return. */