[gnuastro-commits] master 3c3dc91: MakeCatalog: new --areahalfsum option

[Mohammad Akhlaghi]

branch: master
commit 3c3dc91428f4bc9ae1d13c38e9bbd3d3063527e6
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    MakeCatalog: new --areahalfsum option
    
    With this new option, it is now possible to see how many pixels constitute
    half the total brightness (sum of values) within a labeled region (object
    or clump). This can be used to get a zero-th order approximation of the
    effective radius for example. But should be used with care on "objects"
    (which may have multiple peaks), it is perfect for "clumps" though (which
    only have one peak).
    
    This was proposed by Ignacio Trujillo.
    
    This commit completes task #15786.
---
 NEWS                         |  1 +
 bin/mkcatalog/args.h         | 14 ++++++++++++++
 bin/mkcatalog/columns.c      | 23 +++++++++++++++++++++++
 bin/mkcatalog/main.h         |  8 +++++---
 bin/mkcatalog/mkcatalog.c    |  1 +
 bin/mkcatalog/parse.c        | 37 +++++++++++++++++++++++++++++++++++++
 bin/mkcatalog/ui.c           |  4 +++-
 bin/mkcatalog/ui.h           |  1 +
 doc/announce-acknowledge.txt |  1 +
 doc/gnuastro.texi            |  8 ++++++++
 10 files changed, 94 insertions(+), 4 deletions(-)

diff --git a/NEWS b/NEWS
index 32fe210..a33b959 100644
--- a/NEWS
+++ b/NEWS
@@ -42,6 +42,7 @@ See the end of the file for license conditions.
      within the image.
 
   MakeCatalog:
+   --areahalfsum: area containing half the total flux of object or clump.
    - New columns to return the position of pixel with minimum or maximum
      value: '--minvx', '--maxvx', '--minvy', '--maxvy', '--minvz',
      '--maxvz'.
diff --git a/bin/mkcatalog/args.h b/bin/mkcatalog/args.h
index 03018df..d05cdfd 100644
--- a/bin/mkcatalog/args.h
+++ b/bin/mkcatalog/args.h
@@ -1422,6 +1422,20 @@ struct argp_option program_options[] =
       ui_column_codes_ll
     },
     {
+      "areahalfsum",
+      UI_KEY_AREAHALFSUM,
+      0,
+      0,
+      "",
+      UI_GROUP_COLUMNS_MORPHOLOGY,
+      0,
+      GAL_TYPE_INVALID,
+      GAL_OPTIONS_RANGE_ANY,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET,
+      ui_column_codes_ll
+    },
+    {
       "clumpsarea",
       UI_KEY_CLUMPSAREA,
       0,
diff --git a/bin/mkcatalog/columns.c b/bin/mkcatalog/columns.c
index f79d5a5..a201b57 100644
--- a/bin/mkcatalog/columns.c
+++ b/bin/mkcatalog/columns.c
@@ -450,6 +450,21 @@ columns_define_alloc(struct mkcatalogparams *p)
           oiflag[ OCOL_NUMXY ] = ciflag[ CCOL_NUMXY ] = 1;
           break;
 
+        case UI_KEY_AREAHALFSUM:
+          name           = "AREA_HALF_SUM";
+          unit           = "counter";
+          ocomment       = "Number of brightest pixels containing half of 
total sum.";
+          ccomment       = ocomment;
+          otype          = GAL_TYPE_INT32;
+          ctype          = GAL_TYPE_INT32;
+          disp_fmt       = 0;
+          disp_width     = 6;
+          disp_precision = 0;
+          oiflag[ OCOL_NUM        ] = ciflag[ CCOL_NUM        ] = 1;
+          oiflag[ OCOL_SUM        ] = ciflag[ CCOL_SUM        ] = 1;
+          oiflag[ OCOL_NUMHALFSUM ] = ciflag[ CCOL_NUMHALFSUM ] = 1;
+          break;
+
         case UI_KEY_CLUMPSAREA:
           name           = "AREA_CLUMPS";
           unit           = "counter";
@@ -2039,6 +2054,10 @@ columns_fill(struct mkcatalog_passparams *pp)
           ((int32_t *)colarr)[oind] = oi[OCOL_NUMXY];
           break;
 
+        case UI_KEY_AREAHALFSUM:
+          ((int32_t *)colarr)[oind] = oi[OCOL_NUMHALFSUM];
+          break;
+
         case UI_KEY_CLUMPSAREA:
           ((int32_t *)colarr)[oind] = oi[OCOL_C_NUM];
           break;
@@ -2346,6 +2365,10 @@ columns_fill(struct mkcatalog_passparams *pp)
             ((int32_t *)colarr)[cind]=ci[CCOL_NUMXY];
             break;
 
+          case UI_KEY_AREAHALFSUM:
+            ((int32_t *)colarr)[cind]=ci[CCOL_NUMHALFSUM];
+            break;
+
           case UI_KEY_WEIGHTAREA:
             ((int32_t *)colarr)[cind]=ci[CCOL_NUMWHT];
             break;
diff --git a/bin/mkcatalog/main.h b/bin/mkcatalog/main.h
index 195f004..9298a45 100644
--- a/bin/mkcatalog/main.h
+++ b/bin/mkcatalog/main.h
@@ -70,10 +70,11 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
    the FITS standard (fastest dimension is first). */
 enum objectcols
   {
-    OCOL_NUMALL,         /* Number of all pixels with this label.     */
-    OCOL_NUMALLXY,       /* Number of all pixels in first two dims.   */
-    OCOL_NUM,            /* Number of values used in this object.     */
+    OCOL_NUMALL,         /* Area/number of all pixels with this label.*/
+    OCOL_NUMALLXY,       /* Area/Number in first two dimensions.      */
+    OCOL_NUM,            /* Area/Number of values used in this object.*/
     OCOL_NUMXY,          /* Number of values in the first two dims.   */
+    OCOL_NUMHALFSUM,     /* Area/Number containing half of total sum. */
     OCOL_SUM,            /* Sum of (value-sky) in object.             */
     OCOL_SUM_VAR,        /* Variance including values (not just sky). */
     OCOL_MEDIAN,         /* Median of value in object.                */
@@ -130,6 +131,7 @@ enum clumpcols
     CCOL_NUMALLXY,       /* Number of pixels in first two dims.       */
     CCOL_NUM,            /* Number of values used in clump.           */
     CCOL_NUMXY,          /* Number of values only in first two dims.  */
+    CCOL_NUMHALFSUM,     /* Area/Number containing half of total sum. */
     CCOL_SUM,            /* River subtracted brightness.              */
     CCOL_SUM_VAR,        /* Variance including values (not just sky). */
     CCOL_MEDIAN,         /* Median of values in clump.                */
diff --git a/bin/mkcatalog/mkcatalog.c b/bin/mkcatalog/mkcatalog.c
index 1f42887..3204f2c 100644
--- a/bin/mkcatalog/mkcatalog.c
+++ b/bin/mkcatalog/mkcatalog.c
@@ -182,6 +182,7 @@ mkcatalog_single_object(void *in_prm)
       /* If an order-based calculation is requested, another pass is
          necessary. */
       if( p->oiflag[ OCOL_MEDIAN ]
+          || p->oiflag[ OCOL_NUMHALFSUM ]
           || p->oiflag[ OCOL_SIGCLIPNUM ]
           || p->oiflag[ OCOL_SIGCLIPSTD ]
           || p->oiflag[ OCOL_SIGCLIPMEAN ]
diff --git a/bin/mkcatalog/parse.c b/bin/mkcatalog/parse.c
index f1ea95c..ca2b78f 100644
--- a/bin/mkcatalog/parse.c
+++ b/bin/mkcatalog/parse.c
@@ -1182,6 +1182,35 @@ parse_clumps(struct mkcatalog_passparams *pp)
 
 
 
+size_t
+parse_area_of_half_sum(gal_data_t *values, double sum)
+{
+  size_t i;
+  gal_data_t *sorted_d;
+  double sumcheck=0.0f, *sorted;
+
+  /* Allocate the array to use. */
+  sorted_d = ( values->type==GAL_TYPE_FLOAT64
+                ? values
+                : gal_data_copy_to_new_type(values, GAL_TYPE_FLOAT64) );
+
+  /* Sort the desired labels and find the number of elements where we reach
+     half the total sum. */
+  gal_statistics_sort_decreasing(sorted_d);
+  sorted=sorted_d->array;
+  for(i=0;i<values->size;++i)
+    if( (sumcheck+=sorted[i]) > sum/2 )
+      break;
+
+  /* Clean up and return. */
+  if(sorted_d!=values) gal_data_free(sorted_d);
+  return i;
+}
+
+
+
+
+
 void
 parse_order_based(struct mkcatalog_passparams *pp)
 {
@@ -1203,6 +1232,7 @@ parse_order_based(struct mkcatalog_passparams *pp)
   if(tmpsize==0)
     {
       if(p->oiflag[ OCOL_MEDIAN        ]) pp->oi[ OCOL_MEDIAN       ] = NAN;
+      if(p->oiflag[ OCOL_NUMHALFSUM    ]) pp->oi[ OCOL_NUMHALFSUM   ] = 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;
@@ -1212,6 +1242,7 @@ parse_order_based(struct mkcatalog_passparams *pp)
           {
             ci=&pp->ci[ i * CCOL_NUMCOLS ];
             if(p->ciflag[ CCOL_MEDIAN        ]) ci[ CCOL_MEDIAN      ] = NAN;
+            if(p->ciflag[ OCOL_NUMHALFSUM    ]) ci[ CCOL_NUMHALFSUM  ] = 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;
@@ -1321,6 +1352,8 @@ parse_order_based(struct mkcatalog_passparams *pp)
       /* Clean up the sigma-clipped values. */
       gal_data_free(result);
     }
+  if(p->oiflag[ OCOL_NUMHALFSUM ])
+    pp->oi[OCOL_NUMHALFSUM]=parse_area_of_half_sum(objvals, pp->oi[OCOL_SUM]);
 
 
   /* Clean up the object values. */
@@ -1370,6 +1403,10 @@ parse_order_based(struct mkcatalog_passparams *pp)
               gal_data_free(result);
             }
 
+          /* Estimate half of the total sum. */
+          if(p->ciflag[ CCOL_NUMHALFSUM ])
+            ci[CCOL_NUMHALFSUM]=parse_area_of_half_sum(clumpsvals[i], 
ci[CCOL_SUM]);
+
           /* Clean up this clump's values. */
           gal_data_free(clumpsvals[i]);
         }
diff --git a/bin/mkcatalog/ui.c b/bin/mkcatalog/ui.c
index d9f5b86..ae5dd34 100644
--- a/bin/mkcatalog/ui.c
+++ b/bin/mkcatalog/ui.c
@@ -929,6 +929,7 @@ ui_necessary_inputs(struct mkcatalogparams *p, int *values, 
int *sky,
         case OCOL_NUMALLXY:           /* Only object labels. */    break;
         case OCOL_NUM:                *values        = 1;          break;
         case OCOL_NUMXY:              *values        = 1;          break;
+        case CCOL_NUMHALFSUM:         *values        = 1;          break;
         case OCOL_SUM:                *values        = 1;          break;
         case OCOL_SUM_VAR:            *values = *std = 1;          break;
         case OCOL_MEDIAN:             *values        = 1;          break;
@@ -990,10 +991,11 @@ ui_necessary_inputs(struct mkcatalogparams *p, int 
*values, int *sky,
           case CCOL_NUMALLXY:         /* Only clump labels. */     break;
           case CCOL_NUM:              *values        = 1;          break;
           case CCOL_NUMXY:            *values        = 1;          break;
+          case CCOL_NUMHALFSUM:       *values        = 1;          break;
           case CCOL_SUM:              *values        = 1;          break;
           case CCOL_SUM_VAR:          *values = *std = 1;          break;
           case CCOL_MEDIAN:           *values        = 1;          break;
-          case CCOL_SIGCLIPNUM:    *values        = 1;          break;
+          case CCOL_SIGCLIPNUM:       *values        = 1;          break;
           case CCOL_SIGCLIPMEDIAN:    *values        = 1;          break;
           case CCOL_SIGCLIPMEAN:      *values        = 1;          break;
           case CCOL_SIGCLIPSTD:       *values        = 1;          break;
diff --git a/bin/mkcatalog/ui.h b/bin/mkcatalog/ui.h
index 9eb98d7..155801f 100644
--- a/bin/mkcatalog/ui.h
+++ b/bin/mkcatalog/ui.h
@@ -107,6 +107,7 @@ enum option_keys_enum
   UI_KEY_OBJID,                         /* Catalog columns. */
   UI_KEY_IDINHOSTOBJ,
   UI_KEY_AREAXY,
+  UI_KEY_AREAHALFSUM,
   UI_KEY_CLUMPSAREA,
   UI_KEY_WEIGHTAREA,
   UI_KEY_GEOAREA,
diff --git a/doc/announce-acknowledge.txt b/doc/announce-acknowledge.txt
index 40919c8..2b707cb 100644
--- a/doc/announce-acknowledge.txt
+++ b/doc/announce-acknowledge.txt
@@ -5,6 +5,7 @@ Samane Raji
 Joanna Sakowska
 Zahra Sharbaf
 Sachin Kumar Singh
+Ignacio Trujillo
 
 
 
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 6302f3a..7fef1b2 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -16265,6 +16265,14 @@ 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 --areahalfsum
+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.
+
+This option is thus good for clumps (which are defined to have a single peak 
in their morphology), but for objects you should be careful because if they may 
have multiple peaks and the sorting by value doesn't include morphology.
+So if the object includes multiple peaks/clumps at roughly the same level, 
then the area reported by this option will be distributed over all the peaks.
+
+
 @item --clumpsarea
 [Objects] The total area of all the clumps in this object.