[gnuastro-commits] master 3882956 077/113: Narrow-band areas from MakeCatalog in 3D datasets

[Mohammad Akhlaghi]

branch: master
commit 38829569da1978e7e6bfaa19ec1306d2cb6f82d8
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    Narrow-band areas from MakeCatalog in 3D datasets
    
    When working with 3D datasets (especially IFU datasets), it is commonly
    necessary to get narrow-band information for the objects/clumps. With this
    commit the steps have been done for doing this.
    
    Currently it is only for the area, but later, we will be able to do
    elliptical measurements and etc also.
---
 bin/mkcatalog/args.h    |  28 ++++++++++++
 bin/mkcatalog/columns.c |  44 +++++++++++++++++++
 bin/mkcatalog/main.h    |   6 ++-
 bin/mkcatalog/parse.c   | 111 ++++++++++++++++++++++++++++++++++++++++++++----
 bin/mkcatalog/ui.c      |   4 ++
 bin/mkcatalog/ui.h      |   2 +
 doc/gnuastro.texi       |  16 ++++++-
 7 files changed, 200 insertions(+), 11 deletions(-)

diff --git a/bin/mkcatalog/args.h b/bin/mkcatalog/args.h
index c06b3d9..cfdfb5e 100644
--- a/bin/mkcatalog/args.h
+++ b/bin/mkcatalog/args.h
@@ -1131,6 +1131,20 @@ struct argp_option program_options[] =
       ui_column_codes_ll
     },
     {
+      "areaxy",
+      UI_KEY_AREAXY,
+      0,
+      0,
+      "Projected area in first two dimentions.",
+      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,
@@ -1173,6 +1187,20 @@ struct argp_option program_options[] =
       ui_column_codes_ll
     },
     {
+      "geoareaxy",
+      UI_KEY_GEOAREAXY,
+      0,
+      0,
+      "Projected geoarea in first two dimensions.",
+      UI_GROUP_COLUMNS_MORPHOLOGY,
+      0,
+      GAL_TYPE_INVALID,
+      GAL_OPTIONS_RANGE_ANY,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET,
+      ui_column_codes_ll
+    },
+    {
       "semimajor",
       UI_KEY_SEMIMAJOR,
       0,
diff --git a/bin/mkcatalog/columns.c b/bin/mkcatalog/columns.c
index 56f7881..49ee9d3 100644
--- a/bin/mkcatalog/columns.c
+++ b/bin/mkcatalog/columns.c
@@ -272,6 +272,8 @@ columns_sanity_check(struct mkcatalogparams *p)
       for(colcode=p->columnids; colcode!=NULL; colcode=colcode->next)
         switch(colcode->v)
           {
+          case UI_KEY_AREAXY:
+          case UI_KEY_GEOAREAXY:
           case UI_KEY_Z:
           case UI_KEY_GEOZ:
           case UI_KEY_CLUMPSZ:
@@ -416,6 +418,19 @@ columns_define_alloc(struct mkcatalogparams *p)
           oiflag[ OCOL_NUM ] = ciflag[ CCOL_NUM ] = 1;
           break;
 
+        case UI_KEY_AREAXY:
+          name           = "AREAXY";
+          unit           = "counter";
+          ocomment       = "Projected valued pixels in first two dimensions.";
+          ccomment       = ocomment;
+          otype          = GAL_TYPE_INT32;
+          ctype          = GAL_TYPE_INT32;
+          disp_fmt       = 0;
+          disp_width     = 6;
+          disp_precision = 0;
+          oiflag[ OCOL_NUMXY ] = ciflag[ CCOL_NUMXY ] = 1;
+          break;
+
         case UI_KEY_CLUMPSAREA:
           name           = "AREA_CLUMPS";
           unit           = "counter";
@@ -455,6 +470,19 @@ columns_define_alloc(struct mkcatalogparams *p)
           oiflag[ OCOL_NUMALL ] = ciflag[ CCOL_NUMALL ] = 1;
           break;
 
+        case UI_KEY_GEOAREAXY:
+          name           = "AREA_FULL_XY";
+          unit           = "counter";
+          ocomment       = "Project number in first two dimensions.";
+          ccomment       = ocomment;
+          otype          = GAL_TYPE_INT32;
+          ctype          = GAL_TYPE_INT32;
+          disp_fmt       = 0;
+          disp_width     = 6;
+          disp_precision = 0;
+          oiflag[ OCOL_NUMALLXY ] = ciflag[ CCOL_NUMALLXY ] = 1;
+          break;
+
         case UI_KEY_X:
           name           = "X";
           unit           = "pixel";
@@ -1699,6 +1727,10 @@ columns_fill(struct mkcatalog_passparams *pp)
           ((int32_t *)colarr)[oind] = oi[OCOL_NUM];
           break;
 
+        case UI_KEY_AREAXY:
+          ((int32_t *)colarr)[oind] = oi[OCOL_NUMXY];
+          break;
+
         case UI_KEY_CLUMPSAREA:
           ((int32_t *)colarr)[oind] = oi[OCOL_C_NUM];
           break;
@@ -1711,6 +1743,10 @@ columns_fill(struct mkcatalog_passparams *pp)
           ((int32_t *)colarr)[oind] = oi[OCOL_NUMALL];
           break;
 
+        case UI_KEY_GEOAREAXY:
+          ((int32_t *)colarr)[oind] = oi[OCOL_NUMALLXY];
+          break;
+
         case UI_KEY_X:
           ((float *)colarr)[oind] = POS_V_G(oi, OCOL_SUMWHT, OCOL_NUMWHT,
                                             OCOL_VX, OCOL_GX);
@@ -1966,6 +2002,10 @@ columns_fill(struct mkcatalog_passparams *pp)
             ((int32_t *)colarr)[cind]=ci[CCOL_NUM];
             break;
 
+          case UI_KEY_AREAXY:
+            ((int32_t *)colarr)[cind]=ci[CCOL_NUMXY];
+            break;
+
           case UI_KEY_WEIGHTAREA:
             ((int32_t *)colarr)[cind]=ci[CCOL_NUMWHT];
             break;
@@ -1974,6 +2014,10 @@ columns_fill(struct mkcatalog_passparams *pp)
             ((int32_t *)colarr)[cind]=ci[CCOL_NUMALL];
             break;
 
+          case UI_KEY_GEOAREAXY:
+            ((int32_t *)colarr)[cind]=ci[CCOL_NUMALLXY];
+            break;
+
           case UI_KEY_X:
             ((float *)colarr)[cind] = POS_V_G(ci, CCOL_SUMWHT, CCOL_NUMALL,
                                               CCOL_VX, CCOL_GX);
diff --git a/bin/mkcatalog/main.h b/bin/mkcatalog/main.h
index dbf54e3..4519191 100644
--- a/bin/mkcatalog/main.h
+++ b/bin/mkcatalog/main.h
@@ -71,7 +71,9 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 enum objectcols
   {
     OCOL_NUMALL,         /* Number of all pixels with this label.     */
-    OCOL_NUM,            /* Number of pixels with a value.            */
+    OCOL_NUMALLXY,       /* Number of all pixels in first two dims.   */
+    OCOL_NUM,            /* Number of values used in this object.     */
+    OCOL_NUMXY,          /* Number of values in the first two dims.   */
     OCOL_SUM,            /* Sum of (value-sky) in object.             */
     OCOL_SUM_VAR,        /* Varience of sum (for brightness error).   */
     OCOL_MEDIAN,         /* Median of value in object.                */
@@ -113,7 +115,9 @@ enum objectcols
 enum clumpcols
   {
     CCOL_NUMALL,         /* Number of pixels in clump.                */
+    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_SUM,            /* River subtracted brightness.              */
     CCOL_SUM_VAR,        /* Variance of sum (for brightness error).   */
     CCOL_MEDIAN,         /* Median of values in clump.                */
diff --git a/bin/mkcatalog/parse.c b/bin/mkcatalog/parse.c
index c3bdd0d..95ccfc0 100644
--- a/bin/mkcatalog/parse.c
+++ b/bin/mkcatalog/parse.c
@@ -114,8 +114,10 @@ parse_objects(struct mkcatalog_passparams *pp)
   size_t ndim=p->objects->ndim, *dsize=p->objects->dsize;
 
   double *oi=pp->oi;
-  size_t d, increment=0, num_increment=1;
+  gal_data_t *xybin=NULL;
+  uint8_t *xybinarr=NULL, *u, *uf;
   float st, sval, *V=NULL, *SK=NULL, *ST=NULL;
+  size_t d, pind=0, increment=0, num_increment=1;
   int32_t *O, *OO, *C=NULL, *objects=p->objects->array;
   float *std=p->std?p->std->array:NULL, *sky=p->sky?p->sky->array:NULL;
 
@@ -151,6 +153,17 @@ parse_objects(struct mkcatalog_passparams *pp)
                : NULL );
 
 
+  /* If an XY projection area is requested, we'll need to allocate an array
+     to keep the projected space.*/
+  if( oif[    OCOL_NUMALLXY ]
+      || oif[ OCOL_NUMXY    ] )
+    {
+      xybin=gal_data_alloc(NULL, GAL_TYPE_UINT8, 2, &pp->tile->dsize[1], NULL,
+                           1, p->cp.minmapsize, NULL, NULL, NULL);
+      xybinarr=xybin->array;
+    }
+
+
   /* Parse each contiguous patch of memory covered by this object. */
   while( pp->start_end_inc[0] + increment <= pp->start_end_inc[1] )
     {
@@ -176,7 +189,8 @@ parse_objects(struct mkcatalog_passparams *pp)
 
 
               /* Add to the area of this object. */
-              if(oif[ OCOL_NUMALL ]) oi[ OCOL_NUMALL ]++;
+              if(oif[ OCOL_NUMALL   ]) oi[ OCOL_NUMALL ]++;
+              if(oif[ OCOL_NUMALLXY ]) xybinarr[ pind ]=1;
 
 
               /* Geometric coordinate measurements. */
@@ -226,6 +240,7 @@ parse_objects(struct mkcatalog_passparams *pp)
                   /* General flux summations. */
                   if(oif[ OCOL_NUM    ]) oi[ OCOL_NUM     ]++;
                   if(oif[ OCOL_SUM    ]) oi[ OCOL_SUM     ] += *V;
+                  if(oif[ OCOL_NUMXY  ]) xybinarr[ pind ]=2;
 
                   /* Get the necessary clump information. */
                   if(p->clumps && *C>0)
@@ -297,6 +312,7 @@ parse_objects(struct mkcatalog_passparams *pp)
             }
 
           /* Increment the other pointers. */
+          if( xybin                ) ++pind;
           if( p->values            ) ++V;
           if( p->clumps            ) ++C;
           if( p->sky && pp->st_sky ) ++SK;
@@ -307,11 +323,39 @@ parse_objects(struct mkcatalog_passparams *pp)
       /* Increment to the next contiguous region of this tile. */
       increment += ( gal_tile_block_increment(p->objects, dsize,
                                               num_increment++, NULL) );
+
+      /* If a 2D projection is requested, see if we should initialize (set
+         to zero) the projection-index (`pind') not. */
+      if(xybin && (num_increment-1)%p->objects->dsize[1]==0 )
+        pind=0;
+    }
+
+  /* Write the projected area columns. */
+  if(xybin)
+    {
+      /* Any non-zero pixel must be set for NUMALLXY. */
+      uf=(u=xybin->array)+xybin->size;
+      do
+        if(*u)
+          {
+            if(oif[ OCOL_NUMALLXY ]          ) oi[ OCOL_NUMALLXY ]++;
+            if(oif[ OCOL_NUMXY    ] && *u==2 ) oi[ OCOL_NUMXY    ]++;
+          }
+      while(++u<uf);
+
+      /* For a check on the projected 2D areas.
+      if(xybin && pp->object==2)
+        {
+          gal_fits_img_write(xybin, "xybin.fits", NULL, NULL);
+          exit(0);
+        }
+      */
     }
 
   /* Clean up. */
-  if(c)  free(c);
-  if(sc) free(sc);
+  if(c)     free(c);
+  if(sc)    free(sc);
+  if(xybin) gal_data_free(xybin);
 }
 
 
@@ -326,11 +370,12 @@ parse_clumps(struct mkcatalog_passparams *pp)
   size_t ndim=p->objects->ndim, *dsize=p->objects->dsize;
 
   double *ci, *cir;
-  uint8_t *cif=p->ciflag;
+  gal_data_t *xybin=NULL;
   int32_t *O, *OO, *C=NULL, nlab;
+  uint8_t *u, *uf, *cif=p->ciflag;
   float st, 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);
+  size_t i, ii, d, pind=0, increment=0, num_increment=1;
   int32_t *objects=p->objects->array, *clumps=p->clumps->array;
   float *std=p->std?p->std->array:NULL, *sky=p->sky?p->sky->array:NULL;
 
@@ -367,6 +412,18 @@ parse_clumps(struct mkcatalog_passparams *pp)
                      : NULL );
   size_t *dinc = ngblabs ? gal_dimension_increment(ndim, dsize) : NULL;
 
+  /* If an XY projection area is requested, we'll need to allocate an array
+     to keep the projected space.*/
+  if( cif[    CCOL_NUMALLXY ]
+      || cif[ CCOL_NUMXY    ] )
+    {
+      xybin=gal_data_array_calloc(pp->clumpsinobj);
+      for(i=0;i<pp->clumpsinobj;++i)
+        gal_data_initialize(&xybin[i], NULL, GAL_TYPE_UINT8, 2,
+                            &pp->tile->dsize[1], NULL, 1, p->cp.minmapsize,
+                            NULL, NULL, NULL);
+    }
+
 
   /* Parse each contiguous patch of memory covered by this object. */
   while( pp->start_end_inc[0] + increment <= pp->start_end_inc[1] )
@@ -394,7 +451,9 @@ parse_clumps(struct mkcatalog_passparams *pp)
                   ci=&pp->ci[ (*C-1) * CCOL_NUMCOLS ];
 
                   /* Add to the area of this object. */
-                  if(cif[ CCOL_NUMALL ]) ci[ CCOL_NUMALL ]++;
+                  if(cif[ CCOL_NUMALL   ]) ci[ CCOL_NUMALL ]++;
+                  if(cif[ CCOL_NUMALLXY ])
+                    ((uint8_t *)(xybin[*C-1].array))[ pind ] = 1;
 
                   /* Raw-position related measurements. */
                   if(c)
@@ -429,8 +488,10 @@ parse_clumps(struct mkcatalog_passparams *pp)
                   if( p->values && !( p->hasblank && isnan(*V) ) )
                     {
                       /* Fill in the necessary information. */
-                      if(cif[ CCOL_NUM ]) ci[ CCOL_NUM ]++;
-                      if(cif[ CCOL_SUM ]) ci[ CCOL_SUM ] += *V;
+                      if(cif[ CCOL_NUM   ]) ci[ CCOL_NUM ]++;
+                      if(cif[ CCOL_SUM   ]) ci[ CCOL_SUM ] += *V;
+                      if(cif[ CCOL_NUMXY ])
+                        ((uint8_t *)(xybin[*C-1].array))[ pind ] = 2;
                       if( *V > 0.0f )
                         {
                           if(cif[ CCOL_NUMWHT ]) ci[ CCOL_NUMWHT ]++;
@@ -543,6 +604,7 @@ parse_clumps(struct mkcatalog_passparams *pp)
 
           /* Increment the other pointers. */
           ++C;
+          if( xybin                ) ++pind;
           if( p->values            ) ++V;
           if( p->sky && pp->st_sky ) ++SK;
           if( p->std && pp->st_std ) ++ST;
@@ -552,13 +614,44 @@ parse_clumps(struct mkcatalog_passparams *pp)
       /* Increment to the next contiguous region of this tile. */
       increment += ( gal_tile_block_increment(p->objects, dsize,
                                               num_increment++, NULL) );
+
+      /* If a 2D projection is requested, see if we should initialize (set
+         to zero) the projection-index (`pind') not. */
+      if(xybin && (num_increment-1)%p->objects->dsize[1]==0 )
+        pind=0;
     }
 
+
+  /* Write the projected area columns. */
+  if(xybin)
+    for(i=0;i<pp->clumpsinobj;++i)
+      {
+        /* Pointer to make things easier. */
+        ci=&pp->ci[ i * CCOL_NUMCOLS ];
+
+        /* Any non-zero pixel must be set for NUMALLXY. */
+        uf=(u=xybin[i].array)+xybin[i].size;
+        do
+          if(*u)
+            {
+              if(cif[ CCOL_NUMALLXY ]          ) ci[ CCOL_NUMALLXY ]++;
+              if(cif[ CCOL_NUMXY    ] && *u==2 ) ci[ CCOL_NUMXY    ]++;
+            }
+        while(++u<uf);
+
+        /* For a check on the projected 2D areas. */
+        if(xybin && pp->object==2)
+          gal_fits_img_write(&xybin[i], "xybin.fits", NULL, NULL);
+
+      }
+
+
   /* Clean up. */
   if(c)       free(c);
   if(sc)      free(sc);
   if(dinc)    free(dinc);
   if(ngblabs) free(ngblabs);
+  if(xybin)   gal_data_array_free(xybin, pp->clumpsinobj, 1);
 }
 
 
diff --git a/bin/mkcatalog/ui.c b/bin/mkcatalog/ui.c
index be48f60..5554ef8 100644
--- a/bin/mkcatalog/ui.c
+++ b/bin/mkcatalog/ui.c
@@ -698,7 +698,9 @@ ui_necessary_inputs(struct mkcatalogparams *p, int *values, 
int *sky,
       switch(i)
         {
         case OCOL_NUMALL:             /* Only object labels. */    break;
+        case OCOL_NUMALLXY:           /* Only object labels. */    break;
         case OCOL_NUM:                *values        = 1;          break;
+        case OCOL_NUMXY:              *values        = 1;          break;
         case OCOL_SUM:                *values        = 1;          break;
         case OCOL_SUM_VAR:            *values = *std = 1;          break;
         case OCOL_MEDIAN:             *values        = 1;          break;
@@ -747,7 +749,9 @@ ui_necessary_inputs(struct mkcatalogparams *p, int *values, 
int *sky,
         switch(i)
           {
           case CCOL_NUMALL:           /* Only clump labels. */     break;
+          case CCOL_NUMALLXY:         /* Only clump labels. */     break;
           case CCOL_NUM:              *values        = 1;          break;
+          case CCOL_NUMXY:            *values        = 1;          break;
           case CCOL_SUM:              *values        = 1;          break;
           case CCOL_SUM_VAR:          *values = *std = 1;          break;
           case CCOL_MEDIAN:           *values        = 1;          break;
diff --git a/bin/mkcatalog/ui.h b/bin/mkcatalog/ui.h
index 16e63a3..6004fc6 100644
--- a/bin/mkcatalog/ui.h
+++ b/bin/mkcatalog/ui.h
@@ -102,9 +102,11 @@ enum option_keys_enum
 
   UI_KEY_OBJID,                         /* Catalog columns. */
   UI_KEY_IDINHOSTOBJ,
+  UI_KEY_AREAXY,
   UI_KEY_CLUMPSAREA,
   UI_KEY_WEIGHTAREA,
   UI_KEY_GEOAREA,
+  UI_KEY_GEOAREAXY,
   UI_KEY_GEOX,
   UI_KEY_GEOY,
   UI_KEY_GEOZ,
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 33e6412..95ac3b4 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -18010,9 +18010,17 @@ input sky standard deviation image pixels that lie 
over this object.
 
 @item -a
 @itemx --area
-The raw area (number of pixels) in any clump or object independent of what
+The area (number of pixels) in any clump or object independent of what
 pixel it lies over (if it is NaN/blank or unused for example).
 
+@item --areaxy
+@cindex IFU: Integral Field Unit
+@cindex Integral Field Unit
+Similar to @option{--area}, when the clump or object is projected onto the
+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 --clumpsarea
 [Objects] The total area of all the clumps in this object.
 
@@ -18026,6 +18034,12 @@ unlike @option{--area}, pixel values are completely 
ignored in this
 column. For example, if a pixel value is blank, it won't be counted in
 @option{--area}, but will be counted here.
 
+@item --geoareaxy
+Similar to @option{--geoarea}, when the clump or object is projected onto
+the 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 -A
 @itemx --semimajor
 The pixel-value weighted semi-major axis of the profile (assuming it is an