[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 04f4225: Arithmetic: New sigma-clipping coaddi
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 04f4225: Arithmetic: New sigma-clipping coadding operators |
|
Date: |
Mon, 4 Feb 2019 19:43:32 -0500 (EST) |
branch: master
commit 04f4225ff5009c8098b92d637b2f1a198b528b16
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
Arithmetic: New sigma-clipping coadding operators
Until now Arithmetic only had relatively basic operators for coadding
several datasets into one like the mean, median and standard
deviation.
With this commit, sigma-clipping has also been added to reject outliers on
a per-pixel basis.
Also, while doing the checks for this new feature, I noticed that in the
sigma-clipping function, we were mistakenly finding the median before
correcting the size of the dataset, causing some inaccurate results in
special conditions.
---
NEWS | 9 ++++
bin/arithmetic/arithmetic.c | 90 ++++++++++++++++++++++++--------
doc/announce-acknowledge.txt | 1 +
doc/gnuastro.texi | 60 +++++++++++++++++++---
lib/arithmetic.c | 119 ++++++++++++++++++++++++++++++++++++++-----
lib/gnuastro/arithmetic.h | 8 ++-
lib/statistics.c | 24 ++++++---
7 files changed, 260 insertions(+), 51 deletions(-)
diff --git a/NEWS b/NEWS
index 177e78a..0190a3b 100644
--- a/NEWS
+++ b/NEWS
@@ -17,6 +17,11 @@ GNU Astronomy Utilities NEWS -*-
outline -*-
without changing the state of the operand stack (popping the top
element). This can greatly help in debugging/understanding an
Arithmetic command, especially as it gets longer, or more complicated.
+ - Four new operators have beed added to allow coadding multiple datasets
+ into one using sigma-clipping: `sigclip-number', `sigclip-mean',
+ `sigclip-median', and `sigclip-std'. These are very useful when
+ several inputs have strong outliers that affect the median, or the
+ mean is required.
--wcsfile and --wcshdu: these two options can be used to specify a
different file for reading the WCS of the output. This is useful when
the default (theh WCS of the first dataset that is read) is not the
@@ -56,6 +61,10 @@ GNU Astronomy Utilities NEWS -*-
outline -*-
** Changed features
+ Arithmetic:
+ - `num' operator is renamed to `number'.
+ - `numvalue' operator is renamed to `numbervalue'.
+
ConvertType:
--forcemin & --forcemax: until now, `--flminbyte' and `--flmaxbyte' were
used to force the range of conversion to color channels (even if the
diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index 886d61c..986fda2 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -64,21 +64,60 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
/************* Internal functions *************/
/***************************************************************/
#define SET_NUM_OP(CTYPE) { \
- CTYPE a=*(CTYPE *)(data->array); if(a>0) return a; }
+ CTYPE a=*(CTYPE *)(numpop->array); if(a>0) return a; }
static size_t
-pop_number_of_operands(struct arithmeticparams *p, gal_data_t *data,
- char *token_string)
+pop_number_of_operands(struct arithmeticparams *p, int op, char *token_string,
+ gal_data_t **params)
{
+ char *cstring="first";
+ size_t c, numparams=0;
+ gal_data_t *tmp, *numpop;
+
+ /* See if this operator needs any parameters. If so, pop them. */
+ switch(op)
+ {
+ case GAL_ARITHMETIC_OP_SIGCLIP_STD:
+ case GAL_ARITHMETIC_OP_SIGCLIP_MEAN:
+ case GAL_ARITHMETIC_OP_SIGCLIP_MEDIAN:
+ case GAL_ARITHMETIC_OP_SIGCLIP_NUMBER:
+ numparams=2;
+ break;
+ }
+
+ /* If any parameters should be read, read them. */
+ *params=NULL;
+ for(c=0;c<numparams;++c)
+ {
+ /* If it only has one element, save it as floating point and add it
+ to the list. */
+ tmp=operands_pop(p, token_string);
+ if(tmp->size>1)
+ error(EXIT_FAILURE, 0, "the %s popped operand of the \"%s\" "
+ "operator must be a single number", cstring, token_string);
+ tmp=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT32);
+ gal_list_data_add(params, tmp);
+
+ /* A small sanity check (none of the parameters for sigma-clipping
+ can be negative.. */
+ if( ((float *)(tmp->array))[0]<=0.0 )
+ error(EXIT_FAILURE, 0, "the %s popped operand of the \"%s\" "
+ "operator cannot be negative", cstring, token_string);
+
+ /* Increment the counter string. */
+ cstring=c?"third":"second";
+ }
+
/* Check if its a number. */
- if(data->size>1)
- error(EXIT_FAILURE, 0, "the first popped operand to the \"%s\" "
- "operator must be a number, not an array", token_string);
+ numpop=operands_pop(p, token_string);
+ if(numpop->size>1)
+ error(EXIT_FAILURE, 0, "the %s popped operand of the \"%s\" "
+ "operator (number of input datasets) must be a number, not an "
+ "array", cstring, token_string);
/* Check its type and return the value. */
- switch(data->type)
+ switch(numpop->type)
{
-
/* For the integer types, if they are unsigned, then just pass their
value, if they are signed, you have to make sure they are zero or
positive. */
@@ -94,18 +133,20 @@ pop_number_of_operands(struct arithmeticparams *p,
gal_data_t *data,
/* Floating point numbers are not acceptable in this context. */
case GAL_TYPE_FLOAT32:
case GAL_TYPE_FLOAT64:
- error(EXIT_FAILURE, 0, "the first popped operand to the \"%s\" "
- "operator must be an integer type", token_string);
+ error(EXIT_FAILURE, 0, "the %s popped operand of the \"%s\" "
+ "operator (number of input datasets) must be an integer type",
+ cstring, token_string);
default:
error(EXIT_FAILURE, 0, "%s: type code %d not recognized",
- __func__, data->type);
+ __func__, numpop->type);
}
/* If control reaches here, then the number must have been a negative
value, so print an error. */
- error(EXIT_FAILURE, 0, "the first popped operand to the \"%s\" operator "
- "cannot be zero or a negative number", token_string);
+ error(EXIT_FAILURE, 0, "the %s popped operand of the \"%s\" operator "
+ "cannot be zero or a negative number", cstring,
+ token_string);
return 0;
}
@@ -923,8 +964,8 @@ reversepolish(struct arithmeticparams *p)
{ op=GAL_ARITHMETIC_OP_MINVAL; nop=1; }
else if (!strcmp(token->v, "maxvalue"))
{ op=GAL_ARITHMETIC_OP_MAXVAL; nop=1; }
- else if (!strcmp(token->v, "numvalue"))
- { op=GAL_ARITHMETIC_OP_NUMVAL; nop=1; }
+ else if (!strcmp(token->v, "numbervalue"))
+ { op=GAL_ARITHMETIC_OP_NUMBERVAL; nop=1; }
else if (!strcmp(token->v, "sumvalue"))
{ op=GAL_ARITHMETIC_OP_SUMVAL; nop=1; }
else if (!strcmp(token->v, "meanvalue"))
@@ -937,8 +978,8 @@ reversepolish(struct arithmeticparams *p)
{ op=GAL_ARITHMETIC_OP_MIN; nop=-1; }
else if (!strcmp(token->v, "max"))
{ op=GAL_ARITHMETIC_OP_MAX; nop=-1; }
- else if (!strcmp(token->v, "num"))
- { op=GAL_ARITHMETIC_OP_NUM; nop=-1; }
+ else if (!strcmp(token->v, "number"))
+ { op=GAL_ARITHMETIC_OP_NUMBER; nop=-1; }
else if (!strcmp(token->v, "sum"))
{ op=GAL_ARITHMETIC_OP_SUM; nop=-1; }
else if (!strcmp(token->v, "mean"))
@@ -947,6 +988,14 @@ reversepolish(struct arithmeticparams *p)
{ op=GAL_ARITHMETIC_OP_STD; nop=-1; }
else if (!strcmp(token->v, "median"))
{ op=GAL_ARITHMETIC_OP_MEDIAN; nop=-1; }
+ else if (!strcmp(token->v, "sigclip-number"))
+ { op=GAL_ARITHMETIC_OP_SIGCLIP_NUMBER; nop=-1; }
+ else if (!strcmp(token->v, "sigclip-mean"))
+ { op=GAL_ARITHMETIC_OP_SIGCLIP_MEAN; nop=-1; }
+ else if (!strcmp(token->v, "sigclip-median"))
+ { op=GAL_ARITHMETIC_OP_SIGCLIP_MEDIAN; nop=-1; }
+ else if (!strcmp(token->v, "sigclip-std"))
+ { op=GAL_ARITHMETIC_OP_SIGCLIP_STD; nop=-1; }
/* Conditional operators. */
else if (!strcmp(token->v, "lt" ))
@@ -1075,13 +1124,13 @@ reversepolish(struct arithmeticparams *p)
case -1:
/* This case is when the number of operands is itself an
- operand. So the first popped operand must be an
+ operand. So except for sigma-clipping (that has other
+ parameters), the first popped operand must be an
integer number, we will use that to construct a linked
list of any number of operands within the single `d1'
pointer. */
d1=NULL;
- numop=pop_number_of_operands(p, operands_pop(p, token->v),
- token->v);
+ numop=pop_number_of_operands(p, op, token->v, &d2);
for(i=0;i<numop;++i)
gal_list_data_add(&d1, operands_pop(p, token->v));
break;
@@ -1091,7 +1140,6 @@ reversepolish(struct arithmeticparams *p)
"operand(s)", token->v, nop);
}
-
/* Run the arithmetic operation. Note that `gal_arithmetic'
is a variable argument function (like printf). So the
number of arguments it uses depend on the operator. So
diff --git a/doc/announce-acknowledge.txt b/doc/announce-acknowledge.txt
index 5669500..e5c7f79 100644
--- a/doc/announce-acknowledge.txt
+++ b/doc/announce-acknowledge.txt
@@ -2,3 +2,4 @@ Alphabetically ordered list to acknowledge in the next release.
Raúl Infante Sainz
David Valls-Gabaud
+Ignacio Trujillo
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 8143082..1739647 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -12180,7 +12180,7 @@ output of this operand is in the same type as the input.
Maximum (non-blank) value of first operand in the same type, similar to
@command{minvalue}.
address@hidden numvalue
address@hidden numbervalue
Number of non-blank elements in first operand in the @code{uint64} type,
similar to @command{minvalue}.
@@ -12220,8 +12220,7 @@ Important notes:
@itemize
@item
-NaN/blank pixels will be ignored, see @ref{Blank
-pixels}.
+NaN/blank pixels will be ignored, see @ref{Blank pixels}.
@item
The output will have the same type as the inputs. This is natural for the
@@ -12237,7 +12236,7 @@ conversion will be used.
Similar to @command{min}, but the pixels of the output will contain
the maximum of the respective pixels in all operands in the stack.
address@hidden num
address@hidden number
Similar to @command{min}, but the pixels of the output will contain the
number of the respective non-blank pixels in all input operands.
@@ -12257,6 +12256,43 @@ standard deviation of the respective pixels in all
operands in the stack.
Similar to @command{min}, but the pixels of the output will contain
the median of the respective pixels in all operands in the stack.
address@hidden sigclip-number
+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 (paramters for 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 1 it is interpretted 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 third popped operand
+(@command{10}) is number of datasets that will be used (similar to the
+first popped operand to @command{min}).
+
address@hidden
+astarithmetic a.fits b.fits c.fits 3 5 0.2 sigclip-number
address@hidden example
+
address@hidden sigclip-median
+Combine the specified number of inputs into a single output that contains
+the @emph{median} after @mymath{\sigma}-clipping on each element/pixel. For
+more see @command{sigclip-number}.
+
address@hidden sigclip-mean
+Combine the specified number of inputs into a single output that contains
+the @emph{mean} after @mymath{\sigma}-clipping on each element/pixel. For
+more see @command{sigclip-number}.
+
address@hidden sigclip-std
+Combine the specified number of inputs into a single output that contains
+the @emph{standard deviation} after @mymath{\sigma}-clipping on each
+element/pixel. For more see @command{sigclip-number}.
+
@item filter-mean
Apply mean filtering (or @url{https://en.wikipedia.org/wiki/Moving_average,
moving average}) on the input dataset. During mean filtering, each pixel
@@ -26867,7 +26903,7 @@ type, you can convert the input to a floating point
type with
@deffn Macro GAL_ARITHMETIC_OP_MINVAL
@deffnx Macro GAL_ARITHMETIC_OP_MAXVAL
address@hidden Macro GAL_ARITHMETIC_OP_NUMVAL
address@hidden Macro GAL_ARITHMETIC_OP_NUMBERVAL
@deffnx Macro GAL_ARITHMETIC_OP_SUMVAL
@deffnx Macro GAL_ARITHMETIC_OP_MEANVAL
@deffnx Macro GAL_ARITHMETIC_OP_STDVAL
@@ -26886,7 +26922,7 @@ Unary operand absolute-value operator.
@deffn Macro GAL_ARITHMETIC_OP_MIN
@deffnx Macro GAL_ARITHMETIC_OP_MAX
address@hidden Macro GAL_ARITHMETIC_OP_NUM
address@hidden Macro GAL_ARITHMETIC_OP_NUMBER
@deffnx Macro GAL_ARITHMETIC_OP_SUM
@deffnx Macro GAL_ARITHMETIC_OP_MEAN
@deffnx Macro GAL_ARITHMETIC_OP_STD
@@ -26899,6 +26935,18 @@ respective statistical operation on the whole list.
See the discussion
under the @code{min} operator in @ref{Arithmetic operators}.
@end deffn
address@hidden Macro GAL_ARITHMETIC_OP_SIGCLIP_STD
address@hidden Macro GAL_ARITHMETIC_OP_SIGCLIP_MEAN
address@hidden Macro GAL_ARITHMETIC_OP_SIGCLIP_MEDIAN
address@hidden Macro GAL_ARITHMETIC_OP_SIGCLIP_NUMBER
+Similar to the operands above (including @code{GAL_ARITHMETIC_MIN}), except
+that when @code{gal_arithmetic} is called with these operators, it requires
+two arguments. The first is the list of datasets like before, and the
+second is the 2-element list of sigma-clipping parameters. The first
+element in the parameters list is the multiple of sigma and the second is
+the termination criteria (see @ref{Sigma clipping}).
address@hidden deffn
+
@deffn Macro GAL_ARITHMETIC_OP_POW
Binary operator to-power operator. When @code{gal_arithmetic} is called
with any of these operators, it will expect two operands: raising the first
diff --git a/lib/arithmetic.c b/lib/arithmetic.c
index cb975be..d24e960 100644
--- a/lib/arithmetic.c
+++ b/lib/arithmetic.c
@@ -28,6 +28,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <stdlib.h>
#include <stdarg.h>
+#include <gnuastro/list.h>
#include <gnuastro/blank.h>
#include <gnuastro/qsort.h>
#include <gnuastro/pointer.h>
@@ -454,12 +455,12 @@ arithmetic_from_statistics(int operator, int flags,
gal_data_t *input)
switch(operator)
{
- case GAL_ARITHMETIC_OP_MINVAL: out=gal_statistics_minimum(input); break;
- case GAL_ARITHMETIC_OP_MAXVAL: out=gal_statistics_maximum(input); break;
- case GAL_ARITHMETIC_OP_NUMVAL: out=gal_statistics_number(input); break;
- case GAL_ARITHMETIC_OP_SUMVAL: out=gal_statistics_sum(input); break;
- case GAL_ARITHMETIC_OP_MEANVAL: out=gal_statistics_mean(input); break;
- case GAL_ARITHMETIC_OP_STDVAL: out=gal_statistics_std(input); break;
+ case GAL_ARITHMETIC_OP_MINVAL: out=gal_statistics_minimum(input);break;
+ case GAL_ARITHMETIC_OP_MAXVAL: out=gal_statistics_maximum(input);break;
+ case GAL_ARITHMETIC_OP_NUMBERVAL:out=gal_statistics_number(input); break;
+ case GAL_ARITHMETIC_OP_SUMVAL: out=gal_statistics_sum(input); break;
+ case GAL_ARITHMETIC_OP_MEANVAL: out=gal_statistics_mean(input); break;
+ case GAL_ARITHMETIC_OP_STDVAL: out=gal_statistics_std(input); break;
case GAL_ARITHMETIC_OP_MEDIANVAL:
out=gal_statistics_median(input, ip); break;
default:
@@ -836,6 +837,60 @@ arithmetic_where(int flags, gal_data_t *out, gal_data_t
*cond,
+#define MULTIOPERAND_SIGCLIP(TYPE) { \
+ float *sarr; \
+ size_t n, j=0; \
+ gal_data_t *sclip; \
+ TYPE *pixs=gal_pointer_allocate(list->type, dnum, 0, __func__, \
+ "pixs"); \
+ gal_data_t *cont=gal_data_alloc(pixs, list->type, 1, &dnum, NULL, \
+ 0, -1, NULL, NULL, NULL); \
+ \
+ /* Loop over each pixel */ \
+ do \
+ { \
+ /* Initialize. */ \
+ n=0; \
+ \
+ /* Loop over each array. */ \
+ for(i=0;i<dnum;++i) pixs[n++]=a[i][j]; \
+ \
+ /* If there are any elements, measure the */ \
+ if(n) \
+ { \
+ sclip=gal_statistics_sigma_clip(cont, p1, p2, 1, 1); \
+ sarr=sclip->array; \
+ switch(operator) \
+ { \
+ case GAL_ARITHMETIC_OP_SIGCLIP_STD: *o++=sarr[3]; break;\
+ case GAL_ARITHMETIC_OP_SIGCLIP_MEAN: *o++=sarr[2]; break;\
+ case GAL_ARITHMETIC_OP_SIGCLIP_MEDIAN: *o++=sarr[1]; break;\
+ case GAL_ARITHMETIC_OP_SIGCLIP_NUMBER: *o++=sarr[0]; break;\
+ default: \
+ error(EXIT_FAILURE, 0, "%s: a bug! the code %d is not " \
+ "valid for sigma-clipping results", __func__, \
+ operator); \
+ } \
+ \
+ /* Since we are doing sigma-clipping in place, the size, */ \
+ /* and flags need to be reset. */ \
+ cont->flag=0; \
+ cont->size=cont->dsize[0]=dnum; \
+ } \
+ else \
+ *o++=b; \
+ ++j; \
+ } \
+ while(o<of); \
+ \
+ /* Clean up. */ \
+ gal_data_free(cont); \
+ }
+
+
+
+
+
#define MULTIOPERAND_TYPE_SET(TYPE, QSORT_F) { \
TYPE b, **a, *o=out->array, *of=o+out->size; \
size_t i=0; /* Different from the `i' in the main function. */ \
@@ -865,7 +920,7 @@ arithmetic_where(int flags, gal_data_t *out, gal_data_t
*cond,
MULTIOPERAND_MAX(TYPE); \
break; \
\
- case GAL_ARITHMETIC_OP_NUM: \
+ case GAL_ARITHMETIC_OP_NUMBER: \
MULTIOPERAND_NUM; \
break; \
\
@@ -885,6 +940,13 @@ arithmetic_where(int flags, gal_data_t *out, gal_data_t
*cond,
MULTIOPERAND_MEDIAN(TYPE, QSORT_F); \
break; \
\
+ 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); \
+ break; \
+ \
default: \
error(EXIT_FAILURE, 0, "%s: operator code %d not recognized", \
"MULTIOPERAND_TYPE_SET", operator); \
@@ -900,12 +962,14 @@ arithmetic_where(int flags, gal_data_t *out, gal_data_t
*cond,
/* 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.*/
+ 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)
+arithmetic_multioperand(int operator, int flags, gal_data_t *list,
+ gal_data_t *params)
{
uint8_t *hasblank;
size_t i=0, dnum=1;
+ float p1=NAN, p2=NAN;
gal_data_t *out, *tmp, *ttmp;
@@ -914,6 +978,23 @@ arithmetic_multioperand(int operator, int flags,
gal_data_t *list)
if(list==NULL) return NULL;
+ /* If any parameters are given, prepare them. */
+ for(tmp=params; tmp!=NULL; tmp=tmp->next)
+ {
+ /* Basic sanity checks. */
+ if(tmp->size>1)
+ error(EXIT_FAILURE, 0, "%s: parameters must be a single number",
+ __func__);
+ if(tmp->type!=GAL_TYPE_FLOAT32)
+ error(EXIT_FAILURE, 0, "%s: parameters must be float32 type",
+ __func__);
+
+ /* Write them */
+ if(isnan(p1)) p1=((float *)(tmp->array))[0];
+ else p2=((float *)(tmp->array))[0];
+ }
+
+
/* 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)
@@ -1002,6 +1083,7 @@ arithmetic_multioperand(int operator, int flags,
gal_data_t *list)
if(tmp!=out) gal_data_free(tmp);
tmp=ttmp;
}
+ if(params) gal_list_data_free(params);
}
free(hasblank);
return out;
@@ -1369,7 +1451,7 @@ gal_arithmetic_operator_string(int operator)
case GAL_ARITHMETIC_OP_MINVAL: return "minvalue";
case GAL_ARITHMETIC_OP_MAXVAL: return "maxvalue";
- case GAL_ARITHMETIC_OP_NUMVAL: return "numvalue";
+ case GAL_ARITHMETIC_OP_NUMBERVAL: return "numbervalue";
case GAL_ARITHMETIC_OP_SUMVAL: return "sumvalue";
case GAL_ARITHMETIC_OP_MEANVAL: return "meanvalue";
case GAL_ARITHMETIC_OP_STDVAL: return "stdvalue";
@@ -1377,11 +1459,15 @@ gal_arithmetic_operator_string(int operator)
case GAL_ARITHMETIC_OP_MIN: return "min";
case GAL_ARITHMETIC_OP_MAX: return "max";
- case GAL_ARITHMETIC_OP_NUM: return "num";
+ case GAL_ARITHMETIC_OP_NUMBER: return "number";
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_MEDIAN: return "median";
+ 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_TO_UINT8: return "uchar";
case GAL_ARITHMETIC_OP_TO_INT8: return "char";
@@ -1471,7 +1557,7 @@ gal_arithmetic(int operator, int flags, ...)
/* Statistical operators that return one value. */
case GAL_ARITHMETIC_OP_MINVAL:
case GAL_ARITHMETIC_OP_MAXVAL:
- case GAL_ARITHMETIC_OP_NUMVAL:
+ case GAL_ARITHMETIC_OP_NUMBERVAL:
case GAL_ARITHMETIC_OP_SUMVAL:
case GAL_ARITHMETIC_OP_MEANVAL:
case GAL_ARITHMETIC_OP_STDVAL:
@@ -1489,13 +1575,18 @@ gal_arithmetic(int operator, int flags, ...)
/* Multi-operand operators */
case GAL_ARITHMETIC_OP_MIN:
case GAL_ARITHMETIC_OP_MAX:
- case GAL_ARITHMETIC_OP_NUM:
+ case GAL_ARITHMETIC_OP_NUMBER:
case GAL_ARITHMETIC_OP_SUM:
case GAL_ARITHMETIC_OP_MEAN:
case GAL_ARITHMETIC_OP_STD:
case GAL_ARITHMETIC_OP_MEDIAN:
+ case GAL_ARITHMETIC_OP_SIGCLIP_STD:
+ case GAL_ARITHMETIC_OP_SIGCLIP_MEAN:
+ case GAL_ARITHMETIC_OP_SIGCLIP_MEDIAN:
+ case GAL_ARITHMETIC_OP_SIGCLIP_NUMBER:
d1 = va_arg(va, gal_data_t *);
- out=arithmetic_multioperand(operator, flags, d1);
+ d2 = va_arg(va, gal_data_t *);
+ out=arithmetic_multioperand(operator, flags, d1, d2);
break;
diff --git a/lib/gnuastro/arithmetic.h b/lib/gnuastro/arithmetic.h
index 56beedf..db65bc8 100644
--- a/lib/gnuastro/arithmetic.h
+++ b/lib/gnuastro/arithmetic.h
@@ -108,7 +108,7 @@ enum gal_arithmetic_operators
GAL_ARITHMETIC_OP_MINVAL, /* Minimum value of array. */
GAL_ARITHMETIC_OP_MAXVAL, /* Maximum value of array. */
- GAL_ARITHMETIC_OP_NUMVAL, /* Number of (non-blank) elements. */
+ GAL_ARITHMETIC_OP_NUMBERVAL, /* Number of (non-blank) elements. */
GAL_ARITHMETIC_OP_SUMVAL, /* Sum of (non-blank) elements. */
GAL_ARITHMETIC_OP_MEANVAL, /* Mean value of array. */
GAL_ARITHMETIC_OP_STDVAL, /* Standard deviation value of array. */
@@ -116,11 +116,15 @@ enum gal_arithmetic_operators
GAL_ARITHMETIC_OP_MIN, /* Minimum per pixel of multiple arrays. */
GAL_ARITHMETIC_OP_MAX, /* Maximum per pixel of multiple arrays. */
- GAL_ARITHMETIC_OP_NUM, /* Non-blank number of pixels in arrays. */
+ GAL_ARITHMETIC_OP_NUMBER, /* Non-blank number of pixels in arrays. */
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_MEDIAN, /* Median 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_TO_UINT8, /* Convert to uint8_t. */
GAL_ARITHMETIC_OP_TO_INT8, /* Convert to int8_t. */
diff --git a/lib/statistics.c b/lib/statistics.c
index 5ba3c19..4354261 100644
--- a/lib/statistics.c
+++ b/lib/statistics.c
@@ -1397,7 +1397,6 @@ gal_statistics_no_blank_sorted(gal_data_t *input, int
inplace)
}
else noblank=contig;
-
/* 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'. */
@@ -2012,20 +2011,29 @@ gal_statistics_sigma_clip(gal_data_t *input, float
multip, float param,
start=nbs->array;
while(num<maxnum && size)
{
- /* 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;
+
+ /* For a detailed check, just correct the type).
+ if(!quiet)
+ {
+ size_t iii;
+ printf("nbs->size: %zu\n", nbs->size);
+ for(iii=0;iii<nbs->size;++iii)
+ printf("%f\n", ((float *)(nbs->array))[iii]);
+ }
+ */
+
+ /* 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);
- /* Put the three final values in usable (with a type)
- pointers. */
- med = median_d->array;
+ /* Put them in usable (with a type) pointers. */
mean = meanstd->array;
+ med = median_d->array;
std = &((double *)(meanstd->array))[1];
/* If the user wanted to view the steps, show it to them. */
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [gnuastro-commits] master 04f4225: Arithmetic: New sigma-clipping coadding operators,
Mohammad Akhlaghi <=