[gnuastro-commits] master 0ebad2b 1/2: MakeCatalog: Rows with no label are not shown in output

[Mohammad Akhlaghi]

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

    MakeCatalog: Rows with no label are not shown in output
    
    Until now, MakeCatalog would keep a row for every integer from 1 to the
    largest label in the input dataset. For outputs of programs like
    NoiseChisel and Segment, this was fine (the labels are consecutive and
    contiguous), but when the labeled image is created by some other scenario,
    This was not good. In some situations, the labels could start from
    arbitrary large numbers and as a result, there would be many blank rows in
    a very large output table.
    
    With this commit, all the rows/IDs that don't contain any pixels are
    removed from the output catalog, making it much more directly useful.
---
 NEWS                      |  9 +++++++++
 bin/mkcatalog/args.h      | 13 +++++++++++++
 bin/mkcatalog/main.h      |  2 ++
 bin/mkcatalog/mkcatalog.c | 46 ++++++++++++++++++++++++++++++++++++++++++++++
 bin/mkcatalog/ui.c        | 32 +++++++++++++++++++++++++++++---
 bin/mkcatalog/ui.h        |  1 +
 doc/gnuastro.texi         |  7 +++++++
 7 files changed, 107 insertions(+), 3 deletions(-)

diff --git a/NEWS b/NEWS
index c000c11..7375d0a 100644
--- a/NEWS
+++ b/NEWS
@@ -55,6 +55,15 @@ See the end of the file for license conditions.
      --clumpsgeow3: Geometric center of all clumps in 3rd dim.
      --areaxy: Projected area in first two dimentions.
      --geoareaxy: Projected geoarea in first two dimensions.
+   --inbetweenints: output will contain one row for all integers between 1
+     and the largest label in the input (irrespective of their existance in
+     the input image). This was the default/only behavior of MakeCatalog
+     until now. However, there are situations where the labeled input image
+     integers may not be contiguous. For example if the input's only
+     labeled pixel values are 11 and 13 from this release MakeCatalog's
+     output will only have two rows. If you want the old behavior (of one
+     row per integer, even if its not in the image), you can use this
+     option.
 
   MakeProfiles:
    - Can produce mock ellipsoids in a datacube (using X-Z-X Euler angles
diff --git a/bin/mkcatalog/args.h b/bin/mkcatalog/args.h
index 1a2ad60..3bfdc7f 100644
--- a/bin/mkcatalog/args.h
+++ b/bin/mkcatalog/args.h
@@ -244,6 +244,19 @@ struct argp_option program_options[] =
       GAL_OPTIONS_NOT_MANDATORY,
       GAL_OPTIONS_NOT_SET
     },
+    {
+      "inbetweenints",
+      UI_KEY_INBETWEENINTS,
+      0,
+      0,
+      "Keep rows (integer ids) with no labels.",
+      GAL_OPTIONS_GROUP_OUTPUT,
+      &p->inbetweenints,
+      GAL_OPTIONS_NO_ARG_TYPE,
+      GAL_OPTIONS_RANGE_0_OR_1,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
+    },
 
 
 
diff --git a/bin/mkcatalog/main.h b/bin/mkcatalog/main.h
index 00eba86..2721c25 100644
--- a/bin/mkcatalog/main.h
+++ b/bin/mkcatalog/main.h
@@ -184,6 +184,7 @@ struct mkcatalogparams
   float           sfmagnsigma;  /* Surface brightness multiple of sigma.*/
   float             sfmagarea;  /* Surface brightness area (arcsec^2).  */
   uint8_t            spectrum;  /* Object spectrum for 3D datasets.     */
+  uint8_t       inbetweenints;  /* Keep rows (integer ids) with no labels. */
 
   char            *upmaskfile;  /* Name of upper limit mask file.       */
   char             *upmaskhdu;  /* HDU of upper limit mask file.        */
@@ -228,6 +229,7 @@ struct mkcatalogparams
   size_t         *numclumps_c;  /* To sort the clumps table by Obj.ID.  */
   gal_data_t   *specsliceinfo;  /* Slice information for spectra.       */
   gal_data_t         *spectra;  /* Array of datasets containing spectra.*/
+  gal_data_t      *rowsremove;  /* For rows/labels with no pixel.       */
 
   char        *usedvaluesfile;  /* Ptr to final name used for values.   */
   char        *usedclumpsfile;  /* Ptr to final name used for clumps.   */
diff --git a/bin/mkcatalog/mkcatalog.c b/bin/mkcatalog/mkcatalog.c
index 5d12634..dfd24dc 100644
--- a/bin/mkcatalog/mkcatalog.c
+++ b/bin/mkcatalog/mkcatalog.c
@@ -610,6 +610,47 @@ sort_clumps_by_objid(struct mkcatalogparams *p)
 
 
 
+static void
+mkcatalog_remove_labels_with_no_pixel(struct mkcatalogparams *p)
+{
+  size_t i;
+  gal_data_t *tmp, *perm;
+  uint8_t *rarray=p->rowsremove->array;
+  size_t g=0, b=0, numgood=0, *permutation;
+
+  /* Count how many good rows we have. */
+  for(i=0;i<p->numobjects;++i) if(rarray[i]==0) ++numgood;
+
+  /* Define the permutation to keep only the rows that must be shown. */
+  perm=gal_data_alloc(NULL, GAL_TYPE_SIZE_T, 1, &p->numobjects, NULL, 1,
+                      p->cp.minmapsize, p->cp.quietmmap, NULL, NULL, NULL);
+  permutation=perm->array;
+  for(i=0;i<p->numobjects;++i)
+    permutation[i] = rarray[i] ? numgood+b++ : g++;
+
+  /* For a check.
+  for(i=0;i<p->numobjects;++i)
+    printf("%zu: %u (%zu)\n", i+1, rarray[i], permutation[i]);
+  */
+
+  /* Apply the permutation. */
+  for(tmp=p->objectcols;tmp!=NULL;tmp=tmp->next)
+    {
+      /* Apply the permutation. */
+      gal_permutation_apply_inverse(tmp, permutation);
+
+      /* Correct the size. */
+      tmp->size=tmp->dsize[0]=numgood;
+    }
+
+  /* Clean up and return. */
+  gal_data_free(perm);
+}
+
+
+
+
+
 /* Write the produced columns into the output */
 static void
 mkcatalog_write_outputs(struct mkcatalogparams *p)
@@ -619,6 +660,11 @@ mkcatalog_write_outputs(struct mkcatalogparams *p)
   gal_list_str_t *comments;
   int outisfits=gal_fits_name_is_fits(p->objectsout);
 
+  /* If there are object labels without a pixel, then remove them from the
+     objects catalog. */
+  if(p->rowsremove)
+    mkcatalog_remove_labels_with_no_pixel(p);
+
   /* If a catalog is to be generated. */
   if(p->objectcols)
     {
diff --git a/bin/mkcatalog/ui.c b/bin/mkcatalog/ui.c
index 33d1514..780434c 100644
--- a/bin/mkcatalog/ui.c
+++ b/bin/mkcatalog/ui.c
@@ -1354,6 +1354,7 @@ ui_one_tile_per_object(struct mkcatalogparams *p)
 {
   size_t ndim=p->objects->ndim;
 
+  uint8_t *rarray=NULL;
   int32_t *l, *lf, *start;
   size_t i, d, *min, *max, width=2*ndim;
   size_t *minmax=gal_pointer_allocate(GAL_TYPE_SIZE_T,
@@ -1381,7 +1382,8 @@ ui_one_tile_per_object(struct mkcatalogparams *p)
     if(*l>0)
       {
         /* Get the coordinates of this pixel. */
-        gal_dimension_index_to_coord(l-start, ndim, p->objects->dsize, coord);
+        gal_dimension_index_to_coord(l-start, ndim, p->objects->dsize,
+                                     coord);
 
         /* Check to see this coordinate is the smallest/largest found so
            far for this label. Note that labels start from 1, while indexs
@@ -1396,14 +1398,38 @@ ui_one_tile_per_object(struct mkcatalogparams *p)
       }
   while(++l<lf);
 
-  /* If an label doesn't exist in the image, then set the minimum and
+  /* If a label doesn't exist in the image, then set the minimum and
      maximum values to zero. */
   for(i=0;i<p->numobjects;++i)
     for(d=0;d<ndim;++d)
       {
         if( minmax[ i * width + d ] == GAL_BLANK_SIZE_T
             && minmax[ i * width + ndim + d ] == 0 )
-          minmax[ i * width + d ] = minmax[ i * width + ndim + d ] = 0;
+          {
+            /* Set the first and last pixel to 0. */
+            minmax[ i * width + d ] = minmax[ i * width + ndim + d ] = 0;
+
+            /* Flag this label as one that must be removed. The ones that
+               must be removed get a value of `1' in `p->rowsremove'. */
+            if(p->inbetweenints==0)
+              {
+                if(p->rowsremove)
+                  rarray[i]=1;
+                else
+                  {
+                    /* Note that by initializing with zeros, all (the
+                       possibly existing) previous rows that shouldn't be
+                       removed are flagged as zero in this array. */
+                    p->rowsremove=gal_data_alloc(NULL, GAL_TYPE_UINT8, 1,
+                                                 &p->numobjects, NULL, 1,
+                                                 p->cp.minmapsize,
+                                                 p->cp.quietmmap,
+                                                 NULL, NULL, NULL);
+                    rarray=p->rowsremove->array;
+                    rarray[i]=1;
+                  }
+              }
+          }
       }
 
   /* For a check.
diff --git a/bin/mkcatalog/ui.h b/bin/mkcatalog/ui.h
index f4d2f05..c4c860f 100644
--- a/bin/mkcatalog/ui.h
+++ b/bin/mkcatalog/ui.h
@@ -92,6 +92,7 @@ enum option_keys_enum
   UI_KEY_SFMAGNSIGMA,
   UI_KEY_SFMAGAREA,
   UI_KEY_SPECTRUM,
+  UI_KEY_INBETWEENINTS,
   UI_KEY_UPMASKFILE,
   UI_KEY_UPMASKHDU,
   UI_KEY_UPNUM,
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 582e254..bd9d0b4 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -15426,6 +15426,13 @@ pixels that belong to other labels but overlap with 
the 2D projection. This
 can be used to see how reliable the emission line measurement is (on the
 projected spectra) and also if multiple lines (labeled regions) belong to
 the same physical object.
+
+@item --inbetweenints
+Output will contain one row for all integers between 1 and the largest label 
in the input (irrespective of their existance in the input image).
+By default, MakeCatalog's output will only contain rows with integers that 
actually corresponded to atleast one pixel in the input dataset.
+
+For example if the input's only labeled pixel values are 11 and 13, 
MakeCatalog's default output will only have two rows.
+If you use this option, it will have 13 rows and all the columns corresponding 
to integer identifiers that didn't correspond to any pixel will be 0 or NaN 
(depending on context).
 @end table