[gnuastro-commits] master 188a4f4: Clumps can be defined based on local minima in Segment

[Mohammad Akhlaghi]

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

    Clumps can be defined based on local minima in Segment
    
    Until now clumps could only be defined based on local maxima. With this
    commit, `gal_label_oversegment' can also do its oversegmentation starting
    from local minima. As a result, Segment can also do the same. The
    respective parts of the S/N calculation were also corrected to allow this.
    
    To do this, a new `gal_qsort_index_float_increasing' function was also
    defined in the `qsort' library.
---
 NEWS                  | 10 +++++++--
 bin/segment/args.h    | 13 +++++++++++
 bin/segment/clumps.c  |  8 +++----
 bin/segment/main.h    |  1 +
 bin/segment/segment.c |  3 ++-
 bin/segment/ui.h      |  1 +
 doc/gnuastro.texi     | 60 ++++++++++++++++++++++++++++++++-------------------
 lib/gnuastro/label.h  |  3 ++-
 lib/gnuastro/qsort.h  |  3 ++-
 lib/label.c           |  7 ++++--
 lib/qsort.c           |  8 +++++++
 11 files changed, 84 insertions(+), 33 deletions(-)

diff --git a/NEWS b/NEWS
index 129095b..f0014f8 100644
--- a/NEWS
+++ b/NEWS
@@ -36,7 +36,7 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
     --instd: Sky standard deviation as a single value or file (dataset).
     --valuesfile: filename containing the values dataset.
     --valueshdu: HDU/extension containing the values dataset.
-    --clumpscat: Make a clumps catalog also (`WCLUMPS' keyword not used 
anymore).
+    --clumpscat: Make a clumps catalog also.
     --noclumpsort: Don't sort the clumps catalog by host object ID.
     --subtractsky: Subtract the given Sky from the values dataset.
     --variance: input standard deviation image is actually variance.
@@ -66,7 +66,7 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
   Libraries:
     gal_array_read: read array from any of known formats (FITS,TIFF,JPEG,...).
     gal_array_read_to_type: similar to `gal_array_read', but to given type.
-    gal_array_read_one_ch: Read in the data and make sure it is in one channel.
+    gal_array_read_one_ch: Read a dataset, abort if it has multiple channels.
     gal_array_read_one_ch_to_type: Make sure input is in one channel and type.
     gal_binary_label_holes: label the holes within the foreground.
     gal_blank_is: check to see if argument is blank in its type or not.
@@ -83,6 +83,7 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
     gal_pdf_name_is_pdf: Returns 1 if given filename is PDF.
     gal_pdf_suffix_is_pdf: Returns 1 if given suffix is PDF.
     gal_pdf_write: Writes a dataset into an PDF file.
+    gal_qsort_index_float_increasing: Sort indexs in increasing value order.
     gal_tiff_name_is_tiff: check if name contains a TIFF suffix.
     gal_tiff_suffix_is_tiff: check if suffix is a TIFF suffix.
     gal_tiff_dir_string_read: convert a string to a TIFF directory number.
@@ -124,6 +125,10 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
     - With no output name, the output has a `_detected.fits' suffix.
 
   MakeCatalog:
+
+    - The `WCLUMPS' keyword in the objects labeled image is no longer used
+         to see if a clumps catalog should also be made. To build a clumps
+         catalog, you can now use the `--clumpscat' option.
     - Estimation of noise-level is now done per-pixel over the whole
          label. Until now the average noise level was used.
     --objectsfile has been removed. The main input argument is now assumed
@@ -171,6 +176,7 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
 
 
 
+
 * Noteworthy changes in release 0.5 (library 3.0.0) (2017-12-22) [stable]
 
 ** New features
diff --git a/bin/segment/args.h b/bin/segment/args.h
index 280a0cc..6293575 100644
--- a/bin/segment/args.h
+++ b/bin/segment/args.h
@@ -264,6 +264,19 @@ struct argp_option program_options[] =
       GAL_OPTIONS_NOT_SET
     },
     {
+      "minima",
+      UI_KEY_MINIMA,
+      0,
+      0,
+      "Built internal clumps from minima.",
+      UI_GROUP_SEGMENTATION,
+      &p->minima,
+      GAL_OPTIONS_NO_ARG_TYPE,
+      GAL_OPTIONS_RANGE_0_OR_1,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
+    },
+    {
       "snminarea",
       UI_KEY_SNMINAREA,
       "INT",
diff --git a/bin/segment/clumps.c b/bin/segment/clumps.c
index 15074ad..83a0204 100644
--- a/bin/segment/clumps.c
+++ b/bin/segment/clumps.c
@@ -29,7 +29,6 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <string.h>
 
 #include <gnuastro/fits.h>
-#include <gnuastro/qsort.h>
 #include <gnuastro/blank.h>
 #include <gnuastro/label.h>
 #include <gnuastro/threads.h>
@@ -423,7 +422,7 @@ clumps_make_sn_table(struct clumps_thread_params *cltprm)
          noise cases) or the area is smaller than the minimum area to
          calculate signal-to-noise, then set the S/N of this segment to
          zero. */
-      if( I>O && Ni>p->snminarea )
+      if( (p->minima ? O>I : I>O) && Ni>p->snminarea )
         {
           /* For easy reading, define `var' for variance.  */
           var = row[INFO_INSTD] * row[INFO_INSTD];
@@ -435,7 +434,7 @@ clumps_make_sn_table(struct clumps_thread_params *cltprm)
              equal to i. */
           ind = sky0_det1 ? i : counter++;
           if(cltprm->snind) indarr[ind]=i;
-          snarr[ind]=( sqrt(Ni/p->cpscorr) * (I-O)
+          snarr[ind]=( sqrt(Ni/p->cpscorr) * ( p->minima ? O-I : I-O)
                        / sqrt( (I>0?I:-1*I) + (O>0?O:-1*O) + var ) );
         }
       else
@@ -686,7 +685,8 @@ 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);
+                                                     cltprm.topinds,
+                                                     !p->minima);
 
 
           /* Set all river pixels to GAL_LABEL_INIT (to be distinguishable
diff --git a/bin/segment/main.h b/bin/segment/main.h
index 745182f..083045e 100644
--- a/bin/segment/main.h
+++ b/bin/segment/main.h
@@ -60,6 +60,7 @@ struct segmentparams
   char               *stdname;  /* File name of Standard deviation image. */
   char                *stdhdu;  /* HDU of Stanard deviation image.        */
   uint8_t            variance;  /* The input STD is actually variance.    */
+  uint8_t              minima;  /* Build clumps from their minima, maxima.*/
   uint8_t           rawoutput;  /* Output only object and clump labels.   */
 
   float            minskyfrac;  /* Undetected area min. frac. in tile.    */
diff --git a/bin/segment/segment.c b/bin/segment/segment.c
index 212b043..deae6fb 100644
--- a/bin/segment/segment.c
+++ b/bin/segment/segment.c
@@ -542,7 +542,8 @@ 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->clabel, cltprm.topinds,
+                                                 !p->minima);
 
 
       /* Set all the river pixels to zero (we don't need them any more in
diff --git a/bin/segment/ui.h b/bin/segment/ui.h
index 79a2cc5..d591221 100644
--- a/bin/segment/ui.h
+++ b/bin/segment/ui.h
@@ -79,6 +79,7 @@ enum option_keys_enum
   UI_KEY_STD,
   UI_KEY_STDHDU,
   UI_KEY_VARIANCE,
+  UI_KEY_MINIMA,
   UI_KEY_RAWOUTPUT,
   UI_KEY_MINNUMFALSE,
   UI_KEY_GROWNCLUMPS,
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 2d51268..493fb86 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -14731,12 +14731,14 @@ many illustrative figures) is Section 3.2 of
 @url{https://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa [2015]}.
 
 As a summary, Segment first finds true @emph{clump}s over the
-detections. Clumps are associated with local maxima and extend over the
-neighboring pixels until they reach a local minimum (@emph{river}). Segment
-will use the distribution of clump signal-to-noise ratios over the
-undetected regions as reference to find ``true'' clumps over the
-detections. Using the undetected regions can be disabled by directly giving
-a signal-to-noise ratio to @option{--clumpsnthresh}.
+detections. Clumps are associated with local maxima/address@hidden
+default the maximum is used as the first clump pixel, to define clumps
+based on local minima, use the @option{--minima} option.})  and extend over
+the neighboring pixels until they reach a local minimum/maximum
+(@emph{river}). Segment will use the distribution of clump signal-to-noise
+ratios over the undetected regions as reference to find ``true'' clumps
+over the detections. Using the undetected regions can be disabled by
+directly giving a signal-to-noise ratio to @option{--clumpsnthresh}.
 
 The true clumps are then grown to a certain threshold over the
 detections. Based on the strength of the connections between the grown
@@ -15054,6 +15056,13 @@ fields. Operationally, this is almost identical to 
NoiseChisel's
 @option{--minskyfrac} option (@ref{Detection options}). Please see the
 descriptions there for more.
 
address@hidden --minima
+Build the clumps based on the local minima, not maxima. By default, clumps
+are built starting from local maxima (see Figure 8 of
address@hidden://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa
+[2015]}). Therefore, this option can be useful when you are searching for
+true local minima (for example absorption features).
+
 @item -m INT
 @itemx --snminarea=INT
 The minimum area which a clump in the undetected regions should have in
@@ -15103,14 +15112,14 @@ option and visually inspected with 
@option{--checksegmentation}.
 
 @item -v
 @itemx --keepmaxnearriver
-Keep a clump whose maximum flux is 8-connected to a river pixel. By default
-such clumps over detections are considered to be noise and are removed
-irrespective of their brightness (see @ref{Flux Brightness and
-magnitude}). Over large profiles, that sink into the noise very slowly,
-noise can cause part of the profile (which was flat without noise) to
-become a very large and with a very high Signal to noise ratio. In such
-cases, the pixel with the maximum flux in the clump will be immediately
-touching a river pixel.
+Keep a clump whose maximum (minimum if @option{--minima} is called) flux is
+8-connected to a river pixel. By default such clumps over detections are
+considered to be noise and are removed irrespective of their brightness
+(see @ref{Flux Brightness and magnitude}). Over large profiles, that sink
+into the noise very slowly, noise can cause part of the profile (which was
+flat without noise) to become a very large and with a very high Signal to
+noise ratio. In such cases, the pixel with the maximum flux in the clump
+will be immediately touching a river pixel.
 
 @item -s FLT
 @itemx --clumpsnthresh=FLT
@@ -15188,14 +15197,15 @@ definition of clumps and objects, please see Section 
3.2 of
 @url{https://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa [2015]} and
 @ref{Segmentation options}.
 
-The clumps are ``true'' local maxima and their surrounding pixels until a
-local minimum (caused by noise fluctuations, or another ``true''
-clump). Therefore it may happen that some of the input detections aren't
-covered by clumps at all (very diffuse objects without any strong peak),
-while some objects may contain many clumps. Even in those that have clumps,
-there will be regions that are too diffuse. The diffuse regions (within the
-input detected regions) are given a negative label (-1) to help you
-separate them from the undetected regions (with a value of zero).
+The clumps are ``true'' local maxima (minima if @option{--minima} is
+called) and their surrounding pixels until a local minimum/maximum (caused
+by noise fluctuations, or another ``true'' clump). Therefore it may happen
+that some of the input detections aren't covered by clumps at all (very
+diffuse objects without any strong peak), while some objects may contain
+many clumps. Even in those that have clumps, there will be regions that are
+too diffuse. The diffuse regions (within the input detected regions) are
+given a negative label (-1) to help you separate them from the undetected
+regions (with a value of zero).
 
 Each clump is labeled with respect to its host object. Therefore, if an
 object has three clumps for example, the clumps within it have labels 1, 2
@@ -24603,6 +24613,12 @@ main (void)
 The output will be: @code{2, 0, 1, 3}.
 @end deftypefun
 
address@hidden int gal_qsort_index_float_increasing (const void @code{*a}, 
const void @code{*b})
+Similar to @code{gal_qsort_index_float_decreasing}, but will sort the
+indexes such that the values of @code{gal_qsort_index_arr} are in
+increasing order.
address@hidden deftypefun
+
 
 @deftypefun int gal_qsort_TYPE_increasing (const void @code{*a}, const void 
@code{*b})
 When passed to @code{qsort}, this function will sort an @code{TYPE} array
diff --git a/lib/gnuastro/label.h b/lib/gnuastro/label.h
index e73b7d8..d494588 100644
--- a/lib/gnuastro/label.h
+++ b/lib/gnuastro/label.h
@@ -60,7 +60,8 @@ __BEGIN_C_DECLS  /* From C++ preparations */
 /* Functions. */
 size_t
 gal_label_oversegment(gal_data_t *input, gal_data_t *indexs,
-                      gal_data_t *label, size_t *topinds);
+                      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/qsort.h b/lib/gnuastro/qsort.h
index 5daa693..2b20876 100644
--- a/lib/gnuastro/qsort.h
+++ b/lib/gnuastro/qsort.h
@@ -55,7 +55,8 @@ extern float *gal_qsort_index_arr;
 int
 gal_qsort_index_float_decreasing(const void * a, const void * b);
 
-
+int
+gal_qsort_index_float_increasing(const void * a, const void * b);
 
 
 
diff --git a/lib/label.c b/lib/label.c
index 538f704..16178a7 100644
--- a/lib/label.c
+++ b/lib/label.c
@@ -54,7 +54,8 @@ 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)
+                      gal_data_t *label, size_t *topinds,
+                      int min0_max1)
 {
   size_t ndim=input->ndim;
 
@@ -101,7 +102,9 @@ gal_label_oversegment(gal_data_t *input, gal_data_t *indexs,
      defined as static in `gnuastro/qsort.h') */
   gal_qsort_index_arr=input->array;
   qsort(indexs->array, indexs->size, sizeof(size_t),
-        gal_qsort_index_float_decreasing);
+        min0_max1
+        ? gal_qsort_index_float_decreasing
+        : gal_qsort_index_float_increasing );
 
 
   /* Initialize the region we want to over-segment. */
diff --git a/lib/qsort.c b/lib/qsort.c
index 6e6cd70..9800b5c 100644
--- a/lib/qsort.c
+++ b/lib/qsort.c
@@ -40,6 +40,14 @@ gal_qsort_index_float_decreasing(const void * a, const void 
* b)
   return (tb > ta) - (tb < ta);
 }
 
+int
+gal_qsort_index_float_increasing(const void * a, const void * b)
+{
+  float ta=gal_qsort_index_arr[ *(size_t *)a ];
+  float tb=gal_qsort_index_arr[ *(size_t *)b ];
+  return (ta > tb) - (ta < tb);
+}
+