[gnuastro-commits] master 1cb69ef: NoiseChisel segmentation deals with no usable pixels

[Mohammad Akhlaghi]

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

    NoiseChisel segmentation deals with no usable pixels
    
    It may happen that the number of usable pixels to define clumps over
    un-detected pixels in some tiles are zero. However, the sky clump
    identification and S/N calculations weren't ready for that. In such a
    cases, at the start of segmentation, NoiseChisel would crash with a
    "segmentation fault (core dumped)".
    
    The particular case that this happend was when a large area of the
    surroundings of the image were blank. So many tiles were blank. However,
    the particular image and tile sizes used resulted in some tiles only having
    one pixel row of overlap with the non-blank image. Since all pixels on the
    edge of a tile are not considered in finding clumps over the sky region,
    the result was that there were actually no usable pixels for sky
    segmentation within the tile. Such situations won't cause a segmentation
    fault any more.
    
    This bug was reported by Raúl Infante Sainz and Nacho Trujillo.
    
    In parallel some corrections to the book were also made regarding the
    recent updates to MakeProfiles of the last commit.
    
    This fixes bug #51372.
---
 THANKS                   |  1 +
 bin/noisechisel/clumps.c | 38 ++++++++++++++++++++++++++++++--------
 doc/gnuastro.texi        | 18 ++++++++++--------
 3 files changed, 41 insertions(+), 16 deletions(-)

diff --git a/THANKS b/THANKS
index 63651a2..bed7f05 100644
--- a/THANKS
+++ b/THANKS
@@ -22,6 +22,7 @@ support in Gnuastro. The list is ordered alphabetically.
     Rosa Calvi                           address@hidden
     Antonio Diaz Diaz                    address@hidden
     Takashi Ichikawa                     address@hidden
+    Raúl Infante Sainz                   address@hidden
     Brandon Invergo                      address@hidden
     Lee Kelvin                           address@hidden
     Mohammad-Reza Khellat                address@hidden
diff --git a/bin/noisechisel/clumps.c b/bin/noisechisel/clumps.c
index c16b9d0..459b212 100644
--- a/bin/noisechisel/clumps.c
+++ b/bin/noisechisel/clumps.c
@@ -107,6 +107,10 @@ clumps_oversegment(struct clumps_thread_params *cltprm)
   **********************************************/
 
 
+  /* If the size of the indexs is zero, then this function is pointless. */
+  if(indexs->size==0) { cltprm->numinitclumps=0; return; }
+
+
   /* Sort the given indexs based on their flux (`gal_qsort_index_arr' is
      defined as static in `gnuastro/qsort.h') */
   gal_qsort_index_arr=p->conv->array;
@@ -146,10 +150,10 @@ clumps_oversegment(struct clumps_thread_params *cltprm)
 
             /* A small sanity check. */
             if(Q!=NULL || cleanup!=NULL)
-              error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s so we 
"
-                    "can fix this problem. `Q' and `cleanup' should be NULL "
-                    "but while checking the equal flux regions they aren't",
-                    __func__, PACKAGE_BUGREPORT);
+              error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s so "
+                    "we can fix this problem. `Q' and `cleanup' should be "
+                    "NULL but while checking the equal flux regions they "
+                    "aren't", __func__, PACKAGE_BUGREPORT);
 
             /* Add this pixel to a queue. */
             gal_list_sizet_add(&Q, *a);
@@ -808,6 +812,10 @@ clumps_make_sn_table(struct clumps_thread_params *cltprm)
   size_t i, ind, counter=0, infodsize[2]={tablen, INFO_NCOLS};
 
 
+  /* If there were no initial clumps, then ignore this function. */
+  if(cltprm->numinitclumps==0) { cltprm->sn=NULL; return; }
+
+
   /* Allocate the arrays to keep the final S/N table (and possibly S/N
      index) for this object or tile. */
   cltprm->sn        = &cltprm->clprm->sn[ cltprm->id ];
@@ -920,8 +928,8 @@ clumps_correct_sky_labels_for_check(struct 
clumps_thread_params *cltprm,
 
   /* A small sanity check. */
   if(gal_tile_block(tile)!=p->clabel)
-    error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to address the "
-          "problem. `tile->block' must point to the `clabel' dataset",
+    error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to address "
+          "the problem. `tile->block' must point to the `clabel' dataset",
           __func__, PACKAGE_BUGREPORT);
 
 
@@ -1001,6 +1009,7 @@ clumps_find_make_sn_table(void *in_prm)
       cltprm.id = tind  = tprm->indexs[i];
       tile = &p->ltl.tiles[tind];
 
+
       /* Change the tile's pointers to the binary image (which has 1 for
          detected pixels and 0 for un-detected regions). */
       tarray=tile->array;
@@ -1008,6 +1017,7 @@ clumps_find_make_sn_table(void *in_prm)
       tile->array = gal_tile_block_relative_to_other(tile, p->binary);
       tile->block = p->binary;
 
+
       /* Get the number of usable elements in this tile (note that tiles
          can have blank pixels), so we can't simply use `tile->size'. */
       if(tile->flag & GAL_DATA_FLAG_HASBLANK)
@@ -1018,12 +1028,14 @@ clumps_find_make_sn_table(void *in_prm)
         }
       else num=tile->size;
 
+
       /* Find the number of detected pixels over this tile. Since this is
          the binary image, this is just the sum of all the pixels. */
       tmp=gal_statistics_sum(tile);
       numdet=*((double *)(tmp->array));
       gal_data_free(tmp);
 
+
       /* See if this tile should be used or not (has enough undetected
          pixels). Note that it might happen that some tiles are fully
          blank. In such cases, it is important to first check the number of
@@ -1043,6 +1055,7 @@ clumps_find_make_sn_table(void *in_prm)
           tile->array = gal_tile_block_relative_to_other(tile, p->clabel);
           tile->block = p->clabel;
 
+
           /* We need to set all the pixels on the edge of the tile to
              rivers and not include them in the list of indexs to set
              clumps. To do that, we need this tile's starting
@@ -1052,6 +1065,7 @@ clumps_find_make_sn_table(void *in_prm)
                                                          p->clabel->type),
                                        ndim, dsize, scoord);
 
+
           /* Add the index of every sky element to the array of
              indexs. Note that since we know the array is always of type
              `int32_t', we can call the `GAL_TILE_PO_OISET' macro to avoid
@@ -1103,30 +1117,37 @@ clumps_find_make_sn_table(void *in_prm)
                 }
             });
 
+
           /* Correct the number of indexs. */
           cltprm.indexs->size=cltprm.indexs->dsize[0]=c;
 
+
           /* Generate the clumps over this region. */
           clumps_oversegment(&cltprm);
 
+
           /* Set all river pixels to CLUMPS_INIT (to be distinguishable
              from the detected regions). */
           GAL_TILE_PO_OISET( int32_t, int, tile, NULL, 0, 1,
                              {if(*i==CLUMPS_RIVER) *i=CLUMPS_INIT;} );
 
+
           /* For a check, the step variable will be set. */
           if(clprm->step==1)
             { gal_data_free(cltprm.indexs); continue; }
 
+
           /* Make the clump S/N table. */
           clumps_make_sn_table(&cltprm);
 
+
           /* If the user wanted to check the steps, remove the clumps that
              weren't used from the `clabel' image (they have been already
              excluded from the table). */
           if(cltprm.snind)
             clumps_correct_sky_labels_for_check(&cltprm, tile);
 
+
           /* Clean up. */
           gal_data_free(cltprm.indexs);
         }
@@ -1255,10 +1276,11 @@ clumps_true_find_sn_thresh(struct noisechiselparams *p)
   else
     {
       clprm.step=0;
-      gal_threads_spin_off(clumps_find_make_sn_table, &clprm, p->ltl.tottiles,
-                           p->cp.numthreads);
+      gal_threads_spin_off(clumps_find_make_sn_table, &clprm,
+                           p->ltl.tottiles, p->cp.numthreads);
     }
 
+
   /* Destroy the mutex if it was initialized. */
   if( p->cp.numthreads>1 && (p->checksegmentation || p->checkclumpsn) )
     pthread_mutex_destroy(&clprm.labmutex);
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index b901403..46f4c6c 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -1938,8 +1938,8 @@ $ astmkprof -P
 [[[ ... Truncated lines ... ]]]
 
 # Columns, by info (see `--searchin'), or number (starting from 1):
- xcol         2           # Center along first FITS axis (horizontal).
- ycol         3           # Center along second FITS axis (vertical).
+ ccol         2           # Center along first FITS axis (horizontal).
+ ccol         3           # Center along second FITS axis (vertical).
  fcol         4           # sersic (1), moffat (2), gaussian (3),
                           # point (4), flat (5), circumference (6).
  rcol         5           # Effective radius or FWHM in pixels.
@@ -2015,8 +2015,7 @@ necessary parameters and runs MakeProfiles with the 
following command:
 
 @example
 
-$ astmkprof --prepforconv --naxis1=500 --naxis2=500         \
-            --zeropoint=18.0 cat.txt
+$ astmkprof --prepforconv --naxis=500,500 --zeropoint=18.0 cat.txt
 MakeProfiles started on Sat Oct  6 16:26:56 953
   - 6 profiles read from cat.txt
   - Random number generator (RNG) type: mt19937
@@ -2182,8 +2181,7 @@ base=cat
 rm out.fits
 
 # Run MakeProfiles to create an oversampled FITS image.
-astmkprof --prepforconv --naxis1=500 --naxis2=500                \
-          --zeropoint=18.0 "$base".txt
+astmkprof --prepforconv --naxis=500,500 --zeropoint=18.0 "$base".txt
 
 # Convolve the created image with the kernel.
 astconvolve --kernel=0_"$base".fits "$base".fits
@@ -14718,11 +14716,11 @@ $ astmkprof --background=image.fits catalog.txt
 $ astmkprof --kernel=moffat,2.8,5 --oversample=1
 
 ## Make profiles in catalog, using RA and Dec in the given column:
-$ astmkprof --racol=RA_CENTER --ycol=DEC_CENTER catalog.txt
+$ astmkprof --ccol=RA_CENTER --ccol=DEC_CENTER --mode=wcs catalog.txt
 
 ## Make a 1500x1500 merged image (oversampled 500x500) image along
 ## with an individual image for all the profiles in catalog:
-$ astmkprof --individual --oversample 3 -x500 -y500 catalog.txt
+$ astmkprof --individual --oversample 3 --naxis=500,500 catalog.txt
 @end example
 
 @noindent
@@ -15210,6 +15208,10 @@ dimensions, for example @code{CRPIX1} and 
@code{PC2_1}. You can see the
 FITS headers with Gnuastro's @ref{Fits} program using a command like this:
 @command{$ astfits -p image.fits}.
 
+If the values given to all of these options does not correspond to the
+dimensionality of the output dataset, then no WCS information will be
+added.
+
 @table @option
 
 @item --crpix=FLT,FLT