[gnuastro-commits] master 34c1d8e: MakeCatalog: new option to estimate the observed FWHM

[Mohammad Akhlaghi]

branch: master
commit 34c1d8ed27d05355d390c0a29f7a74ea0a0f1f27
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    MakeCatalog: new option to estimate the observed FWHM
    
    The newly added half-sum and frac-sum radius columns of MakeCatalog depend
    on the total flux, which is hard to estimate without deblending. So using
    them had limitations.
    
    However, the FWHM is much more robust to deblending (at least when the
    deblending level is below half the maximum!), so it is now added following
    a similar formalism as the half-sum and frac-sum columns (using the area).
---
 NEWS                    |  1 +
 bin/mkcatalog/args.h    | 14 ++++++++++++++
 bin/mkcatalog/columns.c | 48 ++++++++++++++++++++++++++++++++++--------------
 bin/mkcatalog/main.h    |  2 ++
 bin/mkcatalog/parse.c   | 41 +++++++++++++++++++++++++++++++++--------
 bin/mkcatalog/ui.c      |  2 ++
 bin/mkcatalog/ui.h      |  1 +
 doc/gnuastro.texi       | 10 ++++++++++
 8 files changed, 97 insertions(+), 22 deletions(-)

diff --git a/NEWS b/NEWS
index 3413ef2..9ddfb75 100644
--- a/NEWS
+++ b/NEWS
@@ -42,6 +42,7 @@ See the end of the file for license conditions.
      within the image.
 
   MakeCatalog:
+   --fwhmobs: The observed FWHM in pixels, along the major axis.
    --fracsum: fractions to use in '--fracsumarea1' or '--fracsumarea2'.
    --fracsumarea1: area of given fraction of the summed object or clump values.
    --fracsumarea2: area of given fraction of the summed object or clump values.
diff --git a/bin/mkcatalog/args.h b/bin/mkcatalog/args.h
index 72a8c9d..c534152 100644
--- a/bin/mkcatalog/args.h
+++ b/bin/mkcatalog/args.h
@@ -1617,6 +1617,20 @@ struct argp_option program_options[] =
       ui_column_codes_ll
     },
     {
+      "fwhmobs",
+      UI_KEY_FWHMOBS,
+      0,
+      0,
+      "Full width at half max (non-parametric).",
+      UI_GROUP_COLUMNS_MORPHOLOGY,
+      0,
+      GAL_TYPE_INVALID,
+      GAL_OPTIONS_RANGE_ANY,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET,
+      ui_column_codes_ll
+    },
+    {
       "halfsumarea",
       UI_KEY_HALFSUMAREA,
       0,
diff --git a/bin/mkcatalog/columns.c b/bin/mkcatalog/columns.c
index c9ba603..49f8cfa 100644
--- a/bin/mkcatalog/columns.c
+++ b/bin/mkcatalog/columns.c
@@ -1677,6 +1677,7 @@ columns_define_alloc(struct mkcatalogparams *p)
             oiflag[ OCOL_NUMFRACSUM2 ] = ciflag[ CCOL_NUMFRACSUM2 ] = 1;
           break;
 
+        case UI_KEY_FWHMOBS:
         case UI_KEY_HALFRADIUSOBS:
         case UI_KEY_FRACRADIUSOBS1:
         case UI_KEY_FRACRADIUSOBS2:
@@ -1702,6 +1703,11 @@ columns_define_alloc(struct mkcatalogparams *p)
           oiflag[ OCOL_GXY        ] = ciflag[ CCOL_GXY        ] = 1;
           switch(colcode->v)
             {
+            case UI_KEY_FWHMOBS:
+              name="FWHM_OBS";
+              oiflag[ OCOL_NUMHALFMAX  ] = ciflag[ CCOL_NUMHALFMAX  ] = 1;
+              ocomment = "Full width at half maximum (accounting for 
ellipticity).";
+              break;
             case UI_KEY_HALFRADIUSOBS:
               name="HALF_RADIUS_OBS";
               oiflag[ OCOL_NUMHALFSUM  ] = ciflag[ CCOL_NUMHALFSUM  ] = 1;
@@ -2084,9 +2090,9 @@ columns_fill(struct mkcatalog_passparams *pp)
   int key;
   double tmp;
   void *colarr;
-  size_t tmpind;
   gal_data_t *column;
   double *ci, *oi=pp->oi;
+  size_t tmpind=GAL_BLANK_SIZE_T;
   size_t coord[3]={GAL_BLANK_SIZE_T, GAL_BLANK_SIZE_T, GAL_BLANK_SIZE_T};
 
   size_t i, cind, coind, sr=pp->clumpstartindex, oind=GAL_BLANK_SIZE_T;
@@ -2414,23 +2420,34 @@ columns_fill(struct mkcatalog_passparams *pp)
           break;
 
         case UI_KEY_HALFSUMSB:
+          tmpind = ( key==UI_KEY_HALFRADIUSOBS
+                     ? OCOL_NUMHALFSUM
+                     : ( key==UI_KEY_FRACRADIUSOBS1
+                         ? OCOL_NUMFRACSUM1
+                         : OCOL_NUMFRACSUM2 ) );
           ((float *)colarr)[oind] = oi[OCOL_SUM]/2.0f/oi[OCOL_NUMHALFSUM];
           break;
 
+        case UI_KEY_FWHMOBS:
         case UI_KEY_HALFRADIUSOBS:
         case UI_KEY_FRACRADIUSOBS1:
         case UI_KEY_FRACRADIUSOBS2:
-          /* First derive the axis ratio (as 'tmp'), then use the
-             '--halfsumarea' to estimate the effective radius. */
+          /* First derive the axis ratio (as 'tmp'), then set the index to
+             use and calculate the radius from the area and axis ratio. */
           tmp = ( columns_second_order(pp, oi, UI_KEY_SEMIMINOR, 0)
                   / columns_second_order(pp, oi, UI_KEY_SEMIMAJOR, 0) );
-          tmpind = ( key==UI_KEY_HALFRADIUSOBS
-                     ? OCOL_NUMHALFSUM
-                     : ( key==UI_KEY_FRACRADIUSOBS1
-                         ? OCOL_NUMFRACSUM1
-                         : OCOL_NUMFRACSUM2 ) );
+          switch(key)
+            {
+            case UI_KEY_FWHMOBS:        tmpind=OCOL_NUMHALFMAX;  break;
+            case UI_KEY_HALFRADIUSOBS:  tmpind=OCOL_NUMHALFSUM;  break;
+            case UI_KEY_FRACRADIUSOBS1: tmpind=OCOL_NUMFRACSUM1; break;
+            case UI_KEY_FRACRADIUSOBS2: tmpind=OCOL_NUMFRACSUM2; break;
+            }
           tmp = sqrt( oi[tmpind]/(tmp*M_PI) );
-          ((float *)colarr)[oind] = tmp<1e-6 ? NAN : tmp;
+          if(key==UI_KEY_FWHMOBS)
+            ((float *)colarr)[oind] = (tmp<1e-6 ? NAN : tmp)*2;
+          else
+            ((float *)colarr)[oind] = tmp<1e-6 ? NAN : tmp;
           break;
 
         default:
@@ -2703,16 +2720,19 @@ columns_fill(struct mkcatalog_passparams *pp)
             ((float *)colarr)[cind] = ci[CCOL_SUM]/2.0f/ci[CCOL_NUMHALFSUM];
             break;
 
+          case UI_KEY_FWHMOBS:
           case UI_KEY_HALFRADIUSOBS:
           case UI_KEY_FRACRADIUSOBS1:
           case UI_KEY_FRACRADIUSOBS2:
             tmp = ( columns_second_order(  pp, ci, UI_KEY_SEMIMINOR, 1)
                     / columns_second_order(pp, ci, UI_KEY_SEMIMAJOR, 1) );
-            tmpind = ( key==UI_KEY_HALFRADIUSOBS
-                       ? CCOL_NUMHALFSUM
-                       : ( key==UI_KEY_FRACRADIUSOBS1
-                           ? CCOL_NUMFRACSUM1
-                           : CCOL_NUMFRACSUM2 ) );
+            switch(key)
+              {
+              case UI_KEY_FWHMOBS:        tmpind=CCOL_NUMHALFMAX;  break;
+              case UI_KEY_HALFRADIUSOBS:  tmpind=CCOL_NUMHALFSUM;  break;
+              case UI_KEY_FRACRADIUSOBS1: tmpind=CCOL_NUMFRACSUM1; break;
+              case UI_KEY_FRACRADIUSOBS2: tmpind=CCOL_NUMFRACSUM2; break;
+              }
             tmp = sqrt( ci[tmpind]/(tmp*M_PI) );
             ((float *)colarr)[cind] = tmp<1e-6 ? NAN : tmp;
             break;
diff --git a/bin/mkcatalog/main.h b/bin/mkcatalog/main.h
index d8d1c25..54b5e42 100644
--- a/bin/mkcatalog/main.h
+++ b/bin/mkcatalog/main.h
@@ -109,6 +109,7 @@ enum objectcols
     OCOL_UPPERLIMIT_S,   /* Upper limit one-sigma value.              */
     OCOL_UPPERLIMIT_Q,   /* Quantile of object in random distribution.*/
     OCOL_UPPERLIMIT_SKEW,/* (Mean-Median)/STD of random distribution. */
+    OCOL_NUMHALFMAX,     /* Area/Number of pixels above half of max.  */
     OCOL_NUMHALFSUM,     /* Area/Number containing half of total sum. */
     OCOL_NUMFRACSUM1,    /* Area/Number containing frac of total sum. */
     OCOL_NUMFRACSUM2,    /* Area/Number containing frac of total sum. */
@@ -177,6 +178,7 @@ enum clumpcols
     CCOL_UPPERLIMIT_S,   /* Upper limit one-sigma value.              */
     CCOL_UPPERLIMIT_Q,   /* Quantile of object in random distribution.*/
     CCOL_UPPERLIMIT_SKEW,/* (Mean-Median)/STD of random distribution. */
+    CCOL_NUMHALFMAX,     /* Area/Number of pixels above half of max.  */
     CCOL_NUMHALFSUM,     /* Area/Number containing half of total sum. */
     CCOL_NUMFRACSUM1,    /* Area/Number containing frac of total sum. */
     CCOL_NUMFRACSUM2,    /* Area/Number containing frac of total sum. */
diff --git a/bin/mkcatalog/parse.c b/bin/mkcatalog/parse.c
index 63a9964..5f80da1 100644
--- a/bin/mkcatalog/parse.c
+++ b/bin/mkcatalog/parse.c
@@ -1183,15 +1183,18 @@ parse_clumps(struct mkcatalog_passparams *pp)
 
 
 static double
-parse_frac_sum_find(gal_data_t *sorted_d, double sum, double frac)
+parse_frac_find(gal_data_t *sorted_d, double value, double frac, int dosum)
 {
   size_t i;
-  double sumcheck=0.0f;
+  double check=0.0f;
   double *sorted=sorted_d->array;
 
   /* Parse over the sorted array and find the index. */
   for(i=0;i<sorted_d->size;++i)
-    if( (sumcheck+=sorted[i]) > sum*frac ) break;
+    if(dosum)
+      { if( (check+=sorted[i]) > value*frac ) break; }
+    else
+      { if(         sorted[i]  < value*frac ) break; }
 
   /* Return the final value. Note that if the index is zero, we should
      actually return 1, because we are starting with the maximum. */
@@ -1208,6 +1211,7 @@ parse_area_of_frac_sum(struct mkcatalog_passparams *pp, 
gal_data_t *values,
 {
   struct mkcatalogparams *p=pp->p;
 
+  double max, *sorted;
   gal_data_t *sorted_d;
   uint8_t *flag = o1c0 ? p->oiflag : p->ciflag;
   double sum = o1c0 ? outarr[OCOL_SUM] : outarr[CCOL_SUM];
@@ -1225,15 +1229,28 @@ parse_area_of_frac_sum(struct mkcatalog_passparams *pp, 
gal_data_t *values,
   /* Set the required fractions. */
   if(flag[ o1c0 ? OCOL_NUMHALFSUM : CCOL_NUMHALFSUM ])
     outarr[ o1c0 ? OCOL_NUMHALFSUM : CCOL_NUMHALFSUM ]
-      = parse_frac_sum_find(sorted_d, sum, 0.5f);
+      = parse_frac_find(sorted_d, sum, 0.5f, 1);
 
   if(flag[ o1c0 ? OCOL_NUMFRACSUM1 : CCOL_NUMFRACSUM1 ])
     outarr[ o1c0 ? OCOL_NUMFRACSUM1 : CCOL_NUMFRACSUM1 ]
-      = parse_frac_sum_find(sorted_d, sum, fracsum[0]);
+      = parse_frac_find(sorted_d, sum, fracsum[0], 1);
 
   if(flag[ o1c0 ? OCOL_NUMFRACSUM2 : CCOL_NUMFRACSUM2 ])
     outarr[ o1c0 ? OCOL_NUMFRACSUM2 : CCOL_NUMFRACSUM2 ]
-      = parse_frac_sum_find(sorted_d, sum, fracsum[1]);
+      = parse_frac_find(sorted_d, sum, fracsum[1], 1);
+
+  /* For the FWHM, we'll use the median of the top three pixels for the
+     maximum (to avoid noise). */
+  if(flag[ o1c0 ? OCOL_NUMHALFMAX : CCOL_NUMHALFMAX ])
+    {
+      sorted=sorted_d->array;
+      max = ( sorted_d->size>3
+              ? (sorted[0]+sorted[1]+sorted[2])/3
+              : sorted[0] );
+      outarr[ o1c0 ? OCOL_NUMHALFMAX : CCOL_NUMHALFMAX ]
+        = parse_frac_find(sorted_d, max, 0.5f, 0);
+      printf("area: %f (%f)\n", max, outarr[OCOL_NUMHALFMAX]);
+    }
 
   /* For a check.
   printf("%g, %g, %g\n", outarr[OCOL_NUMFRACSUM1], outarr[OCOL_NUMHALFSUM],
@@ -1270,7 +1287,10 @@ parse_order_based(struct mkcatalog_passparams *pp)
   if(tmpsize==0)
     {
       if(p->oiflag[ OCOL_MEDIAN        ]) pp->oi[ OCOL_MEDIAN       ] = NAN;
+      if(p->oiflag[ OCOL_NUMHALFMAX    ]) pp->oi[ OCOL_NUMHALFMAX   ] = 0;
       if(p->oiflag[ OCOL_NUMHALFSUM    ]) pp->oi[ OCOL_NUMHALFSUM   ] = 0;
+      if(p->oiflag[ OCOL_NUMFRACSUM1   ]) pp->oi[ OCOL_NUMFRACSUM1  ] = 0;
+      if(p->oiflag[ OCOL_NUMFRACSUM2   ]) pp->oi[ OCOL_NUMFRACSUM2  ] = 0;
       if(p->oiflag[ OCOL_SIGCLIPNUM    ]) pp->oi[ OCOL_SIGCLIPNUM   ] = 0;
       if(p->oiflag[ OCOL_SIGCLIPSTD    ]) pp->oi[ OCOL_SIGCLIPSTD   ] = 0;
       if(p->oiflag[ OCOL_SIGCLIPMEAN   ]) pp->oi[ OCOL_SIGCLIPMEAN  ] = NAN;
@@ -1280,7 +1300,10 @@ parse_order_based(struct mkcatalog_passparams *pp)
           {
             ci=&pp->ci[ i * CCOL_NUMCOLS ];
             if(p->ciflag[ CCOL_MEDIAN        ]) ci[ CCOL_MEDIAN      ] = NAN;
+            if(p->oiflag[ OCOL_NUMHALFMAX    ]) ci[ OCOL_NUMHALFMAX  ] = 0;
             if(p->ciflag[ OCOL_NUMHALFSUM    ]) ci[ CCOL_NUMHALFSUM  ] = 0;
+            if(p->oiflag[ OCOL_NUMFRACSUM1   ]) ci[ OCOL_NUMFRACSUM1 ] = 0;
+            if(p->oiflag[ OCOL_NUMFRACSUM2   ]) ci[ OCOL_NUMFRACSUM2 ] = 0;
             if(p->ciflag[ CCOL_SIGCLIPNUM    ]) ci[ CCOL_SIGCLIPNUM  ] = 0;
             if(p->ciflag[ CCOL_SIGCLIPSTD    ]) ci[ CCOL_SIGCLIPSTD  ] = 0;
             if(p->ciflag[ CCOL_SIGCLIPMEAN   ]) ci[ CCOL_SIGCLIPMEAN ] = NAN;
@@ -1392,7 +1415,8 @@ parse_order_based(struct mkcatalog_passparams *pp)
     }
 
   /* Fractional values. */
-  if( p->oiflag[    OCOL_NUMHALFSUM   ]
+  if( p->oiflag[    OCOL_NUMHALFMAX  ]
+      || p->oiflag[ OCOL_NUMHALFSUM  ]
       || p->oiflag[ OCOL_NUMFRACSUM1 ]
       || p->oiflag[ OCOL_NUMFRACSUM2 ] )
     parse_area_of_frac_sum(pp, objvals, pp->oi, 1);
@@ -1445,7 +1469,8 @@ parse_order_based(struct mkcatalog_passparams *pp)
             }
 
           /* Estimate half of the total sum. */
-          if( p->ciflag[ CCOL_NUMHALFSUM ]
+          if( p->ciflag[    CCOL_NUMHALFMAX ]
+              || p->ciflag[ CCOL_NUMHALFSUM ]
               || p->ciflag[ CCOL_NUMFRACSUM1 ]
               || p->ciflag[ CCOL_NUMFRACSUM2 ] )
             parse_area_of_frac_sum(pp, clumpsvals[i], ci, 0);
diff --git a/bin/mkcatalog/ui.c b/bin/mkcatalog/ui.c
index 501d043..c9ed56f 100644
--- a/bin/mkcatalog/ui.c
+++ b/bin/mkcatalog/ui.c
@@ -996,6 +996,7 @@ ui_necessary_inputs(struct mkcatalogparams *p, int *values, 
int *sky,
         case OCOL_UPPERLIMIT_S:       *values        = 1;          break;
         case OCOL_UPPERLIMIT_Q:       *values        = 1;          break;
         case OCOL_UPPERLIMIT_SKEW:    *values        = 1;          break;
+        case OCOL_NUMHALFMAX:         *values        = 1;          break;
         case OCOL_NUMHALFSUM:         *values        = 1;          break;
         case OCOL_NUMFRACSUM1:        *values        = 1;          break;
         case OCOL_NUMFRACSUM2:        *values        = 1;          break;
@@ -1069,6 +1070,7 @@ ui_necessary_inputs(struct mkcatalogparams *p, int 
*values, int *sky,
           case CCOL_UPPERLIMIT_S:     *values        = 1;          break;
           case CCOL_UPPERLIMIT_Q:     *values        = 1;          break;
           case CCOL_UPPERLIMIT_SKEW:  *values        = 1;          break;
+          case CCOL_NUMHALFMAX:       *values        = 1;          break;
           case CCOL_NUMHALFSUM:       *values        = 1;          break;
           case CCOL_NUMFRACSUM1:      *values        = 1;          break;
           case CCOL_NUMFRACSUM2:      *values        = 1;          break;
diff --git a/bin/mkcatalog/ui.h b/bin/mkcatalog/ui.h
index 33b17e3..d7a2303 100644
--- a/bin/mkcatalog/ui.h
+++ b/bin/mkcatalog/ui.h
@@ -169,6 +169,7 @@ enum option_keys_enum
   UI_KEY_GEOSEMIMINOR,
   UI_KEY_GEOAXISRATIO,
   UI_KEY_GEOPOSITIONANGLE,
+  UI_KEY_FWHMOBS,
   UI_KEY_HALFSUMAREA,
   UI_KEY_HALFSUMSB,
   UI_KEY_HALFRADIUSOBS,
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index e142b88..44934e7 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -16270,6 +16270,16 @@ first two dimensions. This is only available for 
3-dimensional
 datasets. When working with Integral Field Unit (IFU) datasets, this
 projection onto the first two dimensions would be a narrow-band image.
 
+@item --fwhmobs
+@cindex FWHM
+The full width at half maximum (in units of pixels, along the semi-major axis) 
of the labeled region (object or clump).
+The maximum value is estimated from the mean of the top-three pixels with the 
highest values (to avoid noise fluctuations which can be large due to Poisson 
noise).
+The number of pixels that have half the value of that maximum are then found 
and a radius is estimated from the area.
+See the description under @option{--halfradiusobs} for more on converting area 
to radius along major axis.
+
+It is therefore important to highlight that this FWHM is not estimated from a 
parametric fit to the labeled region: it does not account for the PSF, and it 
will be strongly affected by noise.
+So when half the maximum value is too close to the noise level, the value 
returned by this function is meaningless (its just noise peaks which are 
randomly distributed over the area).
+
 @item --halfsumarea
 The number of pixels that contain half the object or clump's total sum of 
pixels (half the value in the @option{--brightness} column).
 To count this area, all the non-blank values associated with the given label 
(object or clump) will be sorted and summed in order (starting from the 
maximum), until the sum becomes larger than half the total sum of the label's 
pixels.