gnuastro-commits
[Top][All Lists]
Advanced

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

[gnuastro-commits] master 69bae96: MakeCatalog: new columns for fraction


From: Mohammad Akhlaghi
Subject: [gnuastro-commits] master 69bae96: MakeCatalog: new columns for fractional radius, area, and SB
Date: Wed, 7 Oct 2020 23:51:16 -0400 (EDT)

branch: master
commit 69bae967066a27455b503347f6622af4d6c37f86
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    MakeCatalog: new columns for fractional radius, area, and SB
    
    With the addition of the previously added '--areahalfsum' option, Ignacio
    Trujillo suggested to calculate the approximate radius that corresponds to
    the area, so the new '--halfradiusobs' option was added, and to help in
    sorting the option names with the output of '--help', the original option
    was renamed to '--halfsumarea'. Also, by dividing the flux with the region
    and its area, we now have a new '--halfsumsb' option to report the surface
    brightness.
    
    But to help in identifying the slope of objects, and in particular stars in
    an image, two other options were added for the user to optionally ask for
    areas and radii at custom fractions of total flux: '--fracsumarea1',
    '--fracsumarea2', '--fracradiusobs1' and '--fracradiusobs2'. The custom
    fractions can be set with the newly added '--fracsum' option (which can
    take two floating point values).
    
    These columns were suggested by Ignacio Trujillo.
---
 NEWS                      |   9 ++-
 bin/mkcatalog/args.h      | 142 +++++++++++++++++++++++++++++++----
 bin/mkcatalog/columns.c   | 184 ++++++++++++++++++++++++++++++++++++++++------
 bin/mkcatalog/main.h      |  12 ++-
 bin/mkcatalog/mkcatalog.c |   2 +
 bin/mkcatalog/parse.c     |  67 ++++++++++++++---
 bin/mkcatalog/ui.c        |  44 ++++++++++-
 bin/mkcatalog/ui.h        |  10 ++-
 doc/gnuastro.texi         |  45 +++++++++++-
 9 files changed, 456 insertions(+), 59 deletions(-)

diff --git a/NEWS b/NEWS
index a33b959..3413ef2 100644
--- a/NEWS
+++ b/NEWS
@@ -42,7 +42,14 @@ See the end of the file for license conditions.
      within the image.
 
   MakeCatalog:
-   --areahalfsum: area containing half the total flux of object or clump.
+   --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.
+   --halfsumsb: Surface brightness within area reported by '--halfsumarea'.
+   --halfsumarea: area containing half of the summed object or clump values.
+   --halfradiusobs: radius derived from '--halfsumarea', underestimates r_e.
+   --fracradiusobs1: similar to --halfradiusobs, but for custom fraction 1.
+   --fracradiusobs2: similar to --halfradiusobs, but for custom fraction 2.
    - 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 d05cdfd..72a8c9d 100644
--- a/bin/mkcatalog/args.h
+++ b/bin/mkcatalog/args.h
@@ -287,6 +287,8 @@ struct argp_option program_options[] =
 
 
 
+
+
     /* Upper limit magnitude configurations. */
     {
       0, 0, 0, 0,
@@ -405,6 +407,31 @@ struct argp_option program_options[] =
 
 
 
+    /* Other column configurations. */
+    {
+      0, 0, 0, 0,
+      "Settings for other columns:",
+      UI_GROUP_OTHERSETTINGS
+    },
+    {
+      "fracsum",
+      UI_KEY_FRACSUM,
+      "FLT[,FLT]",
+      0,
+      "Fraction(s) in --fracsumarea1 or --fracsumarea2.",
+      UI_GROUP_OTHERSETTINGS,
+      &p->fracsum,
+      GAL_TYPE_STRING,
+      GAL_OPTIONS_RANGE_ANY,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET,
+      gal_options_parse_csv_float64
+    },
+
+
+
+
+
     /* ID related columns. */
     {
       0, 0, 0, 0,
@@ -1422,20 +1449,6 @@ 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,
@@ -1603,6 +1616,107 @@ struct argp_option program_options[] =
       GAL_OPTIONS_NOT_SET,
       ui_column_codes_ll
     },
+    {
+      "halfsumarea",
+      UI_KEY_HALFSUMAREA,
+      0,
+      0,
+      "Area containing half of --brightness.",
+      UI_GROUP_COLUMNS_MORPHOLOGY,
+      0,
+      GAL_TYPE_INVALID,
+      GAL_OPTIONS_RANGE_ANY,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET,
+      ui_column_codes_ll
+    },
+    {
+      "halfsumsb",
+      UI_KEY_HALFSUMSB,
+      0,
+      0,
+      "Surface brightness within --halfsumarea.",
+      UI_GROUP_COLUMNS_MORPHOLOGY,
+      0,
+      GAL_TYPE_INVALID,
+      GAL_OPTIONS_RANGE_ANY,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET,
+      ui_column_codes_ll
+    },
+    {
+      "halfradiusobs",
+      UI_KEY_HALFRADIUSOBS,
+      0,
+      0,
+      "Radius calculated from --halfsumarea.",
+      UI_GROUP_COLUMNS_MORPHOLOGY,
+      0,
+      GAL_TYPE_INVALID,
+      GAL_OPTIONS_RANGE_ANY,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET,
+      ui_column_codes_ll
+    },
+    {
+      "fracsumarea1",
+      UI_KEY_FRACSUMAREA1,
+      0,
+      0,
+      "Area of 1st --fracsum of --brightness.",
+      UI_GROUP_COLUMNS_MORPHOLOGY,
+      0,
+      GAL_TYPE_INVALID,
+      GAL_OPTIONS_RANGE_ANY,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET,
+      ui_column_codes_ll
+    },
+    {
+      "fracsumarea2",
+      UI_KEY_FRACSUMAREA2,
+      0,
+      0,
+      "Area of 2nd --fracsum of --brightness.",
+      UI_GROUP_COLUMNS_MORPHOLOGY,
+      0,
+      GAL_TYPE_INVALID,
+      GAL_OPTIONS_RANGE_ANY,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET,
+      ui_column_codes_ll
+    },
+    {
+      "fracradiusobs1",
+      UI_KEY_FRACRADIUSOBS1,
+      0,
+      0,
+      "Radius calculated from --fracsumarea1.",
+      UI_GROUP_COLUMNS_MORPHOLOGY,
+      0,
+      GAL_TYPE_INVALID,
+      GAL_OPTIONS_RANGE_ANY,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET,
+      ui_column_codes_ll
+    },
+    {
+      "fracradiusobs2",
+      UI_KEY_FRACRADIUSOBS2,
+      0,
+      0,
+      "Radius calculated from --fracsumarea2.",
+      UI_GROUP_COLUMNS_MORPHOLOGY,
+      0,
+      GAL_TYPE_INVALID,
+      GAL_OPTIONS_RANGE_ANY,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET,
+      ui_column_codes_ll
+    },
+
+
+
 
 
     {0}
diff --git a/bin/mkcatalog/columns.c b/bin/mkcatalog/columns.c
index a201b57..c9ba603 100644
--- a/bin/mkcatalog/columns.c
+++ b/bin/mkcatalog/columns.c
@@ -318,6 +318,9 @@ columns_sanity_check(struct mkcatalogparams *p)
           case UI_KEY_GEOSEMIMAJOR:
           case UI_KEY_GEOSEMIMINOR:
           case UI_KEY_GEOAXISRATIO:
+          case UI_KEY_HALFRADIUSOBS:
+          case UI_KEY_FRACRADIUSOBS1:
+          case UI_KEY_FRACRADIUSOBS2:
           case UI_KEY_GEOPOSITIONANGLE:
             error(EXIT_FAILURE, 0, "columns requiring second moment "
                   "calculations (like semi-major, semi-minor, axis ratio "
@@ -450,21 +453,6 @@ 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";
@@ -1639,6 +1627,99 @@ columns_define_alloc(struct mkcatalogparams *p)
           oiflag[ OCOL_GXY    ] = ciflag[ CCOL_GXY    ] = 1;
           break;
 
+        case UI_KEY_HALFSUMAREA:
+          name           = "HALF_SUM_AREA";
+          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_HALFSUMSB:
+          name           = "HALF_SUM_SB";
+          unit           = MKCATALOG_NO_UNIT;
+          ocomment       = "Surface brightness within area containing half 
total sum.";
+          otype          = GAL_TYPE_FLOAT32;
+          ctype          = GAL_TYPE_FLOAT32;
+          disp_fmt       = GAL_TABLE_DISPLAY_FMT_GENERAL;
+          disp_width     = 10;
+          disp_precision = 5;
+          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_FRACSUMAREA1:
+        case UI_KEY_FRACSUMAREA2:
+          name           = ( colcode->v==UI_KEY_FRACSUMAREA1
+                             ? "FRAC_SUM_AREA_1"
+                             : "FRAC_SUM_AREA_2" );
+          unit           = "counter";
+          ocomment       = "Number of brightest pixels containing given 
fraction 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;
+          if(colcode->v==UI_KEY_FRACSUMAREA1)
+            oiflag[ OCOL_NUMFRACSUM1 ] = ciflag[ CCOL_NUMFRACSUM1 ] = 1;
+          else
+            oiflag[ OCOL_NUMFRACSUM2 ] = ciflag[ CCOL_NUMFRACSUM2 ] = 1;
+          break;
+
+        case UI_KEY_HALFRADIUSOBS:
+        case UI_KEY_FRACRADIUSOBS1:
+        case UI_KEY_FRACRADIUSOBS2:
+          unit           = "pixels";
+          otype          = GAL_TYPE_FLOAT32;
+          ctype          = GAL_TYPE_FLOAT32;
+          disp_fmt       = GAL_TABLE_DISPLAY_FMT_GENERAL;
+          disp_width     = 10;
+          disp_precision = 3;
+          oiflag[ OCOL_NUM        ] = ciflag[ CCOL_NUM        ] = 1; /* 
halfsumarea */
+          oiflag[ OCOL_SUM        ] = ciflag[ CCOL_SUM        ] = 1;
+          oiflag[ OCOL_SUMWHT     ] = ciflag[ CCOL_SUMWHT     ] = 1; /* 
axisratio. */
+          oiflag[ OCOL_VX         ] = ciflag[ CCOL_VX         ] = 1;
+          oiflag[ OCOL_VY         ] = ciflag[ CCOL_VY         ] = 1;
+          oiflag[ OCOL_VXX        ] = ciflag[ CCOL_VXX        ] = 1;
+          oiflag[ OCOL_VYY        ] = ciflag[ CCOL_VYY        ] = 1;
+          oiflag[ OCOL_VXY        ] = ciflag[ CCOL_VXY        ] = 1;
+          oiflag[ OCOL_NUMALL     ] = ciflag[ CCOL_NUMALL     ] = 1;
+          oiflag[ OCOL_GX         ] = ciflag[ CCOL_GX         ] = 1;
+          oiflag[ OCOL_GY         ] = ciflag[ CCOL_GY         ] = 1;
+          oiflag[ OCOL_GXX        ] = ciflag[ CCOL_GXX        ] = 1;
+          oiflag[ OCOL_GYY        ] = ciflag[ CCOL_GYY        ] = 1;
+          oiflag[ OCOL_GXY        ] = ciflag[ CCOL_GXY        ] = 1;
+          switch(colcode->v)
+            {
+            case UI_KEY_HALFRADIUSOBS:
+              name="HALF_RADIUS_OBS";
+              oiflag[ OCOL_NUMHALFSUM  ] = ciflag[ CCOL_NUMHALFSUM  ] = 1;
+              ocomment = "Radius derived from area of half of total sum.";
+              break;
+            case UI_KEY_FRACRADIUSOBS1:
+              name="FRAC_RADIUS_OBS_1";
+              oiflag[ OCOL_NUMFRACSUM1 ] = ciflag[ CCOL_NUMFRACSUM1 ] = 1;
+              ocomment = "Radius derived from area of 1st fraction of total 
sum.";
+              break;
+            case UI_KEY_FRACRADIUSOBS2:
+              name="FRAC_RADIUS_OBS_2";
+              oiflag[ OCOL_NUMFRACSUM2 ] = ciflag[ CCOL_NUMFRACSUM2 ] = 1;
+              ocomment = "Radius derived from area of 2nd fraction of total 
sum.";
+              break;
+            }
+          break;
+
         default:
           error(EXIT_FAILURE, 0, "%s: a bug! please contact us at %s to fix "
                 "the problem. The code %d is not an internally recognized "
@@ -1646,7 +1727,7 @@ columns_define_alloc(struct mkcatalogparams *p)
         }
 
 
-      /* If this is an objects column, add it to the list of columns. We
+      /* If this is an object's column, add it to the list of columns. We
          will be using the 'status' element to keep the MakeCatalog code
          for the columns. */
       if(otype!=GAL_TYPE_INVALID)
@@ -2003,6 +2084,7 @@ 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 coord[3]={GAL_BLANK_SIZE_T, GAL_BLANK_SIZE_T, GAL_BLANK_SIZE_T};
@@ -2054,10 +2136,6 @@ 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;
@@ -2323,6 +2401,38 @@ columns_fill(struct mkcatalog_passparams *pp)
           ((float *)colarr)[oind] = columns_second_order(pp, oi, key, 0);
           break;
 
+        case UI_KEY_HALFSUMAREA:
+          ((int32_t *)colarr)[oind] = oi[OCOL_NUMHALFSUM];
+          break;
+
+        case UI_KEY_FRACSUMAREA1:
+          ((int32_t *)colarr)[oind] = oi[OCOL_NUMFRACSUM1];
+          break;
+
+        case UI_KEY_FRACSUMAREA2:
+          ((int32_t *)colarr)[oind] = oi[OCOL_NUMFRACSUM2];
+          break;
+
+        case UI_KEY_HALFSUMSB:
+          ((float *)colarr)[oind] = oi[OCOL_SUM]/2.0f/oi[OCOL_NUMHALFSUM];
+          break;
+
+        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. */
+          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 ) );
+          tmp = sqrt( oi[tmpind]/(tmp*M_PI) );
+          ((float *)colarr)[oind] = tmp<1e-6 ? NAN : tmp;
+          break;
+
         default:
           error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
                 "solve the problem. the output column code %d not recognized "
@@ -2365,10 +2475,6 @@ 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;
@@ -2581,6 +2687,36 @@ columns_fill(struct mkcatalog_passparams *pp)
             ((float *)colarr)[cind] = columns_second_order(pp, ci, key, 1);
             break;
 
+          case UI_KEY_HALFSUMAREA:
+            ((int32_t *)colarr)[cind] = ci[CCOL_NUMHALFSUM];
+            break;
+
+          case UI_KEY_FRACSUMAREA1:
+            ((int32_t *)colarr)[cind] = ci[CCOL_NUMFRACSUM1];
+            break;
+
+          case UI_KEY_FRACSUMAREA2:
+            ((int32_t *)colarr)[cind] = ci[CCOL_NUMFRACSUM2];
+            break;
+
+          case UI_KEY_HALFSUMSB:
+            ((float *)colarr)[cind] = ci[CCOL_SUM]/2.0f/ci[CCOL_NUMHALFSUM];
+            break;
+
+          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 ) );
+            tmp = sqrt( ci[tmpind]/(tmp*M_PI) );
+            ((float *)colarr)[cind] = tmp<1e-6 ? NAN : tmp;
+            break;
+
           default:
             error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
                   "solve the problem. The output column code %d not "
diff --git a/bin/mkcatalog/main.h b/bin/mkcatalog/main.h
index 9298a45..d8d1c25 100644
--- a/bin/mkcatalog/main.h
+++ b/bin/mkcatalog/main.h
@@ -74,7 +74,6 @@ enum objectcols
     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.                */
@@ -110,6 +109,9 @@ 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_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. */
     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.              */
@@ -131,7 +133,6 @@ 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.                */
@@ -176,6 +177,11 @@ 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_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. */
+    CCOL_FRACSUMAREA1,   /* Area/Number containing frac of total sum. */
+    CCOL_FRACSUMAREA2,   /* Area/Number containing frac of total sum. */
 
     CCOL_NUMCOLS,        /* SHOULD BE LAST: total number of columns.  */
   };
@@ -223,6 +229,8 @@ struct mkcatalogparams
   float              upnsigma;  /* Multiple of sigma to define up-lim.  */
   int32_t       checkuplim[2];  /* Object & clump ID to check dist.     */
 
+  gal_data_t         *fracsum;  /* Fractions to use in --fracsumarea.   */
+
   /* Internal. */
   char           *relabclumps;  /* Name of new file for clump labels.   */
   time_t              rawtime;  /* Starting time of the program.        */
diff --git a/bin/mkcatalog/mkcatalog.c b/bin/mkcatalog/mkcatalog.c
index 3204f2c..d05f80d 100644
--- a/bin/mkcatalog/mkcatalog.c
+++ b/bin/mkcatalog/mkcatalog.c
@@ -186,6 +186,8 @@ mkcatalog_single_object(void *in_prm)
           || p->oiflag[ OCOL_SIGCLIPNUM ]
           || p->oiflag[ OCOL_SIGCLIPSTD ]
           || p->oiflag[ OCOL_SIGCLIPMEAN ]
+          || p->oiflag[ OCOL_NUMFRACSUM1 ]
+          || p->oiflag[ OCOL_NUMFRACSUM2 ]
           || p->oiflag[ OCOL_SIGCLIPMEDIAN ])
         parse_order_based(&pp);
 
diff --git a/bin/mkcatalog/parse.c b/bin/mkcatalog/parse.c
index ca2b78f..63a9964 100644
--- a/bin/mkcatalog/parse.c
+++ b/bin/mkcatalog/parse.c
@@ -1182,12 +1182,36 @@ parse_clumps(struct mkcatalog_passparams *pp)
 
 
 
-size_t
-parse_area_of_half_sum(gal_data_t *values, double sum)
+static double
+parse_frac_sum_find(gal_data_t *sorted_d, double sum, double frac)
 {
   size_t i;
+  double sumcheck=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;
+
+  /* Return the final value. Note that if the index is zero, we should
+     actually return 1, because we are starting with the maximum. */
+  return i==0 ? 1 : i;
+}
+
+
+
+
+
+static void
+parse_area_of_frac_sum(struct mkcatalog_passparams *pp, gal_data_t *values,
+                       double *outarr, int o1c0)
+{
+  struct mkcatalogparams *p=pp->p;
+
   gal_data_t *sorted_d;
-  double sumcheck=0.0f, *sorted;
+  uint8_t *flag = o1c0 ? p->oiflag : p->ciflag;
+  double sum = o1c0 ? outarr[OCOL_SUM] : outarr[CCOL_SUM];
+  double *fracsum = p->fracsum ? p->fracsum->array : NULL;
 
   /* Allocate the array to use. */
   sorted_d = ( values->type==GAL_TYPE_FLOAT64
@@ -1197,14 +1221,28 @@ parse_area_of_half_sum(gal_data_t *values, double sum)
   /* 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;
+
+  /* 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);
+
+  if(flag[ o1c0 ? OCOL_NUMFRACSUM1 : CCOL_NUMFRACSUM1 ])
+    outarr[ o1c0 ? OCOL_NUMFRACSUM1 : CCOL_NUMFRACSUM1 ]
+      = parse_frac_sum_find(sorted_d, sum, fracsum[0]);
+
+  if(flag[ o1c0 ? OCOL_NUMFRACSUM2 : CCOL_NUMFRACSUM2 ])
+    outarr[ o1c0 ? OCOL_NUMFRACSUM2 : CCOL_NUMFRACSUM2 ]
+      = parse_frac_sum_find(sorted_d, sum, fracsum[1]);
+
+  /* For a check.
+  printf("%g, %g, %g\n", outarr[OCOL_NUMFRACSUM1], outarr[OCOL_NUMHALFSUM],
+         outarr[OCOL_NUMFRACSUM2]);
+  exit(1);
+  */
 
   /* Clean up and return. */
   if(sorted_d!=values) gal_data_free(sorted_d);
-  return i;
 }
 
 
@@ -1352,9 +1390,12 @@ 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]);
 
+  /* Fractional values. */
+  if( p->oiflag[    OCOL_NUMHALFSUM   ]
+      || p->oiflag[ OCOL_NUMFRACSUM1 ]
+      || p->oiflag[ OCOL_NUMFRACSUM2 ] )
+    parse_area_of_frac_sum(pp, objvals, pp->oi, 1);
 
   /* Clean up the object values. */
   gal_data_free(objvals);
@@ -1404,8 +1445,10 @@ parse_order_based(struct mkcatalog_passparams *pp)
             }
 
           /* Estimate half of the total sum. */
-          if(p->ciflag[ CCOL_NUMHALFSUM ])
-            ci[CCOL_NUMHALFSUM]=parse_area_of_half_sum(clumpsvals[i], 
ci[CCOL_SUM]);
+          if( p->ciflag[ CCOL_NUMHALFSUM ]
+              || p->ciflag[ CCOL_NUMFRACSUM1 ]
+              || p->ciflag[ CCOL_NUMFRACSUM2 ] )
+            parse_area_of_frac_sum(pp, clumpsvals[i], ci, 0);
 
           /* Clean up this clump's values. */
           gal_data_free(clumpsvals[i]);
diff --git a/bin/mkcatalog/ui.c b/bin/mkcatalog/ui.c
index ae5dd34..501d043 100644
--- a/bin/mkcatalog/ui.c
+++ b/bin/mkcatalog/ui.c
@@ -361,8 +361,10 @@ static void
 ui_read_check_only_options(struct mkcatalogparams *p)
 {
   float tmp;
-  size_t one=1;
+  double *darr;
   char *tailptr;
+  size_t i, one=1;
+  gal_list_i32_t *colcode;
 
   /* If an upper-limit check table is requested with a specific clump, but
      no clump catalog has been requested, then abort and inform the
@@ -407,6 +409,38 @@ ui_read_check_only_options(struct mkcatalogparams *p)
           *((float *)(p->std->array))=tmp;
         }
     }
+
+  /* Make sure that '--fracsum' is given if necessary and that the fracsum
+     values are less than one. */
+  for(colcode=p->columnids; colcode!=NULL; colcode=colcode->next)
+    if( (colcode->v==UI_KEY_FRACSUMAREA1 || colcode->v==UI_KEY_FRACSUMAREA2)
+        && p->fracsum==NULL )
+      error(EXIT_FAILURE, 0, "please specify your required fractions for "
+            "'--fracsumarea' using the '--fracsum' option (for example "
+            "'--fracsum=0.25,0.75' for two columns)");
+  if(p->fracsum)
+    {
+      /* Currently only 2 fracsums are accepted. */
+      if(p->fracsum->size>2)
+        error(EXIT_FAILURE, 0, "%zu values given to '--fracsum', only two "
+              "values are currently supported", p->fracsum->size);
+
+      /* Check the values. */
+      darr=p->fracsum->array;
+      for(i=0;i<p->fracsum->size;++i)
+        if(darr[i]<=0.0f || darr[i]>=1.0f)
+          error(EXIT_FAILURE, 0, "%g is not acceptable for '--fracsum', "
+                "values should be larger than 0 and less than 1", darr[i]);
+
+      /* If a second fracum column is necessary, make sure two values are
+         given to --fracsum. */
+      for(colcode=p->columnids; colcode!=NULL; colcode=colcode->next)
+        if(colcode->v==UI_KEY_FRACSUMAREA2 && p->fracsum->size==1)
+          error(EXIT_FAILURE, 0, "only one value given to '--fracsum', "
+                "but '--fracsumarea2' was requested! You need to give "
+                "the requested fraction as a second value to '--fracsum', "
+                "separated by a comma (,)");
+    }
 }
 
 
@@ -929,7 +963,6 @@ 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;
@@ -963,6 +996,9 @@ 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_NUMHALFSUM:         *values        = 1;          break;
+        case OCOL_NUMFRACSUM1:        *values        = 1;          break;
+        case OCOL_NUMFRACSUM2:        *values        = 1;          break;
         case OCOL_C_NUMALL:           /* Only clump labels.  */    break;
         case OCOL_C_NUM:              *values        = 1;          break;
         case OCOL_C_SUM:              *values        = 1;          break;
@@ -991,7 +1027,6 @@ 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;
@@ -1034,6 +1069,9 @@ 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_NUMHALFSUM:       *values        = 1;          break;
+          case CCOL_NUMFRACSUM1:      *values        = 1;          break;
+          case CCOL_NUMFRACSUM2:      *values        = 1;          break;
           default:
             error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
                   "fix the problem. The code %zu is not a recognized "
diff --git a/bin/mkcatalog/ui.h b/bin/mkcatalog/ui.h
index 155801f..33b17e3 100644
--- a/bin/mkcatalog/ui.h
+++ b/bin/mkcatalog/ui.h
@@ -34,6 +34,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 enum program_args_groups
 {
   UI_GROUP_UPPERLIMIT = GAL_OPTIONS_GROUP_AFTER_COMMON,
+  UI_GROUP_OTHERSETTINGS,
   UI_GROUP_COLUMNS_IDS,
   UI_GROUP_COLUMNS_POSITION_PIXEL,
   UI_GROUP_COLUMNS_POSITION_WCS,
@@ -103,11 +104,11 @@ enum option_keys_enum
   UI_KEY_UPNSIGMA,
   UI_KEY_CHECKUPLIM,
   UI_KEY_NOCLUMPSORT,
+  UI_KEY_FRACSUM,
 
   UI_KEY_OBJID,                         /* Catalog columns. */
   UI_KEY_IDINHOSTOBJ,
   UI_KEY_AREAXY,
-  UI_KEY_AREAHALFSUM,
   UI_KEY_CLUMPSAREA,
   UI_KEY_WEIGHTAREA,
   UI_KEY_GEOAREA,
@@ -168,6 +169,13 @@ enum option_keys_enum
   UI_KEY_GEOSEMIMINOR,
   UI_KEY_GEOAXISRATIO,
   UI_KEY_GEOPOSITIONANGLE,
+  UI_KEY_HALFSUMAREA,
+  UI_KEY_HALFSUMSB,
+  UI_KEY_HALFRADIUSOBS,
+  UI_KEY_FRACSUMAREA1,
+  UI_KEY_FRACSUMAREA2,
+  UI_KEY_FRACRADIUSOBS1,
+  UI_KEY_FRACRADIUSOBS2,
 };
 
 
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 4485070..e142b88 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -15849,6 +15849,11 @@ This option takes two values: the first is the 
multiple of @mymath{\sigma}, and
 If the latter is larger than 1, it is read as an integer number and will be 
the number of times to clip.
 If it is smaller than 1, it is interpreted as the tolerance level to stop 
clipping. See @ref{Sigma clipping} for a complete explanation.
 
+@item --fracsum=FLT[,FLT]
+The fractions of summed flux in objects or clumps to be used in the 
@option{--fracsumarea1} and @option{--fracsumarea2} columns, see 
@ref{MakeCatalog measurements}.
+The values must be larger than 0 and smaller than 1.
+When only @option{--fracsumarea1} is requested, one value must be given to 
this option, but if @option{--fracsumarea2} is also requested, two values must 
be given here.
+The values can be written as simple floating point numbers, or as fractions, 
for example @code{0.25,0.75} and @code{0.25,3/4} are the same.
 @end table
 
 
@@ -16265,13 +16270,49 @@ 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
+@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.
 
 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 --halfsumsb
+Surface brightness (dataset units per pixel) within the area that contains 
half the total sum of the label's pixels (object or clump).
+Effectively, this column is just half the value of the @option{--brightness} 
column, divided by the value of the @option{--halfsumarea} column.
+
+Therefore it is very important to understand the systematics of this column 
and potential biases, please see the description in @option{--halfsumarea}.
+
+@item --halfradiusobs
+Radius (in units of pixels) derived from the area that contains half the total 
sum of the label's pixels (value reported by @option{--halfsumarea}).
+If the area is @mymath{A_h} and the axis ratio is @mymath{q}, then the value 
returned in this column is @mymath{\sqrt{A_h/({\pi}q)}}.
+This option is a good measure of the concentration of the @emph{observed} 
(after PSF convolution and noisy) object or clump,
+But as described below it underestimates the effective radius.
+Also, it should be used in causion with objects, it reliable with clumps, see 
the note under @option{--halfradiusarea}.
+
+@cindex Ellipse area
+@cindex Area, ellipse
+Recall that in general, for an ellipse with semi-major axis @mymath{a}, 
semi-minor axis @mymath{b}, and axis ratio @mymath{q=b/a} the area (@mymath{A}) 
is @mymath{A={\pi}ab={\pi}qa^2}.
+For a circle (where @mymath{q=1}), this simplifies to the familiar 
@mymath{A={\pi}a^2}.
+
+@cindex S@'ersic profile
+@cindex Effective radius
+The suffix @code{obs} is added to the option name to avoid confusing this with 
a profile's effective radius for S@'ersic profiles, commonly written as 
@mymath{r_e}.
+For more on @mymath{r_e}, please see @ref{Galaxies}.
+Therefore, when @mymath{r_e} is meaningful for the target (the target is 
elliptically symmetric and can be parametrized as a S@'ersic profile), 
@mymath{r_e} should be derived from fitting the profile with a S@'ersic 
function which has been convolved with the PSF.
+But from the equation above, you see that this radius is derived from the raw 
image's labeled values (after convolution, with no parametric profile), so this 
column's value will generally be (much) smaller than @mymath{r_e}, depending on 
the PSF, depth of the dataset, the morphology, or if a fraction of the profile 
falls on the edge of the image.
+
+@item --fracsumarea1
+@itemx --fracsumarea2
+Similar to @option{--halfsumarea}, but for user-defined fractions of total sum.
+The fractions are given through the @option{--fracsum} option (that can take 
two values) and is described in @ref{MakeCatalog inputs and basic settings}.
+Recall that in @option{--halfsumarea}, the fraction is fixed to 0.5 so added 
with these two columns, you can sample three parts of the profile.
+
+@item --fracradiusobs1
+@itemx --fracradiusobs2
+Similar to @option{--halfradiusobs}, but for user-defined fractions of total 
sum.
+The fractions are given through the @option{--fracsum} option (that can take 
two values) and is described in @ref{MakeCatalog inputs and basic settings}.
+Recall that in @option{--halfradiusobs}, the fraction is fixed to 0.5 so added 
with these two columns, you can sample three parts of the profile.
 
 @item --clumpsarea
 [Objects] The total area of all the clumps in this object.
@@ -16970,7 +17011,7 @@ Today, most practitioners agree that the flux of 
galaxies can be modeled with on
 @cindex Radius, effective
 @cindex de Vaucouleur profile
 @cindex G@'erard de Vaucouleurs
-G@'erard de Vaucouleurs (1918-1995) was first to show in 1948 that this 
function best fits the galaxy light profiles, with the only difference that he 
held @mymath{n} fixed to a value of 4.
+G@'erard de Vaucouleurs (1918-1995) was first to show in 1948 that this 
function resembles the galaxy light profiles, with the only difference that he 
held @mymath{n} fixed to a value of 4.
 20 years later in 1968, J. L. S@'ersic showed that @mymath{n} can have a 
variety of values and does not necessarily need to be 4.
 This profile depends on the effective radius (@mymath{r_e}) which is defined 
as the radius which contains half of the profile brightness (see @ref{Profile 
magnitude}).
 @mymath{I_e} is the flux at the effective radius.



reply via email to

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