[gnuastro-commits] master d1ce841: NoiseChisel can take image instead of doing convolution

[Mohammad Akhlaghi]

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

    NoiseChisel can take image instead of doing convolution
    
    When testing NoiseChisel's parameters for different datasets/analysis, it
    is really troubling that usually, the single most time (and CPU) consuming
    step is convolution. NoiseChisel now has a new `--convolved' option so the
    user can specify the convolved image, thus telling NoiseChisel to read that
    image instead of doing convolution.
---
 NEWS                                |  6 ++++
 bin/noisechisel/args.h              | 26 +++++++++++++++
 bin/noisechisel/astnoisechisel.conf | 11 ++++---
 bin/noisechisel/main.h              |  2 ++
 bin/noisechisel/noisechisel.c       | 36 ++++++++++++--------
 bin/noisechisel/ui.c                | 48 +++++++++++++++++++++------
 bin/noisechisel/ui.h                |  2 ++
 doc/gnuastro.texi                   | 66 +++++++++++++++++++++++++++++++++----
 8 files changed, 162 insertions(+), 35 deletions(-)

diff --git a/NEWS b/NEWS
index d4a8d20..b80a32c 100644
--- a/NEWS
+++ b/NEWS
@@ -50,6 +50,12 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
   properties of the dataset vary gradually and sampling from the whole
   dataset might produce biased results.
 
+  NoiseChisel: with the new `--convolved' option, NoiseChisel will not
+  convolve the input any more and use the given image. In many cases,
+  convolution is the most time consuming step of NoiseChisel, so with this
+  option, you can greatly speed up your tests (to find the best parameters
+  for a given analysis). See the book for more information and examples.
+
   NoiseChisel: with the new `--widekernel' option it is now possible to use
   a wider kernel to identify which tiles contain signal. The rest of the
   steps (identifying the quantile threshold on the selected tiles and etc)
diff --git a/bin/noisechisel/args.h b/bin/noisechisel/args.h
index a6a68b2..e44efbd 100644
--- a/bin/noisechisel/args.h
+++ b/bin/noisechisel/args.h
@@ -60,6 +60,32 @@ struct argp_option program_options[] =
       GAL_OPTIONS_NOT_SET
     },
     {
+      "convolved",
+      UI_KEY_CONVOLVED,
+      "STR",
+      0,
+      "Convolved image file to avoid convolution.",
+      GAL_OPTIONS_GROUP_INPUT,
+      &p->convolvedname,
+      GAL_TYPE_STRING,
+      GAL_OPTIONS_RANGE_ANY,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
+    },
+    {
+      "convolvedhdu",
+      UI_KEY_CONVOLVEDHDU,
+      "STR",
+      0,
+      "HDU/extension of convolved image in file.",
+      GAL_OPTIONS_GROUP_INPUT,
+      &p->convolvedhdu,
+      GAL_TYPE_STRING,
+      GAL_OPTIONS_RANGE_ANY,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
+    },
+    {
       "widekernel",
       UI_KEY_WIDEKERNEL,
       "STR",
diff --git a/bin/noisechisel/astnoisechisel.conf 
b/bin/noisechisel/astnoisechisel.conf
index bf44ac6..71bdd47 100644
--- a/bin/noisechisel/astnoisechisel.conf
+++ b/bin/noisechisel/astnoisechisel.conf
@@ -18,13 +18,14 @@
 # warranty.
 
 # Input:
- khdu                1
- wkhdu               1
- minskyfrac        0.7
- minnumfalse       100
+ khdu                 1
+ wkhdu                1
+ convolvedhdu         1
+ minskyfrac         0.7
+ minnumfalse        100
 
 # Tessellation
- largetilesize 200,200
+ largetilesize  200,200
 
 # Detection:
  mirrordist         1.5
diff --git a/bin/noisechisel/main.h b/bin/noisechisel/main.h
index c994c4f..2f17f2f 100644
--- a/bin/noisechisel/main.h
+++ b/bin/noisechisel/main.h
@@ -48,6 +48,8 @@ struct noisechiselparams
   char             *inputname;  /* Input filename.                        */
   char            *kernelname;  /* Input kernel filename.                 */
   char        *widekernelname;  /* Name of wider kernel to be used.       */
+  char         *convolvedname;  /* Convolved image (to avoid convolution).*/
+  char          *convolvedhdu;  /* HDU of convolved image.                */
   char                  *khdu;  /* Kernel HDU.                            */
   char                 *wkhdu;  /* Wide kernel HDU.                       */
   uint8_t       skysubtracted;  /* Input has been Sky subtracted before.  */
diff --git a/bin/noisechisel/noisechisel.c b/bin/noisechisel/noisechisel.c
index af2723c..f417454 100644
--- a/bin/noisechisel/noisechisel.c
+++ b/bin/noisechisel/noisechisel.c
@@ -59,20 +59,31 @@ noisechisel_convolve(struct noisechiselparams *p)
   struct gal_tile_two_layer_params *tl=&p->cp.tl;
 
   /* Convovle with sharper kernel. */
-  if(!p->cp.quiet) gettimeofday(&t1, NULL);
-  p->conv=gal_convolve_spatial(tl->tiles, p->kernel, p->cp.numthreads,
-                               1, tl->workoverch);
-  gal_checkset_allocate_copy(p->widekernel?"CONVOLVED-SHARPER":"CONVOLVED",
-                             &p->conv->name);
-
-  /* Report and write check images if necessary. */
-  if(!p->cp.quiet)
+  if(p->conv==NULL)
     {
-      if(p->widekernel)
-        gal_timing_report(&t1, "Convolved with sharper kernel.", 1);
-      else
-        gal_timing_report(&t1, "Convolved with given kernel.", 1);
+      /* Make the convolved image. */
+      if(!p->cp.quiet) gettimeofday(&t1, NULL);
+      p->conv = gal_convolve_spatial(tl->tiles, p->kernel, p->cp.numthreads,
+                                     1, tl->workoverch);
+
+      /* Report and write check images if necessary. */
+      if(!p->cp.quiet)
+        {
+          if(p->widekernel)
+            gal_timing_report(&t1, "Convolved with sharper kernel.", 1);
+          else
+            gal_timing_report(&t1, "Convolved with given kernel.", 1);
+        }
     }
+
+  /* Set a fixed name for the convolved image (since it will be used in
+     many check images). */
+  if(p->conv->name) free(p->conv->name);
+  gal_checkset_allocate_copy( ( p->widekernel
+                                ? "CONVOLVED-SHARPER"
+                                : "CONVOLVED" ), &p->conv->name);
+
+  /* Save the convolution step if necessary. */
   if(p->detectionname)
     {
       gal_fits_img_write(p->input, p->detectionname, NULL, PROGRAM_NAME);
@@ -87,7 +98,6 @@ noisechisel_convolve(struct noisechiselparams *p)
                                     p->cp.numthreads, 1, tl->workoverch);
       gal_checkset_allocate_copy("CONVOLVED-WIDER", &p->wconv->name);
 
-      /* Report the status: */
       if(!p->cp.quiet)
         gal_timing_report(&t1, "Convolved with wider kernel.", 1);
     }
diff --git a/bin/noisechisel/ui.c b/bin/noisechisel/ui.c
index 78ffd59..fb9c423 100644
--- a/bin/noisechisel/ui.c
+++ b/bin/noisechisel/ui.c
@@ -214,6 +214,13 @@ parse_opt(int key, char *arg, struct argp_state *state)
 static void
 ui_read_check_only_options(struct noisechiselparams *p)
 {
+  /* If the convolved option is given, then the convolved HDU is also
+     mandatory. */
+  if(p->convolvedname && p->convolvedhdu==NULL)
+    error(EXIT_FAILURE, 0, "no value given to `--convolvedhdu'. When the "
+          "`--convolved' option is called to specify a convolved image and "
+          "avoid convolution, it is mandatory to also specify a HDU for it");
+
   /* Make sure the connectivities have the correct values. */
   if(p->erodengb!=4 && p->erodengb!=8)
     error(EXIT_FAILURE, 0, "%zu not acceptable for `--erodengb'. It must "
@@ -587,12 +594,27 @@ ui_preparations(struct noisechiselparams *p)
           p->input->ndim);
 
 
-  /* Check for blank values to help later processing.  */
-  gal_blank_present(p->input, 1);
+  /* If a convolved image was given, read it in. Otherwise, read the given
+     kernel. */
+  if(p->convolvedname)
+    {
+      /* Read the input convolved image. */
+      p->conv = gal_fits_img_read_to_type(p->convolvedname, p->convolvedhdu,
+                                          GAL_TYPE_FLOAT32, p->cp.minmapsize,
+                                          0, 0);
+
+      /* Make sure the convolved image is the same size as the input. */
+      if( gal_data_dsize_is_different(p->input, p->conv) )
+        error(EXIT_FAILURE, 0, "%s (hdu %s) which was given to `--convolved' "
+              "is not the same size as the input: %s (hdu: %s)",
+              p->convolvedname, p->convolvedhdu, p->inputname, p->cp.hdu);
+    }
+  else
+    ui_prepare_kernel(p);
 
 
-  /* Read in the kernel for convolution. */
-  ui_prepare_kernel(p);
+  /* Check for blank values to help later processing.  */
+  gal_blank_present(p->input, 1);
 
 
   /* Prepare the tessellation. */
@@ -690,13 +712,19 @@ ui_read_check_inputs_setup(int argc, char *argv[],
       printf("  - Using %zu CPU thread%s\n", p->cp.numthreads,
              p->cp.numthreads==1 ? "." : "s.");
       printf("  - Input: %s (hdu: %s)\n", p->inputname, p->cp.hdu);
-      if(p->kernelname)
-        printf("  - %s: %s (hdu: %s)\n",
-               p->widekernelname ? "Sharp Kernel" : "Kernel",
-               p->kernelname, p->khdu);
+      if(p->convolvedname)
+        printf("  - Convolved input: %s (hdu: %s)\n",
+               p->convolvedname, p->convolvedhdu);
       else
-        printf("  - %s: FWHM=2 pixel Gaussian.\n",
-               p->widekernelname ? "Sharp Kernel" : "Kernel");
+        {
+          if(p->kernelname)
+            printf("  - %s: %s (hdu: %s)\n",
+                   p->widekernelname ? "Sharp Kernel" : "Kernel",
+                   p->kernelname, p->khdu);
+          else
+            printf("  - %s: FWHM=2 pixel Gaussian.\n",
+                   p->widekernelname ? "Sharp Kernel" : "Kernel");
+        }
       if(p->widekernelname)
         printf("  - Wide Kernel: %s (hdu: %s)\n", p->widekernelname,
                p->wkhdu);
diff --git a/bin/noisechisel/ui.h b/bin/noisechisel/ui.h
index 9481981..ec3fe82 100644
--- a/bin/noisechisel/ui.h
+++ b/bin/noisechisel/ui.h
@@ -76,6 +76,8 @@ enum option_keys_enum
   /* Only with long version (start with a value 1000, the rest will be set
      automatically). */
   UI_KEY_KHDU               = 1000,
+  UI_KEY_CONVOLVED,
+  UI_KEY_CONVOLVEDHDU,
   UI_KEY_WKHDU,
   UI_KEY_MINNUMFALSE,
   UI_KEY_ONLYDETECTION,
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 8d073be..34304f7 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -12527,6 +12527,58 @@ Akhlaghi and Ichikawa [2015].
 HDU containing the kernel in the file given to the @option{--kernel}
 option.
 
address@hidden --convolvedname=STR
+Use this file as the convolved image and don't do convolution and ignore
address@hidden NoiseChisel will just check the size of the given
+dataset is the same as the input's size. If a wrong image (with the same
+size) is given to this option, the results (errors, bugs, and etc) are
+un-predictable. So please use this option with care and in a highly
+controlled environment. On such scenario is discussed below.
+
+In almost all situations, as the input gets larger, the single most CPU and
+time consuming step in NoiseChisel is convolution (the first step in its
+processing). The process of testing NoiseChisel for the best parameters in
+a given analysis will involve running NoiseChisel multiple times: to see
+the effect of each new option value. Therefore, once the kernel is
+finalized, re-convolving the input will greatly hinder fast testing of
+higher-level parameters. With this option, you can convolve the input image
+with your chosen kernel once before running NoiseChisel, then feed it to
+NoiseChisel on each test run and thus save valuable time for better/more
+tests.
+
+To build your desired convolution kernel, you can use
address@hidden To convolve the image with a given kernel you can use
address@hidden Spatial domain convolution is mandatory. In the frequency
+domain, blank pixels (if present) will cover the whole image and gradients
+will appear on the edges, see @ref{Spatial vs. Frequency domain}.
+
+Below you can see an example of such a scenario: you want to see how
+variation of the growth level (through the @option{--detgrowquant} option)
+will affect the final result. Recall that you can ignore all the extra
+spaces, new lines, and address@hidden' if you are typing in the terminal (in a
+shell script, remove the @code{$} signs at the start of the lines).
+
address@hidden
+## Make the kernel to convolve with.
+$ astmkprof --oversample=1 --kernel=gaussian,2,5
+
+## Convolve the input with the given kernel.
+$ astconvolve input.fits --kernel=kernel.fits                \
+              --domain=spatial --output=convolved.fits
+
+## Run NoiseChisel with seven growth quantile values.
+$ for g in 60 65 70 75 80 85 90; do                          \
+    astnoisechisel input.fits --convolvedname=convolved.fits \
+                   --detgrowquant=0.$g --output=$g.fits;     \
+  done
address@hidden example
+
+
+
address@hidden --convolvedhdu=STR
+The HDU/extension containing the convolved image in the file given to
address@hidden
+
 @item -w STR
 @itemx --widekernel=STR
 File name of a wider kernel to use in estimating the difference of the mode
@@ -12631,13 +12683,6 @@ catalog). You can then run MakeCatalog on A normally, 
see
 clump labels images to those that NoiseChisel produced for A, see
 @ref{Invoking astmkcatalog}.
 
address@hidden --grownclumps
-In the output (see @ref{NoiseChisel output}) store the grown clumps (or
-full detected region if only one clump was present in that detection). By
-default the original clumps are stored as the third extension of the
-output, if this option is called, it is replaced with the grown clump
-labels.
-
 @item --continueaftercheck
 Continue NoiseChisel after any of the options starting with
 @option{--check}. NoiseChisel involves many steps and as a result, there
@@ -13073,6 +13118,13 @@ NoiseChisel. The `true' detections and clumps can be 
objectively
 identified from the noise characteristics of the image, so you don't
 have to give any hand input Signal to noise ratio.
 
address@hidden --grownclumps
+In the output (see @ref{NoiseChisel output}) store the grown clumps (or
+full detected region if only one clump was present in that detection). By
+default the original clumps are stored as the third extension of the
+output, if this option is called, it is replaced with the grown clump
+labels.
+
 @item --checksegmentation
 A file with the suffix @file{_seg.fits} will be created. This file keeps
 all the relevant steps in finding true clumps and segmenting the detections