[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 8f3b4804: Library (statistics.h): MAD and MAD-
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 8f3b4804: Library (statistics.h): MAD and MAD-clipping functions added |
|
Date: |
Wed, 15 Nov 2023 18:45:16 -0500 (EST) |
branch: master
commit 8f3b48047ce53399f8eb67adafba40d6c001bf49
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Library (statistics.h): MAD and MAD-clipping functions added
Until now, the only type of available clipping of outliers was
sigma-clipping. However, the standard deviation (sigma) is itself the most
strongly affected statistic by an outlier! Therefore, sigma-clipping is not
too effective usually!
With this commit, the statistics library now has functions to find the
median absolute deviation (MAD), and from that, we now also have a
'gal_statistics_clip_mad' funciton. Having that function, many operators
have been added to the Arithmetic program for stacking and
collapsing. Within each a filled re-clipping operator has also been
implemeneted.
---
NEWS | 131 +++-
bin/arithmetic/arithmetic.c | 222 ++++++-
bin/arithmetic/arithmetic.h | 16 +
bin/mkcatalog/parse.c | 63 +-
bin/mkcatalog/upperlimit.c | 38 +-
bin/noisechisel/sky.c | 32 +-
bin/statistics/args.h | 125 ++++
bin/statistics/aststatistics.conf | 1 +
bin/statistics/main.h | 4 +-
bin/statistics/sky.c | 22 +-
bin/statistics/statistics.c | 186 +++++-
bin/statistics/ui.c | 27 +-
bin/statistics/ui.h | 9 +
doc/gnuastro.texi | 1301 ++++++++++++++++++++++++++++++-------
lib/arithmetic.c | 503 +++++++++++---
lib/binary.c | 29 +-
lib/dimension.c | 669 +++++++++++++++++--
lib/gnuastro/arithmetic.h | 17 +
lib/gnuastro/dimension.h | 100 +++
lib/gnuastro/statistics.h | 40 +-
lib/statistics.c | 518 +++++++++++----
21 files changed, 3408 insertions(+), 645 deletions(-)
diff --git a/NEWS b/NEWS
index 4363a4db..550674ca 100644
--- a/NEWS
+++ b/NEWS
@@ -5,6 +5,117 @@ See the end of the file for license conditions.
* Noteworthy changes in release X.XX (library XX.0.0) (YYYY-MM-DD)
** New features
+*** Arithmetic
+ - New operators:
+ - mad: Median Absolute Deviation (MAD) stacking.
+ - madclip-mad: MAD after MAD-clipping stacking.
+ - madclip-std: Standard deviation after MAD-clipping stacking.
+ - madclip-mean: Mean after MAD-clipping stacking.
+ - madclip-median: Median after MAD-clipping stacking.
+ - madclip-number: Number of elements after MAD-clipping stacking.
+ - madclip-fill-mad: MAD after filled MAD-clipping stacking: this
+ involves a two-phase clipping: after the first
+ clipping, Arithmetic will "fill" the unclipped
+ holes in each input and masks them for a second
+ round. This is critical for clipping more
+ diffuse outliers: where a pixel may not be
+ clipped individually and will produce biased
+ results, but due to proximity with many nearby
+ clipped pixels, it can be discarded and thus its
+ effect on the final stack be removed.
+ - madclip-fill-std: Standard deviation after filled MAD-clipping stacking.
+ - madclip-fill-mean: Mean after filled MAD-clipping stacking.
+ - madclip-fill-median: Median after filled MAD-clipping stacking.
+ - madclip-fill-number: Num. of elements after filled MAD-clipping stacking.
+ - sigclip-fill-mad: MAD after filled MAD-clipping stacking.
+ - sigclip-fill-std: Standard deviation after filled MAD-clipping stacking.
+ - sigclip-fill-mean: Mean after filled MAD-clipping stacking.
+ - sigclip-fill-median: Median after filled MAD-clipping stacking.
+ - sigclip-fill-number: Num. of elements after filled MAD-clipping stacking.
+ - sigclip-mad: MAD after sigma-clipping stacking.
+ - collapse-madclip-mad: Collapse dim. by MAD-clipped MAD.
+ - collapse-madclip-std: Collapse dim. by MAD-clipped STD.
+ - collapse-madclip-mean: Collapse dim. by MAD-clipped mean.
+ - collapse-madclip-median: Collapse dim. by MAD-clipped median.
+ - collapse-madclip-number: Collapse dim. by MAD-clipped number.
+ - collapse-madclip-fill-mad: Collapse dim. by filled MAD-clipped MAD.
+ - collapse-madclip-fill-std: Collapse dim. by filled MAD-clipped STD.
+ - collapse-madclip-fill-mean: Collapse dim. by filled MAD-clipped mean.
+ - collapse-madclip-fill-median: Collapse dim. by filled MAD-clipped median.
+ - collapse-madclip-fill-number: Collapse dim. by filled MAD-clipped num.
+ - collapse-sigclip-mad: Collapse dim. by sigma-clipped MAD.
+ - collapse-sigclip-fill-mad: Collapse dim. by filled sigma-clipped MAD.
+ - collapse-sigclip-fill-std: Collapse dim. by filled sigma-clipped STD.
+ - collapse-sigclip-fill-mean: Collapse dim. by filled sigma-clipped mean.
+ - collapse-sigclip-fill-median: Collapse dim. by filled sigma-clip. median.
+ - collapse-sigclip-fill-number: Collapse dim. by filled sigma-clipped num.
+
+*** Statistics
+ --mad: Median Absolute Deviation (MAD) of input dataset.
+ --madclip: MAD-clipping, showing intermediate results.
+ --sigclip-mad: MAD of input after sigma-clipping.
+ --madclip-number: Number of input elements after MAD-clipping.
+ --madclip-median: Median of input elements after MAD-clipping.
+ --madclip-mean: Mean of input elements after MAD-clipping.
+ --madclip-std: Standard deviation of input elements after MAD-clipping.
+ --madclip-mad: MAD of input elements after MAD-clipping.
+ --mclipparams: Parameters of MAD-clipping, similar to sigma-clipping. Note
+ that in a Gaussian distribution, MAD = 0.67 sigma, so if
+ you want something similar to 3-sigma (the default), you
+ should set the MAD multiple to 5 (the default).
+
+*** Library
+**** Functions
+ - gal_dimension_collapse_mclip_mad: MAD-clipped MAD.
+ - gal_dimension_collapse_mclip_fill_mad: filled MAD-clipped MAD.
+ - gal_dimension_collapse_mclip_std: MAD-clipped STD.
+ - gal_dimension_collapse_mclip_fill_std: filled MAD-clipped STD.
+ - gal_dimension_collapse_mclip_mean: MAD-clipped mean.
+ - gal_dimension_collapse_mclip_fill_mean: filled MAD-clipped mean
+ - gal_dimension_collapse_mclip_median: MAD-clipped median.
+ - gal_dimension_collapse_mclip_fill_median: filled MAD-clipped median.
+ - gal_dimension_collapse_mclip_number: MAD-clipped number.
+ - gal_dimension_collapse_mclip_fill_number: filled MAD-clipped number.
+ - gal_dimension_collapse_sclip_mad: sigma-clipped MAD
+ - gal_dimension_collapse_sclip_fill_mad: filled sigma-clipped MAD.
+ - gal_dimension_collapse_sclip_fill_std: filled sigma-clipped STD.
+ - gal_dimension_collapse_sclip_fill_mean: filled sigma-clipped mean
+ - gal_dimension_collapse_sclip_fill_median: filled sigma-clipped median.
+ - gal_dimension_collapse_sclip_fill_number: filled sigma-clipped number.
+ - gal_statistics_clip_mad: MAD clipping of given input.
+ - gal_statistics_mad: return median absolute deviation (MAD).
+ - gal_statistics_median_mad: return median and MAD.
+
+**** Macros
+ - Used by 'gal_arithmetic':
+ - GAL_ARITHMETIC_OP_SIGCLIP_MAD: Sigma-clipped MAD.
+ - GAL_ARITHMETIC_OP_SIGCLIP_MAD: Sigma-clipped STD.
+ - GAL_ARITHMETIC_OP_SIGCLIP_FILL_NUMBER: MAD-clipped num. of arrays.
+ - GAL_ARITHMETIC_OP_SIGCLIP_FILL_MEAN: MAD-clipped mean of arrays.
+ - GAL_ARITHMETIC_OP_SIGCLIP_FILL_MEDIAN: MAD-clipped median of arrays.
+ - GAL_ARITHMETIC_OP_SIGCLIP_FILL_STD: MAD-clipped STD of arrays.
+ - GAL_ARITHMETIC_OP_SIGCLIP_FILL_MAD: MAD-clipped STD of arrays.
+ - GAL_ARITHMETIC_OP_MADCLIP_NUMBER: MAD-clipped number of mult. arrays.
+ - GAL_ARITHMETIC_OP_MADCLIP_MEAN: MAD-clipped mean of multiple arrays.
+ - GAL_ARITHMETIC_OP_MADCLIP_MEDIAN: MAD-clipped median of mult. arrays.
+ - GAL_ARITHMETIC_OP_MADCLIP_STD: MAD-clipped STD of multiple arrays.
+ - GAL_ARITHMETIC_OP_MADCLIP_MAD: MAD-clipped STD of multiple arrays.
+ - GAL_ARITHMETIC_OP_MADCLIP_FILL_NUMBER: MAD-clipped num. of arrays.
+ - GAL_ARITHMETIC_OP_MADCLIP_FILL_MEAN: MAD-clipped mean of arrays.
+ - GAL_ARITHMETIC_OP_MADCLIP_FILL_MEDIAN: MAD-clipped median of arrays.
+ - GAL_ARITHMETIC_OP_MADCLIP_FILL_STD: MAD-clipped STD of arrays.
+ - GAL_ARITHMETIC_OP_MADCLIP_FILL_MAD: MAD-clipped STD of arrays.
+ - Used by 'gal_statistics_clip_sigma' and 'gal_statistics_clip_mad':
+ - GAL_STATISTICS_CLIP_OUTCOL_NUMBER_USED: index of final number in output.
+ - GAL_STATISTICS_CLIP_OUTCOL_MEAN: index of clipped mean in output.
+ - GAL_STATISTICS_CLIP_OUTCOL_STD: index of clipped STD in output.
+ - GAL_STATISTICS_CLIP_OUTCOL_MEDIAN: index of clipped median in output.
+ - GAL_STATISTICS_CLIP_OUTCOL_MAD: index of clipped MAD in output.
+ - GAL_STATISTICS_CLIP_OUTCOL_NUMBER_CLIPS: index of clipped number in
output.
+ - GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_MEAN: bit flag for measuring mean.
+ - GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_STD: bit flag for measuring STD.
+ - GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_MAD: bit flag for measuring MAD.
+
** Removed features
** Changed features
*** Arithmetic
@@ -17,6 +128,12 @@ See the end of the file for license conditions.
Sepideh Eskandarlou and implemented by Faezeh Bidjarchian after a poll
on Gnuastro's Matrix chat channel (#gnuastro:openastronomy.org).
+*** Library
+ - GAL_STATISTICS_CLIP_MAX_CONVERGE: new name for the old
+ 'GAL_STATISTICS_SIG_CLIP_MAX_CONVERGE'. Since we now also have
+ MAD-clipping.
+ - gal_statistics_clip_sigma: new name for 'gal_statistics_sigma_clip'.
+
** Bugs fixed
- bug #64852: astscript-zeropoint: crash when input is in 0-th HDU;
reported by 'danc' in https://savannah.gnu.org/support/?110952; fixed
@@ -131,20 +248,6 @@ See the end of the file for license conditions.
--river-min: minimum river value around a clump.
--river-max: minimum river value around a clump.
-*** Statistics
- --mad: Median Absolute Deviation (MAD) of input dataset.
- --madclip: MAD-clipping, showing intermediate results.
- --sigclip-mad: MAD of input after sigma-clipping.
- --madclip-number: Number of input elements after MAD-clipping.
- --madclip-median: Median of input elements after MAD-clipping.
- --madclip-mean: Mean of input elements after MAD-clipping.
- --madclip-std: Standard deviation of input elements after MAD-clipping.
- --madclip-mad: MAD of input elements after MAD-clipping.
- --mclipparams: Parameters of MAD-clipping, similar to sigma-clipping. Note
- that in a Gaussian distribution, MAD = 0.67 sigma, so if
- you want something similar to 3-sigma (the default), you
- should set the MAD multiple to 5 (the default).
-
*** Table
--info-num-cols: print the number of the input table's columns and abort.
--info-num-rows: print the number of the input table's rows and abort.
diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index 7bcfbd64..5f82fbd5 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -69,8 +69,8 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
if(a>0) return a; }
static size_t
-pop_number_of_operands(struct arithmeticparams *p, int op, char *token_string,
- gal_data_t **params)
+pop_number_of_operands(struct arithmeticparams *p, int op,
+ char *token_string, gal_data_t **params)
{
char *cstring="first";
size_t c, numparams=0;
@@ -83,10 +83,26 @@ pop_number_of_operands(struct arithmeticparams *p, int op,
char *token_string,
case GAL_ARITHMETIC_OP_QUANTILE:
numparams=1;
break;
+ case GAL_ARITHMETIC_OP_MADCLIP_STD:
case GAL_ARITHMETIC_OP_SIGCLIP_STD:
+ case GAL_ARITHMETIC_OP_MADCLIP_MAD:
+ case GAL_ARITHMETIC_OP_SIGCLIP_MAD:
+ case GAL_ARITHMETIC_OP_MADCLIP_MEAN:
case GAL_ARITHMETIC_OP_SIGCLIP_MEAN:
+ case GAL_ARITHMETIC_OP_MADCLIP_MEDIAN:
case GAL_ARITHMETIC_OP_SIGCLIP_MEDIAN:
+ case GAL_ARITHMETIC_OP_MADCLIP_NUMBER:
case GAL_ARITHMETIC_OP_SIGCLIP_NUMBER:
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_MAD:
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_MAD:
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_STD:
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_STD:
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_MEAN:
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_MEAN:
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_MEDIAN:
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_MEDIAN:
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_NUMBER:
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_NUMBER:
numparams=2;
break;
}
@@ -205,6 +221,7 @@ arithmetic_filter(void *in_prm)
gal_data_t *input=afp->input;
size_t sind=-1;
+ int clipflags=0;
size_t ind, index, one=1;
gal_data_t *sigclip, *result=NULL;
size_t *hpfsize=afp->hpfsize, *hnfsize=afp->hnfsize;
@@ -282,15 +299,21 @@ arithmetic_filter(void *in_prm)
case ARITHMETIC_OP_FILTER_SIGCLIP_MEAN:
case ARITHMETIC_OP_FILTER_SIGCLIP_MEDIAN:
- /* Find the sigma-clipped results. */
- sigclip=gal_statistics_sigma_clip(tile, afp->sclip_multip,
- afp->sclip_param, 0, 1);
+ /* The median is always available with a sigma-clip, but the mean
+ needs to be explicitly requested. */
+ if(afp->operator == ARITHMETIC_OP_FILTER_SIGCLIP_MEAN)
+ clipflags = GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_MEAN;
+ sigclip=gal_statistics_clip_sigma(tile, afp->sclip_multip,
+ afp->sclip_param, clipflags,
+ 0, 1);
/* Set the required index. */
switch(afp->operator)
{
- case ARITHMETIC_OP_FILTER_SIGCLIP_MEAN: sind = 2; break;
- case ARITHMETIC_OP_FILTER_SIGCLIP_MEDIAN: sind = 1; break;
+ case ARITHMETIC_OP_FILTER_SIGCLIP_MEAN:
+ sind = GAL_STATISTICS_CLIP_OUTCOL_MEAN; break;
+ case ARITHMETIC_OP_FILTER_SIGCLIP_MEDIAN:
+ sind = GAL_STATISTICS_CLIP_OUTCOL_MEDIAN; break;
default:
error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at "
"%s to fix the problem. The 'afp->operator' value "
@@ -917,8 +940,8 @@ arithmetic_interpolate(struct arithmeticparams *p, int
operator, char *token)
interpop=GAL_INTERPOLATE_NEIGHBORS_FUNC_MEDIAN;
break;
default:
- error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix "
- "the problem. The operator code %d isn't recognized",
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
+ "fix the problem. The operator code %d isn't recognized",
__func__, PACKAGE_BUGREPORT, operator);
}
@@ -943,15 +966,16 @@ arithmetic_collapse_single_value(gal_data_t *input, char
*counter,
int checkint)
{
if( input->ndim!=1 || input->size!=1)
- error(EXIT_FAILURE, 0, "%s popped operand of 'collapse-%s*' operators "
- "(%s) must be a single number (single-element, "
- "one-dimensional dataset). But it has %zu dimension(s) and %zu "
- "element(s).", counter, opstr, description, input->ndim,
input->size);
+ error(EXIT_FAILURE, 0, "%s popped operand of 'collapse-%s*' "
+ "operators (%s) must be a single number (single-element, "
+ "one-dimensional dataset). But it has %zu dimension(s) and "
+ "%zu element(s).", counter, opstr, description, input->ndim,
+ input->size);
if(checkint)
if(input->type==GAL_TYPE_FLOAT32 || input->type==GAL_TYPE_FLOAT64)
- error(EXIT_FAILURE, 0, "%s popped operand of 'collapse-%s*' operators "
- "(%s) must have an integer type, but it has a floating point "
- "type ('%s')", counter, opstr, description,
+ error(EXIT_FAILURE, 0, "%s popped operand of 'collapse-%s*' "
+ "operators (%s) must have an integer type, but it has a "
+ "floating point type ('%s')", counter, opstr, description,
gal_type_name(input->type,1));
}
@@ -976,10 +1000,26 @@ arithmetic_collapse(struct arithmeticparams *p, char
*token, int operator)
/* If we are on any of the sigma-clipping operators, we need to pop two
more operands. */
- if( operator==ARITHMETIC_OP_COLLAPSE_SIGCLIP_STD
+ if( operator==ARITHMETIC_OP_COLLAPSE_MADCLIP_MAD
+ || operator==ARITHMETIC_OP_COLLAPSE_MADCLIP_STD
+ || operator==ARITHMETIC_OP_COLLAPSE_MADCLIP_MEAN
+ || operator==ARITHMETIC_OP_COLLAPSE_MADCLIP_MEDIAN
+ || operator==ARITHMETIC_OP_COLLAPSE_MADCLIP_NUMBER
+ || operator==ARITHMETIC_OP_COLLAPSE_SIGCLIP_MAD
+ || operator==ARITHMETIC_OP_COLLAPSE_SIGCLIP_STD
|| operator==ARITHMETIC_OP_COLLAPSE_SIGCLIP_MEAN
|| operator==ARITHMETIC_OP_COLLAPSE_SIGCLIP_MEDIAN
- || operator==ARITHMETIC_OP_COLLAPSE_SIGCLIP_NUMBER )
+ || operator==ARITHMETIC_OP_COLLAPSE_SIGCLIP_NUMBER
+ || operator==ARITHMETIC_OP_COLLAPSE_MADCLIP_FILL_MAD
+ || operator==ARITHMETIC_OP_COLLAPSE_MADCLIP_FILL_STD
+ || operator==ARITHMETIC_OP_COLLAPSE_MADCLIP_FILL_MEAN
+ || operator==ARITHMETIC_OP_COLLAPSE_MADCLIP_FILL_MEDIAN
+ || operator==ARITHMETIC_OP_COLLAPSE_MADCLIP_FILL_NUMBER
+ || operator==ARITHMETIC_OP_COLLAPSE_SIGCLIP_FILL_MAD
+ || operator==ARITHMETIC_OP_COLLAPSE_SIGCLIP_FILL_STD
+ || operator==ARITHMETIC_OP_COLLAPSE_SIGCLIP_FILL_MEAN
+ || operator==ARITHMETIC_OP_COLLAPSE_SIGCLIP_FILL_MEDIAN
+ || operator==ARITHMETIC_OP_COLLAPSE_SIGCLIP_FILL_NUMBER )
{
param2=operands_pop(p, token); /* Termination criteria. */
param1=operands_pop(p, token); /* Multiple of sigma. */
@@ -996,9 +1036,9 @@ arithmetic_collapse(struct arithmeticparams *p, char
*token, int operator)
dimension=gal_data_copy_to_new_type_free(dimension, GAL_TYPE_LONG);
dim=((long *)(dimension->array))[0];
if(dim<=0)
- error(EXIT_FAILURE, 0, "first popped operand of 'collapse-*' operators "
- "(dimension to collapse) must be positive (larger than zero), it "
- "is %ld", dim);
+ error(EXIT_FAILURE, 0, "first popped operand of 'collapse-*' "
+ "operators (dimension to collapse) must be positive (larger "
+ "than zero), it is %ld", dim);
if(dim > input->ndim)
error(EXIT_FAILURE, 0, "input dataset to '%s' has %zu dimension(s), "
"but you have asked to collapse along dimension %ld", token,
@@ -1006,9 +1046,11 @@ arithmetic_collapse(struct arithmeticparams *p, char
*token, int operator)
if(param1)
{
arithmetic_collapse_single_value(param1, "third", "sigclip-",
- "sigma-clip's multiple of sigma", 0);
+ "sigma-clip's multiple of sigma",
+ 0);
arithmetic_collapse_single_value(param2, "second", "sigclip-",
- "sigma-clip's termination param", 0);
+ "sigma-clip's termination param",
+ 0);
param1=gal_data_copy_to_new_type_free(param1, GAL_TYPE_FLOAT32);
param2=gal_data_copy_to_new_type_free(param2, GAL_TYPE_FLOAT32);
p1=((float *)(param1->array))[0];
@@ -1057,30 +1099,110 @@ arithmetic_collapse(struct arithmeticparams *p, char
*token, int operator)
nt, mms, qmm);
break;
+ case ARITHMETIC_OP_COLLAPSE_MADCLIP_MAD:
+ collapsed=gal_dimension_collapse_mclip_std(input, input->ndim-dim,
+ p1, p2, nt, mms, qmm);
+ break;
+
+ case ARITHMETIC_OP_COLLAPSE_MADCLIP_FILL_MAD:
+ collapsed=gal_dimension_collapse_mclip_fill_std(input,
+ input->ndim-dim, p1, p2, nt, mms, qmm);
+ break;
+
+ case ARITHMETIC_OP_COLLAPSE_MADCLIP_STD:
+ collapsed=gal_dimension_collapse_mclip_std(input, input->ndim-dim,
+ p1, p2, nt, mms, qmm);
+ break;
+
+ case ARITHMETIC_OP_COLLAPSE_MADCLIP_FILL_STD:
+ collapsed=gal_dimension_collapse_mclip_fill_std(input,
+ input->ndim-dim, p1, p2, nt, mms, qmm);
+ break;
+
+ case ARITHMETIC_OP_COLLAPSE_MADCLIP_MEAN:
+ collapsed=gal_dimension_collapse_mclip_mean(input, input->ndim-dim,
+ p1, p2, nt, mms, qmm);
+ break;
+
+ case ARITHMETIC_OP_COLLAPSE_MADCLIP_FILL_MEAN:
+ collapsed=gal_dimension_collapse_mclip_fill_mean(input,
+ input->ndim-dim, p1, p2, nt, mms, qmm);
+ break;
+
+ case ARITHMETIC_OP_COLLAPSE_MADCLIP_MEDIAN:
+ collapsed=gal_dimension_collapse_mclip_median(input, input->ndim-dim,
+ p1, p2, nt, mms, qmm);
+ break;
+
+ case ARITHMETIC_OP_COLLAPSE_MADCLIP_FILL_MEDIAN:
+ collapsed=gal_dimension_collapse_mclip_fill_median(input,
+ input->ndim-dim, p1, p2, nt, mms, qmm);
+ break;
+
+ case ARITHMETIC_OP_COLLAPSE_MADCLIP_NUMBER:
+ collapsed=gal_dimension_collapse_mclip_number(input, input->ndim-dim,
+ p1, p2, nt, mms, qmm);
+ break;
+
+ case ARITHMETIC_OP_COLLAPSE_MADCLIP_FILL_NUMBER:
+ collapsed=gal_dimension_collapse_mclip_fill_number(input,
+ input->ndim-dim, p1, p2, nt, mms, qmm);
+ break;
+
+ case ARITHMETIC_OP_COLLAPSE_SIGCLIP_MAD:
+ collapsed=gal_dimension_collapse_sclip_std(input, input->ndim-dim,
+ p1, p2, nt, mms, qmm);
+ break;
+
+ case ARITHMETIC_OP_COLLAPSE_SIGCLIP_FILL_MAD:
+ collapsed=gal_dimension_collapse_sclip_fill_std(input,
+ input->ndim-dim, p1, p2, nt, mms, qmm);
+ break;
+
case ARITHMETIC_OP_COLLAPSE_SIGCLIP_STD:
collapsed=gal_dimension_collapse_sclip_std(input, input->ndim-dim,
p1, p2, nt, mms, qmm);
break;
+ case ARITHMETIC_OP_COLLAPSE_SIGCLIP_FILL_STD:
+ collapsed=gal_dimension_collapse_sclip_fill_std(input,
+ input->ndim-dim, p1, p2, nt, mms, qmm);
+ break;
+
case ARITHMETIC_OP_COLLAPSE_SIGCLIP_MEAN:
collapsed=gal_dimension_collapse_sclip_mean(input, input->ndim-dim,
p1, p2, nt, mms, qmm);
break;
+ case ARITHMETIC_OP_COLLAPSE_SIGCLIP_FILL_MEAN:
+ collapsed=gal_dimension_collapse_sclip_fill_mean(input,
+ input->ndim-dim, p1, p2, nt, mms, qmm);
+ break;
+
case ARITHMETIC_OP_COLLAPSE_SIGCLIP_MEDIAN:
collapsed=gal_dimension_collapse_sclip_median(input, input->ndim-dim,
p1, p2, nt, mms, qmm);
break;
+ case ARITHMETIC_OP_COLLAPSE_SIGCLIP_FILL_MEDIAN:
+ collapsed=gal_dimension_collapse_sclip_fill_median(input,
+ input->ndim-dim, p1, p2, nt, mms, qmm);
+ break;
+
case ARITHMETIC_OP_COLLAPSE_SIGCLIP_NUMBER:
collapsed=gal_dimension_collapse_sclip_number(input, input->ndim-dim,
p1, p2, nt, mms, qmm);
break;
+ case ARITHMETIC_OP_COLLAPSE_SIGCLIP_FILL_NUMBER:
+ collapsed=gal_dimension_collapse_sclip_fill_number(input,
+ input->ndim-dim, p1, p2, nt, mms, qmm);
+ break;
+
default:
- error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix the "
- "problem. The operator code %d is not recognized", __func__,
- PACKAGE_BUGREPORT, operator);
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
+ "fix the problem. The operator code %d is not recognized",
+ __func__, PACKAGE_BUGREPORT, operator);
}
@@ -1384,6 +1506,18 @@ arithmetic_set_operator(char *string, size_t
*num_operands, int *inlib)
{ op=ARITHMETIC_OP_COLLAPSE_MEDIAN; *num_operands=0; }
else if (!strcmp(string, "collapse-median"))
{ op=ARITHMETIC_OP_COLLAPSE_MEDIAN; *num_operands=0; }
+ else if (!strcmp(string, "collapse-madclip-mad"))
+ { op=ARITHMETIC_OP_COLLAPSE_MADCLIP_MAD; *num_operands=0; }
+ else if (!strcmp(string, "collapse-madclip-std"))
+ { op=ARITHMETIC_OP_COLLAPSE_MADCLIP_STD; *num_operands=0; }
+ else if (!strcmp(string, "collapse-madclip-mean"))
+ { op=ARITHMETIC_OP_COLLAPSE_MADCLIP_MEAN; *num_operands=0; }
+ else if (!strcmp(string, "collapse-madclip-median"))
+ { op=ARITHMETIC_OP_COLLAPSE_MADCLIP_MEDIAN;*num_operands=0;}
+ else if (!strcmp(string, "collapse-madclip-number"))
+ { op=ARITHMETIC_OP_COLLAPSE_MADCLIP_NUMBER;*num_operands=0;}
+ else if (!strcmp(string, "collapse-sigclip-mad"))
+ { op=ARITHMETIC_OP_COLLAPSE_SIGCLIP_MAD; *num_operands=0; }
else if (!strcmp(string, "collapse-sigclip-std"))
{ op=ARITHMETIC_OP_COLLAPSE_SIGCLIP_STD; *num_operands=0; }
else if (!strcmp(string, "collapse-sigclip-mean"))
@@ -1392,6 +1526,26 @@ arithmetic_set_operator(char *string, size_t
*num_operands, int *inlib)
{ op=ARITHMETIC_OP_COLLAPSE_SIGCLIP_MEDIAN;*num_operands=0;}
else if (!strcmp(string, "collapse-sigclip-number"))
{ op=ARITHMETIC_OP_COLLAPSE_SIGCLIP_NUMBER;*num_operands=0;}
+ else if (!strcmp(string, "collapse-madclip-fill-mad"))
+ { op=ARITHMETIC_OP_COLLAPSE_MADCLIP_FILL_MAD; *num_operands=0; }
+ else if (!strcmp(string, "collapse-madclip-fill-std"))
+ { op=ARITHMETIC_OP_COLLAPSE_MADCLIP_FILL_STD; *num_operands=0; }
+ else if (!strcmp(string, "collapse-madclip-fill-mean"))
+ { op=ARITHMETIC_OP_COLLAPSE_MADCLIP_FILL_MEAN; *num_operands=0; }
+ else if (!strcmp(string, "collapse-madclip-fill-median"))
+ { op=ARITHMETIC_OP_COLLAPSE_MADCLIP_FILL_MEDIAN;*num_operands=0;}
+ else if (!strcmp(string, "collapse-madclip-fill-number"))
+ { op=ARITHMETIC_OP_COLLAPSE_MADCLIP_FILL_NUMBER;*num_operands=0;}
+ else if (!strcmp(string, "collapse-sigclip-fill-mad"))
+ { op=ARITHMETIC_OP_COLLAPSE_SIGCLIP_FILL_MAD; *num_operands=0; }
+ else if (!strcmp(string, "collapse-sigclip-fill-std"))
+ { op=ARITHMETIC_OP_COLLAPSE_SIGCLIP_FILL_STD; *num_operands=0; }
+ else if (!strcmp(string, "collapse-sigclip-fill-mean"))
+ { op=ARITHMETIC_OP_COLLAPSE_SIGCLIP_FILL_MEAN; *num_operands=0; }
+ else if (!strcmp(string, "collapse-sigclip-fill-median"))
+ { op=ARITHMETIC_OP_COLLAPSE_SIGCLIP_FILL_MEDIAN;*num_operands=0;}
+ else if (!strcmp(string, "collapse-sigclip-fill-number"))
+ { op=ARITHMETIC_OP_COLLAPSE_SIGCLIP_FILL_NUMBER;*num_operands=0;}
else if (!strcmp(string, "add-dimension-slow"))
{ op=ARITHMETIC_OP_ADD_DIMENSION_SLOW; *num_operands=0; }
else if (!strcmp(string, "add-dimension-fast"))
@@ -1625,10 +1779,26 @@ arithmetic_operator_run(struct arithmeticparams *p, int
operator,
case ARITHMETIC_OP_COLLAPSE_MEAN:
case ARITHMETIC_OP_COLLAPSE_MEDIAN:
case ARITHMETIC_OP_COLLAPSE_NUMBER:
+ case ARITHMETIC_OP_COLLAPSE_MADCLIP_MAD:
+ case ARITHMETIC_OP_COLLAPSE_SIGCLIP_MAD:
+ case ARITHMETIC_OP_COLLAPSE_MADCLIP_STD:
case ARITHMETIC_OP_COLLAPSE_SIGCLIP_STD:
+ case ARITHMETIC_OP_COLLAPSE_MADCLIP_MEAN:
case ARITHMETIC_OP_COLLAPSE_SIGCLIP_MEAN:
+ case ARITHMETIC_OP_COLLAPSE_MADCLIP_MEDIAN:
case ARITHMETIC_OP_COLLAPSE_SIGCLIP_MEDIAN:
+ case ARITHMETIC_OP_COLLAPSE_MADCLIP_NUMBER:
case ARITHMETIC_OP_COLLAPSE_SIGCLIP_NUMBER:
+ case ARITHMETIC_OP_COLLAPSE_MADCLIP_FILL_MAD:
+ case ARITHMETIC_OP_COLLAPSE_SIGCLIP_FILL_MAD:
+ case ARITHMETIC_OP_COLLAPSE_MADCLIP_FILL_STD:
+ case ARITHMETIC_OP_COLLAPSE_SIGCLIP_FILL_STD:
+ case ARITHMETIC_OP_COLLAPSE_MADCLIP_FILL_MEAN:
+ case ARITHMETIC_OP_COLLAPSE_SIGCLIP_FILL_MEAN:
+ case ARITHMETIC_OP_COLLAPSE_MADCLIP_FILL_MEDIAN:
+ case ARITHMETIC_OP_COLLAPSE_SIGCLIP_FILL_MEDIAN:
+ case ARITHMETIC_OP_COLLAPSE_MADCLIP_FILL_NUMBER:
+ case ARITHMETIC_OP_COLLAPSE_SIGCLIP_FILL_NUMBER:
arithmetic_collapse(p, operator_string, operator);
break;
diff --git a/bin/arithmetic/arithmetic.h b/bin/arithmetic/arithmetic.h
index 7abdcfbf..d27ec402 100644
--- a/bin/arithmetic/arithmetic.h
+++ b/bin/arithmetic/arithmetic.h
@@ -53,10 +53,26 @@ enum arithmetic_prog_operators
ARITHMETIC_OP_COLLAPSE_MEAN,
ARITHMETIC_OP_COLLAPSE_MEDIAN,
ARITHMETIC_OP_COLLAPSE_NUMBER,
+ ARITHMETIC_OP_COLLAPSE_SIGCLIP_MAD,
ARITHMETIC_OP_COLLAPSE_SIGCLIP_STD,
ARITHMETIC_OP_COLLAPSE_SIGCLIP_MEAN,
ARITHMETIC_OP_COLLAPSE_SIGCLIP_MEDIAN,
ARITHMETIC_OP_COLLAPSE_SIGCLIP_NUMBER,
+ ARITHMETIC_OP_COLLAPSE_MADCLIP_MAD,
+ ARITHMETIC_OP_COLLAPSE_MADCLIP_STD,
+ ARITHMETIC_OP_COLLAPSE_MADCLIP_MEAN,
+ ARITHMETIC_OP_COLLAPSE_MADCLIP_MEDIAN,
+ ARITHMETIC_OP_COLLAPSE_MADCLIP_NUMBER,
+ ARITHMETIC_OP_COLLAPSE_SIGCLIP_FILL_MAD,
+ ARITHMETIC_OP_COLLAPSE_SIGCLIP_FILL_STD,
+ ARITHMETIC_OP_COLLAPSE_SIGCLIP_FILL_MEAN,
+ ARITHMETIC_OP_COLLAPSE_SIGCLIP_FILL_MEDIAN,
+ ARITHMETIC_OP_COLLAPSE_SIGCLIP_FILL_NUMBER,
+ ARITHMETIC_OP_COLLAPSE_MADCLIP_FILL_MAD,
+ ARITHMETIC_OP_COLLAPSE_MADCLIP_FILL_STD,
+ ARITHMETIC_OP_COLLAPSE_MADCLIP_FILL_MEAN,
+ ARITHMETIC_OP_COLLAPSE_MADCLIP_FILL_MEDIAN,
+ ARITHMETIC_OP_COLLAPSE_MADCLIP_FILL_NUMBER,
ARITHMETIC_OP_ADD_DIMENSION_SLOW,
ARITHMETIC_OP_ADD_DIMENSION_FAST,
ARITHMETIC_OP_REPEAT,
diff --git a/bin/mkcatalog/parse.c b/bin/mkcatalog/parse.c
index a7a6e0d0..eb53564c 100644
--- a/bin/mkcatalog/parse.c
+++ b/bin/mkcatalog/parse.c
@@ -1278,6 +1278,7 @@ parse_order_based(struct mkcatalog_passparams *pp)
double *ci;
float *sigcliparr;
gal_data_t *result;
+ uint8_t clipflags=0;
int32_t *O, *OO, *C=NULL;
size_t i, increment=0, num_increment=1;
gal_data_t *objvals=NULL, **clumpsvals=NULL;
@@ -1408,19 +1409,29 @@ parse_order_based(struct mkcatalog_passparams *pp)
|| p->oiflag[ OCOL_SIGCLIPMEAN ]
|| p->oiflag[ OCOL_SIGCLIPMEDIAN ])
{
- /* Calculate the sigma-clipped results and write them in any
- requested column. */
- result=gal_statistics_sigma_clip(objvals, p->sigmaclip[0],
- p->sigmaclip[1], 1, 1);
+ /* See which optional clipping measurements are necessary and run the
+ clipping. */
+ clipflags=0;
+ if(p->oiflag[ OCOL_SIGCLIPSTD ])
+ clipflags |= GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_STD;
+ if(p->oiflag[ OCOL_SIGCLIPMEAN ])
+ clipflags |= GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_MEAN;
+ result=gal_statistics_clip_sigma(objvals, p->sigmaclip[0],
+ p->sigmaclip[1], clipflags,
+ 1, 1);
sigcliparr=result->array;
if(p->oiflag[ OCOL_SIGCLIPNUM ])
- pp->oi[OCOL_SIGCLIPNUM]=sigcliparr[0];
+ pp->oi[OCOL_SIGCLIPNUM]
+ = sigcliparr[GAL_STATISTICS_CLIP_OUTCOL_NUMBER_USED];
if(p->oiflag[ OCOL_SIGCLIPSTD ])
- pp->oi[OCOL_SIGCLIPSTD]=sigcliparr[3];
+ pp->oi[OCOL_SIGCLIPSTD]
+ = sigcliparr[GAL_STATISTICS_CLIP_OUTCOL_STD];
if(p->oiflag[ OCOL_SIGCLIPMEAN ])
- pp->oi[OCOL_SIGCLIPMEAN]=sigcliparr[2];
+ pp->oi[OCOL_SIGCLIPMEAN]
+ = sigcliparr[GAL_STATISTICS_CLIP_OUTCOL_MEAN];
if(p->oiflag[ OCOL_SIGCLIPMEDIAN ])
- pp->oi[OCOL_SIGCLIPMEDIAN]=sigcliparr[1];
+ pp->oi[OCOL_SIGCLIPMEDIAN]
+ = sigcliparr[GAL_STATISTICS_CLIP_OUTCOL_MEDIAN];
/* Clean up the sigma-clipped values. */
gal_data_free(result);
@@ -1444,7 +1455,7 @@ parse_order_based(struct mkcatalog_passparams *pp)
{
for(i=0;i<pp->clumpsinobj;++i)
{
- /* Set the main row to fill. */
+ /* Set the main row to fill and initialize. */
ci=&pp->ci[ i * CCOL_NUMCOLS ];
/* Median. */
@@ -1471,24 +1482,34 @@ parse_order_based(struct mkcatalog_passparams *pp)
{
if(clumpsvals[i])
{
- result=gal_statistics_sigma_clip(clumpsvals[i],
+ clipflags=0;
+ if(p->oiflag[ OCOL_SIGCLIPSTD ])
+ clipflags |= GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_STD;
+ if(p->oiflag[ OCOL_SIGCLIPMEAN ])
+ clipflags |= GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_MEAN;
+ result=gal_statistics_clip_sigma(clumpsvals[i],
p->sigmaclip[0],
- p->sigmaclip[1], 1, 1);
+ p->sigmaclip[1],
+ clipflags, 1, 1);
sigcliparr=result->array;
if(p->ciflag[ CCOL_SIGCLIPNUM ])
- ci[CCOL_SIGCLIPNUM]=sigcliparr[0];
+ ci[CCOL_SIGCLIPNUM]
+ = sigcliparr[GAL_STATISTICS_CLIP_OUTCOL_NUMBER_USED];
if(p->ciflag[ CCOL_SIGCLIPSTD ])
- ci[CCOL_SIGCLIPSTD]=( sigcliparr[3]
- - ( ci[ CCOL_RIV_SUM ]
- / ci[ CCOL_RIV_NUM ]));
+ ci[CCOL_SIGCLIPSTD]
+ = ( sigcliparr[GAL_STATISTICS_CLIP_OUTCOL_STD]
+ - ( ci[ CCOL_RIV_SUM ]
+ / ci[ CCOL_RIV_NUM ]));
if(p->ciflag[ CCOL_SIGCLIPMEAN ])
- ci[CCOL_SIGCLIPMEAN]=( sigcliparr[2]
- - ( ci[ CCOL_RIV_SUM ]
- / ci[ CCOL_RIV_NUM ]));
+ ci[CCOL_SIGCLIPMEAN]
+ = ( sigcliparr[GAL_STATISTICS_CLIP_OUTCOL_MEAN]
+ - ( ci[ CCOL_RIV_SUM ]
+ / ci[ CCOL_RIV_NUM ]));
if(p->ciflag[ CCOL_SIGCLIPMEDIAN ])
- ci[CCOL_SIGCLIPMEDIAN]=( sigcliparr[1]
- - ( ci[ CCOL_RIV_SUM ]
- / ci[ CCOL_RIV_NUM ]));
+ ci[CCOL_SIGCLIPMEDIAN]
+ = ( sigcliparr[GAL_STATISTICS_CLIP_OUTCOL_MEDIAN]
+ - ( ci[ CCOL_RIV_SUM ]
+ / ci[ CCOL_RIV_NUM ]));
gal_data_free(result);
}
else
diff --git a/bin/mkcatalog/upperlimit.c b/bin/mkcatalog/upperlimit.c
index 0470e162..ec7ddd16 100644
--- a/bin/mkcatalog/upperlimit.c
+++ b/bin/mkcatalog/upperlimit.c
@@ -448,14 +448,17 @@ static void
upperlimit_measure(struct mkcatalog_passparams *pp, int32_t clumplab,
int do_measurement)
{
- float *scarr;
gal_data_t *column;
+ float *mcarr, mcstd;
size_t init_size, col, one=1;
struct mkcatalogparams *p=pp->p;
- gal_data_t *sum, *qfunc=NULL, *sigclip=NULL;
+ gal_data_t *sum, *qfunc=NULL, *madclip=NULL;
double *o = ( clumplab
? &pp->ci[ (clumplab-1) * CCOL_NUMCOLS ]
: pp->oi );
+ uint8_t clipflags = ( GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_MEAN \
+ | GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_STD \
+ | GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_MAD );
/* If the random distribution exsits, then fill it in. */
if(do_measurement)
@@ -476,10 +479,11 @@ upperlimit_measure(struct mkcatalog_passparams *pp,
int32_t clumplab,
/* Similar to the case for sigma-clipping, we'll need to
keep the size here also. */
init_size=pp->up_vals->size;
- sum=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, 1, &one, NULL, 0,
- -1, 1, NULL, NULL, NULL);
+ sum=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, 1, &one, NULL,
+ 0, -1, 1, NULL, NULL, NULL);
((float *)(sum->array))[0]=o[clumplab?CCOL_SUM:OCOL_SUM];
- qfunc=gal_statistics_quantile_function(pp->up_vals, sum, 1);
+ qfunc=gal_statistics_quantile_function(pp->up_vals,
+ sum, 1);
/* Fill in the column. */
col = clumplab ? CCOL_UPPERLIMIT_Q : OCOL_UPPERLIMIT_Q;
@@ -496,37 +500,41 @@ upperlimit_measure(struct mkcatalog_passparams *pp,
int32_t clumplab,
default:
/* We only need to do this once, but the columns can be
requested in any order. */
- if(sigclip==NULL)
+ if(madclip==NULL)
{
/* Calculate the sigma-clipped standard deviation. Since
it is done in place, the size will change, so we'll
keep the size here and put it back after we are
done. */
init_size=pp->up_vals->size;
- sigclip=gal_statistics_sigma_clip(pp->up_vals,
+ madclip=gal_statistics_clip_sigma(pp->up_vals,
p->upsigmaclip[0],
- p->upsigmaclip[1],1,1);
+ p->upsigmaclip[1],
+ clipflags, 1, 1);
pp->up_vals->size=pp->up_vals->dsize[0]=init_size;
- scarr=sigclip->array;
+ mcarr=madclip->array;
/* 1-sigma. */
+ mcstd=mcarr[GAL_STATISTICS_CLIP_OUTCOL_STD];
col = clumplab ? CCOL_UPPERLIMIT_S : OCOL_UPPERLIMIT_S;
- o[col] = scarr[3];
+ o[col] = mcstd;
/* sigma multiplied by 'upnsigma'. */
col = clumplab ? CCOL_UPPERLIMIT_B : OCOL_UPPERLIMIT_B;
- o[col] = scarr[3] * p->upnsigma;
+ o[col] = mcstd * p->upnsigma;
/* Nonparametric skewness [ (Mean-Median)/STD ]. */
col = clumplab?CCOL_UPPERLIMIT_SKEW:OCOL_UPPERLIMIT_SKEW;
- o[col] = ( scarr[2] - scarr[1] ) / scarr[3];
-
- /* Clean up. */
- gal_data_free(sigclip);
+ o[col] = ( ( mcarr[GAL_STATISTICS_CLIP_OUTCOL_MEAN]
+ - mcarr[GAL_STATISTICS_CLIP_OUTCOL_MEDIAN] )
+ / mcstd );
}
break;
}
}
+
+ /* Clean up. */
+ gal_data_free(madclip);
}
else
{
diff --git a/bin/noisechisel/sky.c b/bin/noisechisel/sky.c
index e0383d17..631f8c72 100644
--- a/bin/noisechisel/sky.c
+++ b/bin/noisechisel/sky.c
@@ -61,16 +61,19 @@ sky_mean_std_undetected(void *in_prm)
int setblank, type=GAL_TYPE_FLOAT32;
uint8_t *noskytiles=p->noskytiles->array;
size_t i, tind, numsky, bdsize=2, ndim=p->sky->ndim;
+ gal_data_t *tile, *fusage, *busage, *bintile, *clip;
size_t refarea, twidth=gal_type_sizeof(GAL_TYPE_FLOAT32);
- gal_data_t *tile, *fusage, *busage, *bintile, *sigmaclip;
-
+ uint8_t mclipflags = ( GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_MEAN
+ | GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_STD );
/* 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, p->cp.quietmmap, NULL, NULL, NULL);
+ p->cp.minmapsize, p->cp.quietmmap, NULL, NULL,
+ NULL);
busage=gal_data_alloc(NULL, GAL_TYPE_UINT8, ndim, p->maxtsize, NULL, 0,
- p->cp.minmapsize, p->cp.quietmmap, NULL, NULL, NULL);
+ p->cp.minmapsize, p->cp.quietmmap, NULL, NULL,
+ NULL);
/* An empty dataset to replicate a tile on the binary array. */
@@ -129,9 +132,10 @@ sky_mean_std_undetected(void *in_prm)
gal_blank_flag_apply(fusage, busage);
- /* Do the sigma-clipping. */
- sigmaclip=gal_statistics_sigma_clip(fusage, p->sigmaclip[0],
- p->sigmaclip[1], 1, 1);
+ /* Do the MAD-clipping. */
+ clip=gal_statistics_clip_sigma(fusage, p->sigmaclip[0],
+ p->sigmaclip[1], mclipflags,
+ 1, 1);
/* When there are zero-valued pixels on the edges of the
@@ -140,25 +144,27 @@ sky_mean_std_undetected(void *in_prm)
binary value of 1 and so the Sky and its standard
deviation can become zero. So, we need ignore such
tiles. */
- if( ((float *)(sigmaclip->array))[3]==0.0 )
+ if( ((float *)(clip->array))[3]==0.0 )
setblank=1;
else
{
/* 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
+ sure 'clip' has the same type as the sky and std
arrays. */
- sigmaclip=gal_data_copy_to_new_type_free(sigmaclip, type);
+ clip=gal_data_copy_to_new_type_free(clip, type);
memcpy(gal_pointer_increment(p->sky->array, tind, type),
- gal_pointer_increment(sigmaclip->array, 2, type),
+ gal_pointer_increment(clip->array,
+ GAL_STATISTICS_CLIP_OUTCOL_MEAN, type),
twidth);
memcpy(gal_pointer_increment(p->std->array, tind, type),
- gal_pointer_increment(sigmaclip->array, 3, type),
+ gal_pointer_increment(clip->array,
+ GAL_STATISTICS_CLIP_OUTCOL_STD, type),
twidth);
}
/* Clean up. */
- gal_data_free(sigmaclip);
+ gal_data_free(clip);
}
else
setblank=1;
diff --git a/bin/statistics/args.h b/bin/statistics/args.h
index 135bea4b..2f8aa814 100644
--- a/bin/statistics/args.h
+++ b/bin/statistics/args.h
@@ -196,6 +196,20 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET,
ui_add_to_single_value
},
+ {
+ "mad",
+ UI_KEY_MAD,
+ 0,
+ 0,
+ "Median absolute deviation.",
+ UI_GROUP_SINGLE_VALUE,
+ &p->singlevalue,
+ GAL_OPTIONS_NO_ARG_TYPE,
+ GAL_OPTIONS_RANGE_0_OR_1,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ ui_add_to_single_value
+ },
{
"median",
UI_KEY_MEDIAN,
@@ -364,6 +378,90 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET,
ui_add_to_single_value
},
+ {
+ "sigclip-mad",
+ UI_KEY_SIGCLIPMAD,
+ 0,
+ 0,
+ "Sigma-clipped Median Absolute Deviation.",
+ UI_GROUP_SINGLE_VALUE,
+ &p->singlevalue,
+ GAL_OPTIONS_NO_ARG_TYPE,
+ GAL_OPTIONS_RANGE_0_OR_1,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ ui_add_to_single_value
+ },
+ {
+ "madclip-number",
+ UI_KEY_MADCLIPNUMBER,
+ 0,
+ 0,
+ "Number of elements after MAD-clipping.",
+ UI_GROUP_SINGLE_VALUE,
+ &p->singlevalue,
+ GAL_OPTIONS_NO_ARG_TYPE,
+ GAL_OPTIONS_RANGE_0_OR_1,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ ui_add_to_single_value
+ },
+ {
+ "madclip-median",
+ UI_KEY_MADCLIPMEDIAN,
+ 0,
+ 0,
+ "MAD-clipped median.",
+ UI_GROUP_SINGLE_VALUE,
+ &p->singlevalue,
+ GAL_OPTIONS_NO_ARG_TYPE,
+ GAL_OPTIONS_RANGE_0_OR_1,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ ui_add_to_single_value
+ },
+ {
+ "madclip-mean",
+ UI_KEY_MADCLIPMEAN,
+ 0,
+ 0,
+ "MAD-clipped mean.",
+ UI_GROUP_SINGLE_VALUE,
+ &p->singlevalue,
+ GAL_OPTIONS_NO_ARG_TYPE,
+ GAL_OPTIONS_RANGE_0_OR_1,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ ui_add_to_single_value
+ },
+ {
+ "madclip-std",
+ UI_KEY_MADCLIPSTD,
+ 0,
+ 0,
+ "MAD-clipped standard deviation.",
+ UI_GROUP_SINGLE_VALUE,
+ &p->singlevalue,
+ GAL_OPTIONS_NO_ARG_TYPE,
+ GAL_OPTIONS_RANGE_0_OR_1,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ ui_add_to_single_value
+ },
+ {
+ "madclip-mad",
+ UI_KEY_MADCLIPMAD,
+ 0,
+ 0,
+ "MAD-clipped Median Absolute Deviation.",
+ UI_GROUP_SINGLE_VALUE,
+ &p->singlevalue,
+ GAL_OPTIONS_NO_ARG_TYPE,
+ GAL_OPTIONS_RANGE_0_OR_1,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ ui_add_to_single_value
+ },
@@ -491,6 +589,19 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_MANDATORY,
GAL_OPTIONS_NOT_SET
},
+ {
+ "madclip",
+ UI_KEY_MADCLIP,
+ 0,
+ 0,
+ "Overall MAD-clipping (see '--mclipparams')",
+ UI_GROUP_PARTICULAR_STAT,
+ &p->madclip,
+ GAL_OPTIONS_NO_ARG_TYPE,
+ GAL_OPTIONS_RANGE_0_OR_1,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
{
"contour",
UI_KEY_CONTOUR,
@@ -622,6 +733,20 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET,
gal_options_read_sigma_clip
},
+ {
+ "mclipparams",
+ UI_KEY_MCLIPPARAMS,
+ "FLT,FLT",
+ 0,
+ "MAD clip: Multiple, and tolerance/number.",
+ UI_GROUP_SKY,
+ p->mclipparams,
+ GAL_TYPE_STRING,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ gal_options_read_sigma_clip
+ },
{
"smoothwidth",
UI_KEY_SMOOTHWIDTH,
diff --git a/bin/statistics/aststatistics.conf
b/bin/statistics/aststatistics.conf
index 8d0e807a..49f44988 100644
--- a/bin/statistics/aststatistics.conf
+++ b/bin/statistics/aststatistics.conf
@@ -29,6 +29,7 @@
outliersclip 3,0.2
smoothwidth 3
sclipparams 3,0.1
+ mclipparams 4.48,0.01
# Histogram and CFP settings
numasciibins 70
diff --git a/bin/statistics/main.h b/bin/statistics/main.h
index 0edcc00d..ea2409a5 100644
--- a/bin/statistics/main.h
+++ b/bin/statistics/main.h
@@ -91,7 +91,8 @@ struct statisticsparams
uint8_t cumulative; /* Save cumulative distibution in output. */
double mirror; /* Mirror value for hist and CFP. */
uint8_t sky; /* Find the Sky value over the image. */
- uint8_t sigmaclip; /* So sigma-clipping over all dataset. */
+ uint8_t sigmaclip; /* Do sigma-clipping over all dataset. */
+ uint8_t madclip; /* Do MAD-clipping over all dataset. */
gal_data_t *contour; /* Levels to show contours. */
size_t numbins; /* Number of bins in histogram or CFP. */
@@ -114,6 +115,7 @@ struct statisticsparams
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. */
+ double mclipparams[2]; /* Muliple and parameter of sigma clipping. */
uint8_t ignoreblankintiles;/* Ignore input's blank values. */
diff --git a/bin/statistics/sky.c b/bin/statistics/sky.c
index cc2ac1f0..af83f47a 100644
--- a/bin/statistics/sky.c
+++ b/bin/statistics/sky.c
@@ -55,9 +55,10 @@ sky_on_thread(void *in_prm)
void *tblock=NULL, *tarray=NULL;
int stype=p->sky_t->type, itype=p->input->type;
- gal_data_t *num, *tile, *mean, *meanquant, *sigmaclip;
+ gal_data_t *num, *tile, *mean, *meanquant, *clip;
size_t i, tind, twidth=gal_type_sizeof(p->sky_t->type);
-
+ uint8_t mclipflags = ( GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_MEAN
+ | GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_STD );
/* Find the Sky and its standard deviation on the tiles given to this
thread. */
@@ -97,19 +98,24 @@ sky_on_thread(void *in_prm)
/* Get the sigma-clipped mean and standard deviation. 'inplace'
is irrelevant here because this is a tile and it will be
copied anyway. */
- sigmaclip=gal_statistics_sigma_clip(tile, p->sclipparams[0],
- p->sclipparams[1], 1, 1);
+ clip=gal_statistics_clip_sigma(tile, p->sclipparams[0],
+ p->sclipparams[1], mclipflags,
+ 1, 1);
/* Put the mean and its standard deviation into the respective
place for this tile. */
- sigmaclip=gal_data_copy_to_new_type_free(sigmaclip, stype);
+ clip=gal_data_copy_to_new_type_free(clip, stype);
memcpy(gal_pointer_increment(p->sky_t->array, tind, stype),
- gal_pointer_increment(sigmaclip->array, 2, stype), twidth);
+ gal_pointer_increment(clip->array,
+ GAL_STATISTICS_CLIP_OUTCOL_MEAN,
+ stype), twidth);
memcpy(gal_pointer_increment(p->std_t->array, tind, stype),
- gal_pointer_increment(sigmaclip->array, 3, stype), twidth);
+ gal_pointer_increment(clip->array,
+ GAL_STATISTICS_CLIP_OUTCOL_STD,
+ stype), twidth);
/* Clean up. */
- gal_data_free(sigmaclip);
+ gal_data_free(clip);
}
else
{
diff --git a/bin/statistics/statistics.c b/bin/statistics/statistics.c
index 27f4454d..809e369f 100644
--- a/bin/statistics/statistics.c
+++ b/bin/statistics/statistics.c
@@ -92,16 +92,43 @@ statistics_read_check_args(struct statisticsparams *p)
+/* This function checks if any of the optional clipping measurements are
+ requested (among all the single-valued measuremens) or not. If so, it
+ will add their respective flag. */
+static uint8_t
+statistics_one_row_clip_flags(struct statisticsparams *p, int32_t mean,
+ int32_t std, int32_t mad)
+{
+ uint8_t flags=0;
+ gal_list_i32_t *tmp;
+
+ for(tmp=p->singlevalue; tmp!=NULL; tmp=tmp->next)
+ {
+ if(tmp->v==std) flags |= GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_STD;
+ if(tmp->v==mad) flags |= GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_MAD;
+ if(tmp->v==mean) flags |= GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_MEAN;
+ }
+
+ return flags;
+}
+
+
+
+
+
static void
statistics_print_one_row(struct statisticsparams *p)
{
int mustfree;
char *toprint;
double arg, *d;
+ uint8_t clipflag;
gal_list_i32_t *tmp;
+ gal_data_t *medmad=NULL;
size_t dsize=1, counter;
- gal_data_t *sum=NULL, *med=NULL, *meanstd=NULL, *modearr=NULL;
- gal_data_t *tmpv, *sclip=NULL, *out=NULL, *num=NULL, *min=NULL, *max=NULL;
+ gal_data_t *sclip=NULL, *mclip=NULL;
+ gal_data_t *sum=NULL, *meanstd=NULL, *modearr=NULL;
+ gal_data_t *tmpv, *out=NULL, *num=NULL, *min=NULL, *max=NULL;
/* The user can ask for any of the operators more than once, also some
operators might return more than one usable value (like mode). So we
@@ -121,13 +148,15 @@ statistics_print_one_row(struct statisticsparams *p)
max = max ? max : gal_statistics_maximum(p->input); break;
case UI_KEY_SUM:
sum = sum ? sum : gal_statistics_sum(p->input); break;
- case UI_KEY_MEDIAN:
- med = med ? med : gal_statistics_median(p->sorted, 0); break;
case UI_KEY_STD:
case UI_KEY_MEAN:
case UI_KEY_QUANTOFMEAN:
meanstd = meanstd ? meanstd : gal_statistics_mean_std(p->input);
break;
+ case UI_KEY_MAD:
+ case UI_KEY_MEDIAN:
+ medmad = medmad ? medmad : gal_statistics_median_mad(p->input, 0);
+ break;
case UI_KEY_MODE:
case UI_KEY_MODEQUANT:
case UI_KEY_MODESYM:
@@ -138,17 +167,38 @@ statistics_print_one_row(struct statisticsparams *p)
d=modearr->array;
if(d[2]<GAL_STATISTICS_MODE_GOOD_SYM) d[0]=d[1]=NAN;
break;
+ case UI_KEY_MADCLIPSTD:
+ case UI_KEY_MADCLIPMAD:
+ case UI_KEY_MADCLIPMEAN:
+ case UI_KEY_MADCLIPNUMBER:
+ case UI_KEY_MADCLIPMEDIAN:
+ if(mclip==NULL)
+ {
+ clipflag=statistics_one_row_clip_flags(p, UI_KEY_MADCLIPMEAN,
+ UI_KEY_MADCLIPSTD,
+ UI_KEY_MADCLIPMAD);
+ mclip=gal_statistics_clip_mad(p->sorted, p->mclipparams[0],
+ p->mclipparams[1], clipflag,
+ 0, 1);
+ }
+ break;
case UI_KEY_SIGCLIPSTD:
+ case UI_KEY_SIGCLIPMAD:
case UI_KEY_SIGCLIPMEAN:
case UI_KEY_SIGCLIPNUMBER:
case UI_KEY_SIGCLIPMEDIAN:
- sclip = ( sclip
- ? sclip
- : gal_statistics_sigma_clip(p->sorted, p->sclipparams[0],
- p->sclipparams[1], 0, 1) );
+ if(sclip==NULL)
+ {
+ clipflag=statistics_one_row_clip_flags(p, UI_KEY_SIGCLIPMEAN,
+ UI_KEY_SIGCLIPSTD,
+ UI_KEY_SIGCLIPMAD);
+ sclip=gal_statistics_clip_sigma(p->sorted, p->sclipparams[0],
+ p->sclipparams[1], clipflag,
+ 0, 1);
+ }
break;
- /* Will be calculated as printed. */
+ /* These will be calculated as printed. */
case UI_KEY_QUANTILE:
case UI_KEY_QUANTFUNC:
break;
@@ -176,7 +226,10 @@ statistics_print_one_row(struct statisticsparams *p)
case UI_KEY_MINIMUM: out=min; break;
case UI_KEY_MAXIMUM: out=max; break;
case UI_KEY_SUM: out=sum; break;
- case UI_KEY_MEDIAN: out=med; break;
+ case UI_KEY_MEDIAN:
+ out=statistics_pull_out_element(medmad, 0); mustfree=1; break;
+ case UI_KEY_MAD:
+ out=statistics_pull_out_element(medmad, 1); mustfree=1; break;
case UI_KEY_MEAN:
out=statistics_pull_out_element(meanstd, 0); mustfree=1; break;
case UI_KEY_STD:
@@ -189,14 +242,46 @@ statistics_print_one_row(struct statisticsparams *p)
out=statistics_pull_out_element(modearr, 2); mustfree=1; break;
case UI_KEY_MODESYMVALUE:
out=statistics_pull_out_element(modearr, 3); mustfree=1; break;
+ case UI_KEY_MADCLIPMAD:
+ out=statistics_pull_out_element(mclip,
+ GAL_STATISTICS_CLIP_OUTCOL_MAD);
+ mustfree=1; break;
+ case UI_KEY_MADCLIPSTD:
+ out=statistics_pull_out_element(mclip,
+ GAL_STATISTICS_CLIP_OUTCOL_STD);
+ mustfree=1; break;
+ case UI_KEY_MADCLIPMEAN:
+ out=statistics_pull_out_element(mclip,
+ GAL_STATISTICS_CLIP_OUTCOL_MEAN);
+ mustfree=1; break;
+ case UI_KEY_MADCLIPMEDIAN:
+ out=statistics_pull_out_element(mclip,
+ GAL_STATISTICS_CLIP_OUTCOL_MEDIAN);
+ mustfree=1; break;
+ case UI_KEY_MADCLIPNUMBER:
+ out=statistics_pull_out_element(mclip,
+ GAL_STATISTICS_CLIP_OUTCOL_NUMBER_USED);
+ mustfree=1; break;
+ case UI_KEY_SIGCLIPMAD:
+ out=statistics_pull_out_element(sclip,
+ GAL_STATISTICS_CLIP_OUTCOL_MAD);
+ mustfree=1; break;
case UI_KEY_SIGCLIPSTD:
- out=statistics_pull_out_element(sclip, 3); mustfree=1; break;
+ out=statistics_pull_out_element(sclip,
+ GAL_STATISTICS_CLIP_OUTCOL_STD);
+ mustfree=1; break;
case UI_KEY_SIGCLIPMEAN:
- out=statistics_pull_out_element(sclip, 2); mustfree=1; break;
+ out=statistics_pull_out_element(sclip,
+ GAL_STATISTICS_CLIP_OUTCOL_MEAN);
+ mustfree=1; break;
case UI_KEY_SIGCLIPMEDIAN:
- out=statistics_pull_out_element(sclip, 1); mustfree=1; break;
+ out=statistics_pull_out_element(sclip,
+ GAL_STATISTICS_CLIP_OUTCOL_MEDIAN);
+ mustfree=1; break;
case UI_KEY_SIGCLIPNUMBER:
- out=statistics_pull_out_element(sclip, 0); mustfree=1; break;
+ out=statistics_pull_out_element(sclip,
+ GAL_STATISTICS_CLIP_OUTCOL_NUMBER_USED);
+ mustfree=1; break;
/* Not previously calculated. */
case UI_KEY_QUANTILE:
@@ -253,8 +338,8 @@ statistics_print_one_row(struct statisticsparams *p)
if(min) gal_data_free(min);
if(max) gal_data_free(max);
if(sum) gal_data_free(sum);
- if(med) gal_data_free(med);
if(sclip) gal_data_free(sclip);
+ if(medmad) gal_data_free(medmad);
if(meanstd) gal_data_free(meanstd);
if(modearr) gal_data_free(modearr);
}
@@ -1105,23 +1190,32 @@ print_basics(struct statisticsparams *p)
/************** Sigma clipping ***************/
/*******************************************************************/
void
-print_sigma_clip(struct statisticsparams *p)
+print_clip(struct statisticsparams *p, int sig1_mad0)
{
float *a;
char *mode;
int namewidth=40;
- gal_data_t *sigclip;
+ gal_data_t *result;
+ double clipparams[2];
+ uint8_t flags = ( GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_MEAN
+ | GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_STD
+ | GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_MAD );
+
+ /* Fill in the clipparams. */
+ clipparams[0] = sig1_mad0 ? p->sclipparams[0] : p->mclipparams[0];
+ clipparams[1] = sig1_mad0 ? p->sclipparams[1] : p->mclipparams[1];
/* Set the mode for printing: */
- if( p->sclipparams[1]>=1.0f )
+ if( clipparams[1]>=1.0f )
{
- if( asprintf(&mode, "for %g clips", p->sclipparams[1])<0 )
+ if( asprintf(&mode, "for %g clips", clipparams[1])<0 )
error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
}
else
{
- if( asprintf(&mode, "until relative change in STD is less than %g",
- p->sclipparams[1])<0 )
+ if( asprintf(&mode, "until relative change in %s is less than %g",
+ sig1_mad0 ? "STD" : "MAD (median absolute deviation)",
+ clipparams[1])<0 )
error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
}
@@ -1129,31 +1223,46 @@ print_sigma_clip(struct statisticsparams *p)
if(!p->cp.quiet)
{
print_input_info(p);
- printf("%g-sigma clipping steps %s:\n\n", p->sclipparams[0], mode);
+ printf("%g-%s clipping steps %s:\n\n", clipparams[0],
+ sig1_mad0?"sigma":"MAD", mode);
}
- /* Do the Sigma clipping: */
- sigclip=gal_statistics_sigma_clip(p->sorted, p->sclipparams[0],
- p->sclipparams[1], 0, p->cp.quiet);
- a=sigclip->array;
+ /* Do the Sigma clipping. In-place is not activated because
+ sigma-clipping reduces the total number of points and users may
+ call other operators also. */
+ result = ( sig1_mad0
+ ? gal_statistics_clip_sigma(p->sorted, clipparams[0],
+ clipparams[1], flags, 0, p->cp.quiet)
+ : gal_statistics_clip_mad(p->sorted, clipparams[0],
+ clipparams[1], flags, 0, p->cp.quiet) );
+ a=result->array;
/* Finish the introduction. */
if(!p->cp.quiet)
- printf("-------\nSummary:\n");
+ printf("-------\nStatistics (after clipping):\n");
else
- printf("%g-sigma clipped %s:\n", p->sclipparams[0], mode);
+ printf("%g-%s clipped %s:\n", clipparams[0], sig1_mad0?"sigma":"MAD",
+ mode);
/* Print the final results: */
printf(" %-*s %zu\n", namewidth, "Number of input elements:",
p->input->size);
- if( p->sclipparams[1] < 1.0f )
- printf(" %-*s %d\n", namewidth, "Number of clips:", sigclip->status);
- printf(" %-*s %.0f\n", namewidth, "Final number of elements:", a[0]);
- printf(" %-*s %g\n", namewidth, "Median:", a[1]);
- printf(" %-*s %g\n", namewidth, "Mean:", a[2]);
- printf(" %-*s %g\n", namewidth, "Standard deviation:", a[3]);
+ if( clipparams[1] < 1.0f )
+ printf(" %-*s %.0f\n", namewidth, "Number of clips:",
+ a[GAL_STATISTICS_CLIP_OUTCOL_NUMBER_CLIPS]);
+ printf(" %-*s %.0f\n", namewidth, "Final number of elements:",
+ a[GAL_STATISTICS_CLIP_OUTCOL_NUMBER_USED]);
+ printf(" %-*s %.6e\n", namewidth, "Median:",
+ a[GAL_STATISTICS_CLIP_OUTCOL_MEDIAN]);
+ printf(" %-*s %.6e\n", namewidth, "Mean:",
+ a[GAL_STATISTICS_CLIP_OUTCOL_MEAN]);
+ printf(" %-*s %.6e\n", namewidth, "Standard deviation:",
+ a[GAL_STATISTICS_CLIP_OUTCOL_STD]);
+ printf(" %-*s %.6e\n", namewidth, "Median Absolute Deviation:",
+ a[GAL_STATISTICS_CLIP_OUTCOL_MAD]);
/* Clean up. */
+ gal_data_free(result);
free(mode);
}
@@ -1829,7 +1938,14 @@ statistics(struct statisticsparams *p)
if( p->sigmaclip )
{
print_basic_info=0;
- print_sigma_clip(p);
+ print_clip(p, 1);
+ }
+
+ /* Print the MAD-clipped results. */
+ if( p->madclip )
+ {
+ print_basic_info=0;
+ print_clip(p, 0);
}
/* Make the mirror table. */
diff --git a/bin/statistics/ui.c b/bin/statistics/ui.c
index efd6d815..343a80b3 100644
--- a/bin/statistics/ui.c
+++ b/bin/statistics/ui.c
@@ -409,13 +409,14 @@ ui_read_check_only_options(struct statisticsparams *p)
{
/* The tile or sky mode cannot be called with any other modes. */
if( p->asciihist || p->asciicfp || p->histogram || p->histogram2d
- || p->cumulative || p->sigmaclip || p->fitname
+ || p->cumulative || p->sigmaclip || p->madclip || p->fitname
|| !isnan(p->mirror) )
error(EXIT_FAILURE, 0, "'--ontile' or '--sky' cannot be called "
"with any of the 'particular' calculation options, for "
- "example '--histogram'. This is because the latter work "
- "over the whole dataset and element positions are changed, "
- "but in the former positions are significant");
+ "example '--histogram' or '--madclip'. Because these "
+ "operators change the order of the data and element "
+ "positions are changed, but in the former positions are "
+ "significant");
/* Make sure the tessellation defining options are given. */
if( tl->tilesize==NULL || tl->numchannels==NULL
@@ -469,12 +470,19 @@ ui_read_check_only_options(struct statisticsparams *p)
}
- /* Sigma-clipping needs 'sclipparams'. */
+ /* Sigma-clipping needs 'sclipparams' and MAD-clipping needs
+ 'mclipparams' */
if(p->sigmaclip && isnan(p->sclipparams[0]))
error(EXIT_FAILURE, 0, "'--sclipparams' is necessary with "
"'--sigmaclip'. '--sclipparams' takes two values (separated "
"by a comma) for defining the sigma-clip: the multiple of "
"sigma, and tolerance (<1) or number of clips (>1).");
+ if(p->madclip && isnan(p->mclipparams[0]))
+ error(EXIT_FAILURE, 0, "'--mclipparams' is necessary with "
+ "'--madclip' (median absolute deviation clipping). "
+ "'--mclipparams' takes two values (separated "
+ "by a comma) for defining the MAD-clip: the multiple of "
+ "MAD, and tolerance (<1) or number of clips (>1).");
/* If any of the mode measurements are requested, then 'mirrordist' is
@@ -491,6 +499,7 @@ ui_read_check_only_options(struct statisticsparams *p)
"mode-related single measurements ('--mode', "
"'--modequant', '--modesym', and '--modesymvalue')");
break;
+ case UI_KEY_SIGCLIPMAD:
case UI_KEY_SIGCLIPSTD:
case UI_KEY_SIGCLIPMEAN:
case UI_KEY_SIGCLIPNUMBER:
@@ -800,17 +809,23 @@ ui_make_sorted_if_necessary(struct statisticsparams *p)
case UI_KEY_MEDIAN:
case UI_KEY_QUANTILE:
case UI_KEY_QUANTFUNC:
+ case UI_KEY_MADCLIPMAD:
+ case UI_KEY_SIGCLIPMAD:
+ case UI_KEY_MADCLIPSTD:
case UI_KEY_SIGCLIPSTD:
case UI_KEY_QUANTOFMEAN:
case UI_KEY_SIGCLIPMEAN:
+ case UI_KEY_MADCLIPMEAN:
+ case UI_KEY_MADCLIPNUMBER:
case UI_KEY_SIGCLIPNUMBER:
+ case UI_KEY_MADCLIPMEDIAN:
case UI_KEY_SIGCLIPMEDIAN:
is_necessary=1;
break;
}
/* Check in the rest of the outputs. */
- if( is_necessary==0 && ( p->sigmaclip || !isnan(p->mirror) ) )
+ if( p->sigmaclip || p->madclip || !isnan(p->mirror) )
is_necessary=1;
/* Do the sorting, we will keep the sorted array in a separate space,
diff --git a/bin/statistics/ui.h b/bin/statistics/ui.h
index 11c9f610..f4a38494 100644
--- a/bin/statistics/ui.h
+++ b/bin/statistics/ui.h
@@ -81,6 +81,8 @@ enum option_keys_enum
UI_KEY_MINIMUM,
UI_KEY_MAXIMUM,
UI_KEY_SUM,
+ UI_KEY_MAD,
+ UI_KEY_MADCLIP,
UI_KEY_MODEQUANT,
UI_KEY_MODESYM,
UI_KEY_MODESYMVALUE,
@@ -106,10 +108,17 @@ enum option_keys_enum
UI_KEY_CHECKSKY,
UI_KEY_IGNOREBLANKINTILES,
UI_KEY_SCLIPPARAMS,
+ UI_KEY_MCLIPPARAMS,
UI_KEY_SIGCLIPNUMBER,
UI_KEY_SIGCLIPMEDIAN,
UI_KEY_SIGCLIPMEAN,
UI_KEY_SIGCLIPSTD,
+ UI_KEY_SIGCLIPMAD,
+ UI_KEY_MADCLIPNUMBER,
+ UI_KEY_MADCLIPMEDIAN,
+ UI_KEY_MADCLIPMEAN,
+ UI_KEY_MADCLIPSTD,
+ UI_KEY_MADCLIPMAD,
UI_KEY_GREATEREQUAL2,
UI_KEY_LESSTHAN2,
UI_KEY_ONEBINSTART2,
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index f32c93a9..d4d7aa08 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -20,6 +20,7 @@
@c So dashes and underscores can be used in HTMLs
@allowcodebreaks true
+
@c So function output type is printed on first line
@deftypefnnewline on
@@ -275,6 +276,7 @@ Tutorials
* Zero point of an image:: Estimate the zero point of an image.
* Pointing pattern design:: Optimizing the pointings of your observations.
* Moire pattern in stacking and its correction:: How to avoid this grid-based
artifact.
+* Clipping outliers:: How to avoid outliers in your measurements.
General program usage tutorial
@@ -345,6 +347,13 @@ Pointing pattern design
* Pointings that account for sky curvature:: Sky curvature will cause
problems!
* Accounting for non-exposed pixels:: Parts of the detector do not get
exposed to light.
+Clipping outliers
+
+* Building inputs and analysis without clipping:: Building a dataset for
demonstration below.
+* Sigma clipping:: Standard deviation (STD) clipping.
+* MAD clipping:: Median Absolute Deviation (MAD) clipping.
+* Filled re-clipping:: Two clips with holes filled in the middle.
+
Installation
* Dependencies:: Necessary packages for Gnuastro.
@@ -619,7 +628,6 @@ Statistics
* Histogram and Cumulative Frequency Plot:: Basic definitions.
* 2D Histograms:: Plotting the distribution of two variables.
-* Sigma clipping:: Definition of @mymath{\sigma}-clipping.
* Least squares fitting:: Fitting with various parametric functions.
* Sky value:: Definition and derivation of the Sky value.
* Invoking aststatistics:: Arguments and options to Statistics.
@@ -2044,6 +2052,7 @@ For an explanation of the conventions we use in the
example codes through the bo
* Zero point of an image:: Estimate the zero point of an image.
* Pointing pattern design:: Optimizing the pointings of your observations.
* Moire pattern in stacking and its correction:: How to avoid this grid-based
artifact.
+* Clipping outliers:: How to avoid outliers in your measurements.
@end menu
@@ -9986,7 +9995,7 @@ So we will leave this as an exercise.
-@node Moire pattern in stacking and its correction, , Pointing pattern
design, Tutorials
+@node Moire pattern in stacking and its correction, Clipping outliers,
Pointing pattern design, Tutorials
@section Moir@'e pattern in stacking and its correction
@cindex Moir@'e pattern or fringes
@@ -10212,11 +10221,683 @@ Note that this discussion is on small shifts between
pointings (dithers), not la
+@node Clipping outliers, , Moire pattern in stacking and its correction,
Tutorials
+@section Clipping outliers
+
+@cindex Outlier
+@cindex Clipping of outliers
+Outliers occur often in data sets.
+For example cosmic rays in astronomical imaging: the image of your target
galaxy can be affected by a cosmic ray in one of the five exposures you took in
one night.
+As a result, when you compare the measured magnitude of your target galaxy in
all the exposures, you will get measurements like this (all in magnitudes)
19.8, 20.1, 20.5, 17.0, 19.9 (all fluctuating around magnitude 20, except the
much brighter 17th magnitude measurement).
+
+Normally, you would simply take the mean of these measurements to estimate the
magnitude of your target with more precision.
+However, the 17th magnitude measurement above is clearly wrong and will
significantly affect the mean: without it, the mean magnitude is 20.07, but
with it, the mean is 19.46:
+
+@example
+$ echo " 19.8 20.1 20.5 17 19.9" \
+ | tr ' ' '\n' \
+ | aststatistics --mean
+1.94600000000000e+01
+
+$ echo " 19.8 20.1 20.5 19.9" \
+ | tr ' ' '\n' \
+ | aststatistics --mean
+2.00750000000000e+01
+@end example
+
+This difference of 0.61 magnitudes (or roughly 1.75 times) is significant (for
the definition of magnitudes in astronomy, see @ref{Brightness flux magnitude}).
+In the simple example above, you can visually identify the ``outlier'' and
manually remove it.
+But in most common situations you will not be so lucky!
+For example when you want to stack the five images of the five exposures
above, and each image has @mymath{4000\times4000} (or 16 million!) pixels and
not possible by hand in a reasonable time (an average human's lifetime!).
+
+This tutorial reviews the effect of outliers and different available ways to
remove them.
+In particular, we will be looking at stacking of multiple datasets and
collapsing one dataset along one of its dimensions.
+But the concepts and methods are applicable to any analysis that is affected
by outliers.
+
+@menu
+* Building inputs and analysis without clipping:: Building a dataset for
demonstration below.
+* Sigma clipping:: Standard deviation (STD) clipping.
+* MAD clipping:: Median Absolute Deviation (MAD) clipping.
+* Filled re-clipping:: Two clips with holes filled in the middle.
+@end menu
+
+@node Building inputs and analysis without clipping, Sigma clipping, Clipping
outliers, Clipping outliers
+@subsection Building inputs and analysis without clipping
+
+As described in @ref{Clipping outliers}, the goal of this tutorial is to
demonstrate the effects of outliers and show how to ``clip'' them from basic
statistics measurements.
+This is best done on an actual dataset (rather than pure theory).
+In this section we will build nine noisy images with the script below, such
that one of the images has a circle in the middle.
+We will then stack the 9 images into one final image based on different
statistical measurements: the mean, median, standard deviation (STD), median
absolute deviation (MAD) and number of inputs used in each pixel.
+We will then analyze the resulting stacks to demonstrate the problem with
outliers.
+
+Put the script below into a plain-text file (assuming it is called
@file{script.sh}), and run it with @command{bash ./script.sh}.
+For more on writing and good practices in shell scripts, see @ref{Writing
scripts to automate the steps}.
+The last command of the script above calls DS9 to visualize the five output
stacked images mentioned above.
+
+@verbatim
+# Constants
+list=""
+sigma=10
+number=9
+radius=30
+width=201
+bdir=build
+profsum=3e5
+background=10
+random_seed=1699270427
+
+
+# Clipping parameters (will be set later when we start clipping).
+# clip_multiple: 3 for sigma; 4.48 for MAD
+# clip_tolerance: 0.1 for sigma; 0.01 for MAD
+clip_operator=""
+clip_multiple=""
+clip_tolerance=""
+
+
+# Stop if there is any error.
+set -e
+
+
+# If the build directory does not exist, build it.
+if ! [ -d $bdir ]; then mkdir $bdir; fi
+
+
+# The final image (with largest number) will contain the outlier:
+# we'll put a flat circle in the center of the image as the outlier
+# structure.
+outlier=$bdir/in-$number.fits
+nn=$bdir/$number-no-noise.fits
+export GSL_RNG_SEED=$random_seed
+center=$(echo $width | awk '{print int($1/2)+1}')
+echo "1 $center $center 5 $radius 0 0 1 $profsum 1" \
+ | astmkprof --mode=img --mergedsize=$width,$width \
+ --oversample=1 --output=$nn --mcolissum
+astarithmetic $nn $background + $sigma mknoise-sigma \
+ --envseed -o$outlier
+
+
+# Build pure noise and add elements to the list of images to stack.
+list=$outlier
+numnoise=$(echo $number | awk '{print $1-1}')
+for i in $(seq 1 $numnoise); do
+ img="$bdir/in-$i.fits"
+ if ! [ -f $img ]; then
+ export GSL_RNG_SEED=$(echo $random_seed | awk '{print $1+'$i'}')
+ astarithmetic $width $width 2 makenew float32 $background + \
+ $sigma mknoise-sigma --envseed --output=$img
+ fi
+ list="$list $img"
+done
+
+
+# Stack the images,
+for op in mean median std mad number; do
+ if [ x"$clip_operator" = x ]; then
+ out=$bdir/stack-$op.fits
+ astarithmetic $list $number $op -g1 --output=$out
+ else
+ operator=$clip_operator-$op
+ out=$bdir/stack-$operator.fits
+ astarithmetic $list $number $clip_multiple $clip_tolerance \
+ $operator -g1 --output=$out
+ fi
+done
+
+
+# Collapse the first and last image along the 2nd dimension.
+for i in 1 $number; do
+ if [ x"$clip_operator" = x ]; then
+ out=$bdir/collapsed-$i.fits
+ astarithmetic $bdir/in-$i.fits 2 collapse-median counter \
+ --writeall --output=$out
+ else
+ out=$bdir/collapsed-$clip_operator-$i.fits
+ astarithmetic $bdir/in-$i.fits $clip_multiple $clip_tolerance \
+ 2 collapse-$clip_operator-median counter \
+ --writeall --output=$out
+ fi
+done
+@end verbatim
+
+@noindent
+After the script finishes, you can see the generated input images with the
first command below.
+The second command shows the stacked images.
+
+@example
+$ astscript-fits-view build/in-*.fits --ds9extra="-lock scalelimits yes"
+$ astscript-fits-view build/stack-*.fits
+@end example
+
+Color-blind readers may not clearly see the issue in the opened images with
this color bar.
+In this case, please choose the ``color'' menu at the top of the DS9 and
select ``gray'' or any other color that makes the circle most visible.
+
+The effect of an outlier on the different measurements above can be visually
seen (and quatitatively measured) through the visibility of the circle (that
was only present in one image, of nine).
+Let's look at them one by one (from the one that is most affected to the
least):
+
+@table @file
+@item std.fits
+The standard deviation (third image in DS9) is the most strongly affected
statistic by an outlier.
+This is so strong that the edge of the circle is also clearly visible!
+The standard deviation is calculated by first finding th mean, and estimating
the difference of each element from the mean.
+Those differences are then taken to the power of two and finally the square
root is taken (after a division by the number).
+It is the power-of-two component that amplifies the effect of the single
outlier as you see here.
+
+@item mean.fits
+The mean (first image in DS9) is also affected by the outlier such
+
+@item median.fits
+The median (second image in DS9) is also affected by the outlier; although
much less significantly than the standard deviation or mean.
+At first sight the circle may not be too visible!
+To see it more clearly, click on the ``Analysis'' menu in DS9 and then the
``smooth'' item.
+After smoothing, you will see how the single outlier has leaked into the
median stack.
+
+Intuitively, we would think that since the median is calculated from the
middle element after sorting, the outlier goes to the end and won't affect the
result.
+However, this is not the case as we see here: with 9 elements, the ``central''
element is the 5th (counting from 1; after sorting).
+Since the pixels covered by the circle only have 8 pure noise elements; the
``true'' median should have been the average of the 4th and 5th elements (after
sorting).
+By definition, the 5th element is always larger than the mean of the 4th and
5th (because the 4th element after sorting has a smaller value than the 5th
element).
+Therefore, using the 5th element (after sorting), we are systematically
choosing higher noise values in regions that are covered by the circle!
+
+With larger datasets, the difference between the central elements will be less.
+However, the improved precision (in the elements without an outlier) will also
be more.
+A detailed analysis of the effect of a single outlier on the median based on
the number of inputs can be done as an excersize; but in general, as this
argument shows, the median is not immune to outliers; especially when you care
about low signal-to-noise signal (as we do in astronomy: low surface brightness
science).
+
+@item mad.fits
+The median absolute deviation (fourth image in DS9) is affected by outliers in
a similar fashion to the median.
+
+@item number.fits
+The number image (last image in DS9) shows the number of images that went into
each pixel.
+Since we haven't rejected any outliers (yet!), all the pixels in this image
have a value of 9.
+@end table
+
+The example above included a single outlier.
+But we are not usually that lucky: there are usually more outliers!
+For example, the last command in the script above collapsed @file{1.fits}
(that was pure noise, without the circle) and @file{9.fits} (with the circle)
along their second dimension (the vertical).
+Collapsing was done by taking the median along all the pixels in the vertical
dimension.
+The output of collapsing has one less dimension; in this case, producing a 1D
table (with the same number of rows as the image's horizontal axis).
+To easily plot the output afterwards, we have also used the @code{counter}
operator.
+With the command below, you can open both tables and compare them:
+
+@example
+$ astscript-fits-view build/collapsed-*.fits
+@end example
+
+The last command opens TOPCAT.
+In the ``Graphics'' menu, select plane plot and you will see all the values
fluctuating around zero (with a maximum/minimum around @mymath{\pm2}).
+Afterwards, click on the ``Layers'' menu and click on ``Add position control''.
+In the opened tab at the bottom (where the scroll bar infront of ``Table'' is
empty), select the other table.
+In the regions that there was no circle in any of the vertical axises, the two
match nicely (the noise level is the same).
+However, you see that the image columns that were partly covered by the
outlying circle gradually get more affected as the width of the circle in that
column increases (the full diameter of the circle was in the middle of the
image).
+This shows how the median is biased by outliers as their number increases.
+
+To see the problem more prominently, use the @code{collapse-mean} operator
instead of the median.
+The reason that the mean is more strongly affected by the outlier is exactly
the same as above for the stacking of the input images.
+In the subsections below, we will describe some of the basic ways to reject
the effect of these outliers (and have better stacks or collapses).
+But the methodology is not limited to these types of data and can be
generically applied; unless specified explicitly.
+
+@node Sigma clipping, MAD clipping, Building inputs and analysis without
clipping, Clipping outliers
+@subsection Sigma clipping
+
+Let's assume that you have pure noise (centered on zero) with a clear
@url{https://en.wikipedia.org/wiki/Normal_distribution,Gaussian distribution},
or see @ref{Photon counting noise}.
+Now let's assume you add very bright objects (signal) on the image which have
a very sharp boundary.
+By a sharp boundary, we mean that there is a clear cutoff (from the noise) at
the pixels the objects finish.
+In other words, at their boundaries, the objects do not fade away into the
noise.
+
+In optical astronomical imaging, cosmic rays (when they collide at a near
normal incidence angle) are a very good example of such outliers.
+The tracks they leave behind in the image are perfectly immune to the blurring
caused by the atmosphere on images of stars or galaxies and they have a very
high signal-to-noise ratio.
+They are also very energetic and so their borders are usually clearly
separated from the surrounding noise.
+See Figure 15 in Akhlaghi and Ichikawa,
@url{https://arxiv.org/abs/1505.01664,2015}.
+
+In such a case, when you plot the histogram (see @ref{Histogram and Cumulative
Frequency Plot}) of the distribution, the pixels relating to those objects will
be clearly separate from pixels that belong to parts of the image that did not
have any signal (were just noise).
+In the cumulative frequency plot, after a steady rise (due to the noise), you
would observe a long flat region were for a certain range of data (horizontal
axis), there is no increase in the index (vertical axis).
+
+In the previous section (@ref{Building inputs and analysis without clipping})
we created one such dataset (@file{9.fits}).
+With the command below, let's have a look at its histogram and cumulative
frequency plot (in simple ASCII format; we are decreasing the default number of
bins with @option{--numasciibins} to show them easily within the width of the
print version of this manual; feel free to change this).
+
+@example
+$ aststatistics build/in-9.fits --asciihist --asciicfp \
+ --numasciibins=65
+
+ASCII Histogram:
+Number: 40401
+Y: (linear: 0 to 4191)
+X: (linear: -31.9714 -- 150.323, in 65 bins)
+ | **
+ | ****
+ | ******
+ | ******
+ | ********
+ | ********
+ | **********
+ | ************
+ | **************
+ | ****************** ******
+ |******************************* *************************
+ |-----------------------------------------------------------------
+
+
+ASCII Cumulative frequency plot:
+Y: (linear: 0 to 40401)
+X: (linear: -31.9714 -- 150.323, in 65 bins)
+ | ***************
+ | **********************************************
+ | ***********************************************
+ | *************************************************
+ | **************************************************
+ | ***************************************************
+ | ****************************************************
+ | *****************************************************
+ | ******************************************************
+ | ********************************************************
+ |*****************************************************************
+ |-----------------------------------------------------------------
+@end example
+
+@cindex Cosmic rays
+Outliers like the example above can significantly bias the measurement of the
background noise statistics.
+For example let's compare the median, mean and standard deviation of the image
above with @file{1.fits}:
+
+@example
+$ aststatistics build/in-1.fits --median --mean --std
+9.90529778313248e+00 9.96143102101206e+00 1.00137568561776e+01
+
+$ aststatistics build/in-9.fits --median --mean --std
+1.09305819367634e+01 1.74470443173776e+01 2.88895986970341e+01
+@end example
+
+The effect of the outliers is obvious in all three measures: the median has
become 1.10 times larger, the mean 1.75 times and the standard deviation about
2.88 times!
+The differing effect of outliers in different statistics was already discussed
in @ref{Building inputs and analysis without clipping}; also see
@ref{Quantifying signal in a tile}.
+
+@mymath{\sigma}-clipping is one way to remove/clip the effect of such very
strong outliers in measures like the above.
+@mymath{\sigma}-clipping is defined as the very simple iteration below.
+In each iteration, the range of input data might decrease.
+When the outliers are as strong as above, the outliers will be removed through
this iteration.
+
+@enumerate
+@item
+Calculate the standard deviation (@mymath{\sigma}) and median (@mymath{m}) of
a distribution.
+The median used because, as shown above, the mean is too significantly
affected by the presence of outliers.
+@item
+Remove all points that are smaller or larger than @mymath{m\pm\alpha\sigma}.
+@item
+Go back to step 1, unless the selected exit criteria is reached.
+There are commonly two types of exit criteria (to stop the
@mymath{\sigma}-clipping iteration).
+Within Gnuastro's programs that use sigma-clipping, the exit criteria is the
second value to the @option{--sclipparams} option (the first value is the
@mymath{m} above):
+
+@itemize
+@item
+When a certain number of iterations has taken place (exit criteria is an
integer, larger than 1).
+@item
+When the new measured standard deviation is within a certain tolerance level
of the previous iteration (exit criteria is floating point and less than 1.0).
+The tolerance level is defined by:
+
+@dispmath{\sigma_{old}-\sigma_{new} \over \sigma_{new}}
+
+In each clipping, the dispersion in the distribution is either less or equal.
+So @mymath{\sigma_{old}\geq\sigma_{new}}.
+@end itemize
+@end enumerate
+
+Let's see the algorithm in practice with the @option{--sigmaclip} option of
Gnuastro's Statistics program (using the default configuration of
@mymath{3\sigma} clipping and tolerance of 0.1):
+
+@example
+$ aststatistics build/in-9.fits --sigmaclip
+Statistics (GNU Astronomy Utilities) @value{VERSION}
+-------
+Input: build/in-9.fits (hdu: 1)
+-------
+3-sigma clipping steps until relative change in STD is less than 0.1:
+
+round number median STD
+1 40401 1.09306e+01 2.88896e+01
+2 37660 1.00306e+01 1.07153e+01
+3 37539 1.00080e+01 9.93741e+00
+-------
+Statistics (after clipping):
+ Number of input elements: 40401
+ Number of clips: 2
+ Final number of elements: 37539
+ Median: 1.000803e+01
+ Mean: 1.001822e+01
+ Standard deviation: 9.937410e+00
+ Median Absolute Deviation: 6.772760e+00
+@end example
+
+After the basic information about the input and settings, the Statistics
program has printed the information for each round (iteration) of clipping.
+Initially, there was 40401 elements (the image is @mymath{201\times201}
pixels).
+After the first round of clipping, only 37660 elements remained and because
the difference in standard deviation was larger than the tolerance level, a
third clipping was one.
+But the change in standard deviation after the third clip was smaller than the
tolerance level, so the exit criteria was activated and the clipping finished
with 37539 elements.
+In the end, we see that the final median, mean and standard deviation are very
similar to the data without any outlier (@file{build/1.fits} in the example
above).
+
+The example above provided a single statistic from a single dataset.
+Other scenarios where sigma-clipping becomes necessary are stacking and
collapsing (that was the main goal of the script in @ref{Building inputs and
analysis without clipping}).
+To generate @mymath{\sigma}-clipped stacks and collapsed tables, you just need
to change the values of the three variables of the script (shown below).
+After making this change in your favorite text editor, have a look at the
outputs:
+
+@example
+$ grep ^clip_ script.sh
+clip_operator=sigclip # These variables will be used more
+clip_multiple=3 # effectively with the clipping
+clip_tolerance=0.1 # operators of the next sections.
+
+$ bash ./script.sh
+
+$ astscript-fits-view build/stack-std.fits \
+ build/stack-sigclip-std.fits \
+ build/stack-*mean.fits \
+ build/stack-*median.fits \
+ build/stack-*number.fits \
+ --ds9extra="-tile grid layout 2 4 -scale minmax"
+@end example
+
+It is clear that the @mymath{\sigma}-clipped results have significantly
improved in all four measures (images in the right column in DS9).
+The reason is clear in the @file{stack-sigclip-number.fits} image (which show
how many images were used for each pixel): almost all of the outlying circle
has been accounted for (the pixel values are 8, not 9, showing 8 images went
into those).
+It is the leaked holes in the @file{stack-sigclip-number.fits} image (with
value of 9) that keep the circle in the final stack of the other measures (at
various levels).
+See @ref{Filled re-clipping} for stacking operators that can account for this.
+
+So far, @mymath{\sigma}-clipping has preformed nicely.
+However, there are important caveats to @mymath{\sigma}-clipping that are
listed in the box below and further elaborated (with examples) afterwards.
+
+@cartouche
+@noindent
+@strong{Caveats of @mymath{\sigma}-clipping}: There are some important caveats
to @mymath{\sigma}-clipping:
+@itemize
+@item
+The standard deviation is itself heavily influenced by the presence of
outliers.
+Therefore a sufficiently small number of outliers can expand the standard
deviation such that they stay within the boundaries.
+
+@item
+When the outliers do not consititute a clearly distinct distribution like the
example here, sigma-clipping will not be able to separate them like here.
+@end itemize
+@end cartouche
+
+To demonstrate how weaker outliers will not be clipped in sigma clipping,
let's decrease the total sum of values in the outlying circle, then re-run the
script:
+
+@example
+$ grep ^profsum script.sh
+profsum=1e5
+
+$ bash ./script.sh
+@end example
+
+Let's have a look at the new outlying circle with the first command below.
+With the second command, let's view its pixel value histogram (recall that
previously, the circle had a clearly separate distribution):
+
+@example
+$ astscript-fits-view build/in-9.fits
+
+$ aststatistics build/in-9.fits --asciihist --numasciibins=65
+ASCII Histogram:
+Number: 40401
+Y: (linear: 0 to 2654)
+X: (linear: -31.9714 -- 79.4266, in 65 bins)
+ | **
+ | *****
+ | *********
+ | **********
+ | *************
+ | **************
+ | *****************
+ | *******************
+ | ***********************
+ | ****************************************
+ |*****************************************************************
+ |-----------------------------------------------------------------
+@end example
+
+We see that even tough the circle is still clearly visible in the noise, the
histogram is not longer separate; it has blended into the noise, and just
caused a skewnewss in the otherwise symmetric noise distribution.
+Let's try running the @option{--sigmaclip} option as above:
+
+@example
+$ aststatistics build/in-9.fits --sigmaclip
+Statistics (GNU Astronomy Utilities) @value{VERSION}
+-------
+Input: build/in-9.fits (hdu: 1)
+-------
+3-sigma clipping steps until relative change in STD is less than 0.1:
+
+round number median STD
+1 40401 1.09295e+01 1.34784e+01
+2 39618 1.06762e+01 1.19852e+01
+3 39126 1.05265e+01 1.12983e+01
+-------
+Statistics (after clipping):
+ Number of input elements: 40401
+ Number of clips: 2
+ Final number of elements: 39126
+ Median: 1.052652e+01
+ Mean: 1.114819e+01
+ Standard deviation: 1.129831e+01
+ Median Absolute Deviation: 7.106166e+00
+@end example
+
+We see that the median, mean and standard deviation are over estimated (each
worse than the previous!).
+Let's see how the @mymath{\sigma}-clipping stacking on this outlying flat
circle:
+
+@example
+$ astscript-fits-view build/stack-std.fits \
+ build/stack-sigclip-std.fits \
+ build/stack-*mean.fits \
+ build/stack-*median.fits \
+ build/stack-*number.fits \
+ --ds9extra="-tile grid layout 2 4 -scale minmax"
+@end example
+
+Compared to the previous run (where the outlying circle was brighter), we see
that @mymath{\sigma}-clipping is now less successful in removing the outlying
circle from the stacks; or in the single value measurements.
+This is particularly visible in the @file{stack-sigclip-number.fits} image:
the circle barely visible any more: there is only a very weak clustering of
pixels with a value of 8 over the circle's pixels.
+This has happened because the outliers have biased the standard deviation
itself to a level that includes them with this multiple of the standard
deviation.
+
+To gauge if @mymath{\sigma}-clipping will be useful for your dataset, look at
the histogram (see @ref{Histogram and Cumulative Frequency Plot}).
+The ASCII histogram that is printed on the command-line with
@option{--asciihist} (like above) is good enough in most cases.
+But you can't do this manually in every case (as in the stacking which
involved more than forty thousand pixels)!
+Clipping outliers should be based on a measure of scatter that is less
affected by outliers!
+Therefore, in Gnuastro we also have median absolute deviation (MAD) clipping
which is described in the next section (@ref{MAD clipping}).
+@node MAD clipping, Filled re-clipping, Sigma clipping, Clipping outliers
+@subsection MAD clipping
+@cindex Median absolute deviation (MAD)
+@cindex MAD (median absolute deviation)
+When clipping outliers, it is important that the used measure of dispersion is
itself not strongly affected by the outliers.
+Previously (in @ref{Sigma clipping}), we saw that the standard deviation is
not a good measure of dispersion because of its strong dependency on outliers.
+In this section, we'll introduce clipping operators that are based on the
@url{https://en.wikipedia.org/wiki/Median_absolute_deviation, median absolute
deviation} (MAD).
+The median absolute deviation is defined as the median of the differences of
each element from the median of the elements.
+As mathematically derived in the Wikipedia page above, for a pure Gaussian
distribution, the median absolute deviation will be roughly
@mymath{0.67449\sigma}.
+We can confirm this numerically from the images with pure noise that we
created previously in @ref{Building inputs and analysis without clipping}.
+With the first command below we can see the raw standard deviation and median
absolute deviation values and the second command shows their division:
+@verbatim
+$ aststatistics build/in-1.fits --std --mad
+1.00137568561776e+01 6.74662296703343e+00
+$ aststatistics build/in-1.fits --std --mad | awk '{print $2/$1}'
+0.673735
+@end verbatim
+
+The algorithm of MAD-clipping is identical to @mymath{\sigma}-clipping, except
that instead of @mymath{\sigma}, it uses the median absolute deviation.
+Since the median absolute deviation is smaller than the standard deviation by
roughly 0.67, if you regularly use @mymath{3\sigma} there, you should use
@mymath{(3/0.67)\rm{MAD}=(4.48)\rm{MAD}} when doing MAD-clipping.
+The usual tolerance should also be changed due to the differing nature of the
median absolute deviation (based on sorted differences) in relation to the
standard deviation (based on the sum of squared differences).
+A tolerance of 0.01 is better suited to the termination criteria of
MAD-clipping.
+
+To demonstrate the steps in practice, let's assume you have the original
script in @ref{Building inputs and analysis without clipping} with the changes
shown in the first command below.
+With the second command we'll execute the script, and with the third command
we'll do the iterations of MAD-clipping:
+
+@verbatim
+$ grep '^clip_\|^profsum' script.sh
+profsum=1e5
+clip_operator=madclip
+clip_multiple=4.48
+clip_tolerance=0.01
+
+$ bash ./script.sh
+
+$ aststatistics build/in-9.fits --madclip
+Statistics (GNU Astronomy Utilities) 0.21.6-28a1
+-------
+Input: build/in-9.fits (hdu: 1)
+-------
+4.48-MAD clipping steps until relative change in MAD
+(median absolute deviation) is less than 0.01:
+
+round number median MAD
+1 40401 1.09295e+01 7.38609e+00
+2 38789 1.04261e+01 7.03508e+00
+3 38549 1.03469e+01 6.97927e+00
+-------
+Statistics (after clipping):
+ Number of input elements: 40401
+ Number of clips: 2
+ Final number of elements: 38549
+ Median: 1.034690e+01
+ Mean: 1.068946e+01
+ Standard deviation: 1.062083e+01
+ Median Absolute Deviation: 6.979274e+00
+@end verbatim
+
+We see that the median, mean and standard deviation after MAD-clipping is much
better than the basic @mymath{\sigma}-clipping (see @ref{Sigma clipping}): the
median is now 10.3 (was 10.5 in @mymath{\sigma}-clipping), mean is 10.7 (was
10.11) and the standard deviation is 10.6 (was 10.12).
+
+Let'scompare the MAD-clipped stacks with the results of the previous section.
+Since we want the images shown in a certain order, we'll first construct the
list of images (with a @code{for} loop that will fill the @file{imgs} variable).
+Note that this assumes you have ran and carefully read/understand all the
commands in the previous sections (@ref{Building inputs and analysis without
clipping} and @ref{Sigma clipping}).
+Tip: the three @option{--ds9extra} options ensure that the bottom row (showing
the number of images used in each pixel) has the same scale and limits in all
three columns.
+
+@example
+$ imgs=""
+$ p=build/stack # 'p' is short for "prefix"
+$ for m in std mean median mad number; do \
+ imgs="$imgs $p-$m.fits $p-sigclip-$m.fits $p-madclip-$m.fits"; \
+ done
+$ astscript-fits-view $imgs --ds9extra="-tile grid layout 3 5" \
+ --ds9extra="-scale limits 1 9" \
+ --ds9extra="-frame move back -scale limits 1 9" \
+ --ds9extra="-frame move back -scale limits 1 9"
+@end example
+
+The third column shows the newly created MAD-clipped stacks.
+We see that the outlying circle is much more weaker in the MAD-clipped stacks
than in the @mymath{\sigma}-clipped stacks in all measures (except for the
``number'' measure where the circle should be stronger).
+
+However, unfortunately even MAD-clipping is not perfect and we still see the
circle in all four cases, even with the MAD-clipped median (more clearly: after
smoothing/blocking).
+The reason is similar to what was described in @mymath{\sigma}-clipping (using
the original @code{profsum=3e5}: the leaked holes in the numbers image.
+Because the circle is not too distant from the noise some of its elements do
not get clipped, and their stacked value gets systematically higher than the
rest of the image.
+In Gnuastro, we have a fix for this that is described fully in the next
section (@ref{Filled re-clipping}).
+
+@node Filled re-clipping, , MAD clipping, Clipping outliers
+@subsection Filled re-clipping
+
+When source of the outlier covers more than one element, and its flux is close
to the noise level, not all of its elements will be clipped: because of noise,
some of its elements will remain un-clipped; and thus affecting the output.
+Examples of this were created and thoroughly discussed in previous sections
with @mymath{\sigma}-clipping and MAD-clipping (see @ref{Sigma clipping} and
@ref{MAD clipping}).
+
+To avoid this problem, in Gnuastro we have an extra set of clipping operators
that will do two rounds of clips and with some basic operations in the middle.
+@enumerate
+@item
+The requested clipping is first applied.
+This will the return the median and dispersion measure (MAD or STD).
+@item
+A mask is created for each input image (in stacking) or 1D array (in
collapsing).
+Any pixel that is outside the requested clip range is set to 1; the rest are
set to 0.
+@item
+Isolated regions are found:
+@itemize
+@item
+For 2D images (were each pixel has 8 neighbors) the mask pixels are first
dilated (so the edges of the regions are closed off from the surrounding noise).
+@item
+For 1D arrays (where each element only has two neighbors), the mask is first
eroded.
+This is necessary because the next step (where the holes are filled), two
pixels that have been clipped purely due to noise with a large distance between
them can wrongly mask a very large range of the input data.
+@end itemize
+@item
+Any 0-valued pixel in the masks that are fully surrounded by 1s (or ``holes'')
are filled (given a value of 1).
+@item
+All the pixels that have a value of 1 in the mask are set to NaN in the
respective input data (that the mask corresponds to).
+@item
+The requested clipping is repeated on the newly masked inputs.
+@end enumerate
+
+Through this process, the less significant outliers (which do not get clipped
independently) are clipped based on their surrounding elements.
+The filled re-clipping operators have an extra @code{-fill} in their names.
+For example the filled MAD-clipped mean is called @code{madclip-fill-mean}
(while the simple MAD-clipped mean operator was called @code{madclip-mean}).
+Let's run our script with the filled @mymath{\sigma}-clipping and
@mymath{MAD}-clipping (before each run, make sure the values shown under the
@code{grep} command are correct).
+
+With the last command, we'll view all the outputs generated so far (in this
section and the previous ones (@ref{Building inputs and analysis without
clipping}, @ref{Sigma clipping} and @ref{MAD clipping}):
+
+@verbatim
+$ grep '^clip_\|^profsum' script.sh
+profsum=1e5
+clip_operator=madclip-fill
+clip_multiple=4.48
+clip_tolerance=0.01
+
+$ bash ./script
+
+$ $ grep '^clip_\|^profsum' script.sh
+profsum=1e5
+clip_operator=sigclip-fill
+clip_multiple=3
+clip_tolerance=0.1
+
+$ bash ./script
+
+$ imgs=""
+$ for m in std mean median mad number; do \
+ imgs="$imgs $p-$m.fits $p-sigclip-$m.fits $p-sigclip-fill-$m.fits" \
+ imgs="$p-madclip-$m.fits $p-madclip-fill-$m.fits"; \
+ done
+$ astscript-fits-view $imgs --ds9extra="-tile grid layout 5 6" \
+ --ds9extra="-scale limits 1 9" \
+ --ds9extra="-frame move back -scale limits 1 9" \
+ --ds9extra="-frame move back -scale limits 1 9" \
+ --ds9extra="-frame move back -scale limits 1 9" \
+ --ds9extra="-frame move back -scale limits 1 9"
+@end verbatim
+
+The last column (belonging to the @code{madclip-fill-*} operators) is
@emph{finally} free of the outlying circle (that was present in only one of
nine inputs).
+The filling operation did not affect the @code{sigclip-fill-*} operators
(third column DS9 from the last command above)!
+The reason is clear from the bottom row (showing the number of images used in
each pixel).
+The weak over density of clipped pixels over the circle is barely visible and
was not strong enough for defining ``holes'' (that will be filled).
+On the contrary, when comparing the @file{madclip-number.fits} and
@file{madclip-fill-number.fits}, the filled holes within the circle are clearly
visible.
+
+But the script also generated collapsed columns of @file{build/in-9.fits} (to
a 1D table).
+In this case, for each column, the number of outliers increase as we enter the
circle and reach a maximum in the middle of the image.
+Let's have a look at those outputs:
+
+@example
+$ astscript-fits-view build/collapsed-madclip-9.fits \
+ build/collapsed-madclip-fill-9.fits
+@end example
+
+When comparing the two in TOPCAT (following the same process described in
@ref{Building inputs and analysis without clipping}) you will notice that the
difference is only in the edges of the circle.
+The 4.48 multiple of MAD-clipping (corresponding to 3 sigma), was not
successful in removing the many outlying pixels due to the circle in the
central pixels of the image.
+
+This is a relatively high threshold and was used because for the images, we
only had 9 elements in each clipping for every pixel.
+But for the collapsing, we have many more pixels in each vertical direction fo
the image (201 pixels).
+Let's decrease the threshold to 3 and calculate the collapsed mean after
MAD-clipping, once with filled re-clipping and once without it:
+
+@example
+$ for m in mean number; do \
+ for clip in madclip madclip-fill; do \
+ astarithmetic build/in-9.fits 3 0.01 2 collapse-$clip-$m \
+ counter --writeall -ocollapse-$clip-$m.fits; \
+ done; \
+ done
+@end example
+
+The two loops above created four tables.
+First, with the command below, let's look at the two measured mean values (one
with filling and the other without it):
+
+@example
+$ astscript-fits-view collapse-*-mean.fits
+@end example
+
+In the table without filled re-clipping, you see a small shift in the center
of the image (around 100 in the horizontal axis).
+Let's have a look at the final number of pixels used in each clipping:
+
+@example
+$ astscript-fits-view collapse-*-number.fits
+@end example
+
+The difference is now clearly visible when you plot both in one ``Plane plot''
window.
+In the filled re-clipping case, we see a clear dip in the number of pixels
that very nicely corresponds to the number of pixels associated to the circle.
+But the dip is much more noisy in the simple MAD-clipping.
@@ -10310,8 +10991,8 @@ The @url{http://www.gnu.org/software/gsl/, GNU
Scientific Library}, or GSL, is a
To download and install GSL from source, you can run the following commands.
@example
-$ wget http://ftpmirror.gnu.org/gsl/gsl-latest.tar.gz
-$ tar xf gsl-latest.tar.gz
+$ wget https://ftp.gnu.org/gnu/gsl/gsl-latest.tar.gz
+$ tar -xf gsl-latest.tar.gz
$ cd gsl-X.X # Replace X.X with version number.
$ ./configure CFLAGS="$CFLAGS -g0 -O3"
$ make -j8 # Replace 8 with no. CPU threads.
@@ -10354,7 +11035,7 @@ The commands necessary to download the source,
decompress, build and install CFI
@example
$ urlbase=http://heasarc.gsfc.nasa.gov/FTP/software/fitsio/c
$ wget $urlbase/cfitsio_latest.tar.gz
-$ tar xf cfitsio_latest.tar.gz
+$ tar -xf cfitsio_latest.tar.gz
$ cd cfitsio-X.XX # Replace X.XX with version
$ ./configure --prefix=/usr/local --enable-sse2 --enable-reentrant \
CFLAGS="$CFLAGS -g0 -O3"
@@ -10412,7 +11093,7 @@ To download, configure, build, check and install WCSLIB
from source, you can fol
@example
## Download and unpack the source tarball
$ wget ftp://ftp.atnf.csiro.au/pub/software/wcslib/wcslib.tar.bz2
-$ tar xf wcslib.tar.bz2
+$ tar -xf wcslib.tar.bz2
## In the `cd' command, replace `X.X' with version number.
$ cd wcslib-X.X
@@ -10779,6 +11460,13 @@ This arguably makes Debian-based OSs the largest, and
most used, class of GNU/Li
All of them use Debian's Advanced Packaging Tool (APT, for example,
@command{apt-get}) for managing packages.
@table @asis
+@item Development features (Ubuntu or derivatives)
+By default, a newly installed Ubuntu does not contain the low-level tools that
are necessary for building a software from source.
+Therefore, if you are using Ubuntu, please run the following command.
+@example
+$ sudo apt-get install gcc make zlib1g-dev lzip
+@end example
+
@item Mandatory dependencies
Without these, Gnuastro cannot be built, they are necessary for input/output
and low-level mathematics (see @ref{Mandatory dependencies})!
@example
@@ -10791,7 +11479,7 @@ If present, these libraries can be used in Gnuastro's
build for extra features,
@example
$ sudo apt-get install ghostscript libtool-bin \
libjpeg-dev libtiff-dev \
- libgit2-dev curl lzip
+ libgit2-dev curl
@end example
@item Programs to view FITS images or tables
@@ -20559,29 +21247,12 @@ For example, you have taken ten exposures of your
scientific target, and you wou
When calling these operators you should determine how many operands they
should take in (unlike the rest of the operators that have a fixed number of
input operands).
As described in the first operand below, you do this through their first
popped operand (which should be a single integer number that is larger than
one).
+Below are Some important points for all the stacking operators described in
this section:
-@table @command
-
-@cindex NaN
-@item min
-For each pixel, find the minimum value in all given datasets.
-The output will have the same type as the input.
-
-The first popped operand to this operator must be a positive integer number
which specifies how many further operands should be popped from the stack.
-All the subsequently popped operands must have the same type and size.
-This operator (and all the variable-operand operators similar to it that are
discussed below) will work in multi-threaded mode unless Arithmetic is called
with the @option{--numthreads=1} option, see @ref{Multi-threaded operations}.
-
-Each pixel of the output of the @code{min} operator will be given the minimum
value of the same pixel from all the popped operands/images.
-For example, the following command will produce an image with the same size
and type as the three inputs, but each output pixel value will be the minimum
of the same pixel's values in all three input images.
-
-@example
-$ astarithmetic a.fits b.fits c.fits 3 min --output=min.fits
-@end example
-
-Important notes:
@itemize
@item
+@cindex NaN
NaN/blank pixels will be ignored, see @ref{Blank pixels}.
@item
@@ -20594,91 +21265,152 @@ The operation will be multi-threaded, greatly
speeding up the process if you hav
You can disable multi-threaded operations with the @option{--numthreads=1}
option (see @ref{Multi-threaded operations}).
@end itemize
-@item max
-For each pixel, find the maximum value in all given datasets.
-The output will have the same type as the input.
-This operator is called similar to the @command{min} operator, please see
there for more.
-For example
-@example
-$ astarithmetic a.fits b.fits c.fits 3 max -omax.fits
-@end example
+@table @command
+@item min
+@itemx max
+@itemx sum
+@itemx std
+@itemx mad
+@itemx mean
+@itemx median
+@itemx number
+For each pixel, calculate the respective statistic from in all given datasets.
+For the @code{min} and @code{max} operators, the output will have the same
type as the input. For the @code{number} operator, the output will have an
unsigned 32-bit integer type and the rest will be 32-bit floating point.
+
+The first popped operand to this operator must be a positive integer number
which specifies how many further operands should be popped from the stack.
+All the subsequently popped operands must have the same type and size.
+This operator (and all the variable-operand operators similar to it that are
discussed below) will work in multi-threaded mode unless Arithmetic is called
with the @option{--numthreads=1} option, see @ref{Multi-threaded operations}.
+
+For example, the following command will produce an image with the same size
and type as the three inputs, but each output pixel value will be the minimum
of the same pixel's values in all three input images.
-@item number
-For each pixel count the number of non-blank pixels in all given datasets.
-The output will be an unsigned 32-bit integer datatype (see @ref{Numeric data
types}).
-This operator is called similar to the @command{min} operator, please see
there for more.
-For example
@example
-$ astarithmetic a.fits b.fits c.fits 3 number -onum.fits
+$ astarithmetic a.fits b.fits c.fits 3 min --output=min.fits
@end example
-Some datasets may have blank values (which are also ignored in all similar
operators like @command{min}, @command{sum}, @command{mean} or
@command{median}).
+Regarding the @code{number} operator: some datasets may have blank values
(which are also ignored in all similar operators like @command{min},
@command{sum}, @command{mean} or @command{median}).
Hence, the final pixel values of this operator will not, in general, be equal
to the number of inputs.
This operator is therefore mostly called in parallel with those operators to
know the ``weight'' of each pixel (in case you want to only keep pixels that
had the full exposure for example).
-@item sum
-For each pixel, calculate the sum in all given datasets.
-The output will have the a single-precision (32-bit) floating point type.
-This operator is called similar to the @command{min} operator, please see
there for more.
-For example
+@item quantile
+For each pixel, find the quantile from all given datasets.
+The output will have the same numeric data type and size as the input datasets.
+Besides the input datasets, the quantile operator also needs a single
parameter (the requested quantile).
+The parameter should be the first popped operand, with a value between (and
including) 0 and 1.
+The second popped operand must be the number of datasets to use.
+
+In the example below, the first-popped operand (@command{0.7}) is the
quantile, the second-popped operand (@command{3}) is the number of datasets to
pop.
+
+@example
+astarithmetic a.fits b.fits c.fits 3 0.7 quantile
+@end example
+
+@item madclip-fill-mad
+@itemx madclip-fill-std
+@itemx madclip-fill-mean
+@itemx madclip-fill-median
+@itemx madclip-fill-number
+@cindex Median absolute deviation (MAD)
+@cindex MAD (Median absolute deviation)
+Find the respective statistic after median absolute deviation (MAD) filled
re-clipping of the values of the same element (pixel in an image) of all the
inputs.
+For a complete tutorial on clipping outliers see @ref{Clipping outliers} (if
you haven't read it yet, we encourage you to read through it before continuing).
+For the particular case of filled re-clipping (the @code{madclip-fill-*}
operators here) see @ref{Filled re-clipping}.
+Spoiler alert: this is currently the most robust stacking operator in Gnuastro.
+
+The output data type of all these operators is a 32-bit floating point number,
except for @code{madclip-fill-number}, where an unsigned 32-bit integer is
returned (see @ref{Numeric data types}).
+This operator will combine the given inputs into a single output of the same
dimension as one of the inputs.
+Each pixel of the output contains the requested statistic from the remaining
values after filled MAD re-clipping.
+This operator is very similar to @command{min}, with the exception that it
expects two extra operands (parameters for MAD-clipping) before the total
number of inputs.
+The first popped operand is the termination criteria and the second is the
multiple of the median absolute deviation.
+
+For example, in the command below, the first popped operand (@command{0.01})
is the MAD-clipping termination criteria.
+If the termination criteria is larger than, or equal to, 1 it is interpreted
as a pre-defined the number of clips.
+But if it is between 0 and 1, then it is the tolerance level on the change in
the median absolute devaition (see @ref{MAD clipping}).
+The second popped operand (@command{5}) is the multiple of the median absolute
deviation to use in MAD-clipping.
+The third popped operand (@command{3}) is number of datasets that will be used
(similar to the first popped operand to @command{min}).
+
+@example
+$ astarithmetic a.fits b.fits c.fits 3 5 0.01 madclip-fill-std
+@end example
+
+@item madclip-mad
+@itemx madclip-std
+@itemx madclip-mean
+@itemx madclip-median
+@itemx madclip-number
+Find the number of useful values after median absolute deviation (MAD)
clipping of the values of the same element (pixel in an image) of all the
inputs.
+These operators are called in an identical fashion to the
@code{madclip-fill-*} operators described above; see the description there for
more.
+
+@item sigclip-fill-number
+Find the number of useful values after filled sigma (@mymath{\sigma}, which
stands for the standard deviation) re-clipping of the values of the same
element (pixel in an image) of all the inputs.
+For a complete tutorial on clipping outliers see @ref{Clipping outliers} (if
you haven't read it yet, we encourage you to read through it before continuing).
+For the particular case of filled re-clipping (the @code{sigclip-fill-*}
operators here) see @ref{Filled re-clipping}.
+Spoiler alert: MAD filled re-clipping is currently most robust stacking
operator in Gnuastro (the @code{madclip-fill-*} operators).
+
+This operator will combine the specified number of inputs into a single output
that contains the number of remaining elements after @mymath{\sigma}-clipping
on each element/pixel (for more on @mymath{\sigma}-clipping, see @ref{Sigma
clipping}).
+This operator is very similar to @command{min}, with the exception that it
expects two operands (parameters for @mymath{\sigma}-clipping) before the total
number of inputs.
+The first popped operand is the termination criteria and the second is the
multiple of @mymath{\sigma}.
+
+For example, in the command below, the first popped operand (@command{0.2}) is
the sigma clipping termination criteria.
+If the termination criteria is larger than, or equal to, 1 it is interpreted
as the number of clips to do.
+But if it is between 0 and 1, then it is the tolerance level on the standard
deviation (see @ref{Sigma clipping}).
+The second popped operand (@command{5}) is the multiple of sigma to use in
@mymath{\sigma}-clipping.
+The third popped operand (@command{10}) is number of datasets that will be
used (similar to the first popped operand to @command{min}).
+
@example
-$ astarithmetic a.fits b.fits c.fits 3 sum -ostack-sum.fits
+astarithmetic a.fits b.fits c.fits 3 5 0.2 sigclip-fill-number
@end example
-@item mean
-For each pixel, calculate the mean in all given datasets.
+@item sigclip-fill-median
+For each pixel, find the @mymath{\sigma}-clipped median in all given datasets.
The output will have the a single-precision (32-bit) floating point type.
-This operator is called similar to the @command{min} operator, please see
there for more.
+This operator is called similar to the @command{sigclip-fill-number} operator,
please see there for more.
For example
@example
-$ astarithmetic a.fits b.fits c.fits 3 mean -ocoadd-mean.fits
+astarithmetic a.fits b.fits c.fits 3 5 0.2 sigclip-fill-median
@end example
-@item std
-For each pixel, find the standard deviation in all given datasets.
+@item sigclip-fill-mean
+For each pixel, find the @mymath{\sigma}-clipped mean in all given datasets.
The output will have the a single-precision (32-bit) floating point type.
-This operator is called similar to the @command{min} operator, please see
there for more.
+This operator is called similar to the @command{sigclip-fill-number} operator,
please see there for more.
For example
@example
-$ astarithmetic a.fits b.fits c.fits 3 std -ostd.fits
+astarithmetic a.fits b.fits c.fits 3 5 0.2 sigclip-fill-mean
@end example
-@item median
-For each pixel, find the median in all given datasets.
+@item sigclip-fill-std
+For each pixel, find the @mymath{\sigma}-clipped standard deviation in all
given datasets.
The output will have the a single-precision (32-bit) floating point type.
-This operator is called similar to the @command{min} operator, please see
there for more.
+This operator is called similar to the @command{sigclip-fill-number} operator,
please see there for more.
For example
@example
-$ astarithmetic a.fits b.fits c.fits 3 median \
- --output=stack-median.fits
+astarithmetic a.fits b.fits c.fits 3 5 0.2 sigclip-fill-std
@end example
-@item quantile
-For each pixel, find the quantile from all given datasets.
-The output will have the same numeric data type and size as the input datasets.
-Besides the input datasets, the quantile operator also needs a single
parameter (the requested quantile).
-The parameter should be the first popped operand, with a value between (and
including) 0 and 1.
-The second popped operand must be the number of datasets to use.
-
-In the example below, the first-popped operand (@command{0.7}) is the
quantile, the second-popped operand (@command{3}) is the number of datasets to
pop.
-
+@item sigclip-fill-mad
+For each pixel, find the @mymath{\sigma}-clipped median absolute deviation
(MAD) in all given datasets.
+The output will have the a single-precision (32-bit) floating point type.
+This operator is called similar to the @command{sigclip-fill-number} operator,
please see there for more.
+For example
@example
-astarithmetic a.fits b.fits c.fits 3 0.7 quantile
+astarithmetic a.fits b.fits c.fits 3 5 0.2 sigclip-fill-mad
@end example
@item sigclip-number
-For each pixel, find the sigma-clipped number (after removing outliers) in all
given datasets.
-The output will have the an unsigned 32-bit integer type (see @ref{Numeric
data types}).
+Find the number of useful values after sigma (@mymath{\sigma}, which stands
for the standard deviation) clipping of the values of the same element (pixel
in an image) of all the inputs.
+For a complete tutorial on clipping outliers see @ref{Clipping outliers} (if
you haven't read it yet, we encourage you to read through it before continuing).
+For the particular case of @mymath{\sigma}-clipping (the @code{sigclip-*}
operators here) see @ref{Sigma clipping}.
+Spoiler alert: MAD filled re-clipping is currently most robust stacking
operator in Gnuastro (the @code{madclip-fill-*} operators).
This operator will combine the specified number of inputs into a single output
that contains the number of remaining elements after @mymath{\sigma}-clipping
on each element/pixel (for more on @mymath{\sigma}-clipping, see @ref{Sigma
clipping}).
-This operator is very similar to @command{min}, with the exception that it
expects two operands (parameters for sigma-clipping) before the total number of
inputs.
+This operator is very similar to @command{min}, with the exception that it
expects two operands (parameters for @mymath{\sigma}-clipping) before the total
number of inputs.
The first popped operand is the termination criteria and the second is the
multiple of @mymath{\sigma}.
For example, in the command below, the first popped operand (@command{0.2}) is
the sigma clipping termination criteria.
If the termination criteria is larger than, or equal to, 1 it is interpreted
as the number of clips to do.
But if it is between 0 and 1, then it is the tolerance level on the standard
deviation (see @ref{Sigma clipping}).
-The second popped operand (@command{5}) is the multiple of sigma to use in
sigma-clipping.
+The second popped operand (@command{5}) is the multiple of sigma to use in
@mymath{\sigma}-clipping.
The third popped operand (@command{10}) is number of datasets that will be
used (similar to the first popped operand to @command{min}).
@example
@@ -20686,7 +21418,7 @@ astarithmetic a.fits b.fits c.fits 3 5 0.2
sigclip-number
@end example
@item sigclip-median
-For each pixel, find the sigma-clipped median in all given datasets.
+For each pixel, find the @mymath{\sigma}-clipped median in all given datasets.
The output will have the a single-precision (32-bit) floating point type.
This operator is called similar to the @command{sigclip-number} operator,
please see there for more.
For example
@@ -20695,7 +21427,7 @@ astarithmetic a.fits b.fits c.fits 3 5 0.2
sigclip-median
@end example
@item sigclip-mean
-For each pixel, find the sigma-clipped mean in all given datasets.
+For each pixel, find the @mymath{\sigma}-clipped mean in all given datasets.
The output will have the a single-precision (32-bit) floating point type.
This operator is called similar to the @command{sigclip-number} operator,
please see there for more.
For example
@@ -20704,15 +21436,27 @@ astarithmetic a.fits b.fits c.fits 3 5 0.2
sigclip-mean
@end example
@item sigclip-std
-For each pixel, find the sigma-clipped standard deviation in all given
datasets.
+For each pixel, find the @mymath{\sigma}-clipped standard deviation in all
given datasets.
The output will have the a single-precision (32-bit) floating point type.
This operator is called similar to the @command{sigclip-number} operator,
please see there for more.
For example
@example
astarithmetic a.fits b.fits c.fits 3 5 0.2 sigclip-std
@end example
+
+@item sigclip-mad
+For each pixel, find the @mymath{\sigma}-clipped median absolute deviation
(MAD) in all given datasets.
+The output will have the a single-precision (32-bit) floating point type.
+This operator is called similar to the @command{sigclip-number} operator,
please see there for more.
+For example
+@example
+astarithmetic a.fits b.fits c.fits 3 5 0.2 sigclip-mad
+@end example
@end table
+
+
+
@node Filtering operators, Pooling operators, Stacking operators, Arithmetic
operators
@subsubsection Filtering (smoothing) operators
@@ -21129,6 +21873,33 @@ $ astarithmetic seg.fits -hINPUT-NO-SKY seg.fits
-hOBJECTS \
263 ne nan where trim --output=obj-263.fits
@end example
+@item add-dimension-slow
+Build a higher-dimensional dataset from all the input datasets stacked after
one another (along the slowest dimension).
+The first popped operand has to be a single number.
+It is used by the operator to know how many operands it should pop from the
stack (and the size of the output in the new dimension).
+The rest of the operands must have the same size and numerical data type.
+This operator currently only works for 2D input operands, please contact us if
you want inputs to have different dimensions.
+
+The output's WCS (which should have a different dimensionality compared to the
inputs) can be read from another file with the @option{--wcsfile} option.
+If no file is specified for the WCS, the first dataset's WCS will be used, you
can later add/change the necessary WCS keywords with the FITS keyword
modification features of the Fits program (see @ref{Fits}).
+
+If your datasets do not have the same type, you can use the type
transformation operators of Arithmetic that are discussed below.
+Just beware of overflow if you are transforming to a smaller type, see
@ref{Numeric data types}.
+
+For example, let's assume you have 3 two-dimensional images @file{a.fits},
@file{b.fits} and @file{c.fits} (each with @mymath{200\times100} pixels).
+You can construct a 3D data cube with @mymath{200\times100\times3} voxels
(volume-pixels) using the command below:
+
+@example
+$ astarithmetic a.fits b.fits c.fits 3 add-dimension-slow
+@end example
+
+@item add-dimension-fast
+Similar to @code{add-dimension-slow} but along the fastest dimension.
+This operator currently only works for 1D input operands, please contact us if
you want inputs to have different dimensions.
+
+For example, let's assume you have 3 one-dimensional datasets, each with 100
elements.
+With this operator, you can construct a @mymath{3\times100} pixel FITS image
that has 3 pixels along the horizontal and 5 pixels along the vertical.
+
@item collapse-sum
Collapse the given dataset (second popped operand), by summing all elements
along the first popped operand (a dimension in FITS standard: counting from
one, from fastest dimension).
The returned dataset has one dimension less compared to the input.
@@ -21156,40 +21927,33 @@ The command below will collapse the whole third
dimension into a 2D array the si
$ astarithmetic cube.fits 3 collapse-sum float32
@end example
-@item collapse-mean
-Similar to @option{collapse-sum}, but the returned dataset will be the mean
value along the collapsed dimension, not the sum.
-
-@item collapse-number
-Similar to @option{collapse-sum}, but the returned dataset will be the number
of non-blank values along the collapsed dimension.
-The output will have a 32-bit signed integer type.
-If the input dataset does not have blank values, all the elements in the
returned dataset will have a single value (the length of the collapsed
dimension).
-Therefore this is mostly relevant when there are blank values in the dataset.
-
-@item collapse-min
-Similar to @option{collapse-sum}, but the returned dataset will have the same
numeric type as the input and will contain the minimum value for each pixel
along the collapsed dimension.
+@item collapse-min
+@itemx collapse-max
+@itemx collapse-mean
+@itemx collapse-median
+Similar to @option{collapse-sum}, but the returned dataset will be the desired
statistic along the collapsed dimension, not the sum.
-@item collapse-max
-Similar to @option{collapse-sum}, but the returned dataset will have the same
numeric type as the input and will contain the maximum value for each pixel
along the collapsed dimension.
+@item collapse-madclip-fill-mad
+@itemx collapse-madclip-fill-std
+@itemx collapse-madclip-fill-mean
+@itemx collapse-madclip-fill-median
+@itemx collapse-madclip-fill-number
+Collapse the input dataset (fourth popped operand) along the FITS dimension
given as the first popped operand by calculating the desired statistic after
median absolute deviation (MAD) filled re-cliping.
+The MAD-clipping parameters (namely, the multiple of sigma and termination
criteria) are read as the third and second popped operands respectively.
-@item collapse-median
-Similar to @option{collapse-sum}, but the returned dataset will have the same
numeric type as the input and will contain the median value for each pixel
along the collapsed dimension.
-
-The median involves sorting, therefore @code{collapse-median} will do each
calculation in different CPU threads to speed up the operation.
-By default, Arithmetic will detect and use all available threads, but you can
override this with the @option{--numthreads} (or @option{-N}) option.
-
-@item collapse-sigclip-mean
-Collapse the input dataset (fourth popped operand) along the FITS dimension
given as the first popped operand by calculating the sigma-clipped mean.
-The sigma-clipping parameters (namely, the multiple of sigma and termination
criteria) are read as the third and second popped operands respectively.
-For more on sigma-clipping, see @ref{Sigma clipping}.
+This is the most robust method to reject outliers; for more on filled
re-clipping and its advantages, see @ref{Filled re-clipping}.
+For a more general tutorial on rejecting outliers, see @ref{Clipping outliers}.
+If you have not done this tutorial yet, we recommend you to take an hour or so
and go through that tutorial for optimal understanding and results.
For example, with the command below, the pixels of the input 2 dimensional
@file{image.fits} will be collapsed to a single dimension output.
The first popped operand is @code{2}, so it will collapse all the pixels that
are vertically on top of each other.
Such that the output will have the same number of pixels as the horizontal
axis of the input.
During the collapsing, all pixels that are more than @mymath{3\sigma} (third
popped operand) are rejected, and the clipping will continue until the standard
deviation changes less than @mymath{0.2} between clips.
+Finally the @code{counter} operator is used to have a two-column table with
the first one being a simple counter starting from one (see @ref{Size and
position operators}).
@example
$ astarithmetic image.fits 3 0.2 2 collapse-sigclip-mean \
- --output=collapsed-vertical.fits
+ counter --output=collapsed-vertical.fits
@end example
@cartouche
@@ -21202,44 +21966,29 @@ The FITS format is always superior (since it stores
the value in binary, therefo
But if you are forced to save the output in plain-text, use the @code{float64}
operator after this to change the type to 64-bit floating point (which will
print more decimals).
@end cartouche
-@item collapse-sigclip-std
-Collapse the input dataset along the given FITS dimension by calculating the
sigma-clipped standard deviation.
-Except for returning the standard deviation after clipping, this function is
similar to @code{collapse-sigclip-mean}, see the description of that operator
for more.
-
-@item collapse-sigclip-median
-Collapse the input dataset along the given FITS dimension by calculating the
sigma-clipped median.
-Except for returning the median after clipping, this function is similar to
@code{collapse-sigclip-mean}, see the description of that operator for more.
-
-@item collapse-sigclip-number
-Collapse the input dataset along the given FITS dimension by calculating the
number of elements that remain after sigma-clipped.
-Except for returning the number after clipping, this function is similar to
@code{collapse-sigclip-mean}, see the description of that operator for more.
-
-@item add-dimension-slow
-Build a higher-dimensional dataset from all the input datasets stacked after
one another (along the slowest dimension).
-The first popped operand has to be a single number.
-It is used by the operator to know how many operands it should pop from the
stack (and the size of the output in the new dimension).
-The rest of the operands must have the same size and numerical data type.
-This operator currently only works for 2D input operands, please contact us if
you want inputs to have different dimensions.
-
-The output's WCS (which should have a different dimensionality compared to the
inputs) can be read from another file with the @option{--wcsfile} option.
-If no file is specified for the WCS, the first dataset's WCS will be used, you
can later add/change the necessary WCS keywords with the FITS keyword
modification features of the Fits program (see @ref{Fits}).
-
-If your datasets do not have the same type, you can use the type
transformation operators of Arithmetic that are discussed below.
-Just beware of overflow if you are transforming to a smaller type, see
@ref{Numeric data types}.
-
-For example, let's assume you have 3 two-dimensional images @file{a.fits},
@file{b.fits} and @file{c.fits} (each with @mymath{200\times100} pixels).
-You can construct a 3D data cube with @mymath{200\times100\times3} voxels
(volume-pixels) using the command below:
-
-@example
-$ astarithmetic a.fits b.fits c.fits 3 add-dimension-slow
-@end example
-
-@item add-dimension-fast
-Similar to @code{add-dimension-slow} but along the fastest dimension.
-This operator currently only works for 1D input operands, please contact us if
you want inputs to have different dimensions.
-
-For example, let's assume you have 3 one-dimensional datasets, each with 100
elements.
-With this operator, you can construct a @mymath{3\times100} pixel FITS image
that has 3 pixels along the horizontal and 5 pixels along the vertical.
+@item collapse-madclip-mad
+@itemx collapse-madclip-std
+@itemx collapse-madclip-mean
+@itemx collapse-madclip-median
+@itemx collapse-madclip-number
+Collapse the input dataset (fourth popped operand) along the FITS dimension
given as the first popped operand by calculating the desired statistic after
median absolute deviation (MAD) cliping.
+This operator is called similarly to the @code{collapse-madclip-fill-*}
operators, see the description there for more.
+
+@item collapse-sigclip-fill-mad
+@itemx collapse-sigclip-fill-std
+@itemx collapse-sigclip-fill-mean
+@itemx collapse-sigclip-fill-median
+@itemx collapse-sigclip-fill-number
+Collapse the input dataset (fourth popped operand) along the FITS dimension
given as the first popped operand by calculating the desired statistic after
filled @mymath{\sigma} re-clipping.
+This operator is called similarly to the @code{collapse-madclip-fill-*}
operators, see the description there for more.
+
+@item collapse-sigclip-mad
+@itemx collapse-sigclip-std
+@itemx collapse-sigclip-mean
+@itemx collapse-sigclip-median
+@itemx collapse-sigclip-number
+Collapse the input dataset (fourth popped operand) along the FITS dimension
given as the first popped operand by calculating the desired statistic after
@mymath{\sigma}-cliping.
+This operator is called similarly to the @code{collapse-madclip-fill-*}
operators, see the description there for more.
@end table
@node Conditional operators, Mathematical morphology operators, Dimensionality
changing operators, Arithmetic operators
@@ -24616,7 +25365,6 @@ The Statistics program is designed for such situations.
@menu
* Histogram and Cumulative Frequency Plot:: Basic definitions.
* 2D Histograms:: Plotting the distribution of two variables.
-* Sigma clipping:: Definition of @mymath{\sigma}-clipping.
* Least squares fitting:: Fitting with various parametric functions.
* Sky value:: Definition and derivation of the Sky value.
* Invoking aststatistics:: Arguments and options to Statistics.
@@ -24662,7 +25410,7 @@ Formally, an interval/bin between a and b is
represented by [a, b).
When the over-all range of the dataset is specified (with the
@option{--greaterequal}, @option{--lessthan}, or @option{--qrange} options),
the acceptable values of the dataset are also defined with a similar
inclusive-exclusive manner.
But when the range is determined from the actual dataset (none of these
options is called), the last element in the dataset is included in the last
bin's count.
-@node 2D Histograms, Sigma clipping, Histogram and Cumulative Frequency Plot,
Statistics
+@node 2D Histograms, Least squares fitting, Histogram and Cumulative Frequency
Plot, Statistics
@subsection 2D Histograms
@cindex 2D histogram
@cindex Histogram, 2D
@@ -24917,78 +25665,12 @@ $ pdflatex cmd-report.tex
The improved quality, blending in with the text, vector-graphics resolution
and other features make this plot pleasing to the eye, and let your readers
focus on the main point of your scientific argument.
PGFPlots can also built the PDF of the plot separately from the rest of the
paper/report, see @ref{2D histogram as a table for plotting} for the necessary
changes in the preamble.
-@node Sigma clipping, Least squares fitting, 2D Histograms, Statistics
-@subsection Sigma clipping
-
-Let's assume that you have pure noise (centered on zero) with a clear
@url{https://en.wikipedia.org/wiki/Normal_distribution,Gaussian distribution},
or see @ref{Photon counting noise}.
-Now let's assume you add very bright objects (signal) on the image which have
a very sharp boundary.
-By a sharp boundary, we mean that there is a clear cutoff (from the noise) at
the pixels the objects finish.
-In other words, at their boundaries, the objects do not fade away into the
noise.
-In such a case, when you plot the histogram (see @ref{Histogram and Cumulative
Frequency Plot}) of the distribution, the pixels relating to those objects will
be clearly separate from pixels that belong to parts of the image that did not
have any signal (were just noise).
-In the cumulative frequency plot, after a steady rise (due to the noise), you
would observe a long flat region were for a certain range of data (horizontal
axis), there is no increase in the index (vertical axis).
-@cindex Blurring
-@cindex Cosmic rays
-@cindex Aperture blurring
-@cindex Atmosphere blurring
-Outliers like the example above can significantly bias the measurement of
noise statistics.
-@mymath{\sigma}-clipping is defined as a way to avoid the effect of such
outliers.
-In astronomical applications, cosmic rays (when they collide at a near normal
incidence angle) are a very good example of such outliers.
-The tracks they leave behind in the image are perfectly immune to the blurring
caused by the atmosphere and the aperture.
-They are also very energetic and so their borders are usually clearly
separated from the surrounding noise.
-So @mymath{\sigma}-clipping is very useful in removing their effect on the
data.
-See Figure 15 in Akhlaghi and Ichikawa,
@url{https://arxiv.org/abs/1505.01664,2015}.
-@mymath{\sigma}-clipping is defined as the very simple iteration below.
-In each iteration, the range of input data might decrease and so when the
outliers have the conditions above, the outliers will be removed through this
iteration.
-The exit criteria will be discussed below.
-@enumerate
-@item
-Calculate the standard deviation (@mymath{\sigma}) and median (@mymath{m})
-of a distribution.
-@item
-Remove all points that are smaller or larger than
-@mymath{m\pm\alpha\sigma}.
-@item
-Go back to step 1, unless the selected exit criteria is reached.
-@end enumerate
-@noindent
-The reason the median is used as a reference and not the mean is that the mean
is too significantly affected by the presence of outliers, while the median is
less affected, see @ref{Quantifying signal in a tile}.
-As you can tell from this algorithm, besides the condition above (that the
signal have clear high signal to noise boundaries) @mymath{\sigma}-clipping is
only useful when the signal does not cover more than half of the full data set.
-If they do, then the median will lie over the outliers and
@mymath{\sigma}-clipping might remove the pixels with no signal.
-There are commonly two exit criteria to stop the @mymath{\sigma}-clipping
-iteration:
-
-@itemize
-@item
-When a certain number of iterations has taken place (second value to the
@option{--sclipparams} option is larger than 1).
-@item
-When the new measured standard deviation is within a certain tolerance level
of the old one (second value to the @option{--sclipparams} option is less than
1).
-The tolerance level is defined by:
-
-@dispmath{\sigma_{old}-\sigma_{new} \over \sigma_{new}}
-
-The standard deviation is used because it is heavily influenced by the
presence of outliers.
-Therefore the fact that it stops changing between two iterations is a sign
that we have successfully removed outliers.
-Note that in each clipping, the dispersion in the distribution is either less
or equal.
-So @mymath{\sigma_{old}\geq\sigma_{new}}.
-@end itemize
-
-@cartouche
-@noindent
-When working on astronomical images, objects like galaxies and stars are
blurred by the atmosphere and the telescope aperture, therefore their signal
sinks into the noise very gradually.
-Galaxies in particular do not appear to have a clear high signal to noise
cutoff at all.
-Therefore @mymath{\sigma}-clipping will not be useful in removing their effect
on the data.
-
-To gauge if @mymath{\sigma}-clipping will be useful for your dataset, look at
the histogram (see @ref{Histogram and Cumulative Frequency Plot}).
-The ASCII histogram that is printed on the command-line with
@option{--asciihist} is good enough in most cases.
-@end cartouche
-
-
-@node Least squares fitting, Sky value, Sigma clipping, Statistics
+@node Least squares fitting, Sky value, 2D Histograms, Statistics
@subsection Least squares fitting
@cindex Radial profile
@@ -25692,6 +26374,9 @@ Print the mean (average) of all used elements.
@itemx --std
Print the standard deviation of all used elements.
+@item --mad
+Print the median absolute deviation (MAD) of all used elements.
+
@item -E
@itemx --median
Print the median of all used elements.
@@ -25760,13 +26445,13 @@ distributions are no longer symmetric, see
@option{--mode} and Appendix C
of @url{https://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa 2015} for
more.
-@item --sigclip-number
-Number of elements after applying @mymath{\sigma}-clipping (see @ref{Sigma
clipping}).
-@mymath{\sigma}-clipping configuration is done with the
@option{--sigclipparams} option.
-
-@item --sigclip-median
-Median after applying @mymath{\sigma}-clipping (see @ref{Sigma clipping}).
-@mymath{\sigma}-clipping configuration is done with the
@option{--sigclipparams} option.
+@item --sigclip-std
+@item --sigclip-mad
+@itemx --sigclip-mean
+@itemx --sigclip-number
+@itemx --sigclip-median
+Calculate the desired statistic after applying @mymath{\sigma}-clipping (see
@ref{Sigma clipping}, part of the tutorial @ref{Clipping outliers}).
+@mymath{\sigma}-clipping configuration is done with the @option{--sclipparams}
option.
@cindex Outlier
Here is one scenario where this can be useful: assume you have a table and you
would like to remove the rows that are outliers (not within the
@mymath{\sigma}-clipping range).
@@ -25812,15 +26497,15 @@ $ asttable table.fits --range=COLUMN,$r
To save the resulting table (that is clean of outliers) in another file (for
example, named @file{cleaned.fits}, it can also have a @file{.txt} suffix),
just add @option{--output=cleaned.fits} to the command above.
+@item --madclip-std
+@item --madclip-mad
+@itemx --madclip-mean
+@itemx --madclip-number
+@itemx --madclip-median
+Calculate the desired statistic after applying median absolute deviation (MAD)
clipping (see @ref{MAD clipping}, part of the tutorial @ref{Clipping outliers}).
+MAD-clipping configuration is done with the @option{--mclipparams} option.
-@item --sigclip-mean
-Mean after applying @mymath{\sigma}-clipping (see @ref{Sigma clipping}).
-@mymath{\sigma}-clipping configuration is done with the
@option{--sigclipparams} option.
-
-@item --sigclip-std
-Standard deviation after applying @mymath{\sigma}-clipping (see @ref{Sigma
clipping}).
-@mymath{\sigma}-clipping configuration is done with the
@option{--sigclipparams} option.
-
+This option behaves similarly to @option{--sigclip-*} options, read their
description for usage examples.
@end table
@node Generating histograms and cumulative frequency plots, Fitting options,
Single value measurements, Invoking aststatistics
@@ -25909,10 +26594,15 @@ The third will be number of input points that fall
within that box.
@itemx --cumulative
Save the cumulative frequency plot of the usable values in the input dataset
into a table, similar to @option{--histogram}.
+@item --madclip
+Do median absolute deviation (MAD) clipping on the usable pixels of the input
dataset.
+See @ref{MAD clipping} for a description on MAD-clipping and @ref{Clipping
outliers} for a complete tutorial on clipping of outliers.
+The MAD-clipping parameters can be set through the @option{--mclipparams}
option (see below).
+
@item -s
@itemx --sigmaclip
Do @mymath{\sigma}-clipping on the usable pixels of the input dataset.
-See @ref{Sigma clipping} for a full description on @mymath{\sigma}-clipping
and also to better understand this option.
+See @ref{Sigma clipping} for a full description on @mymath{\sigma}-clipping
and @ref{Clipping outliers} for a complete tutorial on clipping of outliers.
The @mymath{\sigma}-clipping parameters can be set through the
@option{--sclipparams} option (see below).
@item --mirror=FLT
@@ -26242,6 +26932,10 @@ 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 --mclipparams=FLT,FLT
+The MAD-clipping parameters.
+This is very similar to @option{--sclipparams} above, see there for more.
+
@item --outliersclip=FLT,FLT
@mymath{\sigma}-clipping parameters for the outlier rejection of the Sky
value (similar to @option{--sclipparams}).
@@ -36234,6 +36928,21 @@ For more on @code{minmapsize} and @code{quietmmap} see
@ref{Memory management}.
For more on sigma clipping, see @ref{Sigma clipping}.
@end deftypefun
+@deftypefun {gal_data_t *} gal_dimension_collapse_sclip_fill_std (gal_data_t
@code{*in}, size_t @code{c_dim}, float @code{multip}, float @code{param},
size_t @code{numthreads}, size_t @code{minmapsize}, int @code{quietmmap})
+Similar to @code{gal_dimension_collapse_sclip_std}, but with filled
re-clipping (see @ref{Filled re-clipping}).
+@end deftypefun
+
+@deftypefun {gal_data_t *} gal_dimension_collapse_sclip_mad (gal_data_t
@code{*in}, size_t @code{c_dim}, float @code{multip}, float @code{param},
size_t @code{numthreads}, size_t @code{minmapsize}, int @code{quietmmap})
+Collapse the input dataset (@code{in}) along the given dimension
(@code{c_dim}, in C definition: starting from zero, from the slowest
dimension), by finding the median absolute deviation (MAD) of pixels along that
dimension after sigma-clipping.
+Since sigma-clipping involves sorting, this operator benefits from many
threads (which needs to be set with @code{numthreads}).
+For more on @code{minmapsize} and @code{quietmmap} see @ref{Memory management}.
+For more on sigma clipping, see @ref{Sigma clipping}.
+@end deftypefun
+
+@deftypefun {gal_data_t *} gal_dimension_collapse_sclip_fill_mad (gal_data_t
@code{*in}, size_t @code{c_dim}, float @code{multip}, float @code{param},
size_t @code{numthreads}, size_t @code{minmapsize}, int @code{quietmmap})
+Similar to @code{gal_dimension_collapse_sclip_mad}, but with filled
re-clipping (see @ref{Filled re-clipping}).
+@end deftypefun
+
@deftypefun {gal_data_t *} gal_dimension_collapse_sclip_mean (gal_data_t
@code{*in}, size_t @code{c_dim}, float @code{multip}, float @code{param},
size_t @code{numthreads}, size_t @code{minmapsize}, int @code{quietmmap})
Collapse the input dataset (@code{in}) along the given dimension
(@code{c_dim}, in C definition: starting from zero, from the slowest
dimension), by finding the mean of pixels along that dimension after
sigma-clipping.
Since sigma-clipping involves sorting, this operator benefits from many
threads (which needs to be set with @code{numthreads}).
@@ -36241,6 +36950,10 @@ For more on @code{minmapsize} and @code{quietmmap} see
@ref{Memory management}.
For more on sigma clipping, see @ref{Sigma clipping}.
@end deftypefun
+@deftypefun {gal_data_t *} gal_dimension_collapse_sclip_fill_mean (gal_data_t
@code{*in}, size_t @code{c_dim}, float @code{multip}, float @code{param},
size_t @code{numthreads}, size_t @code{minmapsize}, int @code{quietmmap})
+Similar to @code{gal_dimension_collapse_sclip_mean}, but with filled
re-clipping (see @ref{Filled re-clipping}).
+@end deftypefun
+
@deftypefun {gal_data_t *} gal_dimension_collapse_sclip_median (gal_data_t
@code{*in}, size_t @code{c_dim}, float @code{multip}, float @code{param},
size_t @code{numthreads}, size_t @code{minmapsize}, int @code{quietmmap})
Collapse the input dataset (@code{in}) along the given dimension
(@code{c_dim}, in C definition: starting from zero, from the slowest
dimension), by finding the median of pixels along that dimension after
sigma-clipping.
Since sigma-clipping involves sorting, this operator benefits from many
threads (which needs to be set with @code{numthreads}).
@@ -36248,6 +36961,10 @@ For more on @code{minmapsize} and @code{quietmmap} see
@ref{Memory management}.
For more on sigma clipping, see @ref{Sigma clipping}.
@end deftypefun
+@deftypefun {gal_data_t *} gal_dimension_collapse_sclip_fill_median
(gal_data_t @code{*in}, size_t @code{c_dim}, float @code{multip}, float
@code{param}, size_t @code{numthreads}, size_t @code{minmapsize}, int
@code{quietmmap})
+Similar to @code{gal_dimension_collapse_sclip_median}, but with filled
re-clipping (see @ref{Filled re-clipping}).
+@end deftypefun
+
@deftypefun {gal_data_t *} gal_dimension_collapse_sclip_number (gal_data_t
@code{*in}, size_t @code{c_dim}, float @code{multip}, float @code{param},
size_t @code{numthreads}, size_t @code{minmapsize}, int @code{quietmmap})
Collapse the input dataset (@code{in}) along the given dimension
(@code{c_dim}, in C definition: starting from zero, from the slowest
dimension), by finding the number of pixels along that dimension after
sigma-clipping.
Since sigma-clipping involves sorting, this operator benefits from many
threads (which needs to be set with @code{numthreads}).
@@ -36255,6 +36972,65 @@ For more on @code{minmapsize} and @code{quietmmap} see
@ref{Memory management}.
For more on sigma clipping, see @ref{Sigma clipping}.
@end deftypefun
+@deftypefun {gal_data_t *} gal_dimension_collapse_sclip_fill_number
(gal_data_t @code{*in}, size_t @code{c_dim}, float @code{multip}, float
@code{param}, size_t @code{numthreads}, size_t @code{minmapsize}, int
@code{quietmmap})
+Similar to @code{gal_dimension_collapse_sclip_number}, but with filled
re-clipping (see @ref{Filled re-clipping}).
+@end deftypefun
+
+@deftypefun {gal_data_t *} gal_dimension_collapse_mclip_std (gal_data_t
@code{*in}, size_t @code{c_dim}, float @code{multip}, float @code{param},
size_t @code{numthreads}, size_t @code{minmapsize}, int @code{quietmmap})
+Collapse the input dataset (@code{in}) along the given dimension
(@code{c_dim}, in C definition: starting from zero, from the slowest
dimension), by finding the standard deviation of pixels along that dimension
after median absolute deviation (MAD) clipping.
+Since MAD-clipping involves sorting, this operator benefits from many threads
(which needs to be set with @code{numthreads}).
+For more on @code{minmapsize} and @code{quietmmap} see @ref{Memory management}.
+For more on MAD-clipping, see @ref{MAD clipping}.
+@end deftypefun
+
+@deftypefun {gal_data_t *} gal_dimension_collapse_mclip_fill_std (gal_data_t
@code{*in}, size_t @code{c_dim}, float @code{multip}, float @code{param},
size_t @code{numthreads}, size_t @code{minmapsize}, int @code{quietmmap})
+Similar to @code{gal_dimension_collapse_mclip_std}, but with filled
re-clipping (see @ref{Filled re-clipping}).
+@end deftypefun
+
+@deftypefun {gal_data_t *} gal_dimension_collapse_mclip_mad (gal_data_t
@code{*in}, size_t @code{c_dim}, float @code{multip}, float @code{param},
size_t @code{numthreads}, size_t @code{minmapsize}, int @code{quietmmap})
+Collapse the input dataset (@code{in}) along the given dimension
(@code{c_dim}, in C definition: starting from zero, from the slowest
dimension), by finding the median absolute deviation (MAD) of pixels along that
dimension after median absolute deviation (MAD) clipping.
+Since MAD-clipping involves sorting, this operator benefits from many threads
(which needs to be set with @code{numthreads}).
+For more on @code{minmapsize} and @code{quietmmap} see @ref{Memory management}.
+For more on MAD-clipping, see @ref{MAD clipping}.
+@end deftypefun
+
+@deftypefun {gal_data_t *} gal_dimension_collapse_mclip_fill_mad (gal_data_t
@code{*in}, size_t @code{c_dim}, float @code{multip}, float @code{param},
size_t @code{numthreads}, size_t @code{minmapsize}, int @code{quietmmap})
+Similar to @code{gal_dimension_collapse_mclip_mad}, but with filled
re-clipping (see @ref{Filled re-clipping}).
+@end deftypefun
+
+@deftypefun {gal_data_t *} gal_dimension_collapse_mclip_mean (gal_data_t
@code{*in}, size_t @code{c_dim}, float @code{multip}, float @code{param},
size_t @code{numthreads}, size_t @code{minmapsize}, int @code{quietmmap})
+Collapse the input dataset (@code{in}) along the given dimension
(@code{c_dim}, in C definition: starting from zero, from the slowest
dimension), by finding the mean of pixels along that dimension after median
absolute deviation (MAD) clipping.
+Since MAD-clipping involves sorting, this operator benefits from many threads
(which needs to be set with @code{numthreads}).
+For more on @code{minmapsize} and @code{quietmmap} see @ref{Memory management}.
+For more on MAD-clipping, see @ref{MAD clipping}.
+@end deftypefun
+
+@deftypefun {gal_data_t *} gal_dimension_collapse_mclip_fill_mean (gal_data_t
@code{*in}, size_t @code{c_dim}, float @code{multip}, float @code{param},
size_t @code{numthreads}, size_t @code{minmapsize}, int @code{quietmmap})
+Similar to @code{gal_dimension_collapse_mclip_mean}, but with filled
re-clipping (see @ref{Filled re-clipping}).
+@end deftypefun
+
+@deftypefun {gal_data_t *} gal_dimension_collapse_mclip_median (gal_data_t
@code{*in}, size_t @code{c_dim}, float @code{multip}, float @code{param},
size_t @code{numthreads}, size_t @code{minmapsize}, int @code{quietmmap})
+Collapse the input dataset (@code{in}) along the given dimension
(@code{c_dim}, in C definition: starting from zero, from the slowest
dimension), by finding the median of pixels along that dimension after median
absolute deviation (MAD) clipping.
+Since MAD-clipping involves sorting, this operator benefits from many threads
(which needs to be set with @code{numthreads}).
+For more on @code{minmapsize} and @code{quietmmap} see @ref{Memory management}.
+For more on MAD-clipping, see @ref{MAD clipping}.
+@end deftypefun
+
+@deftypefun {gal_data_t *} gal_dimension_collapse_mclip_fill_median
(gal_data_t @code{*in}, size_t @code{c_dim}, float @code{multip}, float
@code{param}, size_t @code{numthreads}, size_t @code{minmapsize}, int
@code{quietmmap})
+Similar to @code{gal_dimension_collapse_mclip_median}, but with filled
re-clipping (see @ref{Filled re-clipping}).
+@end deftypefun
+
+@deftypefun {gal_data_t *} gal_dimension_collapse_mclip_number (gal_data_t
@code{*in}, size_t @code{c_dim}, float @code{multip}, float @code{param},
size_t @code{numthreads}, size_t @code{minmapsize}, int @code{quietmmap})
+Collapse the input dataset (@code{in}) along the given dimension
(@code{c_dim}, in C definition: starting from zero, from the slowest
dimension), by finding the number of pixels along that dimension after median
absolute deviation (MAD) clipping.
+Since MAD-clipping involves sorting, this operator benefits from many threads
(which needs to be set with @code{numthreads}).
+For more on @code{minmapsize} and @code{quietmmap} see @ref{Memory management}.
+For more on MAD-clipping, see @ref{MAD clipping}.
+@end deftypefun
+
+@deftypefun {gal_data_t *} gal_dimension_collapse_mclip_fill_number
(gal_data_t @code{*in}, size_t @code{c_dim}, float @code{multip}, float
@code{param}, size_t @code{numthreads}, size_t @code{minmapsize}, int
@code{quietmmap})
+Similar to @code{gal_dimension_collapse_mclip_number}, but with filled
re-clipping (see @ref{Filled re-clipping}).
+@end deftypefun
+
@deftypefun size_t gal_dimension_remove_extra (size_t @code{ndim}, size_t
@code{*dsize}, struct wcsprm @code{*wcs})
Remove extra dimensions (those that only have a length of 1) from the basic
size information of a dataset.
@code{ndim} is the number of dimensions and @code{dsize} is an array with
@code{ndim} elements containing the size along each dimension in the C
dimension order.
@@ -40530,6 +41306,21 @@ If the symmetricity of the derived mode is less than
this value, all the returne
Macros used to identify if the regularity of the bins when defining bins.
@end deffn
+@deffn Macro GAL_STATISTICS_CLIP_OUTCOL_STD
+@deffnx Macro GAL_STATISTICS_CLIP_OUTCOL_MAD
+@deffnx Macro GAL_STATISTICS_CLIP_OUTCOL_MEAN
+@deffnx Macro GAL_STATISTICS_CLIP_OUTCOL_MEDIAN
+@deffnx Macro GAL_STATISTICS_CLIP_OUTCOL_NUMBER_USED
+@deffnx Macro GAL_STATISTICS_CLIP_OUTCOL_NUMBER_CLIPS
+Macros containing the index of the clipping outputs, see the descriptions of
@code{gal_statistics_clip_sigma} below.
+@end deffn
+
+@deffn Macro GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_STD
+@deffnx Macro GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_MAD
+@deffnx Macro GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_MEAN
+Macros containing bit flags for optional clipping outputs, see the
descriptions of @code{gal_statistics_clip_sigma} below.
+@end deffn
+
@cindex Number
@deftypefun {gal_data_t *} gal_statistics_number (gal_data_t @code{*input})
Return a single-element dataset with type @code{size_t} which contains the
number of non-blank elements in @code{input}.
@@ -40591,6 +41382,24 @@ In this case, the sorting and removal of blank
elements will be done directly on
However, after this function the original dataset may have changed (if it was
not sorted or had blank values).
@end deftypefun
+@cindex Median absolute deviation (MAD)
+@cindex MAD (Median absolute deviation)
+@deftypefun {gal_data_t *} gal_statistics_mad (gal_data_t @code{*input}, int
@code{inplace})
+Return a single-element dataset with same type as input, containing the median
absolute deviation (MAD) of the non-blank values in @code{input}.
+
+If @code{inplace==0}, the input dataset will remain untouched.
+Otherwise, the MAD calculation will be done on the input dataset without
allocating a new one (its values will be changed after this function).
+This is good when you do not need the input after this function and avoid
taking extra RAM and CPU.
+@end deftypefun
+
+@deftypefun {gal_data_t *} gal_statistics_median_mad (gal_data_t
@code{*input}, int @code{inplace})
+Return a two-element dataset with same type as input, containing the median
and median absolute deviation (MAD) of the non-blank values in @code{input}.
+
+If @code{inplace==0}, the input dataset will remain untouched.
+Otherwise, the MAD calculation will be done on the input dataset without
allocating a new one (its values will be changed after this function).
+This is good when you do not need the input after this function and avoid
taking extra RAM and CPU.
+@end deftypefun
+
@cindex Quantile
@deftypefun size_t gal_statistics_quantile_index (size_t @code{size}, double
@code{quantile})
Return the index of the element that has a quantile of @code{quantile}
@@ -40787,7 +41596,7 @@ When a histogram is given and it is normalized, the CFP
will also be normalized
@end deftypefun
-@deftypefun {gal_data_t *} gal_statistics_sigma_clip (gal_data_t
@code{*input}, float @code{multip}, float @code{param}, int @code{inplace}, int
@code{quiet})
+@deftypefun {gal_data_t *} gal_statistics_clip_sigma (gal_data_t
@code{*input}, float @code{multip}, float @code{param}, float
@code{extrastats}, int @code{inplace}, int @code{quiet})
Apply @mymath{\sigma}-clipping on a given dataset and return a dataset that
contains the results.
For a description of @mymath{\sigma}-clipping see @ref{Sigma clipping}.
@code{multip} is the multiple of the standard deviation (or @mymath{\sigma},
that is used to define outliers in each round of clipping).
@@ -40796,21 +41605,43 @@ The role of @code{param} is determined based on its
value.
If @code{param} is larger than @code{1} (one), it must be an integer and will
be interpreted as the number clips to do.
If it is less than @code{1} (one), it is interpreted as the tolerance level to
stop the iteration.
-The returned dataset (let's call it @code{out}) contains a four-element array
with type @code{GAL_TYPE_FLOAT32}.
-The final number of clips is stored in the @code{out->status}.
+The returned dataset (let's call it @code{out}) contains a 6-element array
with type @code{GAL_TYPE_FLOAT32}.
+Through the @code{GAL_STATISTICS_CLIP_OUTCOL_*} macros below, you can access
any particular measurement.
+
@example
+out=gal_statistics_clip_sigma(input, ....);
float *array=out->array;
-array[0]: Number of points used.
-array[1]: Median.
-array[2]: Mean.
-array[3]: Standard deviation.
+
+array[ GAL_STATISTICS_CLIP_OUTCOL_NUMBER_USED ]
+array[ GAL_STATISTICS_CLIP_OUTCOL_MEAN ]
+array[ GAL_STATISTICS_CLIP_OUTCOL_STD ]
+array[ GAL_STATISTICS_CLIP_OUTCOL_MEDIAN ]
+array[ GAL_STATISTICS_CLIP_OUTCOL_MAD ]
+array[ GAL_STATISTICS_CLIP_OUTCOL_NUMBER_CLIPS ]
+@end example
+
+However, note that all are not measured by default!
+Since the mean and MAD are not necessary during sigma-clipping, if you want
them, you have to set the following two bit flags in the @code{extrastats}
argument as below.
+
+@example
+int extrastats=0; /* To initialize all bits */
+
+/* If you want the sigma-clipped MAD. */
+extrastats |= GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_MAD;
+
+/* If you want the sigma-clipped mean. */
+extrastats |= GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_MEAN;
@end example
-If the @mymath{\sigma}-clipping does not converge or all input elements are
-blank, then this function will return NaN values for all the elements
-above.
+If the @mymath{\sigma}-clipping does not converge or all input elements are
blank, then this function will return NaN values for all the elements above.
@end deftypefun
+@deftypefun {gal_data_t *} gal_statistics_clip_mad (gal_data_t @code{*input},
float @code{multip}, float @code{param}, uint8_t @code{extrastats}, int
@code{inplace}, int @code{quiet})
+Similar to @code{gal_statistics_clip_sigma}, but will do median absolute
deviation (MAD) based clipping, see @ref{MAD clipping}.
+
+The only difference is that for this function the MAD is automatically
calculated during clipping.
+It is the mean and standard deviation that will not be calculated unless
requested with the @code{GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_MEAN} and
@code{GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_STD} bit flats respectively.
+@end deftypefun
@deftypefun {gal_data_t *} gal_statistics_outlier_bydistance (int
@code{pos1_neg0}, gal_data_t @code{*input}, size_t @code{window_size}, float
@code{sigma}, float @code{sigclip_multip}, float @code{sigclip_param}, int
@code{inplace}, int @code{quiet})
@@ -45052,7 +45883,7 @@ Let's get back to PGPLOT for the sake of WCSLIB.
Installing it is a little tricky (mainly because it is so old!).
You can download the most recent version from the FTP link in its web
page@footnote{@url{http://www.astro.caltech.edu/~tjp/pgplot/}}.
-You can unpack it with the @command{tar xf} command.
+You can unpack it with the @command{tar -xf} command.
Let's assume the directory you have unpacked it to is @file{PGPLOT}, most
probably it is: @file{/home/username/Downloads/pgplot/}.
Open the @file{drivers.list} file:
@example
diff --git a/lib/arithmetic.c b/lib/arithmetic.c
index 40922aaf..2e6c9fcf 100644
--- a/lib/arithmetic.c
+++ b/lib/arithmetic.c
@@ -45,6 +45,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gnuastro/blank.h>
#include <gnuastro/units.h>
#include <gnuastro/qsort.h>
+#include <gnuastro/binary.h>
#include <gnuastro/pointer.h>
#include <gnuastro/threads.h>
#include <gnuastro/dimension.h>
@@ -1673,6 +1674,9 @@ struct multioperandparams
uint8_t *hasblank; /* Array of 0s or 1s for each input. */
float p1; /* Sigma-cliping parameter 1. */
float p2; /* Sigma-cliping parameter 2. */
+ uint8_t clipflags; /* Extra clipping measurements. */
+ gal_data_t *center; /* Center for dilated clipping. */
+ gal_data_t *spread; /* Spread for dilated clipping. */
};
@@ -1969,11 +1973,58 @@ struct multioperandparams
-#define MULTIOPERAND_SIGCLIP(TYPE) { \
+#define MULTIOPERAND_MAD(TYPE) { \
size_t n, j; \
- gal_data_t *sclip; \
+ gal_data_t *mad; \
+ float *o=p->out->array; \
+ TYPE *pixs=gal_pointer_allocate(p->list->type, p->dnum, 0, \
+ __func__, "pixs"); \
+ gal_data_t *cont=gal_data_alloc(pixs, p->list->type, 1, &p->dnum, \
+ NULL, 0, -1, 1, NULL, NULL, NULL); \
+ \
+ /* Go over all the pixels assigned to this thread. */ \
+ for(tind=0; tprm->indexs[tind] != GAL_BLANK_SIZE_T; ++tind) \
+ { \
+ /* Initialize, 'j' is desired pixel's index. */ \
+ n=0; \
+ j=tprm->indexs[tind]; \
+ \
+ /* Read the necessay values from each input. */ \
+ for(i=0;i<p->dnum;++i) pixs[n++]=a[i][j]; \
+ \
+ /* If there are any elements, do the measurement. */ \
+ if(n) \
+ { \
+ /* Calculate the quantile and put it in the output. */ \
+ mad=gal_statistics_mad(cont, 1); \
+ mad=gal_data_copy_to_new_type_free(mad, GAL_TYPE_FLOAT32); \
+ memcpy(&o[j], mad->array, gal_type_sizeof(GAL_TYPE_FLOAT32)); \
+ gal_data_free(mad); \
+ \
+ /* Since we are sorting in place, the size, and flags */ \
+ /* need to be reset. */ \
+ cont->flag=0; \
+ cont->size=cont->dsize[0]=p->dnum; \
+ } \
+ else \
+ o[j]=NAN; \
+ } \
+ \
+ /* Clean up (note that 'pixs' is inside of 'cont'). */ \
+ gal_data_free(cont); \
+ }
+
+
+
+
+
+#define MULTIOPERAND_CLIP(TYPE, SIG1_MAD0) { \
+ size_t n, j; \
+ gal_data_t *clip; \
uint32_t *N=p->out->array; \
- float *sarr, *o=p->out->array; \
+ float *carr, *o=p->out->array; \
+ float *center=p->center?p->center->array:NULL; \
+ float *spread=p->spread?p->spread->array:NULL; \
TYPE *pixs=gal_pointer_allocate(p->list->type, p->dnum, 0, \
__func__, "pixs"); \
gal_data_t *cont=gal_data_alloc(pixs, p->list->type, 1, &p->dnum, \
@@ -1992,21 +2043,59 @@ struct multioperandparams
/* If there are any usable elements, do the measurement. */ \
if(n) \
{ \
- /* Calculate the sigma-clip and write it in. */ \
- sclip=gal_statistics_sigma_clip(cont, p->p1, p->p2, 1, 1); \
- sarr=sclip->array; \
+ /* Do the clipping and write it in the output. */ \
+ clip = ( SIG1_MAD0 \
+ ? gal_statistics_clip_sigma(cont, p->p1, p->p2, \
+ p->clipflags, 1, 1) \
+ : gal_statistics_clip_mad(cont, p->p1, p->p2, \
+ p->clipflags, 1, 1) ); \
+ carr=clip->array; \
switch(p->operator) \
{ \
- case GAL_ARITHMETIC_OP_SIGCLIP_STD: o[j]=sarr[3]; break;\
- case GAL_ARITHMETIC_OP_SIGCLIP_MEAN: o[j]=sarr[2]; break;\
- case GAL_ARITHMETIC_OP_SIGCLIP_MEDIAN: o[j]=sarr[1]; break;\
- case GAL_ARITHMETIC_OP_SIGCLIP_NUMBER: N[j]=sarr[0]; break;\
+ /* The position of the final value is the same for any */ \
+ /* type of clipping. */ \
+ case GAL_ARITHMETIC_OP_SIGCLIP_MAD: \
+ case GAL_ARITHMETIC_OP_MADCLIP_MAD: \
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_MAD: \
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_MAD: \
+ o[j]=carr[GAL_STATISTICS_CLIP_OUTCOL_MAD]; break; \
+ case GAL_ARITHMETIC_OP_SIGCLIP_STD: \
+ case GAL_ARITHMETIC_OP_MADCLIP_STD: \
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_STD: \
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_STD: \
+ o[j]=carr[GAL_STATISTICS_CLIP_OUTCOL_STD]; break; \
+ case GAL_ARITHMETIC_OP_SIGCLIP_MEAN: \
+ case GAL_ARITHMETIC_OP_MADCLIP_MEAN: \
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_MEAN: \
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_MEAN: \
+ o[j]=carr[GAL_STATISTICS_CLIP_OUTCOL_MEAN]; break; \
+ case GAL_ARITHMETIC_OP_SIGCLIP_MEDIAN: \
+ case GAL_ARITHMETIC_OP_MADCLIP_MEDIAN: \
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_MEDIAN: \
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_MEDIAN: \
+ o[j]=carr[GAL_STATISTICS_CLIP_OUTCOL_MEDIAN]; break; \
+ case GAL_ARITHMETIC_OP_SIGCLIP_NUMBER: \
+ case GAL_ARITHMETIC_OP_MADCLIP_NUMBER: \
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_NUMBER: \
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_NUMBER: \
+ N[j]=carr[GAL_STATISTICS_CLIP_OUTCOL_NUMBER_USED]; break; \
default: \
error(EXIT_FAILURE, 0, "%s: a bug! the code %d is not " \
"valid for sigma-clipping results", __func__, \
p->operator); \
} \
- gal_data_free(sclip); \
+ \
+ /* If we are on clip-dilate operators, keep the values. */ \
+ if(center) \
+ { \
+ center[j]=carr[GAL_STATISTICS_CLIP_OUTCOL_MEDIAN]; \
+ spread[j]=( SIG1_MAD0 \
+ ? carr[GAL_STATISTICS_CLIP_OUTCOL_STD] \
+ : carr[GAL_STATISTICS_CLIP_OUTCOL_MAD] ); \
+ } \
+ \
+ /* Free the output of the clipping. */ \
+ gal_data_free(clip); \
\
/* Since we are doing sigma-clipping in place, the size, */ \
/* and flags need to be reset. */ \
@@ -2073,6 +2162,10 @@ struct multioperandparams
MULTIOPERAND_STD; \
break; \
\
+ case GAL_ARITHMETIC_OP_MAD: \
+ MULTIOPERAND_MAD(TYPE); \
+ break; \
+ \
case GAL_ARITHMETIC_OP_MEDIAN: \
MULTIOPERAND_MEDIAN(TYPE, QSORT_F); \
break; \
@@ -2081,11 +2174,30 @@ struct multioperandparams
MULTIOPERAND_QUANTILE(TYPE); \
break; \
\
+ case GAL_ARITHMETIC_OP_SIGCLIP_MAD: \
case GAL_ARITHMETIC_OP_SIGCLIP_STD: \
case GAL_ARITHMETIC_OP_SIGCLIP_MEAN: \
case GAL_ARITHMETIC_OP_SIGCLIP_MEDIAN: \
case GAL_ARITHMETIC_OP_SIGCLIP_NUMBER: \
- MULTIOPERAND_SIGCLIP(TYPE); \
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_MAD: \
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_STD: \
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_MEAN: \
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_MEDIAN: \
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_NUMBER: \
+ MULTIOPERAND_CLIP(TYPE, 1); \
+ break; \
+ \
+ case GAL_ARITHMETIC_OP_MADCLIP_MAD: \
+ case GAL_ARITHMETIC_OP_MADCLIP_STD: \
+ case GAL_ARITHMETIC_OP_MADCLIP_MEAN: \
+ case GAL_ARITHMETIC_OP_MADCLIP_MEDIAN: \
+ case GAL_ARITHMETIC_OP_MADCLIP_NUMBER: \
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_MAD: \
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_STD: \
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_MEAN: \
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_MEDIAN: \
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_NUMBER: \
+ MULTIOPERAND_CLIP(TYPE, 0); \
break; \
\
default: \
@@ -2156,24 +2268,18 @@ multioperand_on_thread(void *in_prm)
-/* The single operator in this function is assumed to be a linked list. The
- number of operators is determined from the fact that the last node in
- the linked list must have a NULL pointer as its 'next' element. */
-static gal_data_t *
-arithmetic_multioperand(int operator, int flags, gal_data_t *list,
- gal_data_t *params, size_t numthreads)
+static void
+arithmetic_multioperand_prepare(struct multioperandparams *p, int flags,
+ gal_data_t *params)
{
- size_t i=0, dnum=1;
- float p1=NAN, p2=NAN;
- struct multioperandparams p;
- gal_data_t *out, *tmp, *ttmp;
- uint8_t *hasblank, otype=GAL_TYPE_INVALID;
-
-
- /* For generality, 'list' can be a NULL pointer, in that case, this
- function will return a NULL pointer and avoid further processing. */
- if(list==NULL) return NULL;
+ size_t i;
+ gal_data_t *tmp;
+ uint8_t otype=GAL_TYPE_INVALID;
+ /* Initializations. */
+ p->dnum=1;
+ p->clipflags=0;
+ p->p1=p->p2=NAN;
/* If any parameters are given, prepare them. */
for(tmp=params; tmp!=NULL; tmp=tmp->next)
@@ -2187,17 +2293,17 @@ arithmetic_multioperand(int operator, int flags,
gal_data_t *list,
__func__);
/* Write them */
- if(isnan(p1)) p1=((float *)(tmp->array))[0];
- else p2=((float *)(tmp->array))[0];
+ if(isnan(p->p1)) p->p1=((float *)(tmp->array))[0];
+ else p->p2=((float *)(tmp->array))[0];
/* Operator specific, parameter sanity checks. */
- switch(operator)
+ switch(p->operator)
{
case GAL_ARITHMETIC_OP_QUANTILE:
- if(p1<0 || p1>1)
+ if(p->p1<0 || p->p1>1)
error(EXIT_FAILURE, 0, "%s: the parameter given to the "
"'quantile' operator must be between (and including) "
- "0 and 1. The given value is: %g", __func__, p1);
+ "0 and 1. The given value is: %g", __func__, p->p1);
break;
}
}
@@ -2205,16 +2311,16 @@ arithmetic_multioperand(int operator, int flags,
gal_data_t *list,
/* Do a simple sanity check, comparing the operand on top of the list to
the rest of the operands within the list. */
- for(tmp=list->next;tmp!=NULL;tmp=tmp->next)
+ for(tmp=p->list->next;tmp!=NULL;tmp=tmp->next)
{
/* Increment the number of structures. */
- ++dnum;
+ ++p->dnum;
/* Check the types. */
- if(tmp->type!=list->type)
- error(EXIT_FAILURE, 0, "%s: the types of all operands to the '%s' "
- "operator must be same", __func__,
- gal_arithmetic_operator_string(operator));
+ if(tmp->type!=p->list->type)
+ error(EXIT_FAILURE, 0, "%s: the types of all operands to the "
+ "'%s' operator must be same", __func__,
+ gal_arithmetic_operator_string(p->operator));
/* Make sure they actually have any data. */
if(tmp->size==0 || tmp->array==NULL)
@@ -2222,63 +2328,211 @@ arithmetic_multioperand(int operator, int flags,
gal_data_t *list,
"have any data", __func__);
/* Check the sizes. */
- if( gal_dimension_is_different(list, tmp) )
- error(EXIT_FAILURE, 0, "%s: the sizes of all operands to the '%s' "
- "operator must be same", __func__,
- gal_arithmetic_operator_string(operator));
+ if( gal_dimension_is_different(p->list, tmp) )
+ error(EXIT_FAILURE, 0, "%s: the sizes of all operands to the "
+ "'%s' operator must be same", __func__,
+ gal_arithmetic_operator_string(p->operator));
}
- /* Set the output dataset type. */
- switch(operator)
+ /* Set the type of the output and any other operator-specific setting. */
+ switch(p->operator)
{
- case GAL_ARITHMETIC_OP_MIN: otype=list->type; break;
- case GAL_ARITHMETIC_OP_MAX: otype=list->type; break;
+ case GAL_ARITHMETIC_OP_MIN: otype=p->list->type; break;
+ case GAL_ARITHMETIC_OP_MAX: otype=p->list->type; break;
case GAL_ARITHMETIC_OP_NUMBER: otype=GAL_TYPE_UINT32; break;
case GAL_ARITHMETIC_OP_SUM: otype=GAL_TYPE_FLOAT32; break;
case GAL_ARITHMETIC_OP_MEAN: otype=GAL_TYPE_FLOAT32; break;
case GAL_ARITHMETIC_OP_STD: otype=GAL_TYPE_FLOAT32; break;
- case GAL_ARITHMETIC_OP_MEDIAN: otype=GAL_TYPE_FLOAT32; break;
- case GAL_ARITHMETIC_OP_QUANTILE: otype=list->type; break;
+ case GAL_ARITHMETIC_OP_MAD: otype=GAL_TYPE_FLOAT32; break;
+ case GAL_ARITHMETIC_OP_QUANTILE: otype=p->list->type; break;
case GAL_ARITHMETIC_OP_SIGCLIP_STD: otype=GAL_TYPE_FLOAT32; break;
- case GAL_ARITHMETIC_OP_SIGCLIP_MEAN: otype=GAL_TYPE_FLOAT32; break;
- case GAL_ARITHMETIC_OP_SIGCLIP_MEDIAN: otype=GAL_TYPE_FLOAT32; break;
case GAL_ARITHMETIC_OP_SIGCLIP_NUMBER: otype=GAL_TYPE_UINT32; break;
+ case GAL_ARITHMETIC_OP_MADCLIP_NUMBER: otype=GAL_TYPE_UINT32; break;
+ case GAL_ARITHMETIC_OP_MEDIAN:
+ case GAL_ARITHMETIC_OP_SIGCLIP_MEDIAN:
+ case GAL_ARITHMETIC_OP_MADCLIP_MEDIAN:
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_MEDIAN:
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_MEDIAN:
+ otype=GAL_TYPE_FLOAT32; break;
+ case GAL_ARITHMETIC_OP_MADCLIP_MAD:
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_MAD:
+ otype=GAL_TYPE_FLOAT32; break;
+ case GAL_ARITHMETIC_OP_SIGCLIP_MEAN:
+ case GAL_ARITHMETIC_OP_MADCLIP_MEAN:
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_MEAN:
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_MEAN:
+ p->clipflags=GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_MEAN;
+ otype=GAL_TYPE_FLOAT32; break;
+ case GAL_ARITHMETIC_OP_MADCLIP_STD:
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_STD:
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_STD:
+ p->clipflags=GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_STD;
+ otype=GAL_TYPE_FLOAT32; break;
+ case GAL_ARITHMETIC_OP_SIGCLIP_MAD:
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_MAD:
+ p->clipflags=GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_MAD;
+ otype=GAL_TYPE_FLOAT32; break;
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_NUMBER:
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_NUMBER:
+ otype=GAL_TYPE_UINT32; break;
default:
error(EXIT_FAILURE, 0, "%s: operator code %d isn't recognized",
- __func__, operator);
+ __func__, p->operator);
}
-
/* Set the output data structure. */
- if( (flags & GAL_ARITHMETIC_FLAG_INPLACE) && otype==list->type)
- out = list; /* The top element in the list. */
+ if( (flags & GAL_ARITHMETIC_FLAG_INPLACE) && otype==p->list->type)
+ p->out = p->list; /* The top element in the list. */
else
- out = gal_data_alloc(NULL, otype, list->ndim, list->dsize,
- list->wcs, 0, list->minmapsize, list->quietmmap,
- NULL, NULL, NULL);
+ p->out = gal_data_alloc(NULL, otype, p->list->ndim, p->list->dsize,
+ p->list->wcs, 0, p->list->minmapsize,
+ p->list->quietmmap, NULL, NULL, NULL);
+ /* The filled clipping operators need to extra allocations. */
+ switch(p->operator)
+ {
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_MAD:
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_MAD:
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_STD:
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_STD:
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_MEAN:
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_MEAN:
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_MEDIAN:
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_MEDIAN:
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_NUMBER:
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_NUMBER:
+ p->center = gal_data_alloc(NULL, GAL_TYPE_FLOAT32, p->list->ndim,
+ p->list->dsize, p->list->wcs, 0,
+ p->list->minmapsize, p->list->quietmmap,
+ NULL, NULL, NULL);
+ p->spread = gal_data_alloc(NULL, GAL_TYPE_FLOAT32, p->list->ndim,
+ p->list->dsize, p->list->wcs, 0,
+ p->list->minmapsize, p->list->quietmmap,
+ NULL, NULL, NULL);
+ break;
+ default: p->center=p->spread=NULL;
+ }
/* hasblank is used to see if a blank value should be checked for each
list element or not. */
- hasblank=gal_pointer_allocate(GAL_TYPE_UINT8, dnum, 0, __func__,
- "hasblank");
- for(tmp=list;tmp!=NULL;tmp=tmp->next)
- hasblank[i++]=gal_blank_present(tmp, 0);
+ i=0;
+ p->hasblank=gal_pointer_allocate(GAL_TYPE_UINT8, p->dnum, 0, __func__,
+ "hasblank");
+ for(tmp=p->list;tmp!=NULL;tmp=tmp->next)
+ p->hasblank[i++]=gal_blank_present(tmp, 0);
+}
+
+
+
+
+
+static void
+arithmetic_multioperand_clip_fill(struct multioperandparams *p,
+ int operator, gal_data_t *list,
+ size_t numthreads)
+{
+ size_t one=1;
+ int aflags=GAL_ARITHMETIC_FLAG_NUMOK; /* Don't free the inputs. */
+ gal_data_t *tmp, *input, *upper, *lower, *multip;
+
+ /* Find the upper and lower thresholds based on the user's desired
+ multiple. */
+ multip=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, 1, &one, NULL, 0, -1, 1,
+ NULL, NULL, NULL);
+ ((float *)(multip->array))[0]=p->p1;
+ tmp=gal_arithmetic(GAL_ARITHMETIC_OP_MULTIPLY, 1, aflags,
+ p->spread, multip);
+ upper=gal_arithmetic(GAL_ARITHMETIC_OP_PLUS, 1, aflags,
+ p->center, tmp);
+ lower=gal_arithmetic(GAL_ARITHMETIC_OP_MINUS, 1, aflags,
+ p->center, tmp);
+ gal_data_free(tmp);
+
+ /* For a check.
+ gal_fits_img_write(lower, "test.fits", NULL, NULL);
+ gal_fits_img_write(upper, "test.fits", NULL, NULL);
+ printf("%s: GOOD\n", __func__); exit(0);
+ */
+
+ /* Dilate all the sigma-clipped pixels. */
+ for(input=list;input!=NULL;input=input->next)
+ {
+ /* Flag all the pixels that are higher/lower than the limits. */
+ tmp=gal_arithmetic(GAL_ARITHMETIC_OP_OR, 1, aflags,
+ gal_arithmetic(GAL_ARITHMETIC_OP_GT, 1, aflags,
+ input, upper),
+ gal_arithmetic(GAL_ARITHMETIC_OP_LE, 1, aflags,
+ input, lower));
+
+ /* For a 1D array, start with erosion because a two single elements
+ outside the range can mask a very large portion of the input
+ (after "filling" holes). */
+ tmp = ( list->ndim==1
+ ? gal_binary_erode(tmp, 1, 1, 1)
+ : gal_binary_dilate(tmp, 1, 1, 1) );
+ gal_binary_holes_fill(tmp, list->ndim, -1);
+ tmp=gal_binary_erode(tmp, 2, list->ndim, 1);
+ tmp=gal_binary_dilate(tmp, 2, list->ndim, 1);
+
+ /* Set all the 1-vaued pixels in the binary image to NaN in the
+ input. */
+ gal_blank_flag_apply(input, tmp);
+ gal_data_free(tmp);
+
+ /* For a check.
+ gal_fits_img_write(input, "test.fits", NULL, NULL); */
+ }
+ gal_data_free(upper); /* We are freeing these here to avoid*/
+ gal_data_free(lower); /* keeping extra RAM. */
+
+ /* Re-run sigma-clipping on threads. */
+ gal_threads_spin_off(multioperand_on_thread, p, p->out->size,
+ numthreads, list->minmapsize, list->quietmmap);
+}
+
+
+
+/* The single operator in this function is assumed to be a linked list. The
+ number of operators is determined from the fact that the last node in
+ the linked list must have a NULL pointer as its 'next' element. */
+static gal_data_t *
+arithmetic_multioperand(int operator, int flags, gal_data_t *list,
+ gal_data_t *params, size_t numthreads)
+{
+ gal_data_t *tmp, *ttmp;
+ struct multioperandparams p;
- /* Set the parameters necessary for multithreaded operation and spin them
- off to do apply the operator. */
- p.p1=p1;
- p.p2=p2;
- p.out=out;
+ /* For generality, 'list' can be a NULL pointer, in that case, this
+ function will return a NULL pointer and avoid further processing. */
+ if(list==NULL) return NULL;
+
+ /* Sanity checks and preparations. */
p.list=list;
- p.dnum=dnum;
p.operator=operator;
- p.hasblank=hasblank;
- gal_threads_spin_off(multioperand_on_thread, &p, out->size, numthreads,
- list->minmapsize, list->quietmmap);
+ arithmetic_multioperand_prepare(&p, flags, params);
+ /* Spin off the threads apply the operator. */
+ gal_threads_spin_off(multioperand_on_thread, &p, p.out->size,
+ numthreads, list->minmapsize, list->quietmmap);
+
+ /* Do the filled clipping if necessary. */
+ switch(operator)
+ {
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_MAD:
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_MAD:
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_STD:
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_STD:
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_MEAN:
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_MEAN:
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_MEDIAN:
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_MEDIAN:
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_NUMBER:
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_NUMBER:
+ arithmetic_multioperand_clip_fill(&p, operator, list, numthreads);
+ break;
+ }
/* Clean up and return. Note that the operation might have been done in
place. In that case, a list element was used. So we need to check
@@ -2291,13 +2545,13 @@ arithmetic_multioperand(int operator, int flags,
gal_data_t *list,
while(tmp!=NULL)
{
ttmp=tmp->next;
- if(tmp==out) tmp->next=NULL; else gal_data_free(tmp);
+ if(tmp==p.out) tmp->next=NULL; else gal_data_free(tmp);
tmp=ttmp;
}
if(params) gal_list_data_free(params);
}
- free(hasblank);
- return out;
+ free(p.hasblank);
+ return p.out;
}
@@ -2444,8 +2698,15 @@ arithmetic_binary(int operator, int flags, gal_data_t
*l, gal_data_t *r)
if( l->size==0 || l->array==NULL || r->size==0 || r->array==NULL )
{
if(l->array==0 || l->array==NULL)
- { if(flags & GAL_ARITHMETIC_FLAG_FREE) gal_data_free(r); return l;}
- else {if(flags & GAL_ARITHMETIC_FLAG_FREE) gal_data_free(l); return r;}
+ {
+ if(flags & GAL_ARITHMETIC_FLAG_FREE) gal_data_free(r);
+ return l;
+ }
+ else
+ {
+ if(flags & GAL_ARITHMETIC_FLAG_FREE) gal_data_free(l);
+ return r;
+ }
}
@@ -3582,6 +3843,8 @@ gal_arithmetic_set_operator(char *string, size_t
*num_operands)
{ op=GAL_ARITHMETIC_OP_MEAN; *num_operands=-1; }
else if (!strcmp(string, "std"))
{ op=GAL_ARITHMETIC_OP_STD; *num_operands=-1; }
+ else if (!strcmp(string, "mad"))
+ { op=GAL_ARITHMETIC_OP_MAD; *num_operands=-1; }
else if (!strcmp(string, "median"))
{ op=GAL_ARITHMETIC_OP_MEDIAN; *num_operands=-1; }
else if (!strcmp(string, "quantile"))
@@ -3592,8 +3855,40 @@ gal_arithmetic_set_operator(char *string, size_t
*num_operands)
{ op=GAL_ARITHMETIC_OP_SIGCLIP_MEAN; *num_operands=-1; }
else if (!strcmp(string, "sigclip-median"))
{ op=GAL_ARITHMETIC_OP_SIGCLIP_MEDIAN; *num_operands=-1; }
+ else if (!strcmp(string, "sigclip-mad"))
+ { op=GAL_ARITHMETIC_OP_SIGCLIP_MAD; *num_operands=-1; }
else if (!strcmp(string, "sigclip-std"))
{ op=GAL_ARITHMETIC_OP_SIGCLIP_STD; *num_operands=-1; }
+ else if (!strcmp(string, "madclip-number"))
+ { op=GAL_ARITHMETIC_OP_MADCLIP_NUMBER; *num_operands=-1; }
+ else if (!strcmp(string, "madclip-mean"))
+ { op=GAL_ARITHMETIC_OP_MADCLIP_MEAN; *num_operands=-1; }
+ else if (!strcmp(string, "madclip-median"))
+ { op=GAL_ARITHMETIC_OP_MADCLIP_MEDIAN; *num_operands=-1; }
+ else if (!strcmp(string, "madclip-mad"))
+ { op=GAL_ARITHMETIC_OP_MADCLIP_MAD; *num_operands=-1; }
+ else if (!strcmp(string, "madclip-std"))
+ { op=GAL_ARITHMETIC_OP_MADCLIP_STD; *num_operands=-1; }
+ else if (!strcmp(string, "madclip-fill-number"))
+ { op=GAL_ARITHMETIC_OP_MADCLIP_FILL_NUMBER; *num_operands=-1; }
+ else if (!strcmp(string, "madclip-fill-mean"))
+ { op=GAL_ARITHMETIC_OP_MADCLIP_FILL_MEAN; *num_operands=-1; }
+ else if (!strcmp(string, "madclip-fill-median"))
+ { op=GAL_ARITHMETIC_OP_MADCLIP_FILL_MEDIAN; *num_operands=-1; }
+ else if (!strcmp(string, "madclip-fill-mad"))
+ { op=GAL_ARITHMETIC_OP_MADCLIP_FILL_MAD; *num_operands=-1; }
+ else if (!strcmp(string, "madclip-fill-std"))
+ { op=GAL_ARITHMETIC_OP_MADCLIP_FILL_STD; *num_operands=-1; }
+ else if (!strcmp(string, "sigclip-fill-number"))
+ { op=GAL_ARITHMETIC_OP_SIGCLIP_FILL_NUMBER; *num_operands=-1; }
+ else if (!strcmp(string, "sigclip-fill-mean"))
+ { op=GAL_ARITHMETIC_OP_SIGCLIP_FILL_MEAN; *num_operands=-1; }
+ else if (!strcmp(string, "sigclip-fill-median"))
+ { op=GAL_ARITHMETIC_OP_SIGCLIP_FILL_MEDIAN; *num_operands=-1; }
+ else if (!strcmp(string, "sigclip-fill-mad"))
+ { op=GAL_ARITHMETIC_OP_SIGCLIP_FILL_MAD; *num_operands=-1; }
+ else if (!strcmp(string, "sigclip-fill-std"))
+ { op=GAL_ARITHMETIC_OP_SIGCLIP_FILL_STD; *num_operands=-1; }
/* To one-dimension (only based on values). */
else if (!strcmp(string, "unique"))
@@ -3842,19 +4137,48 @@ gal_arithmetic_operator_string(int operator)
case GAL_ARITHMETIC_OP_SUM: return "sum";
case GAL_ARITHMETIC_OP_MEAN: return "mean";
case GAL_ARITHMETIC_OP_STD: return "std";
+ case GAL_ARITHMETIC_OP_MAD: return "mad";
case GAL_ARITHMETIC_OP_MEDIAN: return "median";
case GAL_ARITHMETIC_OP_QUANTILE: return "quantile";
case GAL_ARITHMETIC_OP_SIGCLIP_NUMBER: return "sigclip-number";
case GAL_ARITHMETIC_OP_SIGCLIP_MEDIAN: return "sigclip-median";
case GAL_ARITHMETIC_OP_SIGCLIP_MEAN: return "sigclip-mean";
- case GAL_ARITHMETIC_OP_SIGCLIP_STD: return "sigclip-number";
+ case GAL_ARITHMETIC_OP_SIGCLIP_MAD: return "sigclip-mad";
+ case GAL_ARITHMETIC_OP_SIGCLIP_STD: return "sigclip-std";
+ case GAL_ARITHMETIC_OP_MADCLIP_NUMBER: return "madclip-number";
+ case GAL_ARITHMETIC_OP_MADCLIP_MEDIAN: return "madclip-median";
+ case GAL_ARITHMETIC_OP_MADCLIP_MEAN: return "madclip-mean";
+ case GAL_ARITHMETIC_OP_MADCLIP_MAD: return "madclip-mad";
+ case GAL_ARITHMETIC_OP_MADCLIP_STD: return "madclip-std";
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_NUMBER:
+ return "madclip-dilate-number";
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_MEDIAN:
+ return "madclip-dilate-median";
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_MEAN:
+ return "madclip-dilate-mean";
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_MAD:
+ return "madclip-dilate-mad";
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_STD:
+ return "madclip-dilate-std";
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_NUMBER:
+ return "sigclip-dilate-number";
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_MEDIAN:
+ return "sigclip-dilate-median";
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_MEAN:
+ return "sigclip-dilate-mean";
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_MAD:
+ return "sigclip-dilate-mad";
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_STD:
+ return "sigclip-dilate-std";
case GAL_ARITHMETIC_OP_MKNOISE_SIGMA: return "mknoise-sigma";
- case GAL_ARITHMETIC_OP_MKNOISE_SIGMA_FROM_MEAN: return
"mknoise-sigma-from-mean";
+ case GAL_ARITHMETIC_OP_MKNOISE_SIGMA_FROM_MEAN:
+ return "mknoise-sigma-from-mean";
case GAL_ARITHMETIC_OP_MKNOISE_POISSON: return "mknoise-poisson";
case GAL_ARITHMETIC_OP_MKNOISE_UNIFORM: return "mknoise-uniform";
case GAL_ARITHMETIC_OP_RANDOM_FROM_HIST:return "random-from-hist";
- case GAL_ARITHMETIC_OP_RANDOM_FROM_HIST_RAW:return "random-from-hist-raw";
+ case GAL_ARITHMETIC_OP_RANDOM_FROM_HIST_RAW:
+ return "random-from-hist-raw";
case GAL_ARITHMETIC_OP_STITCH: return "stitch";
@@ -4142,16 +4466,33 @@ gal_arithmetic(int operator, size_t numthreads, int
flags, ...)
/* Multi-operand operators */
case GAL_ARITHMETIC_OP_MIN:
case GAL_ARITHMETIC_OP_MAX:
- case GAL_ARITHMETIC_OP_NUMBER:
case GAL_ARITHMETIC_OP_SUM:
- case GAL_ARITHMETIC_OP_MEAN:
case GAL_ARITHMETIC_OP_STD:
+ case GAL_ARITHMETIC_OP_MAD:
+ case GAL_ARITHMETIC_OP_MEAN:
+ case GAL_ARITHMETIC_OP_NUMBER:
case GAL_ARITHMETIC_OP_MEDIAN:
case GAL_ARITHMETIC_OP_QUANTILE:
+ case GAL_ARITHMETIC_OP_SIGCLIP_MAD:
+ case GAL_ARITHMETIC_OP_MADCLIP_MAD:
case GAL_ARITHMETIC_OP_SIGCLIP_STD:
+ case GAL_ARITHMETIC_OP_MADCLIP_STD:
case GAL_ARITHMETIC_OP_SIGCLIP_MEAN:
+ case GAL_ARITHMETIC_OP_MADCLIP_MEAN:
case GAL_ARITHMETIC_OP_SIGCLIP_MEDIAN:
+ case GAL_ARITHMETIC_OP_MADCLIP_MEDIAN:
case GAL_ARITHMETIC_OP_SIGCLIP_NUMBER:
+ case GAL_ARITHMETIC_OP_MADCLIP_NUMBER:
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_MAD:
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_MAD:
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_STD:
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_STD:
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_MEAN:
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_MEAN:
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_MEDIAN:
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_MEDIAN:
+ case GAL_ARITHMETIC_OP_MADCLIP_FILL_NUMBER:
+ case GAL_ARITHMETIC_OP_SIGCLIP_FILL_NUMBER:
d1 = va_arg(va, gal_data_t *);
d2 = va_arg(va, gal_data_t *);
out=arithmetic_multioperand(operator, flags, d1, d2, numthreads);
diff --git a/lib/binary.c b/lib/binary.c
index 780de4d4..7277d7ce 100644
--- a/lib/binary.c
+++ b/lib/binary.c
@@ -537,8 +537,8 @@ gal_binary_connected_components(gal_data_t *binary,
gal_data_t **out,
array, then give them the blank labeled array. Note that since
their value will not be 0, they will also not be labeled. */
l=lab->array;
- bf=(b=binary->array)+binary->size; /* Library must have no side effect, */
- if( gal_blank_present(binary, 0) ) /* So blank flag should not be changed.*/
+ bf=(b=binary->array)+binary->size; /* Library must have no side effect:*/
+ if( gal_blank_present(binary, 0) ) /* blank flag should not be changed.*/
do *l++ = *b==GAL_BLANK_UINT8 ? GAL_BLANK_INT32 : 0; while(++b<bf);
@@ -548,6 +548,7 @@ gal_binary_connected_components(gal_data_t *binary,
gal_data_t **out,
l=lab->array;
b=binary->array;
for(i=0;i<binary->size;++i)
+
/* Check if this pixel is already labeled. */
if( b[i] && l[i]==0 )
{
@@ -650,7 +651,8 @@ gal_binary_connected_indexs(gal_data_t *binary, int
connectivity)
/* Parsing has finished, put all the indexs into an array. */
onelabarr=gal_list_sizet_to_array(onelab, 1, &onelabnum);
gal_list_data_add_alloc(&lines, onelabarr, GAL_TYPE_SIZE_T, 1,
- &onelabnum, NULL, 0, -1, 1, NULL, NULL, NULL);
+ &onelabnum, NULL, 0, -1, 1, NULL, NULL,
+ NULL);
/* Clean up. */
gal_list_sizet_free(onelab);
@@ -1051,18 +1053,16 @@ gal_binary_holes_fill(gal_data_t *input, int
connectivity, size_t maxsize)
"It has to be between 1 and the number of input's dimensions "
"(%zu)", __func__, connectivity, input->ndim);
-
/* Make the inverse image. */
inv=binary_make_padded_inverse(input, &tile);
-
/* Label the holes */
numholes=gal_binary_connected_components(inv, &holelabs, connectivity);
-
- /* Any pixel with a label larger than 1 is a hole in the input image and
- we should invert the respective pixel. To do it, we'll use the tile
- that was defined before, just change its block and array.*/
+ /* Any pixel with a label that is not touching the edges is a hole in the
+ input image and we should invert the respective pixel. To do it, we'll
+ use the tile that was defined before, just change its block and
+ array. */
in=input->array;
tile->array=gal_tile_block_relative_to_other(tile, holelabs);
tile->block=holelabs; /* has to be after correcting 'tile->array'. */
@@ -1090,9 +1090,16 @@ gal_binary_holes_fill(gal_data_t *input, int
connectivity, size_t maxsize)
/* The type of the tile is already known (it is 'int32_t') and we have no
output, so we'll just put 'int' as a place-holder. In this way we can
avoid the switch statement of GAL_TILE_PARSE_OPERATE, and directly use
- the workhorse macro 'GAL_TILE_PO_OISET'. */
+ the workhorse macro 'GAL_TILE_PO_OISET'.
+
+ For inputs of any dimensionality, the label=1 is not a hole. For 1D
+ data, the last label is also not a hole (because in 1D the start and
+ end of the array don't touch).
+ */
GAL_TILE_PO_OISET(int32_t, int, tile, NULL, 0, 0, {
- *in = *i>1 && *i!=GAL_BLANK_INT32 ? 1 : *in;
+ *in = ( *i==GAL_BLANK_INT32
+ || *i==1
+ || (input->ndim==1 && *i==numholes) ) ? *in : 1;
++in;
});
diff --git a/lib/dimension.c b/lib/dimension.c
index e3a347c6..e9d111d2 100644
--- a/lib/dimension.c
+++ b/lib/dimension.c
@@ -30,10 +30,13 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <stdlib.h>
#include <gnuastro/wcs.h>
+#include <gnuastro/binary.h>
#include <gnuastro/pointer.h>
#include <gnuastro/threads.h>
+#include <gnuastro/convolve.h>
#include <gnuastro/dimension.h>
#include <gnuastro/statistics.h>
+#include <gnuastro/arithmetic.h>
@@ -321,22 +324,22 @@ gal_dimension_dist_elliptical(double *center, double
*pa_deg, double *q,
s3 = sin( pa_deg[2] * DEGREESTORADIANS );
/* Calculate the distance. */
- Xr = x*( c3*c1 - s3*c2*s1 ) + y*( c3*s1 + s3*c2*c1) + z*( s3*s2 );
- Yr = x*( -1*s3*c1 - c3*c2*s1 ) + y*(-1*s3*s1 + c3*c2*c1) + z*( c3*s2 );
- Zr = x*( s1*s2 ) + y*(-1*s2*c1 ) + z*( c2 );
+ Xr = x*( c3*c1 - s3*c2*s1) + y*( c3*s1 + s3*c2*c1) + z*( s3*s2 );
+ Yr = x*(-1*s3*c1 - c3*c2*s1) + y*(-1*s3*s1 + c3*c2*c1) + z*( c3*s2 );
+ Zr = x*( s1*s2 ) + y*(-1*s2*c1 ) + z*( c2 );
return sqrt( Xr*Xr + Yr*Yr/q1/q1 + Zr*Zr/q2/q2 );
break;
default:
error(EXIT_FAILURE, 0, "%s: currently only 2 and 3 dimensional "
- "distances are supported, you have asked for an %zu-dimensional "
- "dataset", __func__, ndim);
+ "distances are supported, you have asked for an "
+ "%zu-dimensional dataset", __func__, ndim);
}
/* Control should never reach here. */
- error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to address the "
- "problem. Control should not reach the end of this function",
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to address "
+ "the problem. Control should not reach the end of this function",
__func__, PACKAGE_BUGREPORT);
return NAN;
}
@@ -372,10 +375,26 @@ enum dimension_collapse_operation
DIMENSION_COLLAPSE_MEAN,
DIMENSION_COLLAPSE_MEDIAN,
DIMENSION_COLLAPSE_NUMBER,
+ DIMENSION_COLLAPSE_MADCLIP_MAD,
+ DIMENSION_COLLAPSE_SIGCLIP_MAD,
+ DIMENSION_COLLAPSE_MADCLIP_STD,
DIMENSION_COLLAPSE_SIGCLIP_STD,
+ DIMENSION_COLLAPSE_MADCLIP_MEAN,
DIMENSION_COLLAPSE_SIGCLIP_MEAN,
+ DIMENSION_COLLAPSE_MADCLIP_MEDIAN,
DIMENSION_COLLAPSE_SIGCLIP_MEDIAN,
+ DIMENSION_COLLAPSE_MADCLIP_NUMBER,
DIMENSION_COLLAPSE_SIGCLIP_NUMBER,
+ DIMENSION_COLLAPSE_MADCLIP_FILL_MAD,
+ DIMENSION_COLLAPSE_SIGCLIP_FILL_MAD,
+ DIMENSION_COLLAPSE_MADCLIP_FILL_STD,
+ DIMENSION_COLLAPSE_SIGCLIP_FILL_STD,
+ DIMENSION_COLLAPSE_MADCLIP_FILL_MEAN,
+ DIMENSION_COLLAPSE_SIGCLIP_FILL_MEAN,
+ DIMENSION_COLLAPSE_MADCLIP_FILL_MEDIAN,
+ DIMENSION_COLLAPSE_SIGCLIP_FILL_MEDIAN,
+ DIMENSION_COLLAPSE_MADCLIP_FILL_NUMBER,
+ DIMENSION_COLLAPSE_SIGCLIP_FILL_NUMBER,
};
@@ -392,8 +411,9 @@ dimension_collapse_sanity_check(gal_data_t *in, gal_data_t
*weight,
/* The requested dimension to collapse cannot be larger than the input's
number of dimensions. */
if( c_dim > (in->ndim-1) )
- error(EXIT_FAILURE, 0, "%s: the input has %zu dimension(s), but you have "
- "asked to collapse dimension %zu", __func__, in->ndim, c_dim);
+ error(EXIT_FAILURE, 0, "%s: the input has %zu dimension(s), but you "
+ "have asked to collapse dimension %zu", __func__, in->ndim,
+ c_dim);
/* If there is no blank value, there is no point in calculating the
number of points in each collapsed dataset (when necessary). In that
@@ -405,8 +425,9 @@ dimension_collapse_sanity_check(gal_data_t *in, gal_data_t
*weight,
if(weight)
{
if( weight->ndim!=1 )
- error(EXIT_FAILURE, 0, "%s: the weight dataset has %zu dimensions, "
- "it must be one-dimensional", __func__, weight->ndim);
+ error(EXIT_FAILURE, 0, "%s: the weight dataset has %zu "
+ "dimensions, it must be one-dimensional", __func__,
+ weight->ndim);
if( in->dsize[c_dim]!=weight->size )
error(EXIT_FAILURE, 0, "%s: the weight dataset has %zu elements, "
"but the input dataset has %zu elements in dimension %zu",
@@ -825,7 +846,8 @@ gal_dimension_collapse_minmax(gal_data_t *in, size_t c_dim,
int max1_min0)
if(hasblank)
{
num=gal_data_alloc(NULL, GAL_TYPE_UINT8, outndim, outdsize, in->wcs,
- 1, in->minmapsize, in->quietmmap, NULL, NULL, NULL);
+ 1, in->minmapsize, in->quietmmap, NULL, NULL,
+ NULL);
iarr=num->array;
}
@@ -881,7 +903,8 @@ struct dimension_sortbased_p
static void
-dimension_csb_copy(gal_data_t *in, size_t from, gal_data_t *work, size_t to)
+dimension_csb_copy(gal_data_t *in, size_t from, gal_data_t *work,
+ size_t to)
{
memcpy(gal_pointer_increment(work->array, to, in->type),
gal_pointer_increment(in->array, from, in->type),
@@ -892,6 +915,235 @@ dimension_csb_copy(gal_data_t *in, size_t from,
gal_data_t *work, size_t to)
+static gal_data_t *
+dimension_collapse_sortbased_operation(struct dimension_sortbased_p *p,
+ gal_data_t *work, uint8_t clipflags,
+ uint8_t isfill)
+{
+ gal_data_t *out=NULL;
+
+ /* Do the operation. */
+ switch(p->operator)
+ {
+ case DIMENSION_COLLAPSE_MEDIAN:
+ out=gal_statistics_median(work, 1);
+ break;
+ case DIMENSION_COLLAPSE_MADCLIP_MAD:
+ case DIMENSION_COLLAPSE_MADCLIP_STD:
+ case DIMENSION_COLLAPSE_MADCLIP_MEAN:
+ case DIMENSION_COLLAPSE_MADCLIP_MEDIAN:
+ case DIMENSION_COLLAPSE_MADCLIP_NUMBER:
+ case DIMENSION_COLLAPSE_MADCLIP_FILL_MAD:
+ case DIMENSION_COLLAPSE_MADCLIP_FILL_STD:
+ case DIMENSION_COLLAPSE_MADCLIP_FILL_MEAN:
+ case DIMENSION_COLLAPSE_MADCLIP_FILL_MEDIAN:
+ case DIMENSION_COLLAPSE_MADCLIP_FILL_NUMBER:
+ /* When we are dealing with a clip-fill operator, the order of
+ the elements matter, so we don't want it to be done inplace.*/
+ out=gal_statistics_clip_mad(work, p->sclipmultip,
+ p->sclipparam, clipflags,
+ isfill?0:1, 1);
+ break;
+
+ case DIMENSION_COLLAPSE_SIGCLIP_MAD:
+ case DIMENSION_COLLAPSE_SIGCLIP_STD:
+ case DIMENSION_COLLAPSE_SIGCLIP_MEAN:
+ case DIMENSION_COLLAPSE_SIGCLIP_MEDIAN:
+ case DIMENSION_COLLAPSE_SIGCLIP_NUMBER:
+ case DIMENSION_COLLAPSE_SIGCLIP_FILL_MAD:
+ case DIMENSION_COLLAPSE_SIGCLIP_FILL_STD:
+ case DIMENSION_COLLAPSE_SIGCLIP_FILL_MEAN:
+ case DIMENSION_COLLAPSE_SIGCLIP_FILL_MEDIAN:
+ case DIMENSION_COLLAPSE_SIGCLIP_FILL_NUMBER:
+ /* When we are dealing with a clip-fill operator, the order of
+ the elements matter, so we don't want it to be done inplace.*/
+ out=gal_statistics_clip_sigma(work, p->sclipmultip,
+ p->sclipparam, clipflags,
+ isfill?0:1, 1);
+ break;
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at '%s' "
+ "to fix the problem. The operator code %d isn't "
+ "recognized", __func__, PACKAGE_BUGREPORT, p->operator);
+ }
+
+ /* Return the output. */
+ return out;
+}
+
+
+
+
+
+static gal_data_t *
+dimension_collapse_sortbased_fill(struct dimension_sortbased_p *p,
+ gal_data_t *fstat, gal_data_t *work,
+ gal_data_t *conv, uint8_t clipflags,
+ int check)
+{
+ size_t one=1;
+ int std1_mad0=0;
+ float *farr=fstat->array;
+ gal_data_t *formask=conv?conv:work, *out=NULL;
+ gal_data_t *tmp, *multip, *upper, *lower, *center, *spread;
+ int aflags=GAL_ARITHMETIC_FLAG_NUMOK; /* Don't free the inputs. */
+
+ /* Basic checks. */
+ switch(p->operator)
+ {
+ case DIMENSION_COLLAPSE_MADCLIP_FILL_MAD:
+ case DIMENSION_COLLAPSE_MADCLIP_FILL_STD:
+ case DIMENSION_COLLAPSE_MADCLIP_FILL_MEAN:
+ case DIMENSION_COLLAPSE_MADCLIP_FILL_MEDIAN:
+ case DIMENSION_COLLAPSE_MADCLIP_FILL_NUMBER:
+ std1_mad0=0; break;
+ case DIMENSION_COLLAPSE_SIGCLIP_FILL_MAD:
+ case DIMENSION_COLLAPSE_SIGCLIP_FILL_STD:
+ case DIMENSION_COLLAPSE_SIGCLIP_FILL_MEAN:
+ case DIMENSION_COLLAPSE_SIGCLIP_FILL_MEDIAN:
+ case DIMENSION_COLLAPSE_SIGCLIP_FILL_NUMBER:
+ std1_mad0=1; break;
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at '%s' to "
+ "fix the problem. The operator code '%d' is not recognized",
+ __func__, PACKAGE_BUGREPORT, p->operator);
+ }
+
+ /* Convert the first clipped statistics ('fstat') into two datasets with
+ the center and spread. */
+ center=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, 1, &one, NULL, 0, -1, 1,
+ NULL, NULL, NULL);
+ spread=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, 1, &one, NULL, 0, -1, 1,
+ NULL, NULL, NULL);
+ ((float *)(center->array))[0]=farr[GAL_STATISTICS_CLIP_OUTCOL_MEDIAN];
+ ((float *)(spread->array))[0]=farr[ std1_mad0
+ ? GAL_STATISTICS_CLIP_OUTCOL_STD
+ : GAL_STATISTICS_CLIP_OUTCOL_MAD ];
+
+ /* Find the upper and lower thresholds based on the user's desired
+ multiple. */
+ multip=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, 1, &one, NULL, 0, -1, 1,
+ NULL, NULL, NULL);
+ ((float *)(multip->array))[0]=p->sclipmultip;
+ tmp=gal_arithmetic(GAL_ARITHMETIC_OP_MULTIPLY, 1, aflags,
+ spread, multip);
+ upper=gal_arithmetic(GAL_ARITHMETIC_OP_PLUS, 1, aflags,
+ center, tmp);
+ lower=gal_arithmetic(GAL_ARITHMETIC_OP_MINUS, 1, aflags,
+ center, tmp);
+ gal_data_free(tmp);
+
+ /* Build a flag with bad elements having a value of 1. For a 1D dataset,
+ the largest possible "hole" size should be smaller than higher
+ dimensionality data, because two elements that are very far from each
+ other can create a very large "hole".*/
+ tmp=gal_arithmetic(GAL_ARITHMETIC_OP_OR, 1, aflags,
+ gal_arithmetic(GAL_ARITHMETIC_OP_GT, 1, aflags,
+ formask, upper),
+ gal_arithmetic(GAL_ARITHMETIC_OP_LE, 1, aflags,
+ formask, lower));
+
+ /* For a check (define 'int check' as an argument, and when calling this
+ function, set 'index==XXX').
+ if(check)
+ {
+ size_t i;
+ double *f=formask->array;
+ uint8_t *u=tmp->array;
+ for(i=0;i<formask->size;++i)
+ printf("%-5zu %-5u %f\n", i, u[i], f[i]);
+ printf("%s: upper:%f, lower: %f\n", __func__,
+ ((float *)(upper->array))[0],
+ ((float *)(lower->array))[0] );
+ }
+ */
+
+ /* For a 1D array, do erosion because after filling holes, two single
+ elements outside the range can mask a very large portion of the
+ input. */
+ tmp = ( formask->ndim==1
+ ? gal_binary_erode(tmp, 1, 1, 1)
+ : gal_binary_dilate(tmp, 1, 1, 1) );
+ gal_binary_holes_fill(tmp, formask->ndim, -1);
+ tmp=gal_binary_erode(tmp, 2, formask->ndim, 1);
+ tmp=gal_binary_dilate(tmp, formask->ndim==1?4:2, formask->ndim, 1);
+
+ /* Apply the flag onto the input (to set the pixels to NaN). */
+ gal_blank_flag_apply(work, tmp);
+
+ /* For a check (define 'int check' as an argument, and when calling this
+ function, set 'index==XXX').
+ if(check)
+ {
+ size_t i;
+ double *f=work->array; // CHECK THE TYPE OF YOUR INPUT
+ uint8_t *u=tmp->array;
+ for(i=0;i<work->size;++i)
+ printf("%-5zu %-5u %f\n", i, u[i], f[i]);
+ printf("%s: GOOD\n", __func__); exit(0);
+ }
+ */
+
+ /* Apply the operation: note that 'isfill' is zero here because we don't
+ care about the order of elements any more ('isfill' is only used to do
+ the operation inplace or not). */
+ out=dimension_collapse_sortbased_operation(p, work, clipflags, 0);
+
+ /* Clean up and return. */
+ gal_data_free(tmp);
+ gal_data_free(fstat);
+ gal_data_free(upper);
+ gal_data_free(lower);
+ gal_data_free(multip);
+ gal_data_free(center);
+ gal_data_free(spread);
+ return out;
+}
+
+
+
+
+
+gal_data_t *
+dimension_collapse_sortbased_conv(gal_data_t *work /*, int check*/)
+{
+ float *karr;
+ gal_data_t *in, *out, *kernel;
+ size_t i, ksize[3]={3,3,3}; /* Extra dimensions will not be used. */
+
+ /* Set the input. */
+ in = ( work->type==GAL_TYPE_FLOAT32
+ ? work
+ : gal_data_copy_to_new_type(work, GAL_TYPE_FLOAT32) );
+
+ /* Build a kernel and fill it with equal values (so their sum is 1). */
+ kernel=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, work->ndim, ksize,
+ NULL, 0, -1, 1, NULL, NULL, NULL);
+ karr=kernel->array;
+ for(i=0;i<kernel->size;++i) karr[i]=1.0/kernel->size;
+
+ /* Convolve the work array with the given kernel. */
+ out=gal_convolve_spatial(in, kernel, 1, 1, 0, 1);
+
+ /* For a check.
+ if(check)
+ {
+ float *ia=in->array, *oa=out->array;
+ for(i=0;i<in->size;++i)
+ printf("%-5zu %-10.3f %-10.3f\n", i, ia[i], oa[i]);
+ }
+ */
+
+ /* Clean up and return. */
+ if(in!=work) gal_data_free(in);
+ gal_data_free(kernel);
+ return out;
+}
+
+
+
+
+
static void *
dimension_collapse_sortbased_worker(void *in_prm)
{
@@ -903,14 +1155,32 @@ dimension_collapse_sortbased_worker(void *in_prm)
gal_data_t *in=p->in;
/* Subsequent definitions. */
- gal_data_t *work, *stat=NULL;
+ uint8_t clipflags=0, isfill=0;
size_t a, b, c, sind=GAL_BLANK_SIZE_T;
+ gal_data_t *work, *conv=NULL, *stat=NULL;
size_t i, j, index, c_dim=p->c_dim, wdsize=in->dsize[c_dim];
/* Allocate the dataset that will be sorted. */
work=gal_data_alloc(NULL, in->type, 1, &wdsize, NULL, 0,
p->minmapsize, p->quietmmap, NULL, NULL, NULL);
+ /* Some things are unique for fill-based operators. */
+ switch(p->operator)
+ {
+ case DIMENSION_COLLAPSE_MADCLIP_FILL_MAD:
+ case DIMENSION_COLLAPSE_SIGCLIP_FILL_MAD:
+ case DIMENSION_COLLAPSE_MADCLIP_FILL_STD:
+ case DIMENSION_COLLAPSE_SIGCLIP_FILL_STD:
+ case DIMENSION_COLLAPSE_MADCLIP_FILL_MEAN:
+ case DIMENSION_COLLAPSE_SIGCLIP_FILL_MEAN:
+ case DIMENSION_COLLAPSE_MADCLIP_FILL_MEDIAN:
+ case DIMENSION_COLLAPSE_SIGCLIP_FILL_MEDIAN:
+ case DIMENSION_COLLAPSE_MADCLIP_FILL_NUMBER:
+ case DIMENSION_COLLAPSE_SIGCLIP_FILL_NUMBER:
+ isfill=1;
+ break;
+ }
+
/* Go over all the actions (pixels in this case) that were assigned to
this thread. */
for(i=0; tprm->indexs[i] != GAL_BLANK_SIZE_T; ++i)
@@ -930,14 +1200,15 @@ dimension_collapse_sortbased_worker(void *in_prm)
{
/* One-dimensional data. */
case 1:
- memcpy(work->array, in->array, in->size*gal_type_sizeof(in->type));
+ memcpy(work->array, in->array,
+ in->size*gal_type_sizeof(in->type));
break;
/* Two dimensional data. */
case 2:
a=in->dsize[0];
b=in->dsize[1];
- if(c_dim) /* c_dim==1 The dim. to collapse is already contiguous. */
+ if(c_dim) /* c_dim==1 dim. to collapse, already contiguous. */
memcpy(work->array,
gal_pointer_increment(in->array, index*b, in->type),
b*gal_type_sizeof(in->type));
@@ -959,7 +1230,8 @@ dimension_collapse_sortbased_worker(void *in_prm)
case 1:
for(j=0;j<b;++j)
- dimension_csb_copy(in, (index/c)*b*c+j*c+(index%c), work, j);
+ dimension_csb_copy(in, (index/c)*b*c+j*c+(index%c),
+ work, j);
break;
case 2: /* Fastest dimension: contiguous in memory. */
@@ -984,45 +1256,88 @@ dimension_collapse_sortbased_worker(void *in_prm)
"dimensions", __func__, PACKAGE_BUGREPORT, in->ndim);
}
- /* For a check.
- if(index==0)
+ /* For a check
+ if(index==250)
{
- float *f=work->array;
+ double *f=work->array;
for(j=0;j<wdsize;++j)
- printf("%zu: %f\n", j, f[j]);
- printf("%s: GOOD\n", __func__); exit(0);
+ printf("%zu %f\n", j, f[j]);
+ printf("%s: GOOD\n", __func__); //exit(0);
}
*/
+ /* Set the necessary flag for extra calculation during sigma-clipping
+ (this is not necessary for some). */
+ switch(p->operator)
+ {
+ case DIMENSION_COLLAPSE_MADCLIP_STD:
+ case DIMENSION_COLLAPSE_MADCLIP_FILL_STD:
+ clipflags=GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_STD; break;
+ case DIMENSION_COLLAPSE_SIGCLIP_MAD:
+ case DIMENSION_COLLAPSE_SIGCLIP_FILL_MAD:
+ clipflags=GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_MAD; break;
+ case DIMENSION_COLLAPSE_MADCLIP_MEAN:
+ case DIMENSION_COLLAPSE_SIGCLIP_MEAN:
+ case DIMENSION_COLLAPSE_MADCLIP_FILL_MEAN:
+ case DIMENSION_COLLAPSE_SIGCLIP_FILL_MEAN:
+ clipflags=GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_MEAN; break;
+ }
+
+ /* Create a convolved image (currently disabled, because 'work' will
+ always be non-NULL, kept for further investigation later!). */
+ conv=work?NULL:dimension_collapse_sortbased_conv(work);
+
/* Do the necessary statistical operation. */
+ stat=dimension_collapse_sortbased_operation(p, conv?conv:work,
+ clipflags, isfill);
+
+ /* If this is a "filling" operation, then repeat the operation with
+ the fill. */
+ if(isfill)
+ stat=dimension_collapse_sortbased_fill(p, stat, work, conv,
+ clipflags, index==-1);
+
+ /* Set the index in the output 'stat' array. These can't be set in
+ the main operation 'switch' because the functions are different,
+ while 'sind' is the same. Note also that this should be done after
+ the actual operation because some operators need to correct the
+ type. */
switch(p->operator)
{
case DIMENSION_COLLAPSE_MEDIAN:
- sind=0;
- stat=gal_statistics_median(work, 1);
- break;
+ stat=gal_data_copy_to_new_type_free(stat, GAL_TYPE_FLOAT32);
+ sind=0; break;
+ case DIMENSION_COLLAPSE_MADCLIP_MAD:
+ case DIMENSION_COLLAPSE_SIGCLIP_MAD:
+ case DIMENSION_COLLAPSE_MADCLIP_FILL_MAD:
+ case DIMENSION_COLLAPSE_SIGCLIP_FILL_MAD:
+ sind=GAL_STATISTICS_CLIP_OUTCOL_MAD; break;
+ case DIMENSION_COLLAPSE_MADCLIP_STD:
case DIMENSION_COLLAPSE_SIGCLIP_STD:
+ case DIMENSION_COLLAPSE_MADCLIP_FILL_STD:
+ case DIMENSION_COLLAPSE_SIGCLIP_FILL_STD:
+ sind=GAL_STATISTICS_CLIP_OUTCOL_STD; break;
+ case DIMENSION_COLLAPSE_MADCLIP_MEAN:
case DIMENSION_COLLAPSE_SIGCLIP_MEAN:
+ case DIMENSION_COLLAPSE_MADCLIP_FILL_MEAN:
+ case DIMENSION_COLLAPSE_SIGCLIP_FILL_MEAN:
+ sind=GAL_STATISTICS_CLIP_OUTCOL_MEAN; break;
+ case DIMENSION_COLLAPSE_MADCLIP_MEDIAN:
case DIMENSION_COLLAPSE_SIGCLIP_MEDIAN:
+ case DIMENSION_COLLAPSE_MADCLIP_FILL_MEDIAN:
+ case DIMENSION_COLLAPSE_SIGCLIP_FILL_MEDIAN:
+ sind=GAL_STATISTICS_CLIP_OUTCOL_MEDIAN; break;
+ case DIMENSION_COLLAPSE_MADCLIP_NUMBER:
case DIMENSION_COLLAPSE_SIGCLIP_NUMBER:
- stat=gal_statistics_sigma_clip(work, p->sclipmultip,
- p->sclipparam, 1, 1);
- switch(p->operator)
- {
- case DIMENSION_COLLAPSE_SIGCLIP_STD: sind=3; break;
- case DIMENSION_COLLAPSE_SIGCLIP_MEAN: sind=2; break;
- case DIMENSION_COLLAPSE_SIGCLIP_MEDIAN:
- stat=gal_data_copy_to_new_type_free(stat, in->type);
- sind=1; break;
- case DIMENSION_COLLAPSE_SIGCLIP_NUMBER:
- stat=gal_data_copy_to_new_type_free(stat, GAL_TYPE_UINT32);
- sind=0; break;
- }
- break;
+ case DIMENSION_COLLAPSE_MADCLIP_FILL_NUMBER:
+ case DIMENSION_COLLAPSE_SIGCLIP_FILL_NUMBER:
+ stat=gal_data_copy_to_new_type_free(stat, GAL_TYPE_UINT32);
+ sind=GAL_STATISTICS_CLIP_OUTCOL_NUMBER_USED; break;
default:
error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at '%s' "
- "to fix the problem. The operator code %d isn't "
- "recognized", __func__, PACKAGE_BUGREPORT, p->operator);
+ "to fix the problem. The operator code '%d' is not "
+ "recognized when writing the output", __func__,
+ PACKAGE_BUGREPORT, p->operator);
}
/* Copy the result from the statistics output into the output array
@@ -1075,11 +1390,29 @@ dimension_collapse_sortbased(gal_data_t *in, size_t
c_dim, int operator,
/* The output array (and its type). */
switch(operator)
{
- case DIMENSION_COLLAPSE_MEDIAN: otype=in->type; break;
- case DIMENSION_COLLAPSE_SIGCLIP_STD: otype=GAL_TYPE_FLOAT32; break;
- case DIMENSION_COLLAPSE_SIGCLIP_MEAN: otype=GAL_TYPE_FLOAT32; break;
- case DIMENSION_COLLAPSE_SIGCLIP_MEDIAN: otype=in->type; break;
- case DIMENSION_COLLAPSE_SIGCLIP_NUMBER: otype=GAL_TYPE_UINT32; break;
+ case DIMENSION_COLLAPSE_MADCLIP_NUMBER:
+ case DIMENSION_COLLAPSE_SIGCLIP_NUMBER:
+ case DIMENSION_COLLAPSE_MADCLIP_FILL_NUMBER:
+ case DIMENSION_COLLAPSE_SIGCLIP_FILL_NUMBER:
+ otype=GAL_TYPE_UINT32; break;
+ case DIMENSION_COLLAPSE_MEDIAN:
+ case DIMENSION_COLLAPSE_MADCLIP_MAD:
+ case DIMENSION_COLLAPSE_SIGCLIP_MAD:
+ case DIMENSION_COLLAPSE_MADCLIP_STD:
+ case DIMENSION_COLLAPSE_SIGCLIP_STD:
+ case DIMENSION_COLLAPSE_MADCLIP_MEAN:
+ case DIMENSION_COLLAPSE_SIGCLIP_MEAN:
+ case DIMENSION_COLLAPSE_MADCLIP_MEDIAN:
+ case DIMENSION_COLLAPSE_SIGCLIP_MEDIAN:
+ case DIMENSION_COLLAPSE_MADCLIP_FILL_MAD:
+ case DIMENSION_COLLAPSE_SIGCLIP_FILL_MAD:
+ case DIMENSION_COLLAPSE_MADCLIP_FILL_STD:
+ case DIMENSION_COLLAPSE_SIGCLIP_FILL_STD:
+ case DIMENSION_COLLAPSE_MADCLIP_FILL_MEAN:
+ case DIMENSION_COLLAPSE_SIGCLIP_FILL_MEAN:
+ case DIMENSION_COLLAPSE_MADCLIP_FILL_MEDIAN:
+ case DIMENSION_COLLAPSE_SIGCLIP_FILL_MEDIAN:
+ otype=GAL_TYPE_FLOAT32; break;
default:
error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at '%s' "
"to fix the problem. The operator code %d is not a "
@@ -1129,6 +1462,185 @@ gal_dimension_collapse_median(gal_data_t *in, size_t
c_dim,
+gal_data_t *
+gal_dimension_collapse_mclip_mad(gal_data_t *in, size_t c_dim,
+ float multip, float param,
+ size_t numthreads, size_t minmapsize,
+ int quietmmap)
+{
+ int op=DIMENSION_COLLAPSE_MADCLIP_MAD;
+ return dimension_collapse_sortbased(in, c_dim, op, multip, param,
+ numthreads, minmapsize, quietmmap);
+}
+
+
+
+
+
+gal_data_t *
+gal_dimension_collapse_mclip_fill_mad(gal_data_t *in, size_t c_dim,
+ float multip, float param,
+ size_t numthreads, size_t minmapsize,
+ int quietmmap)
+{
+ int op=DIMENSION_COLLAPSE_MADCLIP_FILL_MAD;
+ return dimension_collapse_sortbased(in, c_dim, op, multip, param,
+ numthreads, minmapsize, quietmmap);
+}
+
+
+
+
+
+gal_data_t *
+gal_dimension_collapse_mclip_std(gal_data_t *in, size_t c_dim,
+ float multip, float param,
+ size_t numthreads, size_t minmapsize,
+ int quietmmap)
+{
+ int op=DIMENSION_COLLAPSE_MADCLIP_STD;
+ return dimension_collapse_sortbased(in, c_dim, op, multip, param,
+ numthreads, minmapsize, quietmmap);
+}
+
+
+
+
+
+gal_data_t *
+gal_dimension_collapse_mclip_fill_std(gal_data_t *in, size_t c_dim,
+ float multip, float param,
+ size_t numthreads, size_t minmapsize,
+ int quietmmap)
+{
+ int op=DIMENSION_COLLAPSE_MADCLIP_FILL_STD;
+ return dimension_collapse_sortbased(in, c_dim, op, multip, param,
+ numthreads, minmapsize, quietmmap);
+}
+
+
+
+
+
+gal_data_t *
+gal_dimension_collapse_mclip_mean(gal_data_t *in, size_t c_dim,
+ float multip, float param,
+ size_t numthreads, size_t minmapsize,
+ int quietmmap)
+{
+ int op=DIMENSION_COLLAPSE_MADCLIP_MEAN;
+ return dimension_collapse_sortbased(in, c_dim, op, multip, param,
+ numthreads, minmapsize, quietmmap);
+}
+
+
+
+
+
+gal_data_t *
+gal_dimension_collapse_mclip_fill_mean(gal_data_t *in, size_t c_dim,
+ float multip, float param,
+ size_t numthreads, size_t minmapsize,
+ int quietmmap)
+{
+ int op=DIMENSION_COLLAPSE_MADCLIP_FILL_MEAN;
+ return dimension_collapse_sortbased(in, c_dim, op, multip, param,
+ numthreads, minmapsize, quietmmap);
+}
+
+
+
+
+
+gal_data_t *
+gal_dimension_collapse_mclip_median(gal_data_t *in, size_t c_dim,
+ float multip, float param,
+ size_t numthreads, size_t minmapsize,
+ int quietmmap)
+{
+ int op=DIMENSION_COLLAPSE_MADCLIP_MEDIAN;
+ return dimension_collapse_sortbased(in, c_dim, op, multip, param,
+ numthreads, minmapsize, quietmmap);
+}
+
+
+
+
+
+gal_data_t *
+gal_dimension_collapse_mclip_fill_median(gal_data_t *in, size_t c_dim,
+ float multip, float param,
+ size_t numthreads,
+ size_t minmapsize, int quietmmap)
+{
+ int op=DIMENSION_COLLAPSE_MADCLIP_FILL_MEDIAN;
+ return dimension_collapse_sortbased(in, c_dim, op, multip, param,
+ numthreads, minmapsize, quietmmap);
+}
+
+
+
+
+
+gal_data_t *
+gal_dimension_collapse_mclip_number(gal_data_t *in, size_t c_dim,
+ float multip, float param,
+ size_t numthreads, size_t minmapsize,
+ int quietmmap)
+{
+ int op=DIMENSION_COLLAPSE_MADCLIP_NUMBER;
+ return dimension_collapse_sortbased(in, c_dim, op, multip, param,
+ numthreads, minmapsize, quietmmap);
+}
+
+
+
+
+
+gal_data_t *
+gal_dimension_collapse_mclip_fill_number(gal_data_t *in, size_t c_dim,
+ float multip, float param,
+ size_t numthreads,
+ size_t minmapsize, int quietmmap)
+{
+ int op=DIMENSION_COLLAPSE_MADCLIP_FILL_NUMBER;
+ return dimension_collapse_sortbased(in, c_dim, op, multip, param,
+ numthreads, minmapsize, quietmmap);
+}
+
+
+
+
+
+gal_data_t *
+gal_dimension_collapse_sclip_mad(gal_data_t *in, size_t c_dim,
+ float multip, float param,
+ size_t numthreads, size_t minmapsize,
+ int quietmmap)
+{
+ int op=DIMENSION_COLLAPSE_SIGCLIP_MAD;
+ return dimension_collapse_sortbased(in, c_dim, op, multip, param,
+ numthreads, minmapsize, quietmmap);
+}
+
+
+
+
+
+gal_data_t *
+gal_dimension_collapse_sclip_fill_mad(gal_data_t *in, size_t c_dim,
+ float multip, float param,
+ size_t numthreads, size_t minmapsize,
+ int quietmmap)
+{
+ int op=DIMENSION_COLLAPSE_SIGCLIP_FILL_MAD;
+ return dimension_collapse_sortbased(in, c_dim, op, multip, param,
+ numthreads, minmapsize, quietmmap);
+}
+
+
+
+
gal_data_t *
gal_dimension_collapse_sclip_std(gal_data_t *in, size_t c_dim,
@@ -1145,6 +1657,21 @@ gal_dimension_collapse_sclip_std(gal_data_t *in, size_t
c_dim,
+gal_data_t *
+gal_dimension_collapse_sclip_fill_std(gal_data_t *in, size_t c_dim,
+ float multip, float param,
+ size_t numthreads, size_t minmapsize,
+ int quietmmap)
+{
+ int op=DIMENSION_COLLAPSE_SIGCLIP_FILL_STD;
+ return dimension_collapse_sortbased(in, c_dim, op, multip, param,
+ numthreads, minmapsize, quietmmap);
+}
+
+
+
+
+
gal_data_t *
gal_dimension_collapse_sclip_mean(gal_data_t *in, size_t c_dim,
float multip, float param,
@@ -1160,6 +1687,21 @@ gal_dimension_collapse_sclip_mean(gal_data_t *in, size_t
c_dim,
+gal_data_t *
+gal_dimension_collapse_sclip_fill_mean(gal_data_t *in, size_t c_dim,
+ float multip, float param,
+ size_t numthreads, size_t minmapsize,
+ int quietmmap)
+{
+ int op=DIMENSION_COLLAPSE_SIGCLIP_FILL_MEAN;
+ return dimension_collapse_sortbased(in, c_dim, op, multip, param,
+ numthreads, minmapsize, quietmmap);
+}
+
+
+
+
+
gal_data_t *
gal_dimension_collapse_sclip_median(gal_data_t *in, size_t c_dim,
float multip, float param,
@@ -1171,6 +1713,25 @@ gal_dimension_collapse_sclip_median(gal_data_t *in,
size_t c_dim,
numthreads, minmapsize, quietmmap);
}
+
+
+
+
+gal_data_t *
+gal_dimension_collapse_sclip_fill_median(gal_data_t *in, size_t c_dim,
+ float multip, float param,
+ size_t numthreads,
+ size_t minmapsize, int quietmmap)
+{
+ int op=DIMENSION_COLLAPSE_SIGCLIP_FILL_MEDIAN;
+ return dimension_collapse_sortbased(in, c_dim, op, multip, param,
+ numthreads, minmapsize, quietmmap);
+}
+
+
+
+
+
gal_data_t *
gal_dimension_collapse_sclip_number(gal_data_t *in, size_t c_dim,
float multip, float param,
@@ -1186,6 +1747,22 @@ gal_dimension_collapse_sclip_number(gal_data_t *in,
size_t c_dim,
+gal_data_t *
+gal_dimension_collapse_sclip_fill_number(gal_data_t *in, size_t c_dim,
+ float multip, float param,
+ size_t numthreads,
+ size_t minmapsize, int quietmmap)
+{
+ int op=DIMENSION_COLLAPSE_SIGCLIP_FILL_NUMBER;
+ return dimension_collapse_sortbased(in, c_dim, op, multip, param,
+ numthreads, minmapsize, quietmmap);
+}
+
+
+
+
+
+
diff --git a/lib/gnuastro/arithmetic.h b/lib/gnuastro/arithmetic.h
index 31881e8b..66933e22 100644
--- a/lib/gnuastro/arithmetic.h
+++ b/lib/gnuastro/arithmetic.h
@@ -182,12 +182,29 @@ enum gal_arithmetic_operators
GAL_ARITHMETIC_OP_SUM, /* Sum per pixel of multiple arrays. */
GAL_ARITHMETIC_OP_MEAN, /* Mean per pixel of multiple arrays. */
GAL_ARITHMETIC_OP_STD, /* STD per pixel of multiple arrays. */
+ GAL_ARITHMETIC_OP_MAD, /* MAD per pixel of multiple arrays. */
GAL_ARITHMETIC_OP_MEDIAN, /* Median per pixel of multiple arrays. */
GAL_ARITHMETIC_OP_QUANTILE, /* Quantile per pixel of multiple arrays.*/
GAL_ARITHMETIC_OP_SIGCLIP_NUMBER,/* Sigma-clipped number of mult. arrays.*/
GAL_ARITHMETIC_OP_SIGCLIP_MEAN, /* Sigma-clipped mean of multiple arrays.*/
GAL_ARITHMETIC_OP_SIGCLIP_MEDIAN,/* Sigma-clipped median of mult. arrays.*/
GAL_ARITHMETIC_OP_SIGCLIP_STD, /* Sigma-clipped STD of multiple arrays. */
+ GAL_ARITHMETIC_OP_SIGCLIP_MAD, /* Sigma-clipped STD of multiple arrays. */
+ GAL_ARITHMETIC_OP_SIGCLIP_FILL_NUMBER, /* MAD-clipped num. of arrays. */
+ GAL_ARITHMETIC_OP_SIGCLIP_FILL_MEAN, /* MAD-clipped mean of arrays. */
+ GAL_ARITHMETIC_OP_SIGCLIP_FILL_MEDIAN, /* MAD-clipped median of arrays. */
+ GAL_ARITHMETIC_OP_SIGCLIP_FILL_STD, /* MAD-clipped STD of arrays. */
+ GAL_ARITHMETIC_OP_SIGCLIP_FILL_MAD, /* MAD-clipped STD of arrays. */
+ GAL_ARITHMETIC_OP_MADCLIP_NUMBER,/* MAD-clipped number of mult. arrays. */
+ GAL_ARITHMETIC_OP_MADCLIP_MEAN, /* MAD-clipped mean of multiple arrays. */
+ GAL_ARITHMETIC_OP_MADCLIP_MEDIAN,/* MAD-clipped median of mult. arrays. */
+ GAL_ARITHMETIC_OP_MADCLIP_STD, /* MAD-clipped STD of multiple arrays. */
+ GAL_ARITHMETIC_OP_MADCLIP_MAD, /* MAD-clipped STD of multiple arrays. */
+ GAL_ARITHMETIC_OP_MADCLIP_FILL_NUMBER,/* MAD-clipped num. of arrays. */
+ GAL_ARITHMETIC_OP_MADCLIP_FILL_MEAN, /* MAD-clipped mean of arrays. */
+ GAL_ARITHMETIC_OP_MADCLIP_FILL_MEDIAN,/* MAD-clipped median of arrays. */
+ GAL_ARITHMETIC_OP_MADCLIP_FILL_STD, /* MAD-clipped STD of arrays. */
+ GAL_ARITHMETIC_OP_MADCLIP_FILL_MAD, /* MAD-clipped STD of arrays. */
GAL_ARITHMETIC_OP_MKNOISE_SIGMA,/* Fixed-sigma noise to every element. */
GAL_ARITHMETIC_OP_MKNOISE_SIGMA_FROM_MEAN, /* Sigma comes from background. */
diff --git a/lib/gnuastro/dimension.h b/lib/gnuastro/dimension.h
index 6e8426a2..330b9d42 100644
--- a/lib/gnuastro/dimension.h
+++ b/lib/gnuastro/dimension.h
@@ -145,30 +145,128 @@ gal_dimension_collapse_median(gal_data_t *in, size_t
c_dim,
size_t numthreads, size_t minmapsize,
int quietmmap);
+gal_data_t *
+gal_dimension_collapse_mclip_mad(gal_data_t *in, size_t c_dim,
+ float multip, float param,
+ size_t numthreads, size_t minmapsize,
+ int quietmmap);
+
+gal_data_t *
+gal_dimension_collapse_mclip_fill_mad(gal_data_t *in, size_t c_dim,
+ float multip, float param,
+ size_t numthreads, size_t minmapsize,
+ int quietmmap);
+
+gal_data_t *
+gal_dimension_collapse_mclip_std(gal_data_t *in, size_t c_dim,
+ float multip, float param,
+ size_t numthreads, size_t minmapsize,
+ int quietmmap);
+
+gal_data_t *
+gal_dimension_collapse_mclip_fill_std(gal_data_t *in, size_t c_dim,
+ float multip, float param,
+ size_t numthreads, size_t minmapsize,
+ int quietmmap);
+
+gal_data_t *
+gal_dimension_collapse_mclip_mean(gal_data_t *in, size_t c_dim,
+ float multip, float param,
+ size_t numthreads, size_t minmapsize,
+ int quietmmap);
+
+gal_data_t *
+gal_dimension_collapse_mclip_fill_mean(gal_data_t *in, size_t c_dim,
+ float multip, float param,
+ size_t numthreads,
+ size_t minmapsize, int quietmmap);
+
+gal_data_t *
+gal_dimension_collapse_mclip_median(gal_data_t *in, size_t c_dim,
+ float multip, float param,
+ size_t numthreads, size_t minmapsize,
+ int quietmmap);
+
+gal_data_t *
+gal_dimension_collapse_mclip_fill_median(gal_data_t *in, size_t c_dim,
+ float multip, float param,
+ size_t numthreads,
+ size_t minmapsize, int quietmmap);
+
+gal_data_t *
+gal_dimension_collapse_mclip_number(gal_data_t *in, size_t c_dim,
+ float multip, float param,
+ size_t numthreads, size_t minmapsize,
+ int quietmmap);
+
+gal_data_t *
+gal_dimension_collapse_mclip_fill_number(gal_data_t *in, size_t c_dim,
+ float multip, float param,
+ size_t numthreads,
+ size_t minmapsize, int quietmmap);
+
+gal_data_t *
+gal_dimension_collapse_sclip_mad(gal_data_t *in, size_t c_dim,
+ float multip, float param,
+ size_t numthreads, size_t minmapsize,
+ int quietmmap);
+
+gal_data_t *
+gal_dimension_collapse_sclip_fill_mad(gal_data_t *in, size_t c_dim,
+ float multip, float param,
+ size_t numthreads, size_t minmapsize,
+ int quietmmap);
+
gal_data_t *
gal_dimension_collapse_sclip_std(gal_data_t *in, size_t c_dim,
float multip, float param,
size_t numthreads, size_t minmapsize,
int quietmmap);
+gal_data_t *
+gal_dimension_collapse_sclip_fill_std(gal_data_t *in, size_t c_dim,
+ float multip, float param,
+ size_t numthreads,
+ size_t minmapsize, int quietmmap);
+
gal_data_t *
gal_dimension_collapse_sclip_mean(gal_data_t *in, size_t c_dim,
float multip, float param,
size_t numthreads, size_t minmapsize,
int quietmmap);
+gal_data_t *
+gal_dimension_collapse_sclip_fill_mean(gal_data_t *in, size_t c_dim,
+ float multip, float param,
+ size_t numthreads,
+ size_t minmapsize, int quietmmap);
+
gal_data_t *
gal_dimension_collapse_sclip_median(gal_data_t *in, size_t c_dim,
float multip, float param,
size_t numthreads, size_t minmapsize,
int quietmmap);
+gal_data_t *
+gal_dimension_collapse_sclip_fill_median(gal_data_t *in, size_t c_dim,
+ float multip, float param,
+ size_t numthreads,
+ size_t minmapsize, int quietmmap);
+
gal_data_t *
gal_dimension_collapse_sclip_number(gal_data_t *in, size_t c_dim,
float multip, float param,
size_t numthreads, size_t minmapsize,
int quietmmap);
+gal_data_t *
+gal_dimension_collapse_sclip_fill_number(gal_data_t *in, size_t c_dim,
+ float multip, float param,
+ size_t numthreads,
+ size_t minmapsize, int quietmmap);
+
+
+
/************************************************************************/
@@ -179,6 +277,8 @@ gal_dimension_remove_extra(size_t ndim, size_t *dsize,
struct wcsprm *wcs);
+
+
/************************************************************************/
/******************** Neighbors **********************/
/************************************************************************/
diff --git a/lib/gnuastro/statistics.h b/lib/gnuastro/statistics.h
index e9fddf33..4d992a72 100644
--- a/lib/gnuastro/statistics.h
+++ b/lib/gnuastro/statistics.h
@@ -48,7 +48,7 @@ __BEGIN_C_DECLS /* From C++ preparations */
/* Maximum number of tests for sigma-clipping convergence. */
-#define GAL_STATISTICS_SIG_CLIP_MAX_CONVERGE 50
+#define GAL_STATISTICS_CLIP_MAX_CONVERGE 50
/* Least acceptable mode symmetricity.*/
#define GAL_STATISTICS_MODE_GOOD_SYM 0.2f
@@ -57,7 +57,7 @@ __BEGIN_C_DECLS /* From C++ preparations */
-enum bin_status
+enum gal_statistics_bin_status
{
GAL_STATISTICS_BINS_INVALID, /* ==0 by C standard. */
@@ -65,6 +65,23 @@ enum bin_status
GAL_STATISTICS_BINS_IRREGULAR,
};
+enum gal_statistics_clip_outcol
+{
+ GAL_STATISTICS_CLIP_OUTCOL_NUMBER_USED, /* =0 by C standard. */
+ GAL_STATISTICS_CLIP_OUTCOL_MEAN,
+ GAL_STATISTICS_CLIP_OUTCOL_STD,
+ GAL_STATISTICS_CLIP_OUTCOL_MEDIAN,
+ GAL_STATISTICS_CLIP_OUTCOL_MAD,
+ GAL_STATISTICS_CLIP_OUTCOL_NUMBER_CLIPS,
+
+ GAL_STATISTICS_CLIP_OUT_SIZE,
+};
+
+/* The optional measurements to do after sigma-clipping. */
+#define GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_MEAN 0x1
+#define GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_STD 0x2
+#define GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_MAD 0x4
+
/****************************************************************
******** Simple statistics *******
@@ -97,6 +114,12 @@ gal_statistics_std_from_sums(double sum, double sump2,
size_t num);
gal_data_t *
gal_statistics_median(gal_data_t *input, int inplace);
+gal_data_t *
+gal_statistics_mad(gal_data_t *input, int inplace);
+
+gal_data_t *
+gal_statistics_median_mad(gal_data_t *input, int inplace);
+
size_t
gal_statistics_quantile_index(size_t size, double quantile);
@@ -180,14 +203,19 @@ gal_statistics_cfp(gal_data_t *data, gal_data_t *bins,
int normalize);
***************** Outliers ********************
****************************************************************/
gal_data_t *
-gal_statistics_sigma_clip(gal_data_t *input, float multip, float param,
- int inplace, int quiet);
+gal_statistics_clip_sigma(gal_data_t *input, float multip, float param,
+ uint8_t extrastats, int inplace, int quiet);
+
+gal_data_t *
+gal_statistics_clip_mad(gal_data_t *input, float multip, float param,
+ uint8_t extrastats, int inplace, int quiet);
gal_data_t *
gal_statistics_outlier_bydistance(int pos1_neg0, gal_data_t *input,
size_t window_size, float sigma,
- float sigclip_multip, float sigclip_param,
- int inplace, int quiet);
+ float sigclip_multip,
+ float sigclip_param, int inplace,
+ int quiet);
gal_data_t *
gal_statistics_outlier_flat_cfp(gal_data_t *input, size_t numprev,
diff --git a/lib/statistics.c b/lib/statistics.c
index c7913e57..bb4e5ad8 100644
--- a/lib/statistics.c
+++ b/lib/statistics.c
@@ -388,6 +388,139 @@ gal_statistics_median(gal_data_t *input, int inplace)
+
+static void
+statistics_mad_in_sorted_no_blank(gal_data_t *sorted, gal_data_t *med,
+ void *mad_o)
+{
+ uint8_t type;
+ gal_data_t *use, *mad;
+ int flags=GAL_ARITHMETIC_FLAG_INPLACE | GAL_ARITHMETIC_FLAG_NUMOK;
+
+ /* Sanity check. */
+ if(med->type!=sorted->type)
+ error(EXIT_FAILURE, 0, "%s: the input 'sorted' and 'med' arrays "
+ "do not have the same type; they are respectively '%s' and '%s'",
+ __func__, gal_type_name(sorted->type, 1),
+ gal_type_name(med->type, 1));
+
+ /* After subtracting, we will need to sort the array, so a copy is
+ necessary (the input should not be touched). Furthermore, if the input
+ is an un-signed integer, convert it to a signed integer of the next
+ larger size. This is necessary, because half of the values will become
+ negative after subtracting the median. */
+ switch(sorted->type)
+ {
+ case GAL_TYPE_UINT8: type=GAL_TYPE_INT16;
+ case GAL_TYPE_UINT16: type=GAL_TYPE_INT32;
+ case GAL_TYPE_UINT32: type=GAL_TYPE_INT64;
+ case GAL_TYPE_UINT64: type=GAL_TYPE_INT64;
+ default: type=GAL_TYPE_INVALID; /* Not necessary. */
+ }
+ use=gal_data_copy_to_new_type(sorted, ( type!=GAL_TYPE_INVALID
+ ? type : sorted->type) );
+
+ /* Subtract the median from the input. */
+ use=gal_arithmetic(GAL_ARITHMETIC_OP_MINUS, 1, flags, use, med);
+
+ /* Get the absolute value of the differences from the median. The
+ absolute values of the differences can fit into the original input
+ type, so to make things consistent, we'll take it back to the original
+ type. */
+ use=gal_arithmetic(GAL_ARITHMETIC_OP_ABS, 1, flags, use);
+ use=gal_data_copy_to_new_type_free(use, sorted->type);
+ use->flag=0;/* Necessary before new median to re-sort. */
+ mad = gal_statistics_median(use, 1);
+
+ /* For a check:
+ {
+ size_t i;
+ double *u=use->array;
+ double *ma=med->array, *Ma=mad->array, *s=sorted->array;
+ for(i=0;i<sorted->size;++i)
+ printf("%-15g %-15g\n", s[i], u[i]);
+ printf("Median: %g\n", ma[0]);
+ printf("MAD: %g\n", Ma[0]);
+ exit(0);
+ } */
+
+ /* Copy the MAD value into the output pointer. */
+ memcpy(mad_o, mad->array, gal_type_sizeof(mad->type));
+
+ /* Clean up. */
+ gal_data_free(mad);
+ gal_data_free(use);
+}
+
+
+
+
+
+/* Return the median and median absolute deviation. */
+static gal_data_t *
+statistics_median_mad(gal_data_t *input, int inplace, int onlymad)
+{
+ size_t one=1, two=2;
+ gal_data_t *in, *med;
+ gal_data_t *mad, *out;
+
+ /* If the caller only wants the MAD, then the output should only have one
+ element (which is the actual 'mad' that is calculated). */
+ mad = gal_data_alloc(NULL, input->type, 1, &one, NULL, 1,
+ -1, 1, NULL, NULL, NULL);
+ out = ( onlymad
+ ? mad
+ : gal_data_alloc(NULL, input->type, 1, &two, NULL, 1,
+ -1, 1, NULL, NULL, NULL) );
+
+ /* Allocate the input array if we should not work in-place. */
+ in = inplace ? input : gal_data_copy(input);
+
+ /* Calculate the median. */
+ med = gal_statistics_median(in, 1);
+
+ /* Write the MAD into the allocated space. */
+ statistics_mad_in_sorted_no_blank(in, med, mad->array);
+
+ /* If the caller wanted both the median and the MAD, write the median and
+ MAD into the output dataset. */
+ if(onlymad==0)
+ {
+ memcpy(out->array, med->array, gal_type_sizeof(med->type));
+ memcpy(gal_pointer_increment(out->array, 1, out->type), mad->array,
+ gal_type_sizeof(out->type));
+ gal_data_free(mad);
+ }
+
+ /* Clean up and return. */
+ gal_data_free(med);
+ return out;
+}
+
+
+
+
+
+gal_data_t *
+gal_statistics_mad(gal_data_t *input, int inplace)
+{
+ return statistics_median_mad(input, inplace, 1);
+}
+
+
+
+
+
+gal_data_t *
+gal_statistics_median_mad(gal_data_t *input, int inplace)
+{
+ return statistics_median_mad(input, inplace, 0);
+}
+
+
+
+
+
/* For a given size, return the index (starting from zero) that is at the
given quantile. */
size_t
@@ -2266,98 +2399,190 @@ gal_statistics_cfp(gal_data_t *input, gal_data_t
*bins, int normalize)
/****************************************************************
***************** Outliers ********************
****************************************************************/
-/* Sigma-cilp a given distribution:
+static gal_data_t *
+statistics_clip_prepare(gal_data_t *input, gal_data_t *nbs, float multip,
+ float param, int quiet, int sig1_mad0,
+ gal_data_t **center, gal_data_t **spread,
+ char **colnames)
+{
+ float *oa;
+ gal_data_t *out;
+ uint8_t type=gal_tile_block(input)->type;
+ size_t i, one=1, osize=GAL_STATISTICS_CLIP_OUT_SIZE;
- Inputs:
+ /* Some sanity checks. */
+ if( multip<=0 )
+ error(EXIT_FAILURE, 0, "%s: 'multip', must be greater than zero. The "
+ "given value was %g", __func__, multip);
+ if( param<=0 )
+ error(EXIT_FAILURE, 0, "%s: 'param', must be greater than zero. The "
+ "given value was %g", __func__, param);
+ if( param >= 1.0f && ceil(param) != param )
+ error(EXIT_FAILURE, 0, "%s: when 'param' is larger than 1.0, it is "
+ "interpretted as an absolute number of clips. So it must be an "
+ "integer. However, your given value %g", __func__, param);
+ if( (nbs->flag & GAL_DATA_FLAG_SORT_CH)==0 )
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix the "
+ "problem. 'nbs->flag', doesn't have the 'GAL_DATA_FLAG_SORT_CH' "
+ "bit activated", __func__, PACKAGE_BUGREPORT);
+ if( (nbs->flag & GAL_DATA_FLAG_SORTED_I)==0
+ && (nbs->flag & GAL_DATA_FLAG_SORTED_D)==0 )
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix the "
+ "problem. 'nbs' isn't sorted", __func__, PACKAGE_BUGREPORT);
- - 'multip': multiple of the standard deviation,
+ /* Allocate the necessary spaces (spread is only necessary for MAD). */
+ out=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, 1, &osize, NULL, 0,
+ input->minmapsize, input->quietmmap, NULL, NULL, NULL);
+ *center=gal_data_alloc(NULL, type, 1, &one, NULL, 0, input->minmapsize,
+ input->quietmmap, NULL, NULL, NULL);
+ *spread = ( sig1_mad0
+ ? NULL
+ : gal_data_alloc(NULL, type, 1, &one, NULL, 0,
+ input->minmapsize, input->quietmmap,
+ NULL, NULL, NULL) );
+
+ /* Set all the output values to NaN to start with. */
+ oa=out->array;
+ for(i=0;i<GAL_STATISTICS_CLIP_OUT_SIZE;++i) oa[i]=NAN;
- - 'param' must be positive and determines the type of clipping:
+ /* Prepare the column names if the user gave quiet=0. */
+ if(quiet==0)
+ {
+ if(sig1_mad0)
+ {
+ if( asprintf(colnames, "%-5s %-10s %-12s %-12s",
+ "round", "number", "median", "STD")<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation1 error",
+ __func__);
+ }
+ else
+ {
+ if(asprintf(colnames, "%-5s %-10s %-12s %-12s",
+ "round", "number", "median", "MAD")<0)
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation2 error",
+ __func__);
+ }
+ }
- - param<1.0: interpretted as a tolerance level to stop clipping.
+ /* Return the allocated space for the output. */
+ return out;
+}
- - param>=1.0 and an integer: a specific number of times to do the
- clippping.
- Output elements (type FLOAT32):
- - 0: Number of points used.
- - 1: Median.
- - 2: Mean.
- - 3: Standard deviation.
- The way this function works is very simple: first it will sort the input
- (if it isn't sorted). Afterwards, it will recursively change the starting
- point of the array and its size, calcluating the basic statistics in each
- round to define the new starting point and size.
-*/
-#define SIGCLIP(IT) { \
+
+/* Calculate all the extra statistics that are usually useful with
+ sigma-clipping. */
+static void
+statistics_clip_stats_extra(gal_data_t *nbs, float *oa, uint8_t extrastats)
+{
+ gal_data_t *tmp;
+ uint8_t istd = extrastats & GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_STD;
+ uint8_t imad = extrastats & GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_MAD;
+ uint8_t imean = extrastats & GAL_STATISTICS_CLIP_OUTCOL_OPTIONAL_MEAN;
+
+ /* If the "extra" stats are already calculated (for example MAD in
+ MAD-clipping), then there is no need to re-calculate it, so set its
+ conditional variable to 0. Note the '!' at the start of the
+ condition. */
+ if( !(isnan(oa[GAL_STATISTICS_CLIP_OUTCOL_MEAN]) && imean) ) imean=0;
+ if( !(isnan(oa[GAL_STATISTICS_CLIP_OUTCOL_STD]) && istd) ) istd=0;
+ if( !(isnan(oa[GAL_STATISTICS_CLIP_OUTCOL_MAD]) && imad) ) imad=0;
+
+ /* Mean and Standard deviation. */
+ if(imean && istd)
+ {
+ tmp=gal_statistics_mean_std(nbs);
+ oa[ GAL_STATISTICS_CLIP_OUTCOL_STD ] = ((double *)(tmp->array))[1];
+ oa[ GAL_STATISTICS_CLIP_OUTCOL_MEAN ] = ((double *)(tmp->array))[0];
+ gal_data_free(tmp);
+ }
+ else /* Only one of the mean or STD was requested */
+ {
+ if(imean)
+ {
+ tmp=gal_statistics_mean(nbs);
+ oa[ GAL_STATISTICS_CLIP_OUTCOL_MEAN ]
+ = ((double *)(tmp->array))[0];
+ gal_data_free(tmp);
+ }
+ if(istd)
+ {
+ tmp=gal_statistics_std(nbs);
+ oa[ GAL_STATISTICS_CLIP_OUTCOL_STD ]
+ = ((double *)(tmp->array))[0];
+ gal_data_free(tmp);
+ }
+ }
+
+ /* MAD. */
+ if(imad)
+ {
+ tmp=gal_statistics_mad(nbs, 1);
+ tmp=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT32);
+ oa[ GAL_STATISTICS_CLIP_OUTCOL_MAD ] = ((float *)(tmp->array))[0];
+ gal_data_free(tmp);
+ }
+}
+
+
+
+
+
+/* Sigma-cilp a given distribution. The way this function works is very
+ simple: first it will sort the input (if it isn't sorted). Afterwards,
+ it will recursively change the starting point of the array and its size,
+ calcluating the basic statistics in each round to define the new
+ starting point and size. */
+#define CLIPALL(IT) { \
IT *a = nbs->array, *af = a + nbs->size; \
IT *bf = nbs->array, *b = bf + nbs->size - 1; \
\
/* Remove all out-of-range elements from the start of the array. */ \
if( nbs->flag & GAL_DATA_FLAG_SORTED_I ) \
- do if( *a > (*med - (multip * *std)) ) \
+ do if( *a > (center - (multip * spread)) ) \
{ start=a; break; } \
while(++a<af); \
else \
- do if( *a < (*med + (multip * *std)) ) \
+ do if( *a < (center + (multip * spread)) ) \
{ start=a; break; } \
while(++a<af); \
\
/* Remove all out-of-range elements from the end of the array. */ \
if( nbs->flag & GAL_DATA_FLAG_SORTED_I ) \
- do if( *b < (*med + (multip * *std)) ) \
+ do if( *b < (center + (multip * spread)) ) \
{ size=b-a+1; break; } \
while(--b>=bf); \
else \
- do if( *b > (*med - (multip * *std)) ) \
+ do if( *b > (center - (multip * spread)) ) \
{ size=b-a+1; break; } \
while(--b>=bf); \
}
-gal_data_t *
-gal_statistics_sigma_clip(gal_data_t *input, float multip, float param,
- int inplace, int quiet)
+static gal_data_t *
+statistics_clip(gal_data_t *input, float multip, float param,
+ uint8_t extrastats, int inplace, int quiet, int sig1_mad0)
{
float *oa;
+ char *colnames;
+ gal_data_t *spread_d;
void *start, *nbs_array;
- double *med, *mean, *std;
+ size_t i, num=0, size, oldsize;
uint8_t type=gal_tile_block(input)->type;
uint8_t bytolerance = param>=1.0f ? 0 : 1;
- double oldmed=NAN, oldmean=NAN, oldstd=NAN;
- size_t num=0, one=1, four=4, size, oldsize;
- gal_data_t *fcopy, *median_i, *median_d, *out, *meanstd;
+ double center=NAN, spread=NAN, oldspread=NAN;
+ gal_data_t *fcopy, *center_i, *center_d, *spread_i, *out;
gal_data_t *nbs=gal_statistics_no_blank_sorted(input, inplace);
- size_t maxnum = param>=1.0f ? param : GAL_STATISTICS_SIG_CLIP_MAX_CONVERGE;
-
- /* Some sanity checks. */
- if( multip<=0 )
- error(EXIT_FAILURE, 0, "%s: 'multip', must be greater than zero. The "
- "given value was %g", __func__, multip);
- if( param<=0 )
- error(EXIT_FAILURE, 0, "%s: 'param', must be greater than zero. The "
- "given value was %g", __func__, param);
- if( param >= 1.0f && ceil(param) != param )
- error(EXIT_FAILURE, 0, "%s: when 'param' is larger than 1.0, it is "
- "interpretted as an absolute number of clips. So it must be an "
- "integer. However, your given value %g", __func__, param);
- if( (nbs->flag & GAL_DATA_FLAG_SORT_CH)==0 )
- error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix the "
- "problem. 'nbs->flag', doesn't have the 'GAL_DATA_FLAG_SORT_CH' "
- "bit activated", __func__, PACKAGE_BUGREPORT);
- if( (nbs->flag & GAL_DATA_FLAG_SORTED_I)==0
- && (nbs->flag & GAL_DATA_FLAG_SORTED_D)==0 )
- error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix the "
- "problem. 'nbs' isn't sorted", __func__, PACKAGE_BUGREPORT);
+ size_t maxnum = param>=1.0f?param:GAL_STATISTICS_CLIP_MAX_CONVERGE;
+ /* Do sanity checks and allocate space for the output. */
+ out=statistics_clip_prepare(input, nbs, multip, param, quiet, sig1_mad0,
+ ¢er_i, &spread_i, &colnames);
- /* Allocate the necessary spaces. */
- out=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, 1, &four, NULL, 0,
- input->minmapsize, input->quietmmap, NULL, NULL, NULL);
- median_i=gal_data_alloc(NULL, type, 1, &one, NULL, 0, input->minmapsize,
- input->quietmmap, NULL, NULL, NULL);
-
+ /* If we have more than one element, and the user wants to see the
+ progress, then print the column information. */
+ if(!quiet && nbs->size>1) { printf("%s\n", colnames); free(colnames); }
/* Only continue processing if we have non-blank elements. */
oa=out->array;
@@ -2367,41 +2592,37 @@ gal_statistics_sigma_clip(gal_data_t *input, float
multip, float param,
/* There was nothing in the input! */
case 0:
if(!quiet)
- printf("NO SIGMA-CLIPPING: all input elements are blank or input's "
- "size is zero.\n");
- oa[0] = oa[1] = oa[2] = oa[3] = NAN;
+ error(EXIT_SUCCESS, 0, "NO %s-CLIPPING: all input elements "
+ "are blank or input's size is zero",
+ sig1_mad0 ? "SIGMA" : "MAD");
+ for(i=0;i<GAL_STATISTICS_CLIP_OUT_SIZE;++i) oa[i]=NAN;
break;
/* Only one element, convert it to floating point and put it as the
mean and median (the standard deviation will be zero by
definition). */
case 1:
- /* Write the values. */
+
+ /* Write the values in the output array. */
fcopy=gal_data_copy_to_new_type(nbs, GAL_TYPE_FLOAT32);
- oa[0] = 1;
- oa[1] = *((float *)(fcopy->array));
- oa[2] = *((float *)(fcopy->array));
- oa[3] = 0;
+ center=*((float *)(fcopy->array));
gal_data_free(fcopy);
-
- /* Print the comments if requested. */
+ spread=0;
+ size=1;
+ oa[ GAL_STATISTICS_CLIP_OUTCOL_MEDIAN ] = center;
+ oa[ GAL_STATISTICS_CLIP_OUTCOL_NUMBER_USED ] = size;
+ oa[ GAL_STATISTICS_CLIP_OUTCOL_MAD ] = sig1_mad0 ? NAN : spread;
+ oa[ GAL_STATISTICS_CLIP_OUTCOL_STD ] = sig1_mad0 ? spread : NAN;
+
+ /* Print the comments (if requested). */
if(!quiet)
- {
- printf("%-8s %-10s %-15s %-15s %-15s\n",
- "round", "number", "median", "mean", "STD");
- printf("%-8d %-10.0f %-15g %-15g %-15g\n",
- 0, oa[0], oa[1], oa[2], oa[3]);
- }
+ printf("%-5d %-10d %-12.5e %-12.5e\n", 1, 1,
+ oa[ GAL_STATISTICS_CLIP_OUTCOL_MEAN ], 0.0f);
break;
/* More than one element. */
default:
- /* Print the comments if requested. */
- if(!quiet)
- printf("%-8s %-10s %-15s %-15s %-15s\n",
- "round", "number", "median", "mean", "STD");
-
/* Do the clipping, but first initialize the values that will be
changed during the clipping: the start of the array and the
array's size. */
@@ -2424,23 +2645,26 @@ gal_statistics_sigma_clip(gal_data_t *input, float
multip, float param,
}
*/
- /* Find the mean, median and standard deviation. */
- meanstd=gal_statistics_mean_std(nbs);
- statistics_median_in_sorted_no_blank(nbs, median_i->array);
- median_d=gal_data_copy_to_new_type(median_i, GAL_TYPE_FLOAT64);
+ /* Find the center and disperson. */
+ statistics_median_in_sorted_no_blank(nbs, center_i->array);
+ if(sig1_mad0) spread_i=gal_statistics_std(nbs);
+ else statistics_mad_in_sorted_no_blank(nbs, center_i,
+ spread_i->array);
+ center_d=gal_data_copy_to_new_type(center_i, GAL_TYPE_FLOAT64);
+ spread_d=gal_data_copy_to_new_type(spread_i, GAL_TYPE_FLOAT64);
+ if(sig1_mad0) { gal_data_free(spread_i); spread_i=NULL; }
/* Put them in usable (with a type) pointers. */
- mean = meanstd->array;
- med = median_d->array;
- std = &((double *)(meanstd->array))[1];
+ center = ((double *)(center_d->array))[0];
+ spread = ((double *)(spread_d->array))[0];
/* If the user wanted to view the steps, show it to them. */
if(!quiet)
- printf("%-8zu %-10zu %-15g %-15g %-15g\n",
- num+1, size, *med, *mean, *std);
+ printf("%-5zu %-10zu %-12.5e %-12.5e\n", num+1, size, center,
+ spread);
- /* If we are to work by tolerance, then check if we should jump
- out of the loop. Normally, 'oldstd' should be larger than std,
+ /* If we are working by tolerance, check if we should jump out of
+ the loop. Normally, 'oldstd' should be larger than std,
because the possible outliers have been removed. If it is not,
it means that we have clipped too much and must stop anyway,
so we don't need an absolute value on the difference!
@@ -2450,10 +2674,10 @@ gal_statistics_sigma_clip(gal_data_t *input, float
multip, float param,
tolerance (because it will be infinity and thus lager than the
requested tolerance level value). */
if( bytolerance && num>0 )
- if( *std==0 || ((oldstd - *std) / *std) < param )
+ if( spread==0 || ((oldspread - spread) / spread) < param )
{
- if(*std==0) {oldmed=*med; oldstd=*std; oldmean=*mean;}
- gal_data_free(meanstd); gal_data_free(median_d);
+ if(spread==0) oldspread=spread;
+ gal_data_free(spread_d); gal_data_free(center_d);
break;
}
@@ -2462,16 +2686,16 @@ gal_statistics_sigma_clip(gal_data_t *input, float
multip, float param,
pointer and size of the array. */
switch(type)
{
- case GAL_TYPE_UINT8: SIGCLIP( uint8_t ); break;
- case GAL_TYPE_INT8: SIGCLIP( int8_t ); break;
- case GAL_TYPE_UINT16: SIGCLIP( uint16_t ); break;
- case GAL_TYPE_INT16: SIGCLIP( int16_t ); break;
- case GAL_TYPE_UINT32: SIGCLIP( uint32_t ); break;
- case GAL_TYPE_INT32: SIGCLIP( int32_t ); break;
- case GAL_TYPE_UINT64: SIGCLIP( uint64_t ); break;
- case GAL_TYPE_INT64: SIGCLIP( int64_t ); break;
- case GAL_TYPE_FLOAT32: SIGCLIP( float ); break;
- case GAL_TYPE_FLOAT64: SIGCLIP( double ); break;
+ case GAL_TYPE_UINT8: CLIPALL( uint8_t ); break;
+ case GAL_TYPE_INT8: CLIPALL( int8_t ); break;
+ case GAL_TYPE_UINT16: CLIPALL( uint16_t ); break;
+ case GAL_TYPE_INT16: CLIPALL( int16_t ); break;
+ case GAL_TYPE_UINT32: CLIPALL( uint32_t ); break;
+ case GAL_TYPE_INT32: CLIPALL( int32_t ); break;
+ case GAL_TYPE_UINT64: CLIPALL( uint64_t ); break;
+ case GAL_TYPE_INT64: CLIPALL( int64_t ); break;
+ case GAL_TYPE_FLOAT32: CLIPALL( float ); break;
+ case GAL_TYPE_FLOAT64: CLIPALL( double ); break;
default:
error(EXIT_FAILURE, 0, "%s: type code %d not recognized",
__func__, type);
@@ -2479,33 +2703,38 @@ gal_statistics_sigma_clip(gal_data_t *input, float
multip, float param,
/* Set the values from this round in the old elements, so the
next round can compare with, and return then if necessary. */
- oldmed = *med;
- oldstd = *std;
- oldmean = *mean;
+ oldspread = spread;
++num;
- /* Clean up. */
- gal_data_free(meanstd);
- gal_data_free(median_d);
+ /* Clean up: */
+ gal_data_free(spread_d);
+ gal_data_free(center_d);
}
/* If we were in tolerance mode and 'num' and 'maxnum' are equal (the
- loop didn't stop by tolerance), so the outputs should be NaN. */
+ loop didn't stop by tolerance), so the outputs should be NaN. Note
+ that they may have been filled in previous rounds compared to the
+ initialization (where they were all NaN). */
out->status=num;
+ oa[GAL_STATISTICS_CLIP_OUTCOL_NUMBER_CLIPS]=num;
if( size==0 || (bytolerance && num==maxnum) )
- oa[0] = oa[1] = oa[2] = oa[3] = NAN;
+ { for(i=0;i<GAL_STATISTICS_CLIP_OUT_SIZE;++i) oa[i]=NAN; }
else
{
- oa[0] = size;
- oa[1] = oldmed;
- oa[2] = oldmean;
- oa[3] = oldstd;
+ oa[ GAL_STATISTICS_CLIP_OUTCOL_MEDIAN ] = center;
+ oa[ GAL_STATISTICS_CLIP_OUTCOL_NUMBER_USED ] = size;
+ oa[ GAL_STATISTICS_CLIP_OUTCOL_MAD ] = sig1_mad0 ? NAN : spread;
+ oa[ GAL_STATISTICS_CLIP_OUTCOL_STD ] = sig1_mad0 ? spread : NAN;
}
}
+ /* Measure and report the remaining elements if requested. */
+ if(extrastats) statistics_clip_stats_extra(nbs, oa, extrastats);
+
/* Clean up and return. */
nbs->array=nbs_array;
- gal_data_free(median_i);
+ gal_data_free(center_i);
+ gal_data_free(spread_i);
if(nbs!=input) gal_data_free(nbs);
return out;
}
@@ -2514,6 +2743,30 @@ gal_statistics_sigma_clip(gal_data_t *input, float
multip, float param,
+gal_data_t *
+gal_statistics_clip_sigma(gal_data_t *input, float multip, float param,
+ uint8_t extrastats, int inplace, int quiet)
+{
+ return statistics_clip(input, multip, param, extrastats,
+ inplace, quiet, 1);
+}
+
+
+
+
+
+gal_data_t *
+gal_statistics_clip_mad(gal_data_t *input, float multip, float param,
+ uint8_t extrastats, int inplace, int quiet)
+{
+ return statistics_clip(input, multip, param, extrastats,
+ inplace, quiet, 0);
+}
+
+
+
+
+
/* Find the first outlier in a distribution. */
#define OUTLIER_BYTYPE(IT) { \
IT *arr=nbs->array; \
@@ -2528,19 +2781,25 @@ gal_statistics_sigma_clip(gal_data_t *input, float
multip, float param,
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); \
+ sclip=gal_statistics_clip_mad(dist, sigclip_multip, \
+ sigclip_param, clipflags, 0, 1); \
sarr=sclip->array; \
\
/* For a check. */ \
if(quiet==0) \
printf("%f [%zu]: %f (%f, %f) %f\n", (float)(arr[i]), i, \
- (float)(arr[i]-arr[i-1]), sarr[1], sarr[3], \
- (((double)(arr[i]-arr[i-1])) - sarr[1])/sarr[3]); \
+ (float)(arr[i]-arr[i-1]), \
+ sarr[GAL_STATISTICS_CLIP_OUTCOL_NUMBER_USED], \
+ sarr[GAL_STATISTICS_CLIP_OUTCOL_STD], \
+ (((double)(arr[i]-arr[i-1])) \
+ - sarr[GAL_STATISTICS_CLIP_OUTCOL_MEDIAN]) \
+ /sarr[GAL_STATISTICS_CLIP_OUTCOL_STD]); \
\
/* 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] ) \
+ if( (((double)(arr[i]-arr[i-1])) \
+ - sarr[GAL_STATISTICS_CLIP_OUTCOL_MEDIAN]) \
+ > sigma*sarr[GAL_STATISTICS_CLIP_OUTCOL_STD] ) \
{ \
/* Allocate the output dataset. */ \
out=gal_data_alloc(NULL, input->type, 1, &one, NULL, 0, -1, \
@@ -2566,6 +2825,7 @@ gal_statistics_outlier_bydistance(int pos1_neg0,
gal_data_t *input,
double *darr;
size_t i, j, one=1, wtakeone;
gal_data_t *dist, *sclip, *nbs, *out=NULL;
+ uint8_t clipflags=GAL_STATISTICS_CLIP_OUTCOL_STD;
/* Remove all blanks and sort the dataset. */
nbs=gal_statistics_no_blank_sorted(input, inplace);
@@ -2645,19 +2905,22 @@ gal_statistics_outlier_bydistance(int pos1_neg0,
gal_data_t *input,
/* Sigma-clipped median and std for a check. */ \
prev->flag=0; \
prev->size=prev->dsize[0]=numprev; \
- sclip=gal_statistics_sigma_clip(prev, sigclip_multip, \
- sigclip_param, 1, 1); \
+ sclip=gal_statistics_clip_mad(prev, sigclip_multip, \
+ sigclip_param, clipflags, \
+ 1, 1); \
\
sarr=sclip->array; \
- check = (diff - sarr[1]) / sarr[3]; \
+ check = ( (diff - sarr[GAL_STATISTICS_CLIP_OUTCOL_MEDIAN]) \
+ / sarr[GAL_STATISTICS_CLIP_OUTCOL_STD] ); \
\
/* If requested, print the values. */ \
if(!quiet) printf("%-6zu%-15g%-15g%-15g (%g,%g)\n", p-a, \
- (float)(*p), (float)diff, check, sarr[1], \
- sarr[3]); \
+ (float)(*p), (float)diff, check, \
+ sarr[GAL_STATISTICS_CLIP_OUTCOL_MEDIAN], \
+ sarr[GAL_STATISTICS_CLIP_OUTCOL_STD]); \
\
- /* When values are equal, std will be roughly zero. */ \
- if(sarr[3]>1e-6 && check>thresh) \
+ /* When values are equal, std will be roughly zero */ \
+ if(sarr[GAL_STATISTICS_CLIP_OUTCOL_STD]>1e-6 && check>thresh) \
{ \
if(flatind==GAL_BLANK_SIZE_T) \
{ \
@@ -2683,7 +2946,7 @@ gal_statistics_outlier_bydistance(int pos1_neg0,
gal_data_t *input,
} \
} \
while(++p<pp); \
- if(counter+1!=numcontig) flatind=GAL_BLANK_SIZE_T; \
+ if(counter+1!=numcontig) flatind=GAL_BLANK_SIZE_T; \
}
gal_data_t *
@@ -2695,6 +2958,7 @@ gal_statistics_outlier_flat_cfp(gal_data_t *input, size_t
numprev,
float *sarr;
double check;
gal_data_t *nbs, *prev, *out=NULL, *sclip;
+ uint8_t clipflags=GAL_STATISTICS_CLIP_OUTCOL_STD;
size_t d=2, counter=0, one=1, flatind=GAL_BLANK_SIZE_T;
/* Sanity checks. */
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [gnuastro-commits] master 8f3b4804: Library (statistics.h): MAD and MAD-clipping functions added,
Mohammad Akhlaghi <=