gnuastro-commits
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[gnuastro-commits] master 6f4e484 032/113: Merged with recent updates in


From: Mohammad Akhlaghi
Subject: [gnuastro-commits] master 6f4e484 032/113: Merged with recent updates in master
Date: Fri, 16 Apr 2021 10:33:37 -0400 (EDT)

branch: master
commit 6f4e4847d7be3b439897fc4407fd80341d9f39aa
Merge: 6086579 3479c92
Author: Mohammad Akhlaghi <akhlaghi@gnu.org>
Commit: Mohammad Akhlaghi <akhlaghi@gnu.org>

    Merged with recent updates in master
    
    Recent updates are imported into this branch with only a one-line conflict
    resolved.
---
 NEWS                                |  18 ++++
 THANKS                              |   1 +
 bin/fits/main.h                     |  31 +++----
 bin/fits/ui.c                       |  24 ++++-
 bin/noisechisel/args.h              |  13 +++
 bin/noisechisel/astnoisechisel.conf |   1 +
 bin/noisechisel/clumps.c            |   3 +-
 bin/noisechisel/detection.c         |  12 ++-
 bin/noisechisel/main.h              |   1 +
 bin/noisechisel/noisechisel.c       |   2 +-
 bin/noisechisel/segmentation.c      |  15 ++--
 bin/noisechisel/threshold.c         | 153 ++++++++++++++++++++++++++++++-
 bin/noisechisel/ui.h                |   1 +
 doc/announce-acknowledge.txt        |   2 +
 doc/gnuastro.texi                   | 174 +++++++++++++++++++++++++++++-------
 lib/data.c                          |  25 +++++-
 lib/gnuastro-internal/commonopts.h  |   2 +-
 lib/tile.c                          |   2 +-
 lib/wcs.c                           |  60 ++++++++++++-
 tmpfs-config-make                   |   7 +-
 20 files changed, 475 insertions(+), 72 deletions(-)

diff --git a/NEWS b/NEWS
index ec75d66..8ae540f 100644
--- a/NEWS
+++ b/NEWS
@@ -5,12 +5,23 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
 
 ** New features
 
+  Fits: when an extension/HDU is identified on the command-line with the
+  `--hdu' option and no operation is requested, the full list of header
+  keywords in that HDU will be printed (as if only `--printallkeys' was
+  called).
+
   NoiseChisel: with the new `--widekernel' option it is now possible to use
   a wider kernel to identify which tiles contain signal. The rest of the
   steps (identifying the quantile threshold on the selected tiles and etc)
   are done on the dataset convolved with `--kernel' as they were
   before. Since it is time consuming, this is an optional feature.
 
+  NoiseChisel: with the new `--qthreshtilequant' option, it is now possible
+  to discard high-valued (outlier) tiles before estimating qthresh over the
+  whole image. This can be useful in detecting very large diffuse/flat
+  regions that would otherwise be detected as background (and effectively
+  removed).
+
   NoiseChisel: the finally selected true detections are now grown based on
   signal contiguity, not by blind dilation. The growth process is the same
   as the growing of clumps to define objects. Only for true detections, the
@@ -24,6 +35,11 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
 
 ** Removed features
 
+  NoiseChisel: The `--dilate' and `--dilatengb' options have been
+  removed. Growing of true detections is no longer done through dilation
+  but through the `--detgrowquant' and `--detgrowmaxholesize' options (see
+  above).
+
 ** Changed features
 
   `gal_binary_fill_holes' now accepts a `connectivity' and `maxsize'
@@ -41,6 +57,8 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
   Arithmetic not accounting for integer blank pixels in binary operators
   (bug #52014).
 
+  NoiseChisel segfault when memory mapping to a file (bug #52043).
+
 
 
 
diff --git a/THANKS b/THANKS
index 5f34189..1afb512 100644
--- a/THANKS
+++ b/THANKS
@@ -18,6 +18,7 @@ support in Gnuastro. The list is ordered alphabetically (by 
family name).
     Marjan Akbari                        mrjakbari@gmail.com
     Karl Berry                           karl@gnu.org
     Roland Bacon                         roland.bacon@univ-lyon1.fr
+    Leindert Boogaard                    boogaard@strw.leidenuniv.nl
     Nicolas Bouché                       nicolas.bouche@irap.omp.eu
     Fernando Buitrago                    fbuitrago@oal.ul.pt
     Adrian Bunk                          bunk@debian.org
diff --git a/bin/fits/main.h b/bin/fits/main.h
index 46bae7d..782850d 100644
--- a/bin/fits/main.h
+++ b/bin/fits/main.h
@@ -53,21 +53,22 @@ enum fits_mode
 struct fitsparams
 {
   /* From the environment. */
-  struct gal_options_common_params cp;  /* Common parameters.     */
-  char          *filename;     /* Name of input file.             */
-  gal_list_str_t  *remove;     /* Remove extensions from a file.  */
-  gal_list_str_t    *copy;     /* Copy extensions to output.      */
-  gal_list_str_t     *cut;     /* Copy ext. to output and remove. */
-  uint8_t    printallkeys;     /* Print all the header keywords.  */
-  uint8_t            date;     /* Set DATE to current time.       */
-  char           *comment;     /* COMMENT value.                  */
-  char           *history;     /* HISTORY value.                  */
-  gal_list_str_t    *asis;     /* Strings to write asis.          */
-  gal_list_str_t  *delete;     /* Keywords to remove.             */
-  gal_list_str_t  *rename;     /* Rename a keyword.               */
-  gal_list_str_t  *update;     /* For keywords to update.         */
-  gal_list_str_t   *write;     /* Full arg. for keywords to add.  */
-  uint8_t     quitonerror;     /* Quit if an error occurs.        */
+  struct gal_options_common_params cp;  /* Common parameters.       */
+  int  hdu_in_commandline;     /* HDU wasn't given in config. file. */
+  char          *filename;     /* Name of input file.               */
+  gal_list_str_t  *remove;     /* Remove extensions from a file.    */
+  gal_list_str_t    *copy;     /* Copy extensions to output.        */
+  gal_list_str_t     *cut;     /* Copy ext. to output and remove.   */
+  uint8_t    printallkeys;     /* Print all the header keywords.    */
+  uint8_t            date;     /* Set DATE to current time.         */
+  char           *comment;     /* COMMENT value.                    */
+  char           *history;     /* HISTORY value.                    */
+  gal_list_str_t    *asis;     /* Strings to write asis.            */
+  gal_list_str_t  *delete;     /* Keywords to remove.               */
+  gal_list_str_t  *rename;     /* Rename a keyword.                 */
+  gal_list_str_t  *update;     /* For keywords to update.           */
+  gal_list_str_t   *write;     /* Full arg. for keywords to add.    */
+  uint8_t     quitonerror;     /* Quit if an error occurs.          */
 
   /* Internal: */
   int                         mode;  /* Operating on HDUs or keywords.  */
diff --git a/bin/fits/ui.c b/bin/fits/ui.c
index 6f1388e..33576d2 100644
--- a/bin/fits/ui.c
+++ b/bin/fits/ui.c
@@ -266,7 +266,15 @@ ui_read_check_only_options(struct fitsparams *p)
   /* If no options are given, go into HDU mode, which will print the HDU
      information when nothing is asked. */
   if(p->mode==FITS_MODE_INVALID)
-    p->mode=FITS_MODE_HDU;
+    {
+      if(p->hdu_in_commandline)
+        {
+          p->printallkeys=1;
+          p->mode = FITS_MODE_KEY;
+        }
+      else
+        p->mode = FITS_MODE_HDU;
+    }
 }
 
 
@@ -408,8 +416,8 @@ ui_fill_fits_headerll(gal_list_str_t *input, 
gal_fits_list_key_t **output)
               errno=0;
               fvalue=dp=malloc(sizeof *dp);
               if(dp==NULL)
-                error(EXIT_FAILURE, errno, "%s: allocating %zu bytes for `dp'",
-                      __func__, sizeof *dp);
+                error(EXIT_FAILURE, errno, "%s: allocating %zu bytes for "
+                      "`dp'", __func__, sizeof *dp);
               *dp=d;
             }
           else
@@ -461,6 +469,7 @@ ui_preparations(struct fitsparams *p)
 void
 ui_read_check_inputs_setup(int argc, char *argv[], struct fitsparams *p)
 {
+  size_t i;
   struct gal_options_common_params *cp=&p->cp;
 
 
@@ -485,6 +494,15 @@ ui_read_check_inputs_setup(int argc, char *argv[], struct 
fitsparams *p)
     error(EXIT_FAILURE, errno, "parsing arguments");
 
 
+  /* Check if the HDU is specified on the command-line. If so, then later,
+     if no operation is requested, we will print the header of the given
+     HDU.*/
+  for(i=0; !gal_options_is_last(&cp->coptions[i]); ++i)
+    if(cp->coptions[i].key==GAL_OPTIONS_KEY_HDU
+       && cp->coptions[i].set)
+      p->hdu_in_commandline=1;
+
+
   /* Read the configuration files and set the common values. */
   gal_options_read_config_set(&p->cp);
 
diff --git a/bin/noisechisel/args.h b/bin/noisechisel/args.h
index 80de2ce..cf9934b 100644
--- a/bin/noisechisel/args.h
+++ b/bin/noisechisel/args.h
@@ -221,6 +221,19 @@ struct argp_option program_options[] =
       GAL_OPTIONS_NOT_SET
     },
     {
+      "qthreshtilequant",
+      UI_KEY_QTHRESHTILEQUANT,
+      "FLT",
+      0,
+      "Remove tiles at higher quantiles.",
+      ARGS_GROUP_DETECTION,
+      &p->qthreshtilequant,
+      GAL_TYPE_FLOAT32,
+      GAL_OPTIONS_RANGE_GE_0_LE_1,
+      GAL_OPTIONS_MANDATORY,
+      GAL_OPTIONS_NOT_SET
+    },
+    {
       "smoothwidth",
       UI_KEY_SMOOTHWIDTH,
       "INT",
diff --git a/bin/noisechisel/astnoisechisel.conf 
b/bin/noisechisel/astnoisechisel.conf
index a037355..b3ab828 100644
--- a/bin/noisechisel/astnoisechisel.conf
+++ b/bin/noisechisel/astnoisechisel.conf
@@ -30,6 +30,7 @@
  mirrordist         1.5
  modmedqdiff       0.01
  qthresh            0.3
+ qthreshtilequant   1.0
  smoothwidth          3
  erode                2
  erodengb             4
diff --git a/bin/noisechisel/clumps.c b/bin/noisechisel/clumps.c
index 481d44f..8965226 100644
--- a/bin/noisechisel/clumps.c
+++ b/bin/noisechisel/clumps.c
@@ -1324,6 +1324,7 @@ clumps_true_find_sn_thresh(struct noisechiselparams *p)
                            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);
@@ -1441,7 +1442,7 @@ clumps_det_label_indexs(struct noisechiselparams *p)
      to allocate). */
   areas=gal_data_calloc_array(GAL_TYPE_SIZE_T, p->numdetections+1, __func__,
                               "areas");
-  if(p->input->flag & GAL_DATA_FLAG_HASBLANK)
+  if(gal_blank_present(p->input, 1))
     {
       lf=(l=p->olabel->array)+p->olabel->size; /* Blank pixels have a      */
       do if(*l>0) ++areas[*l]; while(++l<lf);  /* negative value in int32. */
diff --git a/bin/noisechisel/detection.c b/bin/noisechisel/detection.c
index b9858f0..405612c 100644
--- a/bin/noisechisel/detection.c
+++ b/bin/noisechisel/detection.c
@@ -387,7 +387,17 @@ detection_pseudo_find(struct noisechiselparams *p, 
gal_data_t *workbin,
 
       /* Clean up: the array in `bin' should just be replaced with that in
          `workbin' because it is used in later steps. */
-      free(workbin->array);
+      if(workbin->mmapname)
+        {
+          /* Delete the memory mapped file and set the filename of `bin'
+             for `workbin'. */
+          remove(workbin->mmapname);
+          free(workbin->mmapname);
+          workbin->mmapname=bin->mmapname;
+          bin->mmapname=NULL;
+        }
+      else
+        free(workbin->array);
       workbin->array=bin->array;
       bin->name=bin->array=NULL;
       gal_data_free(bin);
diff --git a/bin/noisechisel/main.h b/bin/noisechisel/main.h
index 5d08e59..c994c4f 100644
--- a/bin/noisechisel/main.h
+++ b/bin/noisechisel/main.h
@@ -61,6 +61,7 @@ struct noisechiselparams
   float            mirrordist;  /* Maximum distance to check mode sym.    */
   float           modmedqdiff;  /* Difference between mode and median.    */
   float               qthresh;  /* Quantile threshold on convolved image. */
+  float      qthreshtilequant;  /* Remove tiles with lower quantile.      */
   size_t          smoothwidth;  /* Width of flat kernel to smooth.        */
   uint8_t        checkqthresh;  /* Save the quantile threhsold steps.     */
   size_t                erode;  /* Number of erosions after thresholding. */
diff --git a/bin/noisechisel/noisechisel.c b/bin/noisechisel/noisechisel.c
index 58aae42..6ccf12f 100644
--- a/bin/noisechisel/noisechisel.c
+++ b/bin/noisechisel/noisechisel.c
@@ -71,7 +71,7 @@ noisechisel_convolve(struct noisechiselparams *p)
       if(p->widekernel)
         gal_timing_report(&t1, "Convolved with sharper kernel.", 1);
       else
-        gal_timing_report(&t1, "Convolved with kernel.", 1);
+        gal_timing_report(&t1, "Convolved with given kernel.", 1);
     }
   if(p->detectionname)
     {
diff --git a/bin/noisechisel/segmentation.c b/bin/noisechisel/segmentation.c
index 0b6bd65..53be9a7 100644
--- a/bin/noisechisel/segmentation.c
+++ b/bin/noisechisel/segmentation.c
@@ -411,6 +411,7 @@ segmentation_on_threads(void *in_prm)
       /* Find the clumps over this region. */
       clumps_oversegment(&cltprm);
 
+
       /* Make the clump S/N table. This table is made before (possibly)
          stopping the process (if a check is requested). This is because if
          the user has also asked for a check image, we can break out of the
@@ -424,6 +425,7 @@ segmentation_on_threads(void *in_prm)
         clumps_make_sn_table(&cltprm);
       else cltprm.sn=&clprm->sn[ cltprm.id ];
 
+
       /* If the user wanted to check the segmentation steps or the clump
          S/N values in a table, then we have to stop the process at this
          point. */
@@ -435,6 +437,7 @@ segmentation_on_threads(void *in_prm)
       gal_data_free(topinds);
       if(clprm->step==2) continue;
 
+
       /* Set the internal (with the detection) clump and object
          labels. Segmenting a detection into multiple objects is only
          defined when there is more than one true clump over the
@@ -448,12 +451,10 @@ segmentation_on_threads(void *in_prm)
           cltprm.numobjects=1;
           segmentation_relab_noseg(&cltprm);
 
-
           /* If the user wanted a check image, this object doesn't
              change. */
           if( clprm->step >= 3 && clprm->step <= 6) continue;
 
-
           /* 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. In this case, there will always be one
@@ -474,7 +475,6 @@ segmentation_on_threads(void *in_prm)
           if(clprm->step==3)
             { gal_data_free(cltprm.diffuseindexs); continue; }
 
-
           /* If grown clumps are desired instead of the raw clumps, then
              replace all the grown clumps with those in clabel. */
           if(p->grownclumps)
@@ -485,7 +485,6 @@ segmentation_on_threads(void *in_prm)
               while(++s<sf);
             }
 
-
           /* Identify the objects in this detection using the grown clumps
              and correct the grown clump labels into new object labels. */
           segmentation_relab_to_objects(&cltprm);
@@ -496,7 +495,6 @@ segmentation_on_threads(void *in_prm)
               continue;
             }
 
-
           /* Continue the growth and cover the whole area, we don't need
              the diffuse indexs any more, so after filling the detected
              region, free the indexs. */
@@ -514,7 +512,6 @@ segmentation_on_threads(void *in_prm)
           gal_data_free(cltprm.diffuseindexs);
           if(clprm->step==5) { gal_data_free(cltprm.clumptoobj); continue; }
 
-
           /* 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 over the detection. */
@@ -760,8 +757,9 @@ segmentation_detections(struct noisechiselparams *p)
 
             default:
               error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s so "
-                    "we can address the issue. The value %d is not recognized "
-                    "for clprm.step", __func__, PACKAGE_BUGREPORT, clprm.step);
+                    "we can address the issue. The value %d is not "
+                    "recognized for clprm.step", __func__, PACKAGE_BUGREPORT,
+                    clprm.step);
             }
 
           /* Write the demonstration array into the check image. The
@@ -791,6 +789,7 @@ segmentation_detections(struct noisechiselparams *p)
                            p->cp.numthreads);
     }
 
+
   /* Save the final number of objects and clumps. */
   p->numclumps=clprm.totclumps;
   p->numobjects=clprm.totobjects;
diff --git a/bin/noisechisel/threshold.c b/bin/noisechisel/threshold.c
index ad3cce0..a0da6bf 100644
--- a/bin/noisechisel/threshold.c
+++ b/bin/noisechisel/threshold.c
@@ -266,7 +266,6 @@ threshold_interp_smooth(struct noisechiselparams *p, 
gal_data_t **first,
     error(EXIT_FAILURE, 0, "%s: `third' must not have any `next' pointer.",
           __func__);
 
-
   /* Do the interpolation of both arrays. */
   (*first)->next = *second;
   if(third) (*second)->next = *third;
@@ -500,6 +499,149 @@ qthresh_on_tile(void *in_prm)
 
 
 
+/* The main working function for `threshold_qthresh_clean'. The main
+   purpose/problem is this: when we have channels, the qthresh values for
+   each channel should be treated independently. */
+static void
+threshold_qthresh_clean_work(struct noisechiselparams *p, gal_data_t *first,
+                             gal_data_t *second, gal_data_t *third,
+                             size_t start, size_t number)
+{
+  gal_data_t *quantile;
+  size_t i, osize=first->size;
+  float *oa1=NULL, *oa2=NULL, *oa3=NULL;
+  float q, *arr1=NULL, *arr2=NULL, *arr3=NULL;
+
+  /* A small sanity check. */
+  if(first->type!=GAL_TYPE_FLOAT32)
+    error(EXIT_FAILURE, 0, "%s: datatype has to be float32", __func__);
+
+  /* Correct the arrays (if necessary). IMPORTANT: The datasets are
+     multi-dimensional. However, when estimating the quantile, their
+     dimensionality doesn't matter (only the `size' element is checked by
+     `gal_statistics_quantile', not `ndim' or `dsize'). So we just need to
+     correct `size' if channels are to be considered. */
+  if(start || number!=first->size)
+    {
+      /* Keep the original values for re-setting later. */
+      oa1=first->array;
+      oa2=second->array;
+      if(third) oa3=third->array;
+
+      /* Increment the array pointers. */
+      first->array=gal_data_ptr_increment(first->array, start, first->type);
+      second->array=gal_data_ptr_increment(second->array, start,
+                                           second->type);
+      if(third)
+        third->array=gal_data_ptr_increment(third->array, start, third->type);
+
+      /* Correct their sizes. */
+      first->size=number;
+      second->size=number;
+      if(third) third->size=number;
+    }
+
+  /* Find the quantile and remove all tiles that are more than it in the
+     first array. */
+  arr1=first->array;
+  quantile=gal_statistics_quantile(first, p->qthreshtilequant, 0);
+  q=*((float *)(quantile->array));
+  for(i=0;i<first->size;++i)
+    /* Just note that we have blank (NaN) values, so to avoid doing a
+       NaN check with `isnan', we will check if the value is below the
+       quantile, if it succeeds (isn't NaN and is below the quantile),
+       then we'll put it's actual value, otherwise, a NaN. */
+    arr1[i] = arr1[i]<q ? arr1[i] : NAN;
+  gal_data_free(quantile);
+
+  /* Second quantile threshold. */
+  arr2=second->array;
+  quantile=gal_statistics_quantile(second, p->qthreshtilequant, 0);
+  q=*((float *)(quantile->array));
+  for(i=0;i<second->size;++i)
+    arr2[i] = arr2[i]<q ? arr2[i] : NAN;
+  gal_data_free(quantile);
+
+  /* The third (if it exists). */
+  if(third)
+    {
+      arr3=third->array;
+      quantile=gal_statistics_quantile(third, p->qthreshtilequant, 0);
+      q=*((float *)(quantile->array));
+      for(i=0;i<third->size;++i)
+        arr3[i] = arr3[i]<q ? arr3[i] : NAN;
+      gal_data_free(quantile);
+    }
+
+  /* Make sure all three have the same NaN pixels. */
+  for(i=0;i<first->size;++i)
+    if( isnan(arr1[i]) || isnan(arr2[i]) || isnan(arr3[i]) )
+      {
+        arr1[i] = arr2[i] = NAN;
+        if(third) arr3[i] = NAN;
+      }
+
+  /* Correct the values if they were changed. */
+  if(start || number!=osize)
+    {
+      first->array=oa1;
+      second->array=oa2;
+      first->size = second->size = osize;
+      if(third) { third->array=oa3; third->size=osize; }
+    }
+}
+
+
+
+
+
+/* Clean higher valued quantile thresholds: useful when the diffuse (almost
+   flat) structures are much larger than the tile size. */
+static void
+threshold_qthresh_clean(struct noisechiselparams *p, gal_data_t *first,
+                        gal_data_t *second, gal_data_t *third,
+                        char *filename)
+{
+  size_t i;
+  struct gal_tile_two_layer_params *tl=&p->cp.tl;
+
+  /* A small sanity check: */
+  if(first->size!=tl->tottiles)
+    error(EXIT_FAILURE, 0, "%s: `first->size' and `tl->tottiles' must have "
+          "the same value, but they don't: %zu, %zu", __func__, first->size,
+          tl->tottiles);
+
+  /* If the input is from a tile structure and the user has asked to ignore
+     channels, then re-order the values. */
+  for(i=0;i<tl->totchannels;++i)
+    threshold_qthresh_clean_work(p, first, second, third,
+                                 i*tl->tottilesinch, tl->tottilesinch);
+
+  /* If the user wants to see the steps. */
+  if(p->qthreshname)
+    {
+      first->name="QTHRESH_ERODE_CLEAN";
+      second->name="QTHRESH_NOERODE_CLEAN";
+      gal_tile_full_values_write(first, tl, 1, p->qthreshname, NULL,
+                                 PROGRAM_NAME);
+      gal_tile_full_values_write(second, tl, 1, p->qthreshname, NULL,
+                                 PROGRAM_NAME);
+      first->name=second->name=NULL;
+
+      if(third)
+        {
+          third->name="QTHRESH_EXPAND_CLEAN";
+          gal_tile_full_values_write(third, tl, 1, p->qthreshname,
+                                     NULL, PROGRAM_NAME);
+          third->name=NULL;
+        }
+    }
+}
+
+
+
+
+
 void
 threshold_quantile_find_apply(struct noisechiselparams *p)
 {
@@ -583,9 +725,16 @@ threshold_quantile_find_apply(struct noisechiselparams *p)
     }
 
 
+  /* Remove higher thresholds if requested. */
+  if(p->qthreshtilequant!=1.0)
+    threshold_qthresh_clean(p, qprm.erode_th, qprm.noerode_th,
+                            qprm.expand_th ? qprm.expand_th : NULL,
+                            p->qthreshname);
+
+
   /* Interpolate and smooth the derived values. */
   threshold_interp_smooth(p, &qprm.erode_th, &qprm.noerode_th,
-                          qprm.expand_th ? & qprm.expand_th : NULL,
+                          qprm.expand_th ? &qprm.expand_th : NULL,
                           p->qthreshname);
 
 
diff --git a/bin/noisechisel/ui.h b/bin/noisechisel/ui.h
index bfe3b22..c8209ae 100644
--- a/bin/noisechisel/ui.h
+++ b/bin/noisechisel/ui.h
@@ -67,6 +67,7 @@ enum option_keys_enum
   UI_KEY_ONLYDETECTION,
   UI_KEY_GROWNCLUMPS,
   UI_KEY_SMOOTHWIDTH,
+  UI_KEY_QTHRESHTILEQUANT,
   UI_KEY_CHECKQTHRESH,
   UI_KEY_ERODENGB,
   UI_KEY_NOERODEQUANT,
diff --git a/doc/announce-acknowledge.txt b/doc/announce-acknowledge.txt
index b76e5bb..4e2bfea 100644
--- a/doc/announce-acknowledge.txt
+++ b/doc/announce-acknowledge.txt
@@ -1,4 +1,6 @@
 This file is meant to keep the names of the people who's help must be
 acknowledged in the next release.
 
+Leindert Boogaard
 Takashi Ichikawa
+David Valls-Gabaud
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 474640a..55c6d5b 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -4596,6 +4596,37 @@ display in the @option{--help} output of the program.
 
 @table @option
 
+@item --minmapsize=INT
+The minimum size (in bytes) to store the contents of each main processing
+array of a program as a file (on the non-volatile HDD/SSD), not in
+RAM. This can be very useful for large datasets which can be very memory
+intensive such that your RAM will not be sufficient to keep/process them. A
+random filename is assigned to the array (in a @file{.gnuastro} directory
+within the running directory) which will keep the contents of the array as
+long as it is necessary.
+
+When this option has a value of @code{0} (zero), all arrays that use this
+option in a program will actually be in a file (not in RAM). When the value
+is @code{-1} (largest possible number in the unsigned integer types) these
+arrays will be definitely allocated in RAM. However, for complex programs
+like @ref{NoiseChisel}, it is recommended to not set it to @code{0}, but a
+value like @code{10000} so the many small arrays necessary during
+processing are stored in RAM and only larger ones are saved as a file.
+
+Please note that using a non-volatile file (in the HDD/SDD) instead of RAM
+can significantly increase the program's running time, especially on
+HDDs. So it is best to give this option large values by default. You can
+then decrease it for a specific program's invocation on a large input after
+you see memory issues arise (for example an error, or the program not
+aborting and fully consuming your memory).
+
+The random file will be deleted once it is no longer needed by the
+program. The @file{.gnuastro} directory will also be deleted if it has no
+other contents (you may also have configuration files in this directory,
+see @ref{Configuration files}). If you see randomly named files remaining
+in this directory, please send us a bug report so we address the problem,
+see @ref{Report a bug}.
+
 @item -Z INT[,INT[,...]]
 @itemx --tilesize=[,INT[,...]]
 The size of regular tiles for tessellation, see @ref{Tessellation}. For
@@ -6647,10 +6678,13 @@ One line examples:
 $ astfits image.fits
 
 ## Print the header keywords in the second HDU (counting from 0):
-$ astfits image.fits -h1 --printallkeys
+$ astfits image.fits -h1
 
 ## Only print header keywords that contain `NAXIS':
-$ astfits image.fits -h1 --printallkeys | grep NAXIS
+$ astfits image.fits -h1 | grep NAXIS
+
+## Only print the WCS standard PC matrix elements
+$ astfits image.fits -h1 | grep 'PC._.'
 
 ## Copy a HDU from input.fits to out.fits:
 $ astfits input.fits --copy=hdu-name --output=out.fits
@@ -6706,6 +6740,18 @@ HDU (extension) information: `image.fits'.
 3      CATALOG         table_binary    12371x5
 @end example
 
+If a specific HDU is identified on the command-line with the @option{--hdu}
+(or @option{-h} option) and no operation requested, then the full list of
+header keywords in that HDU will be printed (as if the
+@option{--printallkeys} was called, see below). It is important to remember
+that this only occurs when @option{--hdu} is given on the command-line. The
+@option{--hdu} value given in a configuration file will only be used when a
+specific operation on keywords requested. Therefore as described in the
+paragraphs above, when no explicit call to the @option{--hdu} option is
+made on the command-line and no operation is requested (on the command-line
+or configuration files), the basic information of each HDU/extension is
+printed.
+
 The operating mode and input/output options to Fits are similar to the
 other programs and fully described in @ref{Common options}. The options
 particular to Fits can be divided into two groups: 1) those related to
@@ -6736,8 +6782,8 @@ extensions will be copied from the input FITS file to the 
output FITS file
 in the given order (on the command-line and also in configuration files,
 see @ref{Configuration file precedence}). If the separate classes are
 called together in one run of Fits, then first @option{--copy} is run (on
-all specified HDUs), followed @option{--cut} (again on all specified HDUs),
-and then @option{--remove} (on all specified HDUs).
+all specified HDUs), followed by @option{--cut} (again on all specified
+HDUs), and then @option{--remove} (on all specified HDUs).
 
 The @option{--copy} and @option{--cut} options need an output FITS file
 (specified with the @option{--output} option). If the output file exists,
@@ -12177,15 +12223,43 @@ successfully.
 Before using NoiseChisel it is strongly recommended to read
 @url{https://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa [2015]} to
 gain a good understanding of what it does and how each parameter influences
-the output. Thanks to that paper, there is no more need to continue this
-introduction any further and we can just dive into the details of running
-NoiseChisel in @ref{Invoking astnoisechisel}. However, the paper cannot
-undergo any further updates, but NoiseChisel will evolve: better algorithms
-or steps will be found, thus options will be added or removed. So this book
-is the final and definitive guide. To make the transition form the paper to
-this book easier (and encourage reading the paper), below you can see the
-major changes since that paper was published.
+the output. Thanks to that paper, there is no need to go into the details
+of the major processing steps. Hence we can just dive into the details of
+running NoiseChisel in @ref{Invoking astnoisechisel}.
+
+However, the paper cannot undergo any further updates, but NoiseChisel will
+evolve: better algorithms or steps will be found, thus options will be
+added or removed. So this book is the final and definitive guide. For a
+more detailed list of changes in each release, please follow the
+@file{NEWS} file. The @file{NEWS} file is in the released Gnuastro tarball
+(see @ref{Release tarball}). You can also view the most recent @file{NEWS}
+file @url{http://git.savannah.gnu.org/cgit/gnuastro.git/plain/NEWS,
+online}@footnote{The online version of the @file{NEWS} file here may
+contain features that have been implemented, but not yet officially
+released as a tarball (see @ref{Downloading the source}). Therefore, please
+be sure to look at the dates and versions above each group of changed
+features and make sure it corresponds to your installed version. It is
+hence recommended to always stay up to date, see @ref{Announcements}).}.
+
+To make the transition form the paper to this book easier (and encourage
+reading the paper), below you can see the major changes since the paper was
+published. First, the options that have been removed are discussed,
+followed by those that have been added.
+
+@noindent
+Removed options:
+@itemize
+@item
+@option{--dilate}: In the paper, true detections were dilated for a final
+dig into the noise. However, simple 8-connected dilation can produce boxy
+results which are not realistic and could miss diffuse flux. The finalq dig
+into the noise is now done by ``grow''ing the true detections, similar to
+how true clumps were grown, see the description of @option{--detgrowquant}
+below and in @ref{Detection options} for more on the new alternative.
+@end itemize
 
+@noindent
+Added options:
 @itemize
 
 @item
@@ -12213,6 +12287,11 @@ give this option a value of 1 (only the largest valued 
pixel in the input
 will not be eroded).
 
 @item
+@option{--qthreshtilequant}: to manually remove the measured qthresh from
+some tiles. This feature helps in detecting large and extended diffuse
+(almost flat) signal when necessary, see @ref{Detection options}.
+
+@item
 @option{--detgrowquant}: is used to grow the final true detections until a
 given quantile in the same way that clumps are grown during segmentation
 (compare columns 2 and 3 in Figure 10 of the paper). It replaces the old
@@ -12226,6 +12305,11 @@ noise properties and how cleanly it was Sky 
subtracted. The new
 size to fill as part of this growth, see the description in @ref{Detection
 options} for more details.
 
+This new growth process can be much more successful in detecting diffuse
+flux around true detections compared to dilation and give more realistic
+results, but it can also increase the NoiseChisel run time (depending on
+the given value and input size).
+
 @item
 @option{--cleangrowndet}: A process to further clean/remove the possibility
 of false detections, see the descriptions under this option in
@@ -12233,11 +12317,6 @@ of false detections, see the descriptions under this 
option in
 
 @end itemize
 
-For a more detailed list of updates in each release, please follow the
-@file{NEWS} file. The @file{NEWS} file is in the released Gnuastro tarball
-(see @ref{Release tarball}). You can also see it online at
-@url{http://git.savannah.gnu.org/cgit/gnuastro.git/plain/NEWS}.
-
 @node Invoking astnoisechisel,  , NoiseChisel changes after publication, 
NoiseChisel
 @subsection Invoking NoiseChisel
 
@@ -12602,6 +12681,44 @@ is complete, we will have a binary (two valued) image. 
The pixels above the
 threshold are known as foreground pixels (have a value of 1) while those
 which lie below the threshold are known as background (have a value of 0).
 
+@item --qthreshtilequant=FLT
+Only keep tiles which have a q-thresh value above the given quantile of the
+all successful tiles. Hence, when given a value of @code{1}, this option
+will be ignored. When there is more than one channel (and
+@option{--workoverch} is not called), quantile calculation and application
+will be done on each channel independently.
+
+This option is useful when a large and diffuse (almost flat within each
+tile) signal exists with very small regions of Sky. The flat-ness of the
+profile will cause it to successfully pass the tests of @ref{Quantifying
+signal in a tile}. As a result, without this option the flat and diffuse
+signal will be interpretted as sky. In such cases, you can see the status
+of the tiles with the @option{--checkqthresh} option (first image extension
+is enough) and select a quantile through this option to ignore the measured
+values in the higher-valued tiles.
+
+This option can also be useful when there are large bright objects in the
+image with large flat wings which can also pass the tests and give outlier
+values. When there is a sky gradient over the image (mainly due to
+post-processing issues like bad flat fielding), this option must be set to
+@code{1} so it is completely ignored and the sky gradient is accurately
+measured and subtracted.
+
+To get an estimate of the measured qthresh distribution, you can run the
+following commands. The first will create the qthresh check image with one
+value/pixel per tile (see @ref{Processing options}). Open the image in a
+FITS viewer and inspect it. The second command below will print the basic
+information about the distribution of values and the third will print the
+value at the 0.4 quantile. Recall that Gnuastro's Statistics program
+ignores blank values (in this case: tiles with significant signal), see
+@ref{Statistics}.
+
+@example
+$ astnoisechisel image.fits --checkqthresh --oneelempertile
+$ aststatistics image_qthresh.fits
+$ aststatistics image_qthresh.fits --quantile=0.4
+@end example
+
 @item --smoothwidth=INT
 Width of flat kernel used to smooth the interpolated quantile thresholds,
 see @option{--qthresh} for more.
@@ -18261,23 +18378,14 @@ will be deleted.
 The minimum size of an array (in bytes) to store the contents of
 @code{array} as a file (on the non-volatile HDD/SSD), not in RAM. This can
 be very useful for large datasets which can be very memory intensive and
-the user's hardware RAM might not be sufficient to keep/process it. A random
+the user's RAM might not be sufficient to keep/process it. A random
 filename is assigned to the array which is available in the @code{mmapname}
-element of @code{gal_data_t} (above), see there for more.
-
-When this variable has a value of @code{0} (zero), any allocated
-@code{array} will actually be in a file (not in RAM). When the value is
-@code{-1} (largest possible number in the unsigned types including
-@code{size_t}) the array will be definitely allocated in RAM.
-
-Please note that using a non-volatile file instead of RAM will
-significantly increase the programs running time, especially on HDDs. So it
-is best to give this option very large values (depending on how much memory
-you will need for a given input). For example your processing might involve
-a copy of the the input (possibly to a wider data type which takes more
-bytes for each element), so take all such issues into
-consideration. @code{minmapsize} is actually stored in each
-@code{gal_data_t}, so it can be passed on to subsequent/derived datasets.
+element of @code{gal_data_t} (above), see there for more. @code{minmapsize}
+is stored in each @code{gal_data_t}, so it can be passed on to
+subsequent/derived datasets.
+
+See the description of the @option{--minmapsize} option in @ref{Processing
+options} for more on using this value.
 
 @item nwcs
 The number of WCS coordinate representations (for WCSLIB).
diff --git a/lib/data.c b/lib/data.c
index ba7ff40..17f09f4 100644
--- a/lib/data.c
+++ b/lib/data.c
@@ -179,7 +179,7 @@ gal_data_calloc_array(uint8_t type, size_t size, const char 
*funcname,
 
 
 static void
-gal_data_mmap(gal_data_t *data, int clear)
+gal_data_mmap(gal_data_t *data, int clear, size_t minmapsize)
 {
   int filedes;
   uint8_t uc=0;
@@ -222,8 +222,22 @@ gal_data_mmap(gal_data_t *data, int clear)
 
 
   /* Map the memory. */
+  errno=0;
   data->array=mmap(NULL, bsize, PROT_READ | PROT_WRITE, MAP_SHARED,
                    filedes, 0);
+  if(data->array==MAP_FAILED)
+    {
+      if(minmapsize<10000u)
+        fprintf(stderr, "\nIf the processing involves many small mappings "
+                "(along with larger ones), the following error may be "
+                "corrected with a larger value to `minmapsize' (minimum "
+                "number of bytes to use mapping instead of RAM for each "
+                "patch of memory), for example 10000. In this way, mapping "
+                "will only be reserved for larger sizes. The current value is "
+                "%zu.\n\n", minmapsize);
+      error(EXIT_FAILURE, errno, "couldn't map %zu bytes into the file `%s'",
+            bsize, filename);
+    }
 
 
   /* Close the file. */
@@ -318,6 +332,7 @@ gal_data_initialize(gal_data_t *data, void *array, uint8_t 
type,
           data->size *= ( data->dsize[i] = dsize[i] );
         }
 
+
       /* Set the array pointer. If an non-NULL array pointer was given,
          then use it. */
       if(array)
@@ -326,9 +341,9 @@ gal_data_initialize(gal_data_t *data, void *array, uint8_t 
type,
         {
           if(data->size)
             {
-              if( gal_type_sizeof(type)*data->size  > minmapsize )
+              if( gal_type_sizeof(type)*data->size > minmapsize )
                 /* Allocate the space into disk (HDD/SSD). */
-                gal_data_mmap(data, clear);
+                gal_data_mmap(data, clear, minmapsize);
               else
                 /* Allocate the space in RAM. */
                 data->array = ( clear
@@ -432,6 +447,9 @@ gal_data_free_contents(gal_data_t *data)
 
       /* Free the file name space. */
       free(data->mmapname);
+
+      /* Set the name pointer to NULL since it has been freed. */
+      data->mmapname=NULL;
     }
   else
     if(data->array && data->block==NULL) free(data->array);
@@ -505,6 +523,7 @@ gal_data_array_calloc(size_t size)
       out[i].wcs        = NULL;
       out[i].mmapname   = NULL;
       out[i].next       = NULL;
+      out[i].block      = NULL;
       out[i].name = out[i].unit = out[i].comment = NULL;
       out[i].disp_fmt = out[i].disp_width = out[i].disp_precision = -1;
     }
diff --git a/lib/gnuastro-internal/commonopts.h 
b/lib/gnuastro-internal/commonopts.h
index 78ad0d1..be69dee 100644
--- a/lib/gnuastro-internal/commonopts.h
+++ b/lib/gnuastro-internal/commonopts.h
@@ -164,7 +164,7 @@ struct argp_option gal_commonopts_options[] =
       GAL_OPTIONS_KEY_ONEELEMPERTILE,
       0,
       0,
-      "Only display 1 element/tile, not full input res.",
+      "Display 1 element/tile, not full input res.",
       GAL_OPTIONS_GROUP_TESSELLATION,
       &cp->tl.oneelempertile,
       GAL_OPTIONS_NO_ARG_TYPE,
diff --git a/lib/tile.c b/lib/tile.c
index 05c8eac..c972ea5 100644
--- a/lib/tile.c
+++ b/lib/tile.c
@@ -1244,7 +1244,7 @@ gal_tile_full_free_contents(struct 
gal_tile_two_layer_params *tl)
   if(tl->permutation)   free(tl->permutation);
   if(tl->firsttsize)    free(tl->firsttsize);
 
-  /* Free the arrays of `gal_data_c' for each tile and channel. */
+  /* Free the arrays of `gal_data_t' for each tile and channel. */
   if(tl->tiles)    gal_data_array_free(tl->tiles,    tl->tottiles,    0);
   if(tl->channels) gal_data_array_free(tl->channels, tl->totchannels, 0);
 }
diff --git a/lib/wcs.c b/lib/wcs.c
index 1fe18da..80167a1 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -438,10 +438,10 @@ double *
 gal_wcs_pixel_scale(struct wcsprm *wcs)
 {
   gsl_vector S;
-  double *a, *out;
-  int permute_set;
   gsl_matrix A, V;
   size_t i, j, n=wcs->naxis;
+  double *a, *out, maxrow, minrow;
+  int permute_set, warning_printed;
   double *v=gal_data_malloc_array(GAL_TYPE_FLOAT64, n*n, __func__, "v");
   size_t *permutation=gal_data_malloc_array(GAL_TYPE_SIZE_T, n, __func__,
                                             "permutation");
@@ -455,6 +455,62 @@ gal_wcs_pixel_scale(struct wcsprm *wcs)
   a=gal_wcs_warp_matrix(wcs);
 
 
+  /* To avoid confusing issues with floating point errors being written in
+     the non-diagonal elements of the FITS header PC or CD matrices, we
+     need to check if the minimum and maximum values in each row are not
+     several orders of magnitude apart.
+
+     Note that in some cases (for example a spectrum), one axis might be in
+     degrees (value around 1e-5) and the other in angestroms (value around
+     1e-10). So we can't look at the minimum and maximum of the whole
+     matrix. However, in such cases, people will probably not warp/rotate
+     the image to mix the coordinates. So the important thing to check is
+     the minimum and maximum (non-zero) values in each row. */
+  warning_printed=0;
+  for(i=0;i<n;++i)
+    {
+      /* Find the minimum and maximum values in each row. */
+      minrow=FLT_MAX;
+      maxrow=-FLT_MAX;
+      for(j=0;j<n;++j)
+        if(a[i*n+j]!=0.0) /* We aren't concerned with 0 valued elements. */
+          {
+            /* We have to use absolutes because in cases like RA, the
+               diagonal values in different rows may have different signs*/
+            if(fabs(a[i*n+j])<minrow) minrow=fabs(a[i*n+j]);
+            if(fabs(a[i*n+j])>maxrow) maxrow=fabs(a[i*n+j]);
+          }
+
+      /* Do the check, print warning and make correction. */
+      if(maxrow!=minrow && maxrow/minrow>1e4 && warning_printed==0)
+        {
+          fprintf(stderr, "\nWARNING: The input WCS matrix (possibly taken "
+                  "from the FITS header keywords starting with `CD' or `PC') "
+                  "contains values with very different scales (more than "
+                  "10^4 different). This is probably due to floating point "
+                  "errors. These values might bias the pixel scale (and "
+                  "subsequent) calculations.\n\n"
+                  "You can see the respective matrix with one of the "
+                  "following two commands (depending on how the FITS file "
+                  "was written). Recall that if the desired extension/HDU "
+                  "isn't the default, you can choose it with the `--hdu' "
+                  "(or `-h') option before the `|' sign in these commands."
+                  "\n\n"
+                  "    $ astfits file.fits -p | grep 'PC._.'\n"
+                  "    $ astfits file.fits -p | grep 'CD._.'\n\n"
+                  "You can delete the ones with obvious floating point "
+                  "error values using the following command (assuming you "
+                  "want to delete `CD1_2' and `CD2_1'). Afterwards, you can "
+                  "rerun your original command to remove this warning "
+                  "message and possibly correct errors that it might have "
+                  "caused.\n\n"
+                  "    $ astfits file.fits --delete=CD1_2 --delete=CD2_1\n\n"
+                  );
+          warning_printed=1;
+        }
+    }
+
+
   /* Fill in the necessary GSL vector and Matrix structures. */
   S.size=n;     S.stride=1;                S.data=pixscale->array;
   V.size1=n;    V.size2=n;    V.tda=n;     V.data=v;
diff --git a/tmpfs-config-make b/tmpfs-config-make
index 019de25..41cc305 100755
--- a/tmpfs-config-make
+++ b/tmpfs-config-make
@@ -30,6 +30,7 @@
 # Set the variables:
 NUM_THREADS=8
 TMPDIR=/dev/shm
+debug_noshared=0
 
 
 
@@ -123,7 +124,11 @@ cd $build_dir
 # (in `configure.ac'), we have set optimization flags which have to be
 # cancelled in debugging.
 if [ ! -f Makefile ]; then
-    $srcdir/configure --srcdir=$srcdir #CFLAGS="-g -O0" --disable-shared
+    if [ x$debug_noshared = x1 ]; then
+        $srcdir/configure --srcdir=$srcdir CFLAGS="-g -O0" --disable-shared
+    else
+        $srcdir/configure --srcdir=$srcdir
+    fi
 fi
 
 



reply via email to

[Prev in Thread] Current Thread [Next in Thread]