[gnuastro-commits] master 133e277 071/113: Imported recent work from master, conflicts fixed

[Mohammad Akhlaghi]

branch: master
commit 133e2779109ab30c8412f8f1400c328c0d772d7b
Merge: 63b4edd 24757bf
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    Imported recent work from master, conflicts fixed
    
    Some minor conflicts came up and were fixed.
---
 NEWS                  |  15 +-
 bin/segment/clumps.c  |   8 +-
 bin/segment/segment.c |   6 +-
 doc/gnuastro.texi     |  29 ++--
 lib/gnuastro/label.h  |   4 +-
 lib/label.c           | 407 ++++++++++++++++++++++++--------------------------
 tests/during-dev.sh   |  15 +-
 7 files changed, 240 insertions(+), 244 deletions(-)

diff --git a/NEWS b/NEWS
index b1b6b06..4dc4f6c 100644
--- a/NEWS
+++ b/NEWS
@@ -87,7 +87,7 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
     gal_jpeg_read: Reads input JPEG image into `gal_data_t'.
     gal_jpeg_write: Writes a `gal_data_t' into a JPEG file.
     gal_label_grow_indexs: grow known indexs into desired areas.
-    gal_label_oversegment: apply over-segmentation to an input dataset.
+    gal_label_watershed: apply watershed algorithm on desired region.
     gal_label_clump_significance: measure significance of all clumps in region.
     gal_pdf_name_is_pdf: Returns 1 if given filename is PDF.
     gal_pdf_suffix_is_pdf: Returns 1 if given suffix is PDF.
@@ -156,11 +156,14 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
       --detquant     ==> --snquant
     - By default the output detection map is a binary image (values of 0 or 1).
     - With no output name, the output has a `_detected.fits' suffix.
-    - [Now in Segment]: For finding true clumps, the difference in the peak
-      of the clump and the highest valued river pixel, divided by the noise
-      standard deviation are used, not the total signal-to-noise ratio. In
-      initial tests, this algorithm was much more promising in detecting
-      clumps over strong gradients.
+
+  Segment:
+    - [Previously in NoiseChisel]: For finding true clumps, the difference
+      in the peak of the clump and the highest valued river pixel, divided
+      by the noise standard deviation are used. Until now, the total
+      signal-to-noise ratio was used as a criteria. In initial tests, this
+      algorithm was much more promising in detecting clumps over strong
+      gradients and also on flatter gradients.
 
   Table:
     --column: multiple columns (comma separated) can be used in one
diff --git a/bin/segment/clumps.c b/bin/segment/clumps.c
index 0a67803..d9bf72d 100644
--- a/bin/segment/clumps.c
+++ b/bin/segment/clumps.c
@@ -444,10 +444,10 @@ clumps_find_make_sn_table(void *in_prm)
 
 
           /* Generate the clumps over this region. */
-          cltprm.numinitclumps=gal_label_oversegment(p->conv, cltprm.indexs,
-                                                     p->clabel,
-                                                     cltprm.topinds,
-                                                     !p->minima);
+          cltprm.numinitclumps=gal_label_watershed(p->conv, cltprm.indexs,
+                                                   p->clabel,
+                                                   cltprm.topinds,
+                                                   !p->minima);
 
 
           /* Set all river pixels to GAL_LABEL_INIT (to be distinguishable
diff --git a/bin/segment/segment.c b/bin/segment/segment.c
index a456c1c..a0dbd6a 100644
--- a/bin/segment/segment.c
+++ b/bin/segment/segment.c
@@ -564,9 +564,9 @@ segment_on_threads(void *in_prm)
 
 
       /* Find the clumps over this region. */
-      cltprm.numinitclumps=gal_label_oversegment(p->conv, cltprm.indexs,
-                                                 p->clabel, cltprm.topinds,
-                                                 !p->minima);
+      cltprm.numinitclumps=gal_label_watershed(p->conv, cltprm.indexs,
+                                               p->clabel, cltprm.topinds,
+                                               !p->minima);
 
 
       /* Set all the river pixels to zero (we don't need them any more in
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index bcb3825..ef203c2 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -15665,20 +15665,23 @@ detect the diffuse and extended emission, but in 
segmentation, you want to
 detect sharp peaks.
 
 @item
-The criteria to select true from false clumps is the peak signal-to-noise
-ratio. This value is calculated from a clump's peak value (@mymath{C_c})
-and the highest valued river pixel around that clump (@mymath{R_c}). Both
-are calculated on the convolved image (signified by the @mymath{c}
-subscript). To avoid absolute differences, it is then divided by the input
-(not convolved) Sky standard deviation under that clump (@mymath{\sigma})
-as shown below.
+The criteria to select true from false clumps is the peak significance. It
+is defined to be the difference between the clump's peak value
+(@mymath{C_c}) and the highest valued river pixel around that clump
+(@mymath{R_c}). Both are calculated on the convolved image (signified by
+the @mymath{c} subscript). To avoid absolute values (differing from dataset
+to dataset), @mymath{C_c-R_c} is then divided by the Sky standard deviation
+under the river pixel used (@mymath{\sigma_r}) as shown below:
 
-@dispmath{C_c-R_c\over \sigma}
+@dispmath{C_c-R_c\over \sigma_r}
+
+When @option{--minima} is given, the nominator becomes
+@mymath{R_c-C_c}.
 
 The input Sky standard deviation dataset (@option{--std}) is assumed to be
 for the unconvolved image. Therefore a constant factor (related to the
 convolution kernel) is necessary to convert this into an absolute peak
-signal-to-noise ratio@footnote{To get an estimate of the standard deviation
+significance@footnote{To get an estimate of the standard deviation
 correction factor between the input and convolved images, you can take the
 following steps: 1) Mask (set to NaN) all detections on the convolved image
 with the @code{where} operator or @ref{Arithmetic}. 2) Calculate the
@@ -26742,7 +26745,7 @@ dataset's first element/pixel. Therefore it is always 
greater or equal to
 zero and stored in @code{size_t} type.
 @end deftypefun
 
-@deftypefun size_t gal_label_oversegment (gal_data_t @code{*values}, 
gal_data_t @code{*indexs}, gal_data_t @code{*label}, size_t @code{*topinds}, 
int @code{min0_max1})
+@deftypefun size_t gal_label_watershed (gal_data_t @code{*values}, gal_data_t 
@code{*indexs}, gal_data_t @code{*label}, size_t @code{*topinds}, int 
@code{min0_max1})
 @cindex Watershed algorithm
 @cindex Algorithm: watershed
 Use the watershed algorithm@footnote{The watershed algorithm was initially
@@ -26806,7 +26809,7 @@ input->flags &= ~GAL_DATA_FLAG_HASBLANK; /* Set bit to 
0. */
 
 @deftypefun void gal_label_clump_significance (gal_data_t @code{*values}, 
gal_data_t @code{*std}, gal_data_t @code{*label}, gal_data_t @code{*indexs}, 
struct gal_tile_two_layer_params @code{*tl}, size_t @code{numclumps}, size_t 
@code{minarea}, int @code{variance}, int @code{keepsmall}, gal_data_t 
@code{*sig}, gal_data_t @code{*sigind})
 @cindex Clump
-This function is usually called after @code{gal_label_oversegment}, and is
+This function is usually called after @code{gal_label_watershed}, and is
 used as a measure to idenfity which over-segmented ``clumps'' are real and
 which are noise.
 
@@ -26816,7 +26819,7 @@ is only done on pixels which are indexed in 
@code{indexs}. It is expected
 for @code{indexs} to be sorted by their values in @code{values}. If not
 sorted, the measurement may not be reliable. If sorted in a decreasing
 order, then clump building will start from their highest value and
-vice-versa. See the description of @code{gal_label_oversegment} for more on
+vice-versa. See the description of @code{gal_label_watershed} for more on
 @code{indexs}.
 
 Each ``clump'' (identified by a positive integer) is assumed to be
@@ -26878,7 +26881,7 @@ this function. For a demonstration see Columns 2 and 3 
of Figure 10 in
 @url{http://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa [2015]}.
 
 In many aspects, this function is very similar to over-segmentation
-(watershed algorithm, @code{gal_label_oversegment}). The big difference is
+(watershed algorithm, @code{gal_label_watershed}). The big difference is
 that in over-segmentation local maximums (that aren't touching any already
 labeled pixel) get a separate label. However, here the final number of
 labels will not change. All pixels that aren't directly touching a labeled
diff --git a/lib/gnuastro/label.h b/lib/gnuastro/label.h
index 1ee91de..fb21aa5 100644
--- a/lib/gnuastro/label.h
+++ b/lib/gnuastro/label.h
@@ -63,8 +63,8 @@ gal_data_t *
 gal_label_indexs(gal_data_t *labels, size_t numlabs, size_t minmapsize);
 
 size_t
-gal_label_oversegment(gal_data_t *values, gal_data_t *indexs,
-                      gal_data_t *label, size_t *topinds, int min0_max1);
+gal_label_watershed(gal_data_t *values, gal_data_t *indexs,
+                    gal_data_t *label, size_t *topinds, int min0_max1);
 
 void
 gal_label_clump_significance(gal_data_t *values, gal_data_t *std,
diff --git a/lib/label.c b/lib/label.c
index 876f1ec..fb33506 100644
--- a/lib/label.c
+++ b/lib/label.c
@@ -178,8 +178,8 @@ gal_label_indexs(gal_data_t *labels, size_t numlabs, size_t 
minmapsize)
 
 */
 size_t
-gal_label_oversegment(gal_data_t *values, gal_data_t *indexs,
-                      gal_data_t *labels, size_t *topinds, int min0_max1)
+gal_label_watershed(gal_data_t *values, gal_data_t *indexs,
+                    gal_data_t *labels, size_t *topinds, int min0_max1)
 {
   size_t ndim=values->ndim;
 
@@ -507,133 +507,6 @@ gal_label_oversegment(gal_data_t *values, gal_data_t 
*indexs,
 
 
 
-/* Grow the given labels without creating new ones. */
-void
-gal_label_grow_indexs(gal_data_t *labels, gal_data_t *indexs, int withrivers,
-                      int connectivity)
-{
-  int searchngb;
-  size_t *iarray=indexs->array;
-  int32_t n1, nlab, *olabel=labels->array;
-  size_t *s, *sf, thisround, ninds=indexs->size;
-  size_t *dinc=gal_dimension_increment(labels->ndim, labels->dsize);
-
-  /* Some basic sanity checks: */
-  label_check_type(indexs, GAL_TYPE_SIZE_T, "indexs", __func__);
-  label_check_type(labels, GAL_TYPE_INT32,  "labels", __func__);
-  if(indexs->ndim!=1)
-    error(EXIT_FAILURE, 0, "%s: `indexs' has to be a 1D array, but it is "
-          "%zuD", __func__, indexs->ndim);
-
-  /* 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
-     indexed pixels. Note that it will be corrected immediately after the
-     loop has started, it is just important to pass the `while'. */
-  thisround=ninds+1;
-  while( thisround > ninds )
-    {
-      /* `thisround' will keep the number of pixels to be inspected in this
-         round. `ninds' will count the number of pixels left without a
-         label by the end of this round. Since `ninds' comes from the
-         previous loop (or outside, for the first round) it has to be saved
-         in `thisround' to begin counting a fresh. */
-      thisround=ninds;
-      ninds=0;
-
-      /* Go over all the available indexs. NOTE: while the `indexs->array'
-         pointer remains unchanged, `indexs->size' can/will change (get
-         smaller) in every loop. */
-      sf = (s=indexs->array) + indexs->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 neighbors of this pixel. Note that since this
-             macro has multiple loops within it, we can't use
-             break. We'll use the `searchngb' variable instead. */
-          searchngb=1;
-          GAL_DIMENSION_NEIGHBOR_OP(*s, labels->ndim, labels->dsize,
-            connectivity, dinc,
-            {
-              if(searchngb)
-                {
-                  /* For easy reading. */
-                  nlab = olabel[nind];
-
-                  /* This neighbor's label is meaningful. */
-                  if(nlab>0)                   /* This is a real label. */
-                    {
-                      if(n1)       /* A prev. ngb label has been found. */
-                        {
-                          if( n1 != nlab )    /* Different label from   */
-                            {    /* prevously found ngb for this pixel. */
-                              n1=GAL_LABEL_RIVER;
-                              searchngb=0;
-                            }
-                        }
-                      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) searchngb=0;
-                        }
-                    }
-                }
-            } );
-
-          /* The loop over neighbors (above) finishes with three
-             possibilities:
-
-             n1==0                    --> No labeled neighbor was found.
-             n1==GAL_LABEL_RIVER      --> Connecting two labeled regions.
-             n1>0                     --> Only has one neighbouring label.
-
-             The first one means that no neighbors were 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)
-            {
-              /* Set the label. */
-              olabel[*s]=n1;
-
-              /* If this pixel is a river (can only happen when
-                 `withrivers' is zero), keep it in the loop, because we
-                 want the `indexs' dataset to contain all non-positive
-                 (non-labeled) pixels, including rivers. */
-              if(n1==GAL_LABEL_RIVER)
-                iarray[ ninds++ ] = *s;
-            }
-          else
-            iarray[ ninds++ ] = *s;
-
-          /* Correct the size of the `indexs' dataset. */
-          indexs->size = indexs->dsize[0] = ninds;
-        }
-      while(++s<sf);
-    }
-
-  /* Clean up. */
-  free(dinc);
-}
-
-
-
-
-
 
 
 
@@ -650,7 +523,7 @@ gal_label_grow_indexs(gal_data_t *labels, gal_data_t 
*indexs, int withrivers,
 
 
 /**********************************************************************/
-/*************            Clump peak intensity            *************/
+/*************             Clump significance             *************/
 /**********************************************************************/
 static int
 label_clump_significance_sanity(gal_data_t *values, gal_data_t *std,
@@ -719,7 +592,7 @@ label_clump_significance_sanity(gal_data_t *values, 
gal_data_t *std,
           "significance values must have NULL pointers for its `array' "
           "and `dsize' pointers (they will be allocated here)", func);
 
-  /* See if the clumps are to be build starting from local maxima or local
+  /* See if the clumps are to be built starting from local maxima or local
      minima. */
   af=(a=indexs->array)+indexs->size;
   do
@@ -730,11 +603,14 @@ label_clump_significance_sanity(gal_data_t *values, 
gal_data_t *std,
           first=f[*a];
         else
           {
-            /* Note that the elements may have equal values, so for
-               `second', we want the first non-blank AND different
-               value. */
-            if( isnan(second) && f[*a]!=first )
-              second=f[*a];
+            if( isnan(second) )
+              {
+                /* Note that the elements may have equal values, so for
+                   `second', we want the first non-blank AND different
+                   value. */
+                if( f[*a]!=first )
+                  second=f[*a];
+              }
             else
               break;
           }
@@ -763,11 +639,7 @@ label_clump_significance_sanity(gal_data_t *values, 
gal_data_t *std,
    below.*/
 enum infocols
   {
-    INFO_X,              /* Flux weighted X center col, 0 by C std.       */
-    INFO_Y,              /* Flux weighted Y center col.                   */
-    INFO_Z,              /* Flux weighted Z center col.                   */
-    INFO_SFF,            /* Sum of non-negative pixels (for X,Y).         */
-    INFO_INSTD,          /* Standard deviation at clump center.           */
+    INFO_STD,            /* Standard deviation.                           */
     INFO_INAREA,         /* Tatal area within clump.                      */
     INFO_RIVAREA,        /* Tatal area within rivers around clump.        */
     INFO_PEAK_RIVER,     /* Peak (min or max) river value around a clump. */
@@ -779,17 +651,16 @@ static void
 label_clump_significance_raw(gal_data_t *values_d, gal_data_t *std_d,
                              gal_data_t *label_d, gal_data_t *indexs,
                              struct gal_tile_two_layer_params *tl,
-                             double *info, size_t numclumps, size_t minarea,
-                             int variance)
+                             double *info)
 {
   size_t ndim=values_d->ndim, *dsize=values_d->dsize;
 
   double *row;
   size_t i, *a, *af, ii, coord[3];
   size_t nngb=gal_dimension_num_neighbors(ndim);
+  int32_t nlab, *ngblabs, *label=label_d->array;
   float *values=values_d->array, *std=std_d->array;
   size_t *dinc=gal_dimension_increment(ndim, dsize);
-  int32_t lab, nlab, *ngblabs, *label=label_d->array;
 
   /* Allocate the array to keep the neighbor labels of river pixels. */
   ngblabs=gal_pointer_allocate(GAL_TYPE_INT32, nngb, 0, __func__, "ngblabs");
@@ -807,15 +678,6 @@ label_clump_significance_raw(gal_data_t *values_d, 
gal_data_t *std_d,
 
             /* Get the area and flux. */
             ++row[ INFO_INAREA ];
-            if( values[*a]>0.0f )
-              {
-                gal_dimension_index_to_coord(*a, ndim, dsize, coord);
-                row[   INFO_SFF ] += values[*a];
-                row[   INFO_X   ] += values[*a] * coord[0];
-                row[   INFO_Y   ] += values[*a] * coord[1];
-                if(ndim==3)
-                  row[ INFO_Z   ] += values[*a] * coord[2];
-              }
 
             /* In the loop `INFO_INAREA' is just the pixel counter of this
                clump. The pixels are sorted by flux (decreasing for
@@ -862,7 +724,21 @@ label_clump_significance_raw(gal_data_t *values_d, 
gal_data_t *std_d,
 
                         ++row[INFO_RIVAREA];
                         if( row[INFO_RIVAREA]==1.0f )
-                          row[INFO_PEAK_RIVER] = values[*a];
+                          {
+                            /* Get the maximum river value. */
+                            row[INFO_PEAK_RIVER] = values[*a];
+
+                            /* Get the standard deviation. */
+                            if(std_d->size==1 || std_d->size==values_d->size)
+                              row[INFO_STD]=std_d->size==1?std[0]:std[*a];
+                            else
+                              {
+                                gal_dimension_index_to_coord(*a, ndim, dsize,
+                                                             coord);
+                                row[INFO_STD]=
+                                  std[gal_tile_full_id_from_coord(tl,coord)];
+                              }
+                          }
                       }
                   }
               } );
@@ -870,53 +746,6 @@ label_clump_significance_raw(gal_data_t *values_d, 
gal_data_t *std_d,
       }
   while(++a<af);
 
-  /* Based on the position of each clump, find a representative standard
-     deviation. */
-  for(lab=1; lab<=numclumps; ++lab)
-    {
-      /* To help in reading. */
-      row = &info [ lab * INFO_NCOLS ];
-
-      /* The calculations are only necessary for the clumps that satisfy
-         the minimum area. There is no need to waste time on the smaller
-         ones. */
-      if ( row[INFO_INAREA] > minarea )
-        {
-          /* It might happen that none of the pixels were positive
-             (especially over the undetected regions). In that case, set
-             the total area of the clump to zero so it is no longer
-             considered.*/
-          if( row[INFO_SFF]==0.0f )
-            row[INFO_INAREA]=0;
-          else
-            {
-              /* Find the coordinates of the clump's weighted center. */
-              coord[0]=GAL_DIMENSION_FLT_TO_INT(row[INFO_X]/row[INFO_SFF]);
-              coord[1]=GAL_DIMENSION_FLT_TO_INT(row[INFO_Y]/row[INFO_SFF]);
-              if(ndim==3)
-                coord[2]=GAL_DIMENSION_FLT_TO_INT(row[INFO_Z]/row[INFO_SFF]);
-
-              /* Find the corresponding standard deviation. */
-              row[INFO_INSTD]=( std_d->size>1
-                                ? ( std_d->size==values_d->size
-                                    ? std[gal_dimension_coord_to_index(ndim,
-                                                             dsize, coord)]
-                                    : std[gal_tile_full_id_from_coord(tl,
-                                                                    coord)] )
-                                : std[0] );
-              if(variance) row[INFO_INSTD] = sqrt(row[INFO_INSTD]);
-
-              /* For a check
-              printf("---------\n");
-              printf("\t%f --> %zu\n", row[INFO_Y]/row[INFO_SFF], coord[1]);
-              printf("\t%f --> %zu\n", row[INFO_X]/row[INFO_SFF], coord[0]);
-              printf("%u: (%zu, %zu): %.3f\n", lab, coord[1]+1,
-                     coord[0]+1, row[INFO_INSTD]);
-              */
-            }
-        }
-    }
-
   /* Clean up. */
   free(dinc);
   free(ngblabs);
@@ -937,9 +766,9 @@ gal_label_clump_significance(gal_data_t *values, gal_data_t 
*std,
   double *info;
   int max1_min0;
   float *sigarr;
+  double C, R, S, *row;
   int32_t *indarr=NULL;
   size_t i, ind, counter=0;
-  double C, R, S, Ni, *row;
   size_t tablen=numclumps+1;
 
   /* If there were no initial clumps, then ignore this function. */
@@ -974,8 +803,7 @@ gal_label_clump_significance(gal_data_t *values, gal_data_t 
*std,
 
 
   /* First, get the raw information necessary for making the S/N table. */
-  label_clump_significance_raw(values, std, label, indexs, tl, info,
-                               numclumps, minarea, variance);
+  label_clump_significance_raw(values, std, label, indexs, tl, info);
 
 
   /* Calculate the signal to noise ratio for successful clumps */
@@ -988,19 +816,27 @@ gal_label_clump_significance(gal_data_t *values, 
gal_data_t *std,
 
       /* If we have a sufficient area and any rivers were actually found
          for this clump, then do the measurement. */
-      Ni=row[ INFO_INAREA ];
-      if( Ni>minarea && row[ INFO_RIVAREA ])
+      if( row[ INFO_INAREA ]>minarea && row[ INFO_RIVAREA ])
         {
-          /* Calculate the significance ratio, if `keepsmall' is not
+          /* Set the index to write the values. If `keepsmall' is not
              called, we don't care about the IDs of the clumps anymore, so
-             store the Signal to noise ratios contiguously (for easy
-             sorting and etc). Note that counter will always be smaller and
-             equal to i. */
-          S   = row[ INFO_INSTD       ];
-          R   = row[ INFO_PEAK_RIVER  ];
-          C   = row[ INFO_PEAK_CENTER ];
+             store the signal-to-noise ratios contiguously. Note that
+             counter will always be smaller and equal to i. */
           ind = keepsmall ? i : counter++;
 
+          /* For easy reading. */
+          R   = row[ INFO_PEAK_RIVER  ];
+          C   = row[ INFO_PEAK_CENTER ];
+          S   = variance ? sqrt(row[ INFO_STD ]) : row[ INFO_STD ];
+
+          /* NEGATIVE VALUES: Rivers are also defined on the edges of the
+             image and on pixels touching blank pixels. In a strong
+             gradient, such sitations can cause the river to be
+             larger/smaller than the minimum/maximum within the clump. This
+             can only happen in very strong gradients, so for now, I think
+             it is safe to ignore that clump (its negative value will
+             automatically discard it). Later, if we find a problem with
+             this, we'll have to figure out a better solution. */
           if(sigind) indarr[ind]=i;
           sigarr[ind] = ( max1_min0 ? (C-R) : (R-C) ) / S;
         }
@@ -1029,3 +865,148 @@ gal_label_clump_significance(gal_data_t *values, 
gal_data_t *std,
   /* Clean up. */
   free(info);
 }
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/**********************************************************************/
+/*************               Growing labels               *************/
+/**********************************************************************/
+/* Grow the given labels without creating new ones. */
+void
+gal_label_grow_indexs(gal_data_t *labels, gal_data_t *indexs, int withrivers,
+                      int connectivity)
+{
+  int searchngb;
+  size_t *iarray=indexs->array;
+  int32_t n1, nlab, *olabel=labels->array;
+  size_t *s, *sf, thisround, ninds=indexs->size;
+  size_t *dinc=gal_dimension_increment(labels->ndim, labels->dsize);
+
+  /* Some basic sanity checks: */
+  label_check_type(indexs, GAL_TYPE_SIZE_T, "indexs", __func__);
+  label_check_type(labels, GAL_TYPE_INT32,  "labels", __func__);
+  if(indexs->ndim!=1)
+    error(EXIT_FAILURE, 0, "%s: `indexs' has to be a 1D array, but it is "
+          "%zuD", __func__, indexs->ndim);
+
+  /* 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
+     indexed pixels. Note that it will be corrected immediately after the
+     loop has started, it is just important to pass the `while'. */
+  thisround=ninds+1;
+  while( thisround > ninds )
+    {
+      /* `thisround' will keep the number of pixels to be inspected in this
+         round. `ninds' will count the number of pixels left without a
+         label by the end of this round. Since `ninds' comes from the
+         previous loop (or outside, for the first round) it has to be saved
+         in `thisround' to begin counting a fresh. */
+      thisround=ninds;
+      ninds=0;
+
+      /* Go over all the available indexs. NOTE: while the `indexs->array'
+         pointer remains unchanged, `indexs->size' can/will change (get
+         smaller) in every loop. */
+      sf = (s=indexs->array) + indexs->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 neighbors of this pixel. Note that since this
+             macro has multiple loops within it, we can't use
+             break. We'll use the `searchngb' variable instead. */
+          searchngb=1;
+          GAL_DIMENSION_NEIGHBOR_OP(*s, labels->ndim, labels->dsize,
+            connectivity, dinc,
+            {
+              if(searchngb)
+                {
+                  /* For easy reading. */
+                  nlab = olabel[nind];
+
+                  /* This neighbor's label is meaningful. */
+                  if(nlab>0)                   /* This is a real label. */
+                    {
+                      if(n1)       /* A prev. ngb label has been found. */
+                        {
+                          if( n1 != nlab )    /* Different label from   */
+                            {    /* prevously found ngb for this pixel. */
+                              n1=GAL_LABEL_RIVER;
+                              searchngb=0;
+                            }
+                        }
+                      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) searchngb=0;
+                        }
+                    }
+                }
+            } );
+
+          /* The loop over neighbors (above) finishes with three
+             possibilities:
+
+             n1==0                    --> No labeled neighbor was found.
+             n1==GAL_LABEL_RIVER      --> Connecting two labeled regions.
+             n1>0                     --> Only has one neighbouring label.
+
+             The first one means that no neighbors were 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)
+            {
+              /* Set the label. */
+              olabel[*s]=n1;
+
+              /* If this pixel is a river (can only happen when
+                 `withrivers' is zero), keep it in the loop, because we
+                 want the `indexs' dataset to contain all non-positive
+                 (non-labeled) pixels, including rivers. */
+              if(n1==GAL_LABEL_RIVER)
+                iarray[ ninds++ ] = *s;
+            }
+          else
+            iarray[ ninds++ ] = *s;
+
+          /* Correct the size of the `indexs' dataset. */
+          indexs->size = indexs->dsize[0] = ninds;
+        }
+      while(++s<sf);
+    }
+
+  /* Clean up. */
+  free(dinc);
+}
diff --git a/tests/during-dev.sh b/tests/during-dev.sh
index 81cfaf7..e7b9a8f 100755
--- a/tests/during-dev.sh
+++ b/tests/during-dev.sh
@@ -140,13 +140,22 @@ if make -j$numjobs -C "$builddir"; then
     # the last line in the configuration file doesn't actualy end with a
     # new line (in which case the appended string will be added to the end
     # of the last line).
-    cp "$srcdir/bin/gnuastro.conf" "$srcdir/bin/$utilname/"*.conf     \
-       .gnuastro/
+    if [ $utilname = buildprog ]; then
+        extraopts="--la=$builddir/lib/libgnuastro.la"
+        topconfdir="$builddir"
+    else
+        topconfdir="$srcdir"
+    fi
+    cp "$srcdir/bin/gnuastro.conf"                                    \
+       "$topconfdir/bin/$utilname/ast$utilname.conf" .gnuastro/
+
+    # Append `lastconfig' option to `gnuastro.conf', so the program doesn't
+    # go into the system headers.
     echo ""               >> .gnuastro/gnuastro.conf
     echo " lastconfig 1"  >> .gnuastro/gnuastro.conf
 
     # Run the built utility with the given arguments and options.
-    "$utility" $arguments $options
+    "$utility" $arguments $options $extraopts
 
     # Clean up.
     rm -rf .gnuastro