[gnuastro-commits] master 442d052 061/113: Imported work in master, no conflicts

[Mohammad Akhlaghi]

branch: master
commit 442d0525b26155f78fd213b87a9e85e06c289ee1
Merge: 6f85919 f98b018
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    Imported work in master, no conflicts
    
    There were no conflicts in this merge.
---
 NEWS                        |  12 +-
 bin/noisechisel/args.h      |   2 +-
 bin/noisechisel/detection.c |  14 ++-
 bin/segment/args.h          |   2 +-
 bin/segment/clumps.c        |  59 +--------
 bin/segment/segment.c       |  80 ++++++-------
 bin/statistics/ui.c         |   2 +-
 doc/gnuastro.texi           | 282 ++++++++++++++++++++++++++++++++++++++------
 lib/gnuastro/label.h        |   8 +-
 lib/gnuastro/statistics.h   |  13 +-
 lib/label.c                 | 239 +++++++++++++++++++++++++------------
 lib/statistics.c            | 138 ++++++++++++++--------
 12 files changed, 565 insertions(+), 286 deletions(-)

diff --git a/NEWS b/NEWS
index 3ab0703..ae73509 100644
--- a/NEWS
+++ b/NEWS
@@ -95,7 +95,6 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
 ** Removed features
 
   NoiseChisel:
-
     - Segmentation (and thus the options below) moved to the new Segment
       program: --onlydetection, --segsnminarea, --checkclumpsn, --segquant,
       --keepmaxnearriver, --gthresh, --minriverlength, --objbordersn,
@@ -106,6 +105,13 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
   MakeCatalog:
     --skysubtracted: no longer necessary, included in noise measuremnts.
 
+  Library:
+    - The macros `GAL_STATISTICS_SORTED_NOT',
+      `GAL_STATISTICS_SORTED_INVALID', `GAL_STATISTICS_SORTED_INCREASING',
+      `GAL_STATISTICS_SORTED_DECREASING': these macros are removed because
+      we already have the `GAL_DATA_FLAG_SORT*'' bit-flags in `gal_data_t'.
+
+
 ** Changed features
 
   Arithmetic:
@@ -152,6 +158,10 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
         longer needs the last two arguments. A subsequent call to
         `gal_wcs_read' can be used to read the WCS information in the file.
 
+    gal_statistics_is_sorted: can now also update the bit-flags regarding
+        the sorted nature of the input (to optimize future calls to the
+        function).
+
     gal_statistics_quantile_function: returns `inf' or `-inf' if the given
         value is smaller than the minimum or larger than the maximum of the
         input dataset's range. Until now, it would return blank in such
diff --git a/bin/noisechisel/args.h b/bin/noisechisel/args.h
index 76b6cda..3e6e159 100644
--- a/bin/noisechisel/args.h
+++ b/bin/noisechisel/args.h
@@ -416,7 +416,7 @@ struct argp_option program_options[] =
       UI_GROUP_DETECTION,
       &p->snquant,
       GAL_TYPE_FLOAT32,
-      GAL_OPTIONS_RANGE_GT_0_LT_1,
+      GAL_OPTIONS_RANGE_GE_0_LE_1,
       GAL_OPTIONS_MANDATORY,
       GAL_OPTIONS_NOT_SET
     },
diff --git a/bin/noisechisel/detection.c b/bin/noisechisel/detection.c
index 5bb62cc..e0ebf9e 100644
--- a/bin/noisechisel/detection.c
+++ b/bin/noisechisel/detection.c
@@ -1027,11 +1027,11 @@ detection_quantile_expand(struct noisechiselparams *p, 
gal_data_t *workbin)
                                    NULL);
 
       /* Fill in the diffuse indexs and initialize the objects dataset. */
-      b=workbin->array;
-      arr=p->conv->array;
-      d=diffuseindexs->array;
-      e_th=p->exp_thresh_full->array;
-      of=(o=p->olabel->array)+p->olabel->size;
+      b    = workbin->array;
+      arr  = p->conv->array;
+      d    = diffuseindexs->array;
+      e_th = p->exp_thresh_full->array;
+      of   = (o=p->olabel->array) + p->olabel->size;
       do
         {
           /* If the binary value is 1, then we want an initial label of 1
@@ -1048,7 +1048,9 @@ detection_quantile_expand(struct noisechiselparams *p, 
gal_data_t *workbin)
         }
       while(++o<of);
 
-      /* Expand the detections. */
+      /* Expand the detections. Note that because we are only concerned
+         with those regions that are touching a detected region, it is
+         irrelevant to sort the dataset. */
       gal_label_grow_indexs(p->olabel, diffuseindexs, 0, p->olabel->ndim);
 
       /* Only keep the 1 valued pixels in the binary array and fill its
diff --git a/bin/segment/args.h b/bin/segment/args.h
index 78a56c8..3ed23ad 100644
--- a/bin/segment/args.h
+++ b/bin/segment/args.h
@@ -337,7 +337,7 @@ struct argp_option program_options[] =
       UI_GROUP_SEGMENTATION,
       &p->snquant,
       GAL_TYPE_FLOAT32,
-      GAL_OPTIONS_RANGE_GT_0,
+      GAL_OPTIONS_RANGE_GE_0_LE_1,
       GAL_OPTIONS_MANDATORY,
       GAL_OPTIONS_NOT_SET
     },
diff --git a/bin/segment/clumps.c b/bin/segment/clumps.c
index b27fce1..05628d8 100644
--- a/bin/segment/clumps.c
+++ b/bin/segment/clumps.c
@@ -136,7 +136,10 @@ clumps_grow_prepare_initial(struct clumps_thread_params 
*cltprm)
      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. */
+     the exact number.
+
+     Also note that since `indexs' is already sorted, therefore
+     `diffuseindexs' will also be already sorted. */
   cltprm->diffuseindexs=gal_data_alloc(NULL, GAL_TYPE_SIZE_T, 1,
                                        cltprm->indexs->dsize, NULL, 0,
                                        p->cp.minmapsize, NULL, NULL, NULL);
@@ -1050,60 +1053,6 @@ clumps_true_find_sn_thresh(struct segmentparams *p)
 /***********************************************************************/
 /*****************           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 segmentparams *p)
-{
-  size_t i, *areas;
-  int32_t *a, *l, *lf;
-  gal_data_t *labindexs=gal_data_array_calloc(p->numdetections+1);
-
-  /* Find the area in each detected object (to see how much space we need
-     to allocate). If blank values are present, an extra check is
-     necessary, so to get faster results when there aren't any blank
-     values, we'll also do a check. */
-  areas=gal_data_calloc_array(GAL_TYPE_SIZE_T, p->numdetections+1, __func__,
-                              "areas");
-  lf=(l=p->olabel->array)+p->olabel->size;
-  do
-    if(*l>0)  /* Only labeled regions: *l==0 (undetected), *l<0 (blank). */
-      ++areas[*l];
-  while(++l<lf);
-
-  /* For a check.
-  for(i=0;i<p->numdetections+1;++i)
-    printf("detection %zu: %zu\n", i, areas[i]);
-  exit(0);
-  */
-
-  /* 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.*/
-  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>0)  /* No undetected regions (*l==0), or blank (<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)
diff --git a/bin/segment/segment.c b/bin/segment/segment.c
index 312f325..1cf0154 100644
--- a/bin/segment/segment.c
+++ b/bin/segment/segment.c
@@ -588,19 +588,21 @@ segment_on_threads(void *in_prm)
       /* If the user wanted to check the segmentation steps or the clump
          S/N values in a table, then we have to stop the process at this
          point. */
-      if(clprm->step==1 || p->checksn)
+      if( clprm->step==1 || (p->checksn && !p->continueaftercheck ) )
         { gal_data_free(topinds); continue; }
 
 
-      /* Only keep true clumps and abort if the user only wants clumps. */
+      /* Only keep true clumps. */
       clumps_det_keep_true_relabel(&cltprm);
       gal_data_free(topinds);
-      if(clprm->step==2) continue;
 
 
       /* When only clumps are desired ignore the rest of the process. */
       if(!p->onlyclumps)
         {
+          /* Abort the looping here if we don't only want clumps. */
+          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
@@ -799,13 +801,12 @@ static void
 segment_detections(struct segmentparams *p)
 {
   char *msg;
-  int continuecheck=1;
   struct clumps_params clprm;
   gal_data_t *labindexs, *claborig, *demo=NULL;
 
 
   /* Get the indexs of all the pixels in each label. */
-  labindexs=clumps_det_label_indexs(p);
+  labindexs=gal_label_indexs(p->olabel, p->numdetections, p->cp.minmapsize);
 
 
   /* Initialize the necessary thread parameters. Note that since the object
@@ -838,7 +839,16 @@ segment_detections(struct segmentparams *p)
 
 
       /* Do each step. */
-      while(clprm.step<8 && continuecheck)
+      while( clprm.step<8
+
+             /* When the user only wanted clumps, there is no point in
+                continuing beyond step 2. */
+             && !(p->onlyclumps && clprm.step>2)
+
+             /* When the user just wants to check the clump S/N values,
+                then break out of the loop, we don't need the rest of the
+                process any more. */
+             && !( (p->checksn && !p->continueaftercheck) && clprm.step>1 ) )
         {
           /* Reset the temporary copy of clabel back to its original. */
           if(clprm.step>1)
@@ -881,28 +891,18 @@ segment_detections(struct segmentparams *p)
               break;
 
             case 3:
-              /* If the user only wanted clumps, there is no point in
-                 continuing after this point. We DON'T WANT this at the end
-                 of `2', otherwise, the `DET_CLUMPS_TRUE' extension will be
-                 different when `--onlyclumps' is called and when it
-                 isn't.*/
-              if(p->onlyclumps)
-                continuecheck=0;
-              else
+              demo=p->olabel;
+              demo->name = "DET_CLUMPS_GROWN";
+              if(!p->cp.quiet)
                 {
-                  demo=p->olabel;
-                  demo->name = "DET_CLUMPS_GROWN";
-                  if(!p->cp.quiet)
-                    {
-                      gal_timing_report(NULL, "Identify objects...",
-                                        1);
-                      if( asprintf(&msg, "True clumps grown                  "
-                                   "(HDU: `%s').", demo->name)<0 )
-                        error(EXIT_FAILURE, 0, "%s: asprintf allocation",
-                              __func__);
-                      gal_timing_report(NULL, msg, 2);
-                      free(msg);
-                    }
+                  gal_timing_report(NULL, "Identify objects...",
+                                    1);
+                  if( asprintf(&msg, "True clumps grown                  "
+                               "(HDU: `%s').", demo->name)<0 )
+                    error(EXIT_FAILURE, 0, "%s: asprintf allocation",
+                          __func__);
+                  gal_timing_report(NULL, msg, 2);
+                  free(msg);
                 }
               break;
 
@@ -969,18 +969,8 @@ segment_detections(struct segmentparams *p)
                     clprm.step);
             }
 
-          /* 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. */
-          if(continuecheck)
-            gal_fits_img_write(demo, p->segmentationname, NULL, PROGRAM_NAME);
-
-          /* If the user wanted to check the clump S/N values, then break
-             out of the loop, we don't need the rest of the process any
-             more. */
-          if( clprm.step==1
-              && ( p->checksn && !p->continueaftercheck ) ) break;
+          /* Write the demonstration array into the check image.  */
+          gal_fits_img_write(demo, p->segmentationname, NULL, PROGRAM_NAME);
 
           /* Increment the step counter. */
           ++clprm.step;
@@ -998,13 +988,15 @@ segment_detections(struct segmentparams *p)
     }
 
 
-  /* 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 and
+     abort Segment if necessary. */
+  if(p->checksn) segment_save_sn_table(&clprm);
 
 
-  /* If the user wanted to see the S/N table, then make the S/N table. */
-  if(p->checksn) segment_save_sn_table(&clprm);
+  /* Write the final number of objects and clumps to be used beyond this
+     function. */
+  p->numclumps=clprm.totclumps;
+  p->numobjects=clprm.totobjects;
 
 
   /* Clean up allocated structures and destroy the mutex. */
diff --git a/bin/statistics/ui.c b/bin/statistics/ui.c
index 0784419..3819c11 100644
--- a/bin/statistics/ui.c
+++ b/bin/statistics/ui.c
@@ -704,7 +704,7 @@ ui_make_sorted_if_necessary(struct statisticsparams *p)
      point it to the input.*/
   if(is_necessary)
     {
-      if( gal_statistics_is_sorted(p->input) )
+      if( gal_statistics_is_sorted(p->input, 1) )
         p->sorted=p->input;
       else
         {
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 9a3684f..3338cd5 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -581,6 +581,7 @@ Gnuastro library
 * Matching::                    Matching catalogs based on position.
 * Statistical operations::      Functions for basic statistics.
 * Binary datasets::             Datasets that can only have values of 0 or 1.
+* Labeled datasets::            Working with Segmented/labeled datasets.
 * Convolution functions::       Library functions to do convolution.
 * Interpolation::               Interpolate (over blank values possibly).
 * Git wrappers::                Wrappers for functions in libgit2.
@@ -20418,6 +20419,7 @@ documentation will correspond to your installed version.
 * Matching::                    Matching catalogs based on position.
 * Statistical operations::      Functions for basic statistics.
 * Binary datasets::             Datasets that can only have values of 0 or 1.
+* Labeled datasets::            Working with Segmented/labeled datasets.
 * Convolution functions::       Library functions to do convolution.
 * Interpolation::               Interpolate (over blank values possibly).
 * Git wrappers::                Wrappers for functions in libgit2.
@@ -21117,24 +21119,25 @@ below).
 @deftypefun int gal_blank_present (gal_data_t @code{*input}, int 
@code{updateflag})
 Return 1 if the dataset has a blank value and zero if it doesn't. Before
 checking the dataset, this function will look at @code{input}'s flags. If
-the @code{GAL_DATA_FLAG_HASBLANK} or @code{GAL_DATA_FLAG_DONT_CHECK_ZERO}
-bits of @code{input->flag} are set to 1, this function will not do any
-check and will just use the information in the flags. This can greatly
-speed up processing when a dataset needs to be checked multiple times.
+the @code{GAL_DATA_FLAG_BLANK_CH} bit of @code{input->flag} is on, this
+function will not do any check and will just use the information in the
+flags. This can greatly speed up processing when a dataset needs to be
+checked multiple times.
 
-If you want to re-check a dataset which has non-zero flags, then explicitly
-set the appropriate flag to zero before calling this function. When there
-are no other flags, you can just set @code{input->flags} to zero, otherwise
-you can use this expression:
+When the dataset's flags were not used and @code{updateflags} is non-zero,
+this function will set the flags appropriately to avoid having to re-check
+the dataset in future calls. When @code{updateflags==0}, this function has
+no side-effects on the dataset: it will not toggle the flags.
+
+If you want to re-check a dataset with the blank-value-check flag already
+set (for example if you have made changes to it), then explicitly set the
+@code{GAL_DATA_FLAG_BLANK_CH} bit to zero before calling this
+function. When there are no other flags, you can just set the flags to zero
+(@code{input->flags=0}), otherwise you can use this expression:
 
 @example
-input->flags &= ~(GAL_DATA_FLAG_HASBLANK | GAL_DATA_FLAG_USE_ZERO);
+input->flags &= ~GAL_DATA_FLAG_BLANK_CH;
 @end example
-
-When @code{updateflags} is zero, this function has no side-effects on the
-dataset: it will not toggle the flags. When the dataset's flags were not
-used and @code{updateflags} is non-zero, this function will set the flags
-appropriately to avoid having to re-check the dataset in future calls.
 @end deftypefun
 
 
@@ -21146,16 +21149,18 @@ those that aren't.
 
 
 @deftypefun void gal_blank_remove (gal_data_t @code{*input})
-Remove blank elements from a dataset, convert it to a 1D dataset, and
-adjust the size properly (the number of non-blank elements). In practice
-this function doesn't @code{realloc} the input array, it just shifts the
-blank elements to the end and adjusts the size elements of the
-@code{gal_data_t}, see @ref{Generic data container}.
+Remove blank elements from a dataset, convert it to a 1D dataset, adjust
+the size properly (the number of non-blank elements), and toggle the
+blank-value-related bit-flags. In practice this function doesn't
+@code{realloc} the input array, it just shifts the blank elements to the
+end and adjusts the size elements of the @code{gal_data_t}, see
+@ref{Generic data container}.
 
 If all the elements were blank, then @code{input->size} will be zero. This
 is thus a good parameter to check after calling this function to see if
 there actually were any non-blank elements in the input or not and take the
-appropriate measure. This can help avoid strange bugs in later steps.
+appropriate measure. This check is highly recommended because it will avoid
+strange bugs in later steps.
 @end deftypefun
 
 @deftypefun {char *} gal_blank_as_string (uint8_t @code{type}, int 
@code{width})
@@ -25188,14 +25193,6 @@ symmetricity of the derived mode is less than this 
value, all the returned
 values by @code{gal_statistics_mode} will have a value of NaN.
 @end deffn
 
-@deffn  Macro GAL_STATISTICS_SORTED_NOT
-@deffnx Macro GAL_STATISTICS_SORTED_INVALID
-@deffnx Macro GAL_STATISTICS_SORTED_INCREASING
-@deffnx Macro GAL_STATISTICS_SORTED_DECREASING
-Macros used to identify if the dataset is sorted and increasing, sorted and
-decreasing or not sorted.
-@end deffn
-
 @deffn  Macro GAL_STATISTICS_BINS_INVALID
 @deffnx Macro GAL_STATISTICS_BINS_REGULAR
 @deffnx Macro GAL_STATISTICS_BINS_IRREGULAR
@@ -25335,25 +25332,55 @@ plot (with a maximum of one).
 @end deftypefun
 
 
-@deftypefun int gal_statistics_is_sorted (gal_data_t @code{*input})
-Return the respective sort macro (see above) for the @code{input}
-dataset. If the dataset has zero elements, then
-@code{GAL_STATISTICS_SORTED_INVALID} is returned.
+@deftypefun int gal_statistics_is_sorted (gal_data_t @code{*input}, int 
@code{updateflags})
+Return @code{0} if the input is not sorted, if it is sorted, this function
+will return @code{1} and @code{2} if it is increasing or decreasing,
+respectively. This function will abort with an error if @code{input} has
+zero elements. This function will only look into the dataset if the
+@code{GAL_DATA_FLAG_SORT_CH} bit of @code{input->flag} is @code{0}, see
+@ref{Generic data container}.
+
+When the flags don't indicate a previous check @emph{and}
+@code{updateflags} is non-zero, this function will set the flags
+appropriately to avoid having to re-check the dataset in future calls (this
+can be very useful when repeated checks are necessary). When
+@code{updateflags==0}, this function has no side-effects on the dataset: it
+will not toggle the flags.
+
+If you want to re-check a dataset with the blank-value-check flag already
+set (for example if you have made changes to it), then explicitly set the
+@code{GAL_DATA_FLAG_SORT_CH} bit to zero before calling this function. When
+there are no other flags, you can simply set the flags to zero (with
+@code{input->flags=0}), otherwise you can use this expression:
+
+@example
+input->flags &= ~GAL_DATA_FLAG_SORT_CH;
+@end example
+
 @end deftypefun
 
 @deftypefun void gal_statistics_sort_increasing (gal_data_t @code{*input})
-Sort the input dataset (in place) in an increasing order.
+Sort the input dataset (in place) in an increasing order and toggle the
+sort-related bit flags accordingly.
 @end deftypefun
 
 @deftypefun void gal_statistics_sort_decreasing (gal_data_t @code{*input})
-Sort the input dataset (in place) in a decreasing order.
+Sort the input dataset (in place) in a decreasing order and toggle the
+sort-related bit flags accordingly.
 @end deftypefun
 
 @deftypefun {gal_data_t *} gal_statistics_no_blank_sorted (gal_data_t 
@code{*input}, int @code{inplace})
 Remove all the blanks and sort the input dataset. If @code{inplace} is
-non-zero this will happen on the input dataset (and the output will be the
-same as the input). However, if @code{inplace} is zero, this function will
-allocate a new copy of the dataset that is sorted and has no blank values.
+non-zero this will happen on the input dataset (and the output dataset will
+be the input dataset). However, if @code{inplace} is zero, this function
+will allocate a new copy of the dataset that is sorted and has no blank
+values.
+
+This function uses the bit flags of the input, so if you have modified the
+dataset, set @code{input->flags=0} before calling this function. Also note
+that @code{inplace} is only for the dataset elements. Therefore even when
+@code{inplace==0}, if the input is already sorted @emph{and} has no blank
+values, then the flags will be updated to show this.
 
 If all the elements were blank, then the returned dataset's @code{size}
 will be zero. This is thus a good parameter to check after calling this
@@ -25463,7 +25490,7 @@ above.
 
 
 
-@node Binary datasets, Convolution functions, Statistical operations, Gnuastro 
library
+@node Binary datasets, Labeled datasets, Statistical operations, Gnuastro 
library
 @subsection Binary datasets (@file{binary.h})
 
 @cindex Thresholding
@@ -25656,7 +25683,184 @@ can be set with @code{connectivity}. This function 
currently only works on
 a 2D dataset.
 @end deftypefun
 
-@node Convolution functions, Interpolation, Binary datasets, Gnuastro library
+@node Labeled datasets, Convolution functions, Binary datasets, Gnuastro 
library
+@subsection Labeled datasets (@file{label.h})
+
+A labeled dataset is one where each element/pixel has an integer label (or
+counter). The label identifies the group/class that the element belongs
+to. This form of labeling allows the higher-level study of all pixels
+within a certain class.
+
+For example, to detect objects/targets in an image/dataset, you can apply a
+threshold to separate the noise from the signal (to detect diffuse signal,
+a threshold is useless and more advanced methods are necessary, for example
+@ref{NoiseChisel}). But the output of detection is a binary dataset (which
+is just a very low-level labeling of @code{0} for noise and @code{1} for
+signal).
+
+The raw detection map is therefore hardly useful for any kind of analysis
+on objects/targets in the image. One solution is to use a
+connected-components algorithm (see @code{gal_binary_connected_components}
+in @ref{Binary datasets}). It is a simple and useful way to separate/label
+connected patches in the foreground. This higher-level (but still
+elementary) labeling therefore allows you to count how many connected
+patches of signal there are in the dataset and is a major improvement
+compared to the raw detection.
+
+However, when your objects/targets are touching, the simple connected
+components algorithm is not enough and a still higher-level labeling
+mechanism is necessary. This brings us to the necessity of the functions in
+this part of Gnuastro's library. The main inputs to the functions in this
+section are already labeled datasets (for example with the connected
+components algorithm above).
+
+Each of the labeled regions are independent of each other (the labels
+specify different classes of targets). Therefore, especially in large
+datasets, it is often useful to process each label on independent CPU
+threads in parallel rather than in series. Therefore the functions of this
+section actually use an array of pixel/element indexs (belonging to each
+label/class) as the main identifier of a region. Using indexs will also
+allow processing of overlapping labels (for example in deblending
+problems). Just note that overlapping labels are not yet implemented, but
+planned. You can use @code{gal_label_indexs} to generate lists of indexs
+belonging to separate classes from the labeled input.
+
+@deffn  Macro GAL_LABEL_INIT
+@deffnx Macro GAL_LABEL_RIVER
+@deffnx Macro GAL_LABEL_TMPCHECK
+Special negative integer values used internally by some of the functions in
+this section. Recall that meaningful labels are considered to be positive
+integers (@mymath{\geq1}). Zero is conventionally kept for regions with no
+labels, therefore negative integers can be used for any extra
+classification in the labeled datasets.
+@end deffn
+
+@deftypefun {gal_data_t *} gal_label_indexs (gal_data_t @code{*labels}, size_t 
@code{numlabs}, size_t @code{minmapsize})
+
+Return an array of @code{gal_data_t} containers, each containing the pixel
+indexs of the respective label (see @ref{Generic data
+container}). @code{labels} contains the label of each element and has to
+have an @code{GAL_TYPE_INT32} type (see @ref{Library data types}). Only
+positive (greater than zero) values in @code{labels} will be used/indexed,
+other elements will be ignored.
+
+Meaningful labels start from @code{1} and not @code{0}, therefore the
+output array of @code{gal_data_t} will contain @code{numlabs+1} elements.
+The first (zero-th) element of the output (@code{indexs[0]} in the example
+below) will be initialized to a dataset with zero elements. This will allow
+easy (non-confusing) access to the indexs of each (meaningful) label.
+
+@code{numlabs} is the number of labels in the dataset. If it is given a
+value of zero, then the maximum value in the input (largest label) will be
+found and used. Therefore if it is given, but smaller than the actual
+number of labels, this function may/will crash (it will write in
+un-allocated space). @code{numlabs} is therefore useful in a highly
+optimized/checked environment.
+
+For example, if the returned array is called @code{indexs}, then
+@code{indexs[10].size} contains the number of elements that have a label of
+@code{10} in @code{labels} and @code{indexs[10].array} is an array (after
+casting to @code{size_t *}) containing the indexs of each one of those
+elements/pixels.
+
+By @emph{index} we mean the 1D position: the input number of dimentions is
+irrelevant (any dimensionality is supported). In other words, each
+element's index is the number of elements/pixels between it and the
+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})
+@cindex Watershed algorithm
+@cindex Algorithm: watershed
+Use the watershed algorithm@footnote{The watershed algorithm was initially
+introduced by @url{https://doi.org/10.1109/34.87344, Vincent and
+Soille}. It starts from the minima and puts the pixels in, one by one, to
+grow them until the touch (create a watershed). For more, also see the
+Wikipedia article:
+@url{https://en.wikipedia.org/wiki/Watershed_%28image_processing%29}.} to
+``over-segment'' the pixels in the @code{indexs} dataset based on values in
+the @code{values} dataset. Internally, each local extrema (maximum or
+minimum, based on @code{min0_max1}) and its surrounding pixels will be
+given a unique label. For demonstration, see Figures 8 and 9 of
+@url{http://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa [2015]}. If
+@code{topinds!=NULL}, it is assumed to point to an already allocated space
+to write the indexs of the local extrema, otherwise, it is ignored.
+
+The @code{values} dataset must have a 32-bit floating point type
+(@code{GAL_TYPE_FLOAT32}, see @ref{Library data types}) and will only be
+read by this function. @code{indexs} must contain the indexs of the
+elements/pixels that will be over-segmented by this function and have a
+@code{GAL_TYPE_SIZE_T} type, see the description of
+@code{gal_label_indexs}, above. After this function, @code{indexs} will be
+sorted by the values of the respective pixel in @code{values}. The final
+labels will be written in the respective positions of @code{labels}, which
+must have a @code{GAL_TYPE_INT32} type and be the same size as
+@code{values}. All local minima (maxima), when @code{min0_max1} is @code{1}
+(@code{0}), are considered rivers (watersheds) and given a label of
+@code{GAL_LABEL_RIVER} (see above).
+@end deftypefun
+
+@deftypefun void gal_label_grow_indexs (gal_data_t @code{*labels}, gal_data_t 
@code{*indexs}, int @code{withrivers}, int @code{connectivity})
+Grow the (positive) labels of @code{labels} over the pixels in
+@code{indexs} (see description of @code{gal_label_indexs}). The pixels
+(position in @code{indexs}, values in @code{labels}) that must be ``grown''
+must have a value of @code{GAL_LABEL_INIT} in @code{labels} before calling
+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
+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
+pixel just get pushed back to the start of the loop, and the loop iterates
+until its size doesn't change any more. This is because in a generic
+scenario some of the indexed pixels might not be reachable through other
+indexed pixels.
+
+The next major difference with over-segmentation is that when there is only
+one label in growth region(s), it is not mandatory for @code{indexs} to be
+sorted by values. If there are multiple labeled regions in growth
+region(s), then values are important and you can use @code{qsort} with
+@code{gal_qsort_index_float_decreasing} to sort the indexs by values in a
+separate array (see @ref{Qsort functions}).
+
+This function looks for positive-valued neighbors of each pixel in
+@code{indexs} and will label a pixel if it touches one. Therefore, it is
+very important that only pixels/labels that are intended for growth have
+positive values in @code{labels} before calling this function. Any
+non-positive (zero or negative) value will be ignored as a label by this
+function. Thus, it is recommended that while filling in the `indexs' array
+values, you initialize all the pixels that are in @code{indexs} with
+@code{GAL_LABEL_INIT}, and set non-labeled pixels that you don't want to
+grow to @code{0}.
+
+This function will write into both the input datasets. After this function,
+some of the non-positive @code{labels} pixels will have a new positive
+label and the number of useful elements in @code{indexs} will have
+decreased. The index of those pixels that couldn't be labeled will remain
+inside @code{indexs}. If @code{withrivers} is non-zero, then pixels that
+are immediately touching more than one positive value will be given a
+@code{GAL_LABEL_RIVER} label.
+
+@cindex GNU C Library
+Note that the @code{indexs->array} is not re-allocated to its new size at
+the end@footnote{Note that according to the GNU C Library, even a
+@code{realloc} to a smaller size can also cause a re-write of the whole
+array, which is not a cheap operation.}. But since @code{indexs->dsize[0]}
+and @code{indexs->size} have new values after this function is returned,
+the extra elements just won't be used until they are ultimately freed by
+@code{gal_data_free}.
+
+Connectivity is a value between @code{1} (fewest number of neighbors) and
+the number of dimensions in the input (most number of neighbors). For
+example in a 2D dataset, a connectivity of @code{1} and @code{2}
+corresponds to 4-connected and 8-connected neighbors.
+@end deftypefun
+
+
+@node Convolution functions, Interpolation, Labeled datasets, Gnuastro library
 @subsection Convolution functions (@file{convolve.h})
 
 Convolution is a very common operation during data analysis and is
diff --git a/lib/gnuastro/label.h b/lib/gnuastro/label.h
index d494588..679202b 100644
--- a/lib/gnuastro/label.h
+++ b/lib/gnuastro/label.h
@@ -58,10 +58,12 @@ __BEGIN_C_DECLS  /* From C++ preparations */
 
 
 /* Functions. */
+gal_data_t *
+gal_label_indexs(gal_data_t *labels, size_t numlabs, size_t minmapsize);
+
 size_t
-gal_label_oversegment(gal_data_t *input, gal_data_t *indexs,
-                      gal_data_t *label, size_t *topinds,
-                      int min0_max1);
+gal_label_oversegment(gal_data_t *values, gal_data_t *indexs,
+                      gal_data_t *label, size_t *topinds, int min0_max1);
 
 void
 gal_label_grow_indexs(gal_data_t *labels, gal_data_t *indexs, int withrivers,
diff --git a/lib/gnuastro/statistics.h b/lib/gnuastro/statistics.h
index 55a732d..92b7868 100644
--- a/lib/gnuastro/statistics.h
+++ b/lib/gnuastro/statistics.h
@@ -57,17 +57,6 @@ __BEGIN_C_DECLS  /* From C++ preparations */
 
 
 
-/* Enumerators */
-enum is_sorted_outputs
-{
-  GAL_STATISTICS_SORTED_NOT,             /* ==0 by C standard. */
-
-  GAL_STATISTICS_SORTED_INVALID,
-  GAL_STATISTICS_SORTED_INCREASING,
-  GAL_STATISTICS_SORTED_DECREASING,
-};
-
-
 enum bin_status
 {
   GAL_STATISTICS_BINS_INVALID,           /* ==0 by C standard.  */
@@ -144,7 +133,7 @@ gal_statistics_mode_mirror_plots(gal_data_t *input, 
gal_data_t *value,
  ****************************************************************/
 
 int
-gal_statistics_is_sorted(gal_data_t *input);
+gal_statistics_is_sorted(gal_data_t *input, int updateflags);
 
 void
 gal_statistics_sort_increasing(gal_data_t *input);
diff --git a/lib/label.c b/lib/label.c
index 16178a7..6e719cb 100644
--- a/lib/label.c
+++ b/lib/label.c
@@ -25,12 +25,135 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <stdio.h>
 #include <errno.h>
 #include <error.h>
+#include <string.h>
 #include <stdlib.h>
 
 #include <gnuastro/list.h>
 #include <gnuastro/qsort.h>
 #include <gnuastro/label.h>
 #include <gnuastro/dimension.h>
+#include <gnuastro/statistics.h>
+
+
+
+
+
+
+
+
+
+/****************************************************************
+ *****************         Internal          ********************
+ ****************************************************************/
+static void
+label_check_type(gal_data_t *in, uint8_t needed_type, char *variable,
+                 const char *func)
+{
+  if(in->type!=needed_type)
+    error(EXIT_FAILURE, 0, "%s: the `%s' dataset has `%s' type, but it "
+          "must have a `%s' type.\n\n"
+          "You can use `gal_data_copy_to_new_type' or "
+          "`gal_data_copy_to_new_type_free' to convert your input dataset "
+          "to this type before calling this function", func, variable,
+          gal_type_name(in->type, 1), gal_type_name(needed_type, 1));
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/****************************************************************
+ *****************          Indexs           ********************
+ ****************************************************************/
+/* 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 *
+gal_label_indexs(gal_data_t *labels, size_t numlabs, size_t minmapsize)
+{
+  size_t i, *areas;
+  int32_t *a, *l, *lf;
+  gal_data_t *max, *labindexs;
+
+  /* Sanity check. */
+  label_check_type(labels, GAL_TYPE_INT32, "labels", __func__);
+
+  /* If the user hasn't given the number of labels, find it (maximum
+     label). */
+  if(numlabs==0)
+    {
+      max=gal_statistics_maximum(labels);
+      numlabs=*((int32_t *)(max->array));
+      gal_data_free(max);
+    }
+  labindexs=gal_data_array_calloc(numlabs+1);
+
+  /* Find the area in each detected object (to see how much space we need
+     to allocate). If blank values are present, an extra check is
+     necessary, so to get faster results when there aren't any blank
+     values, we'll also do a check. */
+  areas=gal_data_calloc_array(GAL_TYPE_SIZE_T, numlabs+1, __func__, "areas");
+  lf=(l=labels->array)+labels->size;
+  do
+    if(*l>0)  /* Only labeled regions: *l==0 (undetected), *l<0 (blank). */
+      ++areas[*l];
+  while(++l<lf);
+
+  /* For a check.
+  for(i=0;i<numlabs+1;++i)
+    printf("detection %zu: %zu\n", i, areas[i]);
+  exit(0);
+  */
+
+  /* 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.*/
+  for(i=1;i<numlabs+1;++i)
+    gal_data_initialize(&labindexs[i], NULL, GAL_TYPE_SIZE_T, 1,
+                        &areas[i], NULL, 0, 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, (numlabs+1)*sizeof *areas);
+  lf=(a=l=labels->array)+labels->size;
+  do
+    if(*l>0)  /* No undetected regions (*l==0), or blank (<0) */
+      ((size_t *)(labindexs[*l].array))[ areas[*l]++ ] = l-a;
+  while(++l<lf);
+
+  /* Clean up and return. */
+  free(areas);
+  return labindexs;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
 
 
 
@@ -53,17 +176,29 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 */
 size_t
-gal_label_oversegment(gal_data_t *input, gal_data_t *indexs,
-                      gal_data_t *label, size_t *topinds,
-                      int min0_max1)
+gal_label_oversegment(gal_data_t *values, gal_data_t *indexs,
+                      gal_data_t *labels, size_t *topinds, int min0_max1)
 {
-  size_t ndim=input->ndim;
+  size_t ndim=values->ndim;
 
-  float *arr=input->array;
+  float *arr=values->array;
   gal_list_sizet_t *Q=NULL, *cleanup=NULL;
-  size_t *a, *af, ind, *dsize=input->dsize;
+  size_t *a, *af, ind, *dsize=values->dsize;
   size_t *dinc=gal_dimension_increment(ndim, dsize);
-  int32_t n1, nlab, rlab, curlab=1, *clabel=label->array;
+  int32_t n1, nlab, rlab, curlab=1, *labs=labels->array;
+
+
+  /* Sanity checks */
+  label_check_type(values, GAL_TYPE_FLOAT32, "values", __func__);
+  label_check_type(indexs, GAL_TYPE_SIZE_T,  "indexs", __func__);
+  label_check_type(labels, GAL_TYPE_INT32,   "labels", __func__);
+  if( gal_data_dsize_is_different(values, labels) )
+    error(EXIT_FAILURE, 0, "%s: the `values' and `labels' arguments must "
+          "have the same size", __func__);
+  if(indexs->ndim!=1)
+    error(EXIT_FAILURE, 0, "%s: `indexs' has to be a 1D array, but it is "
+          "%zuD", __func__, indexs->ndim);
+
 
   /*********************************************
    For checks and debugging:*
@@ -75,22 +210,19 @@ gal_label_oversegment(gal_data_t *input, gal_data_t 
*indexs,
   char *filename="clumpbuild.fits";
   size_t checkstartind=gal_dimension_coord_to_index(2, dsize, checkstart);
   gal_data_t *tile=gal_data_alloc(gal_data_ptr_increment(arr, checkstartind,
-                                                         input->type),
+                                                         values->type),
                                   GAL_TYPE_INVALID, 2, checkdsize,
                                   NULL, 0, 0, NULL, NULL, NULL);
-  tile->block=input;
+  tile->block=values;
   gal_checkset_writable_remove(filename, 0, 0);
-  if(p->cp.numthreads!=1)
-    error(EXIT_FAILURE, 0, "in the debugging mode of `clumps_oversegment' "
-          "only one thread must be used");
   crop=gal_data_copy(tile);
   gal_fits_img_write(crop, filename, NULL, PROGRAM_NAME);
   gal_data_free(crop);
   printf("blank: %u\nriver: %u\ntmpcheck: %u\ninit: %u\n",
          (int32_t)GAL_BLANK_INT32, (int32_t)GAL_LABEL_RIVER,
          (int32_t)GAL_LABEL_TMPCHECK, (int32_t)GAL_LABEL_INIT);
-  tile->array=gal_tile_block_relative_to_other(tile, clabel);
-  tile->block=clabel;
+  tile->array=gal_tile_block_relative_to_other(tile, labels);
+  tile->block=labels;
   **********************************************/
 
 
@@ -100,7 +232,7 @@ gal_label_oversegment(gal_data_t *input, gal_data_t *indexs,
 
   /* Sort the given indexs based on their flux (`gal_qsort_index_arr' is
      defined as static in `gnuastro/qsort.h') */
-  gal_qsort_index_arr=input->array;
+  gal_qsort_index_arr=values->array;
   qsort(indexs->array, indexs->size, sizeof(size_t),
         min0_max1
         ? gal_qsort_index_float_decreasing
@@ -109,7 +241,7 @@ gal_label_oversegment(gal_data_t *input, gal_data_t *indexs,
 
   /* Initialize the region we want to over-segment. */
   af=(a=indexs->array)+indexs->size;
-  do clabel[*a]=GAL_LABEL_INIT; while(++a<af);
+  do labs[*a]=GAL_LABEL_INIT; while(++a<af);
 
 
   /* Go over all the given indexs and pull out the clumps. */
@@ -118,7 +250,7 @@ gal_label_oversegment(gal_data_t *input, gal_data_t *indexs,
     /* When regions of a constant flux or masked regions exist, some later
        indexs (although they have same flux) will be filled before hand. If
        they are done, there is no need to do them again. */
-    if(clabel[*a]==GAL_LABEL_INIT)
+    if(labs[*a]==GAL_LABEL_INIT)
       {
         /* It might happen where one or multiple regions of the pixels
            under study have the same flux. So two equal valued pixels of
@@ -148,7 +280,7 @@ gal_label_oversegment(gal_data_t *input, gal_data_t *indexs,
             /* Add this pixel to a queue. */
             gal_list_sizet_add(&Q, *a);
             gal_list_sizet_add(&cleanup, *a);
-            clabel[*a] = GAL_LABEL_TMPCHECK;
+            labs[*a] = GAL_LABEL_TMPCHECK;
 
             /* Find all the pixels that have the same flux and are
                connected. */
@@ -166,7 +298,7 @@ gal_label_oversegment(gal_data_t *input, gal_data_t *indexs,
                      if(n1!=GAL_LABEL_RIVER)
                        {
                          /* For easy reading. */
-                         nlab=clabel[ nind ];
+                         nlab=labs[ nind ];
 
                          /* This neighbor's label isn't zero. */
                          if(nlab)
@@ -176,7 +308,7 @@ gal_label_oversegment(gal_data_t *input, gal_data_t *indexs,
                                 to expand the studied region.*/
                              if( nlab==GAL_LABEL_INIT && arr[nind]==arr[*a] )
                                {
-                                 clabel[nind]=GAL_LABEL_TMPCHECK;
+                                 labs[nind]=GAL_LABEL_TMPCHECK;
                                  gal_list_sizet_add(&Q, nind);
                                  gal_list_sizet_add(&cleanup, nind);
                                }
@@ -204,7 +336,7 @@ gal_label_oversegment(gal_data_t *input, gal_data_t *indexs,
                                        outside this loop, for datasets with
                                        no blank element this last step will
                                        be completley ignored. */
-                                    : ( ( (input->flag
+                                    : ( ( (values->flag
                                            & GAL_DATA_FLAG_HASBLANK)
                                           && nlab==GAL_BLANK_INT32 )
                                         ? GAL_LABEL_RIVER : n1 ) );
@@ -214,10 +346,10 @@ gal_label_oversegment(gal_data_t *input, gal_data_t 
*indexs,
                             are on the edge of the indexed region (the
                             neighbor is not in the initial list of pixels
                             to segment). When over-segmenting the noise and
-                            the detections, `clabel' is zero for the parts
+                            the detections, `label' is zero for the parts
                             of the image that we are not interested in
                             here. */
-                         else clabel[*a]=GAL_LABEL_RIVER;
+                         else labs[*a]=GAL_LABEL_RIVER;
                        }
                    } );
               }
@@ -242,7 +374,7 @@ gal_label_oversegment(gal_data_t *input, gal_data_t *indexs,
                 ind=gal_list_sizet_pop(&cleanup);
                 /* If it was on the sides of the image, it has been
                    changed to a river pixel. */
-                if( clabel[ ind ]==GAL_LABEL_TMPCHECK ) clabel[ ind ]=rlab;
+                if( labs[ ind ]==GAL_LABEL_TMPCHECK ) labs[ ind ]=rlab;
               }
           }
 
@@ -268,7 +400,7 @@ gal_label_oversegment(gal_data_t *input, gal_data_t *indexs,
                  if(n1!=GAL_LABEL_RIVER)
                    {
                      /* For easy reading. */
-                     nlab=clabel[ nind ];
+                     nlab=labs[ nind ];
 
                      /* If this neighbor is on a non-processing label, then
                         set the first neighbor accordingly. Note that we
@@ -294,7 +426,7 @@ gal_label_oversegment(gal_data_t *input, gal_data_t *indexs,
                                    doesn't have blank values,
                                    `nlab==GAL_BLANK_INT32' will never be
                                    checked. */
-                                : ( (input->flag & GAL_DATA_FLAG_HASBLANK)
+                                : ( (values->flag & GAL_DATA_FLAG_HASBLANK)
                                     && nlab==GAL_BLANK_INT32
                                     ? GAL_LABEL_RIVER : n1 ) )
 
@@ -322,7 +454,7 @@ gal_label_oversegment(gal_data_t *input, gal_data_t *indexs,
               }
 
             /* Put the found label in the pixel. */
-            clabel[ *a ] = rlab;
+            labs[ *a ] = rlab;
           }
 
         /*********************************************
@@ -334,7 +466,7 @@ gal_label_oversegment(gal_data_t *input, gal_data_t *indexs,
           {
             printf("%zu (%zu: %zu, %zu): %u\n", ++extcount, *a,
                    (*a%dsize[1])-checkstart[1], (*a/dsize[1])-checkstart[0],
-                   clabel[*a]);
+                   labs[*a]);
             crop=gal_data_copy(tile);
             crf=(cr=crop->array)+crop->size;
             do if(*cr==GAL_LABEL_RIVER) *cr=0; while(++cr<crf);
@@ -363,48 +495,7 @@ gal_label_oversegment(gal_data_t *input, gal_data_t 
*indexs,
 
 
 
-/* Grow the given labels over the list of indexs. In over-segmentation
-   (watershed algorithm), the indexs have to be sorted by pixel value and
-   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 pixel just
-   get pushed back to the start of the loop, and the loop iterates until
-   its size doesn't change any more. This is because in a generic scenario
-   some of the indexed pixels might not be reachable. As a result, it is
-   not necessary for the given indexs to be sorted here.
-
-   This function looks for positive-valued neighbors and will label a pixel
-   if it touches one. Therefore, it is is very important that only
-   pixels/labels that are intended for growth have positive values before
-   calling this function. Any non-positive (zero or negative) value will be
-   ignored by this function. Thus, it is recommended that while filling in
-   the `indexs' array values, you initialize all those pixels with
-   `GAL_LABEL_INIT', and set non-labeled pixels that you don't want to grow
-   to 0.
-
-   This function will write into both the input datasets. After this
-   function, some of the non-positive `labels' pixels will have a new
-   positive label and the number of useful elements in `indexs' will have
-   decreased. Those indexs that couldn't be labeled will remain inside
-   `indexs'. If `withrivers' is non-zero, then pixels that are immediately
-   touching more than one positive value will be given a `GAL_LABEL_RIVER'
-   label.
-
-   Note that the `indexs->array' is not re-allocated to its new size at the
-   end. But since `indexs->dsize[0]' and `indexs->size' have new values,
-   the extra elements just won't be used until they are ultimately freed by
-   `gal_data_free'.
-
-   Input:
-
-     labels: The labels array that must be operated on. The pixels that
-             must be "grown" must have the value `GAL_LABEL_INIT' (negative).
-
-     sorted_indexs: The indexs of the pixels that must be grown.
-
-     withrivers: as described above.
-
-     connectivity: connectivity to define neighbors for growth.  */
+/* 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)
@@ -416,12 +507,8 @@ gal_label_grow_indexs(gal_data_t *labels, gal_data_t 
*indexs, int withrivers,
   size_t *dinc=gal_dimension_increment(labels->ndim, labels->dsize);
 
   /* Some basic sanity checks: */
-  if(labels->type!=GAL_TYPE_INT32)
-    error(EXIT_FAILURE, 0, "%s: `labels' has to have type of int32_t",
-          __func__);
-  if(indexs->type!=GAL_TYPE_SIZE_T)
-    error(EXIT_FAILURE, 0, "%s: `indexs' must be `size_t' type but it is "
-          "`%s'", __func__, gal_type_name(indexs->type,1));
+  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);
diff --git a/lib/statistics.c b/lib/statistics.c
index 2d0c0d7..34630bc 100644
--- a/lib/statistics.c
+++ b/lib/statistics.c
@@ -1153,32 +1153,48 @@ gal_statistics_mode_mirror_plots(gal_data_t *input, 
gal_data_t *value,
 /****************************************************************
  ********                      Sort                       *******
  ****************************************************************/
-/* Check if the given dataset is sorted. Output values are:
-
-     - 0: Dataset is not sorted.
-     - 1: Dataset is sorted and increasing or equal.
-     - 2: dataset is sorted and decreasing.                  */
+/* Check if the given dataset is sorted. */
+enum is_sorted_return
+{
+  STATISTICS_IS_SORTED_NOT,                 /* ==0: by C standard. */
+  STATISTICS_IS_SORTED_INCREASING,
+  STATISTICS_IS_SORTED_DECREASING,
+};
 
 #define IS_SORTED(IT) {                                                 \
-    IT *aa=input->array, *a=input->array, *af=a+input->size-1;          \
-    if(a[1]>=a[0]) do if( *(a+1) < *a ) break; while(++a<af);           \
-    else           do if( *(a+1) > *a ) break; while(++a<af);           \
-    return ( a==af                                                      \
-             ? ( aa[1]>=aa[0]                                           \
-                 ? GAL_STATISTICS_SORTED_INCREASING                     \
-                 : GAL_STATISTICS_SORTED_DECREASING )                   \
-             : GAL_STATISTICS_SORTED_NOT );                             \
+  IT *aa=input->array, *a=input->array, *af=a+input->size-1;            \
+  if(a[1]>=a[0]) do if( *(a+1) < *a ) break; while(++a<af);             \
+  else           do if( *(a+1) > *a ) break; while(++a<af);             \
+  out=( a==af                   /* It reached the end of the array. */  \
+          ? ( aa[1]>=aa[0]                                              \
+                ? STATISTICS_IS_SORTED_INCREASING                       \
+                : STATISTICS_IS_SORTED_DECREASING )                     \
+          : STATISTICS_IS_SORTED_NOT );                                 \
   }
 
 int
-gal_statistics_is_sorted(gal_data_t *input)
+gal_statistics_is_sorted(gal_data_t *input, int updateflags)
 {
-  /* A one-element dataset can be considered, sorted, so we'll just return
-     1 (for sorted and increasing). */
+  int out;
+
+  /* If the flags are already set, don't bother going over the dataset. */
+  if( input->flag & GAL_DATA_FLAG_SORT_CH )
+    return ( input->flag & GAL_DATA_FLAG_SORTED_I
+             ? STATISTICS_IS_SORTED_INCREASING
+             : ( input->flag & GAL_DATA_FLAG_SORTED_D
+                 ? STATISTICS_IS_SORTED_DECREASING
+                 : STATISTICS_IS_SORTED_NOT ) );
+
+  /* Parse the array (if necessary). */
   switch(input->size)
     {
-    case 0: return GAL_STATISTICS_SORTED_INVALID;
-    case 1: return GAL_STATISTICS_SORTED_INCREASING;
+    case 0:
+      error(EXIT_FAILURE, 0, "%s: input dataset has 0 elements", __func__);
+
+    /* A one-element dataset can be considered, sorted, so we'll say its
+       increasing. */
+    case 1:
+      out=STATISTICS_IS_SORTED_INCREASING;
 
     /* Do the check. */
     default:
@@ -1200,11 +1216,34 @@ gal_statistics_is_sorted(gal_data_t *input)
         }
     }
 
-  /* Control shouldn't reach this point. */
-  error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s so we can "
-        "fix the problem. Control must not have reached the end of this "
-        "function", __func__, PACKAGE_BUGREPORT);
-  return GAL_STATISTICS_SORTED_INVALID;
+  /* Update the flags, if required. */
+  if(updateflags)
+    {
+      input->flag |= GAL_DATA_FLAG_SORT_CH;
+      switch(out)
+        {
+        case STATISTICS_IS_SORTED_NOT:
+          input->flag &= ~GAL_DATA_FLAG_SORTED_I;
+          input->flag &= ~GAL_DATA_FLAG_SORTED_D;
+          break;
+
+        case STATISTICS_IS_SORTED_INCREASING:
+          input->flag |=  GAL_DATA_FLAG_SORTED_I;
+          input->flag &= ~GAL_DATA_FLAG_SORTED_D;
+          break;
+
+        case STATISTICS_IS_SORTED_DECREASING:
+          input->flag &= ~GAL_DATA_FLAG_SORTED_I;
+          input->flag |=  GAL_DATA_FLAG_SORTED_D;
+          break;
+
+        default:
+          error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix "
+                "the problem. The value %d is not recognized for `out'",
+                __func__, PACKAGE_BUGREPORT, out);
+        }
+    }
+  return out;
 }
 
 
@@ -1219,6 +1258,7 @@ gal_statistics_is_sorted(gal_data_t *input)
 void
 gal_statistics_sort_increasing(gal_data_t *input)
 {
+  /* Do the sorting. */
   if(input->size)
     switch(input->type)
       {
@@ -1246,6 +1286,11 @@ gal_statistics_sort_increasing(gal_data_t *input)
         error(EXIT_FAILURE, 0, "%s: type code %d not recognized",
               __func__, input->type);
       }
+
+  /* Set the flags. */
+  input->flag |=  GAL_DATA_FLAG_SORT_CH;
+  input->flag |=  GAL_DATA_FLAG_SORTED_I;
+  input->flag &= ~GAL_DATA_FLAG_SORTED_D;
 }
 
 
@@ -1256,6 +1301,7 @@ gal_statistics_sort_increasing(gal_data_t *input)
 void
 gal_statistics_sort_decreasing(gal_data_t *input)
 {
+  /* Do the sorting. */
   if(input->size)
     switch(input->type)
       {
@@ -1283,6 +1329,11 @@ gal_statistics_sort_decreasing(gal_data_t *input)
         error(EXIT_FAILURE, 0, "%s: type code %d not recognized",
               __func__, input->type);
       }
+
+  /* Set the flags. */
+  input->flag |=  GAL_DATA_FLAG_SORT_CH;
+  input->flag |=  GAL_DATA_FLAG_SORTED_D;
+  input->flag &= ~GAL_DATA_FLAG_SORTED_I;
 }
 
 
@@ -1301,7 +1352,6 @@ gal_statistics_sort_decreasing(gal_data_t *input)
 gal_data_t *
 gal_statistics_no_blank_sorted(gal_data_t *input, int inplace)
 {
-  int sortstatus;
   gal_data_t *contig, *noblank, *sorted;
 
   /* We need to account for the case that there are no elements in the
@@ -1324,23 +1374,15 @@ gal_statistics_no_blank_sorted(gal_data_t *input, int 
inplace)
       else contig=input;
 
 
-      /* Make sure there is no blanks in the array that will be used. After
-         this step, we won't be dealing with `input' any more, but with
-         `noblank'. */
-      if( gal_blank_present(contig, inplace) )
+      /* Make sure there are no blanks in the array that will be
+         used. After this step, we won't be dealing with `input' any more,
+         but with `noblank'. */
+      if( gal_blank_present(contig, 1) )
         {
           /* See if we should allocate a new dataset to remove blanks or if
              we can use the actual contiguous patch of memory. */
           noblank = inplace ? contig : gal_data_copy(contig);
           gal_blank_remove(noblank);
-
-          /* If we are working in place, then mark that there are no blank
-             pixels. */
-          if(inplace)
-            {
-              noblank->flag |= GAL_DATA_FLAG_BLANK_CH;
-              noblank->flag &= ~GAL_DATA_FLAG_HASBLANK;
-            }
         }
       else noblank=contig;
 
@@ -1350,12 +1392,8 @@ gal_statistics_no_blank_sorted(gal_data_t *input, int 
inplace)
          more but with `sorted'. */
       if(noblank->size)
         {
-          sortstatus=gal_statistics_is_sorted(noblank);
-          if( sortstatus )
-            {
-              sorted=noblank;
-              sorted->status=sortstatus;
-            }
+          if( gal_statistics_is_sorted(noblank, 1) )
+            sorted=noblank;
           else
             {
               if(inplace) sorted=noblank;
@@ -1367,7 +1405,6 @@ gal_statistics_no_blank_sorted(gal_data_t *input, int 
inplace)
                     sorted=gal_data_copy(noblank);
                 }
               gal_statistics_sort_increasing(sorted);
-              sorted->status=GAL_STATISTICS_SORTED_INCREASING;
             }
         }
       else
@@ -1847,7 +1884,7 @@ gal_statistics_cfp(gal_data_t *input, gal_data_t *bins, 
int normalize)
     IT *bf = nbs->array, *b  = bf + nbs->size - 1;                      \
                                                                         \
     /* Remove all out-of-range elements from the start of the array. */ \
-    if(sortstatus==GAL_STATISTICS_SORTED_INCREASING)                    \
+    if( nbs->flag & GAL_DATA_FLAG_SORTED_I )                            \
       do if( *a > (*med - (multip * *std)) )                            \
            { start=a; break; }                                          \
       while(++a<af);                                                    \
@@ -1857,7 +1894,7 @@ gal_statistics_cfp(gal_data_t *input, gal_data_t *bins, 
int normalize)
       while(++a<af);                                                    \
                                                                         \
     /* Remove all out-of-range elements from the end of the array. */   \
-    if(sortstatus==GAL_STATISTICS_SORTED_INCREASING)                    \
+    if( nbs->flag & GAL_DATA_FLAG_SORTED_I )                            \
       do if( *b < (*med + (multip * *std)) )                            \
            { size=b-a+1; break; }                                       \
       while(--b>=bf);                                                   \
@@ -1874,11 +1911,11 @@ gal_statistics_sigma_clip(gal_data_t *input, float 
multip, float param,
   float *oa;
   void *start, *nbs_array;
   double *med, *mean, *std;
+  uint8_t type=gal_tile_block(input)->type;
   uint8_t bytolerance = param>=1.0f ? 0 : 1;
   double oldmed=NAN, oldmean=NAN, oldstd=NAN;
   size_t num=0, one=1, four=4, size, oldsize;
   gal_data_t *median_i, *median_d, *out, *meanstd;
-  int sortstatus, type=gal_tile_block(input)->type;
   gal_data_t *nbs=gal_statistics_no_blank_sorted(input, inplace);
   size_t maxnum = param>=1.0f ? param : GAL_STATISTICS_SIG_CLIP_MAX_CONVERGE;
 
@@ -1894,6 +1931,14 @@ gal_statistics_sigma_clip(gal_data_t *input, float 
multip, float param,
     error(EXIT_FAILURE, 0, "%s: when `param' is larger than 1.0, it is "
           "interpretted as an absolute number of clips. So it must be an "
           "integer. However, your given value %g", __func__, param);
+  if( (nbs->flag & GAL_DATA_FLAG_SORT_CH)==0 )
+    error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix the "
+          "problem. `nbs->flag', doesn't have the `GAL_DATA_FLAG_SORT_CH' "
+          "bit activated", __func__, PACKAGE_BUGREPORT);
+  if( (nbs->flag & GAL_DATA_FLAG_SORTED_I)==0
+      && (nbs->flag & GAL_DATA_FLAG_SORTED_D)==0 )
+    error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix the "
+          "problem. `nbs' isn't sorted", __func__, PACKAGE_BUGREPORT);
 
 
   /* Allocate the necessary spaces. */
@@ -1926,7 +1971,6 @@ gal_statistics_sigma_clip(gal_data_t *input, float 
multip, float param,
          array's size. */
       size=nbs->size;
       start=nbs->array;
-      sortstatus=nbs->status;
       while(num<maxnum && size)
         {
           /* Find the median. */