[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 2dba72f 101/113: Imported work in master, conf
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 2dba72f 101/113: Imported work in master, conflicts fixed, changes made |
|
Date: |
Fri, 16 Apr 2021 10:34:00 -0400 (EDT) |
branch: master
commit 2dba72f627c5c83a8bc094b617f4af91ada8f6e5
Merge: fbbc2fc 97845f9
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Imported work in master, conflicts fixed, changes made
A minor conflict came up with the checking of the new `pseudoconcomp'
option value. It was also necessary to correct the usage of its value which
wasn't a conflict.
---
NEWS | 27 ++++++++
bin/arithmetic/arithmetic.c | 90 ++++++++++++++++++++-------
bin/noisechisel/args.h | 21 +++++--
bin/noisechisel/astnoisechisel.conf | 1 +
bin/noisechisel/detection.c | 3 +-
bin/noisechisel/main.h | 3 +-
bin/noisechisel/noisechisel.c | 4 +-
bin/noisechisel/sky.c | 8 +--
bin/noisechisel/threshold.c | 38 ++++++------
bin/noisechisel/ui.c | 5 +-
bin/noisechisel/ui.h | 3 +-
bin/statistics/args.h | 8 +--
bin/statistics/main.h | 2 +-
bin/statistics/sky.c | 20 +++---
bin/statistics/statistics.c | 3 +-
bin/statistics/ui.h | 2 +-
configure.ac | 32 ++++++----
doc/announce-acknowledge.txt | 2 +
doc/gnuastro.texi | 117 ++++++++++++++++++++++++++++-------
lib/arithmetic.c | 119 +++++++++++++++++++++++++++++++-----
lib/box.c | 12 +++-
lib/gnuastro/arithmetic.h | 8 ++-
lib/statistics.c | 24 +++++---
23 files changed, 424 insertions(+), 128 deletions(-)
diff --git a/NEWS b/NEWS
index 3290258..6bacbd9 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
@@ -45,6 +50,13 @@ GNU Astronomy Utilities NEWS -*-
outline -*-
second input. This greatly simplifies the mergining of different table
columns into one.
+ NoiseChisel:
+ --pseudoconcomp: allows setting the connectivity (4 or 8, in a 2D image)
+ to define separate pseudo-detections. If its stronger,
+ pseudo-detections that are touching on the corner will be identified
+ as one. Until this version, this was hard-written into the code and
+ was the weakest connectivity (4-connected in a 2D image).
+
Library:
GAL_BLANK_LONG: new macro for the `long' type blank value.
GAL_BLANK_ULONG: new macro for the `unsigned long' type blank value.
@@ -56,6 +68,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
@@ -70,6 +86,17 @@ GNU Astronomy Utilities NEWS -*-
outline -*-
on, this option measures the square root of the mean variance, or root
mean square (the correct definition of the Sky standard deviation).
+ NoiseChisel:
+ --ignoreblankintiles: Until now `--ignoreblankinsky', would specify if
+ blank values should also be written into the tiled Sky and Sky
+ standard deviation outputs. But NoiseChisel can optionally produce
+ many more tiled outputs (for example with `--checkqthresh'). So the
+ option was renamed to `--ignoreblankintiles' to highlight that the
+ status of blank elements can be set in all tiled outputs.
+
+ Statistics:
+ --ignoreblankintiles: similar to same option in NoiseChisel.
+
** Bugs fixed
bug #55313: Fits program writing --write values in reverse order
bug #55333: Fits program crash when --write or --update have no value.
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/bin/noisechisel/args.h b/bin/noisechisel/args.h
index 07cfdb0..d665af3 100644
--- a/bin/noisechisel/args.h
+++ b/bin/noisechisel/args.h
@@ -137,13 +137,13 @@ struct argp_option program_options[] =
/* Output options. */
{
- "ignoreblankinsky",
- UI_KEY_IGNOREBLANKINSKY,
+ "ignoreblankintiles",
+ UI_KEY_IGNOREBLANKINTILES,
0,
0,
- "Don't write input's blanks in the Sky output.",
+ "Don't write input's blanks in tiled output.",
GAL_OPTIONS_GROUP_OUTPUT,
- &p->ignoreblankinsky,
+ &p->ignoreblankintiles,
GAL_OPTIONS_NO_ARG_TYPE,
GAL_OPTIONS_RANGE_0_OR_1,
GAL_OPTIONS_NOT_MANDATORY,
@@ -422,6 +422,19 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET
},
{
+ "pseudoconcomp",
+ UI_KEY_PSEUDOCONCOMP,
+ "INT",
+ 0,
+ "4 or 8 neighbors for labeling pseudo-dets.",
+ UI_GROUP_DETECTION,
+ &p->pseudoconcomp,
+ GAL_TYPE_SIZE_T,
+ GAL_OPTIONS_RANGE_GT_0,
+ GAL_OPTIONS_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+ {
"snminarea",
UI_KEY_SNMINAREA,
"INT",
diff --git a/bin/noisechisel/astnoisechisel.conf
b/bin/noisechisel/astnoisechisel.conf
index 2362820..43c1476 100644
--- a/bin/noisechisel/astnoisechisel.conf
+++ b/bin/noisechisel/astnoisechisel.conf
@@ -40,6 +40,7 @@
sigmaclip 3,0.2
dthresh 0.0
holengb 8
+ pseudoconcomp 8
snminarea 10
minnumfalse 100
snquant 0.99
diff --git a/bin/noisechisel/detection.c b/bin/noisechisel/detection.c
index 48c7345..f53973e 100644
--- a/bin/noisechisel/detection.c
+++ b/bin/noisechisel/detection.c
@@ -361,6 +361,7 @@ detection_pseudo_find(struct noisechiselparams *p,
gal_data_t *workbin,
uint8_t *b, *bf;
gal_data_t *bin;
struct fho_params fho_prm={0, NULL, workbin, worklab, p};
+ int con=detection_ngb_to_connectivity(p->input->ndim, p->pseudoconcomp);
/* Set all the initial detected pixels to blank values. */
@@ -474,7 +475,7 @@ detection_pseudo_find(struct noisechiselparams *p,
gal_data_t *workbin,
do if(*b==GAL_BLANK_UINT8) *b = !s0d1; while(++b<bf);
}
*/
- return gal_binary_connected_components(workbin, &worklab, 1);
+ return gal_binary_connected_components(workbin, &worklab, con);
}
diff --git a/bin/noisechisel/main.h b/bin/noisechisel/main.h
index 07867a8..b38686e 100644
--- a/bin/noisechisel/main.h
+++ b/bin/noisechisel/main.h
@@ -52,7 +52,7 @@ struct noisechiselparams
char *whdu; /* Wide kernel HDU. */
uint8_t continueaftercheck; /* Don't abort after the check steps. */
- uint8_t ignoreblankinsky; /* Ignore input's blank values. */
+ uint8_t ignoreblankintiles; /* Ignore input's blank values. */
uint8_t rawoutput; /* Only detection & 1 elem/tile output. */
uint8_t label; /* Label detections that are connected. */
@@ -74,6 +74,7 @@ struct noisechiselparams
uint8_t checkdetsky; /* Check pseudo-detection sky value. */
float dthresh; /* Sigma threshold for Pseudo-detections. */
size_t holengb; /* Connectivity for defining a hole. */
+ size_t pseudoconcomp; /* Connectivity for connected components. */
size_t snminarea; /* Minimum pseudo-detection area for S/N. */
uint8_t checksn; /* Save pseudo-detection S/N values. */
size_t minnumfalse; /* Min No. of det/seg for true quantile. */
diff --git a/bin/noisechisel/noisechisel.c b/bin/noisechisel/noisechisel.c
index 83c3d30..1227e06 100644
--- a/bin/noisechisel/noisechisel.c
+++ b/bin/noisechisel/noisechisel.c
@@ -184,7 +184,7 @@ noisechisel_output(struct noisechiselparams *p)
/* Write the Sky image into the output */
if(p->sky->name) free(p->sky->name);
p->sky->name="SKY";
- gal_tile_full_values_write(p->sky, &p->cp.tl, !p->ignoreblankinsky,
+ gal_tile_full_values_write(p->sky, &p->cp.tl, !p->ignoreblankintiles,
p->cp.output, NULL, PROGRAM_NAME);
p->sky->name=NULL;
@@ -200,7 +200,7 @@ noisechisel_output(struct noisechiselparams *p)
gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "MEDSTD", 0, &p->medstd, 0,
"Median raw tile standard deviation", 0,
p->input->unit);
- gal_tile_full_values_write(p->std, &p->cp.tl, !p->ignoreblankinsky,
+ gal_tile_full_values_write(p->std, &p->cp.tl, !p->ignoreblankintiles,
p->cp.output, keys, PROGRAM_NAME);
p->std->name=NULL;
diff --git a/bin/noisechisel/sky.c b/bin/noisechisel/sky.c
index bb2d091..78abfa2 100644
--- a/bin/noisechisel/sky.c
+++ b/bin/noisechisel/sky.c
@@ -210,10 +210,10 @@ sky_and_std(struct noisechiselparams *p, char *checkname)
{
p->sky->name="SKY";
p->std->name="STD";
- gal_tile_full_values_write(p->sky, tl, 1, checkname, NULL,
- PROGRAM_NAME);
- gal_tile_full_values_write(p->std, tl, 1, checkname, NULL,
- PROGRAM_NAME);
+ gal_tile_full_values_write(p->sky, tl, !p->ignoreblankintiles,
+ checkname, NULL, PROGRAM_NAME);
+ gal_tile_full_values_write(p->std, tl, !p->ignoreblankintiles,
+ checkname, NULL, PROGRAM_NAME);
p->sky->name=p->std->name=NULL;
}
diff --git a/bin/noisechisel/threshold.c b/bin/noisechisel/threshold.c
index 11a56a6..47dcd19 100644
--- a/bin/noisechisel/threshold.c
+++ b/bin/noisechisel/threshold.c
@@ -300,12 +300,13 @@ threshold_interp_smooth(struct noisechiselparams *p,
gal_data_t **first,
(*first)->name="THRESH1_INTERP";
(*second)->name="THRESH2_INTERP";
if(third) (*third)->name="THRESH3_INTERP";
- gal_tile_full_values_write(*first, tl, 1, filename, NULL, PROGRAM_NAME);
- gal_tile_full_values_write(*second, tl, 1, filename, NULL,
- PROGRAM_NAME);
+ gal_tile_full_values_write(*first, tl, !p->ignoreblankintiles,
+ filename, NULL, PROGRAM_NAME);
+ gal_tile_full_values_write(*second, tl, !p->ignoreblankintiles,
+ filename, NULL, PROGRAM_NAME);
if(third)
- gal_tile_full_values_write(*third, tl, 1, filename, NULL,
- PROGRAM_NAME);
+ gal_tile_full_values_write(*third, tl, !p->ignoreblankintiles,
+ filename, NULL, PROGRAM_NAME);
(*first)->name = (*second)->name = NULL;
if(third) (*third)->name=NULL;
}
@@ -340,13 +341,13 @@ threshold_interp_smooth(struct noisechiselparams *p,
gal_data_t **first,
(*first)->name="THRESH1_SMOOTH";
(*second)->name="THRESH2_SMOOTH";
if(third) (*third)->name="THRESH3_SMOOTH";
- gal_tile_full_values_write(*first, tl, 1, filename, NULL,
- PROGRAM_NAME);
- gal_tile_full_values_write(*second, tl, 1, filename, NULL,
- PROGRAM_NAME);
+ gal_tile_full_values_write(*first, tl, !p->ignoreblankintiles,
+ filename, NULL, PROGRAM_NAME);
+ gal_tile_full_values_write(*second, tl, !p->ignoreblankintiles,
+ filename, NULL, PROGRAM_NAME);
if(third)
- gal_tile_full_values_write(*third, tl, 1, filename, NULL,
- PROGRAM_NAME);
+ gal_tile_full_values_write(*third, tl, !p->ignoreblankintiles,
+ filename, NULL, PROGRAM_NAME);
(*first)->name = (*second)->name = NULL;
if(third) (*third)->name=NULL;
}
@@ -645,17 +646,20 @@ threshold_quantile_find_apply(struct noisechiselparams *p)
{
qprm.erode_th->name="QTHRESH_ERODE";
qprm.noerode_th->name="QTHRESH_NOERODE";
- gal_tile_full_values_write(qprm.erode_th, tl, 1, p->qthreshname, NULL,
- PROGRAM_NAME);
- gal_tile_full_values_write(qprm.noerode_th, tl, 1, p->qthreshname, NULL,
- PROGRAM_NAME);
+ gal_tile_full_values_write(qprm.erode_th, tl,
+ !p->ignoreblankintiles,
+ p->qthreshname, NULL, PROGRAM_NAME);
+ gal_tile_full_values_write(qprm.noerode_th, tl,
+ !p->ignoreblankintiles,
+ p->qthreshname, NULL, PROGRAM_NAME);
qprm.erode_th->name=qprm.noerode_th->name=NULL;
if(qprm.expand_th)
{
qprm.expand_th->name="QTHRESH_EXPAND";
- gal_tile_full_values_write(qprm.expand_th, tl, 1, p->qthreshname,
- NULL, PROGRAM_NAME);
+ gal_tile_full_values_write(qprm.expand_th, tl,
+ !p->ignoreblankintiles,
+ p->qthreshname, NULL, PROGRAM_NAME);
qprm.expand_th->name=NULL;
}
}
diff --git a/bin/noisechisel/ui.c b/bin/noisechisel/ui.c
index 1434dff..84f67b2 100644
--- a/bin/noisechisel/ui.c
+++ b/bin/noisechisel/ui.c
@@ -241,9 +241,10 @@ ui_read_check_only_options(struct noisechiselparams *p)
"for it");
/* A general check on the neighbor connectivity values. */
- ui_ngb_check(p->holengb, "holengb");
- ui_ngb_check(p->erodengb, "erodengb");
+ ui_ngb_check(p->holengb, "holengb");
+ ui_ngb_check(p->erodengb, "erodengb");
ui_ngb_check(p->openingngb, "openingngb");
+ ui_ngb_check(p->pseudoconcomp, "pseudoconcomp");
/* Make sure that the no-erode-quantile is not smaller or equal to
qthresh. */
diff --git a/bin/noisechisel/ui.h b/bin/noisechisel/ui.h
index 7cd0bd4..66f32cc 100644
--- a/bin/noisechisel/ui.h
+++ b/bin/noisechisel/ui.h
@@ -93,13 +93,14 @@ enum option_keys_enum
UI_KEY_SKYFRACNOBLANK,
UI_KEY_CHECKDETSKY,
UI_KEY_HOLENGB,
+ UI_KEY_PSEUDOCONCOMP,
UI_KEY_CHECKSN,
UI_KEY_DETGROWMAXHOLESIZE,
UI_KEY_CLEANGROWNDET,
UI_KEY_CHECKDETECTION,
UI_KEY_CHECKSKY,
UI_KEY_RAWOUTPUT,
- UI_KEY_IGNOREBLANKINSKY,
+ UI_KEY_IGNOREBLANKINTILES,
};
diff --git a/bin/statistics/args.h b/bin/statistics/args.h
index 1684c15..7fbd64b 100644
--- a/bin/statistics/args.h
+++ b/bin/statistics/args.h
@@ -548,13 +548,13 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET
},
{
- "ignoreblankinsky",
- UI_KEY_IGNOREBLANKINSKY,
+ "ignoreblankintiles",
+ UI_KEY_IGNOREBLANKINTILES,
0,
0,
- "Don't write input's blanks in the Sky output.",
+ "Don't write input's blanks in the tiled output.",
UI_GROUP_SKY,
- &p->ignoreblankinsky,
+ &p->ignoreblankintiles,
GAL_OPTIONS_NO_ARG_TYPE,
GAL_OPTIONS_RANGE_0_OR_1,
GAL_OPTIONS_NOT_MANDATORY,
diff --git a/bin/statistics/main.h b/bin/statistics/main.h
index 01ae40c..7e835cc 100644
--- a/bin/statistics/main.h
+++ b/bin/statistics/main.h
@@ -92,7 +92,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. */
- uint8_t ignoreblankinsky; /* Ignore input's blank values. */
+ uint8_t ignoreblankintiles;/* Ignore input's blank values. */
/* Internal */
diff --git a/bin/statistics/sky.c b/bin/statistics/sky.c
index 11b6b11..6d6a3fe 100644
--- a/bin/statistics/sky.c
+++ b/bin/statistics/sky.c
@@ -203,10 +203,10 @@ sky(struct statisticsparams *p)
}
if(p->checksky)
{
- gal_tile_full_values_write(p->sky_t, tl, 1, p->checkskyname, NULL,
- PROGRAM_NAME);
- gal_tile_full_values_write(p->std_t, tl, 1, p->checkskyname, NULL,
- PROGRAM_NAME);
+ gal_tile_full_values_write(p->sky_t, tl, !p->ignoreblankintiles,
+ p->checkskyname, NULL, PROGRAM_NAME);
+ gal_tile_full_values_write(p->std_t, tl, !p->ignoreblankintiles,
+ p->checkskyname, NULL, PROGRAM_NAME);
}
@@ -231,9 +231,9 @@ sky(struct statisticsparams *p)
gal_timing_report(&t1, "All blank tiles filled (interplated).", 1);
if(p->checksky)
{
- gal_tile_full_values_write(p->sky_t, tl, !p->ignoreblankinsky,
+ gal_tile_full_values_write(p->sky_t, tl, !p->ignoreblankintiles,
p->checkskyname, NULL, PROGRAM_NAME);
- gal_tile_full_values_write(p->std_t, tl, !p->ignoreblankinsky,
+ gal_tile_full_values_write(p->std_t, tl, !p->ignoreblankintiles,
p->checkskyname, NULL, PROGRAM_NAME);
}
@@ -255,9 +255,9 @@ sky(struct statisticsparams *p)
1);
if(p->checksky)
{
- gal_tile_full_values_write(p->sky_t, tl, !p->ignoreblankinsky,
+ gal_tile_full_values_write(p->sky_t, tl, !p->ignoreblankintiles,
p->checkskyname, NULL, PROGRAM_NAME);
- gal_tile_full_values_write(p->std_t, tl, !p->ignoreblankinsky,
+ gal_tile_full_values_write(p->std_t, tl, !p->ignoreblankintiles,
p->checkskyname, NULL, PROGRAM_NAME);
if(!cp->quiet)
printf(" - Check image written to `%s'.\n", p->checkskyname);
@@ -277,9 +277,9 @@ sky(struct statisticsparams *p)
p->sky_t->name="SKY";
p->std_t->name="SKY_STD";
p->cp.keepinputdir=keepinputdir;
- gal_tile_full_values_write(p->sky_t, tl, !p->ignoreblankinsky, outname,
+ gal_tile_full_values_write(p->sky_t, tl, !p->ignoreblankintiles, outname,
NULL, PROGRAM_NAME);
- gal_tile_full_values_write(p->std_t, tl, !p->ignoreblankinsky, outname,
+ gal_tile_full_values_write(p->std_t, tl, !p->ignoreblankintiles, outname,
NULL, PROGRAM_NAME);
p->sky_t->name = p->std_t->name = NULL;
gal_fits_key_write_filename("input", p->inputname, &p->cp.okeys, 1);
diff --git a/bin/statistics/statistics.c b/bin/statistics/statistics.c
index 606fbb1..c92ca8b 100644
--- a/bin/statistics/statistics.c
+++ b/bin/statistics/statistics.c
@@ -259,7 +259,8 @@ statistics_interpolate_and_write(struct statisticsparams *p,
}
/* Write the values. */
- gal_tile_full_values_write(values, &cp->tl, 1, output, NULL, PROGRAM_NAME);
+ gal_tile_full_values_write(values, &cp->tl, !p->ignoreblankintiles,
+ output, NULL, PROGRAM_NAME);
gal_fits_key_write_filename("input", p->inputname, &p->cp.okeys, 1);
gal_fits_key_write_config(&p->cp.okeys, "Statistics configuration",
"STATISTICS-CONFIG", output, "0");
diff --git a/bin/statistics/ui.h b/bin/statistics/ui.h
index de10940..de4782c 100644
--- a/bin/statistics/ui.h
+++ b/bin/statistics/ui.h
@@ -98,7 +98,7 @@ enum option_keys_enum
UI_KEY_OUTLIERSCLIP,
UI_KEY_SMOOTHWIDTH,
UI_KEY_CHECKSKY,
- UI_KEY_IGNOREBLANKINSKY,
+ UI_KEY_IGNOREBLANKINTILES,
UI_KEY_SCLIPPARAMS,
};
diff --git a/configure.ac b/configure.ac
index 8fea1e7..9318393 100644
--- a/configure.ac
+++ b/configure.ac
@@ -420,16 +420,18 @@ AS_IF([test "x$has_gnulibtool" = "xyes"],
echo "#include <stdio.h>" > $cprog
echo "int main(void){printf(\"success\\n\"); return 0;}" >> $cprog
ltargs="--quiet --tag=CC --mode=link $CC $cprog -O3 -o $outname"
- AS_IF( sh -c "$gnulibtool_exec $ltargs" > /dev/null 2>&1,
- [libtool_shell="sh"],
- [AS_IF([test "x$has_bash" = "xyes"],
- [AS_IF(bash -c "$gnulibtool_exec $ltargs" > /dev/null
2>&1,
- [libtool_shell="bash"]) ],
- [AS_IF([test "x$has_zsh" = "xyes"],
- [AS_IF(zsh -c "$gnulibtool_exec $ltargs" >
/dev/null 2>&1,
- [libtool_shell="zsh"]) ])
- ] )
- ] )
+
+ # Check the shells, starting with known shells and ultimately
+ # trying with `sh' (can be any shell).
+ AS_IF([test "x$has_bash" = "xyes"],
+ [AS_IF(bash -c "$gnulibtool_exec $ltargs" > /dev/null 2>&1,
+ [libtool_shell="bash"],
+ [bash_version=$BASH_VERSION]) ],
+ [AS_IF([test "x$has_zsh" = "xyes"],
+ [AS_IF(zsh -c "$gnulibtool_exec $ltargs" > /dev/null 2>&1,
+ [libtool_shell="zsh"]) ],
+ [AS_IF(sh -c "$gnulibtool_exec $ltargs" > /dev/null 2>&1,
+ [libtool_shell="sh"]) ]) ])
# Clean up: note that no output might have been generated (when no
# proper shell was found). Therefore, for deleting the output file,
@@ -443,7 +445,12 @@ AS_IF([test "x$has_gnulibtool" = "xyes"],
# If a good shell to call Libtool could be found, then Libtool is usable
# within a program (BuildProgram in this case).
AS_IF([test "x$libtool_shell" = "xnone"],
- [usable_libtool=no; anywarnings=yes],
+ [
+ anywarnings=yes;
+ usable_libtool=no;
+ AS_IF([test "x$bash_version" = x], [junk=1],
+ [echo "GNU Bash was found but not a suitable version:
$bash_version"])
+ ],
[
usable_libtool=yes
AC_DEFINE_UNQUOTED([GAL_CONFIG_GNULIBTOOL_SHELL], ["$libtool_shell"],
@@ -455,6 +462,7 @@ AS_IF([test "x$libtool_shell" = "xnone"],
+
# Check Ghostscript: "-dPDFFitPage" option to Ghostscript, used by the
# library to convert from EPS to PDF, has been introduced in Ghostscript
# 9.10. Make sure we have at least that version.
@@ -902,6 +910,8 @@ AS_IF([test x$enable_guide_message = xyes],
AS_ECHO([])
AS_IF([test "x$has_gnulibtool" = "xyes"],
[AS_ECHO([" -- GNU Libtool is present, but couldn't be
run in tested shells."])
+ AS_IF([test "x$bash_version" = x], [junk=1],
+ [echo " (GNU Bash was found but not a
suitable version: $bash_version)"])
AS_ECHO([])],
[AS_IF([test "x$has_libtool" = "xyes"],
[AS_ECHO([" -- A libtool implementation was
found, but it isn't GNU."])
diff --git a/doc/announce-acknowledge.txt b/doc/announce-acknowledge.txt
index d690df8..e5c7f79 100644
--- a/doc/announce-acknowledge.txt
+++ b/doc/announce-acknowledge.txt
@@ -1,3 +1,5 @@
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 f88e3b8..ba09486 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -12178,7 +12178,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}.
-@item numvalue
+@item numbervalue
Number of non-blank elements in first operand in the @code{uint64} type,
similar to @command{minvalue}.
@@ -12218,8 +12218,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
@@ -12235,7 +12234,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.
-@item num
+@item 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.
@@ -12255,6 +12254,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.
+@item 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}).
+
+@example
+astarithmetic a.fits b.fits c.fits 3 5 0.2 sigclip-number
+@end example
+
+@item 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}.
+
+@item 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}.
+
+@item 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
@@ -16179,13 +16215,17 @@ discontinuities do not show up in the final Sky
values. The smoothing is
done through convolution (see @ref{Convolution process}) with a flat
kernel, so the value to this option must be an odd number.
-@item --ignoreblankinsky
-Don't set the input's blank pixels to blank in the output Sky and Sky
-standard deviation datasets. This is only applicable when the output has
-the same size as the input, in other words, when @option{--oneelempertile}
-isn't called. By default, blank values in the input (commonly on the edges
-which are outside the survey/field area) will be set to blank in the output
-Sky and Sky standard deviation also.
+@item --ignoreblankintiles
+Don't set the input's blank pixels to blank in the tiled outputs (for
+example Sky and Sky standard deviation extensions of the output). This is
+only applicable when the tiled output has the same size as the input, in
+other words, when @option{--oneelempertile} isn't called.
+
+By default, blank values in the input (commonly on the edges which are
+outside the survey/field area) will be set to blank in the tiled outputs
+also. But in other scenarios this default behavior is not desired: for
+example if you have masked something in the input, but want the tiled
+output under that also.
@item --checksky
Create a multi-extension FITS file showing the steps that were used to
@@ -16445,6 +16485,17 @@ the image. For more, see the description of this
option in @ref{Detection
options}.
@item
+@option{--holengb}: The connectivity (defined by the number of neighbors)
+to fill holes after applying @option{--dthresh} (above) to find
+pseudo-detections. For more, see the description of this option in
+@ref{Detection options}.
+
+@item
+@option{--pseudoconcomp}: The connectivity (defined by the number of
+neighbors) to find individual pseudo-detections. For more, see the
+description of this option in @ref{Detection options}.
+
+@item
@option{--detgrowquant}: is used to grow the final true detections until a
given quantile in the same way that clumps are grown during segmentation
(compare columns 2 and 3 in Figure 10 of the paper). It replaces the old
@@ -17035,6 +17086,12 @@ walls of the hole are 4-connected. If standard (near
Sky level) values are
given to @option{--dthresh}, setting @option{--holengb=4}, might fill the
complete dataset and thus not create enough pseudo-detections.
+@item --pseudoconcomp=INT
+The connectivity (defined by the number of neighbors) to find individual
+pseudo-detections. If it is a weaker connectivity (4 in a 2D image), then
+pseudo-detections that are connected on the corners will be treated as
+separate.
+
@item -m INT
@itemx --snminarea=INT
The minimum area to calculate the Signal to noise ratio on the
@@ -17235,13 +17292,17 @@ NoiseChisel will abort once its desired extensions
have been written. With
NoiseChisel to continue with the rest of the processing, even after the
requested check files are complete.
-@item --ignoreblankinsky
-Don't set the input's blank pixels to blank in the output Sky and Sky
-standard deviation datasets. This is only applicable when the output has
-the same size as the input, in other words, when @option{--oneelempertile}
-isn't called. By default, blank values in the input (commonly on the edges
-which are outside the survey/field area) will be set to blank in the output
-Sky and Sky standard deviation also.
+@item --ignoreblankintiles
+Don't set the input's blank pixels to blank in the tiled outputs (for
+example Sky and Sky standard deviation extensions of the output). This is
+only applicable when the tiled output has the same size as the input, in
+other words, when @option{--oneelempertile} isn't called.
+
+By default, blank values in the input (commonly on the edges which are
+outside the survey/field area) will be set to blank in the tiled outputs
+also. But in other scenarios this default behavior is not desired: for
+example if you have masked something in the input, but want the tiled
+output under that also.
@item -l
@itemx --label
@@ -27223,7 +27284,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
-@deffnx Macro GAL_ARITHMETIC_OP_NUMVAL
+@deffnx Macro GAL_ARITHMETIC_OP_NUMBERVAL
@deffnx Macro GAL_ARITHMETIC_OP_SUMVAL
@deffnx Macro GAL_ARITHMETIC_OP_MEANVAL
@deffnx Macro GAL_ARITHMETIC_OP_STDVAL
@@ -27242,7 +27303,7 @@ Unary operand absolute-value operator.
@deffn Macro GAL_ARITHMETIC_OP_MIN
@deffnx Macro GAL_ARITHMETIC_OP_MAX
-@deffnx Macro GAL_ARITHMETIC_OP_NUM
+@deffnx Macro GAL_ARITHMETIC_OP_NUMBER
@deffnx Macro GAL_ARITHMETIC_OP_SUM
@deffnx Macro GAL_ARITHMETIC_OP_MEAN
@deffnx Macro GAL_ARITHMETIC_OP_STD
@@ -27255,6 +27316,18 @@ respective statistical operation on the whole list.
See the discussion
under the @code{min} operator in @ref{Arithmetic operators}.
@end deffn
+@deffn Macro GAL_ARITHMETIC_OP_SIGCLIP_STD
+@deffnx Macro GAL_ARITHMETIC_OP_SIGCLIP_MEAN
+@deffnx Macro GAL_ARITHMETIC_OP_SIGCLIP_MEDIAN
+@deffnx 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}).
+@end 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
@@ -29252,8 +29325,8 @@ The next major difference with over-segmentation is
that when there is only
one label in growth region(s), it is not mandatory for @code{indexs} to be
sorted by values. If there are multiple labeled regions in growth
region(s), then values are important and you can use @code{qsort} with
-@code{gal_qsort_index_float_decreasing} to sort the indexs by values in a
-separate array (see @ref{Qsort functions}).
+@code{gal_qsort_index_single_d} to sort the indexs by values in a separate
+array (see @ref{Qsort functions}).
This function looks for positive-valued neighbors of each pixel in
@code{indexs} and will label a pixel if it touches one. Therefore, it is
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/box.c b/lib/box.c
index 035bc14..1a48f71 100644
--- a/lib/box.c
+++ b/lib/box.c
@@ -315,7 +315,7 @@ gal_box_border_from_center(double *center, size_t ndim,
long *width,
input image (fpixel_i[2] and lpixel_i[2]). But those first and last
pixels don't necessarily lie within the image's boundaries. They can be
outside of it or patially overlap with it (see examples below). The job
- of this function is to corret for such situations and find the starting
+ of this function is to correct for such situations and find the starting
and ending points of any overlap.
It is assumed that your output (overlap) image's first pixel lies right
@@ -419,6 +419,11 @@ gal_box_overlap(long *naxes, long *fpixel_i, long
*lpixel_i,
*/
if(fpixel_i[i]<1)
{
+ /* Along any dimension, if `lpixel_i' is also smaller than 1,
+ then there is no overlap. */
+ if(lpixel_i[i]<1) return 0;
+
+ /* Correct the coordinates. */
fpixel_o[i] = -1*fpixel_i[i]+2;
fpixel_i[i] = 1;
}
@@ -437,6 +442,11 @@ gal_box_overlap(long *naxes, long *fpixel_i, long
*lpixel_i,
cropped image we should only fill upto c-n.*/
if(lpixel_i[i]>naxes[i])
{
+ /* Along any dimension, if `fpixel_i' is larger than the image
+ size, there is no overlap. */
+ if(fpixel_i[i]>naxes[i]) return 0;
+
+ /* Correct the coordinates. */
lpixel_o[i] = width - (lpixel_i[i]-naxes[i]);
lpixel_i[i] = naxes[i];
}
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. */
- [gnuastro-commits] master 9e2397a 094/113: Imported recent work in master, no conflicts, (continued)
- [gnuastro-commits] master 9e2397a 094/113: Imported recent work in master, no conflicts, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 9258f68 096/113: Imported recent work in master, minor conflict fixed, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master d0d8d20 109/113: Imported recent work in master, minor conflicts fixed, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 6e11585 111/113: Imported recent work in master, minor conflicts fixed, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master aa80ac4 112/113: Imported recent updates in master branch, minor conflict fixed, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master dadb3f3 107/113: Imported recent work in master, minor conflicts fixed, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 313d522 099/113: Imported recent work in master, minor conflict fixed, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 536b056 110/113: Imported recent changes in master, conflicts in book fixed, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master dd4d43e 113/113: NoiseChisel and Segment: now working on 3D data cubes, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master eb91303 042/113: Merged recent work in master, conflicts fixed, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 2dba72f 101/113: Imported work in master, conflicts fixed, changes made,
Mohammad Akhlaghi <=
- [gnuastro-commits] master b3416de 102/113: Imported recent work in master, conflicts fixed, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 353c0c1 104/113: Imported recent work in master, conflict fixed, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master f21c30f 105/113: Imported recent work in master, no conflicts, Mohammad Akhlaghi, 2021/04/16