[gnuastro-commits] master 6bd397a1: MakeCatalog: same input clump labs when NUMLABS not in clumps image

[Mohammad Akhlaghi]

branch: master
commit 6bd397a19311952f6a327d55f7a12aa1c2b45298
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    MakeCatalog: same input clump labs when NUMLABS not in clumps image
    
    Until now, when the 'NUMLABS' keyword was not in the clumps input of the
    MakeCatalog, it would reset the labels to start from one. This was very
    inconvenient and required the user to map the new labels with the original
    labels to use them.
    
    With this commit, the internal mapping is kept within MakeCatalog and is
    used to put the original labels in the clumps afterwards.
---
 NEWS                    |   9 +++
 bin/mkcatalog/columns.c |   4 +-
 bin/mkcatalog/main.h    |   2 +-
 bin/mkcatalog/ui.c      | 143 ++++++++++++++++++++++++++++++++++++------------
 doc/gnuastro.texi       |  12 +---
 5 files changed, 122 insertions(+), 48 deletions(-)

diff --git a/NEWS b/NEWS
index 25da8461..08cdc0ff 100644
--- a/NEWS
+++ b/NEWS
@@ -133,6 +133,15 @@ See the end of the file for license conditions.
     what is expected by default). Recall that if the sky is not subtracted
     from the values image you can use '--sky' for MakeCatalog to subtract
     it before the measurements. Implemented by Sepideh Eskandarlou.
+  - When the input clump labels do not have the 'NUMLABS' keyword
+    (specifying how many clumps there are in the image, independent of the
+    number of objects), MakeCatalog will no longer re-label the input
+    clumps: the output catalog will have the same clump ids as the
+    labels. Until now, it would generate new labels for the clumps in such
+    cases (along with a new FITS image containing the clumps with the new
+    labels to allow you to match them). The output catalog also had the new
+    label, not the original labels (the old behavior was therefore very
+    inconvenient!).
 
   Segment
   - The default name of the first extension of Segment's output is now
diff --git a/bin/mkcatalog/columns.c b/bin/mkcatalog/columns.c
index d598b047..92d03cd2 100644
--- a/bin/mkcatalog/columns.c
+++ b/bin/mkcatalog/columns.c
@@ -3159,7 +3159,9 @@ columns_fill(struct mkcatalog_passparams *pp)
             break;
 
           case UI_KEY_IDINHOSTOBJ:
-            ((int32_t *)colarr)[cind]=coind+1;
+            ((int32_t *)colarr)[cind] = ( p->origclpid
+                                          ? p->origclpid[pp->object][coind]
+                                          : coind+1 );
             break;
 
           case UI_KEY_AREA:
diff --git a/bin/mkcatalog/main.h b/bin/mkcatalog/main.h
index 14b036a2..4f7f468a 100644
--- a/bin/mkcatalog/main.h
+++ b/bin/mkcatalog/main.h
@@ -281,7 +281,6 @@ struct mkcatalogparams
   float     spatialresolution;  /* Error in area (used in SB error).    */
 
   /* Internal. */
-  char           *relabclumps;  /* Name of new file for clump labels.   */
   time_t              rawtime;  /* Starting time of the program.        */
   gal_data_t          *values;  /* Input.                               */
   gal_data_t         *objects;  /* Object labels.                       */
@@ -293,6 +292,7 @@ struct mkcatalogparams
   float               cpscorr;  /* Counts-per-second correction.        */
   int32_t            *outlabs;  /* Labels in output cat (when necessary)*/
   int32_t         *outlabsinv;  /* Inverse of the 'outlabs' array.      */
+  int32_t         **origclpid;  /* Original clump labels in objects.    */
   size_t           numobjects;  /* Number of object labels in image.    */
   float               clumpsn;  /* Clump S/N threshold.                 */
   size_t            numclumps;  /* Number of clumps in image.           */
diff --git a/bin/mkcatalog/ui.c b/bin/mkcatalog/ui.c
index c9e3136b..6552c4a2 100644
--- a/bin/mkcatalog/ui.c
+++ b/bin/mkcatalog/ui.c
@@ -575,26 +575,35 @@ ui_wcs_info(struct mkcatalogparams *p)
 
 
 
+/* When the number of clumps is not specified a priori (through the
+   'NUMLABS' keyword in the clumps array, we need to map the original
+   labels to new labels that start counting from one and are
+   congtiguous. */
 static size_t
-ui_num_clumps(struct mkcatalogparams *p)
+ui_num_clumps_orig_labs(struct mkcatalogparams *p)
 {
-  int32_t olab;
-  char *basename;
-  int keepinputdir;
-  size_t i, counter, numclumps=0;
-  gal_list_i32_t *tmp, **labsinobj;
-  int32_t *o=p->objects->array, *of=o+p->objects->size, *c=p->clumps->array;
-
-  /* Allocate array of lists to keep the unique labels within each
-     object. */
+  size_t ncinobj;
+  gal_list_i32_t *tmp, **origclumplist;
+  int32_t olab, **origclpid, *c=p->clumps->array;
+  size_t i, j, numclumps=0, nobjpone=p->numobjects+1;
+  int32_t *o=p->objects->array, *of=o+p->objects->size;
+
+  /* Allocate array of lists and arrays to keep the unique labels within
+     each object. The list is used in the first run and the array is later
+     filled based on the list.*/
+  errno=0;
+  origclumplist=calloc(nobjpone, sizeof *origclumplist);
+  if(origclumplist==NULL)
+    error(EXIT_FAILURE, 0, "%s: couldn't allocate %zu bytes for "
+          "'origclumplist'", __func__, nobjpone * sizeof *origclumplist);
   errno=0;
-  labsinobj=calloc(p->numobjects+1, sizeof *labsinobj);
-  if(labsinobj==NULL)
+  origclpid=calloc(nobjpone, sizeof *origclpid);
+  if(origclpid==NULL)
     error(EXIT_FAILURE, 0, "%s: couldn't allocate %zu bytes for "
-          "'labsinobj'", __func__, p->numobjects * sizeof *labsinobj);
+          "'origclpid'", __func__, nobjpone * sizeof *origclpid);
+  p->origclpid=origclpid; /* To avoid 'p->' in the lines below. */
 
-  /* Go over each pixel and find the unique clump labels within each
-     object. */
+  /* Go over each pixel and find the list of original clump labels. */
   do
     {
       /* Do the steps if we are on a clump. */
@@ -602,14 +611,14 @@ ui_num_clumps(struct mkcatalogparams *p)
         {
           /* See if the label has already been found. */
           olab = p->outlabsinv ? p->outlabsinv[*o] : *o;
-          for(tmp=labsinobj[olab];tmp!=NULL;tmp=tmp->next)
+          for(tmp=origclumplist[olab];tmp!=NULL;tmp=tmp->next)
             if(tmp->v==*c) break;
 
           /* When it wasn't found, 'tmp==NULL'. */
           if(tmp==NULL)
             {
               ++numclumps;
-              gal_list_i32_add(&labsinobj[olab], *c);
+              gal_list_i32_add(&origclumplist[olab], *c);
             }
         }
 
@@ -618,8 +627,67 @@ ui_num_clumps(struct mkcatalogparams *p)
     }
   while(++o<of);
 
-  /* Re-write the clump values so their numbering is contiguous, since this
-     is assumed during the later steps. */
+  /* For a check of the list.
+  for(i=1;i<p->numobjects+1;++i)
+    {
+      printf("list%zu: ", i);
+      for(tmp=origclumplist[i];tmp!=NULL;tmp=tmp->next)
+        printf("%d%s", tmp->v, tmp->next?" ":"");
+      printf("\n");
+    }
+  */
+
+  /* Build a separate array of original labels for each objects. */
+  for(i=1;i<p->numobjects+1;++i)
+    {
+      /* Find the number of clumps in this object and allocate. */
+      ncinobj=gal_list_i32_number(origclumplist[i]);
+      origclpid[i]=gal_pointer_allocate(GAL_TYPE_INT32, ncinobj+1, 0,
+                                          __func__, "origclpid");
+
+      /* Put the list into the array in reverse because it was filled as
+         last-in-last-out. We are also setting the last element to blank so
+         we can parse through the labels when necessary. */
+      j=ncinobj;
+      origclpid[i][j--]=GAL_BLANK_INT32;
+      for(tmp=origclumplist[i];tmp!=NULL;tmp=tmp->next)
+        origclpid[i][j--]=tmp->v;
+    }
+
+  /* For a check of the array.
+  for(i=1;i<p->numobjects+1;++i)
+    if(i==9)
+    {
+      printf("arr%zu: ", i);
+      for(j=0; origclpid[i][j]!=GAL_BLANK_INT32; ++j)
+        printf("%d%s", origclpid[i][j],
+               origclpid[i][j+1]==GAL_BLANK_INT32?"":" ");
+      printf("\n");
+    }
+  */
+
+  /* Clean up and return the total number of clumps. */
+  for(i=1;i<p->numobjects;++i) gal_list_i32_free(origclumplist[i]);
+  free(origclumplist);
+  return numclumps;
+}
+
+
+
+
+static size_t
+ui_num_clumps(struct mkcatalogparams *p)
+{
+  size_t i, numclumps=0;
+  int32_t olab, *c=p->clumps->array;
+  int32_t *o=p->objects->array, *of=o+p->objects->size;
+
+  /* Count the total number of clumps and extract the original labels
+     within the clumps into an array ('p->origclpid'). */
+  numclumps=ui_num_clumps_orig_labs(p);
+
+  /* Re-write the clump values so their numbering is contiguous (this is
+     assumed during the later steps). */
   o=p->objects->array;
   c=p->clumps->array;
   do
@@ -627,10 +695,9 @@ ui_num_clumps(struct mkcatalogparams *p)
       /* Do the steps if we are on a clump. */
       if(*o>0 && *c>0)
         {
-          counter=0;
           olab = p->outlabsinv ? p->outlabsinv[*o] : *o;
-          for(tmp=labsinobj[olab];tmp!=NULL;tmp=tmp->next)
-            { counter++; if(tmp->v==*c) {*c=counter; break;} }
+          for(i=0;p->origclpid[olab][i]!=GAL_BLANK_INT32;++i)
+            if(p->origclpid[olab][i]==*c) { *c=i+1; break; }
         }
 
       /* Increment the clumps pointer.*/
@@ -638,19 +705,18 @@ ui_num_clumps(struct mkcatalogparams *p)
     }
   while(++o<of);
 
-  /* Write the created file into a file for the user to inspect. */
+  /* For a check
+  char *basename;
+  char *checkname;
+  int keepinputdir;
   keepinputdir=p->cp.keepinputdir;
   p->cp.keepinputdir = p->cp.output ? 1 : 0;
   basename = p->cp.output ? p->cp.output : p->objectsfile;
-  p->relabclumps=gal_checkset_automatic_output(&p->cp, basename,
-                                               "-clumps-relab.fits");
-  gal_fits_img_write(p->clumps, p->relabclumps, NULL, PROGRAM_STRING);
+  checkname=gal_checkset_automatic_output(&p->cp, basename,
+                                          "-clumps-relab.fits");
+  gal_fits_img_write(p->clumps, checkname, NULL, PROGRAM_STRING);
   p->cp.keepinputdir=keepinputdir;
-
-  /* Clean up. */
-  for(i=0;i<p->numobjects;++i)
-    if(labsinobj[i]) gal_list_i32_free(labsinobj[i]);
-  free(labsinobj);
+  */
 
   /* Return the number of clumps. */
   return numclumps;
@@ -1898,9 +1964,6 @@ ui_read_check_inputs_setup(int argc, char *argv[],
       if(p->clumps)
         printf("  - Clumps:  %s (hdu: %s)\n", p->usedclumpsfile,
                p->clumpshdu);
-      if(p->relabclumps)
-        printf("  - RELABELED CLUMPS (no NUMLABS in original): %s\n",
-               p->relabclumps);
       if(p->values)
         printf("  - Values:  %s (hdu: %s)\n", p->usedvaluesfile,
                p->valueshdu);
@@ -1964,7 +2027,7 @@ ui_read_check_inputs_setup(int argc, char *argv[],
 void
 ui_free_report(struct mkcatalogparams *p, struct timeval *t1)
 {
-  size_t d;
+  size_t i, d;
 
   /* The temporary arrays for WCS coordinates. */
   if(p->wcs_vo ) gal_list_data_free(p->wcs_vo);
@@ -2017,6 +2080,16 @@ ui_free_report(struct mkcatalogparams *p, struct timeval 
*t1)
   if(p->upcheckout) free(p->upcheckout);
   gal_data_array_free(p->tiles, p->numobjects, 0);
 
+  /* The original clump labels. */
+  if(p->origclpid)
+    {
+      for(i=1;i<p->numobjects+1;++i)
+        if(p->origclpid[i])
+          free(p->origclpid[i]);
+      free(p->origclpid);
+    }
+  /**/
+
   /* If the Sky or its STD image were given in tiles, then we defined a
      tile structure to deal with them. The initialization of the tile
      structure is checked with its 'ndim' element. */
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 95ab25dd..bfe56e2c 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -28831,19 +28831,9 @@ When each dataset is in a separate file, these options 
are necessary.
 However, one great advantage of the FITS file format (that is heavily used in 
astronomy) is that it allows the storage of multiple datasets in one file.
 So in most situations (for example, if you are using the outputs of 
@ref{NoiseChisel} or @ref{Segment}), all the necessary input datasets can be in 
one file.
 
-When none of the @option{--*file} options are given, MakeCatalog will assume 
the necessary input datasets are in the file given as its argument (without any 
option).
+When none of the @option{--*file} options are given (for example 
@option{--clumpsfile} or @option{--valuesfile}), MakeCatalog will assume the 
necessary input datasets are available as HDUs in the file given as its 
argument (without any option).
 When the Sky or Sky standard deviation datasets are necessary and the only 
@option{--*file} option called is @option{--valuesfile}, MakeCatalog will 
search for these datasets (with the default/given HDUs) in the file given to 
@option{--valuesfile} (before looking into the main argument file).
 
-When the clumps image (necessary with the @option{--clumpscat} option) is 
used, MakeCatalog looks into the (possibly existing) @code{NUMLABS} keyword for 
the total number of clumps in the image (irrespective of how many objects there 
are).
-If it is not present, it will count them and possibly re-label the clumps so 
the clump labels always start with 1 and finish with the total number of clumps 
in each object.
-The re-labeled clumps image will be stored with the @file{-clumps-relab.fits} 
suffix.
-This can slightly slow-down the run.
-
-Note that @code{NUMLABS} is automatically written by Segment in its outputs, 
so if you are feeding Segment's clump labels, you can benefit from the improved 
speed.
-Otherwise, if you are creating the clumps label dataset manually, it may be 
good to include the @code{NUMLABS} keyword in its header and also be sure that 
there is no gap in the clump labels.
-For example, if an object has three clumps, they are labeled as 1, 2, 3.
-If they are labeled as 1, 3, 4, or any other combination of three positive 
integers that are not an increment of the previous, you might get unknown 
behavior.
-
 It may happen that your labeled objects image was created with a program that 
only outputs floating point files.
 However, you know it only has integer valued pixels that are stored in a 
floating point container.
 In such cases, you can use Gnuastro's Arithmetic program (see 
@ref{Arithmetic}) to change the numerical data type of the image 
(@file{float.fits}) to an integer type image (@file{int.fits}) with a command 
like below: