[gnuastro-commits] master 0bc96c9 122/125: Clumps and their S/N tables over detections

[Mohammad Akhlaghi]

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

    Clumps and their S/N tables over detections
    
    The new implementation of NoiseChisel can now identify all the clumps over
    each detection and also save the S/N table into a file similar to the
    clumps over Sky regions in the new version as described in the previous
    commit. So the users can now actually check specific clump S/N values.
    
    Some other changes with this commit:
    
      - Some parts of the Installation chapter of the book ("CFITSIO" and
        "Installation directory") were edited to be more clear.
    
      - When the input size to `gal_data_initialize' is zero, the allocated
        array pointer will be `NULL'. Previously it would use `mmap', `malloc',
        or `calloc' to allocate a zero-sized array.
    
      - After freeing the allocated spaces, `gal_data_free_contents' will set
        the pointers to `NULL' to avoid later usages (for when the `gal_data_t'
        isn't freed immediately after calling `gal_data_free_contents'.
---
 bin/noisechisel/clumps.c       | 155 +++++++++++++---------
 bin/noisechisel/clumps.h       |  45 ++++++-
 bin/noisechisel/detection.c    |  29 +++--
 bin/noisechisel/main.h         |  12 +-
 bin/noisechisel/segmentation.c | 290 ++++++++++++++++++++++++++++++++++++++++-
 doc/gnuastro.texi              | 140 +++++++++++---------
 lib/data.c                     |  52 +++++---
 7 files changed, 558 insertions(+), 165 deletions(-)

diff --git a/bin/noisechisel/clumps.c b/bin/noisechisel/clumps.c
index ac831d8..688196c 100644
--- a/bin/noisechisel/clumps.c
+++ b/bin/noisechisel/clumps.c
@@ -42,6 +42,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include "main.h"
 
 #include "ui.h"
+#include "clumps.h"
 #include "threshold.h"
 
 
@@ -54,32 +55,8 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 
 /**********************************************************************/
-/*****************     Basic structures and macros    *****************/
+/*****************         Basic macro values         *****************/
 /**********************************************************************/
-/* Parameters for all threads. */
-struct clumps_params
-{
-  int                     step; /* Counter if we want to check steps.      */
-  int                sky0_det1; /* If working on the Sky or Detections.    */
-  gal_data_t               *sn; /* Array of clump S/N tables.              */
-  gal_data_t            *snind; /* Array of clump S/N index (for check).   */
-  struct noisechiselparams  *p; /* Pointer to main NoiseChisel parameters. */
-  pthread_mutex_t     labmutex; /* Mutex to change the total numbers.      */
-};
-
-/* Parameters for one thread. */
-struct clumps_thread_params
-{
-  size_t                    id; /* ID of this detection/tile over tile.    */
-  size_t              *topinds; /* Indexs of all local maxima.             */
-  size_t            numinitial; /* The number of clumps found in this run. */
-  gal_data_t           *indexs; /* Array containing indexs of this object. */
-  gal_data_t             *info; /* Information for all clumps.             */
-  gal_data_t               *sn; /* Signal-to-noise ratio for these clumps. */
-  gal_data_t            *snind; /* Index of S/N for these clumps.          */
-  struct clumps_params  *clprm; /* Pointer to main structure.              */
-};
-
 /* Constants for the clump over-segmentation. */
 #define CLUMPS_RIVER     UINT32_MAX-2
 #define CLUMPS_TMPCHECK  UINT32_MAX-3
@@ -536,7 +513,7 @@ clumps_get_raw_info(struct clumps_thread_params *cltprm)
 
 
 /* Make an S/N table for the clumps in a given region. */
-static void
+void
 clumps_make_sn_table(struct clumps_thread_params *cltprm)
 {
   struct noisechiselparams *p=cltprm->clprm->p;
@@ -557,7 +534,7 @@ clumps_make_sn_table(struct clumps_thread_params *cltprm)
   cltprm->sn->dsize = gal_data_malloc_array(GAL_TYPE_SIZE_T, 1);
   cltprm->sn->array = gal_data_malloc_array(cltprm->sn->type, tablen);
   cltprm->sn->size  = cltprm->sn->dsize[0] = tablen;       /* After dsize. */
-  if(p->checksegmentation || p->checkclumpsn)
+  if( cltprm->clprm->snind )
     {
       cltprm->snind        = &cltprm->clprm->snind [ cltprm->id ];
       cltprm->snind->ndim  = 1;              /* Depends on `cltprm->snind' */
@@ -594,7 +571,7 @@ clumps_make_sn_table(struct clumps_thread_params *cltprm)
          noise cases) or the area is smaller than the minimum area to
          calculate signal-to-noise, then set the S/N of this segment to
          zero. */
-      if( Ni>p->segsnminarea && I>O )   /* This is O, not 0 (zero). */
+      if( I>O && Ni>p->segsnminarea )
         {
           /* Here we have done sky subtraction once. However, if the sky
              was already subtracted (informed by the user), then the
@@ -614,6 +591,7 @@ clumps_make_sn_table(struct clumps_thread_params *cltprm)
         }
       else
         {
+          /* Only over detections, we should put a NaN when the S/N  */
           if(sky0_det1)
             {
               snarr[i]=NAN;
@@ -714,7 +692,7 @@ clumps_find_make_sn_table(void *in_prm)
 
   /* Initialize the parameters for this thread. */
   cltprm.clprm   = clprm;
-  cltprm.topinds = clprm->sky0_det1 ? NULL : NULL;
+  cltprm.topinds = NULL;
 
 
   /* Go over all the tiles given to this thread. */
@@ -800,13 +778,15 @@ clumps_find_make_sn_table(void *in_prm)
                                  NULL, 0, 1);
 
           /* For a check, the step variable will be set. */
-          if(clprm->step==1) continue;
+          if(clprm->step==1)
+            { gal_data_free(cltprm.indexs); continue; }
 
           /* Make the clump S/N table. */
           clumps_make_sn_table(&cltprm);
 
           /* If the user wanted to check the steps, remove the clumps that
-             weren't used. */
+             weren't used from the `clabel' image (they have been already
+             excluded from the table). */
           if(cltprm.snind)
             clumps_correct_sky_labels_for_check(&cltprm, tile);
 
@@ -845,10 +825,69 @@ clumps_find_make_sn_table(void *in_prm)
 
 
 
+/***********************************************************************/
+/*****************           Over detections           *****************/
+/***********************************************************************/
+/* Put the indexs of each labeled region into an array of `gal_data_t's
+   (where each element is a dataset containing the respective label's
+   indexs). */
+gal_data_t *
+clumps_label_indexs(struct noisechiselparams *p)
+{
+  size_t i, *areas;
+  uint32_t *a, *l, *lf;
+  gal_data_t *labindexs=gal_data_array_calloc(p->numinitdets+1);
+
+  /* Find the area in each detected objects (to see how much space we need
+     to allocate). */
+  areas=gal_data_calloc_array(GAL_TYPE_SIZE_T, p->numinitdets+1);
+  lf=(l=p->olabel->array)+p->olabel->size; do ++areas[*l++]; while(l<lf);
+
+  /* Allocate/Initialize the dataset containing the indexs of each
+     object. We don't want the labels of the non-detected regions
+     (areas[0]). So we'll set that to zero.*/
+  areas[0]=0;
+  for(i=1;i<p->numinitdets+1;++i)
+    gal_data_initialize(&labindexs[i], NULL, GAL_TYPE_SIZE_T, 1,
+                        &areas[i], NULL, 0, p->cp.minmapsize, NULL, NULL,
+                        NULL);
+
+  /* Put the indexs into each dataset. We will use the areas array again,
+     but this time, use it as a counter. */
+  memset(areas, 0, (p->numinitdets+1)*sizeof *areas);
+  lf=(a=l=p->olabel->array)+p->olabel->size;
+  do
+    if(*l)  /* We don't want the undetected regions (*l==0) */
+      ((size_t *)(labindexs[*l].array))[ areas[*l]++ ] = l-a;
+  while(++l<lf);
 
-/**********************************************************************/
-/*****************         High level functions        *****************/
-/**********************************************************************/
+  /* Clean up and return. */
+  free(areas);
+  return labindexs;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/***********************************************************************/
+/*****************         High level function         *****************/
+/***********************************************************************/
 /* The job of this function is to find the best signal to noise value to
    use as a threshold to detect real clumps.
 
@@ -864,17 +903,15 @@ clumps_find_make_sn_table(void *in_prm)
    Using these two arrays, after all the threads are finished, we can
    concatenate all the S/N values into one array and send it to the main
    findsnthresh function in thresh.c. */
-#define CLUMPS_SNEXTNAME "CLUMPS_FOR_SN"
-float
-clumps_on_undetected_sn(struct noisechiselparams *p)
+void
+clumps_true_find_sn_thresh(struct noisechiselparams *p)
 {
   char *msg;
-  float snthresh;
   struct timeval t1;
   size_t i, j, c, numsn=0;
   struct clumps_params clprm;
   struct gal_linkedlist_stll *comments=NULL;
-  gal_data_t *clab, *claborig, *sn, *snind, *quant;
+  gal_data_t *claborig, *sn, *snind, *quant;
 
   /* Get starting time for later reporting if necessary. */
   if(!p->cp.quiet) gettimeofday(&t1, NULL);
@@ -903,14 +940,18 @@ clumps_on_undetected_sn(struct noisechiselparams *p)
     }
 
 
-  /* Spin off the threads to start the work. */
+  /* Spin off the threads to start the work. Note that several steps are
+     done on each tile within a thread. So if the user wants to check
+     steps, we need to break out of the processing get an over-all output,
+     then reset the input and call it again. So it will be slower, but its
+     is natural, since the user is testing to find the correct combination
+     of parameters for later use. */
   if(p->segmentationname)
     {
       /* Necessary initializations. */
-      clab=gal_data_copy(p->clabel);
-      claborig=p->clabel;
-      p->clabel=clab;
       clprm.step=1;
+      claborig=p->clabel;
+      p->clabel=gal_data_copy(claborig);
 
       /* Do each step. */
       while(clprm.step<3)
@@ -918,7 +959,7 @@ clumps_on_undetected_sn(struct noisechiselparams *p)
           /* Reset the temporary copy of clabel back to its original. */
           if(clprm.step>1)
             memcpy(p->clabel->array, claborig->array,
-                   clab->size*gal_type_sizeof(clab->type));
+                   claborig->size*gal_type_sizeof(claborig->type));
 
           /* Do this step. */
           gal_threads_spin_off(clumps_find_make_sn_table, &clprm,
@@ -927,12 +968,12 @@ clumps_on_undetected_sn(struct noisechiselparams *p)
           /* Set the extension name. */
           switch(clprm.step)
             {
-            case 1: p->clabel->name = "CLUMPS_ALL";      break;
-            case 2: p->clabel->name = CLUMPS_SNEXTNAME;  break;
+            case 1: p->clabel->name = "CLUMPS_ALL_SKY";  break;
+            case 2: p->clabel->name = "CLUMPS_FOR_SN";   break;
             default:
               error(EXIT_FAILURE, 0, "a bug! the value %d is not recognized "
-                    "in `clumps_on_undetected_sn'. Please contact us at %s "
-                    "so we can address the issue", clprm.step,
+                    "in `clumps_true_find_sn_thresh'. Please contact us at "
+                    "%s so we can address the issue", clprm.step,
                     PACKAGE_BUGREPORT);
             }
 
@@ -969,10 +1010,12 @@ clumps_on_undetected_sn(struct noisechiselparams *p)
 
   /* Allocate the space to keep all the S/N values. */
   sn=gal_data_alloc(NULL, clprm.sn->type, 1, &numsn, NULL, 0,
-                    p->cp.minmapsize, "CLUMP_S/N", "ratio", NULL);
+                    p->cp.minmapsize, "CLUMP_S/N", "ratio",
+                    "Signal-to-noise ratio");
   snind = ( p->checkclumpsn
             ? gal_data_alloc(NULL, clprm.snind->type, 1, &numsn, NULL, 0,
-                             p->cp.minmapsize, "CLUMP_ID", "counter", NULL)
+                             p->cp.minmapsize, "CLUMP_ID", "counter",
+                             "Unique ID for this clump.")
             : NULL );
 
 
@@ -998,9 +1041,8 @@ clumps_on_undetected_sn(struct noisechiselparams *p)
         gal_linkedlist_add_to_stll(&comments, "NOTE: In multi-threaded mode, "
                                    "clump IDs differ in each run and are not "
                                    "sorted.", 1);
-      gal_linkedlist_add_to_stll(&comments, "See also: `"CLUMPS_SNEXTNAME
-                                 "' HDU of output with `--checksegmentation'",
-                                 1);
+      gal_linkedlist_add_to_stll(&comments, "See also: `CLUMPS_FOR_SN' HDU "
+                                 "of output with `--checksegmentation'.", 1);
       gal_linkedlist_add_to_stll(&comments, "S/N of clumps over undetected "
                                  "regions.", 1);
       threshold_write_sn_table(p, sn, snind, p->clumpsn_s_name, comments);
@@ -1010,11 +1052,11 @@ clumps_on_undetected_sn(struct noisechiselparams *p)
 
   /* Find the desired quantile from the full S/N distribution. */
   quant=gal_statistics_quantile(sn, p->segquant, 1);
-  snthresh = *((float *)(quant->array));
+  p->clumpsnthresh = *((float *)(quant->array));
   if(!p->cp.quiet)
     {
       asprintf(&msg, "Clump S/N: %.2f (%.3f quant of %zu).",
-               snthresh, p->segquant, sn->size);
+               p->clumpsnthresh, p->segquant, sn->size);
       gal_timing_report(&t1, msg, 2);
       free(msg);
     }
@@ -1025,7 +1067,4 @@ clumps_on_undetected_sn(struct noisechiselparams *p)
   gal_data_free(quant);
   gal_data_array_free(clprm.sn, p->ltl.tottiles, 1);
   gal_data_array_free(clprm.snind, p->ltl.tottiles, 1);
-
-  /* Return the S/N threshold. */
-  return snthresh;
 }
diff --git a/bin/noisechisel/clumps.h b/bin/noisechisel/clumps.h
index f0b3fb8..466eb62 100644
--- a/bin/noisechisel/clumps.h
+++ b/bin/noisechisel/clumps.h
@@ -23,7 +23,48 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #ifndef CLUMPS_H
 #define CLUMPS_H
 
-float
-clumps_on_undetected_sn(struct noisechiselparams *p);
+
+/* Parameters for all threads. */
+struct clumps_params
+{
+  /* General */
+  int                     step; /* Counter if we want to check steps.      */
+  int                sky0_det1; /* If working on the Sky or Detections.    */
+  struct noisechiselparams  *p; /* Pointer to main NoiseChisel parameters. */
+  pthread_mutex_t     labmutex; /* Mutex to change the total numbers.      */
+
+  /* For Sky region. */
+  gal_data_t               *sn; /* Array of clump S/N tables.              */
+  gal_data_t            *snind; /* Array of clump S/N index (for check).   */
+
+  /* For detections. */
+  gal_data_t        *labindexs; /* Array of `gal_data_t' with obj indexs.  */
+};
+
+
+/* Parameters for one thread. */
+struct clumps_thread_params
+{
+  size_t                    id; /* ID of this detection/tile over tile.    */
+  size_t              *topinds; /* Indexs of all local maxima.             */
+  size_t            numinitial; /* The number of clumps found in this run. */
+  gal_data_t           *indexs; /* Array containing indexs of this object. */
+  gal_data_t             *info; /* Information for all clumps.             */
+  gal_data_t               *sn; /* Signal-to-noise ratio for these clumps. */
+  gal_data_t            *snind; /* Index of S/N for these clumps.          */
+  struct clumps_params  *clprm; /* Pointer to main structure.              */
+};
+
+void
+clumps_oversegment(struct clumps_thread_params *cltprm);
+
+void
+clumps_true_find_sn_thresh(struct noisechiselparams *p);
+
+void
+clumps_make_sn_table(struct clumps_thread_params *cltprm);
+
+gal_data_t *
+clumps_label_indexs(struct noisechiselparams *p);
 
 #endif
diff --git a/bin/noisechisel/detection.c b/bin/noisechisel/detection.c
index a2f8678..fef2c5b 100644
--- a/bin/noisechisel/detection.c
+++ b/bin/noisechisel/detection.c
@@ -110,7 +110,7 @@ detection_initial(struct noisechiselparams *p)
 
 
   /* Label the connected components. */
-  p->numobjects=gal_binary_connected_components(p->binary, &p->olabel, 1);
+  p->numinitdets=gal_binary_connected_components(p->binary, &p->olabel, 1);
   if(p->detectionname)
     {
       p->olabel->name="OPENED-LABELED";
@@ -122,7 +122,7 @@ detection_initial(struct noisechiselparams *p)
   /* Report the ending of initial detection. */
   if(!p->cp.quiet)
     {
-      asprintf(&msg, "%zu initial detections found.", p->numobjects);
+      asprintf(&msg, "%zu initial detections found.", p->numinitdets);
       gal_timing_report(&t0, msg, 1);
       free(msg);
     }
@@ -585,7 +585,7 @@ detection_pseudo_sn(struct noisechiselparams *p, gal_data_t 
*worklab,
 static void
 detection_pseudo_remove_low_sn(struct noisechiselparams *p,
                                gal_data_t *workbin, gal_data_t *worklab,
-                               gal_data_t *sn, float snthresh)
+                               gal_data_t *sn)
 {
   size_t i;
   float *snarr=sn->array;
@@ -598,7 +598,7 @@ detection_pseudo_remove_low_sn(struct noisechiselparams *p,
      (which will fail on any condition), so it acts as if it had an S/N
      lower than the required value. */
   for(i=0;i<sn->size;++i)
-    if( snarr[i] > snthresh ) keep[i]=1;
+    if( snarr[i] > p->detsnthresh ) keep[i]=1;
 
 
   /* Go over the pseudo-detection labels and only keep those that must be
@@ -629,7 +629,6 @@ static gal_data_t *
 detection_pseudo_real(struct noisechiselparams *p)
 {
   char *msg;
-  float snthresh;
   size_t numpseudo;
   struct timeval t1;
   gal_data_t *sn, *quant, *workbin, *worklab;
@@ -649,11 +648,11 @@ detection_pseudo_real(struct noisechiselparams *p)
 
   /* Get the S/N quantile and report it if we are in non-quiet mode. */
   quant=gal_statistics_quantile(sn, p->detquant, 1);
-  snthresh = *((float *)(quant->array));
+  p->detsnthresh = *((float *)(quant->array));
   if(!p->cp.quiet)
     {
       asprintf(&msg, "Pseudo-det S/N: %.2f (%.3f quant of %zu).",
-               snthresh, p->detquant, sn->size);
+               p->detsnthresh, p->detquant, sn->size);
       gal_timing_report(&t1, msg, 2);
       free(msg);
     }
@@ -668,7 +667,7 @@ detection_pseudo_real(struct noisechiselparams *p)
 
 
   /* Remove the pseudo detections with a low S/N. */
-  detection_pseudo_remove_low_sn(p, workbin, worklab, sn, snthresh);
+  detection_pseudo_remove_low_sn(p, workbin, worklab, sn);
 
 
   /* Clean up and return. */
@@ -705,7 +704,8 @@ detection_remove_false_initial(struct noisechiselparams *p,
   size_t i;
   uint8_t *b=workbin->array;
   uint32_t *l=p->olabel->array, *lf=l+p->olabel->size, curlab=1;
-  uint32_t *newlabels=gal_data_calloc_array(GAL_TYPE_UINT32, p->numobjects+1);
+  uint32_t *newlabels=gal_data_calloc_array(GAL_TYPE_UINT32,
+                                            p->numinitdets+1);
 
   /* Find the new labels for all the existing labels. Recall that
      `newlabels' was initialized to zero, so any label that is not given a
@@ -734,7 +734,7 @@ detection_remove_false_initial(struct noisechiselparams *p,
 
   /* Now that we know which labels to keep, set the new labels for the
      detections that must be kept. */
-  for(i=0;i<p->numobjects;++i) if(newlabels[i]) newlabels[i] = curlab++;
+  for(i=0;i<p->numinitdets;++i) if(newlabels[i]) newlabels[i] = curlab++;
 
 
   /* Replace the byt and olab values with their proper values. If the
@@ -817,7 +817,7 @@ detection(struct noisechiselparams *p)
   if(!p->cp.quiet)
     {
       asprintf(&msg, "%zu false initial detections removed.",
-               p->numobjects - num_true_initial);
+               p->numinitdets - num_true_initial);
       gal_timing_report(&t1, msg, 2);
       free(msg);
     }
@@ -827,13 +827,14 @@ detection(struct noisechiselparams *p)
   if(p->dilate)
     {
       gal_binary_dilate(workbin, p->dilate, 8, 1);
-      p->numobjects = gal_binary_connected_components(workbin, &p->olabel, 8);
+      p->numinitdets = gal_binary_connected_components(workbin, &p->olabel,
+                                                       8);
     }
-  else p->numobjects=num_true_initial;
+  else p->numinitdets=num_true_initial;
   if(!p->cp.quiet)
     {
       asprintf(&msg, "%zu detections after %zu dilation%s",
-              p->numobjects, p->dilate, p->dilate>1 ? "s." : ".");
+              p->numinitdets, p->dilate, p->dilate>1 ? "s." : ".");
       gal_timing_report(&t0, msg, 1);
       free(msg);
     }
diff --git a/bin/noisechisel/main.h b/bin/noisechisel/main.h
index 30d0cd3..d1f67c2 100644
--- a/bin/noisechisel/main.h
+++ b/bin/noisechisel/main.h
@@ -104,6 +104,10 @@ struct noisechiselparams
   gal_data_t          *clabel;  /* Labels of clumps in the detection.     */
   gal_data_t             *sky;  /* Mean of undetected pixels, per tile.   */
   gal_data_t             *std;  /* STD of undetected pixels, per tile.    */
+  size_t           maxtcontig;  /* Maximum contiguous space for a tile.   */
+  size_t          maxltcontig;  /* Maximum contiguous space for a tile.   */
+  size_t            *maxtsize;  /* Maximum size of a single small tile.   */
+  size_t           *maxltsize;  /* Maximum size of a single large tile.   */
   time_t              rawtime;  /* Starting time of the program.          */
 
   float                medstd;  /* Median STD before interpolation.       */
@@ -111,12 +115,10 @@ struct noisechiselparams
   float                maxstd;  /* Maximum STD before interpolation.      */
   float               cpscorr;  /* Counts/second correction.              */
 
-  size_t           numobjects;  /* Number of objects detected.            */
+  size_t          numinitdets;  /* Number of objects detected.            */
   size_t            numclumps;  /* Number of true clumps.                 */
-  size_t           maxtcontig;  /* Maximum contiguous space for a tile.   */
-  size_t          maxltcontig;  /* Maximum contiguous space for a tile.   */
-  size_t            *maxtsize;  /* Maximum size of a single small tile.   */
-  size_t           *maxltsize;  /* Maximum size of a single large tile.   */
+  float           detsnthresh;  /* Pseudo-detection S/N threshold.        */
+  float         clumpsnthresh;  /* Clump S/N threshold.                   */
 };
 
 #endif
diff --git a/bin/noisechisel/segmentation.c b/bin/noisechisel/segmentation.c
index c0be1c1..a800ae8 100644
--- a/bin/noisechisel/segmentation.c
+++ b/bin/noisechisel/segmentation.c
@@ -26,9 +26,11 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <errno.h>
 #include <error.h>
 #include <stdlib.h>
+#include <string.h>
 
 #include <gnuastro/fits.h>
 #include <gnuastro/blank.h>
+#include <gnuastro/threads.h>
 
 #include <gnuastro-internal/timing.h>
 
@@ -39,11 +41,284 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include "segmentation.h"
 
 
+
+
+
+/***********************************************************************/
+/*****************            Over detections          *****************/
+/***********************************************************************/
+/* Find the true clumps over each detection. */
+static void *
+segmentation_on_threads(void *in_prm)
+{
+  struct gal_threads_params *tprm=(struct gal_threads_params *)in_prm;
+  struct clumps_params *clprm=(struct clumps_params *)(tprm->params);
+  struct noisechiselparams *p=clprm->p;
+
+  size_t i;
+  gal_data_t *topinds;
+  struct clumps_thread_params cltprm;
+
+  /* Initialize the general parameters for this thread. */
+  cltprm.clprm   = clprm;
+
+  /* Go over all the detections given to this thread (counting from zero.) */
+  for(i=0; tprm->indexs[i] != GAL_THREADS_NON_THRD_INDEX; ++i)
+    {
+
+      /* Set the ID of this detection, note that for the threads, we
+         counted from zero, but the IDs start from 1, so we'll add a 1 to
+         the ID given to this thread. */
+      cltprm.id     = tprm->indexs[i]+1;
+      cltprm.indexs = &clprm->labindexs[ cltprm.id ];
+
+
+      /* The `topinds' array is only necessary when the user wants to
+         ignore true clumps with a peak touching a river. */
+      if(p->keepmaxnearriver==0)
+        {
+          /* Allocate the list of local maxima. For each clump there is
+             going to be one local maxima. But we don't know the number of
+             clumps a-priori, so we'll just allocate the number of pixels
+             given to this detected region. */
+          topinds=gal_data_alloc(NULL, GAL_TYPE_SIZE_T, 1,
+                                 cltprm.indexs->dsize, NULL, 0,
+                                 p->cp.minmapsize, NULL, NULL, NULL);
+          cltprm.topinds=topinds->array;
+        }
+      else { cltprm.topinds=NULL; topinds=NULL; }
+
+
+      /* Find the clumps over this region. */
+      clumps_oversegment(&cltprm);
+
+      /* Make the clump S/N table. This table is made before (possibly)
+         stopping the process (if a check is requested). This is because if
+         the user has also asked for a check image, we can break out of the
+         loop at that point.*/
+      clumps_make_sn_table(&cltprm);
+
+      /* If the user wanted to check the segmentation steps or the clump
+         S/N values in a table, then we have to stop the process at this
+         point. */
+      if(clprm->step==1 || p->checkclumpsn)
+        { gal_data_free(topinds); continue; }
+
+
+
+      /* Clean up. */
+      gal_data_free(topinds);
+    }
+
+  if(tprm->b) pthread_barrier_wait(tprm->b);
+  return NULL;
+}
+
+
+
+
+
+/* If the user wanted to see the S/N table in a file, this function will be
+   called and will do the job. */
+static void
+segmentation_save_sn_table(struct clumps_params *clprm)
+{
+  char *msg;
+  float *sarr;
+  uint32_t *oiarr, *cioarr;
+  size_t i, j, c=0, totclumps=0;
+  gal_data_t *sn, *objind, *clumpinobj;
+  struct noisechiselparams *p=clprm->p;
+  struct gal_linkedlist_stll *comments=NULL;
+
+
+  /* Find the total number of clumps in all the initial detections. Recall
+     that the `size' values were one more than the actual number because
+     the labelings start from 1. */
+  for(i=1;i<p->numinitdets+1;++i) totclumps += clprm->sn[i].size-1;
+
+
+  /* Allocate the columns for the table. */
+  sn=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, 1, &totclumps, NULL, 0,
+                    p->cp.minmapsize, "CLUMP_S/N", "ratio",
+                    "Signal-to-noise ratio.");
+  objind=gal_data_alloc(NULL, GAL_TYPE_UINT32, 1, &totclumps, NULL, 0,
+                        p->cp.minmapsize, "HOST_DET_ID", "counter",
+                        "ID of detection hosting this clump.");
+  clumpinobj=gal_data_alloc(NULL, GAL_TYPE_UINT32, 1, &totclumps, NULL, 0,
+                            p->cp.minmapsize, "CLUMP_ID_IN_OBJ", "counter",
+                            "ID of clump in host detection.");
+
+
+  /* Fill in the columns. */
+  sarr=sn->array;
+  oiarr=objind->array;
+  cioarr=clumpinobj->array;
+  for(i=1;i<p->numinitdets+1;++i)
+    for(j=1;j<clprm->sn[i].size;++j)
+      {
+        oiarr[c]  = i;
+        cioarr[c] = j;
+        sarr[c]   = ((float *)(clprm->sn[i].array))[j];
+        ++c;
+      }
+
+
+  /* Write the comments. */
+  gal_linkedlist_add_to_stll(&comments, "See also: `CLUMPS_ALL_DET' HDU of "
+                             "output with `--checksegmentation'.", 1);
+  asprintf(&msg, "S/N values of `nan': clumps smaller than `--segsnminarea' "
+           "of %zu.", p->segsnminarea);
+  gal_linkedlist_add_to_stll(&comments, msg, 1);
+  free(msg);
+  gal_linkedlist_add_to_stll(&comments, "S/N of clumps over detected "
+                             "regions.", 1);
+  gal_table_comments_add_intro(&comments, PROGRAM_STRING, &p->rawtime);
+
+
+  /* Set the column pointers and write them into a table.. */
+  clumpinobj->next=sn;
+  objind->next=clumpinobj;
+  gal_table_write(objind, comments, p->cp.tableformat, p->clumpsn_d_name, 1);
+
+
+  /* Clean up. */
+  gal_data_free(sn);
+  gal_data_free(objind);
+  gal_data_free(clumpinobj);
+  gal_linkedlist_free_stll(comments, 1);
+
+
+  /* Abort NoiseChisel if necessary. */
+  if(!p->continueaftercheck)
+    ui_abort_after_check(p, p->clumpsn_s_name, p->clumpsn_d_name,
+                         "showing all clump S/N values");
+}
+
+
+
+
+
+/* Find true clumps over the detected regions. */
+void
+segmentation_detections(struct noisechiselparams *p)
+{
+  struct clumps_params clprm;
+  gal_data_t *labindexs, *claborig;
+
+
+  /* Get the indexs of all the pixels in each label. */
+  labindexs=clumps_label_indexs(p);
+
+
+  /* Initialize the necessary thread parameters. Note that since the object
+     labels begin from one, the `sn' array will have one extra element.*/
+  clprm.p=p;
+  clprm.sky0_det1=1;
+  clprm.snind = NULL;
+  clprm.labindexs=labindexs;
+  clprm.sn=gal_data_array_calloc(p->numinitdets+1);
+
+
+  /* Spin off the threads to start the work. Note that several steps are
+     done on each tile within a thread. So if the user wants to check
+     steps, we need to break out of the processing get an over-all output,
+     then reset the input and call it again. So it will be slower, but its
+     is natural, since the user is testing to find the correct combination
+     of parameters for later use. */
+  if(p->segmentationname)
+    {
+      /* Necessary initializations. */
+      clprm.step=1;
+      claborig=p->clabel;
+      p->clabel=gal_data_copy(claborig);
+
+      /* Do each step. */
+      while(clprm.step<3)
+        {
+          /* Reset the temporary copy of clabel back to its original. */
+          if(clprm.step>1)
+            memcpy(p->clabel->array, claborig->array,
+                   claborig->size*gal_type_sizeof(claborig->type));
+
+          /* (Re-)do everything until this step. */
+          gal_threads_spin_off(segmentation_on_threads, &clprm,
+                               p->numinitdets, p->cp.numthreads);
+
+          /* Set the extension name. */
+          switch(clprm.step)
+            {
+            case 1: p->clabel->name = "CLUMPS_ALL_DET";      break;
+            case 2: p->clabel->name = "TRUE_CLUMPS";         break;
+            default:
+              error(EXIT_FAILURE, 0, "a bug! the value %d is not recognized "
+                    "in `segmentation_detections'. Please contact us at "
+                    "%s so we can address the issue", clprm.step,
+                    PACKAGE_BUGREPORT);
+            }
+
+          /* Write the temporary array into the check image. */
+          gal_fits_img_write(p->clabel, p->segmentationname, NULL,
+                             PROGRAM_STRING);
+
+          /* If the user wanted to check the clump S/N values, then break
+             out of the loop, we don't need the rest of the process any
+             more. */
+          if( clprm.step==1
+              && ( p->checkclumpsn && !p->continueaftercheck ) ) break;
+
+          /* Increment the step counter. */
+          ++clprm.step;
+        }
+
+      /* Clean up (we don't need the original any more). */
+      gal_data_free(claborig);
+      p->clabel->name=NULL;
+    }
+  else
+    {
+      clprm.step=0;
+      gal_threads_spin_off(segmentation_on_threads, &clprm, p->numinitdets,
+                           p->cp.numthreads);
+    }
+
+
+  /* If the user wanted to see the S/N table, then make the S/N table. */
+  if(p->checkclumpsn) segmentation_save_sn_table(&clprm);
+
+
+  /* Free the array of indexs. */
+  gal_data_array_free(clprm.sn, p->numinitdets, 1);
+  gal_data_array_free(labindexs, p->numinitdets, 1);
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/***********************************************************************/
+/*****************         High level function         *****************/
+/***********************************************************************/
 void
 segmentation(struct noisechiselparams *p)
 {
+  float *f;
   uint32_t *l, *lf;
-  float *f, snthresh;
 
   /* To keep the user up to date. */
   if(!p->cp.quiet)
@@ -73,8 +348,17 @@ segmentation(struct noisechiselparams *p)
   do if(isnan(*f++)) *l=GAL_BLANK_UINT32; while(++l<lf);
 
 
-  /* Find the clumps over the un-detected regions of the input. */
-  snthresh=clumps_on_undetected_sn(p);
+  /* Find the clump S/N threshold. */
+  clumps_true_find_sn_thresh(p);
+
+
+  /* Reset the clabel array to find true clumps in objects. */
+  f=p->input->array; lf=(l=p->clabel->array)+p->clabel->size;
+  do *l = isnan(*f++) ? GAL_BLANK_UINT32 : 0; while(++l<lf);
+
+
+  /* Find true clumps over the detected regions. */
+  segmentation_detections(p);
 
 
   /* If the user wanted to check the segmentation and hasn't called
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 3851962..a9fa45d 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -2401,15 +2401,15 @@ explains how to check if it has been installed 
correctly.
 CFITSIO comes with two executables called fpack and funpack. From their
 manual: they ``are standalone programs for compressing and uncompressing
 images and tables that are stored in the FITS (Flexible Image Transport
-System) data format.  They are analogous to the gzip and gunzip compression
+System) data format. They are analogous to the gzip and gunzip compression
 programs except that they are optimized for the types of astronomical
-images that are often stored in FITS format''. The instruction below
-explain how to compile and install them on your system along with
-CFITSIO. They are not essential for Gnuastro, since they are just wrappers
-for functions within CFITSIO, but they can come in handy. The @command{make
-utils} command is only available for versions above 3.39, it will build
-these executables along with several other test executables which are
-deleted before the installation (otherwise they will also be installed).
+images that are often stored in FITS format''. The commands below will
+compile and install them on your system along with CFITSIO. They are not
+essential for Gnuastro, since they are just wrappers for functions within
+CFITSIO, but they can come in handy. The @command{make utils} command is
+only available for versions above 3.39, it will build these executables
+along with several other test executables which are deleted before the
+installation (otherwise they will also be installed).
 
 The CFITSIO installation from source process is given below. Let's assume
 you have downloaded
@@ -2446,9 +2446,9 @@ address@hidden E.W., Calabretta M.R. (2002) 
Representation of
 world coordinates in FITS. Astronomy and Astrophysics, 395, 1061-1075.},
 Mark Calabretta. It might be already built and ready in your distribution's
 package management system. However, here the installation from source is
-explained, for the advantages please see @ref{Mandatory dependencies}. To
-install WCSLIB you will need to have CFITSIO already installed, see
address@hidden
+explained, for the advantages of installation from source please see
address@hidden dependencies}. To install WCSLIB you will need to have
+CFITSIO already installed, see @ref{CFITSIO}.
 
 @vindex --without-pgplot
 WCSLIB also has plotting capabilities which use PGPLOT (a plotting library
@@ -3388,15 +3388,16 @@ can use any other directory, but this is the most 
recognized.
 
 The top installation directory will be used to keep all the package's
 components: programs (executables), libraries, include (header) files,
-shared data (like manuals), or configuration files. So it commonly has some
-of the following sub-directories for each class of components respectively:
address@hidden/}, @file{lib/}, @file{include/} @file{man/}, @file{share/},
address@hidden/}. Since the @file{PATH} variable is only used for executables,
-you can add the @file{~/.local/bin} directory to @file{PATH} with the
-following command. As defined below, first the existing value of
address@hidden is used, then your given directory is added to its end and the
-combined value is put back in @file{PATH} (run address@hidden echo $PATH}'
-afterwards to check if it was added).
+shared data (like manuals), or configuration files (see @ref{Review of
+library fundamentals} for a thorough introduction to headers and
+linking). So it commonly has some of the following sub-directories for each
+class of components respectively: @file{bin/}, @file{lib/}, @file{include/}
address@hidden/}, @file{share/}, @file{etc/}. Since the @file{PATH} variable is
+only used for executables, you can add the @file{~/.local/bin} directory to
address@hidden with the following command. As defined below, first the
+existing value of @file{PATH} is used, then your given directory is added
+to its end and the combined value is put back in @file{PATH} (run
address@hidden echo $PATH}' afterwards to check if it was added).
 
 @example
 $ PATH=$PATH:~/.local/bin
@@ -3407,11 +3408,12 @@ $ PATH=$PATH:~/.local/bin
 @cindex Scripts, startup
 Any executable that you installed in this directory will now be usable
 without having to remember and type its full address. However, as soon as
-you leave your current terminal session, this modification will be
-lost. Adding your specified directory to the @file{PATH} environment
-variable each time you start a terminal is also very inconvenient and prone
-to errors. So there are standard `startup files' defined by your
-shell. There is a special startup file for every significant starting step:
+you leave your current terminal session, this modified @file{PATH} variable
+will be forgotten. Adding your specified directory to the @file{PATH}
+environment variable each time you start a terminal is also very
+inconvenient and prone to errors. So there are standard `startup files'
+defined by your shell. There is a special startup file for every
+significant starting step:
 
 @table @asis
 
@@ -3436,13 +3438,13 @@ For security reasons, it is highly recommended to 
directly type in your
 @file{HOME} directory value by hand in startup files instead of using
 variables. So in the following, let's assume your user name is
 address@hidden'. To add @file{~/.local/bin} to your @file{PATH} automatically
-on every startup you have to ``export'' the new value of @command{PATH} in
-the startup file that is most relevant to you with the line
+on any startup file, you have to ``export'' the new value of @command{PATH}
+in the startup file that is most relevant to you with the line
 address@hidden PATH=$PATH:/home/name/.local/bin}'. You can either do it
 manually using a text editor, or by running the following command which
 will add this line as the last line of the file. Let's assume you want to
 add it to @file{~/.bashrc} (afterwards, open your @file{~/.bashrc} with a
-text editor and check it out to see for your self):
+text editor and check it out, to see the result for your self):
 
 @example
 $ echo 'export PATH=$PATH:/home/name/.local/bin' >> ~/.bashrc
@@ -3463,8 +3465,8 @@ install-able files will be put in their respective 
directory under
 @file{~/.local/} (as discussed above). Note that tilde (address@hidden')
 expansion will not happen if you put a address@hidden' between 
@option{--prefix}
 and @file{~/address@hidden you insist on using address@hidden', you can use
address@hidden/.local}.}, so we have avoided it here, see
address@hidden
address@hidden/.local}.}, so we have avoided the @key{=} character
+here which is optional in GNU-style options, see @ref{Options}.
 
 @example
 $ ./configure --prefix ~/.local
@@ -3483,7 +3485,7 @@ CFITSIO, or Gnuastro which needs both.} that depend on 
libraries you
 installed like this, you have to guide your compiler to the necessary
 directories. To do that, you have to define the @command{LDFLAGS} and
 @command{CPPFLAGS} environment variables respectively. This can be done
-while calling @file{./configure}:
+while calling @file{./configure} as shown below:
 
 @example
 $ ./configure LDFLAGS=-L/home/name/.local/lib            \
@@ -3491,17 +3493,19 @@ $ ./configure LDFLAGS=-L/home/name/.local/lib           
 \
               --prefix ~/.local
 @end example
 
-It can be annoying to do this when configuring every package, so you can
-define these two variables in a startup file (discussed above). The
-convention on using these variables doesn't include a colon to separate
-values (as @command{PATH}-like variables do), they use white space
+If you do this for all your libraries, it can be annoying to do this when
+configuring every software that depends on them. Hence, you can define
+these two variables in the most relevant startup file (discussed
+above). The convention on using these variables doesn't include a colon to
+separate values (as @command{PATH}-like variables do), they use white space
 characters and each value is prefixed with a compiler address@hidden
 variables are ultimately used as options while building the programs, so
 every value has be an option name followed be a value as discussed in
 @ref{Options}.} (note the @option{-L} and @option{-I} above, see
address@hidden). Therefore we have to keep the value in double quotation
-signs to keep the white space characters. Run the following two commands to
-do that if you want them in your @file{~/.bashrc}.
address@hidden, @ref{Headers} for @option{-I}, and @ref{Linking} for
address@hidden). Therefore we have to keep the value in double quotation signs
+to keep the white space characters. Run the following two commands to do
+that if you want them in your @file{~/.bashrc}.
 
 @example
 echo 'export LDFLAGS="$LDFLAGS -L/home/name/.local/lib"' >> ~/.bashrc
@@ -3509,23 +3513,35 @@ echo 'export CPPFLAGS="$CPPFLAGS 
-I/home/name/.local/include"' \
      >> ~/.bashrc
 @end example
 
-To be able to link to the locally installed dynamic libraries (which are
-linked at run time) after installation, add @file{~/.local/lib} to your
address@hidden environment variable similar to the example above
-for @command{PATH}. If you also want to access the Info (see @ref{Info})
-and man pages (see @ref{Man pages}) documentations add
address@hidden/.local/share/info} and @file{~/.local/share/man} to your
address@hidden@footnote{Info has the following convention: ``If the
-value of @command{INFOPATH} ends with a colon [or it isn't defined] ...,
-the initial list of directories is constructed by appending the build-time
-default to the value of @command{INFOPATH}.'' So when installing in a
-non-standard directory and if @command{INFOPATH} was not initially defined,
-add a colon to the end of @command{INFOPATH} as shown below, otherwise Info
-will not be able to find system-wide installed
-documentation:@address@hidden 'export
+Dynamic libraries are linked to the executable every time you run a program
+that depends on them (see @ref{Linking} to fully understand this important
+concept). Hence dynamic libraries also require a special path variable
+called @command{LD_LIBRARY_PATH}. To use programs that depend on these
+libraries, you will need to have @file{~/.local/lib} added to your
address@hidden environment variable in the relevant start-up
+file (similar to @command{PATH}) as shown below (the @key{\} is only to fit
+this command within the page width, you can ignore it when typing on the
+terminal):
+
address@hidden
+echo 'export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/home/name/.local/lib' \
+     >> ~/.bashrc
address@hidden example
+
+If you also want to access the Info (see @ref{Info}) and man pages (see
address@hidden pages}) documentations add @file{~/.local/share/info} and
address@hidden/.local/share/man} to your @address@hidden has the
+following convention: ``If the value of @command{INFOPATH} ends with a
+colon [or it isn't defined] ..., the initial list of directories is
+constructed by appending the build-time default to the value of
address@hidden'' So when installing in a non-standard directory and if
address@hidden was not initially defined, add a colon to the end of
address@hidden as shown below, otherwise Info will not be able to find
+system-wide installed documentation:@address@hidden 'export
 INFOPATH=$INFOPATH:/home/name/.local/share/info:' >> ~/address@hidden Note that
 this is only an internal convention of Info, do not use it for other
address@hidden variables.} and @command{MANPATH} environment variables.
address@hidden variables.} and @command{MANPATH} environment variables
+respectively.
 
 @cindex Search directory order
 @cindex Order in search directory
@@ -3546,15 +3562,15 @@ This is good when a library, for example CFITSIO, is 
already present on the
 system, but the system-wide install wasn't configured with the correct
 configuration flags (see @ref{CFITSIO}), or you want to use a newer version
 and you don't have administrator or root access to update it. If you update
address@hidden with the order shown above, the linker will first
-find the CFITSIO you installed for yourself and use it, it will never reach
-the system-wide installation. However there are important security problems
-with this order: all important system-wide programs and libraries can be
-replaced by non-secure versions with the same file names and put in the
-customized directory (@file{./.local/} in this example). So if you choose
-to search in your customized directory first, be sure to keep it clean from
-executables or libraries with the same names as important system programs
-or libraries.
address@hidden with the @file{~/.local/lib} installation first, the
+linker will first find the CFITSIO you installed for yourself and link with
+it. It will never reach the system-wide installation. However there are
+important security problems with this order: all important system-wide
+programs and libraries can be replaced by non-secure versions with the same
+file names and put in the customized directory (@file{./.local/lib} in this
+example). So if you choose to search in your customized directory first, be
+sure to keep it clean from executables or libraries with the same names as
+important system programs or libraries.
 
 
 
diff --git a/lib/data.c b/lib/data.c
index 1c2bd0d..f389042 100644
--- a/lib/data.c
+++ b/lib/data.c
@@ -376,16 +376,20 @@ gal_data_initialize(gal_data_t *data, void *array, 
uint8_t type,
         data->array=array;
       else
         {
-          if( gal_type_sizeof(type)*data->size  > minmapsize )
-            gal_data_mmap(data, clear);
-          else
+          if(data->size)
             {
-              /* Allocate the space for the array. */
-              if(clear)
-                data->array = gal_data_calloc_array(data->type, data->size);
+              if( gal_type_sizeof(type)*data->size  > minmapsize )
+                /* Allocate the space into disk (HDD/SSD). */
+                gal_data_mmap(data, clear);
               else
-                data->array = gal_data_malloc_array(data->type, data->size);
+                /* Allocate the space in RAM. */
+                data->array = ( clear
+                                ? gal_data_calloc_array(data->type,
+                                                        data->size)
+                                : gal_data_malloc_array(data->type,
+                                                        data->size) );
             }
+          else data->array=NULL; /* The given size was zero! */
         }
     }
   else
@@ -484,15 +488,20 @@ gal_data_string_fixed_alloc_size(gal_data_t *data)
 
 
 /* Free the allocated contents of a data structure, not the structure
-   itsself. The reason for this function begin separate from
-   `gal_data_free) is that the data structure might be allocated as an
-   array (statically like `gal_data_t da[20]', or dynamically like
-   `gal_data_t *da; da=malloc(20*sizeof *da);'). In both cases, a loop will
-   be necessary to delete the allocated contents of each element of the
-   data structure array, but not the structure its self. After that loop,
-   if the array of data structures was statically allocated, you don't have
-   to do anything. If it was dynamically allocated, we just have to run
-   `free(da)'.*/
+   itsself. The reason that this function is separate from `gal_data_free'
+   is that the data structure might be allocated as an array (statically
+   like `gal_data_t da[20]', or dynamically like `gal_data_t *da;
+   da=malloc(20*sizeof *da);'). In both cases, a loop will be necessary to
+   delete the allocated contents of each element of the data structure
+   array, but not the structure its self. After that loop, if the array of
+   data structures was statically allocated, you don't have to do
+   anything. If it was dynamically allocated, we just have to run
+   `free(da)'.
+
+   Since we aren't freeing the `gal_data_t' its-self, after the allocated
+   space for each pointer is freed, the pointer is set to NULL for safety
+   (to avoid possible re-calls).
+*/
 void
 gal_data_free_contents(gal_data_t *data)
 {
@@ -501,11 +510,11 @@ gal_data_free_contents(gal_data_t *data)
           "`gal_data_free_contents' was a NULL pointer");
 
   /* Free all the possible allocations. */
-  if(data->name)    free(data->name);
-  if(data->unit)    free(data->unit);
-  if(data->dsize)   free(data->dsize);
-  if(data->wcs)     wcsfree(data->wcs);
-  if(data->comment) free(data->comment);
+  if(data->name)    { free(data->name);    data->name=NULL;    }
+  if(data->unit)    { free(data->unit);    data->unit=NULL;    }
+  if(data->dsize)   { free(data->dsize);   data->dsize=NULL;   }
+  if(data->wcs)     { wcsfree(data->wcs);  data->wcs=NULL;     }
+  if(data->comment) { free(data->comment); data->comment=NULL; }
 
   /* If the data type is string, then each element in the array is actually
      a pointer to the array of characters, so free them before freeing the
@@ -528,6 +537,7 @@ gal_data_free_contents(gal_data_t *data)
     }
   else
     if(data->array) free(data->array);
+  data->array=NULL;
 }