[gnuastro-commits] master d002ec8: MakeCatalog: correct measure Sky standard deviation with --std

[Mohammad Akhlaghi]

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

    MakeCatalog: correct measure Sky standard deviation with --std
    
    Until now, `--std' would simply add the standard deviations and divide them
    by their number (mean standard deviation). But this is not a statistically
    meaningful measure of the Sky standard deviation! The Sky standard
    deviation is the square root of the mean variance.
---
 NEWS                    |  7 +++++++
 bin/mkcatalog/args.h    |  4 ++--
 bin/mkcatalog/columns.c | 11 ++++++-----
 bin/mkcatalog/main.h    |  4 ++--
 bin/mkcatalog/parse.c   | 16 ++++++++--------
 bin/mkcatalog/ui.c      |  4 ++--
 doc/gnuastro.texi       |  6 +++---
 7 files changed, 30 insertions(+), 22 deletions(-)

diff --git a/NEWS b/NEWS
index b110e67..3290258 100644
--- a/NEWS
+++ b/NEWS
@@ -63,6 +63,13 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
      color maps in Gnuastro 0.8, it is also necessary to force a range on
      non-byte datasets. It is thus necessary to use a more generic name.
 
+  MakeCatalog:
+   --std: Until now, this option would measure the mean standard deviation
+     under the label. But this is not a statistically meaningful measure
+     for the Sky standard deviation and could be incorrectly used. From now
+     on, this option measures the square root of the mean variance, or root
+     mean square (the correct definition of the Sky standard deviation).
+
 ** Bugs fixed
   bug #55313: Fits program writing --write values in reverse order
   bug #55333: Fits program crash when --write or --update have no value.
diff --git a/bin/mkcatalog/args.h b/bin/mkcatalog/args.h
index 14b138f..07a24f9 100644
--- a/bin/mkcatalog/args.h
+++ b/bin/mkcatalog/args.h
@@ -1013,7 +1013,7 @@ struct argp_option program_options[] =
       UI_KEY_SKY,
       0,
       0,
-      "Average Sky value under this clump or object.",
+      "Sky value (per pixel).",
       UI_GROUP_COLUMNS_BRIGHTNESS,
       0,
       GAL_TYPE_INVALID,
@@ -1027,7 +1027,7 @@ struct argp_option program_options[] =
       UI_KEY_STD,
       0,
       0,
-      "Average Sky standard deviation under label.",
+      "Sky standard deviation (per pixel).",
       UI_GROUP_COLUMNS_BRIGHTNESS,
       0,
       GAL_TYPE_INVALID,
diff --git a/bin/mkcatalog/columns.c b/bin/mkcatalog/columns.c
index c3edacf..99863a2 100644
--- a/bin/mkcatalog/columns.c
+++ b/bin/mkcatalog/columns.c
@@ -1003,7 +1003,7 @@ columns_define_alloc(struct mkcatalogparams *p)
         case UI_KEY_STD:
           name           = "STD";
           unit           = MKCATALOG_NO_UNIT;
-          ocomment       = "Average of input standard deviation.";
+          ocomment       = "Sky standard deviation under object.";
           ccomment       = ocomment;
           otype          = GAL_TYPE_FLOAT32;
           ctype          = GAL_TYPE_FLOAT32;
@@ -1011,7 +1011,7 @@ columns_define_alloc(struct mkcatalogparams *p)
           disp_width     = 10;
           disp_precision = 4;
           oiflag[ OCOL_NUM    ] = ciflag[ CCOL_NUM    ] = 1;
-          oiflag[ OCOL_SUMSTD ] = ciflag[ CCOL_SUMSTD ] = 1;
+          oiflag[ OCOL_SUMVAR ] = ciflag[ CCOL_SUMVAR ] = 1;
           break;
 
         case UI_KEY_SEMIMAJOR:
@@ -1738,7 +1738,8 @@ columns_fill(struct mkcatalog_passparams *pp)
           break;
 
         case UI_KEY_STD:
-          ((float *)colarr)[oind] = MKC_RATIO(oi[OCOL_SUMSTD], oi[OCOL_NUM]);
+          ((float *)colarr)[oind] = sqrt( MKC_RATIO( oi[ OCOL_SUMVAR ],
+                                                     oi[ OCOL_NUM    ]) );
           break;
 
         case UI_KEY_SEMIMAJOR:
@@ -1954,8 +1955,8 @@ columns_fill(struct mkcatalog_passparams *pp)
             break;
 
           case UI_KEY_STD:
-            ((float *)colarr)[cind] = MKC_RATIO( ci[ CCOL_SUMSTD ],
-                                                 ci[ CCOL_NUM ] );
+            ((float *)colarr)[cind] = sqrt( MKC_RATIO( ci[ CCOL_SUMVAR ],
+                                                       ci[ CCOL_NUM    ] ));
             break;
 
           case UI_KEY_SEMIMAJOR:
diff --git a/bin/mkcatalog/main.h b/bin/mkcatalog/main.h
index 47faa0c..78e3451 100644
--- a/bin/mkcatalog/main.h
+++ b/bin/mkcatalog/main.h
@@ -81,7 +81,7 @@ enum objectcols
     OCOL_VYY,            /* Sum of (value-sky) * y * y.               */
     OCOL_VXY,            /* Sum of (value-sky) * x * y.               */
     OCOL_SUMSKY,         /* Sum of sky value on this object.          */
-    OCOL_SUMSTD,         /* Sum of sky STD value on this object.      */
+    OCOL_SUMVAR,         /* Sum of sky variance value on this object. */
     OCOL_SUMWHT,         /* Sum of positive image pixels.             */
     OCOL_NUMWHT,         /* Number of positive pixels used for wht.   */
     OCOL_GX,             /* Geometric center of object in X.          */
@@ -122,7 +122,7 @@ enum clumpcols
     CCOL_VYY,            /* Sum of flux*y*y of this clump.            */
     CCOL_VXY,            /* Sum of flux*x*y of this clump.            */
     CCOL_SUMSKY,         /* Sum of sky value on this object.          */
-    CCOL_SUMSTD,         /* Sum of sky STD value on this object.      */
+    CCOL_SUMVAR,         /* Sum of sky variance value on this object. */
     CCOL_SUMWHT,         /* Sum of positive image pixels for wht.     */
     CCOL_NUMWHT,         /* Num of positive image pixels for wht.     */
     CCOL_GX,             /* Geometric center of clump in X.           */
diff --git a/bin/mkcatalog/parse.c b/bin/mkcatalog/parse.c
index 052986e..12084c9 100644
--- a/bin/mkcatalog/parse.c
+++ b/bin/mkcatalog/parse.c
@@ -116,7 +116,7 @@ parse_objects(struct mkcatalog_passparams *pp)
   double *oi=pp->oi;
   size_t *tsize=pp->tile->dsize;
   size_t d, increment=0, num_increment=1;
-  float st, sval, *V=NULL, *SK=NULL, *ST=NULL;
+  float var, sval, *V=NULL, *SK=NULL, *ST=NULL;
   int32_t *O, *OO, *C=NULL, *objarr=p->objects->array;
   float *std=p->std?p->std->array:NULL, *sky=p->sky?p->sky->array:NULL;
 
@@ -271,8 +271,8 @@ parse_objects(struct mkcatalog_passparams *pp)
               if(p->std)
                 {
                   sval = pp->st_std ? *ST : (p->std->size>1?std[tid]:std[0]);
-                  st = p->variance ? sqrt(sval) : sval;
-                  if(oif[ OCOL_SUMSTD ]) oi[ OCOL_SUMSTD  ] += st;
+                  var = p->variance ? sval : sval*sval;
+                  if(oif[ OCOL_SUMVAR ]) oi[ OCOL_SUMVAR  ] += var;
                   /* For each pixel, we have a sky contribution to the
                      counts and the signal's contribution. The standard
                      deviation in the sky is simply `sval', but the
@@ -284,7 +284,7 @@ parse_objects(struct mkcatalog_passparams *pp)
                      noisy there will be negative values, and we don't want
                      them to decrease the variance. */
                   if(oif[ OCOL_SUM_VAR ])
-                    oi[ OCOL_SUM_VAR ] += ( (p->variance ? sval : st*st)
+                    oi[ OCOL_SUM_VAR ] += ( (p->variance ? var : sval)
                                             + fabs(*V) );
                 }
             }
@@ -335,7 +335,7 @@ parse_clumps(struct mkcatalog_passparams *pp)
   uint8_t *cif=p->ciflag;
   size_t *tsize=pp->tile->dsize;
   int32_t *O, *OO, *C=NULL, nlab;
-  float st, sval, *V=NULL, *SK=NULL, *ST=NULL;
+  float var, sval, *V=NULL, *SK=NULL, *ST=NULL;
   size_t i, ii, d, increment=0, num_increment=1;
   size_t nngb=gal_dimension_num_neighbors(ndim);
   int32_t *objects=p->objects->array, *clumps=p->clumps->array;
@@ -479,10 +479,10 @@ parse_clumps(struct mkcatalog_passparams *pp)
                       sval = ( pp->st_std
                                ? *ST
                                : (p->std->size>1 ? std[tid] : std[0]) );
-                      st = p->variance ? sqrt(sval) : sval;
-                      if(cif[ CCOL_SUMSTD  ]) ci[ CCOL_SUMSTD  ] += st;
+                      var = p->variance ? sval : sval*sval;
+                      if(cif[ CCOL_SUMVAR  ]) ci[ CCOL_SUMVAR ] += var;
                       if(cif[ CCOL_SUM_VAR ])
-                        ci[ CCOL_SUM_VAR ] += ( (p->variance ? sval : st*st)
+                        ci[ CCOL_SUM_VAR ] += ( (p->variance ? var : sval)
                                                 + fabs(*V) );
                     }
                 }
diff --git a/bin/mkcatalog/ui.c b/bin/mkcatalog/ui.c
index e39a529..5a9ee5c 100644
--- a/bin/mkcatalog/ui.c
+++ b/bin/mkcatalog/ui.c
@@ -709,7 +709,7 @@ ui_necessary_inputs(struct mkcatalogparams *p, int *values, 
int *sky,
         case OCOL_VYY:                *values        = 1;          break;
         case OCOL_VXY:                *values        = 1;          break;
         case OCOL_SUMSKY:             *sky           = 1;          break;
-        case OCOL_SUMSTD:             *std           = 1;          break;
+        case OCOL_SUMVAR:             *std           = 1;          break;
         case OCOL_SUMWHT:             *values        = 1;          break;
         case OCOL_NUMWHT:             *values        = 1;          break;
         case OCOL_GX:                 /* Only object labels. */    break;
@@ -757,7 +757,7 @@ ui_necessary_inputs(struct mkcatalogparams *p, int *values, 
int *sky,
           case CCOL_VYY:              *values        = 1;          break;
           case CCOL_VXY:              *values        = 1;          break;
           case CCOL_SUMSKY:           *sky           = 1;          break;
-          case CCOL_SUMSTD:           *std           = 1;          break;
+          case CCOL_SUMVAR:           *std           = 1;          break;
           case CCOL_SUMWHT:           *values        = 1;          break;
           case CCOL_NUMWHT:           *values        = 1;          break;
           case CCOL_GX:               /* Only clump labels. */     break;
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 3c8e35c..7160a02 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -19397,9 +19397,9 @@ actually the mean value of all the pixels in the sky 
image that lie on
 the same position as the object or clump.
 
 @item --std
-The sky value standard deviation (per pixel) for this clump or
-object. Like @option{--sky}, this is the average of the values in the
-input sky standard deviation image pixels that lie over this object.
+The sky value standard deviation (per pixel) for this clump or object. This
+is the square root of the mean variance under the object, or the root mean
+square.
 
 @item -C
 @itemx --numclumps