[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 800c769 082/113: Imported all the work in mast
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 800c769 082/113: Imported all the work in master, conflicts fixed |
|
Date: |
Fri, 16 Apr 2021 10:33:53 -0400 (EDT) |
branch: master
commit 800c769e666ae9a78c0749fec503a72120ff04a1
Merge: 8d64628 72e46dc
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Imported all the work in master, conflicts fixed
Several conflicts came up during the merge in `bin/mkprof/ui.c',
`bin/noisechisel/detection.c', and `bin/noisechisel/threshold.c'. Besides,
those, the warning messages of commit 62738cc (MakeProfiles only checking
validity of column when relevant) were also added in the 3D profile
checks. Finally, the configuration file of NoiseChisel in 3D was also
updated.
---
NEWS | 61 +-
THANKS | 3 +
bin/arithmetic/args.h | 18 +
bin/arithmetic/arithmetic.c | 6 +-
bin/arithmetic/main.h | 1 +
bin/arithmetic/ui.c | 1 -
bin/arithmetic/ui.h | 3 +-
bin/fits/fits.c | 2 +-
bin/gnuastro.conf | 2 +-
bin/mkprof/args.h | 4 +-
bin/mkprof/main.h | 1 +
bin/mkprof/mkprof.c | 113 +-
bin/mkprof/mkprof.h | 1 +
bin/mkprof/oneprofile.c | 8 +-
bin/mkprof/ui.c | 52 +-
bin/mkprof/ui.h | 3 +
bin/noisechisel/args.h | 73 +-
bin/noisechisel/astnoisechisel-3d.conf | 49 +-
bin/noisechisel/astnoisechisel.conf | 46 +-
bin/noisechisel/detection.c | 40 +-
bin/noisechisel/main.h | 6 +-
bin/noisechisel/noisechisel.c | 10 +-
bin/noisechisel/sky.c | 115 +-
bin/noisechisel/threshold.c | 236 +---
bin/noisechisel/threshold.h | 5 +
bin/noisechisel/ui.c | 1 +
bin/noisechisel/ui.h | 8 +-
bin/statistics/args.h | 33 +-
bin/statistics/aststatistics.conf | 4 +-
bin/statistics/main.h | 4 +-
bin/statistics/sky.c | 47 +-
bin/statistics/ui.c | 18 +-
bin/statistics/ui.h | 5 +-
doc/announce-acknowledge.txt | 4 +
doc/gnuastro.texi | 1362 +++++++++++++-------
lib/Makefile.am | 6 +-
lib/blank.c | 436 ++++---
lib/gnuastro-internal/arithmetic-internal.h | 2 +-
.../{arithmetic-internal.h => tile-internal.h} | 27 +-
lib/gnuastro/blank.h | 3 +
lib/gnuastro/statistics.h | 5 +-
lib/statistics.c | 148 ++-
lib/tile-internal.c | 191 +++
tests/during-dev.sh | 4 +-
44 files changed, 2007 insertions(+), 1160 deletions(-)
diff --git a/NEWS b/NEWS
index 78e8e16..52f7c28 100644
--- a/NEWS
+++ b/NEWS
@@ -9,30 +9,75 @@ GNU Astronomy Utilities NEWS -*-
outline -*-
of output when it is FITS. This is only relevant for some programs,
(for example not the Fits or Table programs).
+ Arithmetic:
+ --onedasimage: write output as an image if it has one dimension, not table.
+
+ NoiseChisel:
+ - New outlier identification algorithm for quantile threshold and Sky
+ estimation. It is useful when there are extended and bright sources
+ in the dataset: the tiles containing very faint signal that pass the
+ general pixel-value distribution test due to the flatness of the
+ extended profiles (outliers), can be identified and removed in
+ comparison with the other passed tiles. The outlier finding algorithm
+ (`gal_statistics_outlier_positive': a new library function) uses the
+ distribution of distances between the sorted elements and is
+ configured with these three options.
+ --outliersclip: Sigma-clipping parameters for the process.
+ --outliersigma: Multiple of sigma to define an outlier.
+
+ Statistics:
+ - Sky estimation: new outlier estimation algorithm similar to NoiseChisel.
+
Library:
- -- gal_fits_key_list_reverse: Reverse the given list of FITS keywords.
- -- gal_fits_key_write_title_in_ptr: Write a two-line title in FITS
keywords.
- -- gal_fits_key_write_in_ptr: New name of gal_fits_key_write.
- -- gal_fits_key_write_version_in_ptr: new name of
gal_fits_key_write_version.
- -- gal_fits_key_write_config: write key list and version as config header.
+ - gal_blank_flag_apply: Set all flagged/masked elements to blank.
+ - gal_fits_key_list_reverse: Reverse the given list of FITS keywords.
+ - gal_fits_key_write_title_in_ptr: Write a two-line title FITS keyword.
+ - gal_fits_key_write_in_ptr: New name of gal_fits_key_write.
+ - gal_fits_key_write_version_in_ptr: old gal_fits_key_write_version.
+ - gal_fits_key_write_config: write key list and version as config header.
+ - gal_statistics_outlier_positive: Find the first positive outlier.
** Removed features
+ NoiseChisel:
+ --mirrordist: mean quantile is now used (not mode), see changes below.
+ --qthreshtilequant: removed due to new outlier rejection algorithm.
+
+ Statistics:
+ --mirrordist: mean quantile is now used (not mode), see changes below.
+
** Changed features
+ Arithmetic:
+ - If the output has one dimension, it will be written as a table, not a
+ FITS image/array. This can be changed with the new `--onedasimage'
+ option.
+
MakeProfiles:
--mergedsize: new name for the old `--naxis' option. Since the option
names and values are now written into the FITS header of the
output, this option's name would get confused with the mandatory
FITS keyword `NAXIS'.
+ NoiseChisel:
+ --meanmedqdiff: new name for `--modmedqdiff'. Until now, the mode's
+ quantile was used to identify tiles with no significant signal. But
+ from this version, the mean's quantile in each tile is used
+ instead. The reason is that the mean is more sensative to outliers
+ (signal).
+
+ Statistics:
+ --meanmedqdiff: new name for `--modmedqdiff'. Similar to NoiseChisel.
+
Library:
- -- gal_fits_key_write: filename and hdu instead of FITS pointer.
- -- gal_fits_key_write_version: filename and hdu instead of FITS pointer.
- -- gal_fits_key_write_filename: write at the top or end of the input list.
+ - gal_fits_key_write: filename and hdu instead of FITS pointer.
+ - gal_fits_key_write_version: filename and hdu instead of FITS pointer.
+ - gal_fits_key_write_filename: write at the top or end of the input list.
** Bugs fixed
bug #54493: Warp crashes when type isn't set.
+ bug #54526: Invalid r, q and truncation of point profiles in MakeProfiles.
+ bug #54579: NoiseChisel pseudo-detection failure when dataset is negative.
diff --git a/THANKS b/THANKS
index 1dde839..f55b913 100644
--- a/THANKS
+++ b/THANKS
@@ -28,6 +28,8 @@ support in Gnuastro. The list is ordered alphabetically (by
family name).
Benjamin Clement benjamin.clement@univ-lyon1.fr
Nima Dehdilani nimadehdilani@gmail.com
Antonio Diaz Diaz antonio@gnu.org
+ Pierre-Alain Duc pierre-alain.duc@astro.unistra.fr
+ Gaspar Galaz ggalaz@astro.puc.cl
Thérèse Godefroy godef.th@free.fr
Madusha Gunawardhana gunawardhana@strw.leidenuniv.nl
Stephen Hamer stephen.hamer@univ-lyon1.fr
@@ -55,6 +57,7 @@ support in Gnuastro. The list is ordered alphabetically (by
family name).
Jenny Sorce jenny.sorce@univ-lyon1.fr
Lee Spitler lee.spitler@mq.edu.au
Richard Stallman rms@gnu.org
+ Michael Stein mstein@astro.rub.de
Ole Streicher olebole@debian.org
Alfred M. Szmidt ams@gnu.org
Michel Tallon mtallon@obs.univ-lyon1.fr
diff --git a/bin/arithmetic/args.h b/bin/arithmetic/args.h
index 8ba927f..ff55ef2 100644
--- a/bin/arithmetic/args.h
+++ b/bin/arithmetic/args.h
@@ -42,6 +42,24 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET
},
+
+
+
+
+ {
+ "onedasimage",
+ UI_KEY_ONEDASIMAGE,
+ 0,
+ 0,
+ "Write 1D outputs as an image, not a table.",
+ GAL_OPTIONS_GROUP_OUTPUT,
+ &p->onedasimage,
+ GAL_OPTIONS_NO_ARG_TYPE,
+ GAL_OPTIONS_RANGE_0_OR_1,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+
{0}
};
diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index e511e69..8afcb7b 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -1169,7 +1169,11 @@ reversepolish(struct arithmeticparams *p)
/* Put a copy of the WCS structure from the reference image, it
will be freed while freeing d1. */
d1->wcs=p->refdata.wcs;
- gal_fits_img_write(d1, p->cp.output, NULL, PROGRAM_NAME);
+ if(d1->ndim==1 && p->onedasimage==0)
+ gal_table_write(d1, NULL, p->cp.tableformat, p->cp.output,
+ "ARITHMETIC", 0);
+ else
+ gal_fits_img_write(d1, p->cp.output, NULL, PROGRAM_NAME);
if(!p->cp.quiet)
printf(" - Output written to %s\n", p->cp.output);
}
diff --git a/bin/arithmetic/main.h b/bin/arithmetic/main.h
index 652f7e8..ce36ebb 100644
--- a/bin/arithmetic/main.h
+++ b/bin/arithmetic/main.h
@@ -74,6 +74,7 @@ struct arithmeticparams
size_t popcounter; /* The number of FITS images popped. */
gal_data_t refdata; /* Container for information of the data. */
char *globalhdu; /* Single HDU for all inputs. */
+ uint8_t onedasimage; /* Write 1D outputs as an image not table. */
gal_data_t *named; /* List containing variables. */
size_t tokencounter; /* Counter for finding place in tokens. */
diff --git a/bin/arithmetic/ui.c b/bin/arithmetic/ui.c
index 3824f3b..9ca9e6e 100644
--- a/bin/arithmetic/ui.c
+++ b/bin/arithmetic/ui.c
@@ -140,7 +140,6 @@ ui_initialize_options(struct arithmeticparams *p,
case GAL_OPTIONS_KEY_TYPE:
case GAL_OPTIONS_KEY_SEARCHIN:
case GAL_OPTIONS_KEY_IGNORECASE:
- case GAL_OPTIONS_KEY_TABLEFORMAT:
cp->coptions[i].flags=OPTION_HIDDEN;
break;
diff --git a/bin/arithmetic/ui.h b/bin/arithmetic/ui.h
index ed36277..c4fed45 100644
--- a/bin/arithmetic/ui.h
+++ b/bin/arithmetic/ui.h
@@ -33,12 +33,13 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
/* Available letters for short options:
a b c d e f i j k l m n p r s t u v w x y z
- A B C E G H J L O Q R W X Y
+ A B C E G H J L Q R W X Y
*/
enum option_keys_enum
{
/* With short-option version. */
UI_KEY_GLOBALHDU = 'g',
+ UI_KEY_ONEDASIMAGE = 'O',
/* Only with long version (start with a value 1000, the rest will be set
automatically). */
diff --git a/bin/fits/fits.c b/bin/fits/fits.c
index 320ec50..55c278b 100644
--- a/bin/fits/fits.c
+++ b/bin/fits/fits.c
@@ -140,7 +140,7 @@ fits_print_extension_info(struct fitsparams *p)
{
case IMAGE_HDU:
gal_fits_img_info(fptr, &type, &ndim, &dsize, NULL, NULL);
- tstr=gal_type_name(type , 1);
+ tstr = ndim==0 ? "no-data" : gal_type_name(type , 1);
break;
case ASCII_TBL:
diff --git a/bin/gnuastro.conf b/bin/gnuastro.conf
index 472879f..c9a6d42 100644
--- a/bin/gnuastro.conf
+++ b/bin/gnuastro.conf
@@ -23,7 +23,7 @@
searchin name
# Tessellation
- tilesize 50,50
+ tilesize 30,30
numchannels 1,1
remainderfrac 0.1
workoverch 0
diff --git a/bin/mkprof/args.h b/bin/mkprof/args.h
index 6d9c224..4ed95f1 100644
--- a/bin/mkprof/args.h
+++ b/bin/mkprof/args.h
@@ -442,7 +442,7 @@ struct argp_option program_options[] =
UI_KEY_QCOL,
"STR/INT",
0,
- "Axis ratio (major/dim2 radius).",
+ "Axis ratio (major/dim2 in 3D).",
UI_GROUP_CATALOG,
&p->qcol,
GAL_TYPE_STRING,
@@ -455,7 +455,7 @@ struct argp_option program_options[] =
UI_KEY_Q2COL,
"STR/INT",
0,
- "Axis ratio (major/dim3 radius).",
+ "Axis ratio (major/dim3 in 3D).",
UI_GROUP_CATALOG,
&p->q2col,
GAL_TYPE_STRING,
diff --git a/bin/mkprof/main.h b/bin/mkprof/main.h
index e71a125..ededc0f 100644
--- a/bin/mkprof/main.h
+++ b/bin/mkprof/main.h
@@ -176,6 +176,7 @@ struct mkprofparams
float *m; /* Magnitude of profile. */
float *t; /* Truncation distance. */
gsl_rng *rng; /* Main instance of random number generator.*/
+ const char *rng_name; /* Name of random number generator. */
time_t rawtime; /* Starting time of the program. */
double *cat; /* Input catalog. */
gal_data_t *log; /* Log data to be printed. */
diff --git a/bin/mkprof/mkprof.c b/bin/mkprof/mkprof.c
index 409c786..fba9ce4 100644
--- a/bin/mkprof/mkprof.c
+++ b/bin/mkprof/mkprof.c
@@ -42,6 +42,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include "main.h"
+#include "ui.h"
#include "mkprof.h" /* Needs main.h astrthreads.h */
#include "oneprofile.h"
@@ -127,8 +128,9 @@ saveindividual(struct mkonthread *mkp)
double *crpix;
long os=p->oversample;
- size_t i, ndim=p->ndim;
+ gal_fits_list_key_t *keys=NULL;
struct builtqueue *ibq=mkp->ibq;
+ size_t i, ndim=p->ndim, id=mkp->ibq->id;
char *filename, *jobname, *outdir=p->outdir;
@@ -145,7 +147,7 @@ saveindividual(struct mkonthread *mkp)
}
- /* Write the array to file (a separately built PSF doesn't need WCS
+ /* Write the array to the file (a separately built PSF doesn't need WCS
coordinates). */
if(ibq->ispsf && p->psfinimg==0)
gal_fits_img_write(ibq->image, filename, NULL, PROGRAM_NAME);
@@ -166,6 +168,104 @@ saveindividual(struct mkonthread *mkp)
ibq->indivcreated=1;
+ /* Write profile settings into the FITS file. */
+ gal_fits_key_list_add(&keys, GAL_TYPE_STRING, "PROFILE", 0,
+ ui_profile_name_write(mkp->func), 0,
+ "Radial function", 0, NULL);
+ gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT64, "XCENTER", 0,
+ &p->x[id], 0, "Center of profile in catalog "
+ "(FITS axis 1)", 0, NULL);
+ gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT64, "YCENTER", 0,
+ &p->y[id], 0, "Center of profile in catalog "
+ "(FITS axis 2)", 0, NULL);
+ if(ndim==3)
+ gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT64, "ZCENTER", 0,
+ &p->z[id], 0, "Center of profile in catalog "
+ "(FITS axis 3)", 0, NULL);
+ gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "RADIUS", 0,
+ &p->r[id], 0, "Radial parameter in catalog",
+ 0, NULL);
+ if( mkp->func==PROFILE_SERSIC || mkp->func==PROFILE_MOFFAT )
+ gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "PINDEX", 0,
+ &p->r[id], 0, "Index (Sersic or Moffat) of profile"
+ "in catalog", 0, NULL);
+ if(ndim==2)
+ {
+ gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "PA_DEG", 0,
+ &p->p1[id], 0, "Position angle of profile in "
+ "catalog", 0, "deg");
+ gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "AXISRATIO", 0,
+ &p->q1[id], 0, "Axis ratio of profile in catalog",
+ 0, NULL);
+ }
+ else
+ {
+ gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "PA1_DEG", 0,
+ &p->p1[id], 0, "First X-Z-X Euler angle in 3D", 0,
+ "deg");
+ gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "PA2_DEG", 0,
+ &p->p2[id], 0, "Second X-Z-X Euler angle in 3D", 0,
+ "deg");
+ gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "PA3_DEG", 0,
+ &p->p3[id], 0, "Third X-Z-X Euler angle in 3D", 0,
+ "deg");
+ gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "AXISRATIO1", 0,
+ &p->q1[id], 0, "Axis ratio along second dim",
+ 0, NULL);
+ gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "AXISRATIO2", 0,
+ &p->q2[id], 0, "Axis ratio along third dim",
+ 0, NULL);
+ }
+ gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "MAGNITUDE", 0,
+ &p->m[id], 0, "Magnitude of profile in catalog",
+ 0, NULL);
+ gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "TRUNCATION", 0,
+ &p->t[id], 0, "Truncation of profile in catalog",
+ 0, NULL);
+ gal_fits_key_list_add(&keys, GAL_TYPE_LONG, "RNGSEED", 0,
+ &mkp->rng_seed, 0, "Seed of random number generator",
+ 0, NULL);
+ gal_fits_key_list_add(&keys, GAL_TYPE_STRING, "RNGNAME", 0,
+ (void *)(p->rng_name), 0,
+ "Name of random number generator", 0, NULL);
+ gal_fits_key_list_add(&keys, GAL_TYPE_SIZE_T, "NUMRANDOM", 0,
+ &p->numrandom, 0,
+ "Number of random points in central pixels", 0, NULL);
+ gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "TOLERANCE", 0,
+ &p->tolerance, 0,
+ "Tolerance level to stop random integration",
+ 0, NULL);
+ gal_fits_key_list_add(&keys, GAL_TYPE_STRING, "MODE", 0,
+ p->mode==MKPROF_MODE_IMG?"img":"wcs", 0,
+ "Coordinates in image or WCS units", 0, NULL);
+ gal_fits_key_list_add(&keys, GAL_TYPE_UINT8, "OVERSAMPLE", 0,
+ &p->oversample, 0, "Oversampling factor", 0, NULL);
+ gal_fits_key_list_add(&keys, GAL_TYPE_UINT8, "TUNITINP", 0,
+ &p->tunitinp, 0, "Truncation is in units of pixels, "
+ "not radius", 0, NULL);
+ if( !isnan(p->zeropoint) )
+ gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "ZEROPOINT", 0,
+ &p->zeropoint, 0, "Zeropoint magnitude", 0, NULL);
+ if( mkp->func==PROFILE_CIRCUMFERENCE )
+ gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "CIRCUMWIDTH", 0,
+ &p->circumwidth, 0, "Width of circumference "
+ "(inward) profiles", 0, NULL);
+ if( mkp->func==PROFILE_FLAT || mkp->func==PROFILE_CIRCUMFERENCE )
+ gal_fits_key_list_add(&keys, GAL_TYPE_UINT8, "MFORFLATPIX", 0,
+ &p->mforflatpix, 0, "Magnitude is flat pixel "
+ "value", 0, NULL);
+ gal_fits_key_list_add(&keys, GAL_TYPE_UINT8, "MCOLISBRIGHTNESS", 0,
+ &p->mcolisbrightness, 0, "Catalog's magnitude is "
+ "actually brightness", 0, NULL);
+ gal_fits_key_list_add(&keys, GAL_TYPE_UINT8, "MAGATPEAK", 0,
+ &p->magatpeak, 0, "Magnitude is for peak pixel, "
+ "not full profile", 0, NULL);
+
+ gal_fits_key_list_reverse(&keys);
+ gal_fits_key_write_config(&keys, "Profile configuration", "PROFILE-CONFIG",
+ filename, "0");
+
+
/* Report if in verbose mode. */
if(!p->cp.quiet)
{
@@ -221,8 +321,13 @@ mkprof_build_single(struct mkonthread *mkp, long
*fpixel_i, long *lpixel_i,
/* Set the seed of the random number generator if the
environment is not to be used. */
- if(mkp->p->envseed==0)
- gsl_rng_set(mkp->rng, gal_timing_time_based_rng_seed());
+ if(mkp->p->envseed)
+ mkp->rng_seed=mkp->p->envseed;
+ else
+ {
+ mkp->rng_seed=gal_timing_time_based_rng_seed();
+ gsl_rng_set(mkp->rng, mkp->rng_seed);
+ }
/* Make the profile */
oneprofile_make(mkp);
diff --git a/bin/mkprof/mkprof.h b/bin/mkprof/mkprof.h
index 4eef34e..c97fcdc 100644
--- a/bin/mkprof/mkprof.h
+++ b/bin/mkprof/mkprof.h
@@ -48,6 +48,7 @@ struct mkonthread
long *onaxes; /* Sides of the unover-sampled image. */
long fpixel_i[3]; /* fpixel_i before running overlap. */
int correction; /* ==1: correct the pixels afterwards. */
+ long rng_seed; /* Seed used to generate this profile. */
/* Random number generator: */
gsl_rng *rng; /* Copy of main random number generator. */
diff --git a/bin/mkprof/oneprofile.c b/bin/mkprof/oneprofile.c
index 06935ab..49d4d69 100644
--- a/bin/mkprof/oneprofile.c
+++ b/bin/mkprof/oneprofile.c
@@ -693,10 +693,12 @@ oneprof_set_prof_params(struct mkonthread *mkp)
mkp->correction = 0;
break;
+
+
default:
- error(EXIT_FAILURE, 0, "%s: a bug! Please contact us so we can correct "
- "this problem. The profile code %u is not recognized.", __func__,
- mkp->func);
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us so we can "
+ "correct this problem. The profile code %u is not recognized.",
+ __func__, mkp->func);
}
}
diff --git a/bin/mkprof/ui.c b/bin/mkprof/ui.c
index 495f61f..d16d90c 100644
--- a/bin/mkprof/ui.c
+++ b/bin/mkprof/ui.c
@@ -144,7 +144,7 @@ ui_profile_name_read(char *string, size_t row)
-static char *
+char *
ui_profile_name_write(int profile_code)
{
switch(profile_code)
@@ -729,10 +729,11 @@ ui_read_cols_2d(struct mkprofparams *p)
/* Check if there is no negative or zero-radius profile. */
for(i=0;i<p->num;++i)
- if(p->r[i]<=0.0f)
+ if(p->f[i]!=PROFILE_POINT && p->r[i]<=0.0f)
error(EXIT_FAILURE, 0, "%s: row %zu, the radius value %g is "
- "not acceptable. It has to be larger than 0", p->catname,
- i+1, p->r[i]);
+ "not acceptable for a `%s' profile. It has to be larger "
+ "than 0", p->catname, i+1, p->r[i],
+ ui_profile_name_write(p->f[i]));
break;
@@ -757,10 +758,11 @@ ui_read_cols_2d(struct mkprofparams *p)
/* Check if there is no negative or >1.0f axis ratio. */
for(i=0;i<p->num;++i)
- if(p->q1[i]<=0.0f || p->q1[i]>1.0f)
+ if( p->f[i]!=PROFILE_POINT && (p->q1[i]<=0.0f || p->q1[i]>1.0f) )
error(EXIT_FAILURE, 0, "%s: row %zu, the axis ratio value %g "
- "is not acceptable. It has to be >0 and <=1", p->catname,
- i+1, p->q1[i]);
+ "is not acceptable for a `%s' profile. It has to be >0 "
+ "and <=1", p->catname, i+1, p->q1[i],
+ ui_profile_name_write(p->f[i]));
break;
@@ -779,10 +781,11 @@ ui_read_cols_2d(struct mkprofparams *p)
/* Check if there is no negative or zero truncation radius. */
for(i=0;i<p->num;++i)
- if(p->t[i]<=0.0f)
+ if(p->f[i]!=PROFILE_POINT && p->t[i]<=0.0f)
error(EXIT_FAILURE, 0, "%s: row %zu, the truncation radius "
- "value %g is not acceptable. It has to be larger than 0",
- p->catname, i+1, p->t[i]);
+ "value %g is not acceptable for a `%s' profile. It has "
+ "to be larger than 0", p->catname, i+1, p->t[i],
+ ui_profile_name_write(p->f[i]));
break;
@@ -936,6 +939,13 @@ ui_read_cols_3d(struct mkprofparams *p)
colname="radius (`rcol')";
corrtype=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT32);
p->r=corrtype->array;
+
+ /* Check if there is no negative or zero-radius profile. */
+ if(p->f[i]!=PROFILE_POINT && p->r[i]<=0.0f)
+ error(EXIT_FAILURE, 0, "%s: row %zu, the radius value %g is "
+ "not acceptable for a `%s' profile. It has to be larger "
+ "than 0", p->catname, i+1, p->r[i],
+ ui_profile_name_write(p->f[i]));
break;
case 6:
@@ -969,10 +979,11 @@ ui_read_cols_3d(struct mkprofparams *p)
/* Check if there is no negative or >1.0f axis ratio. */
for(i=0;i<p->num;++i)
- if(p->q1[i]<=0.0f || p->q1[i]>1.0f)
+ if( p->f[i]!=PROFILE_POINT && (p->q1[i]<=0.0f || p->q1[i]>1.0f) )
error(EXIT_FAILURE, 0, "%s: row %zu, the first axis ratio "
- "(`--qcol') value %g is not acceptable. It has to be "
- ">0 and <=1", p->catname, i+1, p->q1[i]);
+ "value %g is not acceptable for a `%s' profile. It has "
+ "to be >0 and <=1", p->catname, i+1, p->q1[i],
+ ui_profile_name_write(p->f[i]));
break;
case 11:
@@ -982,10 +993,11 @@ ui_read_cols_3d(struct mkprofparams *p)
/* Check if there is no negative or >1.0f axis ratio. */
for(i=0;i<p->num;++i)
- if(p->q2[i]<=0.0f || p->q2[i]>1.0f)
+ if( p->f[i]!=PROFILE_POINT && (p->q2[i]<=0.0f || p->q2[i]>1.0f) )
error(EXIT_FAILURE, 0, "%s: row %zu, the second axis ratio "
- "(`--q2col') value %g is not acceptable. It has to be "
- ">0 and <=1", p->catname, i+1, p->q2[i]);
+ "value %g is not acceptable for a `%s' profile. It has "
+ "to be >0 and <=1", p->catname, i+1, p->q2[i],
+ ui_profile_name_write(p->f[i]));
break;
case 12:
@@ -1002,10 +1014,11 @@ ui_read_cols_3d(struct mkprofparams *p)
/* Check if there is no negative or zero truncation radius. */
for(i=0;i<p->num;++i)
- if(p->t[i]<=0.0f)
+ if(p->f[i]!=PROFILE_POINT && p->t[i]<=0.0f)
error(EXIT_FAILURE, 0, "%s: row %zu, the truncation radius "
- "value %g is not acceptable. It has to be larger than 0",
- p->catname, i+1, p->t[i]);
+ "value %g is not acceptable for a `%s' profile. It has "
+ "to be larger than 0", p->catname, i+1, p->t[i],
+ ui_profile_name_write(p->f[i]));
break;
/* If the index isn't recognized, then it is larger, showing that
@@ -1672,6 +1685,7 @@ ui_preparations(struct mkprofparams *p)
/* Allocate the random number generator: */
gsl_rng_env_setup();
p->rng=gsl_rng_alloc(gsl_rng_ranlxs1);
+ p->rng_name=gsl_rng_name(p->rng);
/* Make the log linked list. */
ui_make_log(p);
diff --git a/bin/mkprof/ui.h b/bin/mkprof/ui.h
index 446071a..7cdd06e 100644
--- a/bin/mkprof/ui.h
+++ b/bin/mkprof/ui.h
@@ -102,6 +102,9 @@ enum option_keys_enum
void
ui_read_check_inputs_setup(int argc, char *argv[], struct mkprofparams *p);
+char *
+ui_profile_name_write(int profile_code);
+
void
ui_free_report(struct mkprofparams *p, struct timeval *t1);
diff --git a/bin/noisechisel/args.h b/bin/noisechisel/args.h
index e96b73d..4236fee 100644
--- a/bin/noisechisel/args.h
+++ b/bin/noisechisel/args.h
@@ -186,26 +186,13 @@ struct argp_option program_options[] =
UI_GROUP_DETECTION
},
{
- "mirrordist",
- UI_KEY_MIRRORDIST,
+ "meanmedqdiff",
+ UI_KEY_MEANMEDQDIFF,
"FLT",
0,
- "Max. dist. (error multip.) to find mode.",
+ "Max. mean and median quant diff. per tile.",
UI_GROUP_DETECTION,
- &p->mirrordist,
- GAL_TYPE_FLOAT32,
- GAL_OPTIONS_RANGE_ANY,
- GAL_OPTIONS_MANDATORY,
- GAL_OPTIONS_NOT_SET
- },
- {
- "modmedqdiff",
- UI_KEY_MODMEDQDIFF,
- "FLT",
- 0,
- "Max. mode and median quant diff. per tile.",
- UI_GROUP_DETECTION,
- &p->modmedqdiff,
+ &p->meanmedqdiff,
GAL_TYPE_FLOAT32,
GAL_OPTIONS_RANGE_GE_0,
GAL_OPTIONS_MANDATORY,
@@ -225,15 +212,29 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET
},
{
- "qthreshtilequant",
- UI_KEY_QTHRESHTILEQUANT,
+ "outliersclip",
+ UI_KEY_OUTLIERSCLIP,
+ "FLT,FLT",
+ 0,
+ "Sigma-clip params for qthresh outliers.",
+ UI_GROUP_DETECTION,
+ &p->outliersclip,
+ GAL_TYPE_STRING,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ gal_options_read_sigma_clip
+ },
+ {
+ "outliersigma",
+ UI_KEY_OUTLIERSIGMA,
"FLT",
0,
- "Remove tiles at higher quantiles.",
+ "Multiple of sigma to define outliers.",
UI_GROUP_DETECTION,
- &p->qthreshtilequant,
+ &p->outliersigma,
GAL_TYPE_FLOAT32,
- GAL_OPTIONS_RANGE_GE_0_LE_1,
+ GAL_OPTIONS_RANGE_GE_0,
GAL_OPTIONS_MANDATORY,
GAL_OPTIONS_NOT_SET
},
@@ -329,20 +330,6 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET
},
{
- "sigmaclip",
- UI_KEY_SIGMACLIP,
- "FLT,FLT",
- 0,
- "Sigma multiple and, tolerance or number.",
- UI_GROUP_DETECTION,
- &p->sigmaclip,
- GAL_TYPE_STRING,
- GAL_OPTIONS_RANGE_ANY,
- GAL_OPTIONS_MANDATORY,
- GAL_OPTIONS_NOT_SET,
- gal_options_read_sigma_clip
- },
- {
"minskyfrac",
UI_KEY_MINSKYFRAC,
"FLT",
@@ -369,6 +356,20 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET
},
{
+ "sigmaclip",
+ UI_KEY_SIGMACLIP,
+ "FLT,FLT",
+ 0,
+ "Sigma multiple and, tolerance or number.",
+ UI_GROUP_DETECTION,
+ &p->sigmaclip,
+ GAL_TYPE_STRING,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ gal_options_read_sigma_clip
+ },
+ {
"dthresh",
UI_KEY_DTHRESH,
"FLT",
diff --git a/bin/noisechisel/astnoisechisel-3d.conf
b/bin/noisechisel/astnoisechisel-3d.conf
index 46c3bbf..5a80eb5 100644
--- a/bin/noisechisel/astnoisechisel-3d.conf
+++ b/bin/noisechisel/astnoisechisel-3d.conf
@@ -18,33 +18,34 @@
# warranty.
# Input:
- khdu 1
- whdu 1
- minskyfrac 0.7
- minnumfalse 100
+ khdu 1
+ whdu 1
+ chdu 1
+ minskyfrac 0.7
+ minnumfalse 100
# Tessellation
- numchannels 1,1,1
- tilesize 15,15,15
- largetilesize 50,50,50
+ numchannels 1,1,1
+ tilesize 15,15,15
+ largetilesize 50,50,50
# Detection:
- mirrordist 1.5
- modmedqdiff 0.01
- qthresh 0.3
- qthreshtilequant 1
- smoothwidth 3
- erode 1
- erodengb 6
- noerodequant 0.9331
- opening 1
- openingngb 18
- sigmaclip 3,0.2
- dthresh -0.1
- snminarea 10
- snquant 0.90
- detgrowquant 0.95
- detgrowmaxholesize 300
+ meanmedqdiff 0.005
+ qthresh 0.3
+ outliersigma 10
+ outliersclip 3,0.2
+ smoothwidth 3
+ erode 1
+ erodengb 6
+ noerodequant 0.99
+ opening 1
+ openingngb 18
+ sigmaclip 3,0.2
+ dthresh 0.0
+ snminarea 20
+ snquant 0.99
+ detgrowquant 0.95
+ detgrowmaxholesize 300
# Operating mode
- continueaftercheck 0
+ continueaftercheck 0
diff --git a/bin/noisechisel/astnoisechisel.conf
b/bin/noisechisel/astnoisechisel.conf
index 5126341..54209cf 100644
--- a/bin/noisechisel/astnoisechisel.conf
+++ b/bin/noisechisel/astnoisechisel.conf
@@ -18,32 +18,32 @@
# warranty.
# Input:
- khdu 1
- whdu 1
- chdu 1
- minskyfrac 0.7
- minnumfalse 100
+ khdu 1
+ whdu 1
+ chdu 1
+ minskyfrac 0.7
+ minnumfalse 100
# Tessellation
- largetilesize 200,200
+ largetilesize 200,200
# Detection:
- mirrordist 1.5
- modmedqdiff 0.01
- qthresh 0.3
- qthreshtilequant 1.0
- smoothwidth 3
- erode 2
- erodengb 4
- noerodequant 0.9331
- opening 1
- openingngb 8
- sigmaclip 3,0.2
- dthresh 0.0
- snminarea 10
- snquant 0.95
- detgrowquant 0.80
- detgrowmaxholesize 100
+ meanmedqdiff 0.005
+ qthresh 0.3
+ outliersigma 10
+ outliersclip 3,0.2
+ smoothwidth 3
+ erode 2
+ erodengb 4
+ noerodequant 0.99
+ opening 1
+ openingngb 8
+ sigmaclip 3,0.2
+ dthresh 0.0
+ snminarea 10
+ snquant 0.99
+ detgrowquant 0.90
+ detgrowmaxholesize 100
# Operating mode
- continueaftercheck 0
+ continueaftercheck 0
diff --git a/bin/noisechisel/detection.c b/bin/noisechisel/detection.c
index 891e41a..45b89c5 100644
--- a/bin/noisechisel/detection.c
+++ b/bin/noisechisel/detection.c
@@ -498,13 +498,12 @@ detection_sn(struct noisechiselparams *p, gal_data_t
*worklab, size_t num,
int32_t *plabend, *indarr=NULL;
double ave, err, *pos, *brightness;
size_t ind, ndim=p->input->ndim, pcols=1+ndim;
+ float s, ss, *f, *ff, *fs, *sky=p->sky->array;
size_t i, j, *area, counter=0, *dsize=p->input->dsize;
- float *img=p->input->array, *f=p->input->array, *ff=f+p->input->size;
int32_t *plab = worklab->array, *dlab = s0d1D2 ? NULL : p->olabel->array;
size_t *coord=gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__,
"coord");
-
/* Sanity check. */
if(p->input->type!=GAL_TYPE_FLOAT32)
error(EXIT_FAILURE, 0, "%s: the input dataset must be float32 type, "
@@ -540,6 +539,8 @@ detection_sn(struct noisechiselparams *p, gal_data_t
*worklab, size_t num,
NULL) );
/* Go over all the pixels and get the necessary information. */
+ fs = f = p->input->array;
+ ff = f + p->input->size;
do
{
/* All this work is only necessary when we are actually on a
@@ -562,22 +563,17 @@ detection_sn(struct noisechiselparams *p, gal_data_t
*worklab, size_t num,
/* Save all the necessary values. */
++area[*plab];
- brightness[*plab] += *f;
- if( *f > 0.0f ) /* For calculatiing the approximate center, */
+ gal_dimension_index_to_coord(f-fs, ndim, dsize, coord);
+ s = sky[ gal_tile_full_id_from_coord(&p->cp.tl, coord) ];
+ ss = *f-s;
+ brightness[*plab] += ss;
+ if( ss > 0.0f ) /* For calculatiing the approximate center, */
{ /* necessary for calculating Sky and STD. */
- pos[*plab*pcols ] += *f;
- switch(ndim)
- {
- case 2:
- pos[*plab*pcols+1] += ( (f-img)/dsize[1] ) * *f;
- pos[*plab*pcols+2] += ( (f-img)%dsize[1] ) * *f;
- break;
- case 3:
- pos[*plab*pcols+1] += ( (f-img) % dsize[1]*dsize[2] ) * *f;
- pos[*plab*pcols+2] += ( (f-img)/dsize[2] ) * *f;
- pos[*plab*pcols+3] += ( (f-img)%dsize[2] ) * *f;
- break;
- }
+ pos[ *plab*pcols ] += ss;
+ pos[ *plab*pcols+1 ] += (double)coord[0] * ss;
+ pos[ *plab*pcols+2 ] += (double)coord[1] * ss;
+ if(ndim==3)
+ pos[*plab*pcols+3 ] += (double)coord[2] * ss;
}
}
@@ -632,9 +628,7 @@ detection_sn(struct noisechiselparams *p, gal_data_t
*worklab, size_t num,
for(j=0;j<ndim;++j)
coord[j]=GAL_DIMENSION_FLT_TO_INT(pos[i*pcols+j+1]/pos[i*pcols]);
- /* Calculate the Sky and Sky standard deviation on this tile. */
- ave -= ((float *)(p->sky->array))[
- gal_tile_full_id_from_coord(&p->cp.tl, coord) ];
+ /* Get the Sky standard deviation on this tile. */
err = ((float *)(p->std->array))[
gal_tile_full_id_from_coord(&p->cp.tl, coord) ];
@@ -659,6 +653,12 @@ detection_sn(struct noisechiselparams *p, gal_data_t
*worklab, size_t num,
}
}
+
+ /* A small sanity check. */
+ if( s0d1D2==0 && counter==0 )
+ error(EXIT_FAILURE, 0, "no sky pseudo-detections.");
+
+
/* If we are in Sky mode, the sizes have to be corrected */
if(s0d1D2==0)
{
diff --git a/bin/noisechisel/main.h b/bin/noisechisel/main.h
index a2590bc..c0d64fd 100644
--- a/bin/noisechisel/main.h
+++ b/bin/noisechisel/main.h
@@ -56,11 +56,11 @@ struct noisechiselparams
uint8_t rawoutput; /* Only detection & 1 elem/tile output. */
uint8_t label; /* Label detections that are connected. */
- float mirrordist; /* Maximum distance to check mode sym. */
- float modmedqdiff; /* Difference between mode and median. */
+ float meanmedqdiff; /* Difference between mode and median. */
float minskyfrac; /* Undetected area min. frac. in tile. */
float qthresh; /* Quantile threshold on convolved image. */
- float qthreshtilequant; /* Remove tiles with lower quantile. */
+ float outliersigma; /* Multiple of sigma to define outlier. */
+ double outliersclip[2]; /* Outlier Sigma-clipping params. */
size_t smoothwidth; /* Interpolation: 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 593552e..a434cbc 100644
--- a/bin/noisechisel/noisechisel.c
+++ b/bin/noisechisel/noisechisel.c
@@ -157,15 +157,17 @@ noisechisel_output(struct noisechiselparams *p)
}
- /* Write the object labels and useful information into it's header. */
+ /* Write the detected pixels and useful information into it's header. */
gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "DETSN", 0, &p->detsnthresh,
0, "Minimum S/N of true pseudo-detections", 0,
"ratio");
if(p->label)
+ gal_fits_key_list_add(&keys, GAL_TYPE_SIZE_T, "NUMLABS", 0,
+ &p->numdetections, 0, "Total number of labels "
+ "(inclusive)", 0, "counter");
+ gal_fits_key_list_reverse(&keys);
+ if(p->label)
{
- gal_fits_key_list_add(&keys, GAL_TYPE_SIZE_T, "NUMLABS", 0,
- &p->numdetections, 0, "Total number of labels "
- "(inclusive)", 0, "counter");
p->olabel->name = "DETECTIONS";
gal_fits_img_write(p->olabel, p->cp.output, keys, PROGRAM_NAME);
p->olabel->name=NULL;
diff --git a/bin/noisechisel/sky.c b/bin/noisechisel/sky.c
index 2d003c9..0e23b44 100644
--- a/bin/noisechisel/sky.c
+++ b/bin/noisechisel/sky.c
@@ -35,6 +35,8 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gnuastro/threads.h>
#include <gnuastro/statistics.h>
+#include <gnuastro-internal/tile-internal.h>
+
#include "main.h"
#include "ui.h"
@@ -56,25 +58,27 @@ sky_mean_std_undetected(void *in_prm)
struct gal_threads_params *tprm=(struct gal_threads_params *)in_prm;
struct noisechiselparams *p=(struct noisechiselparams *)tprm->params;
- double *darr, s, s2;
- int setblank, type=p->sky->type;
- size_t i, tind, numsky, dsize=2;
- gal_data_t *tile, *meanstd_d, *meanstd, *bintile;
+ int setblank, type=GAL_TYPE_FLOAT32;
+ size_t twidth=gal_type_sizeof(GAL_TYPE_FLOAT32);
+ size_t i, tind, numsky, bdsize=2, ndim=p->sky->ndim;
+ gal_data_t *tile, *fusage, *busage, *bintile, *sigmaclip;
- /* A dataset to keep the mean and STD in double type. */
- meanstd_d=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &dsize,
- NULL, 0, -1, NULL, NULL, NULL);
- darr=meanstd_d->array;
+ /* Put the temporary usage space for this thread into a data set for easy
+ processing. */
+ fusage=gal_data_alloc(NULL, type, ndim, p->maxtsize, NULL, 0,
+ p->cp.minmapsize, NULL, NULL, NULL);
+ busage=gal_data_alloc(NULL, GAL_TYPE_UINT8, ndim, p->maxtsize, NULL, 0,
+ p->cp.minmapsize, NULL, NULL, NULL);
/* An empty dataset to replicate a tile on the binary array. */
- bintile=gal_data_alloc(NULL, GAL_TYPE_UINT8, 1, &dsize,
+ bintile=gal_data_alloc(NULL, GAL_TYPE_UINT8, 1, &bdsize,
NULL, 0, -1, NULL, NULL, NULL);
+ bintile->ndim=ndim;
free(bintile->array);
free(bintile->dsize);
bintile->block=p->binary;
- bintile->ndim=p->binary->ndim;
/* Go over all the tiles given to this thread. */
@@ -86,53 +90,61 @@ sky_mean_std_undetected(void *in_prm)
tile = &p->cp.tl.tiles[tind];
/* Correct the fake binary tile's properties to be the same as this
- one, then count the number of zero valued elements in it. */
+ one, then count the number of zero valued elements in it. Note
+ that the `CHECK_BLANK' flag of `GAL_TILE_PARSE_OPERATE' is set to
+ 1. So blank values in the input array are not counted also. */
bintile->size=tile->size;
bintile->dsize=tile->dsize;
bintile->array=gal_tile_block_relative_to_other(tile, p->binary);
GAL_TILE_PARSE_OPERATE(tile, bintile, 1, 1, {if(!*o) numsky++;});
- /* Only continue, if the fraction of Sky values are less than the
+ /* Only continue, if the fraction of Sky values is less than the
requested fraction. */
setblank=0;
if( (float)(numsky)/(float)(tile->size) > p->minskyfrac)
{
- /* Calculate the mean and STD over this tile. */
- s=s2=0.0f;
- GAL_TILE_PARSE_OPERATE(tile, bintile, 1, 1, {
- if(!*o)
- {
- s += *i;
- s2 += *i * *i;
- }
- } );
- darr[0]=s/numsky;
- darr[1]=sqrt( (s2-s*s/numsky)/numsky );
+ /* Re-initialize the usage array's size information (will be
+ corrected to this tile's size by
+ `gal_data_copy_to_allocated'). */
+ busage->ndim = fusage->ndim = ndim;
+ busage->size = fusage->size = p->maxtcontig;
+ gal_data_copy_to_allocated(tile, fusage);
+ gal_data_copy_to_allocated(bintile, busage);
+
+
+ /* Set all the non-zero pixels in `busage' to NaN in `fusage'. */
+ busage->flag = fusage->flag = 0;
+ gal_blank_flag_apply(fusage, busage);
+
+
+ /* Do the sigma-clipping. */
+ sigmaclip=gal_statistics_sigma_clip(fusage, p->sigmaclip[0],
+ p->sigmaclip[1], 1, 1);
/* When there are zero-valued pixels on the edges of the dataset
(that have not been set to NaN/blank), given special
conditions, the whole zero-valued region can get a binary
value of 1 and so the Sky and its standard deviation can
become zero. So, we need ignore such tiles. */
- if(darr[1]==0.0f)
+ if( ((float *)(sigmaclip->array))[3]==0.0 )
setblank=1;
else
{
- /* Convert the mean and std into the same type as the sky and
- std arrays. */
- meanstd=gal_data_copy_to_new_type(meanstd_d, type);
-
- /* Copy the mean and STD to their respective places in the
- tile arrays. */
+ /* Copy the sigma-clipped mean and STD to their respective
+ places in the tile arrays. But first, make sure
+ `sigmaclip' has the same type as the sky and std
+ arrays. */
+ sigmaclip=gal_data_copy_to_new_type_free(sigmaclip, type);
memcpy(gal_pointer_increment(p->sky->array, tind, type),
- meanstd->array, gal_type_sizeof(type));
+ gal_pointer_increment(sigmaclip->array, 2, type),
+ twidth);
memcpy(gal_pointer_increment(p->std->array, tind, type),
- gal_pointer_increment(meanstd->array, 1, type),
- gal_type_sizeof(type));
-
- /* Clean up. */
- gal_data_free(meanstd);
+ gal_pointer_increment(sigmaclip->array, 3, type),
+ twidth);
}
+
+ /* Clean up. */
+ gal_data_free(sigmaclip);
}
else
setblank=1;
@@ -151,8 +163,9 @@ sky_mean_std_undetected(void *in_prm)
/* Clean up and wait for other threads to finish and abort. */
bintile->array=NULL;
bintile->dsize=NULL;
+ gal_data_free(fusage);
+ gal_data_free(busage);
gal_data_free(bintile);
- gal_data_free(meanstd_d);
if(tprm->b) pthread_barrier_wait(tprm->b);
return NULL;
}
@@ -180,10 +193,10 @@ sky_and_std(struct noisechiselparams *p, char *checkname)
/* Allocate space for the mean and standard deviation. */
- p->sky=gal_data_alloc(NULL, p->input->type, p->input->ndim, tl->numtiles,
- NULL, 0, cp->minmapsize, "SKY", p->input->unit, NULL);
- p->std=gal_data_alloc(NULL, p->input->type, p->input->ndim, tl->numtiles,
- NULL, 0, cp->minmapsize, "STD", p->input->unit, NULL);
+ p->sky=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, p->input->ndim, tl->numtiles,
+ NULL, 0, cp->minmapsize, NULL, p->input->unit, NULL);
+ p->std=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, p->input->ndim, tl->numtiles,
+ NULL, 0, cp->minmapsize, NULL, p->input->unit, NULL);
/* Find the Sky and its STD on proper tiles. */
@@ -191,14 +204,29 @@ sky_and_std(struct noisechiselparams *p, char *checkname)
cp->numthreads);
if(checkname)
{
+ p->sky->name="SKY";
+ p->std->name="STD";
gal_tile_full_values_write(p->sky, tl, 1, checkname, NULL,
PROGRAM_NAME);
gal_tile_full_values_write(p->std, tl, 1, checkname, NULL,
PROGRAM_NAME);
+ p->sky->name=p->std->name=NULL;
}
- /* Get the basic information about the standard deviation
- distribution. */
+
+ /* Remove the outliers. */
+ if(p->outliersigma!=0.0)
+ gal_tileinternal_no_outlier(p->sky, p->std, NULL, tl, p->outliersclip,
+ p->outliersigma, checkname);
+
+
+ /* Set the blank checked bit of the ararys to zero, most probably there
+ are tiles with too much signal or outliers. */
+ p->sky->flag &= ~GAL_DATA_FLAG_BLANK_CH;
+ p->std->flag &= ~GAL_DATA_FLAG_BLANK_CH;
+
+
+ /* Basic Sky standard deviation distribution measurements. */
tmp=gal_statistics_median(p->std, 0);
tmp=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT32);
memcpy(&p->medstd, tmp->array, sizeof p->medstd);
@@ -214,6 +242,7 @@ sky_and_std(struct noisechiselparams *p, char *checkname)
memcpy(&p->maxstd, tmp->array, sizeof p->maxstd);
gal_data_free(tmp);
+
/* In case the image is in electrons or counts per second, the standard
deviation of the noise will become smaller than unity, so we need to
correct it in the S/N calculation. So, we'll calculate the correction
diff --git a/bin/noisechisel/threshold.c b/bin/noisechisel/threshold.c
index 70ab66d..71594f0 100644
--- a/bin/noisechisel/threshold.c
+++ b/bin/noisechisel/threshold.c
@@ -37,6 +37,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gnuastro-internal/timing.h>
#include <gnuastro-internal/checkset.h>
+#include <gnuastro-internal/tile-internal.h>
#include "main.h"
@@ -392,12 +393,11 @@ qthresh_on_tile(void *in_prm)
struct qthreshparams *qprm=(struct qthreshparams *)tprm->params;
struct noisechiselparams *p=qprm->p;
- double *darr;
void *tarray=NULL;
int type=qprm->erode_th->type;
- gal_data_t *modeconv = p->wconv ? p->wconv : p->conv;
- gal_data_t *tile, *mode, *qvalue, *usage, *tblock=NULL;
+ gal_data_t *meanconv = p->wconv ? p->wconv : p->conv;
size_t i, tind, twidth=gal_type_sizeof(type), ndim=p->input->ndim;
+ gal_data_t *tile, *mean, *num, *meanquant, *qvalue, *usage, *tblock=NULL;
/* Put the temporary usage space for this thread into a data set for easy
processing. */
@@ -425,31 +425,36 @@ qthresh_on_tile(void *in_prm)
the convolved image, then copy the desired contents into the
already allocated `usage' array. */
tarray=tile->array; tblock=tile->block;
- tile->array=gal_tile_block_relative_to_other(tile, modeconv);
- tile->block=modeconv;
+ tile->array=gal_tile_block_relative_to_other(tile, meanconv);
+ tile->block=meanconv;
gal_data_copy_to_allocated(tile, usage);
- tile->array=tarray; tile->block=tblock;
-
-
- /* Find the mode on this tile, note that we have set the `inplace'
- flag to `1' to avoid extra allocation. */
- mode=gal_statistics_mode(usage, p->mirrordist, 1);
-
-
- /* Check the mode value. Note that if the mode is not accurate, then
- the contents of `darr' will be NaN and all conditions will
- fail. In such cases, the tile will be ignored. */
- darr=mode->array;
- if( fabs(darr[1]-0.5f) < p->modmedqdiff )
+ tile->array=tarray;
+ tile->block=tblock;
+
+
+ /* Find the mean's quantile on this tile, note that we have already
+ copied the tile's dataset to a newly allocated place. So we have
+ set the `inplace' flag to `1' to avoid extra allocation. */
+ mean=gal_statistics_mean(usage);
+ num=gal_statistics_number(usage);
+ mean=gal_data_copy_to_new_type_free(mean, usage->type);
+ meanquant = ( *(size_t *)(num->array)
+ ? gal_statistics_quantile_function(usage, mean, 1)
+ : NULL );
+
+ /* Only continue if the mean's quantile is close enough to the
+ median. */
+ if( meanquant
+ && fabs( *(double *)(meanquant->array)-0.5f) < p->meanmedqdiff )
{
- /* The mode was (possibly) found on the wider convolved image,
- but the qthresh values have to be found on the sharper
- convolved images. This is because the distribution becomes
- more skewed with a wider kernel, helping us find tiles with no
- data more easily. But for the quantile threshold, we want to
- use the sharper convolved image to loose less of the spatial
+ /* The mean was found on the wider convolved image, but the
+ qthresh values have to be found on the sharper convolved
+ images. This is because the distribution becomes more skewed
+ with a wider kernel, helping us find tiles with no data more
+ easily. But for the quantile threshold, we want to use the
+ sharper convolved image to loose less of the spatial
information. */
- if(modeconv!=p->conv)
+ if(meanconv!=p->conv)
{
tarray=tile->array; tblock=tile->block;
tile->array=gal_tile_block_relative_to_other(tile, p->conv);
@@ -474,7 +479,7 @@ qthresh_on_tile(void *in_prm)
gal_data_free(qvalue);
/* Same for the expansion quantile. */
- if(p->detgrowquant!=1.0f)
+ if(qprm->expand_th)
{
qvalue=gal_statistics_quantile(usage, p->detgrowquant, 1);
memcpy(gal_pointer_increment(qprm->expand_th->array, tind,
@@ -489,13 +494,15 @@ qthresh_on_tile(void *in_prm)
tind, type), type);
gal_blank_write(gal_pointer_increment(qprm->noerode_th->array,
tind, type), type);
- if(p->detgrowquant!=1.0f)
+ if(qprm->expand_th)
gal_blank_write(gal_pointer_increment(qprm->expand_th->array,
tind, type), type);
}
/* Clean up and fix the tile's pointers. */
- gal_data_free(mode);
+ gal_data_free(num);
+ gal_data_free(mean);
+ gal_data_free(meanquant);
}
/* Clean up and wait for the other threads to finish, then return. */
@@ -509,149 +516,6 @@ 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_pointer_increment(first->array, start, first->type);
- second->array=gal_pointer_increment(second->array, start,
- second->type);
- if(third)
- third->array=gal_pointer_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)
{
@@ -689,12 +553,11 @@ threshold_quantile_find_apply(struct noisechiselparams *p)
qprm.noerode_th=gal_data_alloc(NULL, p->input->type, p->input->ndim,
tl->numtiles, NULL, 0, cp->minmapsize,
NULL, p->input->unit, NULL);
- if(p->detgrowquant!=1.0f)
- qprm.expand_th=gal_data_alloc(NULL, p->input->type, p->input->ndim,
- tl->numtiles, NULL, 0, cp->minmapsize,
- NULL, p->input->unit, NULL);
- else
- qprm.expand_th=NULL;
+ qprm.expand_th = ( p->detgrowquant!=1.0f
+ ? gal_data_alloc(NULL, p->input->type, p->input->ndim,
+ tl->numtiles, NULL, 0, cp->minmapsize,
+ NULL, p->input->unit, NULL)
+ : NULL );
/* Allocate temporary space for processing in each tile. */
@@ -716,7 +579,7 @@ threshold_quantile_find_apply(struct noisechiselparams *p)
if(qprm.expand_th) qprm.expand_th->flag |= GAL_DATA_FLAG_HASBLANK;
}
qprm.noerode_th->flag |= GAL_DATA_FLAG_BLANK_CH;
- if(p->detgrowquant!=1.0f) qprm.expand_th->flag |= GAL_DATA_FLAG_BLANK_CH;
+ if(qprm.expand_th) qprm.expand_th->flag |= GAL_DATA_FLAG_BLANK_CH;
if(p->qthreshname)
{
qprm.erode_th->name="QTHRESH_ERODE";
@@ -737,11 +600,11 @@ 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);
+ /* Remove outliers if requested. */
+ if(p->outliersigma!=0.0)
+ gal_tileinternal_no_outlier(qprm.erode_th, qprm.noerode_th,
+ qprm.expand_th, &p->cp.tl, p->outliersclip,
+ p->outliersigma, p->qthreshname);
/* Check if the number of acceptable tiles is more than the minimum
@@ -751,7 +614,7 @@ threshold_quantile_find_apply(struct noisechiselparams *p)
num=gal_statistics_number(qprm.erode_th);
nval=((size_t *)(num->array))[0];
if( nval < cp->interpnumngb)
- error(EXIT_FAILURE, 0, "%zu tiles can be used for interpolation of the "
+ error(EXIT_FAILURE, 0, "%zu tile(s) can be used for interpolation of the "
"quantile threshold values over the full dataset. This is smaller "
"than the requested minimum value of %zu (value to the "
"`--interpnumngb' option).\n\n"
@@ -763,9 +626,10 @@ threshold_quantile_find_apply(struct noisechiselparams *p)
"Thus don't loosen them too much. Recall that you can see all the "
"option values to Gnuastro's programs by appending `-P' to the "
"end of your command.\n\n"
- " * Decrease `--interpnumngb' to be smaller than %zu.\n"
" * Slightly decrease `--tilesize' to have more tiles.\n"
- " * Slightly increase `--modmedqdiff' to accept more tiles.\n\n"
+ " * Slightly increase `--meanmedqdiff' to accept more tiles.\n\n"
+ " * Decrease `--outliersigma' to reject less tiles as outliers."
+ " * Decrease `--interpnumngb' to be smaller than %zu.\n"
"Try appending your command with `--checkqthresh' to see the "
"successful tiles (and get a feeling of the cause/solution. Note "
"that the output is a multi-extension FITS file).\n\n"
diff --git a/bin/noisechisel/threshold.h b/bin/noisechisel/threshold.h
index daa6dfc..48680a7 100644
--- a/bin/noisechisel/threshold.h
+++ b/bin/noisechisel/threshold.h
@@ -48,6 +48,11 @@ threshold_interp_smooth(struct noisechiselparams *p,
gal_data_t **first,
char *filename);
void
+threshold_no_outlier(struct noisechiselparams *p, gal_data_t *first,
+ gal_data_t *second, gal_data_t *third,
+ char *filename);
+
+void
threshold_quantile_find_apply(struct noisechiselparams *p);
diff --git a/bin/noisechisel/ui.c b/bin/noisechisel/ui.c
index f48ea6e..c595814 100644
--- a/bin/noisechisel/ui.c
+++ b/bin/noisechisel/ui.c
@@ -115,6 +115,7 @@ ui_initialize_options(struct noisechiselparams *p,
cp->numthreads = gal_threads_number();
cp->coptions = gal_commonopts_options;
+
/* Modify common options. */
for(i=0; !gal_options_is_last(&cp->coptions[i]); ++i)
{
diff --git a/bin/noisechisel/ui.h b/bin/noisechisel/ui.h
index c3e72b2..4c66ef8 100644
--- a/bin/noisechisel/ui.h
+++ b/bin/noisechisel/ui.h
@@ -50,7 +50,7 @@ enum program_args_groups
/* Available letters for short options:
- a b f g i j n u v x y z
+ a b f g i j n r u v x y z
A E G H J O W X Y
*/
enum option_keys_enum
@@ -60,8 +60,7 @@ enum option_keys_enum
UI_KEY_KERNEL = 'k',
UI_KEY_WIDEKERNEL = 'w',
UI_KEY_MINSKYFRAC = 'B',
- UI_KEY_MIRRORDIST = 'r',
- UI_KEY_MODMEDQDIFF = 'Q',
+ UI_KEY_MEANMEDQDIFF = 'Q',
UI_KEY_QTHRESH = 't',
UI_KEY_ERODE = 'e',
UI_KEY_OPENING = 'p',
@@ -83,6 +82,9 @@ enum option_keys_enum
UI_KEY_MINNUMFALSE,
UI_KEY_SMOOTHWIDTH,
UI_KEY_QTHRESHTILEQUANT,
+ UI_KEY_OUTLIERNUM,
+ UI_KEY_OUTLIERSIGMA,
+ UI_KEY_OUTLIERSCLIP,
UI_KEY_CHECKQTHRESH,
UI_KEY_ERODENGB,
UI_KEY_NOERODEQUANT,
diff --git a/bin/statistics/args.h b/bin/statistics/args.h
index 47beddf..e2d37ef 100644
--- a/bin/statistics/args.h
+++ b/bin/statistics/args.h
@@ -468,13 +468,40 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET
},
{
- "modmedqdiff",
- UI_KEY_MODMEDQDIFF,
+ "meanmedqdiff",
+ UI_KEY_MEANMEDQDIFF,
"FLT",
0,
"Max. mode and median quantile diff. per tile.",
UI_GROUP_SKY,
- &p->modmedqdiff,
+ &p->meanmedqdiff,
+ GAL_TYPE_FLOAT32,
+ GAL_OPTIONS_RANGE_GE_0,
+ GAL_OPTIONS_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+ {
+ "outliersclip",
+ UI_KEY_OUTLIERSCLIP,
+ "FLT,FLT",
+ 0,
+ "Sigma-clip params for qthresh outliers.",
+ UI_GROUP_SKY,
+ &p->outliersclip,
+ GAL_TYPE_STRING,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ gal_options_read_sigma_clip
+ },
+ {
+ "outliersigma",
+ UI_KEY_OUTLIERSIGMA,
+ "FLT",
+ 0,
+ "Multiple of sigma to define outliers.",
+ UI_GROUP_SKY,
+ &p->outliersigma,
GAL_TYPE_FLOAT32,
GAL_OPTIONS_RANGE_GE_0,
GAL_OPTIONS_MANDATORY,
diff --git a/bin/statistics/aststatistics.conf
b/bin/statistics/aststatistics.conf
index 07b9c92..0bf3b83 100644
--- a/bin/statistics/aststatistics.conf
+++ b/bin/statistics/aststatistics.conf
@@ -21,7 +21,9 @@
# Sky and its STD settings
khdu 1
- modmedqdiff 0.01
+ meanmedqdiff 0.005
+ outliersigma 10
+ outliersclip 3,0.2
smoothwidth 3
sclipparams 3,0.1
diff --git a/bin/statistics/main.h b/bin/statistics/main.h
index fe58f90..3f7a698 100644
--- a/bin/statistics/main.h
+++ b/bin/statistics/main.h
@@ -86,7 +86,9 @@ struct statisticsparams
char *kernelname; /* File name of kernel to convolve input. */
char *khdu; /* Kernel HDU. */
- float modmedqdiff; /* Mode and median quantile difference. */
+ float meanmedqdiff; /* Mode and median quantile difference. */
+ float outliersigma; /* Multiple of sigma to define outlier. */
+ double outliersclip[2]; /* Outlier Sigma-clipping params. */
size_t smoothwidth; /* Width of flat kernel to smooth interpd. */
uint8_t checksky; /* Save the steps for deriving the Sky. */
double sclipparams[2]; /* Muliple and parameter of sigma clipping. */
diff --git a/bin/statistics/sky.c b/bin/statistics/sky.c
index 8ab6a36..e61c74d 100644
--- a/bin/statistics/sky.c
+++ b/bin/statistics/sky.c
@@ -39,6 +39,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gnuastro-internal/timing.h>
#include <gnuastro-internal/checkset.h>
+#include <gnuastro-internal/tile-internal.h>
#include "main.h"
@@ -52,11 +53,10 @@ sky_on_thread(void *in_prm)
struct gal_threads_params *tprm=(struct gal_threads_params *)in_prm;
struct statisticsparams *p=(struct statisticsparams *)tprm->params;
- double *darr;
- int stype=p->sky_t->type;
void *tblock=NULL, *tarray=NULL;
- gal_data_t *tile, *mode, *sigmaclip;
- size_t i, tind, twidth=gal_type_sizeof(stype);
+ int stype=p->sky_t->type, itype=p->input->type;
+ gal_data_t *num, *tile, *mean, *meanquant, *sigmaclip;
+ size_t i, tind, twidth=gal_type_sizeof(p->sky_t->type);
/* Find the Sky and its standard deviation on the tiles given to this
@@ -67,25 +67,32 @@ sky_on_thread(void *in_prm)
tind = tprm->indexs[i];
tile = &p->cp.tl.tiles[tind];
-
/* If we have a convolved image, temporarily (only for finding the
- mode) change the tile's pointers so we can work on the convolved
- image for the mode. */
+ mean) change the tile's pointers so we can work on the convolved
+ image for the mean. */
if(p->kernel)
{
tarray=tile->array; tblock=tile->block;
tile->array=gal_tile_block_relative_to_other(tile, p->convolved);
tile->block=p->convolved;
}
- mode=gal_statistics_mode(tile, p->mirrordist, 1);
- if(p->kernel) { tile->array=tarray; tile->block=tblock; }
+ /* Calculate the mean's quantile. */
+ mean=gal_statistics_mean(tile);
+ num=gal_statistics_number(tile);
+ mean=gal_data_copy_to_new_type_free(mean, itype);
+ meanquant = ( *(size_t *)(num->array)
+ ? gal_statistics_quantile_function(tile, mean, 1)
+ : NULL );
+
+ /* Reset the pointers of `tile'. */
+ if(p->kernel) { tile->array=tarray; tile->block=tblock; }
- /* Check the mode value. Note that if the mode is in-accurate, then
- the values will be NaN and all conditionals will fail. So, we'll
- go onto finding values for this tile */
- darr=mode->array;
- if( fabs(darr[1]-0.5f) < p->modmedqdiff )
+ /* Check the mean quantile value. Note that if the mode is
+ in-accurate, then the values will be NaN and all conditionals will
+ fail. So, we'll go onto finding values for this tile */
+ if( meanquant
+ && fabs( *(double *)(meanquant->array)-0.5f) < p->meanmedqdiff )
{
/* Get the sigma-clipped mean and standard deviation. `inplace'
is irrelevant here because this is a tile and it will be
@@ -113,7 +120,9 @@ sky_on_thread(void *in_prm)
}
/* Clean up. */
- gal_data_free(mode);
+ gal_data_free(num);
+ gal_data_free(mean);
+ gal_data_free(meanquant);
}
@@ -179,7 +188,6 @@ sky(struct statisticsparams *p)
"SKY STD", p->input->unit, NULL);
-
/* Find the Sky and Sky standard deviation on the tiles. */
if(!cp->quiet) gettimeofday(&t1, NULL);
gal_threads_spin_off(sky_on_thread, p, tl->tottiles, cp->numthreads);
@@ -202,6 +210,13 @@ sky(struct statisticsparams *p)
}
+ /* Remove outliers if requested. */
+ if(p->outliersigma!=0.0)
+ gal_tileinternal_no_outlier(p->sky_t, p->std_t, NULL, &p->cp.tl,
+ p->outliersclip, p->outliersigma,
+ p->checkskyname);
+
+
/* Interpolate the Sky and its standard deviation. */
if(!cp->quiet) gettimeofday(&t1, NULL);
p->sky_t->next=p->std_t;
diff --git a/bin/statistics/ui.c b/bin/statistics/ui.c
index 9ba99b9..4503fb2 100644
--- a/bin/statistics/ui.c
+++ b/bin/statistics/ui.c
@@ -132,7 +132,7 @@ ui_initialize_options(struct statisticsparams *p,
p->quantmax = NAN;
p->mirror = NAN;
p->mirrordist = NAN;
- p->modmedqdiff = NAN;
+ p->meanmedqdiff = NAN;
p->sclipparams[0] = NAN;
p->sclipparams[1] = NAN;
@@ -416,17 +416,17 @@ ui_read_check_only_options(struct statisticsparams *p)
if( p->sky )
{
/* Mandatory options. */
- if ( isnan(p->modmedqdiff) || isnan(p->sclipparams[0])
- || p->cp.interpnumngb==0 || isnan(p->mirrordist) )
- error(EXIT_FAILURE, 0, "`--modmedqdiff', `--sclipparams', "
- "`--mirrordist', and `--interpnumngb' are mandatory with "
- "`--sky'");
+ if( isnan(p->meanmedqdiff) || isnan(p->sclipparams[0])
+ || p->cp.interpnumngb==0 )
+ error(EXIT_FAILURE, 0, "`--meanmedqdiff', `--sclipparams' and "
+ "`--interpnumngb' are mandatory when requesting Sky "
+ "measurement (`--sky')");
/* If mode and median distance is a reasonable value. */
- if(p->modmedqdiff>0.5)
- error(EXIT_FAILURE, 0, "%f not acceptable for `--modmedqdiff'. It "
+ if(p->meanmedqdiff>0.5)
+ error(EXIT_FAILURE, 0, "%f not acceptable for `--meanmedqdiff'. It "
"cannot take values larger than 0.5 (quantile of median)",
- p->modmedqdiff);
+ p->meanmedqdiff);
/* If a kernel name has been given, we need the HDU. */
if(p->kernelname && gal_fits_name_is_fits(p->kernelname)
diff --git a/bin/statistics/ui.h b/bin/statistics/ui.h
index 13efee9..206d1e6 100644
--- a/bin/statistics/ui.h
+++ b/bin/statistics/ui.h
@@ -92,7 +92,10 @@ enum option_keys_enum
UI_KEY_MAXBINONE,
UI_KEY_KHDU,
UI_KEY_MIRRORDIST,
- UI_KEY_MODMEDQDIFF,
+ UI_KEY_MEANMEDQDIFF,
+ UI_KEY_OUTLIERNUM,
+ UI_KEY_OUTLIERSIGMA,
+ UI_KEY_OUTLIERSCLIP,
UI_KEY_SMOOTHWIDTH,
UI_KEY_CHECKSKY,
UI_KEY_IGNOREBLANKINSKY,
diff --git a/doc/announce-acknowledge.txt b/doc/announce-acknowledge.txt
index 3e63c06..e8e934a 100644
--- a/doc/announce-acknowledge.txt
+++ b/doc/announce-acknowledge.txt
@@ -1 +1,5 @@
Alphabetically ordered list to acknowledge in the next release.
+
+Pierre-Alain Duc
+Gaspar Galaz
+Michael Stein
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 38bb4f5..6ca7fbe 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -2515,9 +2515,10 @@ $ ds9
download/hlsp_xdf_hst_wfc3ir-60mas_hudf_f160w_v1_sci.fits \
The first thing you might notice is that the regions with no data have a
value of zero in this image. The next thing might be that the dataset
actually has two ``depth''s (see @ref{Quantifying measurement limits}). The
-exposure time of the inner region is more than 4 times of the outer
+exposure time of the deep inner region is more than 4 times of the outer
parts. Fortunately the XDF survey webpage (above) contains the vertices of
-the deep flat WFC3-IR field. You can use those vertices in @ref{Crop} to
+the deep flat WFC3-IR field. With Gnuastro's Crop program@footnote{To learn
+more about the crop program see @ref{Crop}.}, you can use those vertices to
cutout this deep infra-red region from the larger image. We'll make a
directory called @file{flat-ir} and keep the flat infra-red regions in that
directory (with a `@file{xdf-}' suffix for a shorter and easier filename).
@@ -2552,12 +2553,21 @@ $ for f in f105w f160w; do
\
done
@end example
-Please open these images and inspect them with the same ds9 commands you
-used above. You will see how it is completely flat now and doesn't have
-varying depths. Another important result of this crop is that regions with
-no data now have a NaN (blank) value, not zero. Zero is a meaningful value
-and especially when using NoiseChisel, the input should have NaN values for
-pixels with no data, not zero.
+Please open these images and inspect them with the same @command{ds9}
+command you used above. You will see how it is nicely flat now and doesn't
+have varying depths. Another important result of this crop is that regions
+with no data now have a NaN (Not-a-Number, or a blank value) value, not
+zero. Zero is a number, and thus a meaningful value, especially when you
+later want to NoiseChisel@footnote{As you will see below, unlike most other
+detection algorithms, NoiseChisel detects the objects from their faintest
+parts, it doesn't start with their high signal-to-noise ratio peaks. Since
+the Sky is already subtracted in many images and noise fluctuates around
+zero, zero is commonly higher than the initial threshold applied. Therefore
+not ignoring zero-valued pixels in this image, will cause them to part of
+the detections!}. Generally, when you want to ignore some pixels in a
+dataset, and avoid higher-level ambiguities or complications, it is always
+best to give them blank values (not zero, or some other absurdly large or
+small number).
This is the deepest image we currently have of the sky. The first thing
that comes to mind may be this: ``How large is this field?''. Let's find
@@ -2962,8 +2972,9 @@ $ astwarp flat-ir/xdf-f160w.fits --project=0.001,0.0005
@end example
@noindent
-You can also combine multiple warps in one command. For example to first
-rotate the image, then scale it, run this command:
+If you need to do multiple warps, you can combine them in one call to
+Warp. For example to first rotate the image, then scale it, run this
+command:
@example
$ astwarp flat-ir/xdf-f160w.fits --rotate=20 --scale=0.25
@@ -2995,40 +3006,85 @@ To detect the signal in the image (separate interesting
pixels from noise),
we'll run NoiseChisel (@ref{NoiseChisel}):
@example
-$ astnoisechisel flat-ir/xdf-f160w.fits --output=nc-test.fits
+$ astnoisechisel flat-ir/xdf-f160w.fits
@end example
-NoiseChisel's output is a single file containing multiple extensions. You
-can get basic information about the extensions in a FITS file with
-Gnuastro's Fits program (see @ref{Fits}):
+NoiseChisel's output is a single FITS file containing multiple
+extensions. In the FITS format, each extension contains a separate dataset
+(image in this case). You can get basic information about the extensions in
+a FITS file with Gnuastro's Fits program (see @ref{Fits}):
@example
-$ astfits nc-test.fits
+$ astfits xdf-f160w_detected.fits
@end example
-From the output list, you see NoiseChisel's output contains 5 extensions
-and the first (counted as zero) is blank: it has value of @code{0} in the
-last/size column, showing that it contains no data, it just contains
-meta-data. All of Gnuastro's programs follow this convention of writing no
-data in the first HDU/extension. This allows the first extension to keep
-meta-data about all the extensions and is recommended by the FITS standard,
-see @ref{Fits}.
+From the output list, we see that NoiseChisel's output contains 5
+extensions and the first (counting from zero, with name
+@code{NOISECHISEL-CONFIG}) is empty: it has value of @code{0} in the last
+column (which shows its size). The first extension in all the outputs of
+Gnuastro's programs only contains meta-data: data about/describing the
+datasets within (all) the output's extension(s). This allows the first
+extension to keep meta-data about all the extensions and is recommended by
+the FITS standard, see @ref{Fits} for more. This generic meta-data (for the
+whole file) is very important for being able to reproduce this same result
+later.
-The name of each extension describes its contents: the first
-(@code{INPUT-NO-SKY}) is the Sky-subtracted input that you provided. The
-second (NoiseChisel's main output, @code{DETECTIONS}) has a numeric data
-type of @code{uint8} with only two possible values for all pixels: 0 for
-noise and 1 for signal. The third and fourth (called @code{SKY} and
-@code{SKY_STD}), have the Sky and its standard deviation values of the
-input on a tessellation and were calculated over the undetected regions.
+The second extension of NoiseChisel's output (numbered 1 and named
+@code{INPUT-NO-SKY}) is the Sky-subtracted input that you provided. The
+third (@code{DETECTIONS}) is NoiseChisel's main output which is a binary
+image with only two possible values for all pixels: 0 for noise and 1 for
+signal. Since it only has two values, to avoid taking too much space on
+your computer, its numeric datatype an unsigned 8-bit integer (or
+@code{uint8})@footnote{To learn more about numeric data types see
+@ref{Numeric data types}.}. The fourth and fifth (@code{SKY} and
+@code{SKY_STD}) extensions, have the Sky and its standard deviation values
+for the input on a tile grid and were calculated over the undetected
+regions (for more on the importance of the Sky value, see @ref{Sky value}).
+
+Reproducing your results later (or checking the configuration of the
+program that produced the dataset at a later time during your higher-level
+analysis) is very important in any research. Therefore, Let's first take a
+closer look at the @code{NOISECHISEL-CONFIG} extension. But first, we'll
+run NoiseChisel with @option{-P} to see the option values in a format we
+are already familiar with (to help in the comparison).
+
+@example
+$ astnoisechisel -P
+$ astfits xdf-f160w_detected.fits -h0
+@end example
+
+The first group of FITS header keywords are standard keywords (containing
+the @code{SIMPLE} and @code{BITPIX} keywords the first empty line). They
+are required by the FITS standard and must be present in any FITS
+extension. The second group contain the input file and all the options with
+their values in that run of NoiseChisel. Finally, the last group contain
+the date and version information of Gnuastro and its dependencies. The
+``versions and date'' group of keywords are present in all Gnuastro's FITS
+extension outputs, for more see @ref{Output FITS files}.
+
+Note that if a keyword name is larger than 8 characters, it is preceded by
+a @code{HIERARCH} keyword and that all keyword names are in capital
+letters. Therefore, if you want to see only one keyword's value by feeding
+the output to Grep, you should ask Grep to ignore case with its @option{-i}
+option (short name for @option{--ignore-case}). For example, below we'll
+check the value to the @option{--snminarea} option, note how we don't need
+Grep's @option{-i} option when it is fed with @command{astnoisechisel -P}
+since it is already in small-caps there. The extra white spaces in the
+first command are only to help in readability, you can ignore them when
+typing.
+
+@example
+$ astnoisechisel -P | grep snminarea
+$ astfits xdf-f160w_detected.fits -h0 | grep -i snminarea
+@end example
@cindex DS9
@cindex GNOME
@cindex SAO DS9
-To better understand NoiseChisel's output and the extensions described
-above, let's visually inspect it (here, we'll use SAO DS9). Since the file
-contains multiple related extensions, the easiest way to view all of them
-in DS9 is to open it as a ``Multi-extension data cube'' with the
+Let's continue with the extensions in NoiseChisel's output that contain a
+dataset by visually inspecting them (here, we'll use SAO DS9). Since the
+file contains multiple related extensions, the easiest way to view all of
+them in DS9 is to open the file as a ``Multi-extension data cube'' with the
@option{-mecube} option as shown below@footnote{You can configure your
graphic user interface to open DS9 in multi-extension cube mode by default
when using the GUI (double clicking on the file). If your graphic user
@@ -3037,30 +3093,33 @@ operating systems), a full description is given in
@ref{Viewing
multiextension FITS images}}.
@example
-$ ds9 -mecube nc-test.fits -zscale -zoom to fit
+$ ds9 -mecube xdf-f160w_detected.fits -zscale -zoom to fit
@end example
-The buttons and horizontal scroll bar in the ``cube'' window can be used to
-navigate between the extensions. In this mode, all DS9's settings (for
-example zoom or color-bar) will be identical between the extensions. Try
-zooming into to one part and seeing how the galaxies were detected along
-with the Sky and Sky standard deviation values. Just have in mind that
-NoiseChisel's job is @emph{only} detection (separating signal from noise),
-We'll segment this result later.
+A ``cube'' window opens along with DS9's main window. The buttons and
+horizontal scroll bar in this small new window can be used to navigate
+between the extensions. In this mode, all DS9's settings (for example zoom
+or color-bar) will be identical between the extensions. Try zooming into to
+one part and flipping through the extensions to see how the galaxies were
+detected along with the Sky and Sky standard deviation values for that
+region. Just have in mind that NoiseChisel's job is @emph{only} detection
+(separating signal from noise), We'll do segmentation on this result later
+to find the individual galaxies/peaks over the detected pixels.
-One good way to see if you have missed any signal is to mask all the
-detected pixels and inspect the noise pixels. For this, you can use
-Gnuastro's Arithmetic program (in particular its @code{where} operator, see
-@ref{Arithmetic operators}). With the command below, all detected pixels
-(in the @code{DETECTIONS} extension) will be set to NaN in the output
-(@file{nc-masked.fits}). To make the command easier to read/write, let's
-just put the file name in a shell variable (@code{img}) first. A shell
-variable's value can be retrieved by adding a @code{$} before its name.
+One good way to see if you have missed any signal (small galaxies, or the
+wings of brighter galaxies) is to mask all the detected pixels and inspect
+the noise pixels. For this, you can use Gnuastro's Arithmetic program (in
+particular its @code{where} operator, see @ref{Arithmetic operators}). With
+the command below, all detected pixels (in the @code{DETECTIONS} extension)
+will be set to NaN in the output (@file{nc-masked.fits}). To make the
+command easier to read/write, let's just put the file name in a shell
+variable (@code{img}) first. A shell variable's value can be retrieved by
+adding a @code{$} before its name.
@example
-$ img=nc-test.fits
+$ img=xdf-f160w_detected.fits
$ astarithmetic $img $img nan where -hINPUT-NO-SKY -hDETECTIONS \
- --output=nc-masked.fits
+ --output=mask-det.fits
@end example
@noindent
@@ -3070,84 +3129,232 @@ flip the detected pixel values (from 0 to 1 and
vice-versa) by adding a
@example
$ astarithmetic $img $img not nan where -hINPUT-NO-SKY -hDETECTIONS \
- --output=nc-masked.fits
-@end example
+ --output=mask-sky.fits
+@end example
+
+Looking again at the detected pixels, we see that there are thin
+connections between many of the smaller objects or extending from larger
+objects. This shows that we have dug in too deep, and that we are following
+correlated noise.
+
+Correlated noise is created when we warp datasets from individual exposures
+(that are each slightly offest compared to each other) into the same pixel
+grid, then add them to form the final result. Because it mixes nearby pixel
+values, correlated noise is a form of convolution and it smooths the
+image. In terms of the number of exposures (and thus correlated noise), the
+XDF dataset is by no means an ordinary dataset. It is the result of warping
+and adding roughly 80 separate exposures which can create strong correlated
+noise/smoothing. In common surveys the number of exposures is usually 10 or
+less.
+
+Let's tweak NoiseChisel's configuration a little to get a better result on
+this dataset. Don't forget that ``@emph{Good statistical analysis is not a
+purely routine matter, and generally calls for more than one pass through
+the computer}'' (Anscombe 1973, see @ref{Science and its tools}). A good
+scientist must have a good understanding of her tools to make a meaningful
+analysis. So don't hesitate in playing with the default configuration and
+reviewing the manual when you have a new dataset infront of you. Robust
+data analysis is an art, therefore a good scientist must first be a good
+artist.
NoiseChisel can produce ``Check images'' to help you visualize and inspect
-how each step is completed. You can see all the check images it can produce
-with this command.
+how each step is done. You can see all the check images it can produce with
+this command.
@example
$ astnoisechisel --help | grep check
@end example
-Let's see how the quantile threshold (the first step after convolution) has
-been found and applied in our previous run of NoiseChisel:
+Let's check the overall detection process to get a better feeling of what
+NoiseChisel is doing with the following command. To learn the details of
+NoiseChisel in more detail, please see
+@url{https://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa [2015]}.
+
+@example
+$ astnoisechisel flat-ir/xdf-f160w.fits --checkdetection
+@end example
+
+The check images/tables are also multi-extension FITS files. As you saw
+from the command above, when check datasets are requested, NoiseChisel
+won't go to the end. It will abort as soon as all the extensions of the
+check image are ready. Try listing the extensions of the output with
+@command{astfits} and then opening them with @command{ds9} as we done
+above. In order to understand the parameters and their biases (especially
+as you are starting to use Gnuastro, or running it a new dataset), it is
+@emph{strongly} encouraged to play with the different parameters and use
+the respective check images to see which step is affected by your changes
+and how, for example see @ref{Detecting large extended targets}.
+
+@cindex FWHM
+The @code{OPENED_AND_LABELED} extension shows the initial detection step of
+NoiseChisel. We see these thin connections between smaller points are
+already present here (a relatively early stage in the processing). Such
+connections at the lowest surface brightness limits usually occur when the
+dataset is too smoothed. Because of correlated noise, the dataset is
+already artificially smoothed, therefore futher smoothing it with the
+default kernel may be the problem. Therefore, one solution is to use a
+sharper kernel (NoiseChisel's first step in its processing).
+
+By default NoiseChisel uses a Gaussian with full-width-half-maximum (FWHM)
+of 2 pixels. We can use Gnuastro's MakeProfiles to build a kernel with FWHM
+of 1.5 pixel (truncated at 5 times the FWHM, like the default) using the
+following command. MakeProfiles is a powerful tool to build any number of
+mock profiles on one image or independently, to learn more of its features
+and capabilities, see @ref{MakeProfiles}.
+
+@example
+$ astmkprof --kernel=gaussian,1.5,5 --oversample=1
+@end example
+
+@noindent
+Please open the output @file{kernel.fits} and have a look (it is very small
+and sharp). We can now tell NoiseChisel to use this instead of the default
+kernel with the following command (we'll keep checking the detection steps)
+
+@example
+$ astnoisechisel flat-ir/xdf-f160w.fits --kernel=kernel.fits \
+ --checkdetection
+@end example
+
+Looking at the @code{OPENED_AND_LABELED} extension, we see that the thin
+connections between smaller peaks has now significantly decreased. Going
+two extensions/steps ahead (in the first @code{HOLES-FILLED}), you can see
+that during the process of finding false pseudo-detections, too many holes
+have been filled: see how the many of the brighter galaxies are connected?
+
+Try looking two extensions ahead (in @code{PSEUDOS-FOR-SN}), you can see
+that there aren't too many pseudo-detections because of all those extended
+filled holes. If you look closely, you can see the number of
+pseudo-detections in the result NoiseChisel prints (around 4000). This is
+another side-effect of correlated noise. To address it, we should slightly
+increase the pseudo-detection threhold (@option{--dthresh}, run with
+@option{-P} to see its default value):
+
+@example
+$ astnoisechisel flat-ir/xdf-f160w.fits --kernel=kernel.fits \
+ --dthresh=0.2 --checkdetection
+@end example
+
+Before visually inspecting the check image, you can see the effect of this
+change in NoiseChisel's command-line output: notice how the number of
+pseudos has increased to roughly 5500. Open the check image now and have a
+look, you can see how the pseudo-detections are distributed much more
+evenly in the image. The signal-to-noise ratio of pseudo-detections define
+NoiseChisel's reference for removing false detections, so they are very
+important to get right. Let's have a look at their signal-to-noise
+distribution with @option{--checksn}.
+
+@example
+$ astnoisechisel flat-ir/xdf-f160w.fits --kernel=kernel.fits \
+ --dthresh=0.2 --checkdetection --checksn
+@end example
+
+The output @file{xdf-f160w_detsn_sky.fits} table contains the label and
+signal-to-noise ratio of all the pseudo-detections in
+@code{PSEUDOS-FOR-SN}. You can see the table columns with the first command
+below and get a feeling for its distribution with the second command. We'll
+discuss the two Table and Statistics programs later.
@example
-$ astnoisechisel flat-ir/xdf-f160w.fits --checkqthresh
+$ asttable xdf-f160w_detsn_sky.fits
+$ aststatistics xdf-f160w_detsn_sky.fits -c2
@end example
-The check images/tables are also multi-extension FITS files. As you see,
-when check datasets are requested, NoiseChisel won't go to the end and
-abort as soon as all its extensions are ready. Try listing the extensions
-with @command{astfits} and then opening them with @command{ds9} as we done
-above. It is @emph{strongly} encouraged to play with the different
-parameters and use the respective check images to see which step is
-affected by your change.
+The correlated noise is again visible in this pseudo-detection
+signal-to-noise distribution: it is highly skewed. A small change in the
+quantile will translate into a big change in the S/N value. For example see
+the difference between the three 0.99, 0.95 and 0.90 quantiles with this
+command:
-One major factor determining NoiseChisel's the quantile threshold is
-NoiseChisel's ability to identify signal in a tile (the threshold is only
-found on those tiles with no major signal). Therefore the larger the tiles
-are, number-statistics will cause less scatter, therefore helping
-NoiseChisel find the quantile threshold. However, if the tiles are too
-large, there might not be enough over the dataset or they won't be able to
-deal with possible gradients. Let's see what the default (small) tile size
-was with the following command.
+@example
+$ aststatistics xdf-f160w_detsn_sky.fits -c2 \
+ --quantile=0.99 --quantile=0.95 --quantile=0.90
+@end example
+
+If you run NoiseChisel with @option{-P}, you'll see the default
+signal-to-noise quantile @option{--snquant} is 0.99. In effect with this
+option you specify the purity level you want (contamination by false
+detections). With the @command{aststatistics} command above, you see that a
+small number of extra false detections (impurity) in the final result
+causes a big change in completness (you can detect more lower
+signal-to-noise true detections). So let's loosen-up our desired purity
+level and then mask the detected pixels like before to see if we have
+missed anything.
@example
-$ astnoisechisel -P | grep tilesize
+$ astnoisechisel flat-ir/xdf-f160w.fits --kernel=kernel.fits \
+ --dthresh=0.2 --snquant=0.95
+$ img=xdf-f160w_detected.fits
+$ astarithmetic $img $img nan where -h1 -h2 --output=mask-det.fits
@end example
-Its a 50 by 50 box. Flip back and forth between the @code{CONVOLVED} and
-@code{QTHRESH_ERODE} extensions of the check image to get a feeling of
-which tiles succeeded (have non-blank values)@footnote{The other extensions
-don't matter much for the check here. They show the @option{--qthresh}
-values for the other thresholds (on the same tiles), followed by the
-interpolated values for all the thresholds and afterwards the smoothed
-values that are used for the next steps. The final extension
-(@code{QTHRESH-APPLIED}, shows the result of applying the erode and
-no-erode thresholds.}. Since this is a relatively large image and we don't
-have any gradients, let's increase the tile size to 100 by 100
+Overall it seems good, but if you play a little with the color-bar and look
+closer in the noise, you'll see a few very sharp, but faint, objects that
+have not been detected. This only happens for under-sampled datasets like
+HST (where the pixel size is larger than the point spread function
+FWHM). So this won't happen on ground-based images. Because of this, sharp
+and faint objects will be very small and eroded too easily during
+NoiseChisel's erosion step.
+
+To address this problem of sharp objects, we can use NoiseChisel's
+@option{--noerodequant} option. All pixels above this quantile will not be
+eroded, thus allowing us to preserve faint and sharp objects. Check its
+default value, then run NoiseChisel like below and make the mask again. You
+will see many of those sharp objects are now detected.
@example
-$ astnoisechisel flat-ir/xdf-f160w.fits --tilesize=100,100 \
- --checkqthresh
+$ astnoisechisel flat-ir/xdf-f160w.fits --kernel=kernel.fits \
+ --noerodequant=0.95 --dthresh=0.2 --snquant=0.95
@end example
-Inspecting the check image, you see that there are now only a handful of
-useful tiles in the central parts. This shows the field is too crowded, and
-we should slightly decrease the tile size for a more robust result that
-also covers more of the dataset. Let's set it to a 75 by 75 pixel box:
+This seems to be fine and we can continue with our analysis. Before finally
+running NoiseChisel, let's just see how you can have all the raw outputs of
+NoiseChisel (Detection map and Sky and Sky Standard deviation) in a highly
+compressed format for archivability. For example the Sky-subtracted input
+is a redundant dataset: you can always generate it by subtracting the Sky
+from the input image. With the commands below you can turn the default
+NoiseChisel output that is larger than 100 megabytes in this case into
+about 200 kilobytes by removing all the redundant information in it, then
+compressing it:
@example
-$ astnoisechisel flat-ir/xdf-f160w.fits --tilesize=75,75 \
- --checkqthresh
+$ astnoisechisel flat-ir/xdf-f160w.fits --oneelempertile --rawoutput
+$ gzip --best xdf-f160w_detected.fits
@end example
-The result seems reasonable now: we have a larger tile size than the
-default value, but less scatter, @emph{and} the tiles cover a sufficiently
-wide area of the dataset. So, we'll use this tile size for the next
-steps. But first, let's clean all the temporary files and make a directory
-for the NoiseChisel outputs:
+You can open @file{xdf-f160w_detected.fits.gz} directly in SAO DS9 or feed
+it to any of Gnuastro's programs without having to uncompress
+it. Higher-level programs that take NoiseChisel's output as input can also
+deal with this compressed image where the Sky and its Standard deviation
+are one pixel-per-tile.
+To avoid having to write these options on every call to NoiseChisel, we'll
+just make a configuration file in a visible @file{config} directory. Then
+we'll define the hidden @file{.gnuastro} directory (that all Gnuastro's
+programs will look into for configuration files) as a symbolic link to the
+@file{config} directory. Finally, we'll write the finalized values of the
+options into NoiseChisel's standard configuration file within that
+directory. We'll also put the kernel in a separate directory to keep the
+top directory clean of any files we later need.
+
+@example
+$ mkdir kernel config
+$ ln -s config/ .gnuastro
+$ mv kernel.fits det-kernel.fits
+$ echo "kernel kernel/det-kernel.fits" > config/astnoisechisel.conf
+$ echo "noerodequant 0.95" >> config/astnoisechisel.conf
+$ echo "dthresh 0.2" >> config/astnoisechisel.conf
+$ echo "snquant 0.95" >> config/astnoisechisel.conf
+@end example
+
+@noindent
+We are now ready to finally run NoiseChisel on the two filters and keep the
+output in a dedicated directory (@file{nc}).
@example
$ rm *.fits
$ mkdir nc
-$ astnoisechisel flat-ir/xdf-f160w.fits --tilesize=75,75 \
- --output=nc/xdf-f160w.fits
-$ astnoisechisel flat-ir/xdf-f105w.fits --tilesize=75,75 \
- --output=nc/xdf-f105w.fits
+$ astnoisechisel flat-ir/xdf-f160w.fits --output=nc/xdf-f160w.fits
+$ astnoisechisel flat-ir/xdf-f105w.fits --output=nc/xdf-f105w.fits
@end example
Before continuing with the higher-level processing of this dataset, we'll
@@ -3812,59 +4019,69 @@ $ astnoisechisel r.fits -h0
As described in @ref{NoiseChisel output}, NoiseChisel's default output is a
multi-extension FITS file. A method to view them effectively and easily is
-discussed in @ref{Viewing multiextension FITS images}.
-
-Open the output @file{r_detected.fits} file and you will immediately notice
-how NoiseChisel's default configuration is not suitable for this dataset:
-the Sky estimation has failed so terribly that the tile grid (where the Sky
-was estimated, and subtracted) is visible in the first extension (input
-dataset subtracted by the Sky value). If you look into the third and fourth
-extensions (the Sky and its standard deviation) you will see how they
-exactly map NGC 5195! This is not good! There shouldn't be any signature of
-your extended target on the Sky and its standard deviation images. After
-all, the Sky is suppose to be the average value @emph{in the absence} of
-signal, see @ref{Sky value}.
-
-The fact that signal has been detected as Sky shows that you haven't done a
-good detection. Generally, any time your target is much larger than the
-tile size and the signal is almost flat (like this case), this @emph{will}
-happen, even if it isn't dramatic enough to be seen in the first
-extension. Therefore, @strong{the best place} to check the accuracy of
-your detection is the noise extensions (third and fourth extension) of
-NoiseChisel's output.
+discussed in @ref{Viewing multiextension FITS images}. For more on tweaking
+NoiseChisel and optimizing its output for archiving or sending to
+colleagues, see the NoiseChisel part of the previous tutorial in
+@ref{General program usage tutorial}.
+
+Open the output @file{r_detected.fits} file and have a look at the
+extensions, the first extension is only meta-data and contains
+NoiseChisel's configuration parameters. The rest are the Sky-subtracted
+input, the detection map, Sky values and Sky standard deviation.
+
+Flipping through the extensions in a FITS viewer (for example SAO DS9), you
+will see that the Sky-subtracted image looks reasonable (there are no major
+artifacts due to bad Sky subtraction compared to the input). The second
+extension also seems reasonable with a large detection map that covers the
+whole of NGC5195, but also extends beyond towards the bottom of the
+image. Try going back and forth between the @code{DETECTIONS} and
+@code{SKY} extensions, you will notice that there is still significant
+signal beyond the detected pixels. You can tell that this signal belongs to
+the galaxy because the far-right side of the image is dark and the brighter
+tiles (that weren't interpolated) are surrounding the detected pixels.
+
+The fact that signal from the galaxy remains in the Sky dataset shows that
+you haven't done a good detection. Generally, any time your target is much
+larger than the tile size and the signal is almost flat (like this case),
+this @emph{will} happen. Therefore, when there are large objects in the
+dataset, @strong{the best place} to check the accuracy of your detection is
+the estimated Sky image.
When dominated by the background, noise has a symmetric
distribution. However, signal is not symmetric (we don't have negative
signal). Therefore when non-constant signal is present in a noisy dataset,
the distribution will be positively skewed. This skewness is a good measure
of how much signal we have in the distribution. The skewness can be
-accurately measured by the difference in the mode and median, for more see
-@ref{Quantifying signal in a tile}, and Appendix C
-@url{https://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa [2015]}.
-
-Skewness is only a proxy for signal when the signal has structure (varies
-per pixel). Therefore, when it is approximately constant over a whole tile,
-or sub-set of the image, the signal's effect is just to shift the symmetric
-center of the noise distribution to the positive and there won't be any
-skewness: this positive@footnote{In processed images, where the Sky value
-can be over-estimated, this constant shift can be negative.} shift that
-preserves the symmetric distribution is the Sky value. When there is a
-gradient over the dataset, different tiles will have different constant
+accurately measured by the difference in the mean and median: assuming no
+strong outliers, the more distant they are, the more skewed the dataset
+is. For more see @ref{Quantifying signal in a tile}.
+
+However, skewness is only a proxy for signal when the signal has structure
+(varies per pixel). Therefore, when it is approximately constant over a
+whole tile, or sub-set of the image, the signal's effect is just to shift
+the symmetric center of the noise distribution to the positive and there
+won't be any skewness (major difference between the mean and median): this
+positive@footnote{In processed images, where the Sky value can be
+over-estimated, this constant shift can be negative.} shift that preserves
+the symmetric distribution is the Sky value. When there is a gradient over
+the dataset, different tiles will have different constant
shifts/Sky-values, for example see Figure 11 of
@url{https://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa [2015]}.
-To get less scatter in measuring the mode and median (and thus better
-estimate the skewness), you will need a larger tile. In Gnuastro, you can
-see the option values (@option{--tilesize} in this case) by adding the
-@option{-P} option to your last command. Try it. You can clearly see that
-the default tile size is indeed much smaller than this (huge) dwarf
-galaxy. Therefore NoiseChisel was unable to identify the skewness within
-the tiles under NGC 5159. Recall that NoiseChisel only uses tiles with no
-signal/skewness to define its threshold. Because of this, the threshold has
-been over-estimated on those tiles and further exacerbated the
-non-detection of the diffuse regions. To see which tiles were used for
-estimating the quantile threshold (no skewness was measured), you can use
-NoiseChisel's @option{--checkqthresh} option:
+To get less scatter in measuring the mean and median (and thus better
+estimate the skewness), you will need a larger tile. So let's play with the
+tessellation a little to see how it affects the result. In Gnuastro, you
+can see the option values (@option{--tilesize} in this case) by adding the
+@option{-P} option to your last command. Try running NoiseChisel with
+@option{-P} to see its default tile size.
+
+You can clearly see that the default tile size is indeed much smaller than
+this (huge) galaxy and its tidal features. As a result, NoiseChisel was
+unable to identify the skewness within the tiles under the outer parts of
+M51 and NGC 5159 and the threshold has been over-estimated on those
+tiles. To see which tiles were used for estimating the quantile threshold
+(no skewness was measured), you can use NoiseChisel's
+@option{--checkqthresh} option:
@example
$ astnoisechisel r.fits -h0 --checkqthresh
@@ -3874,29 +4091,50 @@ Notice how this option doesn't allow NoiseChisel to
finish. NoiseChisel
aborted after finding the quantile thresholds. When you call any of
NoiseChisel's @option{--check*} options, by default, it will abort as soon
as all the check steps have been written in the check file (a
-multi-extension FITS file). To optimize the threshold-related settings for
-this image, we'll be playing with this tile for the majority of this
-tutorial. So let's have a closer look at it.
+multi-extension FITS file). This allows you to focus on the problem you
+wanted to check as soon as possible (you can disable this feature with the
+@option{--continueaftercheck} option).
+
+To optimize the threshold-related settings for this image, let's playing
+with this quantile threshold check image a little. Don't forget that
+``@emph{Good statistical analysis is not a purely routine matter, and
+generally calls for more than one pass through the computer}'' (Anscombe
+1973, see @ref{Science and its tools}). A good scientist must have a good
+understanding of her tools to make a meaningful analysis. So don't hesitate
+in playing with the default configuration and reviewing the manual when you
+have a new dataset infront of you. Robust data analysis is an art,
+therefore a good scientist must first be a good artist.
The first extension of @file{r_qthresh.fits} (@code{CONVOLVED}) is the
-convolved input image (where the threshold is defined and applied), for
+convolved input image where the threshold(s) is defined and applied. For
more on the effect of convolution and thresholding, see Sections 3.1.1 and
3.1.2 of @url{https://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa
[2015]}. The second extension (@code{QTHRESH_ERODE}) has a blank value for
all the pixels of any tile that was identified as having significant
-signal. Playing a little with the dynamic range of this extension, you
-clearly see how the non-blank tiles around NGC 5195 have a gradient. You do
-not want this behavior. The ultimate purpose of the next few trials will be
-to remove the gradient from the non-blank tiles.
-
-The next two extensions (@code{QTHRESH_NOERODE} and @code{QTHRESH_EXPAND})
-are for later steps in NoiseChisel. The same tiles are masked, but those
-with a value, have a different value compared to @code{QTHRESH_ERODE}. In
-the subsequent three extensions, you can see how the blank tiles are
-filled/interpolated. The subsequent three extensions show the smoothed tile
-values. Finally in the last extension (@code{QTHRESH-APPLIED}), you can see
-the effect of applying @code{QTHRESH_ERODE} on @code{CONVOLVED} (pixels
-with a value of 0 were below the threshold).
+signal. The next two extensions (@code{QTHRESH_NOERODE} and
+@code{QTHRESH_EXPAND}) are the other two quantile thresholds that are
+necessary in NoiseChisel's later steps. Every step in this file is repeated
+on the three thresholds.
+
+Play a little with the color bar of the @code{QTHRESH_ERODE} extension, you
+clearly see how the non-blank tiles around NGC 5195 have a gradient. As one
+line of attack against discarding too much signal below the threshold,
+NoiseChisel rejects outlier tiles. Go forward by three extensions to
+@code{VALUE1_NO_OUTLIER} and you will see that many of the tiles over the
+galaxy have been removed in this step. For more on the outlier rejection
+algorithm, see the latter half of @ref{Quantifying signal in a tile}.
+
+However, the default outlier rejection parameters weren't enough, and when
+you play with the color-bar, you see that the faintest parts of the galaxy
+outskirts still remain. Therefore have two strategies for approaching this
+problem: 1) Increase the tile size to get more accurate measurements of
+skewness. 2) Strengthen the outlier rejection parameters to discard more of
+the tiles with signal. Fortunately in this image we have a sufficently
+large region on the right of the image that the galaxy doesn't extend
+to. So we can use the more robust first solution. In situations where this
+doesn't happen (for example if the field of view in this image was shifted
+to have more of M51 and less sky) you are limited to a combination of the
+two solutions or just to the second solution.
@cartouche
@noindent
@@ -3910,9 +4148,7 @@ directly feed the convolved image and avoid convolution.
For more on
@option{--convolved}, see @ref{NoiseChisel input}.
@end cartouche
-Fortunately this image is large and has a nice and clean region also
-(filled with very small/distant stars and galaxies). So our first solution
-is to increase the tile size. To identify the skewness caused by NGC 5195
+To identify the skewness caused by the flat NGC 5195 and M51 tidal features
on the tiles under it, we thus have to choose a tile size that is larger
than the gradient of the signal. Let's try a 100 by 100 tile size:
@@ -3920,86 +4156,30 @@ than the gradient of the signal. Let's try a 100 by 100
tile size:
$ astnoisechisel r.fits -h0 --tilesize=100,100 --checkqthresh
@end example
-You can clearly see the effect of this increased tile size: they are much
-larger. Nevertheless the higher valued tiles are still systematically
-surrounding NGC 5195. As a result, when flipping through the interpolated
-and smoothed values, you can still see the signature of the galaxy, and the
-ugly tile signatures are still present in @code{QTHRESH-APPLIED}. So let's
-increase the tile size even further (check the result of the first before
-going to the second):
-
-@example
-$ astnoisechisel r.fits -h0 --tilesize=150,150 --checkqthresh
-$ astnoisechisel r.fits -h0 --tilesize=200,200 --checkqthresh
-@end example
-
-The number of tiles with a gradient does indeed decrease with a larger tile
-size, but you still see a gradient in the raw values. The tile signatures
-in the thresholded image are also still present. These are not a good
-sign. So, let's keep the 200 by 200 tile size and start playing with the
-other constraint that we have: the acceptable distance (in quantile),
-between the mode and median.
-
-The tile size is now very large (16 times the area of the default tile
-size). We thus have much less scatter in the estimation of the mode and
-median and we can afford to decrease the acceptable distance between the
-two. The acceptable distance can be set through the @option{--modmedqdiff}
-option (read as ``mode-median quantile difference''). Before trying the
-next command, run the previous command with a @option{-P} to see what value
-it originally had.
-
-@example
-$ astnoisechisel r.fits -h0 --tilesize=200,200 --modmedqdiff=0.005 \
- --checkqthresh
-@end example
-
-But this command doesn't finish like the previous ones. A
-@file{r_qthresh.fits} file was created, but instead of saying that the
-quantile thresholds have been applied, a long error message is
-printed. Please read the error message before continuing to read here.
-
-The error message fully describes the problem and even proposes
-solutions. As suggested there, the ideal solution would be to use SDSS
-images outside of this field and add them to this one to have a larger
-input image (with a larger area outside the diffuse regions). But we don't
-always have this luxury, so let's keep using to this image for the rest of
-this tutorial.
-
-First, open @file{r_qthresh.fits} and have a look at the successful
-tiles. Unlike the previous @option{--checkqthresh} outputs, this one only
-has four extensions: as the error message explains, the interpolation (to
-give values to blank tiles) has not been done. Therefore its check results
-aren't present.
-
-At the start of the error message, NoiseChisel tells us how many tiles
-passed the test for having no significant signal: six. Looking closely at
-the dataset, we see that outside NGC 5195, there is no background gradient
-(the background is a fixed value). Our tile sizes are also very large (and
-thus don't have much scatter). So a good way to loosen up the parameters
-can be to simply decrease the number of neighboring tiles needed for
-interpolation, with @option{--interpnumngb} (read as ``interpolation number
-of neighbors'').
-
-@example
-$ astnoisechisel r.fits -h0 --tilesize=200,200 --modmedqdiff=0.005 \
- --interpnumngb=6 --checkqthresh
-@end example
-
-There is no longer any significant gradient in the usable tiles and no
-signature of NGC 5195 exists in the interpolated and smoothed values. But
-before finishing the quantile threshold, let's have a closer look at the
-thresholded image in @code{QTHRESH-APPLIED}. Slide the dynamic range in
-your FITS viewer so 0 valued pixels are black and all non-zero pixels are
-white. You will see that the black holes are not evenly distributed. Those
-that follow the tail of NGC 5195 are systematically smaller than those in
-the far-right of the image. This suggests that we can decrease the quantile
-threshold (@code{--qthresh}) even further: there is still signal down
-there!
+You can clearly see the effect of this increased tile size: the tiles are
+much larger (@mymath{\times4} in area) and when you look into
+@code{VALUE1_NO_OUTLIER}, you see that almost all the tiles under the
+galaxy have been discarded and we are only left with tiles in the
+right-most part of the image. The next group of extensions (those ending
+with @code{_INTERP}), give a value to all blank tiles based on the nearest
+tiles with a measurement. The following group of extensions (ending with
+@code{_SMOOTH}) have smoothed the interpolated image to avoid sharp cuts on
+tile edges.
+
+Inspecting @code{THRESH1_SMOOTH}, you can see that there is no longer any
+significant gradient and no major signature of NGC 5195 exists. But before
+finishing the quantile threshold, let's have a closer look at the final
+extension (@code{QTHRESH-APPLIED}) which is thresholded image. Slide the
+dynamic range in your FITS viewer so 0 valued pixels are black and all
+non-zero pixels are white. You will see that the black holes are not evenly
+distributed. Those that follow the tail of NGC 5195 are systematically
+smaller than those in the far-right of the image. This suggests that we can
+decrease the quantile threshold (@code{--qthresh}) even further: there is
+still signal down there!
@example
$ rm r_qthresh.fits
-$ astnoisechisel r.fits -h0 --tilesize=200,200 --modmedqdiff=0.005 \
- --interpnumngb=6 --qthresh=0.2
+$ astnoisechisel r.fits -h0 --tilesize=100,100 --qthresh=0.2
@end example
Since the quantile threshold of the previous command was satisfactory, we
@@ -4007,54 +4187,75 @@ finally removed @option{--checkqthresh} to let
NoiseChisel proceed until
completion. Looking at the @code{DETECTIONS} extension of NoiseChisel's
output, we see the right-ward edges in particular have many holes that are
fully surrounded by signal and the signal stretches out in the noise very
-thinly. This suggests that there is still signal that can be detected. So
-we'll decrease the growth quantile (for larger/deeper growth into the
+thinly. This suggests that there is still signal that can be detected. You
+can confirm this guess by looking at the @code{SKY} extension to see that
+indeed, we still have traces of the galaxy outskirts there. Therefore, we
+should dig deeper into the noise.
+
+Let's decrease the growth quantile (for larger/deeper growth into the
noise, with @option{--detgrowquant}) and increase the size of holes that
can be filled (if they are fully surrounded by signal, with
-@option{--detgrowmaxholesize}). Since we are done with our detection, to
-facilitate later steps, we'll also add the @option{--label} option so the
-connected regions get different labels.
+@option{--detgrowmaxholesize}).
@example
-$ astnoisechisel r.fits -h0 --tilesize=200,200 --modmedqdiff=0.005 \
- --interpnumngb=6 --qthresh=0.2 --detgrowquant=0.6 \
- --detgrowmaxholesize=10000 --label
+$ astnoisechisel r.fits -h0 --tilesize=100,100 --qthresh=0.2 \
+ --detgrowquant=0.65 --detgrowmaxholesize=10000
@end example
Looking into the output, we now clearly see that the tidal features of M51
and NGC 5195 are detected nicely in the same direction as expected (towards
the bottom right side of the image). However, as discussed above, the best
measure of good detection is the noise, not the detections themselves. So
-let's look at the Sky and its Standard deviation. The Sky standard
-deviation no longer has any footprint of NGC 5195. But the Sky still has a
-very faint shadow of the dwarf galaxy (the values on the left are larger
-than on the right). However, this gradient in the Sky (output of first
-command below) is much less (by @mymath{\sim20} times) than the standard
-deviation (output of second command). So we can stop playing with
-NoiseChisel here, and leave the configuration for a more accurate detection
-to you.
+let's look at the Sky and its Standard deviation. The Sky still has a very
+faint shadow of the galaxy outskirts (the values on the left are very
+slightly larger than those on the right).
+
+Let's calculate this gradient as a function of noise. First we'll collapse
+the image along the second (vertical) FITS dimension to have a 1D
+array. Then we'll calculate the top and bottom values of this array and
+define that as the gradient. Then we'll estimate the mean standard
+deviation over the image and divide it by the first value. The first two
+commands are just for demonstration of the collapsed dataset:
@example
-$ aststatistics r_detected.fits -hSKY --maximum --minimum \
- | awk '@{print $1-$2@}'
-$ aststatistics r_detected.fits -hSKY_STD --mean
+$ astarithmetic r_detected.fits 2 collapse-mean -hSKY -ocollapsed.fits
+$ asttable collapsed.fits
+$ skydiff=$(astarithmetic r_detected.fits 2 collapse-mean set-i \
+ i maxvalue i minvalue - -hSKY -q)
+$ echo $skydiff
+$ std=$(aststatistics r_detected.fits -hSKY_STD --mean)
+$ echo $std
+$ echo "$std $skydiff" | awk '@{print $1/$2@}'
@end example
+This gradient in the Sky (output of first command below) is much less (by
+more than 20 times) than the standard deviation (final extension). So we
+can stop configuring NoiseChisel at this point in the tutorial. We leave
+further configuration for a more accurate detection to you as an exercise.
+
Let's see how deeply/successfully we carved out M51 and NGC 5195's tail
from the noise. For this measurement, we'll need to estimate the average
flux on the outer edges of the detection. Fortunately all this can be done
with a few simple commands (and no higher-level language mini-environments)
-using @ref{Arithmetic} and @ref{MakeCatalog}.
+using @ref{Arithmetic} and @ref{MakeCatalog}. First, let's give a separate
+label to all the connected pixels of NoiseChisel's detection map:
+
+@example
+$ astarithmetic r_detected.fits 2 connected-components -hDETECTIONS \
+ -olabeled.fits
+@end example
-The M51 group detection is by far the largest detection in this image. We
-can thus easily find the ID/label that corresponds to it. We'll first run
-MakeCatalog to find the area of all the detections, then we'll use AWK to
-find the ID of the largest object and keep it as a shell variable
-(@code{id}):
+Ofcourse, you can find the the label of the main galaxy visually, but to
+have a little more fun, lets do this automatically. The M51 group detection
+is by far the largest detection in this image. We can thus easily find the
+ID/label that corresponds to it. We'll first run MakeCatalog to find the
+area of all the detections, then we'll use AWK to find the ID of the
+largest object and keep it as a shell variable (@code{id}):
@example
-$ astmkcatalog r_detected.fits --ids --geoarea -hDETECTIONS -ocat.txt
+$ astmkcatalog labeled.fits --ids --geoarea -h1 -ocat.txt
$ id=$(awk '!/^#/@{if($2>max) @{id=$1; max=$2@}@} END@{print id@}' cat.txt)
+$ echo $id
@end example
To separate the outer edges of the detections, we'll need to ``erode'' the
@@ -4072,8 +4273,8 @@ value of 0 in @file{erode.fits}. We'll keep the output in
@file{boundary.fits}.
@example
-$ astarithmetic r_detected.fits $id eq eroded.fits 0 eq and \
- -hDETECTIONS -h1 -oboundary.fits
+$ astarithmetic labeled.fits $id eq eroded.fits 0 eq and \
+ -g1 -oboundary.fits
@end example
Open the image and have a look. You'll see that the detected edge of the
@@ -4086,7 +4287,7 @@ $ astarithmetic r.fits boundary.fits nan where
-ob-masked.fits -h0
@end example
To quantify how deep we have detected the low-surface brightness regions,
-we'll use the command below. In short it just divides all the 1-valued
+we'll use the command below. In short it just divides all the non-zero
pixels of @file{boundary.fits} in the Sky subtracted input (first extension
of NoiseChisel's output) by the pixel standard deviation of the same
pixel. This will give us a signal-to-noise ratio image. The mean value of
@@ -4111,12 +4312,53 @@ $ astarithmetic r_detected.fits boundary.fits not nan
where \
r_detected.fits / \
meanvalue \
-hINPUT-NO-SKY -h1 -hSKY_STD --quiet
---> 0.0511864
@end example
@noindent
The outer wings where therefore non-parametrically detected until
-@mymath{\rm{S/N}\approx0.05}.
+@mymath{\rm{S/N}\approx0.05}!
+
+This is very good, but how much is this in units of standard magnitudes per
+arcseconds squared? To find that, we'll first need to calcaulate how many
+pixels of this image are in one arcsecond-squared. Fortunately the WCS
+headers of Gnuastro's output FITS files (and in particular @code{CDELT})
+give us this information.
+
+@example
+$ n=$(astfits r_detected.fits -h1 \
+ | awk '/CDELT1/ @{p=1/($3*3600); print p*p@}')
+$ echo $n
+@end example
+
+@noindent
+Now, let's calculate the average sky-subtracted flux in the border region
+per pixel.
+
+@example
+$ f=$(astarithmetic r_detected.fits boundary.fits not nan where set-i \
+ i sumvalue i numvalue / -q -hINPUT-NO-SKY)
+$ echo $f
+@end example
+
+@noindent
+We can just multiply the two to get the average flux on this border in one
+arcsecond squared. We also have the r-band SDSS zeropoint
+magnitude@footnote{From
+@url{http://classic.sdss.org/dr7/algorithms/fluxcal.html}} to be
+24.80. Therefore we can get the surface brightness limit (in magnitudes per
+arcsecond squared) using the following command. Just note that @code{log}
+in AWK is in base-2 (not 10) and it doesn't have a @code{log10}
+operator. So we'll do an extra division by @code{log(10)} to correct for
+this.
+
+@example
+$ z=24.80
+$ echo "$n $f $z" | awk '@{print -2.5*log($1*$2)/log(10)+$3@}'
+--> 30.0646
+@end example
+
+This shows that on a single-exposure SDSS image, we have reached a surface
+brightness limit of roughly 30 magnitudes per arcseconds squared!
In interpreting this value, you should just have in mind that NoiseChisel
works based on the contiguity of signal in the pixels. Therefore the larger
@@ -12044,18 +12286,26 @@ $ astarithmetic img1.fits img2.fits img3.fits median
\
-h0 -h1 -h2 --out=median.fits
@end example
-If the output is an image, and the @option{--output} option is not given,
-automatic output will use the name of the first FITS image encountered to
-generate an output file name, see @ref{Automatic output}. Also, output WCS
-information will be taken from the first input image encountered. When the
-output is a single number, that number will be printed in the standard
-output and no output file will be created. Arithmetic's notation for giving
-operands to operators is described in @ref{Reverse polish notation}. To
-ignore certain pixels, set them as blank, see @ref{Blank pixels}, for
-example with the @command{where} operator (see @ref{Arithmetic
-operators}). See @ref{Common options} for a review of the options in all
-Gnuastro programs. Arithmetic just redefines the @option{--hdu} option as
-explained below:
+If the output is not a single number, and the @option{--output} option is
+not given, automatic output will use the name of the first FITS image
+encountered to generate an output file name, see @ref{Automatic
+output}. Also, output WCS information will be taken from the first input
+image encountered. When the output is a single number, that number will be
+printed in the standard output and no output file will be
+created. Arithmetic's notation for giving operands to operators is
+described in @ref{Reverse polish notation}. To ignore certain pixels, set
+them as blank, see @ref{Blank pixels}, for example with the @command{where}
+operator (see @ref{Arithmetic operators}). See @ref{Common options} for a
+review of the options in all Gnuastro programs. Arithmetic just redefines
+the @option{--hdu} option as explained below.
+
+Through operators like those starting with @code{collapse-}, the
+dimensionality of the inputs may not be the same as the outputs. By
+default, when the output is 1D, Arithmetic will write it as a table, not an
+image/array. The format of the output table (plain text or FITS ASCII or
+binary) can be set with the @option{--tableformat} option, see @ref{Input
+output options}). You can disable this feature (write 1D arrays as FITS
+images/arrays) with the @option{--onedasimage} option.
@table @option
@@ -12092,6 +12342,12 @@ option is very convenient when you have many input
files and the dataset of
interest is in the same HDU of all the files. When this option is called,
any values given to the @option{--hdu} option (explained above) are ignored
and will not be used.
+
+@item -O
+@itemx --onedasimage
+When final dataset to write as output only has one dimension, write it as a
+FITS image/array. By default, if the output is 1D, it will be written as a
+table, see above.
@end table
Arithmetic accepts two kinds of input: images and numbers. Images are
@@ -14420,17 +14676,19 @@ One of the most important aspects of a dataset is its
reference value: the
value of the dataset where there is no signal. Without knowing, and thus
removing the effect of, this value it is impossible to compare the derived
results of many high-level analyses over the dataset with other datasets
-(in the attempt to associate our results with the ``real'' world). In
-astronomy, this reference value is known as the ``Sky'' value: the value
-where there is no signal from objects (for example galaxies, stars, planets
-or comets). Depending on the dataset, the Sky value maybe a fixed value
-over the whole dataset, or it may vary based on location. For an example of
-the latter case, see Figure 11 in @url{https://arxiv.org/abs/1505.01664,
-Akhlaghi and Ichikawa (2015)}.
+(in the attempt to associate our results with the ``real'' world).
+
+In astronomy, this reference value is known as the ``Sky'' value: the value
+that noise fluctuates around: where there is no signal from detectable
+objects or artifacts (for example galaxies, stars, planets or comets, star
+spikes or internal optical ghost). Depending on the dataset, the Sky value
+maybe a fixed value over the whole dataset, or it may vary based on
+location. For an example of the latter case, see Figure 11 in
+@url{https://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa (2015)}.
Because of the significance of the Sky value in astronomical data analysis,
we have devoted this subsection to it for a thorough review. We start with
-a thorough discussion on its definition (@ref{Sky value definition}). In
+a thorough discussion on its definition (@ref{Sky value definition}). In
the astronomical literature, researchers use a variety of methods to
estimate the Sky value, so in @ref{Sky value misconceptions}) we review
those and discuss their biases. From the definition of the Sky value, the
@@ -14456,28 +14714,29 @@ This analysis is taken from
@url{https://arxiv.org/abs/1505.01664, Akhlaghi
and Ichikawa (2015)}. Let's assume that all instrument defects -- bias,
dark and flat -- have been corrected and the brightness (see @ref{Flux
Brightness and magnitude}) of a detected object, @mymath{O}, is
-desired. The sources of flux on pixel @mymath{i}@footnote{For this analysis
-the dimension of the data (image) is irrelevant. So if the data is an image
+desired. The sources of flux on pixel@footnote{For this analysis the
+dimension of the data (image) is irrelevant. So if the data is an image
(2D) with width of @mymath{w} pixels, then a pixel located on column
@mymath{x} and row @mymath{y} (where all counting starts from zero and (0,
0) is located on the bottom left corner of the image), would have an index:
-@mymath{i=x+y\times{}w}.} of the image can be written as follows:
+@mymath{i=x+y\times{}w}.} @mymath{i} of the image can be written as
+follows:
@itemize
@item
-Contribution from the target object, (@mymath{O_i}).
+Contribution from the target object (@mymath{O_i}).
@item
-Contribution from other detected objects, (@mymath{D_i}).
+Contribution from other detected objects (@mymath{D_i}).
@item
Undetected objects or the fainter undetected regions of bright
-objects, (@mymath{U_i}).
+objects (@mymath{U_i}).
@item
@cindex Cosmic rays
-A cosmic ray, (@mymath{C_i}).
+A cosmic ray (@mymath{C_i}).
@item
@cindex Background flux
The background flux, which is defined to be the count if none of the
-others exists on that pixel, (@mymath{B_i}).
+others exists on that pixel (@mymath{B_i}).
@end itemize
@noindent
The total flux in this pixel (@mymath{T_i}) can thus be written as:
@@ -14487,9 +14746,9 @@ The total flux in this pixel (@mymath{T_i}) can thus be
written as:
@cindex Cosmic ray removal
@noindent
By definition, @mymath{D_i} is detected and it can be assumed that it is
-correctly estimated (deblended) and subtracted, thus @mymath{D_i=0}. There
-are also methods to detect and remove cosmic rays, for example the method
-described in van Dokkum (2001)@footnote{van Dokkum,
+correctly estimated (deblended) and subtracted, we can thus set
+@mymath{D_i=0}. There are also methods to detect and remove cosmic rays,
+for example the method described in van Dokkum (2001)@footnote{van Dokkum,
P. G. (2001). Publications of the Astronomical Society of the Pacific.
113, 1420.}, or by comparing multiple exposures. This allows us to set
@mymath{C_i=0}. Note that in practice, @mymath{D_i} and @mymath{U_i} are
@@ -14499,37 +14758,39 @@ algorithm is perfect. With these limitations in mind,
the observed Sky
value for this pixel (@mymath{S_i}) can be defined as
@cindex Sky value
-@dispmath{S_i=B_i+U_i.}
+@dispmath{S_i\equiv{}B_i+U_i.}
@noindent
Therefore, as the detection process (algorithm and input parameters)
-becomes more accurate, or @mymath{U_i\to0}, the sky value will tend to
-the background value or @mymath{S_i\to B_i}. Therefore, while
+becomes more accurate, or @mymath{U_i\to0}, the Sky value will tend to the
+background value or @mymath{S_i\to B_i}. Hence, we see that while
@mymath{B_i} is an inherent property of the data (pixel in an image),
-@mymath{S_i} depends on the detection process. Over a group of pixels,
-for example in an image or part of an image, this equation translates
-to the average of undetected pixels. With this definition of sky, the
-object flux in the data can be calculated with
+@mymath{S_i} depends on the detection process. Over a group of pixels, for
+example in an image or part of an image, this equation translates to the
+average of undetected pixels (Sky@mymath{=\sum{S_i}}). With this
+definition of Sky, the object flux in the data can be calculated, per
+pixel, with
@dispmath{ T_{i}=S_{i}+O_{i} \quad\rightarrow\quad
O_{i}=T_{i}-S_{i}.}
-@noindent
@cindex photo-electrons
-Hence, the more accurately @mymath{S_i} is measured, the more accurately
-the brightness (sum of pixel values) of the target object can be measured
-(photometry). Any under-(over-)estimation in the sky will directly
-translate to an over-(under-)estimation of the measured object's
-brightness. In the fainter outskirts of an object a very small fraction of
-the photo-electrons in the pixels actually belong to objects (see Figure 1b
-in @url{https://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa
-(2015)}). Therefore even a small over estimation of the sky value will
+In the fainter outskirts of an object, a very small fraction of the
+photo-electrons in a pixel actually belongs to objects, the rest is caused
+by random factors (noise), see Figure 1b in
+@url{https://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa
+(2015)}. Therefore even a small over estimation of the Sky value will
result in the loss of a very large portion of most galaxies. Besides the
lost area/brightness, this will also cause an over-estimation of the Sky
value and thus even more under-estimation of the object's brightness. It is
thus very important to detect the diffuse flux of a target, even if they
are not your primary target.
+In summary, the more accurately the Sky is measured, the more accurately
+the brightness (sum of pixel values) of the target object can be measured
+(photometry). Any under/over-estimation in the Sky will directly translate
+to an over/under-estimation of the measured object's brightness.
+
@cartouche
@noindent
The @strong{Sky value} is only correctly found when all the detected
@@ -14550,15 +14811,16 @@ signal-based detection is a detection process that
relies heavily on
assumptions about the to-be-detected objects. This method was the most
heavily used technique prior to the introduction of NoiseChisel in that
paper.} use the sky value as a reference to define the detection
-threshold. So these old techniques had to rely on approximations based on
-other assumptions about the data. A review of those other techniques can be
-seen in Appendix A of Akhlaghi and Ichikawa (2015)@footnote{Akhlaghi M.,
-Ichikawa. T. (2015). Astrophysical Journal Supplement Series.}. Since they
-were extensively used in astronomical data analysis for several decades,
-such approximations have given rise to a lot of misconceptions, ambiguities
-and disagreements about the sky value and how to measure it. As a summary,
-the major methods used until now were an approximation of the mode of the
-image pixel distribution and @mymath{\sigma}-clipping.
+threshold. These older techniques therefore had to rely on approximations
+based on other assumptions about the data. A review of those other
+techniques can be seen in Appendix A of
+@url{https://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa (2015)}.
+
+These methods were extensively used in astronomical data analysis for
+several decades, therefore they have given rise to a lot of misconceptions,
+ambiguities and disagreements about the sky value and how to measure it. As
+a summary, the major methods used until now were an approximation of the
+mode of the image pixel distribution and @mymath{\sigma}-clipping.
@itemize
@cindex Histogram
@@ -14571,19 +14833,19 @@ assume (or find) a certain probability density
function (PDF) or use the
histogram. But astronomical datasets can have any distribution, making it
almost impossible to define a generic function. Also, histogram-based
results are very inaccurate (there is a large dispersion) and it depends on
-the histogram bin-widths.
+the histogram bin-widths. Generally, the mode of a distribution also shifts
+as signal is added. Therefore, even if it is accurately measured, the mode
+is a biased measure for the Sky value.
@cindex sigma-clipping
@item
-Another approach was to iteratively clip the brightest pixels in the
-image (which is known as @mymath{\sigma}-clipping, since the reference
-was found from the image mean and its standard deviation or
-@mymath{\sigma}). See @ref{Sigma clipping} for a complete
-explanation. The problem with @mymath{\sigma}-clipping was that real
-astronomical objects have diffuse and faint wings that penetrate
-deeply into the noise. So only removing their brightest parts is
-completely useless in removing the systematic bias an object's fainter
-parts cause in the sky value.
+Another approach was to iteratively clip the brightest pixels in the image
+(which is known as @mymath{\sigma}-clipping). See @ref{Sigma clipping} for
+a complete explanation. @mymath{\sigma}-clipping is useful when there are
+clear outliers (an object with a sharp edge in an image for
+example). However, real astronomical objects have diffuse and faint wings
+that penetrate deeply into the noise, see Figure 1 in
+@url{https://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa (2015)}.
@end itemize
As discussed in @ref{Sky value}, the sky value can only be correctly
@@ -14600,76 +14862,78 @@ are ultimately poor approximations.
@cindex Noise
@cindex Signal
@cindex Gaussian distribution
-Put simply, noise can characterized with a certain spread about a
-characteristic value. In the Gaussian distribution (most commonly used to
-model noise) the spread is defined by the standard deviation about the
-characteristic mean. Before continuing let's clarify some definitions
-first: @emph{Data} is defined as the combination of signal and noise (so a
-noisy image is one @emph{data}-set). @emph{Signal} is defined as the mean
-of the noise on each element (after sky subtraction, see @ref{Sky value
-definition}).
-
-Let's assume that the @emph{background} (see @ref{Sky value definition}) is
-subtracted and is zero. When a data set doesn't have any signal (only
-noise), the mean, median and mode of the distribution are equal within
-statistical errors and approximately equal to the background value. Signal
-always has a positive value and will never become negative, see Figure 1 in
+Put simply, noise can be characterized with a certain spread about the
+measured value. In the Gaussian distribution (most commonly used to model
+noise) the spread is defined by the standard deviation about the
+characteristic mean.
+
+Let's start by clarifying some definitions first: @emph{Data} is defined as
+the combination of signal and noise (so a noisy image is one
+@emph{data}set). @emph{Signal} is defined as the mean of the noise on each
+element. We'll also assume that the @emph{background} (see @ref{Sky value
+definition}) is subtracted and is zero.
+
+When a data set doesn't have any signal (only noise), the mean, median and
+mode of the distribution are equal within statistical errors and
+approximately equal to the background value. Signal always has a positive
+value and will never become negative, see Figure 1 in
@url{https://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa
-(2015)}. Therefore, as more signal is added to the raw noise, the mean,
-median and mode of the dataset (which has both signal and noise) shift to
-the positive. The mean's shift is the largest. The median shifts less,
-since it is defined based on an ordered distribution and so is not affected
-by a small number of outliers. The distribution's mode shifts the least to
-the positive.
+(2015)}. Therefore, as more signal is added, the mean, median and mode of
+the dataset shift to the positive. The mean's shift is the largest. The
+median shifts less, since it is defined based on an ordered distribution
+and so is not affected by a small number of outliers. The distribution's
+mode shifts the least to the positive.
@cindex Mode
@cindex Median
Inverting the argument above gives us a robust method to quantify the
-significance of signal in a dataset. Namely, when the mode and median of a
-distribution are approximately equal, we can argue that there is no
-significant signal. To allow for gradients (which are commonly present in
-ground-based images), we can consider the image to be made of a grid of
-tiles (see @ref{Tessellation}@footnote{The options to customize the
-tessellation are discussed in @ref{Processing options}.}). Hence, from the
-difference of the mode and median on each tile, we can `detect' the
-significance of signal in it. The median of a distribution is defined to be
-the value of the distribution's middle point after sorting (or 0.5
-quantile). Thus, to estimate the presence of signal, we'll compare with the
-quantile of the mode with 0.5, if the difference is larger than the value
-given to the @option{--modmedqdiff} option, this tile will be ignored. You
-can read this option as ``mode-median-quantile-diff''.
-
-This method to use the input's skewness is possible because of a new
-algorithm to find the mode of a distribution that was defined in Appendix C
-of Akhlaghi and Ichikawa (2015). However, the raw dataset's distribution is
-noisy (noise also affects the sorting), so using the argument above on the
+significance of signal in a dataset. Namely, when the mean and median of a
+distribution are approximately equal, or the mean's quantile is around 0.5,
+we can argue that there is no significant signal.
+
+To allow for gradients (which are commonly present in ground-based images),
+we can consider the image to be made of a grid of tiles (see
+@ref{Tessellation}@footnote{The options to customize the tessellation are
+discussed in @ref{Processing options}.}). Hence, from the difference of the
+mean and median on each tile, we can estimate the significance of signal in
+it. The median of a distribution is defined to be the value of the
+distribution's middle point after sorting (or 0.5 quantile). Thus, to
+estimate the presence of signal, we'll compare with the quantile of the
+mean with 0.5. If the absolute difference in a tile is larger than the
+value given to the @option{--meanmedqdiff} option, that tile will be
+ignored. You can read this option as ``mean-median-quantile-difference''.
+
+The raw dataset's distribution is noisy, so using the argument above on the
raw input will give a noisy result. To decrease the noise/error in
estimating the mode, we will use convolution (see @ref{Convolution
process}). Convolution decreases the range of the dataset and enhances its
-skewness, See Section 3.1.1 and Figure 4 in Akhlaghi and Ichikawa
-(2015). This enhanced skewness can be interpreted as an increase in the
-Signal to noise ratio of the objects buried in the noise. Therefore, to
-obtain an even better measure of the presence of signal in a mesh, the
-image can be convolved with a given kernel first.
-
-Note that through the difference of the mode and median we have actually
-`detected' data in the distribution. However this ``detection'' was only
-based on the total distribution of the data in each tile (a much lower
-resolution). This is the main limitation of this technique. The best
-approach is thus to do detection over the dataset, mask all the detected
-pixels and use the undetected regions to estimate the sky and its standard
-deviation.
+skewness, See Section 3.1.1 and Figure 4 in
+@url{https://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa (2015)}. This
+enhanced skewness can be interpreted as an increase in the Signal to noise
+ratio of the objects buried in the noise. Therefore, to obtain an even
+better measure of the presence of signal in a tile, the mean and median
+discussed above are measured on the convolved image.
+
+Through the difference of the mean and median we have actually `detected'
+data in the distribution. However this ``detection'' was only based on the
+total distribution of the data in each tile (a much lower resolution). This
+is the main limitation of this technique. The best approach is thus to do
+detection over the dataset, mask all the detected pixels and use the
+undetected regions to estimate the sky and its standard deviation (possibly
+over a tessellation). This is how NoiseChisel works: it uses the argument
+above to find tiles that are used to find its thresholds. Several
+higher-level steps are done on the thresholded pixels to define the
+higher-level detections (see @ref{NoiseChisel}).
@cindex Cosmic rays
-The mean value of the tiles that have an approximately equal mode and
-median will be the Sky value. However there is one final hurdle:
-astronomical datasets are commonly plagued with Cosmic rays. Images of
-Cosmic rays aren't smoothed by the atmosphere or telescope aperture, so
-they have sharp boundaries. Also, since they don't occupy too many pixels,
-they don't affect the mode and median calculation. But their very high
-values can greatly bias the calculation of the mean (recall how the mean
-shifts the fastest in the presence of outliers), see Figure 15 in Akhlaghi
-and Ichikawa (2015) for one example.
+There is one final hurdle: raw astronomical datasets are commonly peppered
+with Cosmic rays. Images of Cosmic rays aren't smoothed by the atmosphere
+or telescope aperture, so they have sharp boundaries. Also, since they
+don't occupy too many pixels, they don't affect the mode and median
+calculation. But their very high values can greatly bias the calculation of
+the mean (recall how the mean shifts the fastest in the presence of
+outliers), for example see Figure 15 in
+@url{https://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa (2015)}.
The effect of outliers like cosmic rays on the mean and standard deviation
can be removed through @mymath{\sigma}-clipping, see @ref{Sigma clipping}
@@ -14678,8 +14942,56 @@ median are approximately equal in a tile (see
@ref{Tessellation}), the
final Sky value and its standard deviation are determined after
@mymath{\sigma}-clipping with the @option{--sigmaclip} option.
+In the end, some of the tiles will pass the mean and median quantile
+difference test. However, prior to interpolating over the failed tiles,
+another point should be considered: large and extended galaxies, or bright
+stars, have wings which sink into the noise very gradually. In some cases,
+the gradient over these wings can be on scales that is larger than the
+tiles. The mean-median distance test will pass on such tiles and will cause
+a strong peak in the interpolated tile grid, see @ref{Detecting large
+extended targets}.
+The tiles that exist over the wings of large galaxies or bright stars are
+outliers in the distribution of tiles that passed the mean-median quantile
+distance test. Therefore, the final step of ``quantifying signal in a
+tile'' is to look at this distribution and remove the
+outliers. @mymath{\sigma}-clipping is a good solution for removing a few
+outliers, but the problem with outliers of this kind is that there may be
+many such tiles (depending on the large/bright stars/galaxies in the
+image). Therefore a novel outlier rejection algorithm will be used.
+@cindex Outliers
+@cindex Identifying outliers
+To identify the first outlier, we'll use the distribution of distances
+between sorted elements. If there are @mymath{N} successful tiles, for
+every tile, the distance between the adjacent @mymath{N/2} previous
+elements is found, giving a distribution containing @mymath{N/2-1}
+points. The @mymath{\sigma}-clipped median and standard deviation of this
+distribution is then found (@mymath{\sigma}-clipping is configured with
+@option{--outliersclip}). Finally, if the distance between the element and
+its previous element is more than @option{--outliersigma} multiples of the
+@mymath{\sigma}-clipped standard deviation added with the
+@mymath{\sigma}-clipped median, that element is considered an outlier and
+all tiles larger than that value are ignored.
+
+Formally, if we assume there are @mymath{N} elements. They are first
+sorted. Searching for the outlier starts on element @mymath{N/2} (integer
+division). Let's take @mymath{v_i} to be the @mymath{i}-th element of the
+sorted input (with no blank values) and @mymath{m} and @mymath{\sigma} as
+the @mymath{\sigma}-clipped median and standard deviation from the
+distances of the previous @mymath{N/2-1} elements (not including
+@mymath{v_i}). If the value given to @option{--outliersigma} is displayed
+with @mymath{s}, the @mymath{i}-th element is considered as an outlier when
+the condition below is true.
+
+@dispmath{{(v_i-v_{i-1})-m\over \sigma}>s}
+
+Since @mymath{i} begins from the median, the outlier has to be larger than
+the median. You can use the check images (for example @option{--checksky}
+in the Statistics program or @option{--checkqthresh},
+@option{--checkdetsky} and @option{--checksky} options in NoiseChisel for
+any of its steps that uses this outlier rejection) to inspect the steps and
+see which tiles have been discarded as outliers prior to interpolation.
@@ -15230,15 +15542,14 @@ tile, see @ref{Quantifying signal in a tile}.
Kernel HDU to help in estimating the significance of signal in a tile, see
@ref{Quantifying signal in a tile}.
-@item --mirrordist=FLT
-Maximum distance (as a multiple of error) to estimate the difference
-between the input and mirror distributions in finding the mode, see
-Appendix C of @url{https://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa
-2015}, also see @ref{Quantifying signal in a tile}.
-
-@item --modmedqdiff=FLT
-The maximum acceptable distance between the mode and median, see
-@ref{Quantifying signal in a tile}.
+@item --meanmedqdiff=FLT
+The maximum acceptable distance between the quantiles of the mean and
+median, see @ref{Quantifying signal in a tile}. The initial Sky and its
+standard deviation estimates are measured on tiles where the quantiles of
+their mean and median are less distant than the value given to this
+option. For example @option{--meanmedqdiff=0.01} means that only tiles
+where the mean's quantile is between 0.49 and 0.51 (recall that the
+median's quantile is 0.5) will be used.
@item --sclipparams=FLT,FLT
The @mymath{\sigma}-clipping parameters, see @ref{Sigma clipping}. This
@@ -15250,6 +15561,23 @@ section). The second value is the exit criteria. If it
is less than 1, then
it is interpreted as tolerance and if it is larger than one it is a
specific number. Hence, in the latter case the value must be an integer.
+@item --outliersclip=FLT,FLT
+Sigma-clipping parameters for the outlier rejection of the Sky value
+(similar to @option{--sclipparams}).
+
+Outlier rejection is useful when the dataset contains a large and diffuse
+(almost flat within each tile) signal. The flatness of the profile will
+cause it to successfully pass the mean-median quantile difference test, so
+we'll need to use the distribution of successful tiles for removing these
+false positive. For more, see the latter half of @ref{Quantifying signal in
+a tile}.
+
+@item --outliersigma=FLT
+Multiple of sigma to define an outlier in the Sky value estimation. If this
+option is given a value of zero, no outlier rejection will take place. For
+more see @option{--outliersclip} and the latter half of @ref{Quantifying
+signal in a tile}.
+
@item --smoothwidth=INT
Width of a flat kernel to convolve the interpolated tile values. Tile
interpolation is done using the median of the @option{--interpnumngb}
@@ -15468,7 +15796,7 @@ Segmentation has been completely moved to a new
program: @ref{Segment}.
@end itemize
@noindent
-Added options:
+Added features/options:
@itemize
@item
@@ -15490,16 +15818,21 @@ calculated on the dataset that is convolved with the
wider kernel, then the
quantiles are estimated on the image convolved with the sharper kernel.
@item
-@option{--noerodequant}: to specify a quantile threshold where erosion
-will not apply. This is useful to detect sharper point-like sources that
-will be missed due to too much erosion. To completely ignore this features
-give this option a value of 1 (only the largest valued pixel in the input
-will not be eroded).
+The quantile difference to identify tiles with no significant signal is
+measured between the @emph{mean} and median. In the published paper, it was
+between the @emph{mode} and median. The quantile of the mean is more
+sensitive to skewness (the presence of signal), so it is preferable to the
+quantile of the mode. For more see @ref{Quantifying signal in a tile}.
@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}.
+Outlier rejection in quantile thresholds: When there are large galaxies or
+bright stars in the image, their gradient may be on a smaller scale than
+the selected tile size. In such cases, those tiles will be identified as
+tiles with no signal and thus preserved. An outlier identification
+algorithm has been added to NoiseChisel and can be configured with the
+following options: @option{--outliersigma} and @option{--outliersclip}. For
+a more complete description, see the latter half of @ref{Quantifying signal
+in a tile}.
@item
@option{--detgrowquant}: is used to grow the final true detections until a
@@ -15860,26 +16193,52 @@ that is discussed in that section.
@subsubsection Detection options
Detection is the process of separating the pixels in the image into two
-groups: 1) Signal and 2) Noise. Through the parameters below, you can
+groups: 1) Signal, and 2) Noise. Through the parameters below, you can
customize the detection process in NoiseChisel. Recall that you can always
-see the full list of Gnuastro's options with the @option{--help} (see
+see the full list of NoiseChisel's options with the @option{--help} (see
@ref{Getting help}), or @option{--printparams} (or @option{-P}) to see
their values (see @ref{Operating mode options}).
@table @option
-@item -r FLT
-@itemx --mirrordist=FLT
-Maximum distance (as a multiple of error) to estimate the difference
-between the input and mirror distributions in finding the mode, see
-Appendix C of @url{https://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa
-2015}, also see @ref{Quantifying signal in a tile}.
-
@item -Q FLT
-@itemx --modmedqdiff=FLT
-The maximum acceptable distance between the mode and median, see
-@ref{Quantifying signal in a tile}. The quantile threshold will be found on
-tiles that satisfy this mode and median difference.
+@itemx --meanmedqdiff=FLT
+The maximum acceptable distance between the quantiles of the mean and
+median in each tile, see @ref{Quantifying signal in a tile}. The quantile
+threshold estimates are measured on tiles where the quantiles of their mean
+and median are less distant than the value given to this option. For
+example @option{--meanmedqdiff=0.01} means that only tiles where the mean's
+quantile is between 0.49 and 0.51 (recall that the median's quantile is
+0.5) will be used.
+
+@item --outliersclip=FLT,FLT
+Sigma-clipping parameters for the outlier rejection of the quantile
+threshold. The format of the given values is similar to
+@option{--sigmaclip} below. In NoiseChisel, outlier rejection is used in
+three stages:
+@itemize
+@item
+Identifying the quantile thresholds (@option{--qthresh},
+@option{--noerodequant}, and @option{detgrowquant}).
+@item
+Identifying the first estimate of the Sky and its standard deviation values
+for pseudo-detections (@option{--dthresh}).
+@item
+Identifying the final estimate of the Sky and its standard deviation.
+@end itemize
+
+Outlier rejection is useful when the dataset contains a large and diffuse
+(almost flat within each tile) signal. The flatness of the profile will
+cause it to successfully pass the mean-median quantile difference test, so
+we'll need to use the distribution of successful tiles for removing these
+false positive. For more, see the latter half of @ref{Quantifying signal in
+a tile}.
+
+@item --outliersigma=FLT
+Multiple of sigma to define an outlier. If this option is given a value of
+zero, no outlier rejection will take place. For more see
+@option{--outliersclip} and the latter half of @ref{Quantifying signal in a
+tile}.
@item -t FLT
@itemx --qthresh=FLT
@@ -15903,44 +16262,6 @@ 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 flatness 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 interpreted 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.
@@ -16023,18 +16344,6 @@ them. Once opening is complete, we have @emph{initial}
detections.
The class of neighbors used for opening, see @option{--erodengb} for more
information on acceptable values.
-@item -s FLT,FLT
-@itemx --sigmaclip=FLT,FLT
-The @mymath{\sigma}-clipping parameters, see @ref{Sigma clipping}. This
-option takes two values which are separated by a comma (@key{,}). Each
-value can either be written as a single number or as a fraction of two
-numbers (for example @code{3,1/10}). The first value to this option is the
-multiple of @mymath{\sigma} that will be clipped (@mymath{\alpha} in that
-section). The second value is the exit criteria. If it is less than 1, then
-it is interpreted as tolerance and if it is larger than one it is assumed
-to be the fixed number of iterations. Hence, in the latter case the value
-must be an integer.
-
@item -B FLT
@itemx --minskyfrac=FLT
Minimum fraction (value between 0 and 1) of Sky (undetected) areas in a
@@ -16063,10 +16372,24 @@ to the final quantile threshold, this behavior can be
disabled with
pixel size as the input, but with the @option{--oneelempertile} option,
only one pixel will be used for each tile (see @ref{Processing options}).
+@item -s FLT,FLT
+@itemx --sigmaclip=FLT,FLT
+The @mymath{\sigma}-clipping parameters for measuing the initial and final
+Sky values from the undetected pixels, see @ref{Sigma clipping}.
+
+This option takes two values which are separated by a comma (@key{,}). Each
+value can either be written as a single number or as a fraction of two
+numbers (for example @code{3,1/10}). The first value to this option is the
+multiple of @mymath{\sigma} that will be clipped (@mymath{\alpha} in that
+section). The second value is the exit criteria. If it is less than 1, then
+it is interpreted as tolerance and if it is larger than one it is assumed
+to be the fixed number of iterations. Hence, in the latter case the value
+must be an integer.
+
@item -R FLT
@itemx --dthresh=FLT
-The detection threshold: a multiple of the initial sky standard deviation
-added with the initial sky approximation (which you can inspect with
+The detection threshold: a multiple of the initial Sky standard deviation
+added with the initial Sky approximation (which you can inspect with
@option{--checkdetsky}). This flux threshold is applied to the initially
undetected regions on the unconvolved image. The background pixels that are
completely engulfed in a 4-connected foreground region are converted to
@@ -22695,8 +23018,7 @@ function uses C's @code{sizeof} operator to measure the
size of each type.
@end deftypefun
@deftypefun {char *} gal_type_name (uint8_t @code{type}, int @code{long_name})
-Return a string that contains the name of @code{type}. This can be used in
-messages to the users when your function/program accepts many types. It can
+Return a string literal that contains the name of @code{type}. It can
return both short and long formats of the type names (for example
@code{f32} and @code{float32}). If @code{long_name} is non-zero, the long
format will be returned, otherwise the short name will be returned. The
@@ -23026,6 +23348,14 @@ that corresponds to its type. If @code{input} is a
tile over a larger
dataset, only the region that the tile covers will be set to blank.
@end deftypefun
+@deftypefun {char *} gal_blank_as_string (uint8_t @code{type}, int
@code{width})
+Write the blank value for the given data type @code{type} into a string and
+return it. The space for the string is dynamically allocated so it must be
+freed after you are done with it. If @code{width!=0}, then the final string
+will be padded with white space characters to have the requested width if
+it is smaller.
+@end deftypefun
+
@deftypefun int gal_blank_is (void @code{*pointer}, uint8_t @code{type})
Return 1 if the contents of @code{pointer} (assuming a type of @code{type})
is blank. Otherwise, return 0. Note that this function only works on one
@@ -23069,6 +23399,12 @@ Create a dataset of the the same size as the input,
but with an
those that aren't.
@end deftypefun
+@deftypefun void gal_blank_flag_apply (gal_data_t @code{*input}, gal_data_t
@code{*flag})
+Set all non-zero and non-blank elements of @code{flag} to blank in
+@code{input}. @code{flag} has to have an unsigned 8-bit type and be the
+same size as @code{input}.
+@end deftypefun
+
@deftypefun void gal_blank_remove (gal_data_t @code{*input})
Remove blank elements from a dataset, convert it to a 1D dataset, adjust
@@ -23085,11 +23421,10 @@ appropriate measure. This check is highly recommended
because it will avoid
strange bugs in later steps.
@end deftypefun
-@deftypefun {char *} gal_blank_as_string (uint8_t @code{type}, int
@code{width})
-Write the blank value for the given data type @code{type} into a string and
-return it. The space for the string is dynamically allocated so it must be
-freed after you are done with it.
-@end deftypefun
+
+
+
+
@node Library data container, Dimensions, Library blank values, Gnuastro
library
@subsection Data container (@file{data.h})
@@ -27540,6 +27875,55 @@ above.
@end deftypefun
+@deftypefun {gal_data_t *} gal_statistics_outlier_positive (gal_data_t
@code{*input}, float @code{sigma}, float @code{sigclip_multip}, float
@code{sigclip_param}, int @code{inplace}, int @code{quiet})
+Find the first positive outlier in the @code{input} distribution. The
+returned dataset contains a single element: the first positive outlier. It
+is one of the dataset's elements, in the same type as the input. If the
+process fails for any reason (for example no outlier was found), a
+@code{NULL} pointer will be returned.
+
+All (possibly existing) blank elements are first removed from the input
+dataset, then it is sorted. A sliding window equal to half the width of the
+dataset elements is parsed over the dataset. Starting from the middle of
+the dataset (median) in the direction of increasing values. This window is
+used as a reference. The first element where the distance to the previous
+(sorted) element is @code{sigma} units away from the distribution of
+distances in its window is considered an outlier and returned by this
+function.
+
+Formally, if we assume there are @mymath{N} non-blank elements. They are
+first sorted. Searching for the outlier starts on element @mymath{N/2}
+(integer division). Let's take @mymath{v_i} to be the @mymath{i}-th element
+of the sorted input (with no blank values) and @mymath{m} and
+@mymath{\sigma} as the @mymath{\sigma}-clipped median and standard
+deviation from the distances of the previous @mymath{N/2} elements (not
+including @mymath{v_i}). If the value given to @code{sigma} is displayed
+with @mymath{s}, the @mymath{i}-th element is considered as an outlier when
+the condition below is true.
+
+@dispmath{{(v_i-v_{i-1})-m\over \sigma}>s}
+
+The @code{sigclip_multip} and @code{sigclip_param} arguments specify the
+properties of the @mymath{\sigma}-clipping (see @ref{Sigma clipping} for
+more). You see that by this definition, the outlier cannot be any of the
+lower half elements. The advantage of this algorithm compared to
+@mymath{\sigma}-clippign is that it only looks backwards (in the sorted
+array) and parses it in one direction.
+
+If @code{inplace!=0}, the removing of blank elements and sorting will be
+done within the input dataset's allocated space. Otherwise, this function
+will internally allocate (and later free) the necessary space to keep the
+intermediate space that this process requires.
+
+If @code{quiet!=0}, this function will report a warning when the moving
+window size is not much smaller than the number of non-blank elements. In
+such cases, the outlier estimation may be wrong: the outlier point may lie
+before the first checked element.
+
+@end deftypefun
+
+
+
@node Binary datasets, Labeled datasets, Statistical operations, Gnuastro
library
@subsection Binary datasets (@file{binary.h})
diff --git a/lib/Makefile.am b/lib/Makefile.am
index a6ee75a..8a36723 100644
--- a/lib/Makefile.am
+++ b/lib/Makefile.am
@@ -60,8 +60,8 @@ libgnuastro_la_SOURCES = arithmetic.c arithmetic-and.c
arithmetic-bitand.c \
checkset.c convolve.c cosmology.c data.c eps.c fits.c git.c \
interpolate.c jpeg.c label.c list.c match.c options.c pdf.c \
permutation.c pointer.c polygon.c qsort.c dimension.c statistics.c \
- table.c tableintern.c threads.c tiff.c tile.c timing.c txt.c type.c \
- wcs.c
+ table.c tableintern.c threads.c tiff.c tile.c tile-internal.c timing.c \
+ txt.c type.c wcs.c
@@ -107,7 +107,7 @@ EXTRA_DIST = gnuastro.pc.in $(headersdir)/README
$(internaldir)/README \
$(internaldir)/checkset.h $(internaldir)/commonopts.h \
$(internaldir)/config.h.in $(internaldir)/fixedstringmacros.h \
$(internaldir)/options.h $(internaldir)/tableintern.h \
- $(internaldir)/timing.h
+ $(internaldir)/tile-internal.h $(internaldir)/timing.h
diff --git a/lib/blank.c b/lib/blank.c
index b84a5ce..b3c2bb7 100644
--- a/lib/blank.c
+++ b/lib/blank.c
@@ -113,6 +113,184 @@ gal_blank_initialize(gal_data_t *input)
+/* Print the blank value as a string. For the integer types, we'll use the
+ PRIxNN keywords of `inttypes.h' (which is imported into Gnuastro from
+ Gnulib, so we don't necessarily rely on the host system having it). */
+char *
+gal_blank_as_string(uint8_t type, int width)
+{
+ char *blank=NULL, *fmt;
+
+ /* Print the given value. */
+ switch(type)
+ {
+ case GAL_TYPE_BIT:
+ error(EXIT_FAILURE, 0, "%s: bit types are not implemented yet",
+ __func__);
+ break;
+
+ case GAL_TYPE_STRING:
+ if(width)
+ {
+ if( asprintf(&blank, "%*s", width, GAL_BLANK_STRING)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ else
+ {
+ if( asprintf(&blank, "%s", GAL_BLANK_STRING)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ break;
+
+ case GAL_TYPE_UINT8:
+ fmt = width ? "%*"PRIu8 : "%"PRIu8;
+ if(width)
+ {
+ if( asprintf(&blank, fmt, width, (uint8_t)GAL_BLANK_UINT8)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ else
+ {
+ if( asprintf(&blank, fmt, (uint8_t)GAL_BLANK_UINT8)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ break;
+
+ case GAL_TYPE_INT8:
+ fmt = width ? "%*"PRId8 : "%"PRId8;
+ if(width)
+ {
+ if( asprintf(&blank, fmt, width, (int8_t)GAL_BLANK_INT8)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ else
+ {
+ if( asprintf(&blank, fmt, (int8_t)GAL_BLANK_INT8)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ break;
+
+ case GAL_TYPE_UINT16:
+ fmt = width ? "%*"PRIu16 : "%"PRIu16;
+ if(width)
+ {
+ if( asprintf(&blank, fmt, width, (uint16_t)GAL_BLANK_UINT16)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ else
+ {
+ if( asprintf(&blank, fmt, (uint16_t)GAL_BLANK_UINT16)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ break;
+
+ case GAL_TYPE_INT16:
+ fmt = width ? "%*"PRId16 : "%"PRId16;
+ if(width)
+ {
+ if( asprintf(&blank, fmt, width, (int16_t)GAL_BLANK_INT16)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ else
+ {
+ if( asprintf(&blank, fmt, (int16_t)GAL_BLANK_INT16)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ break;
+
+ case GAL_TYPE_UINT32:
+ fmt = width ? "%*"PRIu32 : "%"PRIu32;
+ if(width)
+ {
+ if( asprintf(&blank, fmt, width, (uint32_t)GAL_BLANK_UINT32)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ else
+ {
+ if( asprintf(&blank, fmt, (uint32_t)GAL_BLANK_UINT32)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ break;
+
+ case GAL_TYPE_INT32:
+ fmt = width ? "%*"PRId32 : "%"PRId32;
+ if(width)
+ {
+ if( asprintf(&blank, fmt, width, (int32_t)GAL_BLANK_INT32)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ else
+ {
+ if( asprintf(&blank, fmt, (int32_t)GAL_BLANK_INT32)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ break;
+
+ case GAL_TYPE_UINT64:
+ fmt = width ? "%*"PRIu64 : "%"PRIu64;
+ if(width)
+ {
+ if( asprintf(&blank, fmt, width, (uint64_t)GAL_BLANK_UINT64)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ else
+ {
+ if( asprintf(&blank, fmt, (uint64_t)GAL_BLANK_UINT64)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ break;
+
+ case GAL_TYPE_INT64:
+ fmt = width ? "%*"PRId64 : "%"PRId64;
+ if(width)
+ {
+ if( asprintf(&blank, fmt, width, (int64_t)GAL_BLANK_INT64)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ else
+ {
+ if( asprintf(&blank, fmt, (int64_t)GAL_BLANK_INT64)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ break;
+
+ case GAL_TYPE_FLOAT32:
+ if(width)
+ {
+ if( asprintf(&blank, "%*f", width, (float)GAL_BLANK_FLOAT32)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ else
+ {
+ if( asprintf(&blank, "%f", (float)GAL_BLANK_FLOAT32)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ break;
+
+ case GAL_TYPE_FLOAT64:
+ if(width)
+ {
+ if( asprintf(&blank, "%*f", width, (double)GAL_BLANK_FLOAT64)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ else
+ {
+ if( asprintf(&blank, "%f", (double)GAL_BLANK_FLOAT64)<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ }
+ break;
+
+ default:
+ error(EXIT_FAILURE, 0, "%s: type code %d not recognized", __func__,
+ type);
+ }
+ return blank;
+}
+
+
+
+
+
/* Return 1 if the contents of the pointer (with the given type) is
blank. */
@@ -318,7 +496,7 @@ gal_blank_flag(gal_data_t *input)
/* Go over the pixels and set the output values. */
switch(input->type)
{
- /* Numeric types */
+ /* Numeric types */
case GAL_TYPE_UINT8: FLAG_BLANK( uint8_t ); break;
case GAL_TYPE_INT8: FLAG_BLANK( int8_t ); break;
case GAL_TYPE_UINT16: FLAG_BLANK( uint16_t ); break;
@@ -330,19 +508,19 @@ gal_blank_flag(gal_data_t *input)
case GAL_TYPE_FLOAT32: FLAG_BLANK( float ); break;
case GAL_TYPE_FLOAT64: FLAG_BLANK( double ); break;
- /* String. */
+ /* String. */
case GAL_TYPE_STRING:
do *o++ = !strcmp(*str,GAL_BLANK_STRING); while(++str<strf);
break;
- /* Currently unsupported types. */
+ /* Currently unsupported types. */
case GAL_TYPE_BIT:
case GAL_TYPE_COMPLEX32:
case GAL_TYPE_COMPLEX64:
error(EXIT_FAILURE, 0, "%s: %s type not yet supported",
__func__, gal_type_name(input->type, 1));
- /* Bad input. */
+ /* Bad input. */
default:
error(EXIT_FAILURE, 0, "%s: type value (%d) not recognized",
__func__, input->type);
@@ -361,6 +539,78 @@ gal_blank_flag(gal_data_t *input)
+/* Write a blank value in the input anywhere that the flag dataset is not
+ zero or not blank. */
+#define BLANK_FLAG_APPLY(IT) { \
+ IT *i=input->array, b; \
+ gal_blank_write(&b, input->type); \
+ do {if(*f && *f!=GAL_BLANK_UINT8) *i=b; ++i;} while(++f<ff); \
+ }
+void
+gal_blank_flag_apply(gal_data_t *input, gal_data_t *flag)
+{
+ char **str=input->array;
+ uint8_t *f=flag->array, *ff=f+flag->size;
+
+ /* Sanity check. */
+ if(flag->type!=GAL_TYPE_UINT8)
+ error(EXIT_FAILURE, 0, "%s: the `flag' argument has a `%s' type, it "
+ "must have an unsigned 8-bit type", __func__,
+ gal_type_name(flag->type, 1));
+ if(gal_dimension_is_different(input, flag))
+ error(EXIT_FAILURE, 0, "%s: the `flag' argument doesn't have the same "
+ "size as the `input' argument", __func__);
+
+ /* Write the blank values. */
+ switch(input->type)
+ {
+ /* Basic numeric types. */
+ case GAL_TYPE_UINT8: BLANK_FLAG_APPLY( uint8_t ); break;
+ case GAL_TYPE_INT8: BLANK_FLAG_APPLY( int8_t ); break;
+ case GAL_TYPE_UINT16: BLANK_FLAG_APPLY( uint16_t ); break;
+ case GAL_TYPE_INT16: BLANK_FLAG_APPLY( int16_t ); break;
+ case GAL_TYPE_UINT32: BLANK_FLAG_APPLY( uint32_t ); break;
+ case GAL_TYPE_INT32: BLANK_FLAG_APPLY( int32_t ); break;
+ case GAL_TYPE_UINT64: BLANK_FLAG_APPLY( uint64_t ); break;
+ case GAL_TYPE_INT64: BLANK_FLAG_APPLY( int64_t ); break;
+ case GAL_TYPE_FLOAT32: BLANK_FLAG_APPLY( float ); break;
+ case GAL_TYPE_FLOAT64: BLANK_FLAG_APPLY( double ); break;
+
+ /* Strings. */
+ case GAL_TYPE_STRING:
+ do
+ {
+ if(*f && *f!=GAL_BLANK_UINT8)
+ {
+ free(*str);
+ *str=gal_blank_as_string(GAL_TYPE_STRING, 0);
+ }
+ ++str;
+ }
+ while(++f<ff);
+ break;
+
+ /* Currently unsupported types. */
+ case GAL_TYPE_BIT:
+ case GAL_TYPE_COMPLEX32:
+ case GAL_TYPE_COMPLEX64:
+ error(EXIT_FAILURE, 0, "%s: %s type not yet supported",
+ __func__, gal_type_name(input->type, 1));
+
+ /* Bad input. */
+ default:
+ error(EXIT_FAILURE, 0, "%s: type value (%d) not recognized",
+ __func__, input->type);
+ }
+
+ /* Update the blank flags. */
+ gal_blank_present(input, 1);
+}
+
+
+
+
+
/* Remove blank elements from a dataset, convert it to a 1D dataset and
adjust the size properly. In practice this function doesn't `realloc'
the input array, all it does is to shift the blank eleemnts to the end
@@ -415,181 +665,3 @@ gal_blank_remove(gal_data_t *input)
input->flag |= GAL_DATA_FLAG_BLANK_CH;
input->flag &= ~GAL_DATA_FLAG_HASBLANK;
}
-
-
-
-
-
-/* Print the blank value as a string. For the integer types, we'll use the
- PRIxNN keywords of `inttypes.h' (which is imported into Gnuastro from
- Gnulib, so we don't necessarily rely on the host system having it). */
-char *
-gal_blank_as_string(uint8_t type, int width)
-{
- char *blank=NULL, *fmt;
-
- /* Print the given value. */
- switch(type)
- {
- case GAL_TYPE_BIT:
- error(EXIT_FAILURE, 0, "%s: bit types are not implemented yet",
- __func__);
- break;
-
- case GAL_TYPE_STRING:
- if(width)
- {
- if( asprintf(&blank, "%*s", width, GAL_BLANK_STRING)<0 )
- error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
- }
- else
- {
- if( asprintf(&blank, "%s", GAL_BLANK_STRING)<0 )
- error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
- }
- break;
-
- case GAL_TYPE_UINT8:
- fmt = width ? "%*"PRIu8 : "%"PRIu8;
- if(width)
- {
- if( asprintf(&blank, fmt, width, (uint8_t)GAL_BLANK_UINT8)<0 )
- error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
- }
- else
- {
- if( asprintf(&blank, fmt, (uint8_t)GAL_BLANK_UINT8)<0 )
- error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
- }
- break;
-
- case GAL_TYPE_INT8:
- fmt = width ? "%*"PRId8 : "%"PRId8;
- if(width)
- {
- if( asprintf(&blank, fmt, width, (int8_t)GAL_BLANK_INT8)<0 )
- error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
- }
- else
- {
- if( asprintf(&blank, fmt, (int8_t)GAL_BLANK_INT8)<0 )
- error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
- }
- break;
-
- case GAL_TYPE_UINT16:
- fmt = width ? "%*"PRIu16 : "%"PRIu16;
- if(width)
- {
- if( asprintf(&blank, fmt, width, (uint16_t)GAL_BLANK_UINT16)<0 )
- error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
- }
- else
- {
- if( asprintf(&blank, fmt, (uint16_t)GAL_BLANK_UINT16)<0 )
- error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
- }
- break;
-
- case GAL_TYPE_INT16:
- fmt = width ? "%*"PRId16 : "%"PRId16;
- if(width)
- {
- if( asprintf(&blank, fmt, width, (int16_t)GAL_BLANK_INT16)<0 )
- error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
- }
- else
- {
- if( asprintf(&blank, fmt, (int16_t)GAL_BLANK_INT16)<0 )
- error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
- }
- break;
-
- case GAL_TYPE_UINT32:
- fmt = width ? "%*"PRIu32 : "%"PRIu32;
- if(width)
- {
- if( asprintf(&blank, fmt, width, (uint32_t)GAL_BLANK_UINT32)<0 )
- error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
- }
- else
- {
- if( asprintf(&blank, fmt, (uint32_t)GAL_BLANK_UINT32)<0 )
- error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
- }
- break;
-
- case GAL_TYPE_INT32:
- fmt = width ? "%*"PRId32 : "%"PRId32;
- if(width)
- {
- if( asprintf(&blank, fmt, width, (int32_t)GAL_BLANK_INT32)<0 )
- error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
- }
- else
- {
- if( asprintf(&blank, fmt, (int32_t)GAL_BLANK_INT32)<0 )
- error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
- }
- break;
-
- case GAL_TYPE_UINT64:
- fmt = width ? "%*"PRIu64 : "%"PRIu64;
- if(width)
- {
- if( asprintf(&blank, fmt, width, (uint64_t)GAL_BLANK_UINT64)<0 )
- error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
- }
- else
- {
- if( asprintf(&blank, fmt, (uint64_t)GAL_BLANK_UINT64)<0 )
- error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
- }
- break;
-
- case GAL_TYPE_INT64:
- fmt = width ? "%*"PRId64 : "%"PRId64;
- if(width)
- {
- if( asprintf(&blank, fmt, width, (int64_t)GAL_BLANK_INT64)<0 )
- error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
- }
- else
- {
- if( asprintf(&blank, fmt, (int64_t)GAL_BLANK_INT64)<0 )
- error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
- }
- break;
-
- case GAL_TYPE_FLOAT32:
- if(width)
- {
- if( asprintf(&blank, "%*f", width, (float)GAL_BLANK_FLOAT32)<0 )
- error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
- }
- else
- {
- if( asprintf(&blank, "%f", (float)GAL_BLANK_FLOAT32)<0 )
- error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
- }
- break;
-
- case GAL_TYPE_FLOAT64:
- if(width)
- {
- if( asprintf(&blank, "%*f", width, (double)GAL_BLANK_FLOAT64)<0 )
- error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
- }
- else
- {
- if( asprintf(&blank, "%f", (double)GAL_BLANK_FLOAT64)<0 )
- error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
- }
- break;
-
- default:
- error(EXIT_FAILURE, 0, "%s: type code %d not recognized", __func__,
- type);
- }
- return blank;
-}
diff --git a/lib/gnuastro-internal/arithmetic-internal.h
b/lib/gnuastro-internal/arithmetic-internal.h
index 59a5528..0736465 100644
--- a/lib/gnuastro-internal/arithmetic-internal.h
+++ b/lib/gnuastro-internal/arithmetic-internal.h
@@ -71,4 +71,4 @@ gal_arithmetic_convert_to_compiled_type(gal_data_t *in, int
flags);
__END_C_DECLS /* From C++ preparations */
-#endif /* __GAL_ARITHMETIC_H__ */
+#endif /* __GAL_ARITHMETIC_INTERNAL_H__ */
diff --git a/lib/gnuastro-internal/arithmetic-internal.h
b/lib/gnuastro-internal/tile-internal.h
similarity index 78%
copy from lib/gnuastro-internal/arithmetic-internal.h
copy to lib/gnuastro-internal/tile-internal.h
index 59a5528..83c55ec 100644
--- a/lib/gnuastro-internal/arithmetic-internal.h
+++ b/lib/gnuastro-internal/tile-internal.h
@@ -1,11 +1,12 @@
/*********************************************************************
-Arithmetic -- Preform arithmetic operations on datasets.
+Common tile operations used by some Gnuastro programs, but too specific
+to be in the general library.
This is part of GNU Astronomy Utilities (Gnuastro) package.
Original author:
Mohammad Akhlaghi <mohammad@akhlaghi.org>
Contributing author(s):
-Copyright (C) 2017-2018, Free Software Foundation, Inc.
+Copyright (C) 2018, Free Software Foundation, Inc.
Gnuastro is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
@@ -20,12 +21,11 @@ General Public License for more details.
You should have received a copy of the GNU General Public License
along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
**********************************************************************/
-#ifndef __GAL_ARITHMETIC_INTERNAL_H__
-#define __GAL_ARITHMETIC_INTERNAL_H__
+#ifndef __GAL_TILE_INTERNAL_H__
+#define __GAL_TILE_INTERNAL_H__
/* Include other headers if necessary here. Note that other header files
must be included before the C++ preparations below */
-#include <gnuastro/data.h>
/* When we are within Gnuastro's building process, `IN_GNUASTRO_BUILD' is
@@ -58,17 +58,14 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
/* Actual header contants (the above were for the Pre-processor). */
__BEGIN_C_DECLS /* From C++ preparations */
-int
-gal_arithmetic_binary_checkblank(gal_data_t *l, gal_data_t *r);
-
-char *
-gal_arithmetic_operator_string(int operator);
-
-gal_data_t *
-gal_arithmetic_convert_to_compiled_type(gal_data_t *in, int flags);
-
+void
+gal_tileinternal_no_outlier(gal_data_t *first, gal_data_t *second,
+ gal_data_t *third,
+ struct gal_tile_two_layer_params *tl,
+ double *outliersclip, float outliersigma,
+ char *filename);
__END_C_DECLS /* From C++ preparations */
-#endif /* __GAL_ARITHMETIC_H__ */
+#endif /* __GAL_TILE_INTERNAL_H__ */
diff --git a/lib/gnuastro/blank.h b/lib/gnuastro/blank.h
index eab1fbf..7e36214 100644
--- a/lib/gnuastro/blank.h
+++ b/lib/gnuastro/blank.h
@@ -105,6 +105,9 @@ gal_data_t *
gal_blank_flag(gal_data_t *data);
void
+gal_blank_flag_apply(gal_data_t *input, gal_data_t *flag);
+
+void
gal_blank_remove(gal_data_t *data);
diff --git a/lib/gnuastro/statistics.h b/lib/gnuastro/statistics.h
index ad96d01..2b84756 100644
--- a/lib/gnuastro/statistics.h
+++ b/lib/gnuastro/statistics.h
@@ -173,7 +173,10 @@ gal_data_t *
gal_statistics_sigma_clip(gal_data_t *input, float multip, float param,
int inplace, int quiet);
-
+gal_data_t *
+gal_statistics_outlier_positive(gal_data_t *input, float sigma,
+ float sigclip_multip, float sigclip_param,
+ int inplace, int quiet);
diff --git a/lib/statistics.c b/lib/statistics.c
index ce77950..ea3dee5 100644
--- a/lib/statistics.c
+++ b/lib/statistics.c
@@ -2070,61 +2070,101 @@ gal_statistics_sigma_clip(gal_data_t *input, float
multip, float param,
-#if 0
-/* Using the cumulative distribution function this function will
- remove outliers from a dataset. */
-void
-gal_statistics_remove_outliers_flat_cdf(float *sorted, size_t *outsize)
+/* Find the first outlier in a distribution. */
+#define OUTLIER_BYTYPE(IT) { \
+ IT *arr=nbs->array; \
+ for(i=nbs->size/2;i<nbs->size;++i) \
+ { \
+ /* Fill in the distance array. */ \
+ for(j=0; j<wtakeone; ++j) \
+ darr[j] = arr[i-window_size+j+1] - arr[i-window_size+j]; \
+ \
+ /* Get the sigma-clipped information. */ \
+ sclip=gal_statistics_sigma_clip(dist, sigclip_multip, \
+ sigclip_param, 0, 1); \
+ sarr=sclip->array; \
+ \
+ /* For a check */ \
+ /*printf("%f: %f (%f, %f) %f\n", (float)(arr[i]), */ \
+ /* (float)(arr[i]-arr[i-1]), sarr[1], sarr[3], */ \
+ /* (((double)(arr[i]-arr[i-1])) - sarr[1])/sarr[3]); */ \
+ \
+ /* Terminate the loop if the dist. is larger than requested. */ \
+ /* This shows we have reached the first outlier's position. */ \
+ if( (((double)(arr[i]-arr[i-1])) - sarr[1]) > sigma*sarr[3] ) \
+ { \
+ /* Allocate the output dataset. */ \
+ out=gal_data_alloc(NULL, input->type, 1, &one, NULL, 0, -1, \
+ NULL, NULL, NULL); \
+ \
+ /* Write the outlier, clean up and break. */ \
+ *(IT *)(out->array)=arr[i-1]; \
+ gal_data_free(sclip); \
+ break; \
+ } \
+ \
+ /* Clean up (if we get here). */ \
+ gal_data_free(sclip); \
+ } \
+ }
+gal_data_t *
+gal_statistics_outlier_positive(gal_data_t *input, float sigma,
+ float sigclip_multip, float sigclip_param,
+ int inplace, int quiet)
{
- printf("\n ... in gal_statistics_remove_outliers_flat_cdf ... \n");
- exit(1);
- int firstfound=0;
- size_t size=*outsize, i, maxind;
- float *slopes, minslope, maxslope;
-
- /* Find a slopes array, think of the cumulative frequency plot when
- you want to think about slopes. */
- errno=0; slopes=malloc(size*sizeof *slopes);
- if(slopes==NULL)
- error(EXIT_FAILURE, errno, "%s: %zu bytes for slopes",
- __func__, size*sizeof *slopes);
-
- /* Calcuate the slope of the CDF and put it in the slopes array. */
- for(i=1;i<size-1;++i)
- slopes[i]=2/(sorted[i+1]-sorted[i-1]);
-
- /* Find the position of the maximum slope, note that around the
- distribution mode, the difference between the values varies less,
- so two neighbouring elements have the closest values, hence the
- largest slope (when their difference is in the denominator). */
- gal_statistics_f_max_with_index(slopes+1, size-2, &maxslope, &maxind);
-
- /* Find the minimum slope from the second element (for the first the
- slope is not defined. NOTE; maxind is one smaller than it should
- be because the input array to find it began from the second
- element. */
- gal_statistics_float_second_min(slopes+1, maxind+1, &minslope);
-
- /* Find the second place where the slope falls below `minslope`
- after the maximum position. When found, add it with one to
- account for error. Note that incase there are no outliers by this
- definition, then the final size will be equal to size! */
- for(i=maxind+1;i<size-1;++i)
- if(slopes[i]<minslope)
- {
- if(firstfound)
- break;
- else
- firstfound=1;
- }
- *outsize=i+1;
+ float *sarr;
+ double *darr;
+ size_t i, j, one=1, wtakeone, window_size;
+ gal_data_t *dist, *sclip, *nbs, *out=NULL;
+
+ /* Remove all blanks and sort the dataset. */
+ nbs=gal_statistics_no_blank_sorted(input, inplace);
+ window_size=nbs->size/2;
+
+ /* Only continue if the window size is more than 2 elements (out
+ "outlier" is hard to define on smaller datasets). */
+ if(window_size>2)
+ {
+ /* For a check.
+ if(nbs->type==GAL_TYPE_FLOAT32)
+ {
+ float *n=nbs->array;
+ for(i=0;i<nbs->size;++i)
+ printf("%f\n", n[i]);
+ exit(0);
+ }
+ */
- /*
- for(i=0;i<size;++i)
- printf("%zu\t%.3f\t%.3f\n", i, arr[i], slopes[i]);
- printf("\n\nPlace to cut off for outliers is: %zu\n\n", *outsize);
- */
- free(slopes);
+ /* Allocate space to keep the distances. */
+ wtakeone=window_size-1;
+ dist=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &wtakeone, NULL,
+ 0, -1, NULL, NULL, NULL);
+ darr=dist->array;
+
+ /* Find the outlier based on the type of the input dataset. */
+ switch(input->type)
+ {
+ case GAL_TYPE_UINT8: OUTLIER_BYTYPE( uint8_t ); break;
+ case GAL_TYPE_INT8: OUTLIER_BYTYPE( int8_t ); break;
+ case GAL_TYPE_UINT16: OUTLIER_BYTYPE( uint16_t ); break;
+ case GAL_TYPE_INT16: OUTLIER_BYTYPE( int16_t ); break;
+ case GAL_TYPE_UINT32: OUTLIER_BYTYPE( uint32_t ); break;
+ case GAL_TYPE_INT32: OUTLIER_BYTYPE( int32_t ); break;
+ case GAL_TYPE_UINT64: OUTLIER_BYTYPE( uint64_t ); break;
+ case GAL_TYPE_INT64: OUTLIER_BYTYPE( int64_t ); break;
+ case GAL_TYPE_FLOAT32: OUTLIER_BYTYPE( float ); break;
+ case GAL_TYPE_FLOAT64: OUTLIER_BYTYPE( double ); break;
+ default:
+ error(EXIT_FAILURE, 0, "%s: type code %d not recognized",
+ __func__, input->type);
+ }
+
+ /* Clean up. */
+ gal_data_free(dist);
+ }
+
+ /* Clean up and return. */
+ if(nbs!=input) gal_data_free(nbs);
+ return out;
}
-#endif
diff --git a/lib/tile-internal.c b/lib/tile-internal.c
new file mode 100644
index 0000000..114b02d
--- /dev/null
+++ b/lib/tile-internal.c
@@ -0,0 +1,191 @@
+/*********************************************************************
+Common tile operations used by some Gnuastro programs, but too specific
+to be in the general library.
+This is part of GNU Astronomy Utilities (Gnuastro) package.
+
+Original author:
+ Mohammad Akhlaghi <mohammad@akhlaghi.org>
+Contributing author(s):
+Copyright (C) 2018, Free Software Foundation, Inc.
+
+Gnuastro is free software: you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation, either version 3 of the License, or (at your
+option) any later version.
+
+Gnuastro is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
+**********************************************************************/
+#include <config.h>
+
+#include <stdio.h>
+#include <errno.h>
+#include <error.h>
+#include <stdlib.h>
+
+#include <gnuastro/tile.h>
+#include <gnuastro/pointer.h>
+#include <gnuastro/statistics.h>
+
+#include <gnuastro-internal/tile-internal.h>
+
+
+/* The main working function for `threshold_no_outlier'. The main
+ purpose/problem is this: when we have channels, the qthresh values for
+ each channel should be treated independently. */
+static void
+tileinternal_no_outlier_work(gal_data_t *first, gal_data_t *second,
+ gal_data_t *third, size_t channelid,
+ size_t tottilesinch, double *outliersclip,
+ float outliersigma)
+{
+ gal_data_t *outlier;
+ size_t i, osize=first->size;
+ size_t start=tottilesinch*channelid;
+ float *oa1=NULL, *oa2=NULL, *oa3=NULL;
+ float o, *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 || tottilesinch!=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_pointer_increment(first->array, start, first->type);
+ second->array=gal_pointer_increment(second->array, start,
+ second->type);
+ if(third)
+ third->array=gal_pointer_increment(third->array, start, third->type);
+
+ /* Correct their sizes. */
+ first->size=tottilesinch;
+ second->size=tottilesinch;
+ if(third) third->size=tottilesinch;
+ }
+
+ /* Find the quantile and remove all tiles that are more than it in the
+ first array. */
+ arr1=first->array;
+ outlier=gal_statistics_outlier_positive(first, outliersigma,
+ outliersclip[0], outliersclip[1],
+ 0, 1);
+ if(outlier)
+ {
+ o = *((float *)(outlier->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]<=o ? arr1[i] : NAN;
+ gal_data_free(outlier);
+ }
+
+ /* Second quantile threshold. We are finding the outliers independently
+ on each dataset to later remove any tile that is blank in atleast one
+ of them. */
+ arr2=second->array;
+ outlier=gal_statistics_outlier_positive(second, outliersigma,
+ outliersclip[0], outliersclip[1],
+ 0, 0);
+ if(outlier)
+ {
+ o = *((float *)(outlier->array));
+ for(i=0;i<second->size;++i)
+ arr2[i] = arr2[i]<=o ? arr2[i] : NAN;
+ gal_data_free(outlier);
+ }
+
+ /* The third (if it exists). */
+ if(third)
+ {
+ arr3=third->array;
+ outlier=gal_statistics_outlier_positive(third, outliersigma,
+ outliersclip[0],
+ outliersclip[1], 0, 0);
+ if(outlier)
+ {
+ o = *((float *)(outlier->array));
+ for(i=0;i<third->size;++i)
+ arr3[i] = arr3[i]<=o ? arr3[i] : NAN;
+ gal_data_free(outlier);
+ }
+ }
+
+ /* Make sure all three have the same NaN pixels. */
+ for(i=0;i<first->size;++i)
+ if( isnan(arr1[i]) || isnan(arr2[i]) || (third && isnan(arr3[i])) )
+ {
+ arr1[i] = arr2[i] = NAN;
+ if(third) arr3[i] = NAN;
+ }
+
+
+ /* Correct the values, if they were changed. */
+ if(start || tottilesinch!=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. */
+void
+gal_tileinternal_no_outlier(gal_data_t *first, gal_data_t *second,
+ gal_data_t *third,
+ struct gal_tile_two_layer_params *tl,
+ double *outliersclip, float outliersigma,
+ char *filename)
+{
+ size_t i;
+
+ /* 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);
+
+ /* Do the work. */
+ for(i=0;i<tl->totchannels;++i)
+ tileinternal_no_outlier_work(first, second, third, i, tl->tottilesinch,
+ outliersclip, outliersigma);
+
+ /* If the user wants to see the steps. */
+ if(filename)
+ {
+ first->name="VALUE1_NO_OUTLIER";
+ second->name="VALUE2_NO_OUTLIER";
+ gal_tile_full_values_write(first, tl, 1, filename, NULL, NULL);
+ gal_tile_full_values_write(second, tl, 1, filename, NULL, NULL);
+ first->name=second->name=NULL;
+ if(third)
+ {
+ third->name="VALUE3_NO_OUTLIER";
+ gal_tile_full_values_write(third, tl, 1, filename, NULL, NULL);
+ third->name=NULL;
+ }
+ }
+}
diff --git a/tests/during-dev.sh b/tests/during-dev.sh
index 7a52e82..c3ffb8c 100755
--- a/tests/during-dev.sh
+++ b/tests/during-dev.sh
@@ -1,4 +1,4 @@
-#! /bin/sh
+#! /bin/bash
# Script to rebuild and test a utility during development.
#
# During the development of Gnuastro, you often make changes in a utility
@@ -59,6 +59,7 @@
+
# SET INPUT PARAMETERS
# ====================
@@ -73,6 +74,7 @@ builddir=build
outdir=
+
# Set the utility name, along with its arguments and options. NOTE, for
# multiple arguments and options, please put them all between quotation
# marks so the space characters are included. If there are spaces in the
- [gnuastro-commits] master 5831e9e 052/113: Recent work in master imported, no conflicts, (continued)
- [gnuastro-commits] master 5831e9e 052/113: Recent work in master imported, no conflicts, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 1af4ea9 062/113: Correct checks for kernel dimension in NoiseChisel and Segment, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master d13c715 065/113: Imported work in master, no conflicts, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 91f2d3e 068/113: Imported recent work in master, conflicts fixed, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 63b4edd 070/113: Imported work in master, conflicts fixed, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 0ced9e5 072/113: Imported recent work in master, no conflicts, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 88935ae 073/113: Dataset pointers initialized to NULL in upperlimit_write_check, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master d02c999 079/113: Recent work in master imported, conflicts fixed, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 4021652 078/113: Imported recent work in master, conflicts fixed, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 8d64628 081/113: Recent work in master imported, minor conflicts fixed, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 800c769 082/113: Imported all the work in master, conflicts fixed,
Mohammad Akhlaghi <=
- [gnuastro-commits] master ac91655 084/113: Imported recent work in master, no conflicts, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master d57e0e6 087/113: New spectrum measurement from 3D inputs in MakeCatalog, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master a8428fe 088/113: Imported recent work in master, conflicts fixed, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 175354f 095/113: All MakeCatalog spectrums set to NaN when no measurement was done, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master fe0e594 097/113: Imported recent work in master, no conflicts, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 5d5fd95 106/113: NoiseChisel: independent test of input's dimensions, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master ee4bf1e 069/113: Imported recent work from master, Segment's 3D kernels added, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 2c6c7e0 076/113: Imported recent work, no conflicts, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master e2aba57 083/113: Recent work in master imported, conflicts fixed, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master e92fd9c 039/113: Merged recent work in master, conflicts fixed, Mohammad Akhlaghi, 2021/04/16