[gnuastro-commits] master b9f899c5: Arithmetic: new number-neighbors operator

[Mohammad Akhlaghi]

branch: master
commit b9f899c5b83820750d0546bf466d7ad332f038ee
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    Arithmetic: new number-neighbors operator
    
    Until now, there was no easy way to identify the number of non-zero
    neighbors of each input non-zero pixel.
    
    With this commit, Gnuastro now has the 'number-neighbors' option for doing
    this. To do this, a new 'gal_binary_number_neighbors' was added to the
    'binary.h' library of Gnuastro.
---
 NEWS                        |  3 +++
 bin/arithmetic/arithmetic.c | 38 ++++++++++++++++++++++++---
 bin/arithmetic/arithmetic.h |  1 +
 doc/gnuastro.texi           | 18 +++++++++++++
 lib/binary.c                | 63 +++++++++++++++++++++++++++++++++++++++++++++
 lib/gnuastro/binary.h       |  8 ++++++
 6 files changed, 127 insertions(+), 4 deletions(-)

diff --git a/NEWS b/NEWS
index 401f7bad..5dc63391 100644
--- a/NEWS
+++ b/NEWS
@@ -34,6 +34,8 @@ See the end of the file for license conditions.
        vertices of a rectangle from its center and width along RA and
        Dec. This takes into account the curved nature of the coordinate
        system. Added after discussion with Martin Kuemmel.
+     - number-neighbors: Return the number of non-zero neighbors of each
+       non-zero pixel in a binary image.
    - New operators (only in the Arithmetic program):
      - interpolate-meanngb: interpolate blank values with mean of the
        requested number of nearest neighbors.
@@ -144,6 +146,7 @@ See the end of the file for license conditions.
    - GAL_ARITHMETIC_OP_NANOMAGGY_TO_COUNTS: convert nanomaggy to counts.
    - GAL_ARITHMETIC_OP_BOX_VERTICES_ON_SPHERE: calculate the coordinates of
      vertices of a rectable on a sphere from its center and width/height.
+   - gal_binary_number_neighbors: num. non-zero neighbors of non-zero pixels.
    - gal_data_alloc_empty: Allocate an empty dataset with a given number of
      dimensions.
    - gal_list_f64_to_data: convert list of float64s to a 'gal_data_t'
diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index 97471a8c..76ab80e1 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -578,10 +578,13 @@ arithmetic_binary_sanity_checks(gal_data_t *in, 
gal_data_t *conn,
 
   /* Make sure the array has an unsigned 8-bit type. */
   if(in->type!=GAL_TYPE_UINT8)
-    error(EXIT_FAILURE, 0, "the second popped operand of '%s' doesn't "
-          "have an 8-bit unsigned integer type. It must be a binary "
-          "dataset (only being equal to zero is checked). You can use "
-          "the 'uint8' operator for type conversion", operator);
+    error(EXIT_FAILURE, 0, "the second popped operand of '%s' has a type "
+          "of %s. However, it must be a binary dataset (only being equal "
+          "to zero is checked). You can use the 'uint8' operator for type "
+          "conversion, alternatively, if all your values are positive "
+          "and floating point, you can use '0 gt', if you want non-blank "
+          "values you can use 'isblank not' and many other operators that "
+          "produce a binary output", operator, gal_type_name(in->type, 1));
 
   /* Clean up and return the integer value of 'conn'. */
   gal_data_free(conn);
@@ -623,6 +626,27 @@ arithmetic_erode_dilate(struct arithmeticparams *p, char 
*token, int op)
 
 
 
+static void
+arithmetic_number_neighbors(struct arithmeticparams *p, char *token, int op)
+{
+  int conn_int;
+
+  /* Pop the two necessary operands. */
+  gal_data_t *conn = operands_pop(p, token);
+  gal_data_t *in   = operands_pop(p, token);
+
+  /* Do the sanity checks and do the job. */
+  conn_int=arithmetic_binary_sanity_checks(in, conn, token);
+  in=gal_binary_number_neighbors(in, conn_int, 1);
+
+  /* Push the result onto the stack. */
+  operands_add(p, NULL, in);
+}
+
+
+
+
+
 static void
 arithmetic_connected_components(struct arithmeticparams *p, char *token)
 {
@@ -1321,6 +1345,8 @@ arithmetic_set_operator(char *string, size_t 
*num_operands, int *inlib)
         { op=ARITHMETIC_OP_ERODE;                 *num_operands=0; }
       else if (!strcmp(string, "dilate"))
         { op=ARITHMETIC_OP_DILATE;                *num_operands=0; }
+      else if (!strcmp(string, "number-neighbors"))
+        { op=ARITHMETIC_OP_NUMBER_NEIGHBORS;      *num_operands=0; }
       else if (!strcmp(string, "connected-components"))
         { op=ARITHMETIC_OP_CONNECTED_COMPONENTS;  *num_operands=0; }
       else if (!strcmp(string, "fill-holes"))
@@ -1514,6 +1540,10 @@ arithmetic_operator_run(struct arithmeticparams *p, int 
operator,
           arithmetic_erode_dilate(p, operator_string, operator);
           break;
 
+        case ARITHMETIC_OP_NUMBER_NEIGHBORS:
+          arithmetic_number_neighbors(p, operator_string, operator);
+          break;
+
         case ARITHMETIC_OP_CONNECTED_COMPONENTS:
           arithmetic_connected_components(p, operator_string);
           break;
diff --git a/bin/arithmetic/arithmetic.h b/bin/arithmetic/arithmetic.h
index b82254ba..7abdcfbf 100644
--- a/bin/arithmetic/arithmetic.h
+++ b/bin/arithmetic/arithmetic.h
@@ -37,6 +37,7 @@ enum arithmetic_prog_operators
   ARITHMETIC_OP_FILTER_SIGCLIP_MEDIAN,
   ARITHMETIC_OP_ERODE,
   ARITHMETIC_OP_DILATE,
+  ARITHMETIC_OP_NUMBER_NEIGHBORS,
   ARITHMETIC_OP_CONNECTED_COMPONENTS,
   ARITHMETIC_OP_FILL_HOLES,
   ARITHMETIC_OP_INVERT,
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 157cfd43..5b44afba 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -18122,6 +18122,19 @@ The usage is similar to @code{erode}, for example:
 $ astarithmetic binary.fits 2 dilate -oout.fits
 @end example
 
+@item number-neighbors
+Return a dataset of the same size as the second popped operand, but where each 
non-zero and non-blank input pixel is replaced with the number of its non-zero 
and non-blank neighbors.
+The first popped operand is the connectivity (see above) and must be a 
single-value of an integer type.
+The dataset is assumed to be binary (having an unsigned, 8-bit dataset).
+
+For example with the command below, you can select all pixels above a value of 
100 in your image with the ``greater-than'' or @code{gt} operator (see 
@ref{Conditional operators}).
+Recall that the output of all conditional operators is a binary output (having 
a value of 0 or 1).
+In the same command, we will then find how many neighboring pixels of each 
pixel (that was originally above the threshold) are also above the threshold.
+
+@example
+$ astarithmetic image.fits 100 gt 2 number-neighbors
+@end example
+
 @item connected-components
 @cindex Connected components
 Find the connected components in the input dataset (second popped operand).
@@ -37882,6 +37895,11 @@ foreground. In this implementation, @code{num} 
erosions are going to be
 applied on the dataset, then @code{num} dilations.
 @end deftypefun
 
+@deftypefun {gal_data_t *} gal_binary_number_neighbors (gal_data_t 
@code{*input}, int @code{connectivity}, int @code{inplace})
+Return an image of the same size as the input, but where each non-zero and 
non-blank input pixel is replaced with the number of its non-zero and non-blank 
neighbors.
+The input dataset is assumed to be binary (having an unsigned, 8-bit dataset).
+The neighbors are defined through the @code{connectivity} argument (see above) 
and if @code{inplace!=0}, then the output will be written into the input.
+@end deftypefun
 
 @deftypefun size_t gal_binary_connected_components (gal_data_t @code{*binary}, 
gal_data_t @code{**out}, int @code{connectivity})
 @cindex Breadth first search
diff --git a/lib/binary.c b/lib/binary.c
index 292b0438..780de4d4 100644
--- a/lib/binary.c
+++ b/lib/binary.c
@@ -419,6 +419,69 @@ gal_binary_open(gal_data_t *input, size_t num, int 
connectivity,
 
 
 
+/*********************************************************************/
+/*****************            Neighbors           ********************/
+/*********************************************************************/
+/* This is a general erosion and dilation function. It is less efficient
+   than the more specialized cases above. */
+gal_data_t *
+gal_binary_number_neighbors(gal_data_t *input, int connectivity, int inplace)
+{
+  gal_data_t *out;
+  uint8_t n, *narr, *byt=input->array;
+  size_t i, *dinc=gal_dimension_increment(input->ndim, input->dsize);
+
+  /* Currently this only works on blocks. */
+  if(input->block)
+    error(EXIT_FAILURE, 0, "%s: currently only works on a fully "
+          "allocated block of memory, but the input is a tile (its 'block' "
+          "element is not NULL)", __func__);
+
+  /* The input must have a uint8 datatype. */
+  if(input->type!=GAL_TYPE_UINT8)
+    error(EXIT_FAILURE, 0, "%s: input must have an unsigned 8-bit "
+          "datatype but has a type of %s\n", __func__,
+          gal_type_name(input->type, 1));
+
+  /* Allocate the output dataset. */
+  out = ( inplace
+          ? input
+          : gal_data_alloc(NULL, GAL_TYPE_UINT8, input->ndim, input->dsize,
+                           input->wcs, 1, input->minmapsize, input->quietmmap,
+                           NULL, NULL, NULL) );
+  narr=out->array;
+
+  /* Go over the neighbors of each pixel. */
+  for(i=0;i<input->size;++i)
+    if(byt[i] && byt[i]!=GAL_BLANK_UINT8)
+      {
+        n=0;
+        GAL_DIMENSION_NEIGHBOR_OP(i, input->ndim, input->dsize, connectivity,
+                                  dinc, { n += byt[nind]>0; });
+        narr[i]=n;
+      }
+
+  /* Return the output dataset. */
+  return out;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
 /*********************************************************************/
 /*****************      Connected components      ********************/
 /*********************************************************************/
diff --git a/lib/gnuastro/binary.h b/lib/gnuastro/binary.h
index 39e7b0c2..7758bff7 100644
--- a/lib/gnuastro/binary.h
+++ b/lib/gnuastro/binary.h
@@ -80,6 +80,12 @@ gal_binary_open(gal_data_t *input, size_t num, int 
connectivity,
 
 
 
+/*********************************************************************/
+/*****************            Neighbors           ********************/
+/*********************************************************************/
+gal_data_t *
+gal_binary_number_neighbors(gal_data_t *input, int connectivity, int inplace);
+
 
 
 /*********************************************************************/
@@ -101,6 +107,8 @@ gal_binary_connected_adjacency_list(gal_list_sizet_t 
**listarr,
                                     size_t number, size_t minmapsize,
                                     int quietmmap, size_t *numconnected);
 
+
+
 /*********************************************************************/
 /*****************            Fill holes          ********************/
 /*********************************************************************/