[gnuastro-commits] master de2b395 1/2: NoiseChisel's labeling with large gthresh corrected

[Mohammad Akhlaghi]

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

    NoiseChisel's labeling with large gthresh corrected
    
    With a large `--gthresh' the previous commit would ignore the object and
    clump relabeling, but that was not the way to go, because we still need the
    separate clump labels to be set to 1. So in this commit, when no diffuse
    flux is given, we will manually set `clumptoobj' and `numobjects' and let
    same old algorithms continue.
    
    In the process, some other issues were fixed:
    
     - In MakeCatalog, we were mistakenly setting the `minmapsize' option to
       hidden (so it won't be read and a set to a value of 0). This has been
       corrected.
    
     - NoiseChisel's `--help' documentation for `tableformat' was corrected to
       be shorter to easily display in one line.
    
     - Some typos in the documentation for the discussion on the Adjacency
       matrix were corrected.
    
     - `gal_data_initialize' was mistakenly printing the value of `dsize[i]'
       with `%ld' instead of `%zu'. This remains from the time that `dsize' was
       in long type.
---
 bin/mkcatalog/ui.c             |   1 -
 bin/noisechisel/segmentation.c | 229 ++++++++++++++++++++++-------------------
 bin/noisechisel/ui.c           |   2 +-
 doc/gnuastro.texi              |  21 ++--
 lib/data.c                     |  12 ++-
 5 files changed, 146 insertions(+), 119 deletions(-)

diff --git a/bin/mkcatalog/ui.c b/bin/mkcatalog/ui.c
index c9f982e..71d49da 100644
--- a/bin/mkcatalog/ui.c
+++ b/bin/mkcatalog/ui.c
@@ -140,7 +140,6 @@ ui_initialize_options(struct mkcatalogparams *p,
         case GAL_OPTIONS_KEY_TYPE:
         case GAL_OPTIONS_KEY_SEARCHIN:
         case GAL_OPTIONS_KEY_IGNORECASE:
-        case GAL_OPTIONS_KEY_MINMAPSIZE:
           cp->coptions[i].flags=OPTION_HIDDEN;
           cp->coptions[i].mandatory=GAL_OPTIONS_NOT_MANDATORY;
           break;
diff --git a/bin/noisechisel/segmentation.c b/bin/noisechisel/segmentation.c
index 855f3a3..6c25980 100644
--- a/bin/noisechisel/segmentation.c
+++ b/bin/noisechisel/segmentation.c
@@ -100,14 +100,6 @@ segmentation_relab_to_objects(struct clumps_thread_params 
*cltprm)
   int32_t *ngblabs=gal_data_malloc_array(GAL_TYPE_UINT32, nngb, __func__,
                                          "ngblabs");
 
-  /* If there aren't any diffuse pixels (when a large `--gthresh' was
-     given), then it is just necessary to set two constants. */
-  if(cltprm->diffuseindexs->size==0)
-    {
-      cltprm->clumptoobj=NULL;
-      cltprm->numobjects=cltprm->numtrueclumps;
-      return;
-    }
 
   /* Go over all the still-unlabeled pixels (if they exist) and see which
      labels they touch. In the process, get the average value of the
@@ -115,117 +107,145 @@ segmentation_relab_to_objects(struct 
clumps_thread_params *cltprm)
      matrix. Note that at this point, the rivers are also part of the
      "diffuse" regions. So we don't need to go over all the indexs of this
      object, only its diffuse indexs. */
-  sf=(s=cltprm->diffuseindexs->array)+cltprm->diffuseindexs->size;
-  do
-    /* We only want to work on pixels that have already been identified as
-       touching more than one label: river pixels. */
-    if( olabel[ *s ]==CLUMPS_RIVER )
-      {
-        /* Initialize the values. */
-        i=ii=0;
-        rpnum=1;              /* River-pixel number of points used. */
-        rpsum=imgss[*s];      /* River-pixel sum of values used.    */
-        memset(ngblabs, 0, nngb*sizeof *ngblabs);
-
-        /* Check all the fully-connected neighbors of this pixel and see if
-           it touches a label or not */
-        GAL_DIMENSION_NEIGHBOR_OP(*s, ndim, dsize, ndim, dinc, {
-            if( olabel[nind] > 0 )
+  if(cltprm->diffuseindexs->size)
+    {
+      sf=(s=cltprm->diffuseindexs->array)+cltprm->diffuseindexs->size;
+      do
+        /* We only want to work on pixels that have already been identified
+           as touching more than one label: river pixels. */
+        if( olabel[ *s ]==CLUMPS_RIVER )
+          {
+            /* Initialize the values. */
+            i=ii=0;
+            rpnum=1;              /* River-pixel number of points used. */
+            rpsum=imgss[*s];      /* River-pixel sum of values used.    */
+            memset(ngblabs, 0, nngb*sizeof *ngblabs);
+
+            /* Check all the fully-connected neighbors of this pixel and
+               see if it touches a label or not */
+            GAL_DIMENSION_NEIGHBOR_OP(*s, ndim, dsize, ndim, dinc, {
+                if( olabel[nind] > 0 )
+                  {
+                    /* Add this neighbor's value and increment the number. */
+                    if( !isnan(imgss[nind]) ) { ++rpnum; rpsum+=imgss[nind]; }
+
+                    /* Go over the already found neighbors and see if this
+                       grown clump has already been considered or not. */
+                    for(i=0;i<ii;++i) if(ngblabs[i]==olabel[nind]) break;
+
+                    /* This is the first time we are getting to this
+                       neighbor: */
+                    if(i==ii) ngblabs[ ii++ ] = olabel[nind];
+                  }
+              } );
+
+            /* For a check:
+            if(ii>0)
               {
-                /* Add this neighbor's value and increment the number. */
-                if( !isnan(imgss[nind]) ) { ++rpnum; rpsum+=imgss[nind]; }
-
-                /* Go over the already found neighbors and see if this grown
-                   clump has already been considered or not. */
-                for(i=0;i<ii;++i) if(ngblabs[i]==olabel[nind]) break;
-
-                /* This is the first time we are getting to this neighbor: */
-                if(i==ii) ngblabs[ ii++ ] = olabel[nind];
+                printf("%zu, %zu:\n", *s%dsize[1]+1, *s/dsize[1]+1);
+                for(i=0;i<ii;++i) printf("\t%u\n", ngblabs[i]);
               }
-          } );
-
-        /* For a check:
-        if(ii>0)
-          {
-            printf("%zu, %zu:\n", *s%dsize[1]+1, *s/dsize[1]+1);
-            for(i=0;i<ii;++i) printf("\t%u\n", ngblabs[i]);
+            */
+
+            /* If more than one neighboring label was found, fill in the
+               'sums' and 'nums' adjacency matrixs with the values for this
+               pixel. Recall that ii is the number of neighboring labels to
+               this river pixel. */
+            if(ii>i)
+              for(i=0;i<ii;++i)
+                for(j=0;j<ii;++j)
+                  if(i!=j)
+                    {
+                      /* For safety, we will fill both sides of the
+                         diagonal. */
+                      ++nums[ ngblabs[i] * amwidth + ngblabs[j] ];
+                      ++nums[ ngblabs[j] * amwidth + ngblabs[i] ];
+                      sums[ ngblabs[i] * amwidth + ngblabs[j] ] +=
+                        rpsum/rpnum;
+                      sums[ ngblabs[j] * amwidth + ngblabs[i] ] +=
+                        rpsum/rpnum;
+                    }
           }
-        */
-
-        /* If more than one neighboring label was found, fill in the 'sums'
-           and 'nums' adjacency matrixs with the values for this
-           pixel. Recall that ii is the number of neighboring labels to
-           this river pixel. */
-        if(ii>i)
-          for(i=0;i<ii;++i)
-            for(j=0;j<ii;++j)
-              if(i!=j)
-                {
-                  /* For safety, we will fill both sides of the diagonal. */
-                  ++nums[ ngblabs[i] * amwidth + ngblabs[j] ];
-                  ++nums[ ngblabs[j] * amwidth + ngblabs[i] ];
-                  sums[ ngblabs[i] * amwidth + ngblabs[j] ] += rpsum/rpnum;
-                  sums[ ngblabs[j] * amwidth + ngblabs[i] ] += rpsum/rpnum;
-                }
-      }
-  while(++s<sf);
+      while(++s<sf);
 
 
-  /* We now have the average values and number of all rivers between the
-     grown clumps. We now want to finalize their connection (given the
-     user's criteria). */
-  for(i=1;i<amwidth;++i)
-    for(j=1;j<i;++j)
-      {
-        ii = i * amwidth + j;
-        if(nums[ii]>p->minriverlength)       /* There is a connection. */
+      /* We now have the average values and number of all rivers between
+         the grown clumps. We now want to finalize their connection (given
+         the user's criteria). */
+      for(i=1;i<amwidth;++i)
+        for(j=1;j<i;++j)
           {
-            /* For easy reading. */
-            ave=sums[ii]/nums[ii];
-
-            /* In case the average is negative (only possible if `sums' is
-               negative), don't change the adjacency: it is already
-               initialized to zero. Note that even an area of 1 is
-               acceptable, and we put no area criteria here, because the
-               fact that a river exists between two clumps is important. */
-            if( ave>0.0f && ( c * ave / sqrt(ave+err) ) > p->objbordersn )
+            ii = i * amwidth + j;
+            if(nums[ii]>p->minriverlength)       /* There is a connection. */
               {
-                adjacency[ii]=1;       /* We want to set both sides of the */
-                adjacency[ j * amwidth + i ] = 1; /* matrix diagonal to 1. */
+                /* For easy reading. */
+                ave=sums[ii]/nums[ii];
+
+                /* In case the average is negative (only possible if `sums'
+                   is negative), don't change the adjacency: it is already
+                   initialized to zero. Note that even an area of 1 is
+                   acceptable, and we put no area criteria here, because
+                   the fact that a river exists between two clumps is
+                   important. */
+                if( ave>0.0f && ( c * ave / sqrt(ave+err) ) > p->objbordersn )
+                  {
+                    adjacency[ii]=1;   /* We want to set both sides of the */
+                    adjacency[ j * amwidth + i ] = 1; /* Symmetric matrix. */
+                  }
               }
           }
-      }
 
 
-  /* For a check:
-  if(cltprm->id==XXX)
-    {
-      printf("=====================\n");
-      printf("%zu:\n--------\n", cltprm->id);
-      for(i=1;i<amwidth;++i)
+      /* For a check:
+      if(cltprm->id==XXX)
         {
-          printf(" %zu...\n", i);
-          for(j=1;j<amwidth;++j)
+          printf("=====================\n");
+          printf("%zu:\n--------\n", cltprm->id);
+          for(i=1;i<amwidth;++i)
             {
-              ii=i*amwidth+j;
-              if(nums[ii])
+              printf(" %zu...\n", i);
+              for(j=1;j<amwidth;++j)
                 {
-                  ave=sums[ii]/nums[ii];
-                  printf("    ...%zu: N:%-4zu S:%-10.2f S/N: %-10.2f "
-                         "--> %u\n", j, nums[ii], sums[ii],
-                         c*ave/sqrt(ave+err), adjacency[ii]);
+                  ii=i*amwidth+j;
+                  if(nums[ii])
+                    {
+                      ave=sums[ii]/nums[ii];
+                      printf("    ...%zu: N:%-4zu S:%-10.2f S/N: %-10.2f "
+                             "--> %u\n", j, nums[ii], sums[ii],
+                             c*ave/sqrt(ave+err), adjacency[ii]);
+                    }
                 }
+              printf("\n");
             }
-          printf("\n");
         }
+      */
+
+
+      /* Calculate the new labels for each grown clump. */
+      cltprm->clumptoobj = gal_binary_connected_adjacency_matrix(adjacency_d,
+                                                         &cltprm->numobjects);
+      clumptoobj = cltprm->clumptoobj->array;
     }
-  */
 
+  /* There was no list of diffuse pixels, this happens when the user sets a
+     very high `gthresh' threshold and wants to make sure that each clump
+     is a separate object. So we need to define the number of objects and
+     `clumptoobj' manually. */
+  else
+    {
+      /* Allocate the `clumptoobj' array. */
+      cltprm->clumptoobj = gal_data_alloc(NULL, GAL_TYPE_INT32, 1, &amwidth,
+                                          NULL, 1, p->cp.minmapsize, NULL,
+                                          NULL, NULL);
+      clumptoobj = cltprm->clumptoobj->array;
+
+      /* Fill in the `clumptoobj' array with the indexs of the objects. */
+      for(i=0;i<amwidth;++i) clumptoobj[i]=i;
+
+      /* Set the number of objects. */
+      cltprm->numobjects = cltprm->numtrueclumps;
+    }
 
-  /* Calculate the new labels for each grown clump. */
-  cltprm->clumptoobj = gal_binary_connected_adjacency_matrix(adjacency_d,
-                                                      &cltprm->numobjects);
-  clumptoobj = cltprm->clumptoobj->array;
 
   /* For a check
   if(cltprm->id==XXXX)
@@ -238,6 +258,7 @@ segmentation_relab_to_objects(struct clumps_thread_params 
*cltprm)
     }
   */
 
+
   /* Correct all the labels. */
   sf=(s=cltprm->indexs->array)+cltprm->indexs->size;
   do
@@ -435,11 +456,13 @@ segmentation_on_threads(void *in_prm)
 
           /* If the user has asked for grown clumps in the clumps image
              instead of the raw clumps, then replace the indexs in the
-             `clabel' array is well. */
+             `clabel' array is well. In this case, there will always be one
+             "clump". */
           if(p->grownclumps)
             {
               sf=(s=cltprm.indexs->array)+cltprm.indexs->size;
               do clabel[ *s++ ] = 1; while(s<sf);
+              cltprm.numtrueclumps=1;
             }
         }
       else
@@ -493,10 +516,8 @@ segmentation_on_threads(void *in_prm)
 
           /* Correct the clump labels. Note that this is only necessary
              when there is more than object over the detection or when
-             there were multiple clumps, but no ID conversion was necessary
-             (very high `--gthresh' values). In the latter case,
-             `clumptoobj' will be NULL.*/
-          if(cltprm.numobjects>1 && cltprm.clumptoobj)
+             there were multiple clumps over the detection. */
+          if(cltprm.numobjects>1)
             segmentation_relab_clumps_in_objects(&cltprm);
           gal_data_free(cltprm.clumptoobj);
           if(clprm->step==6) {continue;}
diff --git a/bin/noisechisel/ui.c b/bin/noisechisel/ui.c
index 8601cd0..3165f33 100644
--- a/bin/noisechisel/ui.c
+++ b/bin/noisechisel/ui.c
@@ -142,7 +142,7 @@ ui_initialize_options(struct noisechiselparams *p,
 
         case GAL_OPTIONS_KEY_TABLEFORMAT:
           cp->coptions[i].mandatory=GAL_OPTIONS_MANDATORY;
-          cp->coptions[i].doc="Formats: `txt', `fits-ascii', `fits-binary'.";
+          cp->coptions[i].doc="`txt', `fits-ascii', `fits-binary'.";
           break;
         }
     }
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 8a754e1..60f3d92 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -21261,9 +21261,9 @@ array[3]: Standard deviation.
 Binary datasets only have two (usable) values: 0 (also known as background)
 or 1 (also known as foreground). They are created after some binary
 classificataion is applied to the dataset. The most common is thresholding:
-in an image, pixels with a value above the threshold is commonly assigned a
-value of 1 and those that have a value less than the threshold are assigned
-a value of 0.
+for example in an image, pixels with a value above the threshold are given
+a value of 1 and those with a value less than the threshold are assigned a
+value of 0.
 
 @cindex Connectivity
 @cindex Immediate neighbors
@@ -21392,13 +21392,14 @@ be labeled). Blank pixels in the input will also be 
blank in the output.
 @deftypefun {gal_data_t *} gal_binary_connected_adjacency_matrix (gal_data_t 
@code{*adjacency}, size_t @code{*numconnected})
 @cindex Adjacency matrix
 @cindex Matrix, adjacency
-Find the number of connected labels and new labels based on an an adjacency
-matrix (which must be binary array of type @code{GAL_TYPE_UINT8}). The
-returned dataset is a list of new labels for each old label. In other
-words, this function will find the objects that are connected (possibly
-through a third object) and in the output array, the respective elements
-for all input labels is going to have the same value. The total number of
-connected labels is put into the space that @code{numconnected} points to.
+Find the number of connected labels and new labels based on an adjacency
+matrix, which must be a square binary array (type
address@hidden). The returned dataset is a list of new labels for
+each old label. In other words, this function will find the objects that
+are connected (possibly through a third object) and in the output array,
+the respective elements for all input labels is going to have the same
+value. The total number of connected labels is put into the space that
address@hidden points to.
 
 An adjacency matrix defines connection between two labels. For example,
 let's assume we have 5 labels and we know that labels 1 and 5 are connected
diff --git a/lib/data.c b/lib/data.c
index 3115a1a..97ee961 100644
--- a/lib/data.c
+++ b/lib/data.c
@@ -186,13 +186,16 @@ gal_data_mmap(gal_data_t *data, int clear)
   char *filename;
   size_t bsize=data->size*gal_type_sizeof(data->type);
 
+
   /* Check if the .gnuastro folder exists, write the file there. If it
      doesn't exist, then make the .gnuastro directory.*/
   gal_checkset_mkdir(".gnuastro");
 
+
   /* Set the filename */
   gal_checkset_allocate_copy("./.gnuastro/mmap_XXXXXX", &filename);
 
+
   /* Create a zero-sized file and keep its descriptor.  */
   errno=0;
   /*filedes=open(filename, O_RDWR | O_CREAT | O_EXCL | O_TRUNC );*/
@@ -222,14 +225,17 @@ gal_data_mmap(gal_data_t *data, int clear)
   data->array=mmap(NULL, bsize, PROT_READ | PROT_WRITE, MAP_SHARED,
                    filedes, 0);
 
+
   /* Close the file. */
   if( close(filedes) == -1 )
     error(EXIT_FAILURE, errno, "%s: %s couldn't be closed",
           __func__, filename);
 
+
   /* Keep the filename. */
   data->mmapname=filename;
 
+
   /* If it was supposed to be cleared, then clear the memory. */
   if(clear) memset(data->array, 0, bsize);
 }
@@ -303,9 +309,9 @@ gal_data_initialize(gal_data_t *data, void *array, uint8_t 
type,
         {
           /* Do a small sanity check. */
           if(dsize[i]<=0)
-            error(EXIT_FAILURE, 0, "%s: the size of a dimension cannot be zero 
"
-                  "or negative. dsize[%zu], but has a value of %ld", __func__,
-                  i, dsize[i]);
+            error(EXIT_FAILURE, 0, "%s: the size of a dimension cannot be "
+                  "zero or negative. dsize[%zu], but has a value of %zu",
+                  __func__, i, dsize[i]);
 
           /* Write this dimension's size, also correct the total number of
              elements. */