[gnuastro-commits] master 1e26e7b 1/2: Arithmetic can now find the connected components of input

[Mohammad Akhlaghi]

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

    Arithmetic can now find the connected components of input
    
    Thanks to the library, adding this new operator was very
    straightforward. With it, it is easily possible to produce labeled images
    when a simple threshold is enough for a detection problem. These labeled
    images can then be fed into Segment for segmentation or MakeCatalog for
    higher-level analysis.
---
 NEWS                        |  1 +
 bin/arithmetic/arithmetic.c | 72 +++++++++++++++++++++++++++++++++++++++++++++
 bin/arithmetic/arithmetic.h |  1 +
 doc/gnuastro.texi           | 31 +++++++++++++++++++
 4 files changed, 105 insertions(+)

diff --git a/NEWS b/NEWS
index e34dbe0..e50e6e3 100644
--- a/NEWS
+++ b/NEWS
@@ -13,6 +13,7 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
     - Segment: new program in charge of segmentation over detections.
 
   Arithmetic:
+    - connected-components: label the connected elements of the input.
     - filter-sigclip-mean: sigma-clipped, mean filter operator.
     - filter-sigclip-median: sigma-clipped, median filter operator.
 
diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index 7be17cc..54729e2 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -32,6 +32,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/wcs.h>
 #include <gnuastro/fits.h>
 #include <gnuastro/array.h>
+#include <gnuastro/binary.h>
 #include <gnuastro/threads.h>
 #include <gnuastro/dimension.h>
 #include <gnuastro/statistics.h>
@@ -510,6 +511,71 @@ wrapper_for_filter(struct arithmeticparams *p, char 
*token, int operator)
 
 
 /***************************************************************/
+/*************            Other functions          *************/
+/***************************************************************/
+static void
+arithmetic_connected_components(struct arithmeticparams *p, char *token)
+{
+  int conn_int;
+  gal_data_t *out=NULL;
+
+  /* Pop the two necessary operands. */
+  gal_data_t *conn = operands_pop(p, token);
+  gal_data_t *in   = operands_pop(p, token);
+
+  /* Do proper sanity checks on `conn'. */
+  if(conn->size!=1)
+    error(EXIT_FAILURE, 0, "the first popped operand to "
+          "`connected-components' must be a single number. However, it has "
+          "%zu elements", conn->size);
+  if(conn->type==GAL_TYPE_FLOAT32 || conn->type==GAL_TYPE_FLOAT64)
+    error(EXIT_FAILURE, 0, "the first popped operand to "
+          "`connected-components' is the connectivity (a value between 1 "
+          "and the number of dimensions) therefore, it must NOT be a "
+          "floating point");
+
+  /* Convert the connectivity value to a 32-bit integer and read it in and
+     make sure it is not larger than the number of dimensions. */
+  conn=gal_data_copy_to_new_type_free(conn, GAL_TYPE_INT32);
+  conn_int = *((int32_t *)(conn->array));
+  if(conn_int>in->ndim)
+    error(EXIT_FAILURE, 0, "the first0popped operand to "
+          "`connected-components' (%d) is larger than the number of "
+          "dimensions in the second-popped operand (%zu)", conn_int,
+          in->ndim);
+
+  /* Do the connected components labeling. */
+  gal_binary_connected_components(in, &out, 1);
+
+  /* Push the result onto the stack. */
+  operands_add(p, NULL, out);
+
+  /* Clean up. */
+  gal_data_free(in);
+  gal_data_free(conn);
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/***************************************************************/
 /*************      Reverse Polish algorithm       *************/
 /***************************************************************/
 /* This function implements the reverse polish algorithm as explained
@@ -671,6 +737,8 @@ reversepolish(struct arithmeticparams *p)
             { op=ARITHMETIC_OP_FILTER_SIGCLIP_MEAN;   nop=0;  }
           else if (!strcmp(token->v, "filter-sigclip-median"))
             { op=ARITHMETIC_OP_FILTER_SIGCLIP_MEDIAN; nop=0;  }
+          else if (!strcmp(token->v, "connected-components"))
+            { op=ARITHMETIC_OP_CONNECTED_COMPONENTS;  nop=0;  }
 
 
           /* Finished checks with known operators */
@@ -745,6 +813,10 @@ reversepolish(struct arithmeticparams *p)
                   wrapper_for_filter(p, token->v, op);
                   break;
 
+                case ARITHMETIC_OP_CONNECTED_COMPONENTS:
+                  arithmetic_connected_components(p, token->v);
+                  break;
+
                 default:
                   error(EXIT_FAILURE, 0, "%s: a bug! please contact us at "
                         "%s to fix the problem. The code %d is not "
diff --git a/bin/arithmetic/arithmetic.h b/bin/arithmetic/arithmetic.h
index 5b52da8..da37fa1 100644
--- a/bin/arithmetic/arithmetic.h
+++ b/bin/arithmetic/arithmetic.h
@@ -35,6 +35,7 @@ enum arithmetic_prog_operators
   ARITHMETIC_OP_FILTER_MEAN,
   ARITHMETIC_OP_FILTER_SIGCLIP_MEAN,
   ARITHMETIC_OP_FILTER_SIGCLIP_MEDIAN,
+  ARITHMETIC_OP_CONNECTED_COMPONENTS,
 };
 
 
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index c2514f8..f5b11f7 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -10089,6 +10089,37 @@ dataset. This operator and its necessary operands are 
almost identical to
 median value (which is less affected by outliers than the mean) is added
 back to the stack.
 
address@hidden connected-components
address@hidden Connected components
+Find the connected components in the input dataset (second popped
+operand). The first popped is the connectivity used in the connected
+components algorithm. The second popped operand is the dataset where
+connected components are to be found. It is assumed to be a binary image
+(with values of 0 or 1). It must have an 8-bit unsigned integer type which
+is the format produced by conditional operators. This operator will return
+a labeled dataset where the non-zero pixels in the input will be labeled
+with a counter (starting from 1).
+
+The connectivity is a number between 1 and the number of dimensions in the
+dataset (inclusive). 1 corresponds to the weakest (symmetric) connectivity
+between elements and the number of dimensions the strongest. For example on
+a 2D image, a connectivity of 1 corresponds to 4-connected neighbors and 2
+corresponds to 8-connected neighbors.
+
+One example usage of this operator can be the identification of regions
+above a certain threshold, as in the command below. With this command,
+Arithmetic will first separate all pixels greater than 100 into a binary
+image (where pixels with a value of 1 are above that value). Afterwards, it
+will label all those that are connected.
+
address@hidden
+$ astarithmetic in.fits 100 gt 2 connected-components
address@hidden example
+
+If your input dataset doesn't have a binary type, but you know all its
+values are 0 or 1, you can use the @code{uint8} operator (below) to convert
+it to binary.
+
 @item lt
 Less than: If the second popped (or left operand in infix notation, see
 @ref{Reverse polish notation}) value is smaller than the first popped