gnuastro-commits
[Top][All Lists]
Advanced

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

[gnuastro-commits] master aa55cd2 2/4: MakeCatalog's error is per-pixel,


From: Mohammad Akhlaghi
Subject: [gnuastro-commits] master aa55cd2 2/4: MakeCatalog's error is per-pixel, not average of all area
Date: Wed, 21 Mar 2018 15:11:49 -0400 (EDT)

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

    MakeCatalog's error is per-pixel, not average of all area
    
    Until now, MakeCatalog's error was calculated like this: the average error
    value over the object was found, then using the area and total flux, it was
    converted to the total object's error that would be used in the
    signal-to-noise estimation. This is good when the error (exposure time)
    under the object is flat. However, on surveys, we don't usually have this
    situation and the error can vary strongly on some parts (over one object).
    
    So we now estimate the error on a per-pixel basis and add the errors up (as
    variance). This made the signal-to-noise estimate much more
    cleaner/readable and we MakeCatalog also now accepts a `--brightnesserr'
    option to print the error in measuring the brightness.
    
    While checking, I also noticed that a `break' was missing in the final
    writing of the column values in MakeCatalog's `columns.c'. This is now
    corrected.
---
 NEWS                      |  11 +++--
 bin/mkcatalog/args.h      |  16 ++++++-
 bin/mkcatalog/columns.c   | 109 +++++++++++++++++++++++++++++-----------------
 bin/mkcatalog/main.h      |  11 +++--
 bin/mkcatalog/mkcatalog.c |  86 +++++++++++++++++++++++-------------
 bin/mkcatalog/ui.h        |   1 +
 doc/gnuastro.texi         |  11 ++++-
 7 files changed, 165 insertions(+), 80 deletions(-)

diff --git a/NEWS b/NEWS
index 3fab406..e72b424 100644
--- a/NEWS
+++ b/NEWS
@@ -10,6 +10,7 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
 
   MakeCatalog:
     --variance: input STD image is variance (so take square root before use).
+    --brightnesserr: error in estimating the brightness.
     --mean: calculate the mean pixel value within an object or clump.
     --median: calculate the median pixel value within an object or clump.
     --upperlimitsigma: position in random distribution (in units of sigma).
@@ -20,7 +21,7 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
     --manualbinrange: histogram or CFP range can be outside of distribution.
 
   Libraries:
-    gal_array_read: read from file with all known formats (FITS, TIFF, etc).
+    gal_array_read: read array from any of known formats (FITS,TIFF,JPEG,...).
     gal_array_read_to_type: similar to `gal_array_read', but to given type.
     gal_eps_name_is_eps: Returns 1 if given filename is EPS.
     gal_eps_suffix_is_eps: Returns 1 if given suffix is EPS.
@@ -49,10 +50,14 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
   NoiseChisel:
     --kernel: value `none' will disable convolution.
 
+  MakeCatalog:
+    - Estimation of noise-level is now done per-pixel over the whole
+         label. Until now the average noise level was used.
+
   Table:
     --column: multiple columns (comma separated) can be used in one
-              instance of this option (multiple instances of this option
-              are still acceptable also).
+         instance of this option (multiple instances of this option are
+         still acceptable also).
 
   Libraries:
     gal_fits_img_read: now only reads the data not the WCS, therefore it no
diff --git a/bin/mkcatalog/args.h b/bin/mkcatalog/args.h
index 0389d72..562b08c 100644
--- a/bin/mkcatalog/args.h
+++ b/bin/mkcatalog/args.h
@@ -688,6 +688,20 @@ struct argp_option program_options[] =
       ui_column_codes_ll
     },
     {
+      "brightnesserr",
+      UI_KEY_BRIGHTNESSERR,
+      0,
+      0,
+      "Error (1-sigma) in measuring brightness.",
+      UI_GROUP_COLUMNS_BRIGHTNESS,
+      0,
+      GAL_TYPE_INVALID,
+      GAL_OPTIONS_RANGE_ANY,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET,
+      ui_column_codes_ll
+    },
+    {
       "clumpbrightness",
       UI_KEY_CLUMPSBRIGHTNESS,
       0,
@@ -916,7 +930,7 @@ struct argp_option program_options[] =
       UI_KEY_STD,
       0,
       0,
-      "Average Sky standard deviation.",
+      "Average Sky standard deviation under label.",
       UI_GROUP_COLUMNS_BRIGHTNESS,
       0,
       GAL_TYPE_INVALID,
diff --git a/bin/mkcatalog/columns.c b/bin/mkcatalog/columns.c
index aa354e5..b6ef963 100644
--- a/bin/mkcatalog/columns.c
+++ b/bin/mkcatalog/columns.c
@@ -634,6 +634,21 @@ columns_define_alloc(struct mkcatalogparams *p)
           ciflag[ CCOL_RIV_SUM ] = 1;
           break;
 
+        case UI_KEY_BRIGHTNESSERR:
+          name           = "BRIGHTNESS_ERROR";
+          unit           = p->input->unit ? p->input->unit : "pixelunit";
+          ocomment       = "Error (1-sigma) in measuring brightness.";
+          ccomment       = ocomment;
+          otype          = GAL_TYPE_FLOAT32;
+          ctype          = GAL_TYPE_FLOAT32;
+          disp_fmt       = GAL_TABLE_DISPLAY_FMT_GENERAL;
+          disp_width     = 10;
+          disp_precision = 4;
+          oiflag[ OCOL_SUM_VAR     ] = 1;
+          ciflag[ CCOL_SUM_VAR     ] = 1;
+          ciflag[ CCOL_RIV_SUM_VAR ] = 1;
+          break;
+
         case UI_KEY_CLUMPSBRIGHTNESS:
           name           = "CLUMPS_BRIGHTNESS";
           unit           = p->input->unit ? p->input->unit : "pixelunit";
@@ -712,7 +727,7 @@ columns_define_alloc(struct mkcatalogparams *p)
           break;
 
         case UI_KEY_MAGNITUDEERR:
-          name           = "MAGNITUDE_ERR";
+          name           = "MAGNITUDE_ERROR";
           unit           = "log";
           ocomment       = "Error in measuring magnitude.";
           ccomment       = ocomment;
@@ -721,12 +736,12 @@ columns_define_alloc(struct mkcatalogparams *p)
           disp_fmt       = GAL_TABLE_DISPLAY_FMT_FLOAT;
           disp_width     = 8;
           disp_precision = 3;
-          oiflag[ OCOL_SUMSTD ] = 1;
-          oiflag[ OCOL_NUM    ] = 1;
-          oiflag[ OCOL_SUM    ] = 1;
-          ciflag[ CCOL_SUMSTD ] = 1;
-          ciflag[ CCOL_NUM    ] = 1;
-          ciflag[ CCOL_SUM    ] = 1;
+          oiflag[ OCOL_SUM         ] = 1;
+          oiflag[ OCOL_SUM_VAR     ] = 1;
+          ciflag[ CCOL_SUM         ] = 1;
+          ciflag[ CCOL_SUM_VAR     ] = 1;
+          ciflag[ CCOL_RIV_SUM     ] = 1;
+          ciflag[ CCOL_RIV_SUM_VAR ] = 1;
           break;
 
         case UI_KEY_CLUMPSMAGNITUDE:
@@ -846,12 +861,12 @@ columns_define_alloc(struct mkcatalogparams *p)
           disp_fmt       = GAL_TABLE_DISPLAY_FMT_FLOAT;
           disp_width     = 10;
           disp_precision = 3;
-          oiflag[ OCOL_SUMSTD ] = 1;
-          oiflag[ OCOL_NUM    ] = 1;
-          oiflag[ OCOL_SUM    ] = 1;
-          ciflag[ CCOL_SUMSTD ] = 1;
-          ciflag[ CCOL_NUM    ] = 1;
-          ciflag[ CCOL_SUM    ] = 1;
+          oiflag[ OCOL_SUM         ] = 1;
+          oiflag[ OCOL_SUM_VAR     ] = 1;
+          ciflag[ CCOL_SUM         ] = 1;
+          ciflag[ CCOL_SUM_VAR     ] = 1;
+          ciflag[ CCOL_RIV_SUM     ] = 1;
+          ciflag[ CCOL_RIV_SUM_VAR ] = 1;
           break;
 
         case UI_KEY_SKY:
@@ -1132,38 +1147,33 @@ columns_define_alloc(struct mkcatalogparams *p)
 
 
 
-/* Calculate the Signal to noise ratio for the object. */
+/* Calculate the error in brightness. */
 static double
-columns_sn(struct mkcatalogparams *p, double *row, int o0c1)
+columns_brightness_error(struct mkcatalogparams *p, double *row, int o0c1)
 {
-  double var, sn, std, Ni, I, O;
+  double V = row[ o0c1 ? CCOL_SUM_VAR : OCOL_SUM_VAR ];
+  double OV = (o0c1 && row[ CCOL_RIV_NUM ]) ? row[ CCOL_RIV_SUM_VAR ] : 0.0;
+  return sqrt(V+OV);
+}
+
 
-  /* Get all the values as averages (per pixel). */
-  Ni  = row[ o0c1 ? CCOL_NUM : OCOL_NUM ];
-  I   = MKC_RATIO( row[ o0c1 ? CCOL_SUM    : OCOL_SUM ],    Ni );
-  std = MKC_RATIO( row[ o0c1 ? CCOL_SUMSTD : OCOL_SUMSTD ], Ni );
-  var = (p->skysubtracted ? 2.0f : 1.0f) * std * std;
 
 
-  /* Calculate the S/N. Note that when grown clumps are requested from
-     NoiseChisel, some "clumps" will completely cover their objects and
-     there will be no rivers. So if this is a clump, and the river area is
-     0, we should treat the S/N as a an object. */
-  if( o0c1 && row[ CCOL_RIV_NUM ] )
-    {
-      /* If the Sky is already subtracted, the varience should be counted
-         two times. */
-      O   = row[ CCOL_RIV_SUM ] / row[ CCOL_RIV_NUM ];  /* Outside.  */
-      sn  = ( (I-O)>0
-              ? ( sqrt(Ni/p->cpscorr) * (I-O)
-                  / sqrt( (I>0?I:-1*I) + (O>0?O:-1*O) + var ) )
-              : NAN );
-    }
-  else
-    sn  = I>0 ? sqrt(Ni/p->cpscorr) * I / sqrt( (I>0?I:-1*I) + var ) : NAN;
+
+/* Calculate the Signal to noise ratio for the object. */
+static double
+columns_sn(struct mkcatalogparams *p, double *row, int o0c1)
+{
+  double I = row[ o0c1 ? CCOL_SUM     : OCOL_SUM     ];
+
+  /* When grown clumps are requested from NoiseChisel, some "clumps" will
+     completely cover their objects and there will be no rivers. So if this
+     is a clump, and the river area is 0, we should treat the S/N as a an
+     object. */
+  double O = (o0c1 && row[ CCOL_RIV_NUM ]) ? row[ CCOL_RIV_SUM ] : 0.0 ;
 
   /* Return the derived value. */
-  return sn;
+  return sqrt(1/p->cpscorr) * (I-O) / columns_brightness_error(p, row, o0c1);
 }
 
 
@@ -1312,7 +1322,11 @@ columns_clump_brightness(double *ci)
   `1/S'. So
 
       `DM = 2.5 / ( S * ln(10) )'               */
-#define MAG_ERROR(P,ROW,O0C1) (2.5f/(columns_sn((P),(ROW),(O0C1)) * log(10)))
+#define MAG_ERROR(P,ROW,O0C1) ( 2.5f                                    \
+                                / ( ( columns_sn((P),(ROW),(O0C1)) > 0  \
+                                      ? columns_sn((P),(ROW),(O0C1))    \
+                                      : NAN )                           \
+                                    * log(10) ) )
 
 
 
@@ -1447,6 +1461,12 @@ columns_fill(struct mkcatalog_passparams *pp)
                                       : NAN );
           break;
 
+        case UI_KEY_BRIGHTNESSERR:
+          ((float *)colarr)[oind] = ( oi[ OCOL_NUM ]>0.0f
+                                      ? columns_brightness_error(p, oi, 0)
+                                      : NAN );
+          break;
+
         case UI_KEY_CLUMPSBRIGHTNESS:
           ((float *)colarr)[oind] = ( oi[ OCOL_C_NUM ]>0.0f
                                       ? oi[ OCOL_C_SUM ]
@@ -1622,6 +1642,12 @@ columns_fill(struct mkcatalog_passparams *pp)
             ((float *)colarr)[cind] = columns_clump_brightness(ci);
             break;
 
+          case UI_KEY_BRIGHTNESSERR:
+            ((float *)colarr)[cind] = ( ci[ CCOL_NUM ]>0.0f
+                                        ? columns_brightness_error(p, ci, 1)
+                                        : NAN );
+            break;
+
           case UI_KEY_NORIVERBRIGHTNESS:
             ((float *)colarr)[cind] = ( ci[ CCOL_NUM ]>0.0f
                                         ? ci[ CCOL_SUM ] : NAN );
@@ -1663,6 +1689,7 @@ columns_fill(struct mkcatalog_passparams *pp)
           case UI_KEY_UPPERLIMITSIGMA:
             ((float *)colarr)[cind] = ( columns_clump_brightness(ci)
                                         / ci[ CCOL_UPPERLIMIT_S ] );
+            break;
 
           case UI_KEY_UPPERLIMITQUANTILE:
             ((float *)colarr)[cind] = ci[ CCOL_UPPERLIMIT_Q ];
@@ -1704,7 +1731,7 @@ columns_fill(struct mkcatalog_passparams *pp)
             ((float *)colarr)[cind]
               = ( columns_second_order(pp, ci, UI_KEY_SEMIMINOR, 1)
                   / columns_second_order(pp, ci, UI_KEY_SEMIMAJOR, 1) );
-          break;
+            break;
 
           case UI_KEY_POSITIONANGLE:
             ((float *)colarr)[cind] = columns_second_order(pp, ci, key, 1);
@@ -1722,7 +1749,7 @@ columns_fill(struct mkcatalog_passparams *pp)
             ((float *)colarr)[cind]
               = ( columns_second_order(pp, ci, UI_KEY_GEOSEMIMINOR, 1)
                   / columns_second_order(pp, ci, UI_KEY_GEOSEMIMAJOR, 1) );
-          break;
+            break;
 
           case UI_KEY_GEOPOSITIONANGLE:
             ((float *)colarr)[cind] = columns_second_order(pp, ci, key, 1);
diff --git a/bin/mkcatalog/main.h b/bin/mkcatalog/main.h
index 5e6e6ad..32857ad 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_SUM_VAR,        /* Varience of sum (for brightness error).   */
     OCOL_MEDIAN,         /* Median of value in object.                */
     OCOL_VX,             /* Sum of (value-sky) * x.                   */
     OCOL_VY,             /* Sum of (value-sky) * y.                   */
@@ -106,15 +107,17 @@ enum clumpcols
   {
     CCOL_NUMALL,         /* Area of clump irrespective of threshold.  */
     CCOL_NUM,            /* Area of this clump.                       */
+    CCOL_SUM,            /* River subtracted brightness.              */
+    CCOL_SUM_VAR,        /* Variance of sum (for brightness error).   */
+    CCOL_MEDIAN,         /* Median of values in clump.                */
+    CCOL_RIV_NUM,        /* Num river pixels around this clump.       */
+    CCOL_RIV_SUM,        /* Sum of rivers around clump.               */
+    CCOL_RIV_SUM_VAR,    /* Variance of sum (for error measurements). */
     CCOL_VX,             /* Sum of (value-sky) * x.                   */
     CCOL_VY,             /* Sum of (value-sky) * y.                   */
     CCOL_VXX,            /* Sum of flux*x*x of this clump.            */
     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.          */
     CCOL_SUMSTD,         /* Sum of sky STD value on this object.      */
     CCOL_SUMWHT,         /* Sum of positive image pixels for wht.     */
diff --git a/bin/mkcatalog/mkcatalog.c b/bin/mkcatalog/mkcatalog.c
index 9e224f6..12b7f13 100644
--- a/bin/mkcatalog/mkcatalog.c
+++ b/bin/mkcatalog/mkcatalog.c
@@ -124,7 +124,7 @@ mkcatalog_first_pass(struct mkcatalog_passparams *pp)
   double *oi=pp->oi;
   int32_t *O, *C=NULL;
   size_t d, increment=0, num_increment=1;
-  float ss, *I, *II, *SK, *ST, *input=p->input->array;
+  float ss, st, *I, *II, *SK, *ST, *input=p->input->array;
   size_t *c=gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim, __func__, "c");
   size_t *sc=gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim, __func__, "sc");
 
@@ -196,16 +196,31 @@ mkcatalog_first_pass(struct mkcatalog_passparams *pp)
                  negation will be positive and the calculations below will
                  be done. However, if the user does specify a threhold and
                  the pixel is above the threshold, then (`ss < p->threshold
-                 * *ST') will be false and its logical negation will be
+                 * st') will be false and its logical negation will be
                  positive, so the pixel will be included. */
+              st = p->variance ? sqrt(*ST) : *ST;
               if( !( p->hasblank && isnan(*I) )
-                  && !( (ss = *I - *SK) < p->threshold * *ST ) )
+                  && !( (ss = *I - *SK) < p->threshold * st ) )
                 {
                   /* General flux summations. */
-                  oi[ OCOL_NUM    ]++;
-                  oi[ OCOL_SUM    ] += ss;
-                  oi[ OCOL_SUMSKY ] += *SK;
-                  oi[ OCOL_SUMSTD ] += *ST;
+                  oi[ OCOL_NUM     ]++;
+                  oi[ OCOL_SUM     ] += ss;
+                  oi[ OCOL_SUMSKY  ] += *SK;
+                  oi[ OCOL_SUMSTD  ] += st;
+
+                  /* For each pixel, we have a sky contribution to the
+                     counts and the signal's contribution. The standard
+                     deviation in the sky is simply `*ST', but the standard
+                     deviation of the signal (independent of the sky) is
+                     `sqrt(ss)'. Therefore the total variance of this pixel
+                     is the variance of the sky added with the absolute
+                     value of its sky-subtracted flux. We use the absolute
+                     value, because especially as the signal gets noisy
+                     there will be negative values, and we don't want them
+                     to decrease the variance. */
+                  oi[ OCOL_SUM_VAR ] += (p->variance ? *ST : st*st)+fabs(ss);
+
+                  /* Get the necessary clump information. */
                   if(p->clumps && *C>0)
                     {
                       oi[ OCOL_C_NUM ]++;
@@ -260,12 +275,12 @@ mkcatalog_second_pass(struct mkcatalog_passparams *pp)
   struct mkcatalogparams *p=pp->p;
   size_t ndim=p->input->ndim, *dsize=p->input->dsize;
 
-  double *ci;
+  double *ci, *cir;
   int32_t *O, *C=NULL, nlab, *ngblabs;
   size_t i, ii, d, increment=0, num_increment=1;
   size_t nngb=gal_dimension_num_neighbors(ndim);
   size_t *dinc=gal_dimension_increment(ndim, dsize);
-  float ss, *I, *II, *SK, *ST, *input=p->input->array;
+  float ss, st, *I, *II, *SK, *ST, *input=p->input->array;
   int32_t *objects=p->objects->array, *clumps=p->clumps->array;
   size_t *c=gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim, __func__, "c");
   size_t *sc=gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim, __func__, "sc");
@@ -317,14 +332,17 @@ mkcatalog_second_pass(struct mkcatalog_passparams *pp)
 
                   /* Only use pixels above the threshold, see explanations in
                      first pass for an explanation. */
+                  st = p->variance ? sqrt(*ST) : *ST;
                   if( !( p->hasblank && isnan(*I) )
-                      && !( (ss = *I - *SK) < p->threshold * *ST ) )
+                      && !( (ss = *I - *SK) < p->threshold * st ) )
                     {
                       /* Fill in the necessary information. */
-                      ci[ CCOL_NUM    ]++;
-                      ci[ CCOL_SUM    ] += ss;
-                      ci[ CCOL_SUMSKY ] += *SK;
-                      ci[ CCOL_SUMSTD ] += *ST;
+                      ci[ CCOL_NUM     ]++;
+                      ci[ CCOL_SUM     ] += ss;
+                      ci[ CCOL_SUMSKY  ] += *SK;
+                      ci[ CCOL_SUMSTD  ] += st;
+                      ci[ CCOL_SUM_VAR ] += ( (p->variance ? *ST : st*st)
+                                              + fabs(ss) );
                       if( ss > 0.0f )
                         {
                           ci[ CCOL_NUMWHT ]++;
@@ -373,12 +391,19 @@ mkcatalog_second_pass(struct mkcatalog_passparams *pp)
                            /* It hasn't been considered yet: */
                            if(i==ii)
                              {
+                               /* Make sure it won't be considered any
+                                  more. */
                                ngblabs[ii++] = nlab;
 
-                               ++(pp->ci)[ (nlab-1) * CCOL_NUMCOLS
-                                           + CCOL_RIV_NUM ];
-                               pp->ci[ (nlab-1) * CCOL_NUMCOLS
-                                       + CCOL_RIV_SUM ] += *I-*SK;
+                               /* To help in reading. */
+                               cir=&pp->ci[ (nlab-1) * CCOL_NUMCOLS ];
+
+                               /* Write in the values. */
+                               cir[ CCOL_RIV_NUM ]++;
+                               cir[ CCOL_RIV_SUM ] += *I-*SK;
+                               cir[ CCOL_RIV_SUM_VAR ] +=
+                                 ( (p->variance ? *ST : *ST * *ST)
+                                   + fabs(*I-*SK) );
                              }
                          }
                      });
@@ -416,7 +441,7 @@ mkcatalog_median_pass(struct mkcatalog_passparams *pp)
   gal_data_t *median;
   int32_t *O, *C=NULL;
   gal_data_t **clumpsmed=NULL;
-  float ss, *I, *II, *SK, *ST;
+  float ss, st, *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,
@@ -466,9 +491,10 @@ mkcatalog_median_pass(struct mkcatalog_passparams *pp)
           /* 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. */
+          st = p->variance ? sqrt(*ST) : *ST;
           if( *O==pp->object
               && !( p->hasblank && isnan(*I) )
-              && !( (ss = *I - *SK) < p->threshold * *ST ))
+              && !( (ss = *I - *SK) < p->threshold * st ))
             {
               /* Copy the value for the whole object. */
               ss = *I - *SK;
@@ -1071,6 +1097,16 @@ mkcatalog_write_outputs(struct mkcatalogparams *p)
                       "CLUMPS");
       gal_list_str_free(comments, 1);
     }
+
+  /* Inform the user. */
+  if(p->clumpsout==p->objectsout)
+    printf("  - Output catalog: %s\n", p->objectsout);
+  else
+    {
+      printf("  - Output objects catalog: %s\n", p->objectsout);
+      if(p->clumps)
+        printf("  - Output clumps catalog: %s\n", p->clumpsout);
+    }
 }
 
 
@@ -1098,16 +1134,6 @@ mkcatalog_write_outputs(struct mkcatalogparams *p)
 void
 mkcatalog(struct mkcatalogparams *p)
 {
-  float *f, *fp;
-
-  /* If the given standard deviation array is actually variance, then take
-     its square root. */
-  if(p->variance)
-    {
-      fp = (f=p->std->array) + p->std->size;
-      do *f=sqrt(*f); while(++f<fp);
-    }
-
   /* When more than one thread is to be used, initialize the mutex: we need
      it to assign a column to the clumps in the final catalog. */
   if( p->cp.numthreads > 1 ) pthread_mutex_init(&p->mutex, NULL);
diff --git a/bin/mkcatalog/ui.h b/bin/mkcatalog/ui.h
index b8ef3e5..1f3b480 100644
--- a/bin/mkcatalog/ui.h
+++ b/bin/mkcatalog/ui.h
@@ -114,6 +114,7 @@ enum option_keys_enum
   UI_KEY_CLUMPSW2,
   UI_KEY_CLUMPSGEOW1,
   UI_KEY_CLUMPSGEOW2,
+  UI_KEY_BRIGHTNESSERR,
   UI_KEY_CLUMPSBRIGHTNESS,
   UI_KEY_NORIVERBRIGHTNESS,
   UI_KEY_MEAN,
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 254706b..2df7174 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -15315,10 +15315,19 @@ The HDU of the Sky value image.
 
 @item -t STR
 @itemx --stdfile=STR
-File name of image keeping the Sky value standard deviation for each
+File name of image keeping the @emph{Sky value} standard deviation for each
 pixel. With the @option{--variance} option you can tell MakeCatalog to
 interpret this image as a variance image, not standard deviation.
 
address@hidden note:} This must only be the Sky standard deviation or
+variance, without the signal's contribution. In other words, for each pixel
+that harbors signal, the final standard deviation is more than those
+without signal. MakeCatalog will find the amount of signal within each
+pixel using the Sky image and account for the extra error that it
+needs. Therefore if the input standard deviation (or variance) image also
+contains the contribution of signal to the error, then the final error
+measurements will be over-estimated.
+
 @item --stdhdu=STR
 The HDU of the Sky value standard deviation image.
 



reply via email to

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