[gnuastro-commits] master 48a1059: Function to label holes in a binary image

[Mohammad Akhlaghi]

branch: master
commit 48a1059d66f6d1742019f9b59c3bb1229772a654
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>

    Function to label holes in a binary image
    
    This commit started with a testing a possible new strategy to fill holes in
    NoiseChisel (during the growth of detections, using their S/N instead of
    area). It was thus necessary to label the holes and
    `gal_binary_holes_label' was written.
    
    During the test, I also noticed that the `gal_data_ptr_dist' function name
    was confusing (the "pointer-distance" count be interpretted as the
    byte-base difference). The function is now named `gal_data_num_between' to
    be more clear that it returns the number of elements between the pointers
    in the type of the dataset.
    
    Ultimate, the test didn't give good results: the S/N in the holes of the
    diffuse regions have too much scatter to be meaningful. So the minor
    function name change and new function will remain in the library, but
    NoiseChisel still works as it did before.
---
 NEWS                        |  5 +++
 bin/mkcatalog/upperlimit.c  |  4 +--
 bin/noisechisel/detection.c |  4 +--
 bin/segment/clumps.c        | 14 ++++----
 doc/gnuastro.texi           | 22 ++++++++++---
 lib/binary.c                | 78 ++++++++++++++++++++++++++++++++++++++++++++-
 lib/data.c                  |  6 ++--
 lib/gnuastro/binary.h       |  6 +++-
 lib/gnuastro/data.h         |  2 +-
 lib/tile.c                  | 14 ++++----
 lib/wcs.c                   |  2 +-
 11 files changed, 128 insertions(+), 29 deletions(-)

diff --git a/NEWS b/NEWS
index 96ac855..129095b 100644
--- a/NEWS
+++ b/NEWS
@@ -68,6 +68,7 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
     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_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.
     gal_eps_name_is_eps: Returns 1 if given filename is EPS.
     gal_eps_suffix_is_eps: Returns 1 if given suffix is EPS.
@@ -146,6 +147,10 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
     gal_statistics_number: the output dataset now has a `size_t' type. Until
         now it was `uint64_t'.
 
+  Library:
+    - gal_binary_holes_fill: new name for `gal_binary_fill_holes'.
+    - gal_data_num_between: new name for `gal_data_ptr_dist'.
+
 ** Bug fixes
 
   bug #52979: Many unused result warnings for asprintf in some compilers.
diff --git a/bin/mkcatalog/upperlimit.c b/bin/mkcatalog/upperlimit.c
index c408467..ac34591 100644
--- a/bin/mkcatalog/upperlimit.c
+++ b/bin/mkcatalog/upperlimit.c
@@ -160,8 +160,8 @@ upperlimit_random_range(struct mkcatalog_passparams *pp, 
gal_data_t *tile,
   /* Set the minimum and maximum acceptable value for the range.  */
   if(p->uprange)
     {
-      tstart=gal_data_ptr_dist(tile->block->array, tile->array,
-                               p->objects->type);
+      tstart=gal_data_num_between(tile->block->array, tile->array,
+                                  p->objects->type);
       gal_dimension_index_to_coord(tstart, ndim, dsize, coord);
     }
 
diff --git a/bin/noisechisel/detection.c b/bin/noisechisel/detection.c
index cf1f7d0..fe56734 100644
--- a/bin/noisechisel/detection.c
+++ b/bin/noisechisel/detection.c
@@ -249,7 +249,7 @@ detection_fill_holes_open(void *in_prm)
 
       /* Fill the holes in this tile: holes with maximal connectivity means
          that they are most strongly bounded. */
-      gal_binary_fill_holes(copy, copy->ndim, -1);
+      gal_binary_holes_fill(copy, copy->ndim, -1);
       if(fho_prm->step==1)
         {
           detection_write_in_large(tile, copy);
@@ -1010,7 +1010,7 @@ detection_quantile_expand(struct noisechiselparams *p, 
gal_data_t *workbin)
       bf=(b=workbin->array)+workbin->size;
       do *b = (*o++ == 1); while(++b<bf);
       workbin=gal_binary_dilate(workbin, 1, 1, 1);
-      gal_binary_fill_holes(workbin, 1, p->detgrowmaxholesize);
+      gal_binary_holes_fill(workbin, 1, p->detgrowmaxholesize);
 
       /* Get the labeled image. */
       numexpanded=gal_binary_connected_components(workbin, &p->olabel,
diff --git a/bin/segment/clumps.c b/bin/segment/clumps.c
index 90d1e36..ee21fb3 100644
--- a/bin/segment/clumps.c
+++ b/bin/segment/clumps.c
@@ -621,9 +621,9 @@ clumps_find_make_sn_table(void *in_prm)
              rivers and not include them in the list of indexs to set
              clumps. To do that, we need this tile's starting
              coordinates. */
-          gal_dimension_index_to_coord(gal_data_ptr_dist(p->clabel->array,
-                                                         tile->array,
-                                                         p->clabel->type),
+          gal_dimension_index_to_coord(gal_data_num_between(p->clabel->array,
+                                                            tile->array,
+                                                            p->clabel->type),
                                        ndim, dsize, scoord);
 
 
@@ -661,16 +661,16 @@ clumps_find_make_sn_table(void *in_prm)
                     {
                       if(cltprm.id==282) *i+=2;
                   */
-                      indarr[c++]=gal_data_ptr_dist(p->clabel->array, i,
-                                                    p->clabel->type);
+                      indarr[c++]=gal_data_num_between(p->clabel->array, i,
+                                                       p->clabel->type);
                   /*
                     }
                   else
                     if(cltprm.id==282)
                       {
                         int32_t *clabel=p->clabel->array;
-                        size_t kjd=gal_data_ptr_dist(p->clabel->array, i,
-                                                     p->clabel->type);
+                        size_t kjd=gal_data_num_between(p->clabel->array, i,
+                                                        p->clabel->type);
                         printf("%zu, %zu: %u\n", kjd%dsize[1]+1,
                                kjd/dsize[1]+1, clabel[kjd]);
                       }
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index f8f30a0..2d51268 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -21108,10 +21108,10 @@ function will use the given @code{type} to calculate 
where the incremented
 element is located in memory.
 @end deftypefun
 
address@hidden size_t gal_data_ptr_dist (void @code{*earlier}, void 
@code{*later}, uint8_t @code{type})
-Return the number of elements between @code{earlier} and @code{later}
-assuming each element has a type defined by @code{type} (for the type
-codes, see @ref{Library data types}).
address@hidden size_t gal_data_num_between (void @code{*earlier}, void 
@code{*later}, uint8_t @code{type})
+Return the number of elements (in the given @code{type}) between
address@hidden and @code{later}. For the type codes, see @ref{Library data
+types}).
 @end deftypefun
 
 @deftypefun {void *} gal_data_malloc_array (uint8_t @code{type}, size_t 
@code{size}, const char @code{*funcname}, const char @code{*varname})
@@ -25234,7 +25234,19 @@ Although the adjacency matrix as used here is 
symmetric, currently this
 function assumes that it is filled on both sides of the diagonal.
 @end deftypefun
 
address@hidden void gal_binary_fill_holes (gal_data_t @code{*input}, int 
@code{connectivity})
address@hidden {gal_data_t *} gal_binary_holes_label (gal_data_t @code{*input}, 
int @code{connectivity}, size_t @code{*numholes})
+Label all the holes in the foreground (non-zero elements in input) as
+independent regions. Holes are background regions (zero-valued in input)
+that are fully surrounded by the foreground, as defined by
address@hidden The returned dataset has a 32-bit signed integer type
+with the size of the input. All holes in the input will have
+labels/counters greater or equal to @code{1}. The rest of the background
+regions will still have a value of @code{0} and the initial foreground
+pixels will have a value of @code{-1}. The total number of holes will be
+written where @code{numholes} points to.
address@hidden deftypefun
+
address@hidden void gal_binary_holes_fill (gal_data_t @code{*input}, int 
@code{connectivity})
 Fill all the holes (0 valued pixels surrounded by 1 valued pixels) within a
 region of the binary @code{input} dataset. The connectivity of the holes
 can be set with @code{connectivity}. This function currently only works on
diff --git a/lib/binary.c b/lib/binary.c
index be916cf..14cf026 100644
--- a/lib/binary.c
+++ b/lib/binary.c
@@ -618,6 +618,82 @@ binary_make_padded_inverse(gal_data_t *input, gal_data_t 
**outtile)
 
 
 
+gal_data_t *
+gal_binary_holes_label(gal_data_t *input, int connectivity,
+                       size_t *numholes)
+{
+  size_t d;
+  int32_t *lab;
+  gal_data_t *inv, *tile, *holelabs=NULL;
+
+  /* A small sanity check. */
+  if( input->type != GAL_TYPE_UINT8 )
+    error(EXIT_FAILURE, 0, "%s: input must have `uint8' type, but its "
+          "input dataset has `%s' type", __func__,
+          gal_type_name(input->type, 1));
+
+
+  /* Make the inverse image. */
+  inv=binary_make_padded_inverse(input, &tile);
+
+
+  /* Label the holes. Recall that the first label is just the undetected
+     regions, so we should subtract that from the total number.*/
+  *numholes=gal_binary_connected_components(inv, &holelabs, connectivity);
+  *numholes -= 1;
+
+
+  /* Any pixel with a label larger than 1 is a hole in the input image and
+     we should invert the respective pixel. To do it, we'll use the tile
+     that was defined before, just change its block and array.*/
+  tile->array=gal_tile_block_relative_to_other(tile, holelabs);
+  tile->block=holelabs; /* has to be after correcting `tile->array'. */
+
+
+  /* The type of the tile is already known (it is `int32_t') and we have no
+     output/other, so we'll just put `int' as a place-holder. In this way
+     we can avoid the switch statement of GAL_TILE_PARSE_OPERATE, and
+     directly use the workhorse macro `GAL_TILE_PO_OISET'. */
+  lab=(holelabs)->array;
+  GAL_TILE_PO_OISET(int32_t, int, tile, NULL, 0, 0, {
+      *lab++ = ( *i
+                 ? ( *i==1
+                     ? 0        /* Originally, was background. */
+                     : *i-1 )   /* Real label: -1 (background has 1). */
+                 : -1 );        /* Originally, was foreground. */
+    });
+
+
+  /* Clean up */
+  tile->array=NULL;
+  gal_data_free(inv);
+  gal_data_free(tile);
+
+
+  /* Correct the sizes of the hole labels array. We have already filled the
+     from the start, effectively removing the paddings. Therefore, ee will
+     just correct the sizes and we won't bother actually re-allocating the
+     array size in memory. According to the GNU C library's description
+     after `realloc': "In several allocation implementations, making a
+     block smaller sometimes necessitates copying it, so it can fail if no
+     other space is available.". The extra padding is only 2 pixels wide,
+     thus, the extra space is negligible compared to the actual array. So
+     it isn't worth possibly having to copy the whole array to another
+     location. Later, when we free the space, the kernel knows how much the
+     allocated size is. */
+  for(d=0;d<input->ndim;++d)
+    holelabs->dsize[d] = input->dsize[d];
+  holelabs->size=input->size;
+
+
+  /* Return the number of holes.  */
+  return holelabs;
+}
+
+
+
+
+
 /* Fill all the holes in an input unsigned char array.
 
    The basic method is this:
@@ -641,7 +717,7 @@ binary_make_padded_inverse(gal_data_t *input, gal_data_t 
**outtile)
       Any pixel with a label larger than 1, is therefore a bounded
       hole that is not 8-connected to the rest of the holes.  */
 void
-gal_binary_fill_holes(gal_data_t *input, int connectivity, size_t maxsize)
+gal_binary_holes_fill(gal_data_t *input, int connectivity, size_t maxsize)
 {
   uint8_t *in;
   uint32_t *i, *fi;
diff --git a/lib/data.c b/lib/data.c
index 9317b62..58f03d8 100644
--- a/lib/data.c
+++ b/lib/data.c
@@ -109,10 +109,10 @@ gal_data_ptr_increment(void *pointer, size_t increment, 
uint8_t type)
 
 
 
-/* Find the distance between two void pointers with a given type. See the
-   explanations before `gal_data_ptr_increment'. */
+/* Find the number of values between two void pointers with a given
+   type. See the explanations before `gal_data_ptr_increment'. */
 size_t
-gal_data_ptr_dist(void *earlier, void *later, uint8_t type)
+gal_data_num_between(void *earlier, void *later, uint8_t type)
 {
   char *e=(char *)earlier, *l=(char *)later;
   return (l-e)/gal_type_sizeof(type);
diff --git a/lib/gnuastro/binary.h b/lib/gnuastro/binary.h
index 4ae9357..5766123 100644
--- a/lib/gnuastro/binary.h
+++ b/lib/gnuastro/binary.h
@@ -97,8 +97,12 @@ gal_binary_connected_adjacency_matrix(gal_data_t *adjacency,
 /*********************************************************************/
 /*****************            Fill holes          ********************/
 /*********************************************************************/
+gal_data_t *
+gal_binary_holes_label(gal_data_t *input, int connectivity,
+                       size_t *numholes);
+
 void
-gal_binary_fill_holes(gal_data_t *input, int connectivity, size_t maxsize);
+gal_binary_holes_fill(gal_data_t *input, int connectivity, size_t maxsize);
 
 
 
diff --git a/lib/gnuastro/data.h b/lib/gnuastro/data.h
index c6332de..2e0aeee 100644
--- a/lib/gnuastro/data.h
+++ b/lib/gnuastro/data.h
@@ -239,7 +239,7 @@ void *
 gal_data_ptr_increment(void *pointer, size_t increment, uint8_t type);
 
 size_t
-gal_data_ptr_dist(void *earlier, void *later, uint8_t type);
+gal_data_num_between(void *earlier, void *later, uint8_t type);
 
 void *
 gal_data_malloc_array(uint8_t type, size_t size, const char *funcname,
diff --git a/lib/tile.c b/lib/tile.c
index b46c002..71c036b 100644
--- a/lib/tile.c
+++ b/lib/tile.c
@@ -65,7 +65,7 @@ gal_tile_start_coord(gal_data_t *tile, size_t *start_coord)
   else
     {
       /* Calculate the coordinates of the first pixel of the tile. */
-      ind = gal_data_ptr_dist(block->array, tile->array, block->type);
+      ind = gal_data_num_between(block->array, tile->array, block->type);
       gal_dimension_index_to_coord(ind, ndim, block->dsize, start_coord);
     }
 }
@@ -97,7 +97,7 @@ gal_tile_start_end_coord(gal_data_t *tile, size_t *start_end, 
int rel_block)
 
   /* Get the starting index. Note that for the type we need the allocated
      block dataset and can't rely on the tiles. */
-  start_ind=gal_data_ptr_dist(block->array, tile->array, block->type);
+  start_ind=gal_data_num_between(block->array, tile->array, block->type);
 
   /* Get the coordinates of the starting point relative to the allocated
      block. */
@@ -108,7 +108,7 @@ gal_tile_start_end_coord(gal_data_t *tile, size_t 
*start_end, int rel_block)
   if(host!=block)
     {
       /* Get the host's starting coordinates. */
-      start_ind=gal_data_ptr_dist(block->array, host->array, block->type);
+      start_ind=gal_data_num_between(block->array, host->array, block->type);
 
       /* Temporarily put the host's coordinates in the place held for the
          ending coordinates. */
@@ -145,7 +145,8 @@ gal_tile_start_end_ind_inclusive(gal_data_t *tile, 
gal_data_t *work,
   /* The starting index can be found from the distance of the `tile->array'
      pointer and `block->array' pointer. IMPORTANT: with the type of the
      block array.  */
-  start_end_inc[0]=gal_data_ptr_dist(block->array, tile->array, block->type);
+  start_end_inc[0]=gal_data_num_between(block->array, tile->array,
+                                        block->type);
 
 
   /* To find the end index, we need to know the coordinates of the starting
@@ -515,8 +516,9 @@ gal_tile_block_relative_to_other(gal_data_t *tile, 
gal_data_t *other)
 {
   gal_data_t *block=gal_tile_block(tile);
   return gal_data_ptr_increment(other->array,
-                                gal_data_ptr_dist(block->array,
-                                                  tile->array, block->type),
+                                gal_data_num_between(block->array,
+                                                     tile->array,
+                                                     block->type),
                                 other->type);
 }
 
diff --git a/lib/wcs.c b/lib/wcs.c
index a41393f..14d01bb 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -280,7 +280,7 @@ gal_wcs_on_tile(gal_data_t *tile)
       tile->wcs=gal_wcs_copy(block->wcs);
 
       /* Find the coordinates of the tile's starting index. */
-      start_ind=gal_data_ptr_dist(block->array, tile->array, block->type);
+      start_ind=gal_data_num_between(block->array, tile->array, block->type);
       gal_dimension_index_to_coord(start_ind, ndim, block->dsize, coord);
 
       /* Correct the copied WCS structure. Note that crpix is indexed in