[gnuastro-commits] master 280fee2: Threshod for no erosion added to NoiseChisel

[Mohammad Akhlaghi]

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

    Threshod for no erosion added to NoiseChisel
    
    Pure erosion is going to carve off sharp and small objects completely out
    of the image. Therefore to avoid missing such sharp and small objects
    (which have significant pixels, but not over a large area) NoiseChisel now
    has a new `--noerodequant' option. Through it, the level of significance
    can be specified. All pixels with a value larger than this significance
    level will not be eroded.
    
    In order to do this cleanly, the `gal_arraymanip_uchar_replace' was added
    to the arraymanip library and some corrections were made in the comments of
    the mesh library.
    
    This finishes task #14139.
---
 doc/gnuastro.texi                   |   15 ++++++++
 lib/arraymanip.c                    |   17 +++++++--
 lib/gnuastro/arraymanip.h           |    4 +++
 lib/gnuastro/mesh.h                 |    2 +-
 lib/mesh.c                          |   12 ++++---
 src/noisechisel/args.h              |   15 +++++++-
 src/noisechisel/astnoisechisel.conf |    1 +
 src/noisechisel/binary.h            |   17 +++++++--
 src/noisechisel/detection.c         |   11 ++++--
 src/noisechisel/main.h              |    2 ++
 src/noisechisel/thresh.c            |   68 +++++++++++++++++++++++------------
 src/noisechisel/ui.c                |   18 ++++++++++
 12 files changed, 145 insertions(+), 37 deletions(-)

diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 7192a7a..91fb63b 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -10707,6 +10707,21 @@ edge with it. The 8-connected neighbors on the other 
hand include the
 corner with this pixel. See Figure 6 (a) and (b) in Akhlaghi and
 Ichikawa (2015) for a demonstration.
 
address@hidden --noerodequant
+Pure erosion is going to carve off sharp and small objects completely out
+of the detected regions. This option can be used to avoid missing such
+sharp and small objects (which have significant pixels, but not over a
+large area). All pixels with a value larger than the significance level
+specified by this option will not be eroded during the erosion step
+above. However, they will undergo the erosion and dilation of the opening
+step below.
+
+Like the @option{--qthresh} option, the significance level is
+determined using the quantile (a value between 0 and 1). Just as a
+reminder, in the normal distribution, @mymath{1\sigma}, @mymath{1.5\sigma},
+and @mymath{2\sigma} are approximately on the 0.84, 0.93, and 0.98
+quantiles.
+
 @item -p
 @itemx --opening
 (@option{=INT}) Depth of opening to be applied to the eroded binary
diff --git a/lib/arraymanip.c b/lib/arraymanip.c
index 419f880..5e8583e 100644
--- a/lib/arraymanip.c
+++ b/lib/arraymanip.c
@@ -210,8 +210,7 @@ gal_arraymanip_freplace_value(float *in, size_t size, float 
from, float to)
 void
 gal_arraymanip_freplace_nonnans(float *in, size_t size, float to)
 {
-  float *fpt;
-  fpt=in+size;
+  float *fpt=in+size;
   do
     *in = isnan(*in) ? *in : to;
   while(++in<fpt);
@@ -262,6 +261,20 @@ gal_arraymanip_no_nans_double(double *in, size_t *size)
 
 
 
+void
+gal_arraymanip_uchar_replace(unsigned char *in, size_t size,
+                             unsigned char from, unsigned char to)
+{
+  unsigned char *fpt=in+size;
+  do
+    *in = *in==from ? to : *in;
+  while(++in<fpt);
+}
+
+
+
+
+
 
 
 
diff --git a/lib/gnuastro/arraymanip.h b/lib/gnuastro/arraymanip.h
index 69765f4..aa33cfc 100644
--- a/lib/gnuastro/arraymanip.h
+++ b/lib/gnuastro/arraymanip.h
@@ -60,6 +60,10 @@ void
 gal_arraymanip_no_nans_double(double *in, size_t *size);
 
 void
+gal_arraymanip_uchar_replace(unsigned char *in, size_t size,
+                             unsigned char from, unsigned char to);
+
+void
 gal_arraymanip_fmultip_const(float *in, size_t size, float a);
 
 void
diff --git a/lib/gnuastro/mesh.h b/lib/gnuastro/mesh.h
index ec4e9bb..c533d39 100644
--- a/lib/gnuastro/mesh.h
+++ b/lib/gnuastro/mesh.h
@@ -72,7 +72,7 @@ struct gal_mesh_thread_params
 
    In short, the meshs in each channel have to be contiguous to
    facilitate the neighbor analysis in interpolation and other channel
-   specific jobs.g
+   specific jobs.
 
    The operations on the meshs might need more than one output, for
    example the mean and the standard deviation. So we have two garrays
diff --git a/lib/mesh.c b/lib/mesh.c
index 59ff44d..1aaa3a8 100644
--- a/lib/mesh.c
+++ b/lib/mesh.c
@@ -860,10 +860,10 @@ gal_mesh_free_mesh(struct gal_mesh_params *mp)
    1. A pointer to the gal_mesh_params structure that keeps all the
       information.
 
-   2. A pointer to a function that returns and gets a `void *' as its
-      only argument. This function will be directly given to
-      pthread_create. Through this function, you can any function that
-      you wish to operate on the mesh grid with.
+   2. A pointer to a function that returns and gets a `void *' as its only
+      argument. This function will be directly given to
+      pthread_create. Through this argument, you can choose the function to
+      operate on the mesh grid.
 
    3. The size of each element to copy the mesh grid into, this has to
       be type size of the same type that constitutes `img' in
@@ -873,7 +873,9 @@ gal_mesh_free_mesh(struct gal_mesh_params *mp)
       example sort them) in each mesh.
 
    4. If the value of this argument is 1, then a second garray will be
-      allocated in case your operation needs one.
+      allocated in case your operation needs it.
+
+   5. Wether the allocated garrays should be initialized or not.
 */
 void
 gal_mesh_operate_on_mesh(struct gal_mesh_params *mp,
diff --git a/src/noisechisel/args.h b/src/noisechisel/args.h
index 832f8ca..c016078 100644
--- a/src/noisechisel/args.h
+++ b/src/noisechisel/args.h
@@ -75,7 +75,7 @@ const char doc[] =
    f j w x z
    A J K W X Y Z
 
-   Used numbers: <=518 (except 503, and 515)
+   Used numbers: <=518 (except 515)
 
    Options with keys (second structure element) larger than 500 do not
    have a short version.
@@ -328,6 +328,14 @@ static struct argp_option options[] =
       4
     },
     {
+      "noerodequant",
+      503,
+      "FLT",
+      0,
+      "Quantile for no erosion.",
+      4
+    },
+    {
       "opening",
       'p',
       "INT",
@@ -681,6 +689,11 @@ parse_opt(int key, char *arg, struct argp_state *state)
                               NULL, 0);
       p->up.erodengbset=1;
       break;
+    case 503:
+      gal_checkset_float_l_0_s_1(arg, &p->qthresh, "noerodequant", key,
+                                 SPACK, NULL, 0);
+      p->up.noerodequantset=1;
+      break;
     case 'p':
       gal_checkset_sizet_el_zero(arg, &p->opening, "opening", key, SPACK,
                                  NULL, 0);
diff --git a/src/noisechisel/astnoisechisel.conf 
b/src/noisechisel/astnoisechisel.conf
index 32869be..4836354 100644
--- a/src/noisechisel/astnoisechisel.conf
+++ b/src/noisechisel/astnoisechisel.conf
@@ -45,6 +45,7 @@
  qthresh          0.3
  erode              2
  erodengb           4
+ noerodequant  0.9331
  opening            1
  openingngb         8
  sigclipmultip      3
diff --git a/src/noisechisel/binary.h b/src/noisechisel/binary.h
index a6d42bd..21db117 100644
--- a/src/noisechisel/binary.h
+++ b/src/noisechisel/binary.h
@@ -34,8 +34,21 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 
 
-/* Special values: */
-#define BINARYTMP 2
+/* Special values:
+
+     BINARYNOOP:   Value that no binary operation should be preformed on.
+     BINARYTMP:    Temporary value to use within one function.
+
+   Note that through the `setbytblank' function, practically we are using
+   the value `GAL_FITS_BYTE_BLANK' (from `gnuastro/fits.h') for blank
+   binary values. Recall that due to the nature of the CPU (which operates
+   on 8-bits), in practice it is much more efficient to work on a byte (or
+   8-bits) rather than each bit. So in practice we can use 256 values for
+   meta-data analyais (like blank values, or temporary values at etc),
+   eventhough the main values we are working with are 0 and 1.
+*/
+#define BINARYNOOP  2
+#define BINARYTMP   3
 
 
 
diff --git a/src/noisechisel/detection.c b/src/noisechisel/detection.c
index 891a021..2fedd47 100644
--- a/src/noisechisel/detection.c
+++ b/src/noisechisel/detection.c
@@ -66,14 +66,13 @@ void
 initialdetection(struct noisechiselparams *p)
 {
   int verb=p->cp.verb;
-  char report[GAL_TIMING_VERB_MSG_LENGTHS_2_V];
   size_t i, s0=p->smp.s0, s1=p->smp.s1;
   char *detectionname=p->detectionname;
+  char report[GAL_TIMING_VERB_MSG_LENGTHS_2_V];
 
   /* Find the threshold and apply it, if there are blank pixels in the
      image, set those pixels to the binary blank value too: */
   findapplyqthreshold(p);
-  if(p->anyblank)  setbytblank(p->img, p->byt, s0*s1);
   if(detectionname)
     gal_fits_array_to_file(detectionname, "Thresholded", BYTE_IMG,
                            p->byt, s0, s1, p->anyblank, p->wcs, NULL,
@@ -87,13 +86,17 @@ initialdetection(struct noisechiselparams *p)
 
 
 
-  /* Erode the thresholded image: */
+
+  /* Erode the thresholded image, then convert the pixels that should not
+     have been eroded back to 1 so the next step (opening) can be
+     completely ignorant of them. */
   if(p->erodengb==4)
     for(i=0;i<p->erode;++i)
       dilate0_erode1_4con(p->byt, s0, s1, 1);
   else
     for(i=0;i<p->erode;++i)
       dilate0_erode1_8con(p->byt, s0, s1, 1);
+  gal_arraymanip_uchar_replace(p->byt, s0*s1, BINARYNOOP, 1);
   if(detectionname)
     gal_fits_array_to_file(detectionname, "Eroded", BYTE_IMG, p->byt,
                            s0, s1, p->anyblank, p->wcs, NULL,
@@ -107,6 +110,7 @@ initialdetection(struct noisechiselparams *p)
 
 
 
+
   /* Do the opening: */
   opening(p->byt, s0, s1, p->opening, p->openingngb);
   if(detectionname)
@@ -122,6 +126,7 @@ initialdetection(struct noisechiselparams *p)
 
 
 
+
   /* Label the connected regions. Note that p->olab was allocated in
      ui.c and will be freed there. */
   p->numobjects=BF_concmp(p->byt, p->olab, s0, s1, p->anyblank, 4);
diff --git a/src/noisechisel/main.h b/src/noisechisel/main.h
index 2043a9d..3b58877 100644
--- a/src/noisechisel/main.h
+++ b/src/noisechisel/main.h
@@ -66,6 +66,7 @@ struct uiparams
   int           qthreshset;
   int             erodeset;
   int          erodengbset;
+  int      noerodequantset;
   int           openingset;
   int        openingngbset;
   int          minbfracset;
@@ -122,6 +123,7 @@ struct noisechiselparams
   /* Detection: */
   float              *conv; /* Convolved image.                            */
   float            qthresh; /* Quantile threshold on convolved img.        */
+  float       noerodequant; /* Quantile for no erosion.                    */
   size_t             erode; /* Number of times to erode thresholded image. */
   int             erodengb; /* Use 4 or 8 connectivity in erosion.         */
   size_t           opening; /* Depth of opening to apply to eroded image.  */
diff --git a/src/noisechisel/thresh.c b/src/noisechisel/thresh.c
index 950b78c..48105ed 100644
--- a/src/noisechisel/thresh.c
+++ b/src/noisechisel/thresh.c
@@ -35,6 +35,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/spatialconvolve.h>
 
 #include "main.h"
+#include "binary.h"
 #include "noisechisel.h"
 
 
@@ -59,7 +60,6 @@ qthreshonmesh(void *inparam)
   float *mponeforall=mp->oneforall;
   float *oneforall=&mponeforall[mtp->id*mp->maxs0*mp->maxs1];
 
-  float qthresh=p->qthresh;
   size_t modeindex=(size_t)(-1);
   size_t s0, s1, ind, row, num, start, is1=mp->s1;
   size_t i, *indexs=&mp->indexs[mtp->id*mp->thrdcols];
@@ -97,12 +97,18 @@ qthreshonmesh(void *inparam)
       /* Do the desired operation on the mesh: */
       if(num)
         {
-          qsort(oneforall, num, sizeof *oneforall, gal_qsort_float_increasing);
-          gal_mode_index_in_sorted(oneforall, num, mirrordist, &modeindex,
-                                   &modesym);
-          if( modesym>GAL_MODE_SYM_GOOD && 
(float)modeindex/(float)num>minmodeq)
-            mp->garray1[ind]=
-              oneforall[gal_statistics_index_from_quantile(num, qthresh)];
+          qsort(oneforall, num, sizeof *oneforall,
+                gal_qsort_float_increasing);
+          gal_mode_index_in_sorted(oneforall, num, mirrordist,
+                                   &modeindex, &modesym);
+          if( modesym>GAL_MODE_SYM_GOOD
+              && (float)modeindex/(float)num>minmodeq)
+            {
+              mp->garray1[ind] = oneforall[
+                 gal_statistics_index_from_quantile(num, p->qthresh) ];
+              mp->garray2[ind] = oneforall[
+                 gal_statistics_index_from_quantile(num, p->noerodequant) ];
+            }
         }
     }
 
@@ -127,17 +133,18 @@ qthreshonmesh(void *inparam)
 void
 applythreshold(struct noisechiselparams *p)
 {
-  struct gal_mesh_params *mp=&p->smp; /* `mp' instead of `smp' so you can try 
*/
-                                      /* with p->lmp if you like.             
*/
+  struct gal_mesh_params *mp=&p->smp; /* `mp' instead of `smp' so you     */
+                                      /* can try with p->lmp if you like. */
+  float qthreshv, nethreshv;
   unsigned char *b, *byt=p->byt;
   size_t gs0=mp->gs0, gs1=mp->gs1, nch1=mp->nch1;
   size_t s0, s1, fs1=mp->gs1*mp->nch1, chid, inchid;
+  float *f, *fp, *garray1=mp->garray1, *array=p->conv;
   size_t *types=mp->types, *ts0=mp->ts0, *ts1=mp->ts1;
   size_t i, f0, f1, nmeshi=mp->nmeshi, nmeshc=mp->nmeshc;
   size_t row, startind, *start=mp->start, meshid, is1=mp->s1;
-  float *f, *fp, thresh, *garray1=mp->garray1, *array=p->conv;
 
-    /* Fill the array: */
+  /* Fill the array: */
   for(i=0;i<nmeshi;++i)
     {
       /* Set the proper meshid. */
@@ -153,19 +160,32 @@ applythreshold(struct noisechiselparams *p)
 
       /* Fill the output array with the value in this mesh: */
       row=0;
-      thresh=garray1[i];
+      qthreshv=garray1[i];
       s0=ts0[types[meshid]];
       s1=ts1[types[meshid]];
       startind=start[meshid];
+      nethreshv=mp->garray2[i];
       do
         {
-          b = byt + startind + row * is1;
-          fp= ( f = array + startind + row * is1 ) + s1;
-          do *b++ = isnan(*f) || *f>thresh ? 1 : 0; while(++f<fp);
+          b  = byt + startind + row * is1;
+          fp = ( f = array + startind + row * is1 ) + s1;
+          do
+            /* The image might have NaN pixels too, in which case any
+               comparison will fail and the value in `b' will be 0, This is
+               no problem, after this loop, all NaN pixels will be set to a
+               fixed blank value. */
+            *b++ = ( *f>qthreshv
+                     ? ( *f>nethreshv ? BINARYNOOP : 1 )
+                     : 0 );
+          while(++f<fp);
           ++row;
         }
       while(row<s0);
     }
+
+  /* If there are any NaN pixels (in the float image), set them to the
+     blank value for the unsigned character type. */
+  if(p->anyblank) setbytblank(p->img, p->byt, p->smp.s0*p->smp.s1);
 }
 
 
@@ -177,21 +197,23 @@ findapplyqthreshold(struct noisechiselparams *p)
 {
   struct gal_mesh_params *mp=&p->smp;
 
-  /* Find the threshold on each mesh: */
-  gal_mesh_operate_on_mesh(mp, qthreshonmesh, sizeof(float), 0, 1);
+  /* Find the threshold on each mesh that has the proper conditions. */
+  gal_mesh_operate_on_mesh(mp, qthreshonmesh, sizeof(float), 1, 1);
   if(p->threshname)
-    gal_mesh_value_file(mp, p->threshname, "Quantile values", NULL,
-                        p->wcs, SPACK_STRING);
+    gal_mesh_value_file(mp, p->threshname, "Quantile threshold value",
+                        "No erode threshold value", p->wcs, SPACK_STRING);
 
+  /* Interpolate to have a value for each mesh. */
   gal_mesh_interpolate(mp, "Interpolating quantile threshold");
   if(p->threshname)
-    gal_mesh_value_file(mp, p->threshname, "Interpolated", NULL,
-                        p->wcs, SPACK_STRING);
+    gal_mesh_value_file(mp, p->threshname, "qthreshv interpolated",
+                        "nethreshv interpolated", p->wcs, SPACK_STRING);
 
+  /* Smooth the interpolation  */
   gal_mesh_smooth(mp);
   if(p->threshname)
-    gal_mesh_value_file(mp, p->threshname, "smoothed", NULL,
-                        p->wcs, SPACK_STRING);
+    gal_mesh_value_file(mp, p->threshname, "qthreshv smoothed",
+                        "nethreshv smoothed", p->wcs, SPACK_STRING);
 
   /* Apply the threshold on all the pixels: */
   applythreshold(p);
diff --git a/src/noisechisel/ui.c b/src/noisechisel/ui.c
index 3df6e97..47793a7 100644
--- a/src/noisechisel/ui.c
+++ b/src/noisechisel/ui.c
@@ -264,6 +264,13 @@ readconfig(char *filename, struct noisechiselparams *p)
                                   SPACK, filename, lineno);
           up->erodengbset=1;
         }
+      else if(strcmp(name, "noerodequant")==0)
+        {
+          if(up->noerodequantset) continue;
+          gal_checkset_float_l_0_s_1(value, &p->noerodequant, name, key,
+                                     SPACK, filename, lineno);
+          up->noerodequantset=1;
+        }
       else if(strcmp(name, "opening")==0)
         {
           if(up->openingset) continue;
@@ -470,6 +477,8 @@ printvalues(FILE *fp, struct noisechiselparams *p)
     fprintf(fp, CONF_SHOWFMT"%lu\n", "erode", p->erode);
   if(up->erodengbset)
     fprintf(fp, CONF_SHOWFMT"%d\n", "erodengb", p->erodengb);
+  if(up->noerodequantset)
+    fprintf(fp, CONF_SHOWFMT"%.3f\n", "noerodequant", p->noerodequant);
   if(up->openingset)
     fprintf(fp, CONF_SHOWFMT"%lu\n", "opening", p->opening);
   if(up->openingngbset)
@@ -577,6 +586,8 @@ checkifset(struct noisechiselparams *p)
     GAL_CONFIGFILES_REPORT_NOTSET("erode");
   if(up->erodengbset==0)
     GAL_CONFIGFILES_REPORT_NOTSET("erodengb");
+  if(up->noerodequantset==0)
+    GAL_CONFIGFILES_REPORT_NOTSET("noerodequant");
   if(up->openingset==0)
     GAL_CONFIGFILES_REPORT_NOTSET("opening");
   if(up->openingngbset==0)
@@ -645,6 +656,13 @@ sanitycheck(struct noisechiselparams *p)
   /* Make sure the input file exists. */
   gal_checkset_check_file(p->up.inputname);
 
+  /* Make sure that the noerode quantile is larger than qthresh. */
+  if( p->noerodequant <= p->qthresh)
+    error(EXIT_FAILURE, 0, "The quantile for no erosion (`--noerodequant') "
+          "must be larger than the base quantile threshold (`--qthresh', "
+          "or `-t'). You have provided %.4f and %.4f for the former and "
+          "latter, respectively.", p->noerodequant, p->qthresh);
+
   /* Set the maskname and mask hdu accordingly: */
   gal_fits_file_or_ext_name(p->up.inputname, p->cp.hdu, p->up.masknameset,
                                  &p->up.maskname, p->up.mhdu, p->up.mhduset,