[gnuastro-commits] master c250a6d 123/125: NoiseChisel's clumps grown and objects defined

[Mohammad Akhlaghi]

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

    NoiseChisel's clumps grown and objects defined
    
    The new implementation of NoiseChisel will now fully grow the clumps,
    determine the connectivity of the clumps and identify which clumps should
    be merged into one object. It will finally expand the grown regions to
    completely label the diffuse detected area.
    
    The old function to find connected objects through an adjacency matrix
    which was only usable within NoiseChisel has now been moved to a separate
    library function: `gal_binary_connected_adjacency_matrix' (in
    `gnuastro/binary.h'). All the old NoiseChisel functions for the steps above
    have also been re-written in a much more cleaner way using the new tools,
    but they are specific to NoiseChisel and weren't be taken to the libraries.
---
 bin/noisechisel/clumps.c       | 530 ++++++++++++++++++++++++++++++---------
 bin/noisechisel/clumps.h       |  35 ++-
 bin/noisechisel/detection.c    |  16 +-
 bin/noisechisel/main.h         |   4 +-
 bin/noisechisel/segmentation.c | 544 +++++++++++++++++++++++++++++++++++++++--
 bin/noisechisel/ui.c           |   2 +-
 lib/binary.c                   | 111 +++++++++
 lib/convolve.c                 |   2 +-
 lib/gnuastro/binary.h          |   5 +-
 9 files changed, 1091 insertions(+), 158 deletions(-)

diff --git a/bin/noisechisel/clumps.c b/bin/noisechisel/clumps.c
index 688196c..ad20722 100644
--- a/bin/noisechisel/clumps.c
+++ b/bin/noisechisel/clumps.c
@@ -54,31 +54,6 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 
 
-/**********************************************************************/
-/*****************         Basic macro values         *****************/
-/**********************************************************************/
-/* Constants for the clump over-segmentation. */
-#define CLUMPS_RIVER     UINT32_MAX-2
-#define CLUMPS_TMPCHECK  UINT32_MAX-3
-#define CLUMPS_INIT      UINT32_MAX-4
-#define CLUMPS_MAXLAB    UINT32_MAX-5   /* Largest clump label (unsigned). */
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
 /****************************************************************
  *****************   Over segmentation       ********************
  ****************************************************************/
@@ -205,8 +180,8 @@ clumps_oversegment(struct clumps_thread_params *cltprm)
                          if( nlab==CLUMPS_INIT && arr[nind]==arr[*a] )
                            {
                              clabel[nind]=CLUMPS_TMPCHECK;
-                             gal_linkedlist_add_to_sll(&cleanup, nind);
                              gal_linkedlist_add_to_sll(&Q, nind);
+                             gal_linkedlist_add_to_sll(&cleanup, nind);
                            }
 
                          /* If this neighbour has a positive nlab, it
@@ -242,7 +217,7 @@ clumps_oversegment(struct clumps_thread_params *cltprm)
             else
               {
                 rlab = curlab++;
-                if( cltprm->topinds)        /* This is a local maximum of  */
+                if( cltprm->topinds )       /* This is a local maximum of  */
                   cltprm->topinds[rlab]=*a; /* this region, save its index.*/
               }
 
@@ -332,7 +307,7 @@ clumps_oversegment(struct clumps_thread_params *cltprm)
   while(++a<af);
 
   /* Save the total number of clumps. */
-  cltprm->numinitial=curlab-1;
+  cltprm->numinitclumps=curlab-1;
 
   /* Set all the river pixels to zero. Note that this is only necessary for
      the detected clumps. When finding clumps over the Sky, we will be
@@ -373,6 +348,282 @@ clumps_oversegment(struct clumps_thread_params *cltprm)
 
 
 
+
+
+/**********************************************************************/
+/*****************              Grow clumps           *****************/
+/**********************************************************************/
+/* Make the preparations for the intiial growing the clumps to identify
+   objects: a single standard deviation for the whole object and preparing
+   the labels (because the growth is going to happen on the `olabel'
+   image. */
+void
+clumps_grow_prepare_initial(struct clumps_thread_params *cltprm)
+{
+  gal_data_t *indexs=cltprm->indexs;
+  gal_data_t *input=cltprm->clprm->p->input;
+  struct noisechiselparams *p=cltprm->clprm->p;
+
+  size_t *s, *sf, *dsize=input->dsize;
+  size_t ndiffuse=0, coord[2], *dindexs;
+  double wcoord[2]={0.0f,0.0f}, brightness=0.0f;
+  float glimit, *imgss=input->array, *std=p->std->array;
+  uint32_t *olabel=p->olabel->array, *clabel=p->clabel->array;
+
+
+  /* Find the flux weighted center (meaningful only for positive valued
+     pixels). */
+  sf=(s=indexs->array)+indexs->size;
+  do
+    if( imgss[ *s ] > 0.0f )
+      {
+        brightness += imgss[ *s ];
+        wcoord[0]  += imgss[ *s ] * (*s/dsize[1]);
+        wcoord[1]  += imgss[ *s ] * (*s%dsize[1]);
+      }
+  while(++s<sf);
+
+
+  /* Calculate the center, if no pixels were positive, use the
+     geometric center (irrespective of flux). */
+  if(brightness==0.0f)
+    {
+      sf=(s=indexs->array)+indexs->size;
+      do
+        {
+          wcoord[0] += *s / dsize[1];
+          wcoord[1] += *s % dsize[1];
+        }
+      while(++s<sf);
+      brightness = indexs->size;
+    }
+
+
+  /* Convert floatint point coordinates to FITS integers. */
+  coord[0] = GAL_DIMENSION_FLT_TO_INT(wcoord[0]);
+  coord[1] = GAL_DIMENSION_FLT_TO_INT(wcoord[1]);
+
+
+  /* Find the growth limit  */
+  cltprm->std = std[ gal_tile_full_id_from_coord(&p->cp.tl, coord) ];
+  glimit = p->gthresh * cltprm->std;
+
+
+  /* Allocate space to keep the diffuse indexs over this detection. We need
+     to keep the actual indexs since it is our only connection to the
+     object at this stage: we are also going to re-label the pixels to
+     grow. For most astronomical objects, the major part of the detection
+     area is going to be diffuse flux, so we will just allocate the same
+     size as `indexs' array (the `dsize' will be corrected after getting
+     the exact number. */
+  cltprm->diffuseindexs=gal_data_alloc(NULL, GAL_TYPE_SIZE_T, 1,
+                                       cltprm->indexs->dsize, NULL, 0,
+                                       p->cp.minmapsize, NULL, NULL, NULL);
+  dindexs=cltprm->diffuseindexs->array;
+  sf=(s=indexs->array)+indexs->size;
+  do
+    {
+      if( clabel[*s] ) olabel[*s] = clabel[*s];
+      else
+        {
+          olabel[*s] = CLUMPS_INIT;
+          if( imgss[*s]>glimit ) dindexs[ ndiffuse++ ] = *s;
+        }
+    }
+  while(++s<sf);
+
+
+  /* When a user just wants clumps and doesn't want multiple clumps to be
+     considered part of one object, they are going to set an impossibly
+     high `--gthresh' so no growth actually happens. In that case, we don't
+     need the allocated array of indexs. So just free it. */
+  if(ndiffuse==0)
+    {
+      free(cltprm->diffuseindexs->array);
+      cltprm->diffuseindexs->array=NULL;
+    }
+
+
+  /* Correct the sizes of the `diffuseindexs' data structure. */
+  cltprm->diffuseindexs->size = cltprm->diffuseindexs->dsize[0] = ndiffuse;
+}
+
+
+
+
+
+/* Add all the remaining pixels in the detection (below the growth
+   threshold, or those that were not touching). Note that initially
+   `diffuseindexs' was filled with the pixels that are above the growth
+   threshold. That was necessary for identifying the objects. Now that we
+   have identified the objects and labeled them, we want to add the
+   remaining diffuse pixels to it too before doing the final growth.
+
+   Note that the most efficient way is just to re-fill the `diffuseindexs'
+   array instead of adding the pixels below the threshold and sorting them
+   afterwards.*/
+void
+clumps_grow_prepare_final(struct clumps_thread_params *cltprm)
+{
+  size_t ndiffuse=0;
+  size_t *dindexs=cltprm->diffuseindexs->array;
+  uint32_t *olabel=cltprm->clprm->p->olabel->array;
+  size_t *s=cltprm->indexs->array, *sf=s+cltprm->indexs->size;
+
+  /* Recall that we initially allocated `diffuseindexs' to have the same
+     size as the indexs. So there is no problem if there are more pixels in
+     this final round compared to the initial round. */
+  do
+    if( olabel[*s] > CLUMPS_MAXLAB )
+      dindexs[ ndiffuse++ ] = *s;
+  while(++s<sf);
+
+  /* Correct the sizes of the `diffuseindexs' data structure. */
+  cltprm->diffuseindexs->size = cltprm->diffuseindexs->dsize[0] = ndiffuse;
+}
+
+
+
+
+
+/* Grow the true clumps over the diffuse regions of a detection. Note that
+   unlike before, were river pixels would get a separate label for them
+   selves, here, they don't, they just get set back to SEGMENTINIT. This is
+   because some of the pixels that lie immediately between two labeled
+   regions might not be in the blankinds array (they were below the
+   threshold). So we have to find river pixels later on, after the growth
+   is done independently.
+
+   This function is going to be used before identifying objects and also
+   after it (to completely fill in the diffuse area). The distinguishing
+   point between these two steps is the presence of rivers, so you can use
+   the `withrivers' argument. */
+void
+clumps_grow(struct clumps_thread_params *cltprm, int withrivers)
+{
+  struct noisechiselparams *p=cltprm->clprm->p;
+  gal_data_t *diffuseindexs=cltprm->diffuseindexs;
+  size_t ndim=p->input->ndim, *dsize=p->input->dsize;
+
+  int stopngbsearch;
+  size_t *diarray=cltprm->diffuseindexs->array;
+  uint32_t n1, nlab, *olabel=p->olabel->array;
+  size_t *dinc=gal_dimension_increment(ndim, dsize);
+  size_t *s, *sf, thisround, ndiffuse=diffuseindexs->size;
+
+  /* The basic idea is this: after growing, not all the blank pixels are
+     necessarily filled, for example the pixels might belong to two regions
+     above the growth threshold. So the pixels in between them (which are
+     below the threshold will not ever be able to get a label). Therefore,
+     the safest way we can terminate the loop of growing the objects is to
+     stop it when the number of pixels left to fill in this round
+     (thisround) equals the number of blanks.
+
+     To start the loop, we set `thisround' to one more than the number of
+     diffuse pixels. Note that it will be corrected immediately after the
+     loop has started, it is just important to pass the `while'. */
+  thisround=ndiffuse+1;
+  while( thisround > ndiffuse )
+    {
+      /* `thisround' will keep the number of pixels to be inspected in this
+         round. `ndiffuse' will count the number of pixels left without a
+         label by the end of this round. Since `ndiffuse' comes from the
+         previous loop (or outside, for the first round) it has to be saved
+         in `thisround' to begin counting a fresh. */
+      thisround=ndiffuse;
+      ndiffuse=0;
+
+      /* Go over all the available indexs. */
+      sf=(s=diffuseindexs->array)+diffuseindexs->size;
+      do
+        {
+          /* We'll begin by assuming the nearest neighbor of this pixel has
+             no label (has a value of 0). */
+          n1=0;
+
+          /* Check the very closest neighbors of each pixel (4-connectivity
+             in a 2D image). Note that since this macro has multiple loops
+             within it, we can't use break. We'll use a variable instead. */
+          stopngbsearch=0;
+          GAL_DIMENSION_NEIGHBOR_OP(*s, ndim, dsize, 1, dinc,
+            {
+              if(!stopngbsearch)
+                {
+                  nlab = olabel[nind];
+                  if(nlab && nlab<CLUMPS_MAXLAB)  /* This is a real label. */
+                    {
+                      if(n1)  /* A prev. neighboring label has been found. */
+                        {
+                          if( n1 != nlab )       /* Different label from   */
+                            {  /* prevously found neighbor for this pixel. */
+                              n1=CLUMPS_RIVER;
+                              stopngbsearch=1;
+                            }
+                        }
+                      else
+                        {      /* This is the first labeld neighbor found. */
+                          n1=nlab;
+
+                          /* If we want to completely fill in the region
+                             (`withrivers==0'), then there is no point in
+                             looking in other neighbors, the first neighbor
+                             we find is the one we'll use. */
+                          if(!withrivers) stopngbsearch=1;
+                        }
+                    }
+                }
+            } );
+
+          /* The loop above over neighbors finishes with three
+             possibilities:
+
+               n1==0                    --> No labeled neighbor was found.
+               n1==CLUMPS_RIVER         --> Connecting two labeled regions.
+               n1>0 && n1<CLUMPS_MAXLAB --> Only has one neighbouring label.
+
+             The first one means that no neighbors was found and this pixel
+             should be kept for the next loop (we'll be growing the objects
+             pixel-layer by pixel-layer). In the other two cases, we just
+             need to write in the value of `n1'. */
+          if(n1)
+            {
+              olabel[*s]=n1;
+              if(n1>CLUMPS_MAXLAB)             /* We want to keep river   */
+                {
+                  if(withrivers==0) printf("%zu\n", *s);
+                diarray[ ndiffuse++ ] = *s;    /* pixels the diffuse list.*/
+                }
+            }
+          else
+            diarray[ ndiffuse++ ] = *s;
+
+          /* Correct the size of the `diffuseindexs' dataset. */
+          diffuseindexs->size = diffuseindexs->dsize[0] = ndiffuse;
+        }
+      while(++s<sf);
+    }
+
+  /* Clean up. */
+  free(dinc);
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
 /**********************************************************************/
 /*****************             S/N threshold          *****************/
 /**********************************************************************/
@@ -478,7 +729,7 @@ clumps_get_raw_info(struct clumps_thread_params *cltprm)
   /* Do the final preparations. All the calculations are only necessary for
      the clumps that satisfy the minimum area. So there is no need to waste
      time on the smaller ones. */
-  for(lab=1; lab<=cltprm->numinitial; ++lab)
+  for(lab=1; lab<=cltprm->numinitclumps; ++lab)
     {
       row = &info [ lab * INFO_NCOLS ];
       if ( row[INFO_INAREA] > p->segsnminarea )
@@ -517,7 +768,7 @@ void
 clumps_make_sn_table(struct clumps_thread_params *cltprm)
 {
   struct noisechiselparams *p=cltprm->clprm->p;
-  size_t tablen=cltprm->numinitial+1;
+  size_t tablen=cltprm->numinitclumps+1;
 
   float *snarr;
   uint32_t *indarr=NULL;
@@ -624,7 +875,7 @@ clumps_correct_sky_labels_for_check(struct 
clumps_thread_params *cltprm,
                                     gal_data_t *tile)
 {
   gal_data_t *newinds;
-  size_t len=cltprm->numinitial+1;
+  size_t len=cltprm->numinitclumps+1;
   struct noisechiselparams *p=cltprm->clprm->p;
   uint32_t *ninds, curlab, *l=cltprm->snind->array, *lf=l+cltprm->snind->size;
 
@@ -812,82 +1063,6 @@ 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);
-
-  /* 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.
 
@@ -907,11 +1082,13 @@ void
 clumps_true_find_sn_thresh(struct noisechiselparams *p)
 {
   char *msg;
+  uint8_t *b;
+  uint32_t *l, *lf;
   struct timeval t1;
   size_t i, j, c, numsn=0;
   struct clumps_params clprm;
   struct gal_linkedlist_stll *comments=NULL;
-  gal_data_t *claborig, *sn, *snind, *quant;
+  gal_data_t *sn, *tmp, *snind, *quant, *claborig;
 
   /* Get starting time for later reporting if necessary. */
   if(!p->cp.quiet) gettimeofday(&t1, NULL);
@@ -968,8 +1145,8 @@ clumps_true_find_sn_thresh(struct noisechiselparams *p)
           /* Set the extension name. */
           switch(clprm.step)
             {
-            case 1: p->clabel->name = "CLUMPS_ALL_SKY";  break;
-            case 2: p->clabel->name = "CLUMPS_FOR_SN";   break;
+            case 1: p->clabel->name = "SKY_CLUMPS_ALL";  break;
+            case 2: p->clabel->name = "SKY_CLUMPS_FOR_SN";   break;
             default:
               error(EXIT_FAILURE, 0, "a bug! the value %d is not recognized "
                     "in `clumps_true_find_sn_thresh'. Please contact us at "
@@ -977,9 +1154,17 @@ clumps_true_find_sn_thresh(struct noisechiselparams *p)
                     PACKAGE_BUGREPORT);
             }
 
-          /* Write the temporary array into the check image. */
-          gal_fits_img_write(p->clabel, p->segmentationname, NULL,
-                             PROGRAM_STRING);
+          /* Write the demonstration array into the check image. The
+             default values are hard to view, so we'll make a copy of the
+             demo, set all Sky regions to blank and all clump macro values
+             to zero. */
+          tmp=gal_data_copy(p->clabel);
+          b=p->binary->array; lf=(l=tmp->array)+tmp->size;
+          do *l = ( *b++
+                    ? GAL_BLANK_UINT32
+                    : ( *l > CLUMPS_MAXLAB ? 0 : *l ) );  while(++l<lf);
+          gal_fits_img_write(tmp, p->segmentationname, NULL, PROGRAM_STRING);
+          gal_data_free(tmp);
 
           /* Increment the step counter. */
           ++clprm.step;
@@ -1041,8 +1226,9 @@ clumps_true_find_sn_thresh(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_FOR_SN' HDU "
-                                 "of output with `--checksegmentation'.", 1);
+      gal_linkedlist_add_to_stll(&comments, "See also: `SKY_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);
@@ -1061,6 +1247,7 @@ clumps_true_find_sn_thresh(struct noisechiselparams *p)
       free(msg);
     }
 
+
   /* Clean up. */
   gal_data_free(sn);
   gal_data_free(snind);
@@ -1068,3 +1255,116 @@ clumps_true_find_sn_thresh(struct noisechiselparams *p)
   gal_data_array_free(clprm.sn, p->ltl.tottiles, 1);
   gal_data_array_free(clprm.snind, p->ltl.tottiles, 1);
 }
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/***********************************************************************/
+/*****************           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_det_label_indexs(struct noisechiselparams *p)
+{
+  size_t i, *areas;
+  uint32_t *a, *l, *lf;
+  gal_data_t *labindexs=gal_data_array_calloc(p->numdetections+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->numdetections+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->numdetections+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->numdetections+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);
+
+  /* Clean up and return. */
+  free(areas);
+  return labindexs;
+}
+
+
+
+
+
+/* Only keep true clumps over detections. */
+void
+clumps_det_keep_true_relabel(struct clumps_thread_params *cltprm)
+{
+  struct noisechiselparams *p=cltprm->clprm->p;
+  size_t ndim=p->input->ndim, *dsize=p->input->dsize;
+
+  int istouching;
+  size_t *s, *sf;
+  float *sn=cltprm->sn->array;
+  uint32_t *clabel=p->clabel->array;
+  size_t i, *dinc=gal_dimension_increment(ndim, dsize);
+  uint32_t curlab=1, *newlabs=gal_data_calloc_array(GAL_TYPE_UINT32,
+                                                    cltprm->numinitclumps+1);
+
+  /* Set the new labels. Here we will also be removing clumps with a peak
+     that touches a river pixel. */
+  if(p->keepmaxnearriver)
+    {
+      for(i=1;i<cltprm->numinitclumps+1;++i)
+        if( sn[i] > p->clumpsnthresh ) newlabs[i]=curlab++;
+    }
+  else
+    {
+      for(i=1;i<cltprm->numinitclumps+1;++i)
+        {
+          /* Check if all the neighbors of this top element are touching a
+             river or not. */
+          istouching=0;
+          GAL_DIMENSION_NEIGHBOR_OP(cltprm->topinds[i], ndim, dsize, ndim,
+                                    dinc,
+                                    { if(clabel[nind]==0) istouching=1; });
+
+          /* If the peak isn't touching a river, then check its S/N and if
+             that is also good, give it a new label. */
+          if( !istouching && sn[i] > p->clumpsnthresh ) newlabs[i]=curlab++;
+        }
+    }
+
+  /* Save the total number of true clumps in this detection. */
+  cltprm->numtrueclumps=curlab-1;
+
+  /* Correct the clump labels. */
+  sf=(s=cltprm->indexs->array)+cltprm->indexs->size;
+  do if(clabel[*s]) clabel[*s]=newlabs[ clabel[*s] ]; while(++s<sf);
+
+  /* Clean up. */
+  free(dinc);
+  free(newlabs);
+}
diff --git a/bin/noisechisel/clumps.h b/bin/noisechisel/clumps.h
index 466eb62..7e17c90 100644
--- a/bin/noisechisel/clumps.h
+++ b/bin/noisechisel/clumps.h
@@ -24,6 +24,13 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #define CLUMPS_H
 
 
+/* Constants for the clump over-segmentation. */
+#define CLUMPS_RIVER     UINT32_MAX-2
+#define CLUMPS_TMPCHECK  UINT32_MAX-3
+#define CLUMPS_INIT      UINT32_MAX-4
+#define CLUMPS_MAXLAB    UINT32_MAX-5   /* Largest clump label (unsigned). */
+
+
 /* Parameters for all threads. */
 struct clumps_params
 {
@@ -39,32 +46,52 @@ struct clumps_params
 
   /* For detections. */
   gal_data_t        *labindexs; /* Array of `gal_data_t' with obj indexs.  */
+  size_t            totobjects; /* Total number of objects at any point.   */
+  size_t             totclumps; /* Total number of clumps at any point.    */
 };
 
 
-/* Parameters for one thread. */
+/* Parameters for one thread (a tile or a detected region). */
 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. */
+  size_t         numinitclumps; /* Number of clumps in tile/detection.     */
+  size_t         numtrueclumps; /* Number of true clumps in tile/detection.*/
+  size_t            numobjects; /* Number of objects over this clump.      */
+  float                    std; /* Standard deviation of noise on center.  */
+  gal_data_t           *indexs; /* Array containing indexs of this det.    */
+  gal_data_t    *diffuseindexs; /* Diffuse region (after finding clumps).  */
   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.          */
+  gal_data_t       *clumptoobj; /* Index of object that a clump belongs to.*/
   struct clumps_params  *clprm; /* Pointer to main structure.              */
 };
 
+
 void
 clumps_oversegment(struct clumps_thread_params *cltprm);
 
 void
+clumps_grow_prepare_initial(struct clumps_thread_params *cltprm);
+
+void
+clumps_grow_prepare_final(struct clumps_thread_params *cltprm);
+
+void
+clumps_grow(struct clumps_thread_params *cltprm, int withrivers);
+
+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);
+clumps_det_label_indexs(struct noisechiselparams *p);
+
+void
+clumps_det_keep_true_relabel(struct clumps_thread_params *cltprm);
 
 #endif
diff --git a/bin/noisechisel/detection.c b/bin/noisechisel/detection.c
index fef2c5b..1aebc48 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->numinitdets=gal_binary_connected_components(p->binary, &p->olabel, 1);
+  p->numinitialdets=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->numinitdets);
+      asprintf(&msg, "%zu initial detections found.", p->numinitialdets);
       gal_timing_report(&t0, msg, 1);
       free(msg);
     }
@@ -705,7 +705,7 @@ detection_remove_false_initial(struct noisechiselparams *p,
   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->numinitdets+1);
+                                            p->numinitialdets+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->numinitdets;++i) if(newlabels[i]) newlabels[i] = curlab++;
+  for(i=0;i<p->numinitialdets;++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->numinitdets - num_true_initial);
+               p->numinitialdets - num_true_initial);
       gal_timing_report(&t1, msg, 2);
       free(msg);
     }
@@ -827,14 +827,14 @@ detection(struct noisechiselparams *p)
   if(p->dilate)
     {
       gal_binary_dilate(workbin, p->dilate, 8, 1);
-      p->numinitdets = gal_binary_connected_components(workbin, &p->olabel,
+      p->numdetections = gal_binary_connected_components(workbin, &p->olabel,
                                                        8);
     }
-  else p->numinitdets=num_true_initial;
+  else p->numdetections=num_true_initial;
   if(!p->cp.quiet)
     {
       asprintf(&msg, "%zu detections after %zu dilation%s",
-              p->numinitdets, p->dilate, p->dilate>1 ? "s." : ".");
+              p->numdetections, 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 d1f67c2..367048f 100644
--- a/bin/noisechisel/main.h
+++ b/bin/noisechisel/main.h
@@ -115,8 +115,10 @@ struct noisechiselparams
   float                maxstd;  /* Maximum STD before interpolation.      */
   float               cpscorr;  /* Counts/second correction.              */
 
-  size_t          numinitdets;  /* Number of objects detected.            */
+  size_t       numinitialdets;  /* Number of initial detections.          */
+  size_t        numdetections;  /* Number of final detections.            */
   size_t            numclumps;  /* Number of true clumps.                 */
+  size_t           numobjects;  /* Number of objects.                     */
   float           detsnthresh;  /* Pseudo-detection S/N threshold.        */
   float         clumpsnthresh;  /* Clump S/N threshold.                   */
 };
diff --git a/bin/noisechisel/segmentation.c b/bin/noisechisel/segmentation.c
index a800ae8..6856742 100644
--- a/bin/noisechisel/segmentation.c
+++ b/bin/noisechisel/segmentation.c
@@ -30,7 +30,9 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 #include <gnuastro/fits.h>
 #include <gnuastro/blank.h>
+#include <gnuastro/binary.h>
 #include <gnuastro/threads.h>
+#include <gnuastro/dimension.h>
 
 #include <gnuastro-internal/timing.h>
 
@@ -43,6 +45,278 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 
 
+/***********************************************************************/
+/*****************      Relabeling (grown) clumps      *****************/
+/***********************************************************************/
+/* Correct the label of an detection when it doesn't need segmentation (it
+   is fully one object). The final labels for the object(s) with a detected
+   region will be set later (don't forget that we have detections that are
+   composed of multiple objects). So the labels within each detection start
+   from 1.*/
+static void
+segmentation_relab_noseg(struct clumps_thread_params *cltprm)
+{
+  uint32_t *olabel=cltprm->clprm->p->olabel->array;
+  size_t *s=cltprm->indexs->array, *sf=s+cltprm->indexs->size;
+  do olabel[ *s ] = 1; while(++s<sf);
+}
+
+
+
+
+
+/* Find the adjacency matrixs (number, sum and signal to noise) for the
+   rivers between potentially separate objects in a detection region. They
+   have to be allocated prior to entering this funciton.
+
+   The way to find connected objects is through an adjacency matrix. It is
+   a square matrix with a side equal to numobjs. So to see if regions `a`
+   and `b` are connected. All we have to do is to look at element
+   a*numobjs+b or b*numobjs+a and get the answer. Since the number of
+   objects in a given region will not be too high, this is efficient. */
+static void
+segmentation_relab_to_objects(struct clumps_thread_params *cltprm)
+{
+  size_t amwidth=cltprm->numtrueclumps+1;
+  struct noisechiselparams *p=cltprm->clprm->p;
+  size_t ndim=p->input->ndim, *dsize=p->input->dsize;
+
+  size_t mdsize[2]={amwidth, amwidth};
+  gal_data_t *nums_d=gal_data_alloc(NULL, GAL_TYPE_SIZE_T, 2, mdsize, NULL,
+                                    1, p->cp.minmapsize, NULL, NULL, NULL);
+  gal_data_t *sums_d=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 2, mdsize, NULL,
+                                    1, p->cp.minmapsize, NULL, NULL, NULL);
+  gal_data_t *adjacency_d=gal_data_alloc(NULL, GAL_TYPE_UINT8, 2, mdsize,
+                                         NULL, 1, p->cp.minmapsize, NULL,
+                                         NULL, NULL);
+  float *imgss=p->input->array;
+  uint8_t *adjacency=adjacency_d->array;
+  size_t nngb=gal_dimension_num_neighbors(ndim);
+  uint32_t *clumptoobj, *olabel=p->olabel->array;
+  size_t *dinc=gal_dimension_increment(ndim, dsize);
+  size_t *s, *sf, i, j, ii, rpnum, *nums=nums_d->array;
+  double ave, rpsum, c=sqrt(1/p->cpscorr), *sums=sums_d->array;
+  uint32_t *ngblabs=gal_data_malloc_array(GAL_TYPE_UINT32, nngb);
+  double err=cltprm->std*cltprm->std*(p->skysubtracted?1.0f:2.0f);
+
+
+  /* Go over all the still-unlabeled pixels and see which labels they
+     touch. In the process, get the average value of the river-pixel values
+     and put them in the respective adjacency matrix. Note that at this
+     point, the rivers are also part of the "diffuse" regions. So we don't
+     need to go over all the indexs of this object, only its diffuse
+     indexs. */
+  sf=(s=cltprm->diffuseindexs->array)+cltprm->diffuseindexs->size;
+  do
+    /* We only want to work on pixels that have already been identified as
+       touching more than one label: river pixels. */
+    if( olabel[ *s ]==CLUMPS_RIVER )
+      {
+        /* Initialize the values. */
+        ii=0;
+        rpnum=1;              /* River-pixel number of points used. */
+        rpsum=imgss[*s];      /* River-pixel sum of values used.    */
+        memset(ngblabs, 0, nngb*sizeof *ngblabs);
+
+        /* Check all the fully-connected neighbors of this pixel and see if
+           it touches a label or not */
+        GAL_DIMENSION_NEIGHBOR_OP(*s, ndim, dsize, ndim, dinc, {
+            if( olabel[nind] && olabel[nind] < CLUMPS_MAXLAB )
+              {
+                /* Add this neighbor's value and increment the number. */
+                if( !isnan(imgss[nind]) ) { ++rpnum; rpsum+=imgss[nind]; }
+
+                /* Go over the already found neighbors and see if this grown
+                   clump has already been considered or not. */
+                for(i=0;i<ii;++i) if(ngblabs[i]==olabel[nind]) break;
+
+                /* This is the first time we are getting to this neighbor: */
+                if(i==ii) ngblabs[ ii++ ] = olabel[nind];
+              }
+          } );
+
+        /* For a check:
+        if(ii>0)
+          {
+            printf("%zu, %zu:\n", *s%dsize[1]+1, *s/dsize[1]+1);
+            for(i=0;i<ii;++i) printf("\t%u\n", ngblabs[i]);
+          }
+        */
+
+        /* If more than one neighboring label was found, fill in the 'sums'
+           and 'nums' adjacency matrixs with the values for this
+           pixel. Recall that ii is the number of neighboring labels to
+           this river pixel. */
+        if(ii>i)
+          for(i=0;i<ii;++i)
+            for(j=0;j<ii;++j)
+              if(i!=j)
+                {
+                  /* For safety, we will fill both sides of the diagonal. */
+                  ++nums[ ngblabs[i] * amwidth + ngblabs[j] ];
+                  ++nums[ ngblabs[j] * amwidth + ngblabs[i] ];
+                  sums[ ngblabs[i] * amwidth + ngblabs[j] ] += rpsum/rpnum;
+                  sums[ ngblabs[j] * amwidth + ngblabs[i] ] += rpsum/rpnum;
+                }
+      }
+  while(++s<sf);
+
+
+  /* We now have the average values and number of all rivers between the
+     grown clumps. We now want to finalize their connection (given the
+     user's criteria). */
+  for(i=1;i<amwidth;++i)
+    for(j=1;j<i;++j)
+      {
+        ii = i * amwidth + j;
+        if(nums[ii]>p->minriverlength)       /* There is a connection. */
+          {
+            /* For easy reading. */
+            ave=sums[ii]/nums[ii];
+
+            /* In case the average is negative (only possible if `sums' is
+               negative), don't change the adjacency: it is already
+               initialized to zero. Note that even an area of 1 is
+               acceptable, and we put no area criteria here, because the
+               fact that a river exists between two clumps is important. */
+            if( ave>0.0f && ( c * ave / sqrt(ave+err) ) > p->objbordersn )
+              {
+                adjacency[ii]=1;       /* We want to set both sides of the */
+                adjacency[ j * amwidth + i ] = 1; /* matrix diagonal to 1. */
+              }
+          }
+      }
+
+
+  /* For a check:
+  printf("=====================\n");
+  printf("%zu:\n--------\n", cltprm->id);
+  for(i=1;i<amwidth;++i)
+    {
+      printf(" %zu...\n", i);
+      for(j=1;j<amwidth;++j)
+        {
+          ii=i*amwidth+j;
+          if(nums[ii])
+            {
+              ave=sums[ii]/nums[ii];
+              printf("    ...%zu: N:%-4zu S:%-10.2f S/N: %-10.2f "
+                     "--> %u\n", j, nums[ii], sums[ii], c*ave/sqrt(ave+err),
+                       adjacency[ii]);
+            }
+        }
+      printf("\n");
+    }
+  */
+
+
+  /* Calculate the new labels for each grown clump. */
+  cltprm->clumptoobj = gal_binary_connected_adjacency_matrix(adjacency_d,
+                                                      &cltprm->numobjects);
+  clumptoobj = cltprm->clumptoobj->array;
+
+
+  /* Correct all the labels. */
+  sf=(s=cltprm->indexs->array)+cltprm->indexs->size;
+  do
+    if( olabel[*s]<CLUMPS_MAXLAB )
+      olabel[*s] = clumptoobj[ olabel[*s] ];
+  while(++s<sf);
+
+
+  /* Clean up and return. */
+  free(dinc);
+  free(ngblabs);
+  gal_data_free(nums_d);
+  gal_data_free(sums_d);
+  gal_data_free(adjacency_d);
+}
+
+
+
+
+/* The correspondance between the clumps and objects has been found. With
+   this function, we want to correct the clump labels so the clump IDs in
+   each object start from 1 and are contiguous. */
+static void
+segmentation_relab_clumps_in_objects(struct clumps_thread_params *cltprm)
+{
+  size_t numobjects=cltprm->numobjects, numtrueclumps=cltprm->numtrueclumps;
+
+  uint32_t *clumptoobj=cltprm->clumptoobj->array;
+  uint32_t *clabel=cltprm->clprm->p->clabel->array;
+  size_t i, *s=cltprm->indexs->array, *sf=s+cltprm->indexs->size;
+  size_t *nclumpsinobj=gal_data_calloc_array(GAL_TYPE_SIZE_T, numobjects+1);
+  uint32_t *newlabs=gal_data_calloc_array(GAL_TYPE_UINT32, numtrueclumps+1);
+
+  /* Fill both arrays. */
+  for(i=1;i<numtrueclumps+1;++i)
+    newlabs[i] = ++nclumpsinobj[ clumptoobj[i] ];
+
+  /* Reset the clump labels over the detection region. */
+  do if(clabel[*s]) clabel[*s] = newlabs[ clabel[*s] ]; while(++s<sf);
+
+  /* Clean up. */
+  free(newlabs);
+  free(nclumpsinobj);
+}
+
+
+
+
+
+/* Prior to this function, the objects have labels that are unique and
+   contiguous (the labels are contiguous, not the objects!) within each
+   detection and start from 1. However, for the final output, it is
+   necessary that each object over the whole dataset have a unique
+   ID. Since multiple threads are working on separate objects at every
+   instance, this function will use a mutex to limit the reading and
+   writing to the variable keeping the total number of objects counter. */
+static void
+segmentation_relab_overall(struct clumps_thread_params *cltprm)
+{
+  struct clumps_params *clprm=cltprm->clprm;
+  uint32_t startinglab, *olabel=clprm->p->olabel->array;
+  size_t *s=cltprm->indexs->array, *sf=s+cltprm->indexs->size;
+
+  /* Lock the mutex if we are working on more than one thread. NOTE: it is
+     very important to keep the number of operations within the mutex to a
+     minimum so other threads don't get delayed. */
+  if(clprm->p->cp.numthreads>1)
+    pthread_mutex_lock(&clprm->labmutex);
+
+  /* Save the total number of objects found so far into `startinglab', then
+     increment the total number of objects and clumps. */
+  startinglab        = clprm->totobjects;
+  clprm->totobjects += cltprm->numobjects;
+  clprm->totclumps  += cltprm->numtrueclumps;
+
+  /* Unlock it (if it was locked). */
+  if(clprm->p->cp.numthreads>1)
+    pthread_mutex_destroy(&clprm->labmutex);
+
+  /* Increase all the object labels by `startinglab'. */
+  do olabel[*s] += startinglab; while(++s<sf);
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
 
 /***********************************************************************/
 /*****************            Over detections          *****************/
@@ -55,9 +329,10 @@ segmentation_on_threads(void *in_prm)
   struct clumps_params *clprm=(struct clumps_params *)(tprm->params);
   struct noisechiselparams *p=clprm->p;
 
-  size_t i;
+  size_t i, *s, *sf;
   gal_data_t *topinds;
   struct clumps_thread_params cltprm;
+  uint32_t *clabel=p->clabel->array, *olabel=p->olabel->array;
 
   /* Initialize the general parameters for this thread. */
   cltprm.clprm   = clprm;
@@ -65,7 +340,6 @@ segmentation_on_threads(void *in_prm)
   /* 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. */
@@ -95,8 +369,15 @@ segmentation_on_threads(void *in_prm)
       /* 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);
+         loop at that point.
+
+         Note that the array of `gal_data_t' that keeps the S/N table for
+         each detection is allocated before threading starts. However, when
+         the user wants to inspect the steps, this function is called
+         multiple times. So we need to avoid over-writing the allocations. */
+      if( clprm->sn[ cltprm.id ].dsize==NULL )
+        clumps_make_sn_table(&cltprm);
+      else cltprm.sn=&clprm->sn[ cltprm.id ];
 
       /* 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
@@ -104,12 +385,103 @@ segmentation_on_threads(void *in_prm)
       if(clprm->step==1 || p->checkclumpsn)
         { gal_data_free(topinds); continue; }
 
+      /* Only keep true clumps. */
+      clumps_det_keep_true_relabel(&cltprm);
+      gal_data_free(topinds);
+      if(clprm->step==2) continue;
+
+      /* Set the internal (with the detection) clump and object
+         labels. Segmenting a detection into multiple objects is only
+         defined when there is more than one true clump over the
+         detection. When there is only one true clump
+         (cltprm->numtrueclumps==1) or none (p->numtrueclumps==0), then
+         just set the required preliminaries to make the next steps be
+         generic for all cases. */
+      if(cltprm.numtrueclumps<=1)
+        {
+          /* Set the basics. */
+          cltprm.numobjects=1;
+          segmentation_relab_noseg(&cltprm);
 
 
-      /* Clean up. */
-      gal_data_free(topinds);
+          /* If the user wanted a check image, this object doesn't
+             change. */
+          if( clprm->step >= 3 && clprm->step <= 6) continue;
+
+
+          /* If the user has asked for grown clumps in the clumps image
+             instead of the raw clumps, then replace the indexs in the
+             `clabel' array is well. */
+          if(p->grownclumps)
+            {
+              sf=(s=cltprm.indexs->array)+cltprm.indexs->size;
+              do clabel[ *s++ ] = 1; while(s<sf);
+            }
+        }
+      else
+        {
+          /* Grow the true clumps over the detection. */
+          clumps_grow_prepare_initial(&cltprm);
+          if(cltprm.diffuseindexs->size) clumps_grow(&cltprm, 1);
+          if(clprm->step==3)
+            { gal_data_free(cltprm.diffuseindexs); continue; }
+
+
+          /* If grown clumps are desired instead of the raw clumps, then
+             replace all the grown clumps with those in clabel. */
+          if(p->grownclumps)
+            {
+              sf=(s=cltprm.indexs->array)+cltprm.indexs->size;
+              do
+                if(olabel[*s]<CLUMPS_MAXLAB) clabel[*s]=olabel[*s];
+              while(++s<sf);
+            }
+
+
+          /* Identify the objects within the grown clumps and correct the
+             grown clump labels into new object labels. */
+          segmentation_relab_to_objects(&cltprm);
+          if(clprm->step==4)
+            {
+              gal_data_free(cltprm.clumptoobj);
+              gal_data_free(cltprm.diffuseindexs);
+              continue;
+            }
+
+
+          /* Continue the growth and cover the whole area, we don't need
+             the diffuse indexs any more, so after filling the detected
+             region, free the indexs. */
+          if( cltprm.numobjects == 1 )
+            segmentation_relab_noseg(&cltprm);
+          else
+            {
+              /* Correct the labels so every non-labeled pixel can be
+                 grown. */
+              clumps_grow_prepare_final(&cltprm);
+
+              /* Cover the whole area. */
+              clumps_grow(&cltprm, 0);
+            }
+          gal_data_free(cltprm.diffuseindexs);
+          if(clprm->step==5) { gal_data_free(cltprm.clumptoobj); continue; }
+
+
+          /* Correct the clump labels. Note that this is only necessary
+             when there is more than object over the detection. When there
+             is only one object over the full detection or if there is only
+             one clump, the existing clump labels are fine.  */
+          if(cltprm.numobjects>1)
+            segmentation_relab_clumps_in_objects(&cltprm);
+          gal_data_free(cltprm.clumptoobj);
+          if(clprm->step==6) continue;
+        }
+
+      /* Convert the object labels to their final value */
+      segmentation_relab_overall(&cltprm);
     }
 
+  /* Wait until all the threads finish then return. */
   if(tprm->b) pthread_barrier_wait(tprm->b);
   return NULL;
 }
@@ -135,7 +507,7 @@ segmentation_save_sn_table(struct clumps_params *clprm)
   /* 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;
+  for(i=1;i<p->numdetections+1;++i) totclumps += clprm->sn[i].size-1;
 
 
   /* Allocate the columns for the table. */
@@ -154,7 +526,7 @@ segmentation_save_sn_table(struct clumps_params *clprm)
   sarr=sn->array;
   oiarr=objind->array;
   cioarr=clumpinobj->array;
-  for(i=1;i<p->numinitdets+1;++i)
+  for(i=1;i<p->numdetections+1;++i)
     for(j=1;j<clprm->sn[i].size;++j)
       {
         oiarr[c]  = i;
@@ -200,24 +572,33 @@ segmentation_save_sn_table(struct clumps_params *clprm)
 
 
 /* Find true clumps over the detected regions. */
-void
+static void
 segmentation_detections(struct noisechiselparams *p)
 {
+  char *msg;
+  uint8_t *b;
+  uint32_t *l, *lf;
   struct clumps_params clprm;
-  gal_data_t *labindexs, *claborig;
+  gal_data_t *tmp, *demo, *labindexs, *claborig;
 
 
   /* Get the indexs of all the pixels in each label. */
-  labindexs=clumps_label_indexs(p);
+  labindexs=clumps_det_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.totclumps=0;
+  clprm.totobjects=0;
   clprm.snind = NULL;
   clprm.labindexs=labindexs;
-  clprm.sn=gal_data_array_calloc(p->numinitdets+1);
+  clprm.sn=gal_data_array_calloc(p->numdetections+1);
+
+
+  /* When more than one thread is to be used, initialize the mutex. */
+  if( p->cp.numthreads > 1 ) pthread_mutex_init(&clprm.labmutex, NULL);
 
 
   /* Spin off the threads to start the work. Note that several steps are
@@ -233,8 +614,9 @@ segmentation_detections(struct noisechiselparams *p)
       claborig=p->clabel;
       p->clabel=gal_data_copy(claborig);
 
+
       /* Do each step. */
-      while(clprm.step<3)
+      while(clprm.step<8)
         {
           /* Reset the temporary copy of clabel back to its original. */
           if(clprm.step>1)
@@ -243,13 +625,96 @@ segmentation_detections(struct noisechiselparams *p)
 
           /* (Re-)do everything until this step. */
           gal_threads_spin_off(segmentation_on_threads, &clprm,
-                               p->numinitdets, p->cp.numthreads);
+                               p->numdetections, 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;
+            case 1:
+              demo=p->clabel;
+              demo->name = "DET_CLUMPS_ALL";
+              if(!p->cp.quiet)
+                {
+                  asprintf(&msg, "Identified clumps over detections  "
+                           "(HDU: `%s').", demo->name);
+                  gal_timing_report(NULL, msg, 2);
+                  free(msg);
+                }
+              break;
+
+            case 2:
+              demo=p->clabel;
+              demo->name = "DET_CLUMPS_TRUE";
+              if(!p->cp.quiet)
+                {
+                  asprintf(&msg, "True clumps found                  "
+                           "(HDU: `%s').", demo->name);
+                  gal_timing_report(NULL, msg, 2);
+                  free(msg);
+                }
+              break;
+
+            case 3:
+              demo=p->olabel;
+              demo->name = "DET_CLUMPS_GROWN";
+              if(!p->cp.quiet)
+                {
+                  gal_timing_report(NULL, "Starting to identify objects.", 1);
+                  asprintf(&msg, "True clumps grown                  "
+                           "(HDU: `%s').", demo->name);
+                  gal_timing_report(NULL, msg, 2);
+                  free(msg);
+                }
+              break;
+
+            case 4:
+              demo=p->olabel;
+              demo->name = "DET_OBJ_IDENTIFIED";
+              if(!p->cp.quiet)
+                {
+                  asprintf(&msg, "Identified objects over detections "
+                           "(HDU: `%s').", demo->name);
+                  gal_timing_report(NULL, msg, 2);
+                  free(msg);
+                }
+              break;
+
+            case 5:
+              demo=p->olabel;
+              demo->name = "DET_OBJECTS_FULL";
+              if(!p->cp.quiet)
+                {
+                  asprintf(&msg, "Objects grown to cover full area   "
+                           "(HDU: `%s').", demo->name);
+                  gal_timing_report(NULL, msg, 2);
+                  free(msg);
+                }
+              break;
+
+            case 6:
+              demo=p->clabel;
+              demo->name = "CLUMPS_FINAL";
+              if(!p->cp.quiet)
+                {
+                  asprintf(&msg, "Clumps given their final label     "
+                           "(HDU: `%s').", demo->name);
+                  gal_timing_report(NULL, msg, 2);
+                  free(msg);
+                }
+              break;
+
+            case 7:
+              demo=p->olabel;
+              demo->name = "OBJECTS_FINAL";
+              if(!p->cp.quiet)
+                {
+                  asprintf(&msg, "Objects given their final label    "
+                           "(HDU: `%s').", demo->name);
+                  gal_timing_report(NULL, msg, 2);
+                  free(msg);
+                }
+              break;
+
             default:
               error(EXIT_FAILURE, 0, "a bug! the value %d is not recognized "
                     "in `segmentation_detections'. Please contact us at "
@@ -257,9 +722,17 @@ segmentation_detections(struct noisechiselparams *p)
                     PACKAGE_BUGREPORT);
             }
 
-          /* Write the temporary array into the check image. */
-          gal_fits_img_write(p->clabel, p->segmentationname, NULL,
-                             PROGRAM_STRING);
+          /* Write the demonstration array into the check image. The
+             default values are hard to view, so we'll make a copy of the
+             demo, set all Sky regions to blank and all clump macro values
+             to zero. */
+          tmp=gal_data_copy(demo);
+          b=p->binary->array; lf=(l=tmp->array)+tmp->size;
+          do *l = ( *b++
+                    ? ( *l > CLUMPS_MAXLAB ? 0 : *l )
+                    : GAL_BLANK_UINT32 );  while(++l<lf);
+          gal_fits_img_write(tmp, p->segmentationname, NULL, PROGRAM_STRING);
+          gal_data_free(tmp);
 
           /* 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
@@ -273,23 +746,27 @@ segmentation_detections(struct noisechiselparams *p)
 
       /* Clean up (we don't need the original any more). */
       gal_data_free(claborig);
-      p->clabel->name=NULL;
+      p->olabel->name = p->clabel->name = NULL;
     }
   else
     {
       clprm.step=0;
-      gal_threads_spin_off(segmentation_on_threads, &clprm, p->numinitdets,
+      gal_threads_spin_off(segmentation_on_threads, &clprm, p->numdetections,
                            p->cp.numthreads);
     }
 
+  /* Save the final number of objects and clumps. */
+  p->numclumps=clprm.totclumps;
+  p->numobjects=clprm.totobjects;
 
   /* 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);
+  /* Clean up allocated structures and destroy the mutex. */
+  gal_data_array_free(clprm.sn, p->numdetections+1, 1);
+  gal_data_array_free(labindexs, p->numdetections+1, 1);
+  if( p->cp.numthreads>1 ) pthread_mutex_destroy(&clprm.labmutex);
 }
 
 
@@ -318,12 +795,17 @@ void
 segmentation(struct noisechiselparams *p)
 {
   float *f;
+  char *msg;
   uint32_t *l, *lf;
+  struct timeval t1;
 
   /* To keep the user up to date. */
   if(!p->cp.quiet)
-    gal_timing_report(NULL, "Starting over-segmentation (finding clumps).",
-                      1);
+    {
+      if(!p->cp.quiet) gettimeofday(&t1, NULL);
+      gal_timing_report(NULL, "Starting segmentation.",
+                        1);
+    }
 
 
   /* If a check segmentation image was requested, then put in the
@@ -332,8 +814,10 @@ segmentation(struct noisechiselparams *p)
     {
       gal_fits_img_write(p->input, p->segmentationname, NULL, PROGRAM_STRING);
       gal_fits_img_write(p->conv, p->segmentationname, NULL, PROGRAM_STRING);
+      p->olabel->name="DETECTION_LABELS";
       gal_fits_img_write(p->olabel, p->segmentationname, NULL,
                          PROGRAM_STRING);
+      p->olabel->name=NULL;
     }
 
 
@@ -360,6 +844,14 @@ segmentation(struct noisechiselparams *p)
   /* Find true clumps over the detected regions. */
   segmentation_detections(p);
 
+  /* Report the results and timing to the user. */
+  if(!p->cp.quiet)
+    {
+      asprintf(&msg, "%zu object%s""containing %zu clump%sfound.",
+               p->numobjects-1, p->numobjects==2 ? " " : "s ",
+               p->numclumps-1,  p->numclumps ==2 ? " " : "s ");
+      gal_timing_report(&t1, msg, 1);
+    }
 
   /* If the user wanted to check the segmentation and hasn't called
      `continueaftercheck', then stop NoiseChisel. */
diff --git a/bin/noisechisel/ui.c b/bin/noisechisel/ui.c
index dd562aa..25f739c 100644
--- a/bin/noisechisel/ui.c
+++ b/bin/noisechisel/ui.c
@@ -694,7 +694,7 @@ ui_abort_after_check(struct noisechiselparams *p, char 
*filename,
           "------------------------------------------------\n"
           "%s (%s) has been created.\n\n"
           "If you want %s to continue its processing AND save any "
-          "requested check(s), please run it again with "
+          "requested check outputs, please run it again with "
           "`--continueaftercheck'.\n"
           "------------------------------------------------\n",
           PROGRAM_NAME, name, description, PROGRAM_NAME);
diff --git a/lib/binary.c b/lib/binary.c
index a5bd201..f283450 100644
--- a/lib/binary.c
+++ b/lib/binary.c
@@ -328,6 +328,16 @@ gal_binary_open(gal_data_t *input, size_t num, int 
connectivity,
 
 
 
+
+
+
+
+
+
+
+
+
+
 /*********************************************************************/
 /*****************      Connected components      ********************/
 /*********************************************************************/
@@ -446,6 +456,107 @@ gal_binary_connected_components(gal_data_t *binary, 
gal_data_t **out,
 
 
 
+/* Given an adjacency matrix (which should be binary), find the number of
+   connected objects and return an array of new labels for each old
+   label. In other words, this function will find the objects that are
+   connected (possibly through a third object) and in the output array, the
+   respective elements for all these objects is going to have the same
+   value. The total number of connected labels is put into the place
+   pointed to by `numconnected'.
+
+   Labels begin from 1 (0 is kept for non-labeled regions usually). So if
+   you have 3 initial objects/labels, the input matrix to this function
+   should have a size of 4x4. The first (label 0) row and column are not
+   going to be parsed/checked.
+
+   The adjacency matrix needs to be completely filled (on both sides of the
+   diagonal) for this function.
+
+   If the input adjacency matrix has a size of `amsize * amsize', the
+   output will have a size of `amsize' with each index having a new label
+   in its place. */
+gal_data_t *
+gal_binary_connected_adjacency_matrix(gal_data_t *adjacency,
+                                      size_t *numconnected)
+{
+  uint32_t *newlabs;
+  gal_data_t *newlabs_d;
+  struct gal_linkedlist_sll *Q=NULL;
+  size_t i, j, p, num=adjacency->dsize[0];
+  uint8_t *adj=adjacency->array, curlab=1;
+
+  /* Some small sanity checks. */
+  if(adjacency->type != GAL_TYPE_UINT8)
+    error(EXIT_FAILURE, 0, "`gal_binary_connected_adjacency_matrix' only "
+          "accepts a `uint8' numeric type. However, the input dataset has "
+          "type of `%s'", gal_type_to_string(adjacency->type, 1));
+
+  if(adjacency->ndim != 2)
+    error(EXIT_FAILURE, 0, "`gal_binary_connected_adjacency_matrix' only "
+          "accepts a 2-dimensional array. However, the input dataset has "
+          "%zu dimensions", adjacency->ndim);
+
+  if(adjacency->dsize[0] != adjacency->dsize[1])
+    error(EXIT_FAILURE, 0, "`gal_binary_connected_adjacency_matrix' only "
+          "accepts a square matrix. However, the input dataset has a size "
+          "of %zu x %zu", adjacency->dsize[0], adjacency->dsize[1]);
+
+
+  /* Allocate (and clear) the output datastructure. */
+  newlabs_d=gal_data_alloc(NULL, GAL_TYPE_UINT32, 1, &num, NULL, 1,
+                           adjacency->minmapsize, NULL, NULL, NULL);
+  newlabs=newlabs_d->array;
+
+
+  /* Go over the input matrix and apply the same principle as we used to
+     identify connected components in an image: through a queue, find those
+     elements that are connected. */
+  for(i=1;i<num;++i)
+    if(newlabs[i]==0)
+      {
+        /* Add this old label to the list that must be corrected. */
+        gal_linkedlist_add_to_sll(&Q, i);
+
+        /* Continue while the list has elements. */
+        while(Q!=NULL)
+          {
+            /* Pop the top old-label from the list. */
+            gal_linkedlist_pop_from_sll(&Q, &p);
+
+            /* If it has already been labeled then ignore it. */
+            if( newlabs[p]!=curlab )
+              {
+                /* Give it the new label. */
+                newlabs[p]=curlab;
+
+                /* Go over the adjacecny matrix row for this touching
+                   object and see if there are any not-yet-labeled objects
+                   that are touching it. */
+                for(j=1;j<num;++j)
+                  if( adj[ p*num+j ] && newlabs[j]==0 )
+                    gal_linkedlist_add_to_sll(&Q, j);
+              }
+          }
+
+        /* Increment the current label. */
+        ++curlab;
+      }
+
+
+  /* For a check.
+  printf("=== Old labels --> new labels ===\n");
+  for(i=1;i<num;++i) printf("%zu: %u\n", i, newlabs[i]);
+  */
+
+  /* Return the output. */
+  *numconnected = curlab-1;
+  return newlabs_d;
+}
+
+
+
+
+
 
 
 
diff --git a/lib/convolve.c b/lib/convolve.c
index ca0a685..556eefe 100644
--- a/lib/convolve.c
+++ b/lib/convolve.c
@@ -536,7 +536,7 @@ gal_convolve_spatial_general(gal_data_t *tiles, gal_data_t 
*kernel,
   else
     {
       name = ( block->name
-               ? gal_checkset_malloc_cat("CONVL_", block->name) : NULL );
+               ? gal_checkset_malloc_cat(block->name, "_CONVOLVED") : NULL );
       out=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, block->ndim, block->dsize,
                          block->wcs, 0, block->minmapsize, name,
                          block->unit, NULL);
diff --git a/lib/gnuastro/binary.h b/lib/gnuastro/binary.h
index b06cb5e..cb19650 100644
--- a/lib/gnuastro/binary.h
+++ b/lib/gnuastro/binary.h
@@ -89,8 +89,9 @@ size_t
 gal_binary_connected_components(gal_data_t *binary, gal_data_t **out,
                                 int connectivity);
 
-
-
+gal_data_t *
+gal_binary_connected_adjacency_matrix(gal_data_t *adjacency,
+                                      size_t *numnewlabs);
 
 
 /*********************************************************************/