[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master ce3ebeb: New sigma-clipped mean and median fil
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master ce3ebeb: New sigma-clipped mean and median filters for Arithmetic |
|
Date: |
Thu, 22 Mar 2018 10:22:20 -0400 (EDT) |
branch: master
commit ce3ebeb284bdea1311c2bb3a3633286a7d478a96
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
New sigma-clipped mean and median filters for Arithmetic
The Arithmetic program now has two new filters: `filter-sigclip-mean' and
`filter-sigclip-median'. They can be used to remove strong outliers in the
filter before estimating the mean or median.
Arithmetic's filtering operators were also measuring the filter on blank
elements (in the input). This is not desired: a blank element in the input
should also be blank in the output. So a check was added to each pixel
before the filter is measured for it.
While trying to ignore blank pixels, I noticed that we don't have a simple
function to check a single value to see if its blank or not. So a new
function was defined for this job: `gal_blank_is'.
In the process, I also noticed that when using `gal_blank_remove' or
`gal_statistics_no_blank_sorted', we weren't checking for the case that all
the input pixels are blank. All occurances of these two functions now check
for this.
---
NEWS | 7 +-
bin/arithmetic/arithmetic.c | 364 +++++++++++++++++++++++++++-----------------
bin/arithmetic/arithmetic.h | 4 +-
bin/noisechisel/threshold.c | 7 +
bin/statistics/ui.c | 6 +
doc/gnuastro.texi | 71 ++++++++-
lib/blank.c | 51 +++++++
lib/gnuastro/blank.h | 3 +
lib/statistics.c | 335 ++++++++++++++++++++++------------------
9 files changed, 551 insertions(+), 297 deletions(-)
diff --git a/NEWS b/NEWS
index 2cfe321..9e1c36d 100644
--- a/NEWS
+++ b/NEWS
@@ -5,11 +5,15 @@ GNU Astronomy Utilities NEWS -*-
outline -*-
** New features
+ Arithmetic:
+ - filter-sigclip-mean: sigma-clipped, mean filter operator.
+ - filter-sigclip-median: sigma-clipped, median filter operator.
+
ConvertType:
- TIFF images can also be used as input.
MakeCatalog:
- --variance: input STD image is variance (so take square root before use).
+ --variance: input standard deviation image is actually variance.
--brightnesserr: error in estimating the brightness.
--mean: calculate the mean pixel value within an object or clump.
--median: calculate the median pixel value within an object or clump.
@@ -24,6 +28,7 @@ GNU Astronomy Utilities NEWS -*-
outline -*-
Libraries:
gal_array_read: read array from any of known formats (FITS,TIFF,JPEG,...).
gal_array_read_to_type: similar to `gal_array_read', but to given type.
+ gal_blank_is: check to see if argument is blank in its type or not.
gal_eps_name_is_eps: Returns 1 if given filename is EPS.
gal_eps_suffix_is_eps: Returns 1 if given suffix is EPS.
gal_eps_to_pt: Converts dataset size to PostScript points.
diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index 9bc4786..7b2b067 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -131,14 +131,16 @@ pop_number_of_operands(struct arithmeticparams *p,
gal_data_t *data,
struct arithmetic_filter_p
{
- int operator; /* The type of filtering. */
- size_t *fsize; /* Filter size. */
- size_t *hpfsize; /* Positive Half-filter size. */
- size_t *hnfsize; /* Negative Half-filter size. */
- gal_data_t *input; /* Input dataset. */
- gal_data_t *out; /* Output dataset. */
-
- int hasblank; /* If the dataset has blank values.*/
+ int operator; /* The type of filtering. */
+ size_t *fsize; /* Filter size. */
+ size_t *hpfsize; /* Positive Half-filter size. */
+ size_t *hnfsize; /* Negative Half-filter size. */
+ float sclip_multip; /* Sigma multiple in sigma-clipping. */
+ float sclip_param; /* Termination critera in sigma-cliping. */
+ gal_data_t *input; /* Input dataset. */
+ gal_data_t *out; /* Output dataset. */
+
+ int hasblank; /* If the dataset has blank values. */
};
@@ -153,9 +155,9 @@ arithmetic_filter(void *in_prm)
struct arithmetic_filter_p *afp=(struct arithmetic_filter_p *)tprm->params;
gal_data_t *input=afp->input;
- size_t ind,index;
- int out_type_checked=0;
- gal_data_t *result=NULL;
+ size_t sind=-1;
+ size_t ind, index, one=1;
+ gal_data_t *sigclip, *result=NULL;
size_t *hpfsize=afp->hpfsize, *hnfsize=afp->hnfsize;
size_t *tsize, *dsize=input->dsize, *fsize=afp->fsize;
size_t i, j, coord[ARITHMETIC_FILTER_DIM], ndim=input->ndim;
@@ -172,85 +174,128 @@ arithmetic_filter(void *in_prm)
/* Go over all the pixels that were assigned to this thread. */
for(i=0; tprm->indexs[i] != GAL_BLANK_SIZE_T; ++i)
{
- /* Get the coordinate of the pixel. */
+ /* For easy reading, put the index in `ind'. */
ind=tprm->indexs[i];
- gal_dimension_index_to_coord(ind, ndim, dsize, coord);
- /* See which dimensions need trimming. */
- tile->size=1;
- for(j=0;j<ndim;++j)
+ /* If we are on a blank element, then just set the output to blank
+ also. */
+ if( afp->hasblank
+ && gal_blank_is(gal_data_ptr_increment(input->array, ind,
+ input->type), input->type) )
+ gal_blank_write(gal_data_ptr_increment(afp->out->array, ind,
+ afp->out->type),
+ afp->out->type);
+ else
{
- /* Estimate the coordinate of the filter's starting point. Note
- that we are dealing with size_t (unsigned int) type here, so
- there are no negatives. A negative result will produce an
- extremely large number, so instead of checking for negative,
- we can just see if the result of a subtraction is less than
- the width of the input. */
- if( (coord[j] - hnfsize[j] > dsize[j])
- || (coord[j] + hpfsize[j] >= dsize[j]) )
+ /* Get the coordinate of the pixel. */
+ gal_dimension_index_to_coord(ind, ndim, dsize, coord);
+
+ /* See which dimensions need trimming. */
+ tile->size=1;
+ for(j=0;j<ndim;++j)
{
- start[j] = ( (coord[j] - hnfsize[j] > dsize[j])
- ? 0 : coord[j] - hnfsize[j] );
- end[j] = ( (coord[j] + hpfsize[j] >= dsize[j])
- ? dsize[j]
- : coord[j] + hpfsize[j] + 1);
- tsize[j] = end[j] - start[j];
+ /* Estimate the coordinate of the filter's starting
+ point. Note that we are dealing with size_t (unsigned int)
+ type here, so there are no negatives. A negative result
+ will produce an extremely large number, so instead of
+ checking for negative, we can just see if the result of a
+ subtraction is less than the width of the input. */
+ if( (coord[j] - hnfsize[j] > dsize[j])
+ || (coord[j] + hpfsize[j] >= dsize[j]) )
+ {
+ start[j] = ( (coord[j] - hnfsize[j] > dsize[j])
+ ? 0 : coord[j] - hnfsize[j] );
+ end[j] = ( (coord[j] + hpfsize[j] >= dsize[j])
+ ? dsize[j]
+ : coord[j] + hpfsize[j] + 1);
+ tsize[j] = end[j] - start[j];
+ }
+ else /* NOT on the edge (given requested filter width). */
+ {
+ tsize[j] = fsize[j];
+ start[j] = coord[j] - hnfsize[j];
+ }
+ tile->size *= tsize[j];
}
- else /* We are NOT on the edge (given requested filter width). */
+
+ /* For a test.
+ printf("coord: %zu, %zu\n", coord[1]+1, coord[0]+1);
+ printf("\tstart: %zu, %zu\n", start[1]+1, start[0]+1);
+ printf("\ttsize: %zu, %zu\n", tsize[1], tsize[0]);
+ */
+
+ /* Set the tile's starting pointer. */
+ index=gal_dimension_coord_to_index(ndim, dsize, start);
+ tile->array=gal_data_ptr_increment(input->array, index,
+ input->type);
+
+ /* Do the necessary calculation. */
+ switch(afp->operator)
{
- tsize[j] = fsize[j];
- start[j] = coord[j] - hnfsize[j];
- }
- tile->size *= tsize[j];
- }
+ case ARITHMETIC_OP_FILTER_MEDIAN:
+ result=gal_statistics_median(tile, 0);
+ break;
- /* Set the tile's starting pointer. */
- index=gal_dimension_coord_to_index(ndim, dsize, start);
- tile->array=gal_data_ptr_increment(input->array, index, input->type);
- /* Do the necessary calculation. */
- switch(afp->operator)
- {
- case ARITHMETIC_OP_FILTER_MEDIAN:
- result=gal_statistics_median(tile, 0);
- break;
+ case ARITHMETIC_OP_FILTER_MEAN:
+ result=gal_statistics_mean(tile);
+ break;
- case ARITHMETIC_OP_FILTER_MEAN:
- result=gal_statistics_mean(tile);
- break;
- default:
- error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
- "fix the problem. `afp->operator' code %d is not "
- "recognized", PACKAGE_BUGREPORT, __func__, afp->operator);
- }
+ 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);
+
+ /* 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;
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at "
+ "%s to fix the problem. The `afp->operator' value "
+ "%d is not recognized as sigma-clipped median or "
+ "mean", __func__, PACKAGE_BUGREPORT, afp->operator);
+ }
+
+ /* Allocate the output and write the value into it. */
+ result=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, 1, &one, NULL,
+ 0, -1, NULL, NULL, NULL);
+ ((float *)(result->array))[0] =
+ ((float *)(sigclip->array))[sind];
+
+ /* Clean up. */
+ gal_data_free(sigclip);
+ break;
+
+
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s "
+ "to fix the problem. `afp->operator' code %d is not "
+ "recognized", PACKAGE_BUGREPORT, __func__,
+ afp->operator);
+ }
+
+ /* Make sure the output array type and result's type are the
+ same. */
+ if(result->type!=afp->out->type)
+ result=gal_data_copy_to_new_type_free(result, afp->out->type);
- /* Make sure the output array type and result's type are the same. We
- only need to do this once, but we'll suffice to once for each
- thread for simplicify of the code, it is too negligible to have
- any real affect. */
- if( out_type_checked==0)
- {
- if(result->type!=afp->out->type )
- error(EXIT_FAILURE, 0, "%s: a bug! please contact us at %s so "
- "we can address the problem. The tyes of `result' and "
- "`out' aren't the same, they are respectively: `%s' and "
- "`%s'", __func__, PACKAGE_BUGREPORT,
- gal_type_name(result->type, 1),
- gal_type_name(afp->out->type, 1));
- out_type_checked=1;
- }
- /* Copy the result into the output array. */
- memcpy(gal_data_ptr_increment(afp->out->array, ind, afp->out->type),
- result->array, gal_type_sizeof(afp->out->type));
+ /* Copy the result into the output array. */
+ memcpy(gal_data_ptr_increment(afp->out->array, ind,
+ afp->out->type),
+ result->array, gal_type_sizeof(afp->out->type));
- /* Clean up. */
- gal_data_free(result);
+ /* Clean up for this pixel. */
+ gal_data_free(result);
+ }
}
- /* Clean up. */
+ /* Clean up for this thread. */
tile->array=NULL;
tile->block=NULL;
gal_data_free(tile);
@@ -268,12 +313,14 @@ arithmetic_filter(void *in_prm)
static void
wrapper_for_filter(struct arithmeticparams *p, char *token, int operator)
{
- size_t i=0, ndim, one=1;
int type=GAL_TYPE_INVALID;
+ size_t i=0, ndim, nparams, one=1;
struct arithmetic_filter_p afp={0};
size_t fsize[ARITHMETIC_FILTER_DIM];
- gal_data_t *tmp, *tmp2, *zero, *comp, *fsize_list=NULL;
+ gal_data_t *tmp, *tmp2, *zero, *comp, *params_list=NULL;
size_t hnfsize[ARITHMETIC_FILTER_DIM], hpfsize[ARITHMETIC_FILTER_DIM];
+ int issigclip=(operator==ARITHMETIC_OP_FILTER_SIGCLIP_MEAN
+ || operator==ARITHMETIC_OP_FILTER_SIGCLIP_MEDIAN);
/* Get the input's number of dimensions. */
@@ -297,21 +344,40 @@ wrapper_for_filter(struct arithmeticparams *p, char
*token, int operator)
NULL, NULL);
- /* Based on the dimensions of the popped operand, pop the necessary
- number of operands. */
- for(i=0;i<ndim;++i)
- gal_list_data_add(&fsize_list, operands_pop(p, token));
+ /* Based on the first popped operand's dimensions and the operator, of
+ pop the necessary number of operands. */
+ nparams = ndim + (issigclip ? 2 : 0 );
+ for(i=0;i<nparams;++i)
+ gal_list_data_add(¶ms_list, operands_pop(p, token));
- /* Make sure the filter size only has single values. */
+ /* Make sure the parameters only have single values. */
i=0;
- for(tmp=fsize_list; tmp!=NULL; tmp=tmp->next)
+ for(tmp=params_list; tmp!=NULL; tmp=tmp->next)
{
++i;
if(tmp->size!=1)
- error(EXIT_FAILURE, 0, "the filter length values given to the "
- "filter operators can only be numbers. Value number %zu has "
- "%zu elements, so its an array", ndim-i-1, tmp->size);
+ error(EXIT_FAILURE, 0, "the parameters given to the filtering "
+ "operators can only be numbers. Value number %zu has %zu "
+ "elements, so its an array", i, tmp->size);
+ }
+
+
+ /* If this is a sigma-clipping filter, the top two operands are the
+ sigma-clipping parameters. */
+ if(issigclip)
+ {
+ /* Read the sigma-clipping multiple (first element in the list). */
+ tmp=gal_list_data_pop(¶ms_list);
+ tmp=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT32);
+ afp.sclip_multip=*(float *)(tmp->array);
+ gal_data_free(tmp);
+
+ /* Read the sigma-clipping termination parameter. */
+ tmp=gal_list_data_pop(¶ms_list);
+ tmp=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT32);
+ afp.sclip_param=*(float *)(tmp->array);
+ gal_data_free(tmp);
}
@@ -324,7 +390,7 @@ wrapper_for_filter(struct arithmeticparams *p, char *token,
int operator)
values must be written in the inverse order since the user gives
dimensions with the FITS standard. */
i=ndim-1;
- for(tmp=fsize_list; tmp!=NULL; tmp=tmp->next)
+ for(tmp=params_list; tmp!=NULL; tmp=tmp->next)
{
/* Make sure the user has given an integer type. */
if(tmp->type==GAL_TYPE_FLOAT32 || tmp->type==GAL_TYPE_FLOAT64)
@@ -370,18 +436,26 @@ wrapper_for_filter(struct arithmeticparams *p, char
*token, int operator)
{ hnfsize[i]=fsize[i]/2; hpfsize[i]=fsize[i]/2-1; }
}
+ /* For a test.
+ printf("fsize: %zu, %zu\n", fsize[0], fsize[1]);
+ printf("hnfsize: %zu, %zu\n", hnfsize[0], hnfsize[1]);
+ printf("hpfsize: %zu, %zu\n", hpfsize[0], hpfsize[1]);
+ */
/* See if the input has blank pixels. */
afp.hasblank=gal_blank_present(afp.input, 1);
+
/* Set the type of the output dataset. */
switch(operator)
{
case ARITHMETIC_OP_FILTER_MEDIAN:
+ case ARITHMETIC_OP_FILTER_SIGCLIP_MEDIAN:
type=afp.input->type;
break;
case ARITHMETIC_OP_FILTER_MEAN:
+ case ARITHMETIC_OP_FILTER_SIGCLIP_MEAN:
type=GAL_TYPE_FLOAT64;
break;
@@ -391,12 +465,14 @@ wrapper_for_filter(struct arithmeticparams *p, char
*token, int operator)
PACKAGE_BUGREPORT, __func__, operator);
}
+
/* Allocate the output dataset. Note that filtering doesn't change
the units of the dataset. */
afp.out=gal_data_alloc(NULL, type, ndim, afp.input->dsize,
afp.input->wcs, 0, afp.input->minmapsize,
NULL, afp.input->unit, NULL);
+
/* Spin off threads for each pixel. */
gal_threads_spin_off(arithmetic_filter, &afp, afp.input->size,
p->cp.numthreads);
@@ -409,7 +485,7 @@ wrapper_for_filter(struct arithmeticparams *p, char *token,
int operator)
/* Clean up and add the output on top of the stack */
gal_data_free(afp.input);
- gal_list_data_free(fsize_list);
+ gal_list_data_free(params_list);
}
@@ -472,123 +548,127 @@ reversepolish(struct arithmeticparams *p)
/* Order is the same as in the manual. */
/* Simple arithmetic operators. */
if (!strcmp(token->v, "+" ))
- { op=GAL_ARITHMETIC_OP_PLUS; nop=2; }
+ { op=GAL_ARITHMETIC_OP_PLUS; nop=2; }
else if (!strcmp(token->v, "-" ))
- { op=GAL_ARITHMETIC_OP_MINUS; nop=2; }
+ { op=GAL_ARITHMETIC_OP_MINUS; nop=2; }
else if (!strcmp(token->v, "x" ))
- { op=GAL_ARITHMETIC_OP_MULTIPLY; nop=2; }
+ { op=GAL_ARITHMETIC_OP_MULTIPLY; nop=2; }
else if (!strcmp(token->v, "/" ))
- { op=GAL_ARITHMETIC_OP_DIVIDE; nop=2; }
+ { op=GAL_ARITHMETIC_OP_DIVIDE; nop=2; }
else if (!strcmp(token->v, "%" ))
- { op=GAL_ARITHMETIC_OP_MODULO; nop=2; }
+ { op=GAL_ARITHMETIC_OP_MODULO; nop=2; }
/* Mathematical Operators. */
else if (!strcmp(token->v, "abs"))
- { op=GAL_ARITHMETIC_OP_ABS; nop=1; }
+ { op=GAL_ARITHMETIC_OP_ABS; nop=1; }
else if (!strcmp(token->v, "pow"))
- { op=GAL_ARITHMETIC_OP_POW; nop=2; }
+ { op=GAL_ARITHMETIC_OP_POW; nop=2; }
else if (!strcmp(token->v, "sqrt"))
- { op=GAL_ARITHMETIC_OP_SQRT; nop=1; }
+ { op=GAL_ARITHMETIC_OP_SQRT; nop=1; }
else if (!strcmp(token->v, "log"))
- { op=GAL_ARITHMETIC_OP_LOG; nop=1; }
+ { op=GAL_ARITHMETIC_OP_LOG; nop=1; }
else if (!strcmp(token->v, "log10"))
- { op=GAL_ARITHMETIC_OP_LOG10; nop=1; }
+ { op=GAL_ARITHMETIC_OP_LOG10; nop=1; }
/* Statistical/higher-level operators. */
else if (!strcmp(token->v, "minvalue"))
- { op=GAL_ARITHMETIC_OP_MINVAL; nop=1; }
+ { op=GAL_ARITHMETIC_OP_MINVAL; nop=1; }
else if (!strcmp(token->v, "maxvalue"))
- { op=GAL_ARITHMETIC_OP_MAXVAL; nop=1; }
+ { op=GAL_ARITHMETIC_OP_MAXVAL; nop=1; }
else if (!strcmp(token->v, "numvalue"))
- { op=GAL_ARITHMETIC_OP_NUMVAL; nop=1; }
+ { op=GAL_ARITHMETIC_OP_NUMVAL; nop=1; }
else if (!strcmp(token->v, "sumvalue"))
- { op=GAL_ARITHMETIC_OP_SUMVAL; nop=1; }
+ { op=GAL_ARITHMETIC_OP_SUMVAL; nop=1; }
else if (!strcmp(token->v, "meanvalue"))
- { op=GAL_ARITHMETIC_OP_MEANVAL; nop=1; }
+ { op=GAL_ARITHMETIC_OP_MEANVAL; nop=1; }
else if (!strcmp(token->v, "stdvalue"))
- { op=GAL_ARITHMETIC_OP_STDVAL; nop=1; }
+ { op=GAL_ARITHMETIC_OP_STDVAL; nop=1; }
else if (!strcmp(token->v, "medianvalue"))
- { op=GAL_ARITHMETIC_OP_MEDIANVAL; nop=1; }
+ { op=GAL_ARITHMETIC_OP_MEDIANVAL; nop=1; }
else if (!strcmp(token->v, "min"))
- { op=GAL_ARITHMETIC_OP_MIN; nop=-1; }
+ { op=GAL_ARITHMETIC_OP_MIN; nop=-1; }
else if (!strcmp(token->v, "max"))
- { op=GAL_ARITHMETIC_OP_MAX; nop=-1; }
+ { op=GAL_ARITHMETIC_OP_MAX; nop=-1; }
else if (!strcmp(token->v, "num"))
- { op=GAL_ARITHMETIC_OP_NUM; nop=-1; }
+ { op=GAL_ARITHMETIC_OP_NUM; nop=-1; }
else if (!strcmp(token->v, "sum"))
- { op=GAL_ARITHMETIC_OP_SUM; nop=-1; }
+ { op=GAL_ARITHMETIC_OP_SUM; nop=-1; }
else if (!strcmp(token->v, "mean"))
- { op=GAL_ARITHMETIC_OP_MEAN; nop=-1; }
+ { op=GAL_ARITHMETIC_OP_MEAN; nop=-1; }
else if (!strcmp(token->v, "std"))
- { op=GAL_ARITHMETIC_OP_STD; nop=-1; }
+ { op=GAL_ARITHMETIC_OP_STD; nop=-1; }
else if (!strcmp(token->v, "median"))
- { op=GAL_ARITHMETIC_OP_MEDIAN; nop=-1; }
+ { op=GAL_ARITHMETIC_OP_MEDIAN; nop=-1; }
/* Conditional operators. */
else if (!strcmp(token->v, "lt" ))
- { op=GAL_ARITHMETIC_OP_LT; nop=2; }
+ { op=GAL_ARITHMETIC_OP_LT; nop=2; }
else if (!strcmp(token->v, "le"))
- { op=GAL_ARITHMETIC_OP_LE; nop=2; }
+ { op=GAL_ARITHMETIC_OP_LE; nop=2; }
else if (!strcmp(token->v, "gt" ))
- { op=GAL_ARITHMETIC_OP_GT; nop=2; }
+ { op=GAL_ARITHMETIC_OP_GT; nop=2; }
else if (!strcmp(token->v, "ge"))
- { op=GAL_ARITHMETIC_OP_GE; nop=2; }
+ { op=GAL_ARITHMETIC_OP_GE; nop=2; }
else if (!strcmp(token->v, "eq"))
- { op=GAL_ARITHMETIC_OP_EQ; nop=2; }
+ { op=GAL_ARITHMETIC_OP_EQ; nop=2; }
else if (!strcmp(token->v, "ne"))
- { op=GAL_ARITHMETIC_OP_NE; nop=2; }
+ { op=GAL_ARITHMETIC_OP_NE; nop=2; }
else if (!strcmp(token->v, "and"))
- { op=GAL_ARITHMETIC_OP_AND; nop=2; }
+ { op=GAL_ARITHMETIC_OP_AND; nop=2; }
else if (!strcmp(token->v, "or"))
- { op=GAL_ARITHMETIC_OP_OR; nop=2; }
+ { op=GAL_ARITHMETIC_OP_OR; nop=2; }
else if (!strcmp(token->v, "not"))
- { op=GAL_ARITHMETIC_OP_NOT; nop=1; }
+ { op=GAL_ARITHMETIC_OP_NOT; nop=1; }
else if (!strcmp(token->v, "isblank"))
- { op=GAL_ARITHMETIC_OP_ISBLANK; nop=1; }
+ { op=GAL_ARITHMETIC_OP_ISBLANK; nop=1; }
else if (!strcmp(token->v, "where"))
- { op=GAL_ARITHMETIC_OP_WHERE; nop=3; }
+ { op=GAL_ARITHMETIC_OP_WHERE; nop=3; }
/* Bitwise operators. */
else if (!strcmp(token->v, "bitand"))
- { op=GAL_ARITHMETIC_OP_BITAND; nop=2; }
+ { op=GAL_ARITHMETIC_OP_BITAND; nop=2; }
else if (!strcmp(token->v, "bitor"))
- { op=GAL_ARITHMETIC_OP_BITOR; nop=2; }
+ { op=GAL_ARITHMETIC_OP_BITOR; nop=2; }
else if (!strcmp(token->v, "bitxor"))
- { op=GAL_ARITHMETIC_OP_BITXOR; nop=2; }
+ { op=GAL_ARITHMETIC_OP_BITXOR; nop=2; }
else if (!strcmp(token->v, "lshift"))
- { op=GAL_ARITHMETIC_OP_BITLSH; nop=2; }
+ { op=GAL_ARITHMETIC_OP_BITLSH; nop=2; }
else if (!strcmp(token->v, "rshift"))
- { op=GAL_ARITHMETIC_OP_BITRSH; nop=2; }
+ { op=GAL_ARITHMETIC_OP_BITRSH; nop=2; }
else if (!strcmp(token->v, "bitnot"))
- { op=GAL_ARITHMETIC_OP_BITNOT; nop=1; }
+ { op=GAL_ARITHMETIC_OP_BITNOT; nop=1; }
/* Type conversion. */
else if (!strcmp(token->v, "uint8"))
- { op=GAL_ARITHMETIC_OP_TO_UINT8; nop=1; }
+ { op=GAL_ARITHMETIC_OP_TO_UINT8; nop=1; }
else if (!strcmp(token->v, "int8"))
- { op=GAL_ARITHMETIC_OP_TO_INT8; nop=1; }
+ { op=GAL_ARITHMETIC_OP_TO_INT8; nop=1; }
else if (!strcmp(token->v, "uint16"))
- { op=GAL_ARITHMETIC_OP_TO_UINT16; nop=1; }
+ { op=GAL_ARITHMETIC_OP_TO_UINT16; nop=1; }
else if (!strcmp(token->v, "int16"))
- { op=GAL_ARITHMETIC_OP_TO_INT16; nop=1; }
+ { op=GAL_ARITHMETIC_OP_TO_INT16; nop=1; }
else if (!strcmp(token->v, "uint32"))
- { op=GAL_ARITHMETIC_OP_TO_UINT32; nop=1; }
+ { op=GAL_ARITHMETIC_OP_TO_UINT32; nop=1; }
else if (!strcmp(token->v, "int32"))
- { op=GAL_ARITHMETIC_OP_TO_INT32; nop=1; }
+ { op=GAL_ARITHMETIC_OP_TO_INT32; nop=1; }
else if (!strcmp(token->v, "uint64"))
- { op=GAL_ARITHMETIC_OP_TO_UINT64; nop=1; }
+ { op=GAL_ARITHMETIC_OP_TO_UINT64; nop=1; }
else if (!strcmp(token->v, "int64"))
- { op=GAL_ARITHMETIC_OP_TO_INT64; nop=1; }
+ { op=GAL_ARITHMETIC_OP_TO_INT64; nop=1; }
else if (!strcmp(token->v, "float32"))
- { op=GAL_ARITHMETIC_OP_TO_FLOAT32; nop=1; }
+ { op=GAL_ARITHMETIC_OP_TO_FLOAT32; nop=1; }
else if (!strcmp(token->v, "float64"))
- { op=GAL_ARITHMETIC_OP_TO_FLOAT64; nop=1; }
+ { op=GAL_ARITHMETIC_OP_TO_FLOAT64; nop=1; }
/* Filters. */
- else if (!strcmp(token->v, "filter-median"))
- { op=ARITHMETIC_OP_FILTER_MEDIAN; nop=0; }
else if (!strcmp(token->v, "filter-mean"))
- { op=ARITHMETIC_OP_FILTER_MEAN; nop=0; }
+ { op=ARITHMETIC_OP_FILTER_MEAN; nop=0; }
+ else if (!strcmp(token->v, "filter-median"))
+ { op=ARITHMETIC_OP_FILTER_MEDIAN; nop=0; }
+ else if (!strcmp(token->v, "filter-sigclip-mean"))
+ { op=ARITHMETIC_OP_FILTER_SIGCLIP_MEAN; nop=0; }
+ else if (!strcmp(token->v, "filter-sigclip-median"))
+ { op=ARITHMETIC_OP_FILTER_SIGCLIP_MEDIAN; nop=0; }
/* Finished checks with known operators */
@@ -658,6 +738,8 @@ reversepolish(struct arithmeticparams *p)
{
case ARITHMETIC_OP_FILTER_MEAN:
case ARITHMETIC_OP_FILTER_MEDIAN:
+ case ARITHMETIC_OP_FILTER_SIGCLIP_MEAN:
+ case ARITHMETIC_OP_FILTER_SIGCLIP_MEDIAN:
wrapper_for_filter(p, token->v, op);
break;
diff --git a/bin/arithmetic/arithmetic.h b/bin/arithmetic/arithmetic.h
index fc5b65c..a9ab1fa 100644
--- a/bin/arithmetic/arithmetic.h
+++ b/bin/arithmetic/arithmetic.h
@@ -31,8 +31,10 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
library. */
enum arithmetic_prog_operators
{
- ARITHMETIC_OP_FILTER_MEDIAN=GAL_ARITHMETIC_OP_LAST_CODE,
+ ARITHMETIC_OP_FILTER_MEDIAN = GAL_ARITHMETIC_OP_LAST_CODE,
ARITHMETIC_OP_FILTER_MEAN,
+ ARITHMETIC_OP_FILTER_SIGCLIP_MEAN,
+ ARITHMETIC_OP_FILTER_SIGCLIP_MEDIAN,
};
diff --git a/bin/noisechisel/threshold.c b/bin/noisechisel/threshold.c
index 098f69e..09f66a9 100644
--- a/bin/noisechisel/threshold.c
+++ b/bin/noisechisel/threshold.c
@@ -196,10 +196,17 @@ threshold_write_sn_table(struct noisechiselparams *p,
gal_data_t *insn,
set of blank elements, but checking on the integer array is faster. */
if( gal_blank_present(inind, 1) )
{
+ /* Remove blank elements. */
ind=gal_data_copy(inind);
sn=gal_data_copy(insn);
gal_blank_remove(ind);
gal_blank_remove(sn);
+
+ /* A small sanity check. */
+ if(ind->size==0 || sn->size==0)
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix "
+ "the problem. For some reason, all the elements in `ind' or "
+ "`sn' are blank", __func__, PACKAGE_BUGREPORT);
}
else
{
diff --git a/bin/statistics/ui.c b/bin/statistics/ui.c
index 36c8810..57a2aaf 100644
--- a/bin/statistics/ui.c
+++ b/bin/statistics/ui.c
@@ -852,6 +852,12 @@ ui_preparations(struct statisticsparams *p)
{
/* Only keep the elements we want. */
gal_blank_remove(p->input);
+
+ /* Make sure there actually are any (non-blank) elements left. */
+ if(p->input->size==0)
+ error(EXIT_FAILURE, 0, "%s: all elements are blank",
+ gal_fits_name_save_as_string(p->inputname, cp->hdu));
+
p->input->flag &= ~GAL_DATA_FLAG_HASBLANK ;
p->input->flag |= GAL_DATA_FLAG_BLANK_CH ;
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index d674010..df88856 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -10017,7 +10017,7 @@ viewed in ds9) and 4 pixels along the second FITS
dimension (vertical).
Each pixel will be placed in the center of the box that the mean is
calculated on. If the given width along a dimension is even, then the
center is assumed to be between the pixels (not in the center of a
-pixel). When the pixel is close to the center, the pixels of the box that
+pixel). When the pixel is close to the edge, the pixels of the box that
fall outside the image are ignored. Therefore, on the edge, less points
will be used in calculating the mean.
@@ -10036,6 +10036,39 @@ The median is less susceptible to outliers compared to
the mean. As a
result, after median filtering, the pixel values will be more discontinuous
than mean filtering.
address@hidden filter-sigclip-mean
+Apply a @mymath{\sigma}-clipped mean filtering onto the input dataset. This
+is very similar to @code{filter-mean}, except that all outliers (identified
+by the @mymath{\sigma}-clipping algorithm) have been removed, see
address@hidden clipping} for more on the basics of this algorithm. As described
+there, two extra input parameters are necessary for
address@hidden: the multiple of @mymath{\sigma} and the
+termination critera. @code{filter-sigclip-mean} therefore needs to pop two
+other operands from the stack after the dimensions of the box.
+
+For example the line below uses the same box size as the example of
address@hidden However, all elements in the box that are iteratively
+beyond @mymath{3\sigma} of the distribution's median are removed from the
+final calculation of the mean until the change in @mymath{\sigma} is less
+than @mymath{0.2}.
+
address@hidden
+$ astarithmetic 3 0.2 5 4 image.fits filter-sigclip-mean
address@hidden example
+
+The median (which needs a sorted dataset) is necessary for
address@hidden, therefore @code{filter-sigclip-mean} can be
+significantly slower than @code{filter-mean}. However, if there are strong
+outliers in the dataset that you want to ignore (for example emission lines
+on a spectrum when finding the continuum), this is a much better solution.
+
address@hidden filter-sigclip-median
+Apply a @mymath{\sigma}-clipped median filtering onto the input
+dataset. This operator and its necessary operands are almost identical to
address@hidden, except that after @mymath{\sigma}-clipping, the
+median value (which is less affected by outliers than the mean) is added
+back to the stack.
+
@item lt
Less than: If the second popped (or left operand in infix notation, see
@ref{Reverse polish notation}) value is smaller than the first popped
@@ -19844,6 +19877,17 @@ that corresponds to its type. If @code{input} is a
tile over a larger
dataset, only the region that the tile covers will be set to blank.
@end deftypefun
address@hidden int gal_blank_is (void @code{*pointer}, uint8_t @code{type})
+Return 1 if the contents of @code{pointer} (assuming a type of @code{type})
+is blank. Otherwise, return 0. Note that this function only works on one
+element of the given type. So if @code{pointer} is an array, only its first
+element will be checked. Therefore for strings, the type of @code{pointer}
+is assumed to be @code{char *}. To check if an array/dataset has blank
+elements or to find which elements in an array are blank, you can use
address@hidden or @code{gal_blank_flag} respectively (described
+below).
address@hidden deftypefun
+
@deftypefun int gal_blank_present (gal_data_t @code{*input}, int
@code{updateflag})
Return 1 if the dataset has a blank value and zero if it doesn't. Before
@@ -19882,6 +19926,11 @@ adjust the size properly (the number of non-blank
elements). In practice
this function doesn't @code{realloc} the input array, it just shifts the
blank elements to the end and adjusts the size elements of the
@code{gal_data_t}, see @ref{Generic data container}.
+
+If all the elements were blank, then @code{input->size} will be zero. This
+is thus a good parameter to check after calling this function to see if
+there actually were any non-blank elements in the input or not and take the
+appropriate measure. This can help avoid strange bugs in later steps.
@end deftypefun
@deftypefun {char *} gal_blank_as_string (uint8_t @code{type}, int
@code{width})
@@ -23983,8 +24032,10 @@ array[3]: value at the end of symmetricity.
@deftypefun {gal_data_t *} gal_statistics_mode_mirror_plots (gal_data_t
@code{*input}, gal_data_t @code{*value}, size_t @code{numbins}, int
@code{inplace}, double @code{*mirror_val})
Make a mirrored histogram and cumulative frequency plot (with
address@hidden) with the mirror distribution of the @code{input} with a
-value at @code{value}.
address@hidden) with the mirror distribution of the @code{input} having a
+value in @code{value}. If all the input elements are blank, or the mirror
+value is outside the range of the input, this function will return a
address@hidden pointer.
The output is a list of data structures (see @ref{List of gal_data_t}): the
first is the bins with one bin at the mirror point, the second is the
@@ -24010,6 +24061,12 @@ Remove all the blanks and sort the input dataset. If
@code{inplace} is
non-zero this will happen on the input dataset (and the output will be the
same as the input). However, if @code{inplace} is zero, this function will
allocate a new copy of the dataset that is sorted and has no blank values.
+
+If all the elements were blank, then the returned dataset's @code{size}
+will be zero. This is thus a good parameter to check after calling this
+function to see if there actually were any non-blank elements in the input
+or not and take the appropriate measure. This can help avoid strange bugs
+in later steps.
@end deftypefun
@deftypefun {gal_data_t *} gal_statistics_regular_bins (gal_data_t
@code{*input}, gal_data_t @code{*inrange}, size_t @code{numbins}, double
@code{onebinstart})
@@ -24089,8 +24146,8 @@ normalized CFP's maximum value is 1.
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 (@mymath{\sigma} that is used to define outliers in each round of
-clipping).
+deviation (or @mymath{\sigma}, that is used to define outliers in each
+round of clipping).
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
@@ -24105,6 +24162,10 @@ array[1]: Median.
array[2]: Mean.
array[3]: Standard deviation.
@end example
+
+If the @mymath{\sigma}-clipping doesn't converge or all input elements are
+blank, then this function will return NaN values for all the elements
+above.
@end deftypefun
diff --git a/lib/blank.c b/lib/blank.c
index 98b69f3..11cbd6b 100644
--- a/lib/blank.c
+++ b/lib/blank.c
@@ -112,6 +112,57 @@ gal_blank_initialize(gal_data_t *input)
+
+/* Return 1 if the contents of the pointer (with the given type) is
+ blank. */
+int
+gal_blank_is(void *pointer, uint8_t type)
+{
+ switch(type)
+ {
+ /* Numeric types */
+ case GAL_TYPE_UINT8: return *(uint8_t *)pointer==GAL_BLANK_UINT8;
+ case GAL_TYPE_INT8: return *(int8_t *)pointer==GAL_BLANK_INT8;
+ case GAL_TYPE_UINT16: return *(uint16_t *)pointer==GAL_BLANK_UINT16;
+ case GAL_TYPE_INT16: return *(int16_t *)pointer==GAL_BLANK_INT16;
+ case GAL_TYPE_UINT32: return *(uint32_t *)pointer==GAL_BLANK_UINT32;
+ case GAL_TYPE_INT32: return *(int32_t *)pointer==GAL_BLANK_INT32;
+ case GAL_TYPE_UINT64: return *(uint64_t *)pointer==GAL_BLANK_UINT64;
+ case GAL_TYPE_INT64: return *(int64_t *)pointer==GAL_BLANK_INT64;
+ case GAL_TYPE_FLOAT32: return isnan( *(float *)(pointer) );
+ case GAL_TYPE_FLOAT64: return isnan( *(double *)(pointer) );
+
+ /* String. */
+ case GAL_TYPE_STRING: if(!strcmp(pointer,GAL_BLANK_STRING)) return 1;
+
+ /* Complex types */
+ case GAL_TYPE_COMPLEX32:
+ case GAL_TYPE_COMPLEX64:
+ error(EXIT_FAILURE, 0, "%s: complex types are not yet supported",
+ __func__);
+
+ /* Bit. */
+ case GAL_TYPE_BIT:
+ error(EXIT_FAILURE, 0, "%s: bit type datasets are not yet supported",
+ __func__);
+
+ default:
+ error(EXIT_FAILURE, 0, "%s: type value (%d) not recognized",
+ __func__, type);
+ }
+
+ /* Control should not reach here, so print an error if it does, then
+ return a 0 (just to avoid compiler warnings). */
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to address the "
+ "problem. Control should not reach the end of this funciton",
+ __func__, PACKAGE_BUGREPORT);
+ return 0;
+}
+
+
+
+
+
/* Return 1 if the dataset has a blank value and zero if it doesn't. Before
checking the dataset, this function will look at its flags. If the
`GAL_DATA_FLAG_HASBLANK' or `GAL_DATA_FLAG_DONT_CHECK_ZERO' bits of
diff --git a/lib/gnuastro/blank.h b/lib/gnuastro/blank.h
index 674496f..eab1fbf 100644
--- a/lib/gnuastro/blank.h
+++ b/lib/gnuastro/blank.h
@@ -96,6 +96,9 @@ char *
gal_blank_as_string(uint8_t type, int width);
int
+gal_blank_is(void *pointer, uint8_t type);
+
+int
gal_blank_present(gal_data_t *input, int updateflag);
gal_data_t *
diff --git a/lib/statistics.c b/lib/statistics.c
index 405e535..36259ce 100644
--- a/lib/statistics.c
+++ b/lib/statistics.c
@@ -269,7 +269,10 @@ gal_statistics_median(gal_data_t *input, int inplace)
NULL, NULL, NULL);
/* Write the median. */
- statistics_median_in_sorted_no_blank(nbs, out->array);
+ if(nbs->size)
+ statistics_median_in_sorted_no_blank(nbs, out->array);
+ else
+ gal_blank_write(out->array, out->type);
/* Clean up (if necessary), then return the output */
if(nbs!=input) gal_data_free(nbs);
@@ -326,19 +329,26 @@ gal_statistics_quantile(gal_data_t *input, double
quantile, int inplace)
gal_data_t *out=gal_data_alloc(NULL, nbs->type, 1, &dsize,
NULL, 1, -1, NULL, NULL, NULL);
- /* Find the index of the quantile. */
- index=gal_statistics_quantile_index(nbs->size, quantile);
-
- /* Write the value at this index into the output. */
- if(index==GAL_BLANK_SIZE_T)
+ /* Only continue processing if there are non-blank elements. */
+ if(nbs->size)
{
- blank=gal_data_malloc_array(nbs->type, 1, __func__, "blank");
- memcpy(out->array, blank, gal_type_sizeof(nbs->type));
- free(blank);
+ /* Find the index of the quantile. */
+ index=gal_statistics_quantile_index(nbs->size, quantile);
+
+ /* Write the value at this index into the output. */
+ if(index==GAL_BLANK_SIZE_T)
+ {
+ blank=gal_data_malloc_array(nbs->type, 1, __func__, "blank");
+ memcpy(out->array, blank, gal_type_sizeof(nbs->type));
+ free(blank);
+ }
+ else
+ memcpy(out->array,
+ gal_data_ptr_increment(nbs->array, index, nbs->type),
+ gal_type_sizeof(nbs->type));
}
else
- memcpy(out->array, gal_data_ptr_increment(nbs->array, index, nbs->type),
- gal_type_sizeof(nbs->type));
+ gal_blank_write(out->array, out->type);
/* Clean up and return. */
if(nbs!=input) gal_data_free(nbs);
@@ -397,23 +407,27 @@ gal_statistics_quantile_function_index(gal_data_t *input,
gal_data_t *value,
error(EXIT_FAILURE, 0, "%s: the types of the input dataset and requested "
"value have to be the same", __func__);
- /* Find the result: */
- switch(nbs->type)
- {
- case GAL_TYPE_UINT8: STATS_QFUNC_IND( uint8_t ); break;
- case GAL_TYPE_INT8: STATS_QFUNC_IND( int8_t ); break;
- case GAL_TYPE_UINT16: STATS_QFUNC_IND( uint16_t ); break;
- case GAL_TYPE_INT16: STATS_QFUNC_IND( int16_t ); break;
- case GAL_TYPE_UINT32: STATS_QFUNC_IND( uint32_t ); break;
- case GAL_TYPE_INT32: STATS_QFUNC_IND( int32_t ); break;
- case GAL_TYPE_UINT64: STATS_QFUNC_IND( uint64_t ); break;
- case GAL_TYPE_INT64: STATS_QFUNC_IND( int64_t ); break;
- case GAL_TYPE_FLOAT32: STATS_QFUNC_IND( float ); break;
- case GAL_TYPE_FLOAT64: STATS_QFUNC_IND( double ); break;
- default:
- error(EXIT_FAILURE, 0, "%s: type code %d not recognized",
- __func__, nbs->type);
- }
+ /* Only continue processing if we have non-blank elements. */
+ if(nbs->size)
+ /* Find the result: */
+ switch(nbs->type)
+ {
+ case GAL_TYPE_UINT8: STATS_QFUNC_IND( uint8_t ); break;
+ case GAL_TYPE_INT8: STATS_QFUNC_IND( int8_t ); break;
+ case GAL_TYPE_UINT16: STATS_QFUNC_IND( uint16_t ); break;
+ case GAL_TYPE_INT16: STATS_QFUNC_IND( int16_t ); break;
+ case GAL_TYPE_UINT32: STATS_QFUNC_IND( uint32_t ); break;
+ case GAL_TYPE_INT32: STATS_QFUNC_IND( int32_t ); break;
+ case GAL_TYPE_UINT64: STATS_QFUNC_IND( uint64_t ); break;
+ case GAL_TYPE_INT64: STATS_QFUNC_IND( int64_t ); break;
+ case GAL_TYPE_FLOAT32: STATS_QFUNC_IND( float ); break;
+ case GAL_TYPE_FLOAT64: STATS_QFUNC_IND( double ); break;
+ default:
+ error(EXIT_FAILURE, 0, "%s: type code %d not recognized",
+ __func__, nbs->type);
+ }
+ else
+ index=GAL_BLANK_SIZE_T;
/* Clean up and return. */
if(nbs!=input) gal_data_free(nbs);
@@ -447,32 +461,38 @@ gal_statistics_quantile_function(gal_data_t *input,
gal_data_t *value,
NULL, 1, -1, NULL, NULL, NULL);
size_t ind=gal_statistics_quantile_function_index(input, value, inplace);
- /* Note that counting of the index starts from 0, so for the quantile we
- should divided by (size - 1). */
- d=out->array;
- if(ind==GAL_BLANK_SIZE_T)
+ /* Only continue processing if there are non-blank values. */
+ if(nbs->size)
{
- /* See if the value is larger or smaller than the input's minimum or
- maximum. */
- switch(nbs->type)
+ /* Note that counting of the index starts from 0, so for the quantile
+ we should divided by (size - 1). */
+ d=out->array;
+ if(ind==GAL_BLANK_SIZE_T)
{
- case GAL_TYPE_UINT8: STATS_QFUNC( uint8_t ); break;
- case GAL_TYPE_INT8: STATS_QFUNC( int8_t ); break;
- case GAL_TYPE_UINT16: STATS_QFUNC( uint16_t ); break;
- case GAL_TYPE_INT16: STATS_QFUNC( int16_t ); break;
- case GAL_TYPE_UINT32: STATS_QFUNC( uint32_t ); break;
- case GAL_TYPE_INT32: STATS_QFUNC( int32_t ); break;
- case GAL_TYPE_UINT64: STATS_QFUNC( uint64_t ); break;
- case GAL_TYPE_INT64: STATS_QFUNC( int64_t ); break;
- case GAL_TYPE_FLOAT32: STATS_QFUNC( float ); break;
- case GAL_TYPE_FLOAT64: STATS_QFUNC( double ); break;
- default:
- error(EXIT_FAILURE, 0, "%s: type code %d not recognized",
- __func__, nbs->type);
+ /* See if the value is larger or smaller than the input's minimum
+ or maximum. */
+ switch(nbs->type)
+ {
+ case GAL_TYPE_UINT8: STATS_QFUNC( uint8_t ); break;
+ case GAL_TYPE_INT8: STATS_QFUNC( int8_t ); break;
+ case GAL_TYPE_UINT16: STATS_QFUNC( uint16_t ); break;
+ case GAL_TYPE_INT16: STATS_QFUNC( int16_t ); break;
+ case GAL_TYPE_UINT32: STATS_QFUNC( uint32_t ); break;
+ case GAL_TYPE_INT32: STATS_QFUNC( int32_t ); break;
+ case GAL_TYPE_UINT64: STATS_QFUNC( uint64_t ); break;
+ case GAL_TYPE_INT64: STATS_QFUNC( int64_t ); break;
+ case GAL_TYPE_FLOAT32: STATS_QFUNC( float ); break;
+ case GAL_TYPE_FLOAT64: STATS_QFUNC( double ); break;
+ default:
+ error(EXIT_FAILURE, 0, "%s: type code %d not recognized",
+ __func__, nbs->type);
+ }
}
+ else
+ d[0] = (double)ind / ((double)(nbs->size - 1));
}
else
- d[0] = (double)ind / ((double)(nbs->size - 1));
+ gal_blank_write(out->array, out->type);
/* Clean up and return. */
if(nbs!=input) gal_data_free(nbs);
@@ -1055,6 +1075,8 @@ gal_statistics_mode_mirror_plots(gal_data_t *input,
gal_data_t *value,
gal_data_t *nbs=gal_statistics_no_blank_sorted(input, inplace);
size_t ind=gal_statistics_quantile_function_index(nbs, value, inplace);
+ /* Only continue if we actually have non-blank elements. */
+ if(nbs->size==0) return NULL;
/* If the given mirror was outside the range of the input, then index
will be 0 (below the range) or -1 (above the range), in that case, we
@@ -1244,8 +1266,7 @@ gal_statistics_sort_decreasing(gal_data_t *input)
This function can also work on tiles, in that case, `inplace' is
useless, because a tile doesn't own its dataset and the dataset is not
- contiguous.
-*/
+ contiguous. */
gal_data_t *
gal_statistics_no_blank_sorted(gal_data_t *input, int inplace)
{
@@ -1290,27 +1311,33 @@ gal_statistics_no_blank_sorted(gal_data_t *input, int
inplace)
else noblank=contig;
- /* Make sure the array is sorted. After this step, we won't be dealing
- with `noblank' any more but with `sorted'. */
- sortstatus=gal_statistics_is_sorted(noblank);
- if( sortstatus )
+ /* If there are any non-blank elements, make sure the array is
+ sorted. After this step, we won't be dealing with `noblank' any more
+ but with `sorted'. */
+ if(noblank->size)
{
- sorted=noblank;
- sorted->status=sortstatus;
- }
- else
- {
- if(inplace) sorted=noblank;
+ sortstatus=gal_statistics_is_sorted(noblank);
+ if( sortstatus )
+ {
+ sorted=noblank;
+ sorted->status=sortstatus;
+ }
else
{
- if(noblank!=input) /* no-blank has already been allocated. */
- sorted=noblank;
+ if(inplace) sorted=noblank;
else
- sorted=gal_data_copy(noblank);
+ {
+ if(noblank!=input) /* no-blank has already been allocated. */
+ sorted=noblank;
+ else
+ sorted=gal_data_copy(noblank);
+ }
+ gal_statistics_sort_increasing(sorted);
+ sorted->status=GAL_STATISTICS_SORTED_INCREASING;
}
- gal_statistics_sort_increasing(sorted);
- sorted->status=GAL_STATISTICS_SORTED_INCREASING;
}
+ else
+ sorted=noblank;
/* Return final array. */
@@ -1799,6 +1826,7 @@ gal_data_t *
gal_statistics_sigma_clip(gal_data_t *input, float multip, float param,
int inplace, int quiet)
{
+ float *oa;
void *start, *nbs_array;
double *med, *mean, *std;
uint8_t bytolerance = param>=1.0f ? 0 : 1;
@@ -1830,95 +1858,104 @@ gal_statistics_sigma_clip(gal_data_t *input, float
multip, float param,
NULL, NULL, NULL);
- /* Print the comments. */
- 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. */
- size=nbs->size;
- sortstatus=nbs->status;
- nbs_array=start=nbs->array;
- while(num<maxnum)
+ /* Only continue processing if we have non-blank elements. */
+ oa=out->array;
+ nbs_array=nbs->array;
+ if(nbs->size==0)
{
- /* Find the median. */
- 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 average and Standard deviation, note that both `start'
- and `size' will be different in the next round. */
- nbs->array = start;
- nbs->size = oldsize = size;
- meanstd=gal_statistics_mean_std(nbs);
-
- /* Put the three final values in usable (with a type) pointers. */
- med = median_d->array;
- mean = meanstd->array;
- std = &((double *)(meanstd->array))[1];
-
- /* 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);
-
- /* If we are to work by tolerance, then 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! */
- if( bytolerance && num>0 && ((oldstd - *std) / *std) < param )
- break;
-
- /* Clip all the elements outside of the desired range: since the
- array is sorted, this means to just change the starting pointer
- and size of the array. */
- switch(type)
+ printf("NO SIGMA-CLIPPING: all input elements are blank.");
+ oa[0] = oa[1] = oa[2] = oa[3] = NAN;
+ }
+ else
+ {
+ /* Print the comments. */
+ 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. */
+ size=nbs->size;
+ start=nbs->array;
+ sortstatus=nbs->status;
+ while(num<maxnum)
{
- 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;
- default:
- error(EXIT_FAILURE, 0, "%s: type code %d not recognized",
- __func__, type);
- }
+ /* Find the median. */
+ 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 average and Standard deviation, note that both
+ `start' and `size' will be different in the next round. */
+ nbs->array = start;
+ nbs->size = oldsize = size;
+ meanstd=gal_statistics_mean_std(nbs);
+
+ /* Put the three final values in usable (with a type)
+ pointers. */
+ med = median_d->array;
+ mean = meanstd->array;
+ std = &((double *)(meanstd->array))[1];
+
+ /* 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);
+
+ /* If we are to work by tolerance, then 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! */
+ if( bytolerance && num>0 && ((oldstd - *std) / *std) < param )
+ break;
+
+ /* Clip all the elements outside of the desired range: since the
+ array is sorted, this means to just change the starting
+ 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;
+ default:
+ error(EXIT_FAILURE, 0, "%s: type code %d not recognized",
+ __func__, type);
+ }
- /* 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;
- ++num;
+ /* 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;
+ ++num;
- /* Clean up: */
- gal_data_free(meanstd);
- gal_data_free(median_d);
- }
+ /* Clean up: */
+ gal_data_free(meanstd);
+ gal_data_free(median_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. */
- out->status=num;
- if( bytolerance && num==maxnum )
- {
- ((float *)(out->array))[0] = NAN;
- ((float *)(out->array))[1] = NAN;
- ((float *)(out->array))[2] = NAN;
- ((float *)(out->array))[3] = NAN;
- }
- else
- {
- ((float *)(out->array))[0] = size;
- ((float *)(out->array))[1] = oldmed;
- ((float *)(out->array))[2] = oldmean;
- ((float *)(out->array))[3] = oldstd;
+ /* 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. */
+ out->status=num;
+ if( bytolerance && num==maxnum )
+ oa[0] = oa[1] = oa[2] = oa[3] = NAN;
+ else
+ {
+ oa[0] = size;
+ oa[1] = oldmed;
+ oa[2] = oldmean;
+ oa[3] = oldstd;
+ }
}
/* Clean up and return. */
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [gnuastro-commits] master ce3ebeb: New sigma-clipped mean and median filters for Arithmetic,
Mohammad Akhlaghi <=