[gnuastro-commits] master 924b90e: NoiseChisel: New options to remove hard-coded configuration steps

[Mohammad Akhlaghi]

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

    NoiseChisel: New options to remove hard-coded configuration steps
    
    Some steps were hardcoded in Gnuastro's code and thus couldn't be modified
    by the user. With this commit, several of these configuration steps have
    been opened up with the new `--dopening', `--dopeningngb'.
    
    Also, the `--snthresh' option was added to allow the user to skip the step
    to identify false pseudo-detections.
---
 NEWS                                | 20 ++++++++---
 bin/noisechisel/args.h              | 41 +++++++++++++++++++++-
 bin/noisechisel/astnoisechisel.conf |  2 ++
 bin/noisechisel/detection.c         | 68 +++++++++++++++++++++----------------
 bin/noisechisel/main.h              |  5 ++-
 bin/noisechisel/ui.c                |  5 +++
 bin/noisechisel/ui.h                |  3 ++
 doc/gnuastro.texi                   | 34 +++++++++++++++++++
 8 files changed, 141 insertions(+), 37 deletions(-)

diff --git a/NEWS b/NEWS
index b6c2caf..ab729be 100644
--- a/NEWS
+++ b/NEWS
@@ -51,11 +51,21 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
      columns into one.
 
   NoiseChisel:
-   --pseudoconcomp: allows setting the connectivity (4 or 8, in a 2D image)
-     to define separate pseudo-detections. If its stronger,
-     pseudo-detections that are touching on the corner will be identified
-     as one. Until this version, this was hard-written into the code and
-     was the weakest connectivity (4-connected in a 2D image).
+
+   --snthresh: Manually set the signal-to-noise ratio of true
+     pseudo-detections. With this option, NoiseChisel will not attempt to
+     find pseudo-detections over the noisy regions of the dataset, but will
+     directly go onto applying the manually input value.
+   - Several new options have been added to remove NoiseChisel's dependence
+     on values that were hard-coded in its source and thus not modifiable
+     at run-time by the user. To allow full configurability, these steps
+     can also be configured by the user.
+     --pseudoconcomp: allows setting the connectivity (4 or 8, in a 2D
+       image) to define separate pseudo-detections. If its stronger,
+       pseudo-detections that are touching on the corner will be identified
+       as one.
+     --dopening: The number of openings to do after applying `--dthresh'.
+     --dopeningngb: The connectivity (4 or 8) to use for `--dopening'.
 
   Library:
     GAL_BLANK_LONG: new macro for the `long' type blank value.
diff --git a/bin/noisechisel/args.h b/bin/noisechisel/args.h
index 19bb426..8b85805 100644
--- a/bin/noisechisel/args.h
+++ b/bin/noisechisel/args.h
@@ -325,7 +325,7 @@ struct argp_option program_options[] =
       UI_GROUP_DETECTION,
       &p->opening,
       GAL_TYPE_SIZE_T,
-      GAL_OPTIONS_RANGE_GT_0,
+      GAL_OPTIONS_RANGE_GE_0,
       GAL_OPTIONS_MANDATORY,
       GAL_OPTIONS_NOT_SET
     },
@@ -409,6 +409,32 @@ struct argp_option program_options[] =
       GAL_OPTIONS_NOT_SET
     },
     {
+      "dopening",
+      UI_KEY_DOPENING,
+      "INT",
+      0,
+      "Depth of opening after dthresh.",
+      UI_GROUP_DETECTION,
+      &p->dopening,
+      GAL_TYPE_SIZE_T,
+      GAL_OPTIONS_RANGE_GE_0,
+      GAL_OPTIONS_MANDATORY,
+      GAL_OPTIONS_NOT_SET
+    },
+    {
+      "dopeningngb",
+      UI_KEY_DOPENINGNGB,
+      "INT",
+      0,
+      "4 or 8 connectivity for dthresh opening.",
+      UI_GROUP_DETECTION,
+      &p->dopeningngb,
+      GAL_TYPE_SIZE_T,
+      GAL_OPTIONS_RANGE_GT_0,
+      GAL_OPTIONS_MANDATORY,
+      GAL_OPTIONS_NOT_SET
+    },
+    {
       "holengb",
       UI_KEY_HOLENGB,
       "INT",
@@ -487,6 +513,19 @@ struct argp_option program_options[] =
       GAL_OPTIONS_NOT_SET
     },
     {
+      "snthresh",
+      UI_KEY_SNTHRESH,
+      "FLT",
+      0,
+      "Manually input pseudo-det S/N threshold.",
+      UI_GROUP_DETECTION,
+      &p->snthresh,
+      GAL_TYPE_FLOAT32,
+      GAL_OPTIONS_RANGE_GE_0,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
+    },
+    {
       "detgrowquant",
       UI_KEY_DETGROWQUANT,
       "FLT",
diff --git a/bin/noisechisel/astnoisechisel.conf 
b/bin/noisechisel/astnoisechisel.conf
index b872b76..cbb4596 100644
--- a/bin/noisechisel/astnoisechisel.conf
+++ b/bin/noisechisel/astnoisechisel.conf
@@ -39,6 +39,8 @@
  minskyfrac            0.7
  sigmaclip           3,0.2
  dthresh               0.0
+ dopening                1
+ dopeningngb             4
  holengb                 8
  pseudoconcomp           4
  snminarea              10
diff --git a/bin/noisechisel/detection.c b/bin/noisechisel/detection.c
index 1af2355..4552c7a 100644
--- a/bin/noisechisel/detection.c
+++ b/bin/noisechisel/detection.c
@@ -96,7 +96,7 @@ detection_initial(struct noisechiselparams *p)
   if(!p->cp.quiet)
     {
       if( asprintf(&msg, "Eroded %zu time%s (%zu-connected).", p->erode,
-                   p->erode>1?"s":"", p->erodengb)<0 )
+                   p->erode!=1?"s":"", p->erodengb)<0 )
         error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
       gal_timing_report(&t1, msg, 2);
       free(msg);
@@ -293,7 +293,7 @@ detection_fill_holes_open(void *in_prm)
         }
 
       /* Open all the regions. */
-      gal_binary_open(copy, 1, 1, 1);
+      gal_binary_open(copy, p->dopening, p->dopeningngb==4 ? 1 : 2, 1);
 
       /* Write the copied region back into the large input and AFTERWARDS,
          correct the tile's pointers, the pointers must not be corrected
@@ -762,36 +762,44 @@ detection_pseudo_real(struct noisechiselparams *p)
   workbin->flag=p->input->flag;
 
 
-  /* Over the Sky: find the pseudo-detections and make the S/N table. */
-  if(!p->cp.quiet) gettimeofday(&t1, NULL);
-  numpseudo=detection_pseudo_find(p, workbin, worklab, 0);
-  sn=detection_sn(p, worklab, numpseudo, 0, "PSEUDOS-FOR-SN");
-
-
-  /* A small sanity check */
-  if( sn->size < p->minnumfalse)
-    error(EXIT_FAILURE, 0, "only %zu pseudo-detections could be found over "
-          "the sky region to estimate an S/N. This is less than %zu (value "
-          "to `--minnumfalse' option). Please adjust parameters like "
-          "`--dthresh', `--snminarea', or make sure that there actually "
-          "is sufficient sky area after initial detection. You can use "
-          "`--checkdetection' to see every step until this point", sn->size,
-          p->minnumfalse);
-
-
-  /* Get the S/N quantile and report it if we are in non-quiet mode. */
-  quant=gal_statistics_quantile(sn, p->snquant, 1);
-  p->detsnthresh = *((float *)(quant->array));
-  if(!p->cp.quiet)
+  /* If a manual S/N is given, then don't bother with using
+     pseudo-detections. */
+  if( isnan(p->snthresh) )
     {
-      if( asprintf(&msg, "Pseudo-det S/N: %.3f (%g quant of %zu).",
-                   p->detsnthresh, p->snquant, sn->size)<0 )
-        error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
-      gal_timing_report(&t1, msg, 2);
-      free(msg);
+      /* Over the Sky: find the pseudo-detections and make the S/N
+         table. */
+      if(!p->cp.quiet) gettimeofday(&t1, NULL);
+      numpseudo=detection_pseudo_find(p, workbin, worklab, 0);
+      sn=detection_sn(p, worklab, numpseudo, 0, "PSEUDOS-FOR-SN");
+
+
+      /* A small sanity check */
+      if( sn->size < p->minnumfalse)
+        error(EXIT_FAILURE, 0, "only %zu pseudo-detections could be found "
+              "over the sky region to estimate an S/N. This is less than "
+              "%zu (value to `--minnumfalse' option). Please adjust "
+              "parameters like `--dthresh', `--snminarea', or make sure "
+              "that there actually is sufficient sky area after initial "
+              "detection. You can use `--checkdetection' to see every step "
+              "until this point", sn->size, p->minnumfalse);
+
+
+      /* Get the S/N quantile and report it if we are in non-quiet mode. */
+      quant=gal_statistics_quantile(sn, p->snquant, 1);
+      p->detsnthresh = *((float *)(quant->array));
+      if(!p->cp.quiet)
+        {
+          if( asprintf(&msg, "Pseudo-det S/N: %.3f (%g quant of %zu).",
+                       p->detsnthresh, p->snquant, sn->size)<0 )
+            error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+          gal_timing_report(&t1, msg, 2);
+          free(msg);
+        }
+      gal_data_free(sn);
+      gal_data_free(quant);
     }
-  gal_data_free(sn);
-  gal_data_free(quant);
+  else
+    p->detsnthresh=p->snthresh;
 
 
   /* Over the detections: find pseudo-detections and make S/N table. */
diff --git a/bin/noisechisel/main.h b/bin/noisechisel/main.h
index b38686e..f0d62c5 100644
--- a/bin/noisechisel/main.h
+++ b/bin/noisechisel/main.h
@@ -73,12 +73,15 @@ struct noisechiselparams
   double         sigmaclip[2];  /* Sigma-clipping parameters.             */
   uint8_t         checkdetsky;  /* Check pseudo-detection sky value.      */
   float               dthresh;  /* Sigma threshold for Pseudo-detections. */
+  size_t             dopening;  /* Depth of opening after dthresh.        */
+  size_t          dopeningngb;  /* Connectivity for opening after dthresh.*/
   size_t              holengb;  /* Connectivity for defining a hole.      */
   size_t        pseudoconcomp;  /* Connectivity for connected components. */
   size_t            snminarea;  /* Minimum pseudo-detection area for S/N. */
   uint8_t             checksn;  /* Save pseudo-detection S/N values.      */
   size_t          minnumfalse;  /* Min No. of det/seg for true quantile.  */
-  float              snquant;  /* True detection quantile.               */
+  float               snquant;  /* True detection quantile.               */
+  float              snthresh;  /* Manually input S/N value.              */
   float          detgrowquant;  /* Quantile to grow true detections.      */
   size_t   detgrowmaxholesize;  /* Max. size of holes to fill in growth.  */
   uint8_t       cleangrowndet;  /* Remove grown objects with small S/N.   */
diff --git a/bin/noisechisel/ui.c b/bin/noisechisel/ui.c
index e66a8d2..318963a 100644
--- a/bin/noisechisel/ui.c
+++ b/bin/noisechisel/ui.c
@@ -115,6 +115,8 @@ ui_initialize_options(struct noisechiselparams *p,
   cp->numthreads         = gal_threads_number();
   cp->coptions           = gal_commonopts_options;
 
+  /* Program specific initialization. */
+  p->snthresh=NAN;
 
   /* Modify common options. */
   for(i=0; !gal_options_is_last(&cp->coptions[i]); ++i)
@@ -234,6 +236,9 @@ ui_read_check_only_options(struct noisechiselparams *p)
   if(p->openingngb!=4 && p->openingngb!=8)
     error(EXIT_FAILURE, 0, "%zu not acceptable for `--openingngb'. It must "
           "be 4 or 8 (specifying the type of connectivity)", p->openingngb);
+  if(p->dopeningngb!=4 && p->dopeningngb!=8)
+    error(EXIT_FAILURE, 0, "%zu not acceptable for `--dopeningngb'. It must "
+          "be 4 or 8 (specifying the type of connectivity)", p->dopeningngb);
   if(p->holengb!=4 && p->holengb!=8)
     error(EXIT_FAILURE, 0, "%zu not acceptable for `--holengb'. It must "
           "be 4 or 8 (specifying the type of connectivity)", p->holengb);
diff --git a/bin/noisechisel/ui.h b/bin/noisechisel/ui.h
index 66f32cc..c42df7d 100644
--- a/bin/noisechisel/ui.h
+++ b/bin/noisechisel/ui.h
@@ -92,9 +92,12 @@ enum option_keys_enum
   UI_KEY_OPENINGNGB,
   UI_KEY_SKYFRACNOBLANK,
   UI_KEY_CHECKDETSKY,
+  UI_KEY_DOPENING,
+  UI_KEY_DOPENINGNGB,
   UI_KEY_HOLENGB,
   UI_KEY_PSEUDOCONCOMP,
   UI_KEY_CHECKSN,
+  UI_KEY_SNTHRESH,
   UI_KEY_DETGROWMAXHOLESIZE,
   UI_KEY_CLEANGROWNDET,
   UI_KEY_CHECKDETECTION,
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 202384e..7b18c00 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -16516,6 +16516,16 @@ the image. For more, see the description of this 
option in @ref{Detection
 options}.
 
 @item
address@hidden: Number of openings after applying
address@hidden For more, see the description of this option in
address@hidden options}.
+
address@hidden
address@hidden: Number of openings after applying
address@hidden For more, see the description of this option in
address@hidden options}.
+
address@hidden
 @option{--holengb}: The connectivity (defined by the number of neighbors)
 to fill holes after applying @option{--dthresh} (above) to find
 pseudo-detections. For more, see the description of this option in
@@ -16527,6 +16537,11 @@ neighbors) to find individual pseudo-detections. For 
more, see the
 description of this option in @ref{Detection options}.
 
 @item
address@hidden: Manually set the S/N of true pseudo-detections and
+thus avoid the need to manually identify this value. For more, see the
+description of this option in @ref{Detection options}.
+
address@hidden
 @option{--detgrowquant}: is used to grow the final true detections until a
 given quantile in the same way that clumps are grown during segmentation
 (compare columns 2 and 3 in Figure 10 of the paper). It replaces the old
@@ -17032,6 +17047,14 @@ ratio of the resulting `pseudo-detections' are used to 
identify true
 vs. false detections. See Section 3.1.5 and Figure 7 in Akhlaghi and
 Ichikawa (2015) for a very complete explanation.
 
address@hidden --dopening=INT
+The number of openings to do after applying @option{--dthresh}.
+
address@hidden --dopeningngb=INT
+The connectivity used in the opening of @option{--dopening}. In a 2D image
+this must be either 4 or 8. The stronger the connectivity, the more smaller
+regions will be discarded.
+
 @item --holengb=INT
 The connectivity (defined by the number of neighbors) to fill holes after
 applying @option{--dthresh} (above) to find pseudo-detections. For example
@@ -17092,6 +17115,17 @@ that this is only calculated for the large mesh grids 
that satisfy the
 minimum fraction of undetected pixels (value of @option{--minbfrac}) and
 minimum number of pseudo-detections (value of @option{--minnumfalse}).
 
address@hidden --snthresh=FLT
+Manually set the signal-to-noise ratio of true pseudo-detections. With this
+option, NoiseChisel will not attempt to find pseudo-detections over the
+noisy regions of the dataset, but will directly go onto applying the
+manually input value.
+
+This option is useful in crowded images where there is no blank sky to find
+the sky pseudo-detections. You can get this value on a similarly reduced
+dataset (from another region of the Sky with more undetected regions
+spaces).
+
 @item -d FLT
 @itemx --detgrowquant=FLT
 Quantile limit to ``grow'' the final detections. As discussed in the