[gnuastro-commits] master e74f2b6 4/7: MakeCatalog makes clump catalog only when asked

[Mohammad Akhlaghi]

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

    MakeCatalog makes clump catalog only when asked
    
    MakeCatalogs will now read the `WCLUMPS' keyword in the header of the
    objects HDU to see if it should make a clump catalog or not. If a
    clump catalog isn't needed, then no input clump HDU will also be
    necessary.
    
    The value for `WCLUMPS' in the objects HDU should be a "yes" (not
    case-sensitive). If it is anything else, or doesn't exist, then it is
    assumed that no clump catalog is needed. To do the case insensitive check,
    it was necessary to use the `strcasecmp' function which is not fully
    portable in some systems, Gnulib's `strcase' module is now also imported
    during the bootstrapping.
    
    MakeCatalog (and more generally `gal_fits_read_keywords') can now deal
    with the non-existence of a header keyword. So the procedures to read
    other keywords were also updated: if they are not given, they are
    calculated (where possible). To do this in an effective and clean
    manner, a new `readkeywords' function was defined in MakeCatalog's
    `ui.c'.
    
    The detection and clump limiting S/N values (that can only be obtained
    by NoiseChisel) are now given a NaN value if they aren't present in
    the headers and will not be printed in the catalog's comments.
---
 bootstrap.conf            |    1 +
 src/mkcatalog/Makefile.am |    8 +--
 src/mkcatalog/mkcatalog.c |   43 ++++++-----
 src/mkcatalog/ui.c        |  175 +++++++++++++++++++++++++++++++++------------
 4 files changed, 162 insertions(+), 65 deletions(-)

diff --git a/bootstrap.conf b/bootstrap.conf
index 5908d99..88ca654 100644
--- a/bootstrap.conf
+++ b/bootstrap.conf
@@ -163,6 +163,7 @@ gnulib_modules="
     argp
     error
     nproc
+    strcase
     gendocs
     progname
     git-version-gen
diff --git a/src/mkcatalog/Makefile.am b/src/mkcatalog/Makefile.am
index 2f96074..2077c3c 100644
--- a/src/mkcatalog/Makefile.am
+++ b/src/mkcatalog/Makefile.am
@@ -25,12 +25,12 @@
 ## Utility and its sources
 bin_PROGRAMS = astmkcatalog
 
-astmkcatalog_SOURCES = main.c main.h cite.h ui.c ui.h args.h   \
+astmkcatalog_SOURCES = main.c main.h cite.h ui.c ui.h args.h          \
 mkcatalog.c mkcatalog.h columns.c columns.h
 
-astmkcatalog_LDADD = $(top_builddir)/bootstrapped/lib/libgnu.la        \
--lgalconfigfiles -lgalfits -lgalwcs -lgalcheckset -lgaltiming  \
--lgallinkedlist -lgaltxtarray
+astmkcatalog_LDADD = $(top_builddir)/bootstrapped/lib/libgnu.la       \
+-lgalconfigfiles -lgalstatistics -lgalarraymanip -lgalqsort -lgalfits \
+-lgalwcs -lgalcheckset -lgaltiming -lgallinkedlist -lgaltxtarray
 
 
 
diff --git a/src/mkcatalog/mkcatalog.c b/src/mkcatalog/mkcatalog.c
index af178c3..80352d9 100644
--- a/src/mkcatalog/mkcatalog.c
+++ b/src/mkcatalog/mkcatalog.c
@@ -154,7 +154,7 @@ firstpass(struct mkcatalogparams *p)
               }
           }
 
-        if(clumps[i]>0)
+        if(clumps && clumps[i]>0)
           {
             /* The largest clump ID over each object is the number of
                clumps that object has. */
@@ -191,11 +191,10 @@ firstpass(struct mkcatalogparams *p)
 
 
 
-/* In the second pass, we have the number of clumps so we can find
-   store the total values for the clumps. In this second round we can
-   also find second order moments of the objects if we want to. */
+/* In the second pass, we have the number of clumps so we can find the
+   total values for the clumps. */
 void
-secondpass(struct mkcatalogparams *p)
+clumppass(struct mkcatalogparams *p)
 {
   float imgss;
   double x, y, sx, sy;
@@ -412,6 +411,11 @@ makeoutput(struct mkcatalogparams *p)
   for(p->obj0clump1=0;p->obj0clump1<2;++p->obj0clump1)
     {
 
+      /* If no clumps image was provided, then ignore the clumps
+         catalog. */
+      if(p->obj0clump1==1 && p->clumps==NULL)
+        continue;
+
 
       /* Do the preparations for this round: */
       p->intcounter=p->accucounter=p->curcol=0;
@@ -483,13 +487,16 @@ makeoutput(struct mkcatalogparams *p)
       strcat(comment, p->line);
 
       sn = p->obj0clump1 ? p->clumpsn : p->detsn;
-      sprintf(tline, "%s limiting Signal to noise ratio: ",
-              p->obj0clump1 ? "Clump" : "Detection");
-      sprintf(p->line, "# "CATDESCRIPTLENGTH"%.3f\n", tline, sn);
-      strcat(comment, p->line);
-      if(p->obj0clump1==0)
-          strcat(comment, "# (NOTE: limits above are for detections, not "
-                 "objects)\n");
+      if(!isnan(sn))
+        {
+          sprintf(tline, "%s limiting Signal to noise ratio: ",
+                  p->obj0clump1 ? "Clump" : "Detection");
+          sprintf(p->line, "# "CATDESCRIPTLENGTH"%.3f\n", tline, sn);
+          strcat(comment, p->line);
+          if(p->obj0clump1==0)
+            strcat(comment, "# (NOTE: limits above are for detections, not "
+                   "objects)\n");
+        }
 
 
       /* If cpscorr was used, report it: */
@@ -518,8 +525,9 @@ makeoutput(struct mkcatalogparams *p)
       strcat(comment, "#\n# Columns:\n# --------\n");
 
 
-      /* Fill the catalog array, in the end set the last elements in intcols 
and
-         accucols to -1, so gal_txtarray_array_to_txt knows when to stop. */
+      /* Fill the catalog array, in the end set the last elements in
+         intcols and accucols to -1, so gal_txtarray_array_to_txt knows
+         when to stop. */
       for(p->curcol=0;p->curcol<p->numcols;++p->curcol)
         {
           col=cols[p->curcol];
@@ -687,8 +695,9 @@ makeoutput(struct mkcatalogparams *p)
 
 
       /* Write the catalog to file: */
-      gal_txtarray_array_to_txt(p->cat, p->num, p->numcols, comment, 
p->intcols,
-                                p->accucols, space, prec, 'f', p->filename);
+      gal_txtarray_array_to_txt(p->cat, p->num, p->numcols, comment,
+                                p->intcols, p->accucols, space, prec,
+                                'f', p->filename);
 
       /* Clean up: */
       free(p->intcols);
@@ -723,7 +732,7 @@ mkcatalog(struct mkcatalogparams *p)
 {
   /* Run through the data for the first time: */
   firstpass(p);
-  secondpass(p);
+  if(p->clumps) clumppass(p);
 
   /* Write the output: */
   makeoutput(p);
diff --git a/src/mkcatalog/ui.c b/src/mkcatalog/ui.c
index f982081..8035f7d 100644
--- a/src/mkcatalog/ui.c
+++ b/src/mkcatalog/ui.c
@@ -36,6 +36,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/timing.h>     /* Includes time.h and sys/time.h */
 #include <gnuastro/checkset.h>
 #include <gnuastro/txtarray.h>
+#include <gnuastro/statistics.h>
 #include <gnuastro/commonargs.h>
 #include <gnuastro/configfiles.h>
 
@@ -836,51 +837,43 @@ checkifset(struct mkcatalogparams *p)
 void
 sanitycheck(struct mkcatalogparams *p)
 {
-  struct gal_fits_read_header_keys keys[2];
+  struct uiparams *up=&p->up;
+  struct gal_fits_read_header_keys keys[1];
 
   /* Make sure the input file exists. */
-  gal_checkset_check_file(p->up.inputname);
+  gal_checkset_check_file(up->inputname);
 
   /* Set the names of the files. */
-  gal_fits_file_or_ext_name(p->up.inputname, p->cp.hdu, p->up.masknameset,
-                            &p->up.maskname, p->up.mhdu, p->up.mhduset,
+  gal_fits_file_or_ext_name(up->inputname, p->cp.hdu, up->masknameset,
+                            &up->maskname, up->mhdu, up->mhduset,
                             "mask");
-  gal_fits_file_or_ext_name(p->up.inputname, p->cp.hdu,
-                            p->up.objlabsnameset, &p->up.objlabsname,
-                            p->up.objhdu, p->up.objhduset,
+  gal_fits_file_or_ext_name(up->inputname, p->cp.hdu,
+                            up->objlabsnameset, &up->objlabsname,
+                            up->objhdu, up->objhduset,
                             "object labels");
-  gal_fits_file_or_ext_name(p->up.inputname, p->cp.hdu,
-                            p->up.clumplabsnameset, &p->up.clumplabsname,
-                            p->up.clumphdu, p->up.clumphduset,
-                            "clump labels");
-  gal_fits_file_or_ext_name(p->up.inputname, p->cp.hdu, p->up.skynameset,
-                            &p->up.skyname, p->up.skyhdu, p->up.skyhduset,
+  gal_fits_file_or_ext_name(up->inputname, p->cp.hdu, up->skynameset,
+                            &up->skyname, up->skyhdu, up->skyhduset,
                             "sky value image");
-  gal_fits_file_or_ext_name(p->up.inputname, p->cp.hdu, p->up.stdnameset,
-                            &p->up.stdname, p->up.stdhdu, p->up.stdhduset,
+  gal_fits_file_or_ext_name(up->inputname, p->cp.hdu, up->stdnameset,
+                            &up->stdname, up->stdhdu, up->stdhduset,
                             "sky standard deviation");
 
-  /* Read the number of labels for the objects:  */
-  keys[0].keyname="DETSN";        keys[0].datatype=TDOUBLE;
-  keys[1].keyname="NOBJS";        keys[1].datatype=TLONG;
-  gal_fits_read_keywords(p->up.objlabsname, p->up.objhdu, keys, 2);
-  p->detsn=keys[0].d;
-  p->numobjects=keys[1].l;
-
-  /* Read the clumps information. Note that the datatypes don't change. */
-  keys[0].keyname="CLUMPSN";
-  keys[1].keyname="NCLUMPS";
-  gal_fits_read_keywords(p->up.clumplabsname, p->up.clumphdu, keys, 2);
-  p->clumpsn=keys[0].d;
-  p->numclumps=keys[1].l;
-
-  /* Read the minimum and maximum standard deviation values. */
-  keys[0].keyname="MINSTD";       keys[0].datatype=TFLOAT;
-  keys[1].keyname="MEDSTD";       keys[1].datatype=TFLOAT;
-  gal_fits_read_keywords(p->up.stdname, p->up.stdhdu, keys, 2);
-  p->minstd=keys[0].f;
-  p->medstd=keys[1].f;
-  p->cpscorr = p->minstd>1 ? 1.0f : p->minstd;
+  /* The WCLUMPS (with-clumps) keyword in the object HDU says that there is
+     also a clumps image accompaning the object image. If it exists and its
+     value is `yes' (not case sensitive), then set the image name to the
+     proper string. Otherwise, set the image name to NULL, so future
+     functions can check. */
+  keys[0].str[0]='\0';
+  keys[0].datatype=TSTRING;
+  keys[0].keyname="WCLUMPS";
+  gal_fits_read_keywords(up->objlabsname, up->objhdu, keys, 1);
+  if(strcasecmp(keys[0].str, "yes"))
+    up->clumplabsname=NULL;
+  else
+    gal_fits_file_or_ext_name(up->inputname, p->cp.hdu,
+                              up->clumplabsnameset, &up->clumplabsname,
+                              up->clumphdu, up->clumphduset,
+                              "clump labels");
 
   /* When the RA and Dec are needed, make sure that the X and Y
      columns and the RA and Dec columns in the information array are
@@ -891,7 +884,7 @@ sanitycheck(struct mkcatalogparams *p)
 
      NOTE: the information array is separate from the output array
   */
-  if(p->up.raset || p->up.decset)
+  if(up->raset || up->decset)
     {
       if( OFlxWhtX!=OFlxWhtY-1 || OFlxWhtRA!=OFlxWhtDec-1 )
         error(EXIT_FAILURE, 0, "a bug! Please contact us at %s so we can "
@@ -909,10 +902,10 @@ sanitycheck(struct mkcatalogparams *p)
     }
   else
     {
-      gal_checkset_automatic_output(p->up.inputname, "_o.txt",
+      gal_checkset_automatic_output(up->inputname, "_o.txt",
                                     p->cp.removedirinfo, p->cp.dontdelete,
                                     &p->ocatname);
-      gal_checkset_automatic_output(p->up.inputname, "_c.txt",
+      gal_checkset_automatic_output(up->inputname, "_c.txt",
                                     p->cp.removedirinfo, p->cp.dontdelete,
                                     &p->ccatname);
     }
@@ -998,6 +991,92 @@ checksetfloat(struct mkcatalogparams *p, char *filename, 
char *hdu,
 
 
 
+/* Read the necessary keywords and if they aren't present do the
+   appropriate action. */
+void
+readkeywords(struct mkcatalogparams *p)
+{
+  long numobjects;
+  size_t size=p->s0*p->s1;
+  struct uiparams *up=&p->up;
+  struct gal_fits_read_header_keys keys[2];
+
+  /* Read keywords from the standard deviation image:  */
+  keys[0].keyname="MINSTD";   keys[0].datatype=TFLOAT;
+  keys[1].keyname="MEDSTD";   keys[1].datatype=TFLOAT;
+  gal_fits_read_keywords(up->stdname, up->stdhdu, keys, 2);
+
+  /* The minimum standard deviation value. */
+  if(keys[0].status)
+    gal_statistics_float_min(p->std, size, &p->minstd);
+  else
+    p->minstd=keys[0].f;
+  p->cpscorr = p->minstd>1 ? 1.0f : p->minstd;
+
+  /* Median standard deviation value (only necessary in catalog
+     comments). */
+  if(keys[1].status)
+    {
+      p->medstd=gal_statistics_median(p->std, size);
+      fprintf(stderr, "Warning: Could not find the MEDSTD keyword in %s "
+              "(hdu: %s). The median standard deviation is thus found on"
+              "the (interpolated) standard deviation image. NoiseChisel "
+              "finds the median before interpolation, so the reported value "
+              "in the final catalog will not be accurate (will depend on "
+              "how many meshs were blank and their spatial position "
+              "relative to the non-blank ones.", up->stdname, up->stdhdu);
+    }
+  else
+    p->medstd=keys[1].f;
+
+  /* Read the keywords from the objects image:  */
+  keys[0].keyname="DETSN";     keys[0].datatype=TDOUBLE;
+  keys[1].keyname="NOBJS";     keys[1].datatype=TLONG;
+  gal_fits_read_keywords(up->objlabsname, up->objhdu, keys, 2);
+
+  /* If DETSN is not given, there is no way we can calculate it here, so we
+     will just set it to NaN to check and not report in the end. */
+  p->detsn = keys[0].status ? NAN : keys[0].d;
+
+  /* Read the total number of objects. */
+  if(keys[1].status)
+    {
+      gal_statistics_long_non_blank_max(p->objects, size, &numobjects,
+                                        GAL_FITS_LONG_BLANK);
+      p->numobjects=numobjects;
+    }
+  else
+    p->numobjects=keys[1].l;
+
+  /* Read the clumps information if it is necessary.
+
+     Note that unlike finding the number of objects, finding the number of
+     clumps is not an easy process, since the clumps of each object start
+     with a label of one. So if the number of clumps is not given, we have
+     to abort. For now, is up to the program that built the clumps to give
+     the total number, later we can take a procedure to find them (for
+     example first only taking the positive labels (that are clumps) and
+     making a binary image, then running the connected component algorithm
+     to find the number of clumps in the image.*/
+  if(up->clumplabsname)
+    {
+      keys[0].keyname="CLUMPSN";
+      keys[1].keyname="NCLUMPS";
+      gal_fits_read_keywords(up->clumplabsname, up->clumphdu, keys, 2);
+      p->clumpsn = keys[0].status ? NAN : keys[0].d;
+      if(keys[1].status)
+        error(EXIT_FAILURE, 0, "couldn't find NCLUMPS in the header of "
+              "%s (hdu: %s).", up->clumplabsname, up->clumphdu);
+      else
+        p->numclumps=keys[1].l;
+    }
+
+}
+
+
+
+
+
 
 void
 preparearrays(struct mkcatalogparams *p)
@@ -1185,12 +1264,19 @@ preparearrays(struct mkcatalogparams *p)
                         &p->wcs);
 
 
-      /* Read and check the other arrays: */
-      checksetlong(p, p->up.objlabsname, p->up.objhdu, &p->objects);
-      checksetlong(p, p->up.clumplabsname, p->up.clumphdu, &p->clumps);
+      /* Read and check the object, sky and skystd HDUs. Note that the
+         clumps image is only used when the objects image says a clumps
+         image exists. */
       checksetfloat(p, p->up.skyname, p->up.skyhdu, &p->sky);
       checksetfloat(p, p->up.stdname, p->up.stdhdu, &p->std);
+      checksetlong(p, p->up.objlabsname, p->up.objhdu, &p->objects);
+      if(p->up.clumplabsname)
+        checksetlong(p, p->up.clumplabsname, p->up.clumphdu, &p->clumps);
+      else
+        p->clumps=NULL;
 
+      /* Read the necessary keywords. */
+      readkeywords(p);
 
       /* Allocate the catalog arrays: */
       if(p->objncols>0 && p->numobjects>0)
@@ -1312,8 +1398,9 @@ setparams(int argc, char *argv[], struct mkcatalogparams 
*p)
         printf("  - Mask   %s (hdu: %s)\n", p->up.maskname, p->up.mhdu);
       printf("  - Objects %s (hdu: %s)\n", p->up.objlabsname,
              p->up.objhdu);
-      printf("  - Clumps  %s (hdu: %s)\n", p->up.clumplabsname,
-             p->up.clumphdu);
+      if(p->up.clumplabsname)
+        printf("  - Clumps  %s (hdu: %s)\n", p->up.clumplabsname,
+               p->up.clumphdu);
       printf("  - Sky     %s (hdu: %s)\n", p->up.skyname, p->up.skyhdu);
       printf("  - Sky STD %s (hdu: %s)\n", p->up.stdname, p->up.stdhdu);
     }