[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 3b26d5c 114/125: SubtractSky program removed,
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 3b26d5c 114/125: SubtractSky program removed, Statistics estimates Sky |
|
Date: |
Sun, 23 Apr 2017 22:36:50 -0400 (EDT) |
branch: master
commit 3b26d5cec7522523cd8fc40283f7575d28210320
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
SubtractSky program removed, Statistics estimates Sky
Thanks to Arithmetic, the `Subtract' in the SubtractSky program was useless
for a long time, but it still remained. So after the old SubtractSky was
re-written to work with all the new infra-structure, it was decided that
ultimately, finding the Sky value (and its standard deviation) over an
image is a statistical operation. As a result, the estimation of the Sky
value has been added as another operation for the Statistics program. The
book has also been corrected to accommodate all these changes.
In the process of merging the new infra-structure into the sky finding
function, several other issues were addressed
- The default tile size was changed to 50 from 40. This was done because
50 is a more common number (people might think there is a special reason
for choosing 40 and get confused). Also, it allows for more pixels, thus
getting better estimates of the measurements on the tiles.
- Tile interpolation options have been moved to the common arguments,
since they are also needed in NoiseChisel.
- The new `gal_options_read_sigma_clip' function will parse the values
given to sigmaclipping parameters.
- When the mode symmetricity is not good (mode is not accurate),
`gal_statistics_mode' will return an array of all NaN values.
- `gal_statistics_sigma_clip' now also uses
`gal_statistics_no_blank_sorted'. Before it had its own implementation
which didn't account for tiles.
- Several new functions in `lib/tile.c'.
- The main Gnuastro general configuration file `gnuastro.conf' was moved
to the `bin' directory from the `lib' directory. The files in the `lib'
directory are meant to be for the library. So this file was out of place
there. It is now managed by the top level `Makefile.am'.
---
Makefile.am | 13 +-
bin/convolve/ui.c | 2 +
{lib => bin}/gnuastro.conf | 10 +-
bin/statistics/Makefile.am | 4 +-
bin/statistics/args.h | 212 ++-
bin/statistics/aststatistics.conf | 8 +-
bin/statistics/main.h | 19 +-
bin/statistics/sky.c | 272 ++++
.../subtractsky.h => statistics/sky.h} | 10 +-
bin/statistics/statistics.c | 89 +-
bin/statistics/ui.c | 268 ++--
bin/statistics/ui.h | 11 +-
bin/subtractsky/Makefile.am | 40 -
bin/subtractsky/args.h | 505 -------
bin/subtractsky/astsubtractsky.conf | 40 -
bin/subtractsky/authors-cite.h | 38 -
bin/subtractsky/main.c | 56 -
bin/subtractsky/main.h | 100 --
bin/subtractsky/subtractsky.c | 247 ----
bin/subtractsky/ui.c | 627 ---------
bin/subtractsky/ui.h | 32 -
configure.ac | 10 -
doc/Makefile.am | 9 +-
doc/gnuastro.texi | 1400 +++++++++-----------
lib/Makefile.am | 20 +-
lib/commonopts.h | 28 +-
lib/convolve.c | 15 +-
lib/data.c | 76 +-
lib/fits.c | 49 +-
lib/gnuastro/data.h | 13 +-
lib/gnuastro/interpolate.h | 2 +-
lib/gnuastro/linkedlist.h | 25 +
lib/gnuastro/statistics.h | 2 +-
lib/gnuastro/tile.h | 15 +-
lib/interpolate.c | 182 ++-
lib/linkedlist.c | 90 ++
lib/options.c | 75 ++
lib/options.h | 8 +-
lib/statistics.c | 123 +-
lib/tile.c | 115 +-
tests/Makefile.am | 10 +-
tests/during-dev.sh | 4 +-
tests/prepconf.sh | 4 +-
.../subtractsky.sh => statistics/estimate_sky.sh} | 11 +-
tmpfs-config-make | 12 +-
45 files changed, 1917 insertions(+), 2984 deletions(-)
diff --git a/Makefile.am b/Makefile.am
index 6f3338a..65f78db 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -80,9 +80,6 @@ endif
if COND_STATISTICS
MAYBE_STATISTICS = bin/statistics
endif
-if COND_SUBTRACTSKY
- MAYBE_SUBTRACTSKY = bin/subtractsky
-endif
if COND_TABLE
MAYBE_TABLE = bin/table
endif
@@ -110,8 +107,8 @@ endif
SUBDIRS = bootstrapped/lib $(MAYBE_GNULIBCHECK) lib $(MAYBE_ARITHMETIC) \
$(MAYBE_CONVERTT) $(MAYBE_CONVOLVE) $(MAYBE_COSMICCAL) $(MAYBE_CROP) \
$(MAYBE_FITS) $(MAYBE_MKCATALOG) $(MAYBE_MKNOISE) $(MAYBE_MKPROF) \
- $(MAYBE_NOISECHISEL) $(MAYBE_STATISTICS) $(MAYBE_SUBTRACTSKY) \
- $(MAYBE_TABLE) $(MAYBE_TEMPLATE) $(MAYBE_WARP) doc tests
+ $(MAYBE_NOISECHISEL) $(MAYBE_STATISTICS) $(MAYBE_TABLE) \
+ $(MAYBE_TEMPLATE) $(MAYBE_WARP) doc tests
@@ -127,6 +124,12 @@ dist_doc_DATA = README
+# Installed system configuration files
+# ====================================
+dist_sysconf_DATA = bin/gnuastro.conf
+
+
+
## Files that are only distributed
## ===============================
EXTRA_DIST = bootstrap bootstrap.conf genauthors .dir-locals.el .version \
diff --git a/bin/convolve/ui.c b/bin/convolve/ui.c
index 4b1225d..866e2c9 100644
--- a/bin/convolve/ui.c
+++ b/bin/convolve/ui.c
@@ -145,6 +145,8 @@ ui_initialize_options(struct convolveparams *p,
case GAL_OPTIONS_KEY_SEARCHIN:
case GAL_OPTIONS_KEY_IGNORECASE:
case GAL_OPTIONS_KEY_TABLEFORMAT:
+ case GAL_OPTIONS_KEY_INTERPNUMNGB:
+ case GAL_OPTIONS_KEY_INTERPONLYBLANK:
cp->coptions[i].flags=OPTION_HIDDEN;
break;
}
diff --git a/lib/gnuastro.conf b/bin/gnuastro.conf
similarity index 84%
rename from lib/gnuastro.conf
rename to bin/gnuastro.conf
index ab2e5c6..d64529f 100644
--- a/lib/gnuastro.conf
+++ b/bin/gnuastro.conf
@@ -23,10 +23,12 @@
searchin name
# Tessellation
- tilesize 40,40
- numchannels 1,1
- remainderfrac 0.1
- workoverch 0
+ tilesize 50,50
+ numchannels 1,1
+ remainderfrac 0.1
+ workoverch 0
+ interpnumngb 9
+ interponlyblank 0
# Output:
tableformat fits-binary
diff --git a/bin/statistics/Makefile.am b/bin/statistics/Makefile.am
index e19778a..6454a24 100644
--- a/bin/statistics/Makefile.am
+++ b/bin/statistics/Makefile.am
@@ -30,9 +30,9 @@ bin_PROGRAMS = aststatistics
aststatistics_LDADD = -lgnuastro
-aststatistics_SOURCES = main.c ui.c statistics.c
+aststatistics_SOURCES = main.c ui.c sky.c statistics.c
-EXTRA_DIST = main.h authors-cite.h args.h ui.h statistics.h
+EXTRA_DIST = main.h authors-cite.h args.h ui.h sky.h statistics.h
diff --git a/bin/statistics/args.h b/bin/statistics/args.h
index 01b8b72..9617ddc 100644
--- a/bin/statistics/args.h
+++ b/bin/statistics/args.h
@@ -95,7 +95,7 @@ struct argp_option program_options[] =
GAL_OPTIONS_RANGE_ANY,
GAL_OPTIONS_NOT_MANDATORY,
GAL_OPTIONS_NOT_SET,
- ui_parse_numbers
+ ui_read_quantile_range
},
@@ -114,32 +114,6 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_MANDATORY,
GAL_OPTIONS_NOT_SET
},
- {
- "interponlyblank",
- ARGS_OPTION_KEY_INTERPONLYBLANK,
- 0,
- 0,
- "Only interpolate over the blank values.",
- GAL_OPTIONS_GROUP_TESSELLATION,
- &p->interponlyblank,
- GAL_OPTIONS_NO_ARG_TYPE,
- GAL_OPTIONS_RANGE_0_OR_1,
- GAL_OPTIONS_NOT_MANDATORY,
- GAL_OPTIONS_NOT_SET
- },
- {
- "interpnumngb",
- ARGS_OPTION_KEY_INTERPNUMNGB,
- 0,
- 0,
- "Number of neighbors to use for interpolation.",
- GAL_OPTIONS_GROUP_TESSELLATION,
- &p->interpnumngb,
- GAL_OPTIONS_NO_ARG_TYPE,
- GAL_OPTIONS_RANGE_GT_0,
- GAL_OPTIONS_NOT_MANDATORY,
- GAL_OPTIONS_NOT_SET
- },
@@ -280,7 +254,7 @@ struct argp_option program_options[] =
ARGS_OPTION_KEY_MODE,
0,
0,
- "Mode (Appendix C of arXiv:1505.011664).",
+ "Mode (Appendix C of arXiv:1505.01664).",
ARGS_GROUP_SINGLE_VALUE,
&p->singlevalue,
GAL_OPTIONS_NO_ARG_TYPE,
@@ -331,19 +305,6 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET,
ui_add_to_single_value
},
- {
- "ontile",
- ARGS_OPTION_KEY_ONTILE,
- 0,
- 0,
- "Do the measurements on separate tiles, not all.",
- ARGS_GROUP_SINGLE_VALUE,
- &p->ontile,
- GAL_OPTIONS_NO_ARG_TYPE,
- GAL_OPTIONS_RANGE_0_OR_1,
- GAL_OPTIONS_NOT_MANDATORY,
- GAL_OPTIONS_NOT_SET
- },
@@ -406,40 +367,166 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET
},
{
+ "mirror",
+ ARGS_OPTION_KEY_MIRROR,
+ "FLT",
+ 0,
+ "Save the histogram and CFP of the mirror dist.",
+ ARGS_GROUP_PARTICULAR_STAT,
+ &p->mirror,
+ GAL_DATA_TYPE_FLOAT64,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+ {
+ "ontile",
+ ARGS_OPTION_KEY_ONTILE,
+ 0,
+ 0,
+ "Single values on separate tiles, not full input.",
+ ARGS_GROUP_PARTICULAR_STAT,
+ &p->ontile,
+ GAL_OPTIONS_NO_ARG_TYPE,
+ GAL_OPTIONS_RANGE_0_OR_1,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+ {
+ "sky",
+ ARGS_OPTION_KEY_SKY,
+ 0,
+ 0,
+ "Find the Sky and its STD over the tessellation.",
+ ARGS_GROUP_PARTICULAR_STAT,
+ &p->sky,
+ GAL_OPTIONS_NO_ARG_TYPE,
+ GAL_OPTIONS_RANGE_0_OR_1,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+ {
"sigmaclip",
ARGS_OPTION_KEY_SIGMACLIP,
- "FLT,FLT",
0,
- "Multiple of sigma and tolerance/number.",
+ 0,
+ "Overall sigma-clipping (see `--sclipparams')",
ARGS_GROUP_PARTICULAR_STAT,
- NULL,
+ &p->sigmaclip,
+ GAL_OPTIONS_NO_ARG_TYPE,
+ GAL_OPTIONS_RANGE_0_OR_1,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+
+
+
+
+
+ {
+ 0, 0, 0, 0,
+ "Sky and Sky STD settings",
+ ARGS_GROUP_SKY
+ },
+ {
+ "kernel",
+ ARGS_OPTION_KEY_KERNEL,
+ "STR",
+ 0,
+ "File name of kernel to convolve input.",
+ ARGS_GROUP_SKY,
+ &p->kernelname,
GAL_DATA_TYPE_STRING,
GAL_OPTIONS_RANGE_ANY,
GAL_OPTIONS_NOT_MANDATORY,
- GAL_OPTIONS_NOT_SET,
- ui_parse_numbers
+ GAL_OPTIONS_NOT_SET
},
{
- "mirror",
- ARGS_OPTION_KEY_MIRROR,
+ "khdu",
+ ARGS_OPTION_KEY_KHDU,
+ "STR",
+ 0,
+ "HDU/extension name or number of kernel.",
+ ARGS_GROUP_SKY,
+ &p->khdu,
+ GAL_DATA_TYPE_STRING,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+ {
+ "mirrordist",
+ ARGS_OPTION_KEY_MIRRORDIST,
"FLT",
0,
- "Save the histogram and CFP of the mirror dist.",
- ARGS_GROUP_PARTICULAR_STAT,
- &p->mirror,
- GAL_DATA_TYPE_FLOAT64,
+ "Max. distance (error multip.) to find mode.",
+ ARGS_GROUP_SKY,
+ &p->mirrordist,
+ GAL_DATA_TYPE_FLOAT32,
GAL_OPTIONS_RANGE_ANY,
GAL_OPTIONS_NOT_MANDATORY,
GAL_OPTIONS_NOT_SET
},
-
+ {
+ "modmedqdiff",
+ ARGS_OPTION_KEY_MODMEDQDIFF,
+ "FLT",
+ 0,
+ "Max. mode and median quantile diff. per tile.",
+ ARGS_GROUP_SKY,
+ &p->modmedqdiff,
+ GAL_DATA_TYPE_FLOAT32,
+ GAL_OPTIONS_RANGE_GE_0,
+ GAL_OPTIONS_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+ {
+ "sclipparams",
+ ARGS_OPTION_KEY_SCLIPPARAMS,
+ "FLT,FLT",
+ 0,
+ "Multiple of sigma and tolerance/number.",
+ ARGS_GROUP_SKY,
+ p->sclipparams,
+ GAL_DATA_TYPE_STRING,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ gal_options_read_sigma_clip
+ },
+ {
+ "smoothwidth",
+ ARGS_OPTION_KEY_SMOOTHWIDTH,
+ "INT",
+ 0,
+ "Sky: flat kernel width to smooth interpolated.",
+ ARGS_GROUP_SKY,
+ &p->smoothwidth,
+ GAL_DATA_TYPE_SIZE_T,
+ GAL_OPTIONS_RANGE_0_OR_ODD,
+ GAL_OPTIONS_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+ {
+ "checksky",
+ ARGS_OPTION_KEY_CHECKSKY,
+ 0,
+ 0,
+ "Store steps in `_sky_steps.fits' file.",
+ ARGS_GROUP_SKY,
+ &p->checksky,
+ GAL_OPTIONS_NO_ARG_TYPE,
+ GAL_OPTIONS_RANGE_0_OR_1,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
{
0, 0, 0, 0,
- "Settings",
+ "Histogram and CFP settings",
ARGS_GROUP_HIST_CFP
},
{
@@ -447,7 +534,7 @@ struct argp_option program_options[] =
ARGS_OPTION_KEY_NUMBINS,
"INT",
0,
- "Number of bins in histogram or CFP.",
+ "No. of bins in histogram or CFP tables.",
ARGS_GROUP_HIST_CFP,
&p->numbins,
GAL_DATA_TYPE_SIZE_T,
@@ -460,7 +547,7 @@ struct argp_option program_options[] =
ARGS_OPTION_KEY_NUMASCIIBINS,
"INT",
0,
- "Number of bins in ASCII histogram or CFP plots.",
+ "No. of bins in ASCII histogram or CFP plots.",
ARGS_GROUP_HIST_CFP,
&p->numasciibins,
GAL_DATA_TYPE_SIZE_T,
@@ -520,19 +607,6 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_MANDATORY,
GAL_OPTIONS_NOT_SET
},
- {
- "mirrordist",
- ARGS_OPTION_KEY_MIRRORDIST,
- "FLT",
- 0,
- "Maximum dist. (in multip. of error) to find mode.",
- ARGS_GROUP_HIST_CFP,
- &p->mirrordist,
- GAL_DATA_TYPE_FLOAT32,
- GAL_OPTIONS_RANGE_ANY,
- GAL_OPTIONS_NOT_MANDATORY,
- GAL_OPTIONS_NOT_SET
- },
{0}
diff --git a/bin/statistics/aststatistics.conf
b/bin/statistics/aststatistics.conf
index 32db56a..d2ad1bb 100644
--- a/bin/statistics/aststatistics.conf
+++ b/bin/statistics/aststatistics.conf
@@ -19,9 +19,11 @@
# Input image:
-# Tessellation
- interponlyblank 1
- interpnumngb 9
+# Sky and its STD settings
+ khdu 1
+ modmedqdiff 0.01
+ smoothwidth 3
+ sclipparams 3,0.1
# Histogram and CFP settings
numasciibins 70
diff --git a/bin/statistics/main.h b/bin/statistics/main.h
index 6e8cfd8..057541b 100644
--- a/bin/statistics/main.h
+++ b/bin/statistics/main.h
@@ -66,16 +66,14 @@ struct statisticsparams
float quantmax; /* Quantile maximum. */
uint8_t ontile; /* Do single value calculations on tiles. */
uint8_t interpolate; /* Use interpolation to fill blank tiles. */
- uint8_t interponlyblank; /* Only interpolate over blank values. */
- size_t interpnumngb; /* Number of neighbors for interpolation. */
uint8_t asciihist; /* Print an ASCII histogram. */
uint8_t asciicfp; /* Print an ASCII cumulative frequency plot.*/
uint8_t histogram; /* Save histogram in output. */
uint8_t cumulative; /* Save cumulative distibution in output. */
- float sigclipmultip; /* Multiple of sigma in sigma clipping. */
- float sigclipparam; /* Tolerance to stop or number of clips. */
double mirror; /* Mirror value for hist and CFP. */
+ uint8_t sky; /* Find the Sky value over the image. */
+ uint8_t sigmaclip; /* So sigma-clipping over all dataset. */
size_t numbins; /* Number of bins in histogram or CFP. */
size_t numasciibins; /* Number of bins in ASCII plots. */
@@ -85,6 +83,14 @@ struct statisticsparams
uint8_t maxbinone; /* Set the maximum bin to 1. */
float mirrordist; /* Maximum distance after mirror for mode. */
+ char *kernelname; /* File name of kernel to convolve input. */
+ char *khdu; /* Kernel HDU. */
+ float modmedqdiff; /* Mode and median quantile difference. */
+ 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. */
+
+
/* Internal */
uint8_t inputformat; /* Format of input dataset. */
int numoutfiles; /* Number of output files made in this run. */
@@ -94,6 +100,11 @@ struct statisticsparams
gal_data_t *reference; /* Reference data structure. */
int isfits; /* Input is a FITS file. */
int hdu_type; /* Type of HDU (image or table). */
+ gal_data_t *kernel; /* Kernel for convolution of input for Sky. */
+ gal_data_t *convolved; /* Convolved input. */
+ gal_data_t *sky_t; /* Sky on each tile. */
+ gal_data_t *std_t; /* Sky standard deviation on each tile. */
+ char *checkskyname; /* Name of file for Sky calculation steps. */
time_t rawtime; /* Starting time of the program. */
};
diff --git a/bin/statistics/sky.c b/bin/statistics/sky.c
new file mode 100644
index 0000000..79c56b0
--- /dev/null
+++ b/bin/statistics/sky.c
@@ -0,0 +1,272 @@
+/*********************************************************************
+Statistics - Statistical analysis on input dataset.
+Statistics is part of GNU Astronomy Utilities (Gnuastro) package.
+
+Original author:
+ Mohammad Akhlaghi <address@hidden>
+Contributing author(s):
+Copyright (C) 2015, Free Software Foundation, Inc.
+
+Gnuastro is free software: you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation, either version 3 of the License, or (at your
+option) any later version.
+
+Gnuastro is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
+**********************************************************************/
+#include <config.h>
+
+#include <stdio.h>
+#include <errno.h>
+#include <error.h>
+#include <string.h>
+#include <stdlib.h>
+
+#include <gnuastro/fits.h>
+#include <gnuastro/mesh.h>
+#include <gnuastro/qsort.h>
+#include <gnuastro/array.h>
+#include <gnuastro/blank.h>
+#include <gnuastro/convolve.h>
+#include <gnuastro/statistics.h>
+#include <gnuastro/interpolate.h>
+
+#include <timing.h>
+#include <checkset.h>
+
+#include "main.h"
+
+
+
+
+
+static void *
+sky_on_thread(void *in_prm)
+{
+ struct gal_threads_params *tprm=(struct gal_threads_params *)in_prm;
+ struct statisticsparams *p=(struct statisticsparams *)tprm->params;
+ int type=p->input->type, stype=GAL_DATA_TYPE_FLOAT32;
+
+ double *darr;
+ void *tblock=NULL, *tarray=NULL;
+ gal_data_t *tile, *mode, *sigmaclip;
+ size_t i, tind, twidth=gal_data_sizeof(stype);
+
+
+ /* Find the Sky and its standard deviation on the tiles given to this
+ thread. */
+ for(i=0; tprm->indexs[i] != GAL_THREADS_NON_THRD_INDEX; ++i)
+ {
+ /* Set the tile and copy its values into the array we'll be using. */
+ tind = tprm->indexs[i];
+ tile = &p->cp.tl.tiles[tind];
+
+
+ /* If we have a convolved image, temporarily (only for finding the
+ mode) change the tile's pointers so we can work on the convolved
+ image for the mode. */
+ if(p->kernel)
+ {
+ tarray=tile->array; tblock=tile->block;
+ tile->array=gal_tile_block_relative_to_other(tile, p->convolved);
+ tile->block=p->convolved;
+ }
+ mode=gal_statistics_mode(tile, p->mirrordist, 1);
+ if(p->kernel) { tile->array=tarray; tile->block=tblock; }
+
+
+ /* Check the mode value. Note that if the mode is in-accurate, then
+ the values will be NaN and all conditionals will fail. So, we'll
+ go onto finding values for this tile */
+ darr=mode->array;
+ if( fabs(darr[1]-0.5f) < p->modmedqdiff )
+ {
+ /* Get the sigma-clipped mean and standard deviation. `inplace'
+ is irrelevant here because this is a tile and it will be
+ copied anyway. */
+ sigmaclip=gal_statistics_sigma_clip(tile, p->sclipparams[0],
+ p->sclipparams[1], 1, 1);
+
+ /* Put the mean and its standard deviation into the respective
+ place for this tile. */
+ sigmaclip=gal_data_copy_to_new_type_free(sigmaclip, stype);
+ memcpy(gal_data_ptr_increment(p->sky_t->array, tind, stype),
+ gal_data_ptr_increment(sigmaclip->array, 2, stype), twidth);
+ memcpy(gal_data_ptr_increment(p->std_t->array, tind, stype),
+ gal_data_ptr_increment(sigmaclip->array, 3, stype), twidth);
+
+ /* Clean up. */
+ gal_data_free(sigmaclip);
+ }
+ else
+ {
+ gal_blank_write(gal_data_ptr_increment(p->sky_t->array, tind, type),
+ type);
+ gal_blank_write(gal_data_ptr_increment(p->std_t->array, tind, type),
+ type);
+ }
+
+ /* Clean up. */
+ gal_data_free(mode);
+ }
+
+
+ /* Wait for all threads to finish and return. */
+ if(tprm->b) pthread_barrier_wait(tprm->b);
+ return NULL;
+}
+
+
+
+
+
+void
+sky(struct statisticsparams *p)
+{
+ char *msg, *outname;
+ struct timeval t0, t1;
+ gal_data_t *num, *tmp;
+ struct gal_options_common_params *cp=&p->cp;
+ struct gal_tile_two_layer_params *tl=&cp->tl;
+
+
+ /* Print basic information */
+ if(!cp->quiet)
+ {
+ gettimeofday(&t0, NULL);
+ printf("%s\n", PROGRAM_STRING);
+ printf("Estimating Sky (reference value) and its STD.\n");
+ printf("-----------\n");
+ printf(" - Using %zu CPU thread%s.\n", cp->numthreads,
+ cp->numthreads==1 ? "" : "s");
+ printf(" - Input: %s (hdu: %s)\n", p->inputname, cp->hdu);
+ if(p->kernelname)
+ printf(" - Kernel: %s (hdu: %s)\n", p->kernelname, p->khdu);
+ }
+
+ /* When checking steps, the input image is the first extension. */
+ if(p->checksky)
+ gal_fits_img_write(p->input, p->checkskyname, NULL, PROGRAM_STRING);
+
+
+ /* Convolve the image (if desired). */
+ if(p->kernel)
+ {
+ if(!cp->quiet) gettimeofday(&t1, NULL);
+ p->convolved=gal_convolve_spatial(tl->tiles, p->kernel,
+ cp->numthreads, 1, tl->workoverch);
+ if(p->checksky)
+ gal_fits_img_write(p->convolved, p->checkskyname, NULL,
+ PROGRAM_STRING);
+ if(!cp->quiet)
+ gal_timing_report(&t1, "Input convolved with kernel.", 1);
+ }
+
+
+ /* Make the arrays keeping the Sky and Sky standard deviation values. */
+ p->sky_t=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT32, p->input->ndim,
+ tl->numtiles, NULL, 0, p->input->minmapsize, "SKY",
+ p->input->unit, NULL);
+ p->std_t=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT32, p->input->ndim,
+ tl->numtiles, NULL, 0, p->input->minmapsize,
+ "SKY STD", p->input->unit, NULL);
+
+
+
+ /* Find the Sky and Sky standard deviation on the tiles. */
+ if(!cp->quiet) gettimeofday(&t1, NULL);
+ gal_threads_spin_off(sky_on_thread, p, tl->tottiles, cp->numthreads);
+ if(!cp->quiet)
+ {
+ num=gal_statistics_number(p->sky_t);
+ asprintf(&msg, "Sky and its STD found on %lu/%zu tiles.",
+ *((uint64_t *)(num->array)), tl->tottiles );
+ gal_timing_report(&t1, msg, 1);
+ gal_data_free(num);
+ free(msg);
+ }
+ if(p->checksky)
+ {
+ gal_tile_full_values_write(p->sky_t, tl, p->checkskyname,
+ PROGRAM_STRING);
+ gal_tile_full_values_write(p->std_t, tl, p->checkskyname,
+ PROGRAM_STRING);
+ }
+
+
+ /* Interpolate the Sky and its standard deviation. */
+ if(!cp->quiet) gettimeofday(&t1, NULL);
+ p->sky_t->next=p->std_t;
+ tmp=gal_interpolate_close_neighbors(p->sky_t, tl, cp->interpnumngb,
+ cp->numthreads, cp->interponlyblank, 1);
+ gal_data_free(p->sky_t);
+ gal_data_free(p->std_t);
+ p->sky_t=tmp;
+ p->std_t=tmp->next;
+ p->sky_t->next=p->std_t->next=NULL;
+ if(!cp->quiet)
+ gal_timing_report(&t1, "All blank tiles filled (interplated).", 1);
+ if(p->checksky)
+ {
+ gal_tile_full_values_write(p->sky_t, tl, p->checkskyname,
+ PROGRAM_STRING);
+ gal_tile_full_values_write(p->std_t, tl, p->checkskyname,
+ PROGRAM_STRING);
+ }
+
+
+ /* Smooth the Sky and Sky STD arrays. */
+ if(p->smoothwidth>1)
+ {
+ if(!cp->quiet) gettimeofday(&t1, NULL);
+ tmp=gal_tile_full_values_smooth(p->sky_t, tl, p->smoothwidth,
+ p->cp.numthreads);
+ gal_data_free(p->sky_t);
+ p->sky_t=tmp;
+ tmp=gal_tile_full_values_smooth(p->std_t, tl, p->smoothwidth,
+ p->cp.numthreads);
+ gal_data_free(p->std_t);
+ p->std_t=tmp;
+ if(!cp->quiet)
+ gal_timing_report(&t1, "Smoothed Sky and Sky STD values on tiles.",
+ 1);
+ if(p->checksky)
+ {
+ gal_tile_full_values_write(p->sky_t, tl, p->checkskyname,
+ PROGRAM_STRING);
+ gal_tile_full_values_write(p->std_t, tl, p->checkskyname,
+ PROGRAM_STRING);
+ }
+ }
+
+
+ /* Save the Sky and its standard deviation */
+ outname=gal_checkset_automatic_output(&p->cp,
+ ( p->cp.output
+ ? p->cp.output
+ : p->inputname ), "_sky.fits");
+ gal_tile_full_values_write(p->sky_t, tl, outname, PROGRAM_STRING);
+ gal_tile_full_values_write(p->std_t, tl, outname, PROGRAM_STRING);
+ if(!cp->quiet)
+ printf(" - Written to `%s'.\n", outname);
+
+
+ /* Clean up and return. */
+ free(outname);
+ gal_data_free(p->sky_t);
+ gal_data_free(p->std_t);
+ gal_data_free(p->convolved);
+
+ if(!cp->quiet)
+ {
+ printf("-----------\n");
+ gal_timing_report(&t0, "Completed in:", 0);
+ printf("-----------\n");
+ }
+}
diff --git a/bin/subtractsky/subtractsky.h b/bin/statistics/sky.h
similarity index 80%
rename from bin/subtractsky/subtractsky.h
rename to bin/statistics/sky.h
index 8fe8c49..3cdb2e1 100644
--- a/bin/subtractsky/subtractsky.h
+++ b/bin/statistics/sky.h
@@ -1,6 +1,6 @@
/*********************************************************************
-SubtractSky - Find and subtract the sky value from an image.
-SubtractSky is part of GNU Astronomy Utilities (Gnuastro) package.
+Statistics - Statistical analysis on input dataset.
+Statistics is part of GNU Astronomy Utilities (Gnuastro) package.
Original author:
Mohammad Akhlaghi <address@hidden>
@@ -20,10 +20,10 @@ General Public License for more details.
You should have received a copy of the GNU General Public License
along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
**********************************************************************/
-#ifndef SUBTRACTSKY_H
-#define SUBTRACTSKY_H
+#ifndef SKY_H
+#define SKY_H
void
-subtractsky(struct subtractskyparams *p);
+sky(struct statisticsparams *p);
#endif
diff --git a/bin/statistics/statistics.c b/bin/statistics/statistics.c
index a31c6e7..d9d141d 100644
--- a/bin/statistics/statistics.c
+++ b/bin/statistics/statistics.c
@@ -45,6 +45,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include "main.h"
#include "ui.h"
+#include "sky.h"
#include "statistics.h"
@@ -235,29 +236,25 @@ statistics_print_one_row(struct statisticsparams *p)
/*******************************************************************/
static void
statistics_interpolate_and_write(struct statisticsparams *p,
- gal_data_t *values)
+ gal_data_t *values, char *output)
{
gal_data_t *interpd;
- char *output=p->cp.output;
- struct gal_tile_two_layer_params *tl=&p->cp.tl;
-
+ struct gal_options_common_params *cp=&p->cp;
/* Do the interpolation (if necessary). */
if( p->interpolate
- && !(p->interponlyblank && gal_blank_present(values)==0) )
+ && !(p->cp.interponlyblank && gal_blank_present(values)==0) )
{
- interpd=gal_interpolate_close_neighbors(values, tl, p->interpnumngb,
- p->cp.numthreads,
- p->interponlyblank);
+ interpd=gal_interpolate_close_neighbors(values, &cp->tl,
+ cp->interpnumngb,
+ cp->numthreads,
+ cp->interponlyblank, 0);
gal_data_free(values);
values=interpd;
}
/* Write the values. */
- gal_tile_full_values_write(values, tl, output, PROGRAM_STRING);
-
- /* Clean up. */
- gal_data_free(values);
+ gal_tile_full_values_write(values, &cp->tl, output, PROGRAM_STRING);
}
@@ -267,36 +264,19 @@ statistics_interpolate_and_write(struct statisticsparams
*p,
static void
statistics_on_tile(struct statisticsparams *p)
{
- double arg;
- uint8_t type;
- size_t tind, dsize=1, mind;
- gal_data_t *tmp, *tmpv, *ttmp;
- gal_data_t *check, *tile, *values;
+ double arg=0;
+ gal_data_t *tile, *values;
+ size_t tind, dsize=1, mind=-1;
+ uint8_t type=GAL_DATA_TYPE_INVALID;
struct gal_linkedlist_ill *operation;
+ gal_data_t *tmp=NULL, *tmpv=NULL, *ttmp;
struct gal_options_common_params *cp=&p->cp;
struct gal_tile_two_layer_params *tl=&p->cp.tl;
-
- /* Set the output name. */
- if(cp->output==NULL)
- cp->output=gal_checkset_automatic_output(cp, p->inputname,
+ char *output=gal_checkset_automatic_output(cp, cp->output
+ ? cp->output
+ : p->inputname,
"_ontile.fits");
- /* Make the tile structure. */
- gal_tile_full_two_layers(p->input, tl);
-
-
- /* Save the tile IDs if they are requested. */
- if(tl->tilecheckname)
- {
- check=gal_tile_block_check_tiles(tl->tiles);
- if(p->inputformat==INPUT_FORMAT_IMAGE)
- gal_fits_img_write(check, tl->tilecheckname, NULL, PROGRAM_NAME);
- else
- gal_table_write(check, NULL, cp->tableformat, tl->tilecheckname,
- cp->dontdelete);
- gal_data_free(check);
- }
-
/* Do the operation on each tile. */
for(operation=p->singlevalue; operation!=NULL; operation=operation->next)
{
@@ -406,22 +386,22 @@ statistics_on_tile(struct statisticsparams *p)
/* Put the output value into the `values' array and clean up. */
tmp=gal_data_copy_to_new_type_free(tmp, type);
- memcpy(gal_data_ptr_increment(values->array, tind++,
- values->type),
+ memcpy(gal_data_ptr_increment(values->array, tind++, values->type),
tmp->array, gal_data_sizeof(type));
gal_data_free(tmp);
}
/* Do the interpolation (if necessary) and write the array into the
output. */
- statistics_interpolate_and_write(p, values);
+ statistics_interpolate_and_write(p, values, output);
/* Clean up. */
+ gal_data_free(values);
if(operation->v==ARGS_OPTION_KEY_QUANTFUNC) gal_data_free(tmpv);
}
/* Clean up. */
- gal_tile_full_free_contents(tl);
+ free(output);
}
@@ -880,36 +860,36 @@ print_sigma_clip(struct statisticsparams *p)
gal_data_t *sigclip;
/* Set the mode for printing: */
- if( p->sigclipparam>=1.0f )
- asprintf(&mode, "for %g clips", p->sigclipparam);
+ if( p->sclipparams[1]>=1.0f )
+ asprintf(&mode, "for %g clips", p->sclipparams[1]);
else
asprintf(&mode, "until relative change in STD is less than %g",
- p->sigclipparam);
+ p->sclipparams[1]);
/* Report the status */
if(!p->cp.quiet)
{
print_input_info(p);
- printf("%g-sigma clipping steps %s:\n\n", p->sigclipmultip, mode);
+ printf("%g-sigma clipping steps %s:\n\n", p->sclipparams[0], mode);
}
/* Do the Sigma clipping: */
- sigclip=gal_statistics_sigma_clip(p->sorted, p->sigclipmultip,
- p->sigclipparam, p->cp.quiet);
+ sigclip=gal_statistics_sigma_clip(p->sorted, p->sclipparams[0],
+ p->sclipparams[1], 0, p->cp.quiet);
a=sigclip->array;
/* Finish the introduction. */
if(!p->cp.quiet)
printf("-------\nSummary:\n");
else
- printf("%g-sigma clipped %s:\n", p->sigclipmultip, mode);
+ printf("%g-sigma clipped %s:\n", p->sclipparams[0], mode);
/* Print the final results: */
printf(" %-*s %zu\n", namewidth, "Number of input elements:",
p->input->size);
- if( p->sigclipparam < 1.0f )
+ if( p->sclipparams[1] < 1.0f )
printf(" %-*s %d\n", namewidth, "Number of clips:", sigclip->status);
- printf(" %-*s %g\n", namewidth, "Final number of elements:", a[0]);
+ printf(" %-*s %.0f\n", namewidth, "Final number of elements:", a[0]);
printf(" %-*s %g\n", namewidth, "Median:", a[1]);
printf(" %-*s %g\n", namewidth, "Mean:", a[2]);
printf(" %-*s %g\n", namewidth, "Standard deviation:", a[3]);
@@ -951,6 +931,13 @@ statistics(struct statisticsparams *p)
else statistics_print_one_row(p);
}
+ /* Find the Sky value if called. */
+ if(p->sky)
+ {
+ sky(p);
+ print_basic_info=0;
+ }
+
/* Print the ASCII plots if requested. */
if(p->asciihist || p->asciicfp)
{
@@ -966,7 +953,7 @@ statistics(struct statisticsparams *p)
}
/* Print the sigma-clipped results. */
- if( !isnan(p->sigclipmultip) )
+ if( p->sigmaclip )
{
print_basic_info=0;
print_sigma_clip(p);
diff --git a/bin/statistics/ui.c b/bin/statistics/ui.c
index 282c56f..b06b026 100644
--- a/bin/statistics/ui.c
+++ b/bin/statistics/ui.c
@@ -90,6 +90,7 @@ enum program_args_groups
{
ARGS_GROUP_SINGLE_VALUE = GAL_OPTIONS_GROUP_AFTER_COMMON,
ARGS_GROUP_PARTICULAR_STAT,
+ ARGS_GROUP_SKY,
ARGS_GROUP_HIST_CFP,
};
@@ -128,6 +129,7 @@ ui_initialize_options(struct statisticsparams *p,
cp->program_authors = PROGRAM_AUTHORS;
cp->coptions = gal_commonopts_options;
cp->numthreads = gal_threads_number();
+ cp->tl.remainderfrac = NAN;
/* Program-specific initializers */
p->lessthan = NAN;
@@ -137,8 +139,9 @@ ui_initialize_options(struct statisticsparams *p,
p->quantmax = NAN;
p->mirror = NAN;
p->mirrordist = NAN;
- p->sigclipparam = NAN;
- p->sigclipmultip = NAN;
+ p->modmedqdiff = NAN;
+ p->sclipparams[0] = NAN;
+ p->sclipparams[1] = NAN;
/* Set the mandatory common options. */
for(i=0; !gal_options_is_last(&cp->coptions[i]); ++i)
@@ -290,8 +293,8 @@ ui_add_to_single_value(struct argp_option *option, char
*arg,
static void *
-ui_parse_numbers(struct argp_option *option, char *arg,
- char *filename, size_t lineno, void *params)
+ui_read_quantile_range(struct argp_option *option, char *arg,
+ char *filename, size_t lineno, void *params)
{
char *str;
gal_data_t *in;
@@ -300,117 +303,47 @@ ui_parse_numbers(struct argp_option *option, char *arg,
/* For the `--printparams' (`-P') option:*/
if(lineno==-1)
{
- switch(option->key)
- {
- case ARGS_OPTION_KEY_SIGMACLIP:
- asprintf(&str, "%g,%g", p->sigclipmultip, p->sigclipparam);
- break;
- case ARGS_OPTION_KEY_QRANGE:
- if( isnan(p->quantmax) ) asprintf(&str, "%g", p->quantmin);
- else asprintf(&str, "%g,%g", p->quantmin, p->quantmax);
- break;
- default:
- error(EXIT_FAILURE, 0, "a bug! option `%s' not recognized "
- "in `ui_parse_numbers' (when called for printing). Please "
- "contact us at %s to fix the problem", option->name,
- PACKAGE_BUGREPORT);
- }
+ if( isnan(p->quantmax) ) asprintf(&str, "%g", p->quantmin);
+ else asprintf(&str, "%g,%g", p->quantmin, p->quantmax);
return str;
}
/* Parse the inputs. */
in=gal_options_parse_list_of_numbers(arg, filename, lineno);
- /* Checks depending on the option. */
- switch(option->key)
- {
- case ARGS_OPTION_KEY_SIGMACLIP:
- /* Check if there was only two numbers. */
- if(in->size!=2)
- error_at_line(EXIT_FAILURE, 0, filename, lineno, "the `--%s' "
- "option takes two values (separated by a comma) for "
- "defining the sigma-clip. However, %zu numbers were "
- "read in the string `%s' (value to this option).\n\n"
- "The first number is the multiple of sigma, and the "
- "second is either the tolerance (if its is less than "
- "1.0), or a specific number of times to clip (if it "
- "is equal or larger than 1.0).", option->name, in->size,
- arg);
-
- /* Read the values in. */
- p->sigclipparam = ((double *)(in->array))[1];
- p->sigclipmultip = ((double *)(in->array))[0];
-
- /* Some sanity checks. */
- if( p->sigclipmultip<=0 )
- error_at_line(EXIT_FAILURE, 0, filename, lineno, "the first value to "
- "the `--%s' option (multiple of sigma), must be "
- "greater than zero. From the string `%s' (value to "
- "this option), you have given a value of %g for the "
- "first value", option->name, arg, p->sigclipmultip);
-
- /* If the second value must also be positive. */
- if( p->sigclipparam<=0 )
- error_at_line(EXIT_FAILURE, 0, filename, lineno, "the second value "
- "to the `--%s' option (tolerance to stop clipping or "
- "number of clips), must be greater than zero. From "
- "the string `%s' (value to this option), you have "
- "given a value of %g for the second value",
- option->name, arg, p->sigclipparam);
-
- /* if the second value is larger or equal to 1.0, it must be an
- integer. */
- if( p->sigclipparam >= 1.0f && ceil(p->sigclipparam) != p->sigclipparam)
- error_at_line(EXIT_FAILURE, 0, filename, lineno, "when the second "
- "value to the `--%s' option is >=1, it is interpretted "
- "as an absolute number of clips. So it must be an "
- "integer. However, your second value is a floating "
- "point number: %g (parsed from `%s')", option->name,
- p->sigclipparam, arg);
- break;
-
- case ARGS_OPTION_KEY_QRANGE:
- /* Check if there was only two numbers. */
- if(in->size!=1 && in->size!=2)
- error_at_line(EXIT_FAILURE, 0, filename, lineno, "the `--%s' "
- "option takes one or two values values (separated by "
- "a comma) to define the range of used values with "
- "quantiles. However, %zu numbers were read in the "
- "string `%s' (value to this option).\n\n"
- "If there is only one number as input, it will be "
- "interpretted as the lower quantile (Q) range. The "
- "higher range will be set to the quantile (1-Q). "
- "When two numbers are given, they will be used as the "
- "lower and higher quantile range respectively",
- option->name, in->size, arg);
-
- /* Read the values in. */
- p->quantmin = ((double *)(in->array))[0];
- if(in->size==2) p->quantmax = ((double *)(in->array))[1];
-
- /* Make sure the values are between 0 and 1. */
- if( (p->quantmin<0 || p->quantmin>1)
- || ( !isnan(p->quantmax) && (p->quantmax<0 || p->quantmax>1) ) )
- error_at_line(EXIT_FAILURE, 0, filename, lineno, "values to the "
- "`--quantrange' option must be between 0 and 1 "
- "(inclusive). Your input was: `%s'", arg);
- break;
-
- /* When only one value is given, make sure it is less than 0.5. */
- if( !isnan(p->quantmax) && p->quantmin>0.5 )
- error(EXIT_FAILURE, 0, "%g>=0.5! When only one value is given to the "
- "`--%s' option, the range is defined as Q and 1-Q. Thus, the "
- "value must be less than 0.5", p->quantmin, option->name);
-
- default:
- error(EXIT_FAILURE, 0, "a bug! option `%s' not recognized "
- "in `ui_parse_numbers' (when called for printing). Please "
- "contact us at %s to fix the problem", option->name,
- PACKAGE_BUGREPORT);
- }
-
-
- /* No return value necessary for this function. */
+ /* Check if there was only two numbers. */
+ if(in->size!=1 && in->size!=2)
+ error_at_line(EXIT_FAILURE, 0, filename, lineno, "the `--%s' "
+ "option takes one or two values values (separated by "
+ "a comma) to define the range of used values with "
+ "quantiles. However, %zu numbers were read in the "
+ "string `%s' (value to this option).\n\n"
+ "If there is only one number as input, it will be "
+ "interpretted as the lower quantile (Q) range. The "
+ "higher range will be set to the quantile (1-Q). "
+ "When two numbers are given, they will be used as the "
+ "lower and higher quantile range respectively",
+ option->name, in->size, arg);
+
+ /* Read the values in. */
+ p->quantmin = ((double *)(in->array))[0];
+ if(in->size==2) p->quantmax = ((double *)(in->array))[1];
+
+ /* Make sure the values are between 0 and 1. */
+ if( (p->quantmin<0 || p->quantmin>1)
+ || ( !isnan(p->quantmax) && (p->quantmax<0 || p->quantmax>1) ) )
+ error_at_line(EXIT_FAILURE, 0, filename, lineno, "values to the "
+ "`--quantrange' option must be between 0 and 1 "
+ "(inclusive). Your input was: `%s'", arg);
+
+ /* When only one value is given, make sure it is less than 0.5. */
+ if( !isnan(p->quantmax) && p->quantmin>0.5 )
+ error(EXIT_FAILURE, 0, "%g>=0.5! When only one value is given to the "
+ "`--%s' option, the range is defined as Q and 1-Q. Thus, the "
+ "value must be less than 0.5", p->quantmin, option->name);
+
+ /* Clean up and return. */
+ gal_data_free(in);
return NULL;
}
@@ -442,6 +375,7 @@ static void
ui_read_check_only_options(struct statisticsparams *p)
{
struct gal_linkedlist_ill *tmp;
+ struct gal_tile_two_layer_params *tl=&p->cp.tl;
/* Check if the format of the output table is valid, given the type of
the output. */
@@ -454,15 +388,61 @@ ui_read_check_only_options(struct statisticsparams *p)
"(for example `--median') must be requested with the `--ontile' "
"option: there is no value to put in each tile");
- /* The tile mode cannot be called with any other modes. */
- if(p->ontile && (p->asciihist || p->asciicfp || p->histogram
- || p->cumulative || !isnan(p->sigclipmultip)
- || !isnan(p->mirror) ) )
- error(EXIT_FAILURE, 0, "`--ontile' cannot be called with any of the "
- "`particular' calculation options, for example `--histogram'. "
- "This is because these options deal with the whole dataset, "
- "but in tessellation mode, the input dataset is broken up into "
- "individual tiles");
+ /* Tessellation related options. */
+ if( p->ontile || p->sky )
+ {
+ /* The tile or sky mode cannot be called with any other modes. */
+ if(p->asciihist || p->asciicfp || p->histogram || p->cumulative
+ || p->sigmaclip || !isnan(p->mirror) )
+ error(EXIT_FAILURE, 0, "`--ontile' or `--sky' cannot be called with "
+ "any of the `particular' calculation options, for example "
+ "`--histogram'. This is because the latter work over the whole "
+ "dataset and element positions are changed, but in the former "
+ "positions are significant");
+
+ /* Make sure the tessellation defining options are given. */
+ if( tl->tilesize==NULL || tl->numchannels==NULL
+ || isnan(tl->remainderfrac) )
+ error(EXIT_FAILURE, 0, "`--tilesize', `--numchannels', and "
+ "`--remainderfrac' are mandatory options when dealing with "
+ "a tessellation (in `--ontile' or `--sky' mode). Atleast "
+ "one of these options wasn't given a value.");
+ }
+
+
+ /* In Sky mode, several options are mandatory. */
+ if( p->sky )
+ {
+ /* Mandatory options. */
+ if ( isnan(p->modmedqdiff) || isnan(p->sclipparams[0])
+ || p->cp.interpnumngb==0 || isnan(p->mirrordist) )
+ error(EXIT_FAILURE, 0, "`--modmedqdiff', `--sclipparams', "
+ "`--mirrordist', and `--interpnumngb' are mandatory with "
+ "`--sky'");
+
+ /* If mode and median distance is a reasonable value. */
+ if(p->modmedqdiff>0.5)
+ error(EXIT_FAILURE, 0, "%f not acceptable for `--modmedqdiff'. It "
+ "cannot take values larger than 0.5 (quantile of median)",
+ p->modmedqdiff);
+
+ /* If a kernel name has been given, we need the HDU. */
+ if(p->kernelname && gal_fits_name_is_fits(p->kernelname)
+ && p->khdu==NULL )
+ error(EXIT_FAILURE, 0, "no HDU specified for the kernel image. When "
+ "A HDU is necessary for FITS files. You can use the `--khdu' "
+ "(`-u') option and give it the HDU number (starting from "
+ "zero), extension name, or anything acceptable by CFITSIO");
+ }
+
+
+ /* Sigma-clipping needs `sclipparams'. */
+ if(p->sigmaclip && isnan(p->sclipparams[0]))
+ error(EXIT_FAILURE, 0, "`--sclipparams' is necessary with `--sigmaclip'. "
+ "`--sclipparams' takes two values (separated by a comma) for "
+ "defining the sigma-clip: the multiple of sigma, and tolerance "
+ "(<1) or number of clips (>1).");
+
/* If any of the mode measurements are requested, then `mirrordist' is
mandatory. */
@@ -707,8 +687,7 @@ ui_make_sorted_if_necessary(struct statisticsparams *p)
}
/* Check in the rest of the outputs. */
- if( is_necessary==0
- && ( !isnan(p->sigclipmultip) || !isnan(p->mirror) ) )
+ if( is_necessary==0 && ( p->sigmaclip || !isnan(p->mirror) ) )
is_necessary=1;
/* Do the sorting, we will keep the sorted array in a separate space,
@@ -802,7 +781,9 @@ ui_read_columns(struct statisticsparams *p)
void
ui_preparations(struct statisticsparams *p)
{
+ gal_data_t *check;
struct gal_options_common_params *cp=&p->cp;
+ struct gal_tile_two_layer_params *tl=&cp->tl;
/* Read the input. */
if(p->isfits && p->hdu_type==IMAGE_HDU)
@@ -819,33 +800,49 @@ ui_preparations(struct statisticsparams *p)
p->inputformat=INPUT_FORMAT_TABLE;
}
- /* Tile and channel sanity checks. */
- if(p->ontile)
+ /* Read the convolution kernel if necessary. */
+ if(p->sky && p->kernelname)
+ p->kernel=gal_fits_img_read_kernel(p->kernelname, p->khdu,
+ cp->minmapsize);
+
+ /* Tile and channel sanity checks and preparations. */
+ if(p->ontile || p->sky)
{
- /* Check the tiles. */
- gal_tile_full_sanity_check(p->inputname, p->cp.hdu, p->input, &cp->tl);
+ /* Check the tiles and make the tile structure. */
+ gal_tile_full_sanity_check(p->inputname, p->cp.hdu, p->input, tl);
+ gal_tile_full_two_layers(p->input, tl);
+ gal_tile_full_permutation(tl);
- /* Set the tile file name if the user wants to check the tiles. */
- if(cp->tl.checktiles)
+ /* Make the tile check image if requested. */
+ if(tl->checktiles)
{
- cp->tl.tilecheckname=gal_checkset_automatic_output(&p->cp,
- p->inputname,
- "_tiled.fits");
- gal_checkset_check_remove_file(cp->tl.tilecheckname, 0,
- cp->dontdelete);
+ tl->tilecheckname=gal_checkset_automatic_output(cp, p->inputname,
+ "_tiled.fits");
+ check=gal_tile_block_check_tiles(tl->tiles);
+ if(p->inputformat==INPUT_FORMAT_IMAGE)
+ gal_fits_img_write(check, tl->tilecheckname, NULL, PROGRAM_NAME);
+ else
+ gal_table_write(check, NULL, cp->tableformat, tl->tilecheckname,
+ cp->dontdelete);
+ gal_data_free(check);
}
+
+ /* Set the steps image name. */
+ if(p->sky && p->checksky)
+ p->checkskyname=gal_checkset_automatic_output(cp, p->inputname,
+ "_sky_steps.fits");
}
/* Set the out-of-range values in the input to blank. */
ui_out_of_range_to_blank(p);
- /* If we are not to work on tiles: */
- if(!p->ontile)
+ /* If we are not to work on tiles, then re-order and change the input. */
+ if(p->ontile==0 && p->sky==0)
{
- /* Only keep the numbers we want. */
+ /* Only keep the elements we want. */
gal_blank_remove(p->input);
- /* Make sure there is any data remaining: */
+ /* Make sure there is data remaining: */
if(p->input->size==0)
error(EXIT_FAILURE, 0, "%s: no data, maybe the `--greaterequal' or "
"`--lessthan' options need to be adjusted",
@@ -967,5 +964,6 @@ ui_free_report(struct statisticsparams *p)
gal_data_free(p->input);
gal_data_free(p->reference);
gal_linkedlist_free_dll(p->tp_args);
+ gal_tile_full_free_contents(&p->cp.tl);
gal_linkedlist_free_ill(p->singlevalue);
}
diff --git a/bin/statistics/ui.h b/bin/statistics/ui.h
index 2ab02b8..6d50231 100644
--- a/bin/statistics/ui.h
+++ b/bin/statistics/ui.h
@@ -29,7 +29,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
/* Available letters for short options:
- a b e f j k p v w x y z
+ a b e f j p v w x z
B G J L R W X Y
*/
enum option_keys_enum
@@ -52,6 +52,8 @@ enum option_keys_enum
ARGS_OPTION_KEY_NORMALIZE = 'n',
ARGS_OPTION_KEY_ONTILE = 't',
ARGS_OPTION_KEY_INTERPOLATE = 'i',
+ ARGS_OPTION_KEY_SKY = 'y',
+ ARGS_OPTION_KEY_KERNEL = 'k',
/* Only with long version (start with a value 1000, the rest will be set
automatically). */
@@ -71,9 +73,12 @@ enum option_keys_enum
ARGS_OPTION_KEY_LOWERBIN,
ARGS_OPTION_KEY_ONEBINSTART,
ARGS_OPTION_KEY_MAXBINONE,
+ ARGS_OPTION_KEY_KHDU,
ARGS_OPTION_KEY_MIRRORDIST,
- ARGS_OPTION_KEY_INTERPONLYBLANK,
- ARGS_OPTION_KEY_INTERPNUMNGB,
+ ARGS_OPTION_KEY_MODMEDQDIFF,
+ ARGS_OPTION_KEY_SMOOTHWIDTH,
+ ARGS_OPTION_KEY_CHECKSKY,
+ ARGS_OPTION_KEY_SCLIPPARAMS,
};
diff --git a/bin/subtractsky/Makefile.am b/bin/subtractsky/Makefile.am
deleted file mode 100644
index f56a27d..0000000
--- a/bin/subtractsky/Makefile.am
+++ /dev/null
@@ -1,40 +0,0 @@
-## Process this file with automake to produce Makefile.inx
-##
-## Original author:
-## Mohammad Akhlaghi <address@hidden>
-## Contributing author(s):
-## Copyright (C) 2015, Free Software Foundation, Inc.
-##
-## Gnuastro is free software: you can redistribute it and/or modify it
-## under the terms of the GNU General Public License as published by
-## the Free Software Foundation, either version 3 of the License, or
-## (at your option) any later version.
-##
-## Gnuastro is distributed in the hope that it will be useful, but
-## WITHOUT ANY WARRANTY; without even the implied warranty of
-## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-## General Public License for more details.
-##
-## You should have received a copy of the GNU General Public License
-## along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
-
-
-## Necessary flags. NOTE: $(top_srcdir)/bootstrapped/lib is only necessary
-## for internally compiled utilities and libraries. It must not be included
-## during the tests since the bootstrapped libraries are not installed.
-AM_CPPFLAGS = -I\$(top_srcdir)/bootstrapped/lib
-
-
-## Program definition (name, linking, sources and headers)
-bin_PROGRAMS = astsubtractsky
-
-astsubtractsky_LDADD = -lgnuastro
-
-astsubtractsky_SOURCES = main.c ui.c subtractsky.c
-
-EXTRA_DIST = main.h authors-cite.h args.h ui.h subtractsky.h
-
-
-## The configuration file (distribute and install).
-## NOTE: the man page is created in doc/Makefile.am
-dist_sysconf_DATA = astsubtractsky.conf
diff --git a/bin/subtractsky/args.h b/bin/subtractsky/args.h
deleted file mode 100644
index 411d647..0000000
--- a/bin/subtractsky/args.h
+++ /dev/null
@@ -1,505 +0,0 @@
-/*********************************************************************
-SubtractSky - Find and subtract the sky value from an image.
-SubtractSky is part of GNU Astronomy Utilities (Gnuastro) package.
-
-Original author:
- Mohammad Akhlaghi <address@hidden>
-Contributing author(s):
-Copyright (C) 2015, Free Software Foundation, Inc.
-
-Gnuastro is free software: you can redistribute it and/or modify it
-under the terms of the GNU General Public License as published by the
-Free Software Foundation, either version 3 of the License, or (at your
-option) any later version.
-
-Gnuastro is distributed in the hope that it will be useful, but
-WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-General Public License for more details.
-
-You should have received a copy of the GNU General Public License
-along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
-**********************************************************************/
-#ifndef ARGS_H
-#define ARGS_H
-
-#include <argp.h>
-
-#include <commonargs.h>
-#include <fixedstringmacros.h>
-
-
-
-
-
-
-
-
-
-
-/**************************************************************/
-/************** argp.h definitions ***************/
-/**************************************************************/
-
-
-
-
-/* Definition parameters for the argp: */
-const char *argp_program_version=SPACK_STRING"\n"GAL_STRINGS_COPYRIGHT
- "\n\nWritten by Mohammad Akhlaghi";
-const char *argp_program_bug_address=PACKAGE_BUGREPORT;
-static char args_doc[] = "ASTRdata";
-
-
-
-
-
-const char doc[] =
- /* Before the list of options: */
- GAL_STRINGS_TOP_HELP_INFO
- SPACK_NAME" Finds the sky value over a grid on the input and subtracts "
- "it from the image to give a clear and uniform output. \n"
- GAL_STRINGS_MORE_HELP_INFO
- /* After the list of options: */
- "\v"
- PACKAGE_NAME" home page: "PACKAGE_URL;
-
-
-
-
-
-/* Available letters for short options:
-
- c e f g i j l m p r v w x y z
- A B C E F G I J O R W X Y Z
-
- Number keys free: >=510
-
- Options with keys (second structure element) larger than 500 do not
- have a short version.
- */
-static struct argp_option options[] =
- {
- {
- 0, 0, 0, 0,
- "Input:",
- 1
- },
- {
- "mask",
- 'M',
- "STR",
- 0,
- "Mask image file name.",
- 1
- },
- {
- "mhdu",
- 'H',
- "STR",
- 0,
- "Mask image header name.",
- 1
- },
- {
- "kernel",
- 'k',
- "STR",
- 0,
- "Kernel image file name for convolution.",
- 1
- },
- {
- "khdu",
- 'U',
- "STR",
- 0,
- "Kernel image header name for convolution.",
- 1
- },
-
-
- {
- 0, 0, 0, 0,
- "Output:",
- 2
- },
- {
- "checksky",
- 502,
- 0,
- 0,
- "Store final sky and its STD in `_sky.fits' file.",
- 2
- },
- {
- "checkskystd",
- 505,
- 0,
- 0,
- "Include sky standard deviation in all checks too.",
- 2
- },
- {
- "checkconvolution",
- 507,
- 0,
- 0,
- "Store convolved image in `_conv.fits' file.",
- 2
- },
-
-
-
-
- {
- 0, 0, 0, 0,
- "Mesh grid:",
- 3
- },
- {
- "meshsize",
- 's',
- "INT",
- 0,
- "Size of each mesh (tile) in the grid.",
- 3
- },
- {
- "nch1",
- 'a',
- "INT",
- 0,
- "Number of channels along first FITS axis.",
- 3
- },
- {
- "nch2",
- 'b',
- "INT",
- 0,
- "Number of channels along second FITS axis.",
- 3
- },
- {
- "lastmeshfrac",
- 'L',
- "INT",
- 0,
- "Fraction of last mesh area to add new.",
- 3
- },
- {
- "mirrordist",
- 'd',
- "FLT",
- 0,
- "Distance beyond mirror point. Multiple of std.",
- 3
- },
- {
- "minmodeq",
- 'Q',
- "FLT",
- 0,
- "Minimum acceptable quantile for the mode.",
- 3
- },
- {
- "interponlyblank",
- 508,
- 0,
- 0,
- "Only interpolate over the blank pixels.",
- 3
- },
- {
- "numnearest",
- 'n',
- "INT",
- 0,
- "Number of nearest neighbors to interpolate.",
- 3
- },
- {
- "smoothwidth",
- 'T',
- "INT",
- 0,
- "Width of smoothing kernel (odd number).",
- 3
- },
- {
- "fullconvolution",
- 506,
- 0,
- 0,
- "Ignore channels in imageconvolution.",
- 3
- },
- {
- "fullinterpolation",
- 503,
- 0,
- 0,
- "Ignore channels in interpolation.",
- 3
- },
- {
- "fullsmooth",
- 504,
- 0,
- 0,
- "Ignore channels in smoothing.",
- 3
- },
- {
- "checkmesh",
- 500,
- 0,
- 0,
- "Store mesh IDs in `_mesh.fits' file.",
- 3
- },
- {
- "meshbasedcheck",
- 501,
- 0,
- 0,
- "Each mesh in one pixel in mesh check images.",
- 3
- },
-
-
-
- {
- 0, 0, 0, 0,
- "Statistics:",
- 4
- },
- {
- "sigclipmultip",
- 'u',
- "FLT",
- 0,
- "Multiple of standard deviation in sigma-clipping.",
- 4
- },
- {
- "sigcliptolerance",
- 't',
- "FLT",
- 0,
- "Difference in STD tolerance to halt iteration.",
- 4
- },
-
-
-
- {
- 0, 0, 0, 0,
- "Operating modes:",
- -1
- },
-
-
-
- {0}
- };
-
-
-
-
-
-/* Parse a single option: */
-static error_t
-parse_opt(int key, char *arg, struct argp_state *state)
-{
- /* Save the arguments structure: */
- struct subtractskyparams *p = state->input;
-
- /* Set the pointer to the common parameters for all programs
- here: */
- state->child_inputs[0]=&p->cp;
-
- /* In case the user incorrectly uses the equal sign (for example
- with a short format or with space in the long format, then `arg`
- start with (if the short version was called) or be (if the long
- version was called with a space) the equal sign. So, here we
- check if the first character of arg is the equal sign, then the
- user is warned and the program is stopped: */
- if(arg && arg[0]=='=')
- argp_error(state, "incorrect use of the equal sign (`=`). For short "
- "options, `=` should not be used and for long options, "
- "there should be no space between the option, equal sign "
- "and value");
-
- switch(key)
- {
-
-
- /* Input: */
- case 'M':
- gal_checkset_allocate_copy_set(arg, &p->up.maskname, &p->up.masknameset);
- break;
- case 'H':
- gal_checkset_allocate_copy_set(arg, &p->up.mhdu, &p->up.mhduset);
- break;
- case 'K':
- gal_checkset_allocate_copy_set(arg, &p->up.kernelname,
- &p->up.kernelnameset);
- break;
- case 'U':
- gal_checkset_allocate_copy_set(arg, &p->up.khdu, &p->up.khduset);
- break;
-
- /* Output: */
- case 502:
- p->skyname="a";
- break;
- case 505:
- p->checkstd=1;
- break;
-
- /* Mesh grid: */
- case 's':
- gal_checkset_sizet_l_zero(arg, &p->mp.meshsize, "meshsize", key, SPACK,
- NULL, 0);
- p->up.meshsizeset=1;
- break;
- case 'a':
- gal_checkset_sizet_l_zero(arg, &p->mp.nch1, "nch1", key, SPACK, NULL, 0);
- p->up.nch1set=1;
- break;
- case 'b':
- gal_checkset_sizet_l_zero(arg, &p->mp.nch2, "nch2", key, SPACK, NULL, 0);
- p->up.nch2set=1;
- break;
- case 'L':
- gal_checkset_float_l_0_s_1(arg, &p->mp.lastmeshfrac, "lastmeshfrac", key,
- SPACK, NULL, 0);
- p->up.lastmeshfracset=1;
- break;
- case 'd':
- gal_checkset_float_l_0(arg, &p->mp.mirrordist, "mirrordist", key, SPACK,
- NULL, 0);
- p->up.mirrordistset=1;
- break;
- case 'Q':
- gal_checkset_float_l_0_s_1(arg, &p->mp.minmodeq, "minmodeq", key, SPACK,
- NULL, 0);
- p->up.minmodeqset=1;
- break;
- case 508:
- p->mp.interponlyblank=1;
- break;
- case 'n':
- gal_checkset_sizet_l_zero(arg, &p->mp.numnearest, "numnearest", key,
- SPACK, NULL, 0);
- p->up.numnearestset=1;
- break;
- case 'T':
- gal_checkset_sizet_p_odd(arg, &p->mp.smoothwidth, "smoothwidth", key,
- SPACK, NULL, 0);
- p->up.smoothwidthset=1;
- break;
- case 506:
- p->mp.fullconvolution=1;
- p->up.fullconvolutionset=1;
- break;
- case 503:
- p->mp.fullinterpolation=1;
- p->up.fullinterpolationset=1;
- break;
- case 504:
- p->mp.fullsmooth=1;
- p->up.fullsmoothset=1;
- break;
- case 500:
- p->meshname="a"; /* Just a placeholder! It will be corrected later */
- break;
- case 507:
- p->convname="a";
- break;
- case 501:
- p->mp.meshbasedcheck=1;
- break;
-
-
- /* Statistics */
- case 'u':
- gal_checkset_float_l_0(arg, &p->sigclipmultip, "sigclipmultip", key,
SPACK,
- NULL, 0);
- p->up.sigclipmultipset=1;
- break;
- case 't':
- gal_checkset_float_l_0_s_1(arg, &p->sigcliptolerance, "sigcliptolerance",
- key, SPACK, NULL, 0);
- p->up.sigcliptoleranceset=1;
- break;
-
-
- /* Operating modes: */
-
-
- /* Read the non-option arguments: */
- case ARGP_KEY_ARG:
-
- /* See what type of input value it is and put it in. */
- if( gal_fits_name_is_fits(arg) )
- {
- if(p->up.inputname)
- argp_error(state, "only one input image should be given");
- else
- p->up.inputname=arg;
- }
- else
- argp_error(state, "%s is not a valid file type", arg);
- break;
-
-
-
-
-
- /* The command line options and arguments are finished. */
- case ARGP_KEY_END:
- if(p->cp.setdirconf==0 && p->cp.setusrconf==0
- && p->cp.printparams==0)
- {
- if(state->arg_num==0)
- argp_error(state, "no argument given");
- if(p->up.inputname==NULL)
- argp_error(state, "no input FITS image(s) provided");
- }
- break;
-
-
-
-
-
- default:
- return ARGP_ERR_UNKNOWN;
- }
- return 0;
-}
-
-
-
-
-
-/* Specify the children parsers: */
-struct argp_child children[]=
- {
- {&commonargp, 0, NULL, 0},
- {0, 0, 0, 0}
- };
-
-
-
-
-
-/* Basic structure defining the whole argument reading process. */
-static struct argp thisargp = {options, parse_opt, args_doc,
- doc, children, NULL, NULL};
-
-#endif
diff --git a/bin/subtractsky/astsubtractsky.conf
b/bin/subtractsky/astsubtractsky.conf
deleted file mode 100644
index 14f8aa7..0000000
--- a/bin/subtractsky/astsubtractsky.conf
+++ /dev/null
@@ -1,40 +0,0 @@
-# Default parameters (System) for SubtractSky.
-# SubtractSky is part of GNU Astronomy Utitlies.
-#
-# Use the long option name of each paramter followed by
-# a value. The name and value should be separated by
-# atleast one of the following charaacters:
-# space, `,`, `=` or `:`
-#
-# Run with `--help` option or read the manual for a full
-# explanation of what each option means.
-#
-# NOTE I: All counting is from zero, not one.
-# NOTE II: Lines starting with `#` are ignored.
-#
-# Copying and distribution of this file, with or without modification,
-# are permitted in any medium without royalty provided the copyright
-# notice and this notice are preserved. This file is offered as-is,
-# without any warranty.
-
-# Input:
- khdu 1
-
-# Output:
-
-# Mesh grid:
- meshsize 32
- nch1 1
- nch2 1
- lastmeshfrac 0.6
- numnearest 9
- smoothwidth 5
- fullconvolution 0
- fullinterpolation 0
- fullsmooth 0
-
-# Statistics:
- mirrordist 1.5
- minmodeq 0.49
- sigclipmultip 3
- sigcliptolerance 0.2
\ No newline at end of file
diff --git a/bin/subtractsky/authors-cite.h b/bin/subtractsky/authors-cite.h
deleted file mode 100644
index ad08147..0000000
--- a/bin/subtractsky/authors-cite.h
+++ /dev/null
@@ -1,38 +0,0 @@
-/*********************************************************************
-SubtractSky - Find and subtract the sky value from an image.
-SubtractSky is part of GNU Astronomy Utilities (Gnuastro) package.
-
-Original author:
- Mohammad Akhlaghi <address@hidden>
-Contributing author(s):
-Copyright (C) 2017, Free Software Foundation, Inc.
-
-Gnuastro is free software: you can redistribute it and/or modify it
-under the terms of the GNU General Public License as published by the
-Free Software Foundation, either version 3 of the License, or (at your
-option) any later version.
-
-Gnuastro is distributed in the hope that it will be useful, but
-WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-General Public License for more details.
-
-You should have received a copy of the GNU General Public License
-along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
-**********************************************************************/
-#ifndef AUTHORS_CITE_H
-#define AUTHORS_CITE_H
-
-/* When any specific citation is necessary, please add its BibTeX (from ADS
- hopefully) to this variable along with a title decribing what this
- paper/book does for the progarm in a short line. In the following line
- put a row of `-' with the same length and then put the BibTeX.
-
- See the `gnuastro_bibtex' variable in `lib/options' (from the top
- Gnuastro source code directory) as an example.*/
-
-#define PROGRAM_BIBTEX "";
-
-#define PROGRAM_AUTHORS "Mohammad Akhlaghi";
-
-#endif
diff --git a/bin/subtractsky/main.c b/bin/subtractsky/main.c
deleted file mode 100644
index 7a9729f..0000000
--- a/bin/subtractsky/main.c
+++ /dev/null
@@ -1,56 +0,0 @@
-/*********************************************************************
-SubtractSky - Find and subtract the sky value from an image.
-SubtractSky is part of GNU Astronomy Utilities (Gnuastro) package.
-
-Original author:
- Mohammad Akhlaghi <address@hidden>
-Contributing author(s):
-Copyright (C) 2015, Free Software Foundation, Inc.
-
-Gnuastro is free software: you can redistribute it and/or modify it
-under the terms of the GNU General Public License as published by the
-Free Software Foundation, either version 3 of the License, or (at your
-option) any later version.
-
-Gnuastro is distributed in the hope that it will be useful, but
-WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-General Public License for more details.
-
-You should have received a copy of the GNU General Public License
-along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
-**********************************************************************/
-#include <config.h>
-
-#include <stdio.h>
-#include <stdlib.h>
-
-#include <timing.h> /* Includes time.h and sys/time.h */
-
-#include "main.h"
-
-#include "ui.h" /* needs main.h. */
-#include "subtractsky.h" /* needs main.h. */
-
-int
-main (int argc, char *argv[])
-{
- struct timeval t1;
- struct subtractskyparams p={{{0},0},0};
-
- /* Set the starting time. */
- time(&p.rawtime);
- gettimeofday(&t1, NULL);
-
- /* Read the input parameters. */
- setparams(argc, argv, &p);
-
- /* Run MakeProfiles */
- subtractsky(&p);
-
- /* Free all non-freed allocations. */
- freeandreport(&p, &t1);
-
- /* Return successfully.*/
- return EXIT_SUCCESS;
-}
diff --git a/bin/subtractsky/main.h b/bin/subtractsky/main.h
deleted file mode 100644
index 8205be7..0000000
--- a/bin/subtractsky/main.h
+++ /dev/null
@@ -1,100 +0,0 @@
-/*********************************************************************
-SubtractSky - Find and subtract the sky value from an image.
-SubtractSky is part of GNU Astronomy Utilities (Gnuastro) package.
-
-Original author:
- Mohammad Akhlaghi <address@hidden>
-Contributing author(s):
-Copyright (C) 2015, Free Software Foundation, Inc.
-
-Gnuastro is free software: you can redistribute it and/or modify it
-under the terms of the GNU General Public License as published by the
-Free Software Foundation, either version 3 of the License, or (at your
-option) any later version.
-
-Gnuastro is distributed in the hope that it will be useful, but
-WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-General Public License for more details.
-
-You should have received a copy of the GNU General Public License
-along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
-**********************************************************************/
-#ifndef MAIN_H
-#define MAIN_H
-
-#include <gnuastro/mesh.h>
-#include <gnuastro/fits.h>
-
-#include <commonparams.h>
-
-/* Progarm name macros: */
-#define SPACK "astsubtractsky" /* Subpackage executable name. */
-#define SPACK_NAME "SubtractSky" /* Subpackage full name. */
-#define SPACK_STRING SPACK_NAME" ("PACKAGE_NAME") "PACKAGE_VERSION
-
-
-
-
-
-struct uiparams
-{
- char *inputname; /* Name of input file. */
- char *maskname; /* Name of mask image file. */
- char *mhdu; /* Name of mask image header name. */
- char *kernelname;
- char *khdu;
- int masknameset;
- int mhduset;
- int kernelnameset;
- int khduset;
- int numnearestset;
- int smoothwidthset;
- int mirrordistset;
- int minmodeqset;
- int fullconvolutionset;
- int fullinterpolationset;
- int fullsmoothset;
- int sigclipmultipset;
- int sigcliptoleranceset;
- int meshsizeset;
- int nch1set;
- int nch2set;
- int lastmeshfracset;
-};
-
-
-
-
-
-struct subtractskyparams
-{
- /* Other structures: */
- struct uiparams up; /* User interface parameters. */
- struct gal_commonparams cp; /* Common parameters. */
- struct gal_mesh_params mp; /* Mesh grid of input image. */
-
- /* Input: */
- int nwcs; /* Number of WCS structures. */
- struct wcsprm *wcs; /* Pointer to WCS structures. */
- int bitpix; /* Input image bitpix value. */
- int anyblank; /* ==1: thereare blank pixels in input image. */
-
- /* output: */
- int checkstd; /* ==1: include the sky STD in checks. */
- char *meshname; /* Name of --checkmesh output. */
- char *convname; /* Name of --checkconvolution output. */
- char *skyname; /* Name of sky and its STD image. */
-
- /* Statistics: */
- float sigclipmultip; /* Multiple of standard deviation, sigma clip. */
- float sigcliptolerance; /* Tolerance in sigma clip. */
-
- /* Operating mode: */
-
- /* Internal: */
- float *conv; /* Convolved input image. */
- time_t rawtime; /* Starting time of the program. */
-};
-
-#endif
diff --git a/bin/subtractsky/subtractsky.c b/bin/subtractsky/subtractsky.c
deleted file mode 100644
index a6e91f2..0000000
--- a/bin/subtractsky/subtractsky.c
+++ /dev/null
@@ -1,247 +0,0 @@
-/*********************************************************************
-SubtractSky - Find and subtract the sky value from an image.
-SubtractSky is part of GNU Astronomy Utilities (Gnuastro) package.
-
-Original author:
- Mohammad Akhlaghi <address@hidden>
-Contributing author(s):
-Copyright (C) 2015, Free Software Foundation, Inc.
-
-Gnuastro is free software: you can redistribute it and/or modify it
-under the terms of the GNU General Public License as published by the
-Free Software Foundation, either version 3 of the License, or (at your
-option) any later version.
-
-Gnuastro is distributed in the hope that it will be useful, but
-WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-General Public License for more details.
-
-You should have received a copy of the GNU General Public License
-along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
-**********************************************************************/
-#include <config.h>
-
-#include <stdio.h>
-#include <errno.h>
-#include <error.h>
-#include <stdlib.h>
-
-#include <gnuastro/fits.h>
-#include <gnuastro/mesh.h>
-#include <gnuastro/qsort.h>
-#include <gnuastro/array.h>
-#include <gnuastro/statistics.h>
-
-#include <timing.h>
-
-#include "main.h"
-
-
-void *
-avestdonthread(void *inparam)
-{
- struct gal_mesh_thread_params *mtp=(struct gal_mesh_thread_params *)inparam;
- struct gal_mesh_params *mp=mtp->mp;
- struct subtractskyparams *p=(struct subtractskyparams *)mp->params;
-
- float *mponeforall=mp->oneforall;
- float *oneforall=&mponeforall[mtp->id*mp->maxs0*mp->maxs1];
-
- float *f, *cofa, *c=NULL; /* cofa: convolved-oneforall */
- size_t modeindex=(size_t)(-1);
- size_t s0, s1, ind, row, num, start, is1=mp->s1;
- size_t i, *indexs=&mp->indexs[mtp->id*mp->thrdcols];
- float mirrordist=mp->mirrordist, minmodeq=mp->minmodeq;
- float ave, med, std, *sorted, sigclipmultip=p->sigclipmultip;
- float *imgend, *inimg=mp->img, *inconv=p->conv, modesym=0.0f;
- float *conv=NULL, *img, sigcliptolerance=p->sigcliptolerance;
-
- /* Allocate the oneforall array for the convolved image, since we
- only want to use those meshs whose convolved image mode is larger
- than minmodeq. */
- if(p->conv!=mp->img)
- {
- errno=0; cofa=malloc(mp->maxs0*mp->maxs1*sizeof *oneforall);
- if(cofa==NULL)
- error(EXIT_FAILURE, errno, "unable to allocate %zu bytes for"
- "cofa in avestdonthread of subtractsky.c",
- mp->maxs0*mp->maxs1*sizeof *mp->img);
- }
- else cofa=NULL;
-
- /* Start this thread's work: */
- for(i=0;indexs[i]!=GAL_THREADS_NON_THRD_INDEX;++i)
- {
- /* Prepare the values: */
- num=row=0;
- f=oneforall;
- ind=indexs[i];
- start=mp->start[ind];
- s0=mp->ts0[mp->types[ind]];
- s1=mp->ts1[mp->types[ind]];
- if(cofa) c=sorted=cofa; else sorted=oneforall;
-
- /* Copy all the non-NaN pixels images pixels of this mesh into
- the mesh array. Note that currently, the spatial positioning
- of the pixels is irrelevant, so we only keep those that are
- non-NaN. Recall that both the convolved an unconvolved image
- have the same NaN pixels.*/
- do
- {
- if(cofa) conv = inconv + start + row * is1;
- imgend=(img = inimg + start + row++ * is1 ) + s1;
- do
- {
- if(!isnan(*img))
- {
- ++num;
- *f++ = *img;
- if(cofa) *c++=*conv;
- }
- if(cofa) ++conv;
- }
- while(++img<imgend);
- }
- while(row<s0);
-
- /* Do the desired operation on the mesh: */
- qsort(sorted, num, sizeof *oneforall, gal_qsort_float_increasing);
- gal_statistics_mode_index_in_sorted(sorted, num, mirrordist,
- &modeindex, &modesym);
- if( modesym>GAL_STATISTICS_MODE_SYM_GOOD
- && (float)modeindex/(float)num>minmodeq )
- {
- /* If cofa was defined, then oneforall was not sorted. */
- if(cofa)
- qsort(oneforall, num, sizeof *oneforall,
- gal_qsort_float_increasing);
-
- /* Do sigma-clipping and save the result if it is
- accurate. Note that all meshs were initialized to NaN, so
- if they don't fit the criteria, they can simply be
- ignored. */
- if(gal_statistics_sigma_clip_converge(oneforall, 1, num,
- sigclipmultip,
sigcliptolerance,
- &ave, &med, &std, 0))
- {
- mp->cgarray1[ind]=ave;
- if(mp->ngarrays==2) mp->cgarray2[ind]=std;
- }
- }
- }
-
- /* Free any allocated space and if multiple threads were used, wait
- until all other threads finish. */
- if(p->conv!=mp->img) free(cofa);
- if(mp->numthreads>1)
- pthread_barrier_wait(&mp->b);
- return NULL;
-}
-
-
-
-
-
-void
-subtractsky(struct subtractskyparams *p)
-{
- struct gal_mesh_params *mp=&p->mp;
-
- long *meshindexs;
- struct timeval t1;
- int checkstd=p->checkstd;
- size_t s0=mp->s0, s1=mp->s1;
- float *sky, *std, *skysubtracted;
-
-
- /* Prepare the mesh array. */
- gettimeofday(&t1, NULL);
- gal_mesh_make_mesh(mp);
- if(p->meshname)
- {
- gal_mesh_check_mesh_id(mp, &meshindexs);
- gal_fits_array_to_file(p->meshname, "Input", FLOAT_IMG,
- p->mp.img, s0, s1, p->anyblank, p->wcs,
- NULL, SPACK_STRING);
- gal_fits_array_to_file(p->meshname, "MeshIndexs", LONG_IMG,
- meshindexs, s0, s1, 0, p->wcs,
- NULL, SPACK_STRING);
- free(meshindexs);
- }
- if(p->cp.verb) gal_timing_report(&t1, "Mesh grid ready.", 1);
-
-
-
- /* Convolve the image if the user has asked for it: */
- if(p->up.kernelnameset)
- {
- gal_mesh_spatial_convolve_on_mesh(mp, &p->conv);
- if(p->convname)
- {
- gal_fits_array_to_file(p->convname, "Input", FLOAT_IMG,
- p->mp.img, s0, s1, p->anyblank,
- p->wcs, NULL, SPACK_STRING);
- gal_fits_array_to_file(p->convname, "Input", FLOAT_IMG,
- p->conv, s0, s1, p->anyblank, p->wcs,
- NULL, SPACK_STRING);
- }
- }
- else p->conv=p->mp.img;
- if(p->cp.verb)
- gal_timing_report(&t1, "Input image convolved with kernel.", 1);
-
-
-
- /* Find the sky value and its standard deviation on each mesh. */
- gal_mesh_operate_on_mesh(mp, avestdonthread, sizeof(float), checkstd, 1);
- if(p->skyname)
- gal_mesh_value_file(mp, p->skyname, "Sky value", "Sky STD", p->wcs,
- SPACK_STRING);
- if(p->cp.verb)
- gal_timing_report(&t1, "Sky and its STD found on some meshes.", 1);
-
-
-
- /* Interpolate over the meshs to fill all the blank ones in both the
- sky and the standard deviation arrays: */
- gal_mesh_interpolate(mp, "Interpolating the sky and its standard deviation");
- if(p->skyname)
- gal_mesh_value_file(mp, p->skyname, "Sky Interpolated",
- "Sky STD interpolated", p->wcs, SPACK_STRING);
- if(p->cp.verb)
- gal_timing_report(&t1, "All blank meshs filled (interplated).", 1);
-
-
-
- /* Smooth the interpolated array: */
- if(mp->smoothwidth>1)
- {
- gal_mesh_smooth(mp);
- if(p->cp.verb)
- gal_timing_report(&t1, "Mesh grid smoothed.", 1);
- }
-
-
- /* Make the sky array and save it if the user has asked for it: */
- gal_mesh_check_garray(mp, &sky, &std);
- if(p->skyname)
- gal_mesh_value_file(mp, p->skyname, "Sky Smoothed", "Sky STD smoothed",
- p->wcs, SPACK_STRING);
-
-
- /* Subtract the sky value */
- gal_array_fmultip_const(sky, s0*s1, -1.0f);
- skysubtracted=gal_array_fsum_arrays(mp->img, sky, s0*s1);
- gal_fits_array_to_file(p->cp.output ,"SkySubtracted", FLOAT_IMG,
- skysubtracted, s0, s1, p->anyblank, p->wcs,
- NULL, SPACK_STRING);
-
-
- /* Clean up: */
- free(sky);
- gal_mesh_free_mesh(mp);
- free(skysubtracted);
- if(checkstd) free(std);
- if(p->up.kernelnameset) free(p->conv);
-}
diff --git a/bin/subtractsky/ui.c b/bin/subtractsky/ui.c
deleted file mode 100644
index d170092..0000000
--- a/bin/subtractsky/ui.c
+++ /dev/null
@@ -1,627 +0,0 @@
-/*********************************************************************
-SubtractSky - Find and subtract the sky value from an image.
-SubtractSky is part of GNU Astronomy Utilities (Gnuastro) package.
-
-Original author:
- Mohammad Akhlaghi <address@hidden>
-Contributing author(s):
-Copyright (C) 2015, Free Software Foundation, Inc.
-
-Gnuastro is free software: you can redistribute it and/or modify it
-under the terms of the GNU General Public License as published by the
-Free Software Foundation, either version 3 of the License, or (at your
-option) any later version.
-
-Gnuastro is distributed in the hope that it will be useful, but
-WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-General Public License for more details.
-
-You should have received a copy of the GNU General Public License
-along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
-**********************************************************************/
-#include <config.h>
-
-#include <math.h>
-#include <stdio.h>
-#include <errno.h>
-#include <error.h>
-#include <stdlib.h>
-#include <string.h>
-#include <fitsio.h>
-
-#include <gnuastro/fits.h>
-#include <gnuastro/array.h>
-#include <gnuastro/txtarray.h>
-#include <gnuastro/statistics.h>
-
-#include <nproc.h> /* From Gnulib. */
-#include <timing.h> /* Includes time.h and sys/time.h */
-#include <checkset.h>
-#include <commonargs.h>
-#include <configfiles.h>
-
-#include "main.h"
-
-#include "ui.h" /* Needs main.h */
-#include "args.h" /* Needs main.h, includes argp.h. */
-
-
-/* Set the file names of the places where the default parameters are
- put. */
-#define CONFIG_FILE SPACK CONF_POSTFIX
-#define SYSCONFIG_FILE SYSCONFIG_DIR "/" CONFIG_FILE
-#define USERCONFIG_FILEEND USERCONFIG_DIR CONFIG_FILE
-#define CURDIRCONFIG_FILE CURDIRCONFIG_DIR CONFIG_FILE
-
-
-
-
-
-
-
-
-
-
-/**************************************************************/
-/************** Options and parameters ***************/
-/**************************************************************/
-void
-readconfig(char *filename, struct subtractskyparams *p)
-{
- FILE *fp;
- size_t lineno=0, len=200;
- char *line, *name, *value;
- struct uiparams *up=&p->up;
- struct gal_commonparams *cp=&p->cp;
- char key='a'; /* Not used, just a place holder. */
-
- /* When the file doesn't exist or can't be opened, it is ignored. It
- might be intentional, so there is no error. If a parameter is
- missing, it will be reported after all defaults are read. */
- fp=fopen(filename, "r");
- if (fp==NULL) return;
-
-
- /* Allocate some space for `line` with `len` elements so it can
- easily be freed later on. The value of `len` is arbitarary at
- this point, during the run, getline will change it along with the
- pointer to line. */
- errno=0;
- line=malloc(len*sizeof *line);
- if(line==NULL)
- error(EXIT_FAILURE, errno, "ui.c: %zu bytes in readdefaults",
- len * sizeof *line);
-
- /* Read the tokens in the file: */
- while(getline(&line, &len, fp) != -1)
- {
- /* Prepare the "name" and "value" strings, also set lineno. */
- GAL_CONFIGFILES_START_READING_LINE;
-
-
- /* Inputs: */
- if(strcmp(name, "hdu")==0)
- gal_checkset_allocate_copy_set(value, &cp->hdu, &cp->hduset);
-
- else if(strcmp(name, "mask")==0)
- gal_checkset_allocate_copy_set(value, &up->maskname,
- &up->masknameset);
-
- else if(strcmp(name, "mhdu")==0)
- gal_checkset_allocate_copy_set(value, &up->mhdu, &up->mhduset);
-
- else if(strcmp(name, "kernel")==0)
- gal_checkset_allocate_copy_set(value, &up->kernelname,
- &up->kernelnameset);
-
- else if(strcmp(name, "khdu")==0)
- gal_checkset_allocate_copy_set(value, &up->khdu, &up->khduset);
-
-
-
- /* Outputs */
- else if(strcmp(name, "output")==0)
- gal_checkset_allocate_copy_set(value, &cp->output, &cp->outputset);
-
-
- /* Mesh grid: */
- else if(strcmp(name, "meshsize")==0)
- {
- if(up->meshsizeset) continue;
- gal_checkset_sizet_l_zero(value, &p->mp.meshsize, name,
- key, SPACK, filename, lineno);
- up->meshsizeset=1;
- }
- else if(strcmp(name, "nch1")==0)
- {
- if(up->nch1set) continue;
- gal_checkset_sizet_l_zero(value, &p->mp.nch1, name, key, SPACK,
- filename, lineno);
- up->nch1set=1;
- }
- else if(strcmp(name, "nch2")==0)
- {
- if(up->nch2set) continue;
- gal_checkset_sizet_l_zero(value, &p->mp.nch2, name, key, SPACK,
- filename, lineno);
- up->nch2set=1;
- }
- else if(strcmp(name, "lastmeshfrac")==0)
- {
- if(up->lastmeshfracset) continue;
- gal_checkset_float_l_0_s_1(value, &p->mp.lastmeshfrac, name,
- key, SPACK, filename, lineno);
- up->lastmeshfracset=1;
- }
- else if(strcmp(name, "mirrordist")==0)
- {
- if(up->mirrordistset) continue;
- gal_checkset_float_l_0(value, &p->mp.mirrordist, name, key,
- SPACK, filename, lineno);
- up->mirrordistset=1;
- }
- else if(strcmp(name, "minmodeq")==0)
- {
- if(up->minmodeqset) continue;
- gal_checkset_float_l_0_s_1(value, &p->mp.minmodeq, name, key,
- SPACK, filename, lineno);
- up->minmodeqset=1;
- }
- else if(strcmp(name, "numnearest")==0)
- {
- if(up->numnearestset) continue;
- gal_checkset_sizet_l_zero(value, &p->mp.numnearest, name,
- key, SPACK, filename, lineno);
- up->numnearestset=1;
- }
- else if(strcmp(name, "smoothwidth")==0)
- {
- if(up->smoothwidthset) continue;
- gal_checkset_sizet_p_odd(value, &p->mp.smoothwidth, name,
- key, SPACK, filename, lineno);
- up->smoothwidthset=1;
- }
- else if(strcmp(name, "fullconvolution")==0)
- {
- if(up->fullconvolutionset) continue;
- gal_checkset_int_zero_or_one(value, &p->mp.fullconvolution,
- name, key, SPACK, filename, lineno);
- up->fullconvolutionset=1;
- }
- else if(strcmp(name, "fullinterpolation")==0)
- {
- if(up->fullinterpolationset) continue;
- gal_checkset_int_zero_or_one(value, &p->mp.fullinterpolation,
- name, key, SPACK, filename, lineno);
- up->fullinterpolationset=1;
- }
- else if(strcmp(name, "fullsmooth")==0)
- {
- if(up->fullsmoothset) continue;
- gal_checkset_int_zero_or_one(value, &p->mp.fullsmooth, name, key,
- SPACK, filename, lineno);
- up->fullsmoothset=1;
- }
-
-
- /* Statistics: */
- else if(strcmp(name, "sigclipmultip")==0)
- {
- if(up->sigclipmultipset) continue;
- gal_checkset_float_l_0(value, &p->sigclipmultip, name, key,
- SPACK, filename, lineno);
- up->sigclipmultipset=1;
- }
- else if(strcmp(name, "sigcliptolerance")==0)
- {
- if(up->sigcliptoleranceset) continue;
- gal_checkset_float_l_0_s_1(value, &p->sigcliptolerance, name,
- key, SPACK, filename, lineno);
- up->sigcliptoleranceset=1;
- }
-
-
- /* Operating modes: */
- /* Read options common to all programs */
- GAL_CONFIGFILES_READ_COMMONOPTIONS_FROM_CONF
-
-
- else
- error_at_line(EXIT_FAILURE, 0, filename, lineno,
- "`%s` not recognized.\n", name);
- }
-
- free(line);
- fclose(fp);
-}
-
-
-
-
-
-void
-printvalues(FILE *fp, struct subtractskyparams *p)
-{
- struct uiparams *up=&p->up;
- struct gal_mesh_params *mp=&p->mp;
- struct gal_commonparams *cp=&p->cp;
-
- /* Print all the options that are set. Separate each group with a
- commented line explaining the options in that group. */
- fprintf(fp, "\n# Input:\n");
- if(cp->hduset)
- {
- if(gal_checkset_string_has_space(cp->hdu))
- fprintf(fp, CONF_SHOWFMT"\"%s\"\n", "hdu", cp->hdu);
- else
- fprintf(fp, CONF_SHOWFMT"%s\n", "hdu", cp->hdu);
- }
- if(up->masknameset)
- {
- if(gal_checkset_string_has_space(up->maskname))
- fprintf(fp, CONF_SHOWFMT"\"%s\"\n", "mask", up->maskname);
- else
- fprintf(fp, CONF_SHOWFMT"%s\n", "mask", up->maskname);
- }
- if(up->mhdu)
- {
- if(gal_checkset_string_has_space(up->mhdu))
- fprintf(fp, CONF_SHOWFMT"\"%s\"\n", "mhdu", up->mhdu);
- else
- fprintf(fp, CONF_SHOWFMT"%s\n", "mhdu", up->mhdu);
- }
- if(up->kernelnameset)
- {
- if(gal_checkset_string_has_space(up->kernelname))
- fprintf(fp, CONF_SHOWFMT"\"%s\"\n", "kernel", up->kernelname);
- else
- fprintf(fp, CONF_SHOWFMT"%s\n", "kernel", up->kernelname);
- }
- if(up->khdu)
- {
- if(gal_checkset_string_has_space(up->khdu))
- fprintf(fp, CONF_SHOWFMT"\"%s\"\n", "khdu", up->khdu);
- else
- fprintf(fp, CONF_SHOWFMT"%s\n", "khdu", up->khdu);
- }
-
-
- fprintf(fp, "\n# Output:\n");
- if(cp->outputset)
- fprintf(fp, CONF_SHOWFMT"%s\n", "output", cp->output);
-
-
- fprintf(fp, "\n# Mesh grid:\n");
- if(up->meshsizeset)
- fprintf(fp, CONF_SHOWFMT"%zu\n", "meshsize", mp->meshsize);
- if(up->nch1set)
- fprintf(fp, CONF_SHOWFMT"%zu\n", "nch1", mp->nch1);
- if(up->nch2set)
- fprintf(fp, CONF_SHOWFMT"%zu\n", "nch2", mp->nch2);
- if(up->lastmeshfracset)
- fprintf(fp, CONF_SHOWFMT"%.3f\n", "lastmeshfrac", mp->lastmeshfrac);
- if(up->mirrordistset)
- fprintf(fp, CONF_SHOWFMT"%.3f\n", "mirrordist", mp->mirrordist);
- if(up->minmodeqset)
- fprintf(fp, CONF_SHOWFMT"%.3f\n", "minmodeq", mp->minmodeq);
- if(up->numnearestset)
- fprintf(fp, CONF_SHOWFMT"%zu\n", "numnearest", mp->numnearest);
- if(up->smoothwidthset)
- fprintf(fp, CONF_SHOWFMT"%zu\n", "smoothwidth", mp->smoothwidth);
- if(up->fullconvolutionset)
- fprintf(fp, CONF_SHOWFMT"%d\n", "fullconvolution",
- mp->fullconvolution);
- if(up->fullinterpolationset)
- fprintf(fp, CONF_SHOWFMT"%d\n", "fullinterpolation",
- mp->fullinterpolation);
- if(up->fullsmoothset)
- fprintf(fp, CONF_SHOWFMT"%d\n", "fullsmooth", mp->fullsmooth);
-
-
- fprintf(fp, "\n# Statistics:\n");
- if(up->sigclipmultipset)
- fprintf(fp, CONF_SHOWFMT"%.3f\n", "sigclipmultip", p->sigclipmultip);
- if(up->sigcliptoleranceset)
- fprintf(fp, CONF_SHOWFMT"%.3f\n", "sigcliptolerance",
- p->sigcliptolerance);
-
- /* For the operating mode, first put the macro to print the common
- options, then the (possible options particular to this
- program). */
- fprintf(fp, "\n# Operating mode:\n");
- GAL_CONFIGFILES_PRINT_COMMONOPTIONS;
-}
-
-
-
-
-
-
-/* Note that numthreads will be used automatically based on the
- configure time. */
-void
-checkifset(struct subtractskyparams *p)
-{
- struct uiparams *up=&p->up;
- struct gal_commonparams *cp=&p->cp;
-
- int intro=0;
- if(cp->hduset==0)
- GAL_CONFIGFILES_REPORT_NOTSET("hdu");
- if(up->khduset==0)
- GAL_CONFIGFILES_REPORT_NOTSET("khdu");
-
- /* Mesh grid: */
- if(up->meshsizeset==0)
- GAL_CONFIGFILES_REPORT_NOTSET("meshsize");
- if(up->nch1set==0)
- GAL_CONFIGFILES_REPORT_NOTSET("nch1");
- if(up->nch2set==0)
- GAL_CONFIGFILES_REPORT_NOTSET("nch2");
- if(up->lastmeshfracset==0)
- GAL_CONFIGFILES_REPORT_NOTSET("lastmeshfrac");
- if(up->mirrordistset==0)
- GAL_CONFIGFILES_REPORT_NOTSET("mirrordist");
- if(up->minmodeqset==0)
- GAL_CONFIGFILES_REPORT_NOTSET("minmodeq");
- if(up->numnearestset==0)
- GAL_CONFIGFILES_REPORT_NOTSET("numnearest");
- if(up->smoothwidthset==0)
- GAL_CONFIGFILES_REPORT_NOTSET("smoothwidth");
- if(up->fullconvolutionset==0)
- GAL_CONFIGFILES_REPORT_NOTSET("fullconvolution");
- if(up->fullinterpolationset==0)
- GAL_CONFIGFILES_REPORT_NOTSET("fullinterpolation");
- if(up->fullsmoothset==0)
- GAL_CONFIGFILES_REPORT_NOTSET("fullsmooth");
-
- /* Statistics: */
- if(up->sigclipmultipset==0)
- GAL_CONFIGFILES_REPORT_NOTSET("sigclipmultip");
- if(up->sigcliptoleranceset==0)
- GAL_CONFIGFILES_REPORT_NOTSET("sigcliptolerance");
-
- GAL_CONFIGFILES_END_OF_NOTSET_REPORT;
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-/**************************************************************/
-/*************** Sanity Check *******************/
-/**************************************************************/
-void
-sanitycheck(struct subtractskyparams *p)
-{
-
- /* Make sure the input file exists. */
- gal_checkset_check_file(p->up.inputname);
-
- /* Set the maskname and mask hdu accordingly: */
- gal_fits_file_or_ext_name(p->up.inputname, p->cp.hdu, p->up.masknameset,
- &p->up.maskname, p->up.mhdu, p->up.mhduset,
- "mask");
-
- /* Set the output name: */
- if(p->cp.output)
- gal_checkset_check_remove_file(p->cp.output, p->cp.dontdelete);
- else
- gal_checkset_automatic_output(p->up.inputname, "_skysubed.fits",
- p->cp.removedirinfo, p->cp.dontdelete,
- &p->cp.output);
-
- /* Set the sky image name: */
-
- /* Set the check image names: */
- if(p->meshname)
- {
- p->meshname=NULL; /* Was not allocated before! */
- gal_checkset_automatic_output(p->up.inputname, "_mesh.fits",
- p->cp.removedirinfo, p->cp.dontdelete,
- &p->meshname);
- }
- if(p->convname)
- {
- p->convname=NULL; /* Was not allocated before! */
- gal_checkset_automatic_output(p->up.inputname, "_conv.fits",
- p->cp.removedirinfo, p->cp.dontdelete,
- &p->convname);
- }
- if(p->skyname)
- {
- p->skyname=NULL; /* Was not allocated before! */
- gal_checkset_automatic_output(p->up.inputname, "_sky.fits",
- p->cp.removedirinfo, p->cp.dontdelete,
- &p->skyname);
- }
-
-
- /* Other checks: */
- if(p->mp.numnearest<GAL_MESH_MIN_ACCEPTABLE_NEAREST)
- error(EXIT_FAILURE, 0, "the smallest possible number for "
- "`--numnearest' (`-n') is %d. You have asked for: %zu",
- GAL_MESH_MIN_ACCEPTABLE_NEAREST, p->mp.numnearest);
-
- /* Set the constants in the gal_mesh_params structure. */
- p->mp.params=p;
- p->mp.numthreads=p->cp.numthreads;
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-/**************************************************************/
-/*************** Preparations *******************/
-/**************************************************************/
-void
-preparearrays(struct subtractskyparams *p)
-{
- struct gal_mesh_params *mp=&p->mp;
-
- /* Read the input image. */
- gal_fits_file_to_float(p->up.inputname, p->up.maskname, p->cp.hdu,
- p->up.mhdu, (float **)&p->mp.img, &p->bitpix,
- &p->anyblank, &mp->s0, &mp->s1);
- gal_fits_read_wcs(p->up.inputname, p->cp.hdu, 0, 0, &p->nwcs, &p->wcs);
-
- /* Read the kernel: */
- if(p->up.kernelnameset)
- gal_fits_prep_float_kernel(p->up.kernelname, p->up.khdu, &mp->kernel,
- &mp->ks0, &mp->ks1);
-
- /* Check if the input sizes and channel sizes are exact
- multiples. */
- if( mp->s0%mp->nch2 || mp->s1%mp->nch1 )
- error(EXIT_FAILURE, 0, "the input image size (%zu x %zu) is not an "
- "exact multiple of the number of the given channels (%zu, %zu) "
- "in the respective axis", mp->s1, mp->s0, mp->nch1, mp->nch2);
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-/**************************************************************/
-/************ Set the parameters *************/
-/**************************************************************/
-void
-setparams(int argc, char *argv[], struct subtractskyparams *p)
-{
- struct gal_commonparams *cp=&p->cp;
-
- /* Set the non-zero initial values, the structure was initialized to
- have a zero value for all elements. */
- cp->spack = SPACK;
- cp->verb = 1;
- cp->numthreads = num_processors(NPROC_CURRENT);
- cp->removedirinfo = 1;
-
- /* Read the arguments. */
- errno=0;
- if(argp_parse(&thisargp, argc, argv, 0, 0, p))
- error(EXIT_FAILURE, errno, "parsing arguments");
-
- /* Add the user default values and save them if asked. */
- GAL_CONFIGFILES_CHECK_SET_CONFIG;
-
- /* Check if all the required parameters are set. */
- checkifset(p);
-
- /* Print the values for each parameter. */
- if(cp->printparams)
- GAL_CONFIGFILES_REPORT_PARAMETERS_SET;
-
- /* Do a sanity check. */
- sanitycheck(p);
-
- /* Make the array of input images. */
- preparearrays(p);
-
- /* Everything is ready, notify the user of the program starting. */
- if(cp->verb)
- {
- printf(SPACK_NAME" started on %s", ctime(&p->rawtime));
- printf(" - Using %zu CPU threads.\n", p->cp.numthreads);
- printf(" - Input: %s (hdu: %s)\n", p->up.inputname, p->cp.hdu);
- if(p->up.maskname)
- printf(" - Mask: %s (hdu: %s)\n", p->up.maskname, p->up.mhdu);
- if(p->up.kernelnameset)
- printf(" - Kernel: %s (hdu: %s)\n", p->up.kernelname, p->up.khdu);
- }
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-/**************************************************************/
-/************ Free allocated, report *************/
-/**************************************************************/
-void
-freeandreport(struct subtractskyparams *p, struct timeval *t1)
-{
- /* Free the allocated arrays: */
- free(p->mp.img);
- free(p->cp.hdu);
- free(p->up.khdu);
- free(p->up.mhdu);
- free(p->cp.output);
-
- /* Free all the allocated names: */
- if(p->meshname) free(p->meshname);
- if(p->up.kernelnameset) free(p->up.kernelname);
-
- /* Free the mask image name. Note that p->up.inputname was not
- allocated, but given to the program by the operating system. */
- if(p->up.maskname && p->up.maskname!=p->up.inputname)
- free(p->up.maskname);
-
- /* Free the WCS structure: */
- if(p->wcs)
- wcsvfree(&p->nwcs, &p->wcs);
-
- /* Print the final message. */
- if(p->cp.verb)
- gal_timing_report(t1, SPACK_NAME" finished in: ", 0);
-}
diff --git a/bin/subtractsky/ui.h b/bin/subtractsky/ui.h
deleted file mode 100644
index 7cab110..0000000
--- a/bin/subtractsky/ui.h
+++ /dev/null
@@ -1,32 +0,0 @@
-/*********************************************************************
-SubtractSky - Find and subtract the sky value from an image.
-SubtractSky is part of GNU Astronomy Utilities (Gnuastro) package.
-
-Original author:
- Mohammad Akhlaghi <address@hidden>
-Contributing author(s):
-Copyright (C) 2015, Free Software Foundation, Inc.
-
-Gnuastro is free software: you can redistribute it and/or modify it
-under the terms of the GNU General Public License as published by the
-Free Software Foundation, either version 3 of the License, or (at your
-option) any later version.
-
-Gnuastro is distributed in the hope that it will be useful, but
-WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-General Public License for more details.
-
-You should have received a copy of the GNU General Public License
-along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
-**********************************************************************/
-#ifndef IMCROPUI_H
-#define IMCROPUI_H
-
-void
-setparams(int argc, char *argv[], struct subtractskyparams *p);
-
-void
-freeandreport(struct subtractskyparams *p, struct timeval *t1);
-
-#endif
diff --git a/configure.ac b/configure.ac
index 373d990..d642ba8 100644
--- a/configure.ac
+++ b/configure.ac
@@ -508,12 +508,6 @@ AC_ARG_ENABLE([statistics],
[AS_IF([test "x$enable_statistics" != xno],
[enable_statistics=yes; ayes=true])],
[enable_statistics=notset])
-AC_ARG_ENABLE([subtractsky],
- [AS_HELP_STRING([--enable-subtractsky],
- [Install SubtractSky and other enabled programs.])],
- [AS_IF([test "x$enable_subtractsky" != xno],
- [enable_subtractsky=yes; ayes=true])],
- [enable_subtractsky=notset])
AC_ARG_ENABLE([table],
[AS_HELP_STRING([--enable-table],
[Install Table and other enabled programs.])],
@@ -558,7 +552,6 @@ AS_IF([test $ayes = true ],
AS_IF([test $enable_mkprof = notset], [enable_mkprof=no])
AS_IF([test $enable_noisechisel = notset], [enable_noisechisel=no])
AS_IF([test $enable_statistics = notset], [enable_statistics=no])
- AS_IF([test $enable_subtractsky = notset], [enable_subtractsky=no])
AS_IF([test $enable_table = notset], [enable_table=no])
# AS_IF([test $enable_TEMPLATE = notset], [enable_TEMPLATE=no])
AS_IF([test $enable_warp = notset], [enable_warp=no])
@@ -575,7 +568,6 @@ AS_IF([test $ayes = true ],
AS_IF([test $enable_mkprof = notset], [enable_mkprof=yes])
AS_IF([test $enable_noisechisel = notset], [enable_noisechisel=yes])
AS_IF([test $enable_statistics = notset], [enable_statistics=yes])
- AS_IF([test $enable_subtractsky = notset], [enable_subtractsky=yes])
AS_IF([test $enable_table = notset], [enable_table=yes])
# AS_IF([test $enable_TEMPLATE = notset], [enable_TEMPLATE=yes])
AS_IF([test $enable_warp = notset], [enable_warp=yes])
@@ -598,7 +590,6 @@ AM_CONDITIONAL([COND_MKNOISE], [test $enable_mknoise = yes])
AM_CONDITIONAL([COND_MKPROF], [test $enable_mkprof = yes])
AM_CONDITIONAL([COND_NOISECHISEL], [test $enable_noisechisel = yes])
AM_CONDITIONAL([COND_STATISTICS], [test $enable_statistics = yes])
-AM_CONDITIONAL([COND_SUBTRACTSKY], [test $enable_subtractsky = yes])
AM_CONDITIONAL([COND_TABLE], [test $enable_table = yes])
#AM_CONDITIONAL([COND_TEMPLATE], [test $enable_TEMPLATE = yes])
AM_CONDITIONAL([COND_WARP], [test $enable_warp = yes])
@@ -627,7 +618,6 @@ AC_CONFIG_FILES([Makefile
bin/arithmetic/Makefile
bin/statistics/Makefile
bin/noisechisel/Makefile
- bin/subtractsky/Makefile
bootstrapped/lib/Makefile
bootstrapped/tests/Makefile
])
diff --git a/doc/Makefile.am b/doc/Makefile.am
index f75c829..f764c77 100644
--- a/doc/Makefile.am
+++ b/doc/Makefile.am
@@ -114,9 +114,6 @@ endif
if COND_STATISTICS
MAYBE_STATISTICS_MAN = man/aststatistics.1
endif
-if COND_SUBTRACTSKY
- MAYBE_SUBTRACTSKY_MAN = man/astsubtractsky.1
-endif
if COND_TABLE
MAYBE_TABLE_MAN = man/asttable.1
endif
@@ -130,7 +127,7 @@ dist_man_MANS = $(MAYBE_ARITHMETIC_MAN)
$(MAYBE_CONVERTT_MAN) \
$(MAYBE_CONVOLVE_MAN) $(MAYBE_COSMICCAL_MAN) $(MAYBE_CROP_MAN) \
$(MAYBE_FITS_MAN) $(MAYBE_WARP_MAN) $(MAYBE_MKCATALOG_MAN) \
$(MAYBE_MKNOISE_MAN) $(MAYBE_MKPROF_MAN) $(MAYBE_NOISECHISEL_MAN) \
- $(MAYBE_STATISTICS_MAN) $(MAYBE_SUBTRACTSKY_MAN) $(MAYBE_TABLE_MAN)
+ $(MAYBE_STATISTICS_MAN) $(MAYBE_TABLE_MAN)
## See if help2man is present or not. When help2man doesn't exist, we don't
@@ -192,10 +189,6 @@ man/aststatistics.1: $(top_srcdir)/bin/statistics/args.h
$(ALLMANSDEP)
$(MAYBE_HELP2MAN) -n "calculate statistics of a dataset" \
--libtool $(toputildir)/statistics/aststatistics
-man/astsubtractsky.1: $(top_srcdir)/bin/subtractsky/args.h $(ALLMANSDEP)
- $(MAYBE_HELP2MAN) -n "subtract the sky value in an image" \
- --libtool $(toputildir)/subtractsky/astsubtractsky
-
man/asttable.1: $(top_srcdir)/bin/table/args.h $(ALLMANSDEP)
$(MAYBE_HELP2MAN) -n "read and write ASCII and FITS tables" \
--libtool $(toputildir)/table/asttable
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 9717e7f..2c68e90 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -364,7 +364,6 @@ Data manipulation
* Arithmetic:: Arithmetic on input data.
* Convolve:: Convolve an image with a kernel.
* Warp:: Warp/Transform an image to a different grid.
-* SubtractSky:: Find the sky and subtract it from image.
Crop
@@ -417,24 +416,6 @@ Warp
* Resampling:: Warping an image is re-sampling it.
* Invoking astwarp:: Arguments and options for Warp.
-SubtractSky
-
-* Sky value:: Definition of the sky value.
-* Tiling an image:: Defining a mesh grid on the image.
-* Invoking astsubtractsky:: Options and arguments to SubtractSky.
-
-Sky value
-
-* Finding the sky value:: How SubtractSky finds the sky value.
-* Sky value misconceptions:: Sky value misconceptions and wrong approaches.
-
-Tiling an image
-
-* Quantifying signal in a mesh:: Good meshs for gradients.
-* Grid interpolation and smoothing:: Filling the blank grid elements.
-* Checking grid values:: Check the values on the grid.
-* Mesh grid options:: Common options related to the mesh grid.
-
Data analysis
* Statistics:: Calculate dataset statistics.
@@ -445,8 +426,15 @@ Statistics
* Histogram and Cumulative Frequency Plot:: Basic definitions.
* Sigma clipping:: Definition of @mymath{\sigma}-clipping
+* Sky value::
* Invoking aststatistics:: Arguments and options to Statistics.
+Sky value
+
+* Sky value definition:: Definition of the Sky/reference value.
+* Sky value misconceptions:: Wrong methods to estimate the Sky value.
+* Quantifying signal in a tile:: Method to estimate the presence of signal.
+
NoiseChisel
* Invoking astnoisechisel:: Options and arguments for NoiseChisel.
@@ -2214,16 +2202,40 @@ night's measurements on the ecliptic.
@node Installation, Common program behavior, Tutorials, Top
@chapter Installation
address@hidden This link is put here because the `Quick start' section of the
first
address@hidden chapter is not the most eye-catching part of the manual and some
users
address@hidden were seen to follow this ``Installation'' chapter title in
search of the
address@hidden tarball and fast instructions.
@cindex Installation
-To successfully install and thus use Gnuastro, first please make sure you
-have the dependencies installed (see @ref{Dependencies}). Only the first
-group are required when you are building Gnuastro using a tarball (see
address@hidden tarball}), they are very basic and low-level tools used in
-most astronomical software, so you might already have them installed, if
-not they are very easy to install as described for each. @ref{Downloading
-the source} discusses the two methods you can obtain the source code: as a
-tarball (a significant snapshot in Gnuastro's history), or the full
-history.
+The latest released version of Gnuastro source code is always available at
+the following URL:
+
address@hidden://ftpmirror.gnu.org/gnuastro/gnuastro-latest.tar.gz}
+
address@hidden
address@hidden start} discribes the commands necessary to configure, build, and
+install Gnuastro on your system. This chapter will be useful in cases where
+the simple procedure above is not sufficient, for example your system lacks
+a mandatory/optional dependency (in other words, you can't pass the
address@hidden ./configure} step), or you want greater customization, or you
+want to build and install Gnuastro from other random points in its history,
+or you want a higher level of control on the installation. Thus if you were
+happy with downloading the tarball and following @ref{Quick start}, then
+you can safely ignore this chapter and come back to it in the future if you
+need more customization.
+
address@hidden describes the mandatory, optional and bootstrapping
+dependencies of Gnuastro. Only the first group are required/mandatory when
+you are building Gnuastro using a tarball (see @ref{Release tarball}), they
+are very basic and low-level tools used in most astronomical software, so
+you might already have them installed, if not they are very easy to install
+as described for each. @ref{Downloading the source} discusses the two
+methods you can obtain the source code: as a tarball (a significant
+snapshot in Gnuastro's history), or the full
address@hidden@ref{Bootstrapping dependencies} are required if you clone
+the full history.}. The latter allows you to build Gnuastro at any random
+point in its history (for example to get bug fixes or new features that are
+not released as a tarball yet).
The building and installation of Gnuastro is heavily customizable, to learn
more about them, see @ref{Build and install}. This section is essentially a
@@ -4486,11 +4498,29 @@ though. You can see how many pixels were within each
tile (for example to
weight the values or discard some for later processing) with Gnuastro's
Statistics (see @ref{Statistics}) as shown below. The output FITS file is
going to have two extensions, one with the median calculated on each tile
-and one with the number of elements that each tile covers.
+and one with the number of elements that each tile covers. You can then use
+the @code{where} operator in @ref{Arithmetic} to set the values of all
+tiles that don't have the regular area to a blank value.
@example
-$ aststatistics --median --number --ontile input.fits
+$ aststatistics --median --number --ontile input.fits \
+ --oneelempertile --output=o.fits
+$ REGULAR_AREA=1600 # Check second extension of `o.fits'.
+$ astarithmetic o.fits o.fits $REGULAR_AREA ne nan where \
+ -h1 -h2
@end example
+
+Note that if @file{input.fits} also has blank values, then the median on
+tiles with blank values will also be ignored with the command above (which
+is desirable).
+
+
address@hidden --inteponlyblank
+When values are to be interpolated, only change the values of the blank
+elements, keep the non-blank elements untouched.
+
address@hidden --interpnumngb=INT
+The number of nearby non-blank neighbors to use for interpolation.
@end table
@node Operating mode options, , Processing options, Common options
@@ -5265,7 +5295,7 @@ speicy a type. For example, as a value to the common
option @option{--type}
@item i32
@itemx int32
-32-bit signed integers, range:@* @mymath{[-2^{32}\rm{\ to\ }2^{32}-1]} or
+32-bit signed integers, range:@* @mymath{[-2^{31}\rm{\ to\ }2^{31}-1]} or
@mymath{[-2147483648\rm{\ to\ }2147483647]}.
@item u64
@@ -5275,7 +5305,7 @@ speicy a type. For example, as a value to the common
option @option{--type}
@item i64
@itemx int64
-64-bit signed integers, range:@* @mymath{[-2^{64}\rm{\ to\ }2^{64}-1]} or
+64-bit signed integers, range:@* @mymath{[-2^{63}\rm{\ to\ }2^{63}-1]} or
@mymath{[-9223372036854775808\rm{\ to\ }9223372036854775807]}.
@item f32
@@ -5283,9 +5313,9 @@ speicy a type. For example, as a value to the common
option @option{--type}
32-bit (single-precision) floating point types. The maximum (minimum is its
negative) possible value is
@mymath{3.402823\times10^{38}}. Single-precision floating points can
-accurately represent a floating point number @mymath{\sim7.2} significant
-decimals. Given the heavy noise in astronomical data, this is usually more
-than sufficient for storing results.
+accurately represent a floating point number upto @mymath{\sim7.2}
+significant decimals. Given the heavy noise in astronomical data, this is
+usually more than sufficient for storing results.
@item f64
@itemx float64
@@ -5300,12 +5330,13 @@ floating points can accurately represent a floating
point number
@cartouche
@noindent
@strong{Some file formats don't recognize all types.} Some file formats
-don't recognize all the types, for example the FITS binary table standard
-does not define the @code{unsigned int} and @code{unsigned long} types. The
-FITS standard also only defines the following types for images:
address@hidden char}, @code{short}, @code{int}, @code{long}, @code{float},
-and @code{double}. Although CFITSIO allows writing/reading to/from a wider
-range of types.
+don't recognize all the types, for example the FITS standard (see
address@hidden) does not define @code{uint64} in binary tables or images. When
+a type is not acceptable for output into a given file format, the
+respective Gnuastro program or library will let you know and abort. On the
+command-line, you can use the @ref{Arithmetic} program to convert the
+numerical type of a dataset, in the libraries, you can call
address@hidden
@end cartouche
@@ -5724,15 +5755,16 @@ you can define a tile grid over the image. When the
tile sizes are set
properly, the background's variation over each tile will be negligble,
allowing you to measure (and subtract) it. In other cases (for example
spatial domain convolution in Gnuastro, see @ref{Convolve}), it might
-simply be for speed of processing, since each tile can be processed
-independently on a separate CPU thread. In the arts and mathematics, this
-process is formally known as
address@hidden://en.wikipedia.org/wiki/Tessellation, tessellation}.
+simply be for speed of processing: each tile can be processed independently
+on a separate CPU thread. In the arts and mathematics, this process is
+formally known as @url{https://en.wikipedia.org/wiki/Tessellation,
+tessellation}.
The size of the regular tiles (in units of data-elements, or pixels in an
image) can be defined with the @option{--tilesize} option. It takes
multiple numbers (separated by a comma) which will be the length along the
-respective dimension (in Fortran/FITS order). For example
+respective dimension (in Fortran/FITS dimension order). Divisions are also
+acceptable, but must result in an integer. For example
@option{--tilesize=30,40} can be used for an image (a 2D dataset). The
regular tile size along the first FITS axis (horizontal when viewed in SAO
ds9) will be 30 pixels and along the second it will be 40 pixels. Ideally,
@@ -5750,16 +5782,19 @@ less than the size along that dimension), the remainder
will be added to
one regular tile's size and the large tile will be cut in half and put in
the two ends of the grid/tessellation. In this way, all the tiles in the
central regions of the dataset will have the regular tile sizes and the
-tiles on the edge will be slightly smaller. You can set the significance
-fraction with the @option{--remainderfrac} option.
+tiles on the edge will be slightly larger/smaller depending on the
+remainder significance. The fraction which defines the remainder
+significance along all dimensions can be set through
address@hidden
The best tile size is directly related to the spatial properties of the
-gradient on the image. In practice we assume that the gradient is not
-present over each tile. So if there is a strong gradient (for example in
-long wavelength ground based images) or the image is of a crowded area
-where there isn't too much blank area, you have to choose a smaller tile
-size. A larger mesh will give more pixels and so the scatter in the results
-will be less.
+property you want to study (for example, gradient on the image). In
+practice we assume that the gradient is not present over each tile. So if
+there is a strong gradient (for example in long wavelength ground based
+images) or the image is of a crowded area where there isn't too much blank
+area, you have to choose a smaller tile size. A larger mesh will give more
+pixels and and so the scatter in the results will be less (better
+statistics).
@cindex CCD
@cindex Amplifier
@@ -5788,24 +5823,26 @@ systematic biases will then propagate to all subsequent
measurements we do
on the data (for example photometry and subsequent stellar mass and star
formation rate measurements in the case of galaxies).
-Therefore an accurate analysis requires a lower-level tesselation: with
-larger tiles, each large tile covering one amplifier channel. For clarity
-we'll call these larger tiles ``channels''. The number of channels along
-each dimension is defined through the @option{--numchannels}. Each channel
-will then be covered by its own individual smaller tessellation (with tile
-sizes determined by the @option{--tilesize} option). This will allow two
-adjacent pixels from different channels to never be considered together in
-the higher-level analysis. If the image is processed or the detector only
-has one amplifier, you can set the number of channels in both axes to 1.
-
-The final tesselation can be seen on the image with the
+Therefore an accurate analysis requires a two layer tesselation: the top
+layer contains larger tiles, each covering one amplifier channel. For
+clarity we'll call these larger tiles ``channels''. The number of channels
+along each dimension is defined through the @option{--numchannels}. Each
+channel is then covered by its own individual smaller tessellation (with
+tile sizes determined by the @option{--tilesize} option). This will allow
+independent analyis of two adjacent pixels from different channels if
+necessary. If the image is processed or the detector only has one
+amplifier, you can set the number of channels in both dimension to 1.
+
+The final tesselation can be inspected on the image with the
@option{--checktiles} option that is available to all programs which use
tessellation for localized operations. When this option is called, a FITS
-file with a @file{_mesh.fits} suffix will be created along with the
-outputs, see @ref{Automatic output}. Each pixel in this image has the ID of
-the tile that covers it. If the number of channels in any dimension are
-larger than unity, you will notice that the tile IDs are defined such that
-the first channels is covered first, then the second and so on.
+file with a @file{_tiled.fits} suffix will be created along with the
+outputs, see @ref{Automatic output}. Each pixel in this image has the
+number of the tile that covers it. If the number of channels in any
+dimension are larger than unity, you will notice that the tile IDs are
+defined such that the first channels is covered first, then the second and
+so on. For the full list of processing-related common options (including
+tessellation options), please see @ref{Processing options}.
@@ -6295,11 +6332,13 @@ END
The most low-level and basic property of a dataset is how it is stored. To
process, archive and transmit the data, you need a container to store it
first. From the start of the computer age, different formats have been
-difined to store data, optimized for pariticular applications, one
-format/container can never be useful for all applications. In astronomy,
-FITS is the most common format of data storage and transmission. It has
-many useful features, for example multiple sub-containers (also known as
-extensions or header data units, HDUs) within one file. Each HDU can store
+defined to store data, optimized for pariticular applications. One
+format/container can never be useful for all applications: the storage
+defines the application and vice-versa. In astronomy, the Flexible Image
+Transport System (FITS) standard has become the most common format of data
+storage and transmission. It has many useful features, for example multiple
+sub-containers (also known as extensions or header data units, HDUs) within
+one file, or support for tables as well as images. Each HDU can store an
independent dataset and its corresponding meta-data. Therefore, Gnuastro
has one program (see @ref{Fits}) specifically designed to manipulate FITS
HDUs and the meta-data (header keywords) in each HDU.
@@ -7440,7 +7479,6 @@ transformation to it.
* Arithmetic:: Arithmetic on input data.
* Convolve:: Convolve an image with a kernel.
* Warp:: Warp/Transform an image to a different grid.
-* SubtractSky:: Find the sky and subtract it from image.
@end menu
@node Crop, Arithmetic, Data manipulation, Data manipulation
@@ -9169,7 +9207,7 @@ of Principia) showed that:
@mymath{ix=\ln(\cos{x}+i\sin{x})}.} (1707 --
@mymath{e^{iv+2\pi}=e^{iv}}. Later, Caspar Wessel (mathematician and
cartographer 1745 -- 1818 A.D.) showed how complex numbers can be
displayed as vectors on a plane. Euler's identity might seem counter
-intuitive at first, so we will try to explain it geometrically (for a more
+intuitive at first, so we will try to explain it geometrically (for deeper
physical insight). On the real-imaginary 2D plane (like the left hand plot
in each box of @ref{iandtime}), multiplying a number by @mymath{i} can be
interpretted as rotating the point by @mymath{90} degrees (for example the
@@ -10204,7 +10242,7 @@ explanations under the @option{--makekernel} option for
more information.
address@hidden Warp, SubtractSky, Convolve, Data manipulation
address@hidden Warp, , Convolve, Data manipulation
@section Warp
Image warping is the process of mapping the pixels of one image onto a new
pixel grid. This process is sometimes known as transformation, however
@@ -10790,579 +10828,11 @@ of the pixels in the input and output images will be
the same).
address@hidden SubtractSky, , Warp, Data manipulation
address@hidden SubtractSky
-
address@hidden Atmosphere
address@hidden Flat field
address@hidden Stray light
address@hidden Bias subtraction
-Raw astronomical images (and even poorly processed images) don't usually
-have a uniform `sky' value over their surface prior to processing, see
address@hidden value} for a complete definition of the sky value. However, a
-uniform sky value over the image is vital for further processing. For
-ground based images (particularly at longer wavelengths) this can be due to
-actual variations in the atmosphere. Another cause might be systematic
-biases in the instrument or prior processing. For example stray light in
-the telescope/camera or bad flat-fielding or bias subtraction. The latter
-is a major issue in space based images where the atmosphere is no longer a
-problem.
-
address@hidden Grid
address@hidden Mesh
address@hidden Gradient
-SubtractSky is a tool to find the sky value and its standard deviation
-on a grid over the image. The size of the grid will determine how
-accurately it can account for gradients in the sky value. Such that if
-the gradient (change of sky value) is too sharp, a smaller grid size
-has to be chosen. However the results will be most accurate with a
-larger grid size.
-
-
address@hidden
-* Sky value:: Definition of the sky value.
-* Tiling an image:: Defining a mesh grid on the image.
-* Invoking astsubtractsky:: Options and arguments to SubtractSky.
address@hidden menu
-
address@hidden Sky value, Tiling an image, SubtractSky, SubtractSky
address@hidden Sky value
-
address@hidden Sky value
-The discussion here is taken from Akhlaghi and Ichikawa
-(2015)@footnote{See the section on sky in Akhlaghi M.,
-Ichikawa. T. (2015). Astrophysical Journal Supplement Series.}. Let's
-assume that all instrument defects -- bias, dark and flat -- have been
-corrected and the brightness (see @ref{Flux Brightness and magnitude})
-of a detected object, @mymath{O}, is desired. The sources of flux on
-pixel @address@hidden this analysis the dimension of the
-data (image) is irrelevant. So if the data is an image (2D) with width
-of @mymath{w} pixels, then a pixel located on column @mymath{x} and
-row @mymath{y} (where all counting starts from zero and (0, 0) is
-located on the bottom left corner of the image), would have an index:
address@hidden of the image can be written as follows:
-
address@hidden
address@hidden
-Contribution from the target object, (@mymath{O_i}).
address@hidden
-Contribution from other detected objects, (@mymath{D_i}).
address@hidden
-Undetected objects or the fainter undetected regions of bright
-objects, (@mymath{U_i}).
address@hidden
address@hidden Cosmic rays
-A cosmic ray, (@mymath{C_i}).
address@hidden
address@hidden Background flux
-The background flux, which is defined to be the count if none of the
-others exists on that pixel, (@mymath{B_i}).
address@hidden itemize
address@hidden
-The total flux in this pixel, @mymath{T_i}, can thus be written as:
-
address@hidden
-
address@hidden Cosmic ray removal
address@hidden
-By definition, @mymath{D_i} is detected and it can be assumed that it
-is correctly subtracted, so that @mymath{D_i} can be set to
-zero. There are also methods to detect and remove cosmic rays (for
-example, van Dokkum (2001)@footnote{van Dokkum,
-P. G. (2001). Publications of the Astronomical Society of the Pacific.
-113, 1420.}) enabling us to set @mymath{C_i=0}. Note that in practice,
address@hidden and @mymath{U_i} are correlated, because they both
-directly depend on the detection algorithm and its input
-parameters. Also note that no detection or cosmic ray removal
-algorithm is perfect. With these limitations in mind, the observed sky
-value for this pixel (@mymath{S_i}) can be defined as
-
address@hidden Sky value
address@hidden
-
address@hidden
-Therefore, as the detection process (algorithm and input parameters)
-becomes more accurate, or @mymath{U_i\to0}, the sky value will tend to
-the background value or @mymath{S_i\to B_i}. Therefore, while
address@hidden is an inherent property of the data (pixel in an image),
address@hidden depends on the detection process. Over a group of pixels,
-for example in an image or part of an image, this equation translates
-to the average of undetected pixels. With this definition of sky, the
-object flux in the data can be calculated with
-
address@hidden T_{i}=S_{i}+O_{i} \quad\rightarrow\quad
- O_{i}=T_{i}-S_{i}.}
-
address@hidden
address@hidden photo-electrons
-Hence, the more accurately @mymath{S_i} is measured, the more
-accurately the brightness (sum of pixel values) of the target object
-can be measured (photometry). Similarly, any under-(over-)estimation
-in the sky will directly translate to an over(under) estimation of the
-measured object's brightness. In the fainter outskirts of an object a
-very small fraction of the photo-electrons in the pixels actually
-belong to objects. Therefore even a small over estimation of the sky
-value will result in the loss of a very large portion of most
-galaxies. Based on the definition above, the sky value is only
-correctly found when all the detected objects (@mymath{D_i} and
address@hidden) have been removed from the data.
-
-
address@hidden
-* Finding the sky value:: How SubtractSky finds the sky value.
-* Sky value misconceptions:: Sky value misconceptions and wrong approaches.
address@hidden menu
-
-
-
address@hidden Finding the sky value, Sky value misconceptions, Sky value, Sky
value
address@hidden Finding the sky value
-
address@hidden Data
address@hidden Distribution mode
address@hidden Distribution mean
address@hidden Distribution median
address@hidden Mode of distribution
address@hidden Mean of distribution
address@hidden Median of distribution
-This technique to find the sky value in a distribution was initially
-proposed in Akhlaghi and Ichikawa address@hidden M.,
-Ichikawa. T. (2015). Astrophysical Journal Supplement Series.}.
-
-Note that through the difference of the mode and median we have
-actually `detected' data in the distribution. However this detection
-was only based on the total distribution of the data, not its spatial
-position in each mesh. So we adhere to the definition of Sky value in
address@hidden value}. Finding the median is very easy, the main problem is
-in finding the mode through a robust method. In Appendix C of Akhlaghi
-and Ichikawa (2015) a new approach to finding the mode in any
-astronomically relevant distribution is introduced.
-
address@hidden Cosmic rays
-Even when the mode and median are approximately equal, Cosmic rays can
-significantly bias the calculation of the average. Even if they are
-very few. However, usually, Cosmic rays have very sharp boundaries and
-do not fade away into the noise. Therefore, when the histogram of the
-distribution is plotted, they are clearly separate from the rest of
-the data. For example see Figure 15 in Akhlaghi and Ichikawa
-(2015). In such cases, @mymath{\sigma}-clipping is a perfect tool to
-remove the effect of such objects in the average and standard
-deviation. See @ref{Sigma clipping} for a complete explanation. So
-after asserting that the mode and median are approximately equal in a
-mesh (see @ref{Tiling an image}), convergence-based
address@hidden is also applied before getting the final sky
-value and its standard deviation for a mesh.
-
-
address@hidden Sky value misconceptions, , Finding the sky value, Sky value
address@hidden Sky value misconceptions
-
-As defined in @ref{Sky value}, the sky value is only accurately
-defined when the detection algorithm is not significantly reliant on
-the sky value. In particular its detection threshold. However, most
-signal-based detection tools @footnote{According to Akhlaghi and
-Ichikawa (2015), signal-based detection is a detection process that
-relies heavily on assumptions about the to-be-detected objects. This
-method was the most heavily used technique prior to the introduction
-of NoiseChisel in that paper.} used the sky value as a reference to
-define the detection threshold. So these old techniques had to rely on
-approximations based on other assumptions about the data. A review of
-those other techniques can be seen in Appendix A of Akhlaghi and
-Ichikawa (2015)@footnote{Akhlaghi M.,
-Ichikawa. T. (2015). Astrophysical Journal Supplement Series.}. Since
-they were extensively used in astronomical data analysis for several
-decades, such approximations have given rise to a lot of
-misconceptions, ambiguities and disagreements about the sky value and
-how to measure it. As a summary, the major methods used until now were
-an approximation of the mode of the image pixel distribution and
address@hidden
-
address@hidden
address@hidden Histogram
address@hidden Distribution mode
address@hidden Mode of a distribution
address@hidden Probability density function
address@hidden
-To find the mode of a distribution those methods would either have to
-assume (or find) a certain probability density function (PDF) or use
-the histogram. But the image pixels can have any distribution, and the
-histogram results are very inaccurate (there is a large dispersion)
-and depend on bin-widths.
-
address@hidden sigma-clipping
address@hidden
-Another approach was to iteratively clip the brightest pixels in the
-image (which is known as @mymath{\sigma}-clipping, since the reference
-was found from the image mean and its standard deviation or
address@hidden). See @ref{Sigma clipping} for a complete
-explanation. The problem with @mymath{\sigma}-clipping was that real
-astronomical objects have diffuse and faint wings that penetrate
-deeply into the noise. So only removing their brightest parts is
-completely useless in removing the systematic bias an object's fainter
-parts cause in the sky value.
address@hidden itemize
-
-As discussed in @ref{Sky value}, the sky value can only be correctly
-defined as the average of undetected pixels. Therefore all such
-approaches that try to approximate the sky value prior to detection
-are ultimately poor approximations.
-
-
-
-
address@hidden Tiling an image, Invoking astsubtractsky, Sky value, SubtractSky
address@hidden Tiling an image
-
address@hidden Grid
address@hidden Tile
address@hidden Gradient
address@hidden Mesh grid
address@hidden Tile grid
address@hidden Tessellation
-
-
address@hidden
-* Quantifying signal in a mesh:: Good meshs for gradients.
-* Grid interpolation and smoothing:: Filling the blank grid elements.
-* Checking grid values:: Check the values on the grid.
-* Mesh grid options:: Common options related to the mesh grid.
address@hidden menu
-
-
address@hidden Quantifying signal in a mesh, Grid interpolation and smoothing,
Tiling an image, Tiling an image
address@hidden Quantifying signal in a mesh
-
address@hidden Noise
address@hidden Gaussian distribution
-Noise is characterized with a fixed background value and a certain
-distribution. For example, for the Gaussian distribution these two are
-the mean and standard deviation. When we have absolutely no signal and
-only noise in a data set, the mean, median and mode of the distribution
-are equal within statistical errors and approximately equal to the
-background value. For the next paragraph, let's assume that the
-background is subtracted and is zero.
-
-Data always has a positive value and will never become negative, see
-Figure 1 in Akhlaghi and Ichikawa (2015). Therefore, as data is buried
-into the noise, the mean, median and mode shift to the positive. The
-mean is the fastest in this shift. The median is slower since it is
-defined based on an ordered distribution and so is not affected by a
-small (less than half) number of outliers. Finally, the mode is the
-slowest to shift to the positive.
-
-Inverting the argument above provides us with the basis of Gnuastro's
-algorithm to quantify the presence of signal in a mesh. Namely, when
-the mode and median of a distribution are approximately equal, we can
-argue that there is no significant signal in that mesh. So we can
-consider the image to be made of a grid and use this argument to
-`detect' signal in each grid element. The median is defined to be the
-value of the 0.5 quantile in the image. So the only necessary
-parameter is the minimum acceptable quantile (smaller than 0.5) for
-the mode in a mesh that we deem accurate. See @ref{Mesh grid options}
-for an explanation of the options used to customize this behavior.
-
-Since there is sufficient signal in the mesh to bias the analysis on
-that mesh, any grid element whose mode quantile is smaller than the
-minimum acceptable quantile is usually kept with a blank value and no
-value is given to it. Finally, when all the grid elements have been
-checked, we can interpolate over all the empty elements and smooth the
-final result to find the sky value over the full image. See @ref{Grid
-interpolation and smoothing}.
-
-Convolving a data set (that contains signal and noise), creates a
-positive skewness in it depending on the fraction of data present in
-the distribution and also the convolution kernel. See Section 3.1.1 in
-Akhlaghi and Ichikawa (2015) and @ref{Convolution process}. This
-skewness can be interpreted as an increase in the Signal to noise
-ratio of the objects buried in the noise. Therefore, to obtain an
-even better measure of the presence of signal in a mesh, the image can
-be convolved with a given PSF first. This positive skew will result in
-more distance between the mode an median thereby enabling a more
-accurate detection of fainter signal, for example the faint wings of
-bright stars and galaxies. When convolving over a mesh grid, the
-pixels in each channel will be treated independently. This can be
-disabled with the @option{--fullconvolution} option. See
address@hidden kernel} for the respective options.
-
address@hidden Grid interpolation and smoothing, Checking grid values,
Quantifying signal in a mesh, Tiling an image
address@hidden Grid interpolation and smoothing
-
-On some of the grid elements, the desired value will not be found, for
-example the Sky value in @ref{Finding the sky value}. This can happen
-for lots of reasons depending on the job that was to be done on a
-mesh, for example when a large galaxy is present in the image, no Sky
-value will be found for the grid elements that lie over it. However,
-the Sky value should be `guessed' over that part of the image. We
-cannot just ignore those regions! To fill in such blank grid elements,
-we use interpolation.
-
address@hidden Spline interpolation
address@hidden Linear interpolation
address@hidden Bicubic interpolation
address@hidden Interpolation, spline
address@hidden Interpolation, bicubic
address@hidden Interpolation, bi-linear
-Parametric interpolations like bi-linear, bicubic or spline
-interpolations are not used because they fail terribly on the edges of
-the image. For example see Figure 16 in Akhlaghi and Ichikawa
-(2015). They are also prone to cause significant gradients in the
-image. So to find the interpolated value for each grid element,
-Gnuastro will look at a certain number of the nearest neighbors. The
-exact number can be specified through @option{--numnearest}. The
-median value of those grid elements will be taken as the final value
-for each mesh. The median is chosen because on the fainter wings of
-bright objects, the mean can easily become biased. If the number of
-meshes with a good value is less than the value given to
address@hidden, then the program will abort and notify you. In
-such cases you can either decrease the value to this option or set
-less restrictive requirements (for example a smaller
address@hidden, or larger/smaller meshs) at the expense of less
-accurate results.
-
address@hidden PSF
-By default the process above will occur on all the grid elements, not
-just the ones that are blank. This is done to avoid biased results on
-the faint wings of bright galaxies and stars (the PSF). Because their
-flux penetrates into the noise very slowly, it might not be possible
-to completely identify that flux and ignore that mesh. Therefore, if
-interpolation is only done on blank pixels, such false positives can
-cause a bias in the vicinity of bright objects, particularly after
-smoothing in the next step. In order to only interpolate over blank
-meshs and leave the values of the successful meshs untouched, the
address@hidden option can be used.
-
address@hidden Average filter
-Once all the grid elements are filled, the values given to all the
-meshs should be smoothed. This is because the median was used in the
-interpolation. The median is robust in the face of outliers, but there
-might be strong differences from one grid element to the next. By
-smoothing the grid, the variation between grid element to grid element
-will be far less. To smooth the mesh values, Gnuastro uses an average
-filter. An average filter is just a convolution of the mesh grid
-(spatial convolution with edge correction) by a kernel with all
-elements having an equal weight, see @ref{Convolution process}. The
-kernel is a square. The length of the kernel edge can be set in units
-of meshs through the @option{--smoothwidth} option (which has to be an
-odd number as with any kernel). If a width of 1 is given for the
-kernel width, then no smoothing will take place.
-
-By default the interpolation and smoothing explained above are done
-independently on each channel. In some circumstances, it might be
-preferable to do either one of these two steps on the whole image,
-independent of the channels. For example when there are gradients over
-the image and their variation over the image is stronger than the
-variation caused by the channels. Through the two options
address@hidden and @option{--fullsmooth} you can ask for
-interpolation or smoothing to use all the meshs in the image, not just
-those in the same channel of a mesh. Note that it is still very
-important that no mesh contain pixels from two channels. Since the
-pixels have been assigned to a mesh prior to these steps (see
address@hidden an image}), there is no problem in this regard.
-
-Even after smoothing, a simple visual examination of the values given
-to the pixels in each mesh over the image will give a very `boxy' or
-pixelated impression. To our eyes it feels that it would be much
-better if the values could be smoothed in sub-mesh (pixel) scales to
-obtain a smooth variation. An example is the results of bi-linear and
-bi-cubic interpolations, see Figure 16 in Akhlaghi and Ichikawa (2015)
-for several examples. The reasons we are not doing that level of
-smoothing are two fold:
-
address@hidden
-
address@hidden
-The difference in value between neighboring meshs is not statistically
-significant. After smoothing, the variation between the neighboring
-mesh values will be very small. It does visually appear to be a lot
-when viewing in image viewers like SAO ds9. This is because such
-viewers use the minimum and maximum range of all the pixel values as a
-reference to choose the color given to each pixel. However compared to
-the standard deviation of the noise in the image, the different values
-of each mesh are completely negligible. Please confirm this for your
-self on your own datasets to clearly understand.
-
address@hidden
-Such pixel-based (not mesh-based) interpolation will be very time
-consuming. Since the difference between the neighboring meshs is
-statistically negligible, it is simply not worth the time investment.
-
address@hidden itemize
-
address@hidden Checking grid values, Mesh grid options, Grid interpolation and
smoothing, Tiling an image
address@hidden Checking grid values
-
-Programs that use the mesh grid to calculate some value, also have
-checking options (that start with @option{--check}). When such options
-are called, the program will create a specific output file depending
-on what you want to check. When the operation is done on the meshes,
-this file shows the values assigned to each mesh grid element. It has
-three extensions:
-
address@hidden
-
address@hidden
-The value that was initially calculated for each grid element. If a
-mesh was not successful in providing a value, a NaN value is stored
-for that mesh.
-
address@hidden
-The interpolated values, see @ref{Grid interpolation and smoothing}.
-
address@hidden
-The smoothed values, see @ref{Grid interpolation and smoothing}.
-
address@hidden enumerate
-If more than one value is calculated for each mesh (for example the
-mean and its standard deviation), then each step has that many
-extensions. By default, these check image outputs have the same size
-(pixels) as the input image. So all the pixels within a mesh are given
-the same value. This is useful if you want to later apply these values
-to the image through another program for example. There is a one to
-one correspondence between the input image pixels and the grid showing
-you the mesh values for that pixel.
-
-In other situations, the input pixel resolution might not be important
-for you and you just want to see the relative mesh values over the
-image. In such cases, you can call the @option{--meshbasedcheck}
-option so the check image only has one pixel for each mesh. This image
-will only be as big as the full mesh grid and there will be no world
-coordinate system. When the input images are really large, this can
-make a difference in both the processing of the programs and in
-viewing the images.
-
-Another case when only one pixel for each mesh will be useful is when
-you might want to display the mesh values in a document and you don't
-want the volume of the document (in bytes) to get too high. For
-example if the mesh sizes are 30 pixels by 30 pixels, then the mesh
-grid created through this option will take @mymath{30\times30=900}
-times less space! You can resize the standard output, but the borders
-between the meshs will be blurred. The best case is to call that
-program with this option.
-
address@hidden Mesh grid options, , Checking grid values, Tiling an image
address@hidden Mesh grid options
-
address@hidden @option
-
address@hidden -d FLT
address@hidden --mirrordist=FLT
-The distance beyond the mirror point (in units of the error in the mirror
-point) to check for finding the mode in each mesh. This is part of the
-process to quantify the presence of signal in a mesh, see @ref{Quantifying
-signal in a mesh}. See appendix C in Akhlaghi and Ichikawa (2015) for a
-complete explanation of the mode-finding algorithm. The value to this
-option is shown as @mymath{\alpha} in that appendix.
-
address@hidden -Q FLT
address@hidden --minmodeq=FLT
-The minimum acceptable quantile for the mode of each mesh. The median is on
-the 0.5 quantile of the image and as long as we have positive signal (all
-astronomically relevant observations), the mode will be less than the
-median. The sky value is only found on meshes where the median and mode are
-approximately the same, see @ref{Sky value}.
-
address@hidden --interponlyblank
-Only interpolate the blank pixels. By default, interpolation will
-happen on all the mesh grids, not just the blank ones. See @ref{Grid
-interpolation and smoothing}.
-
address@hidden -n INT
address@hidden --numnearest=INT
-The number of nearest grid elements with a successful sky value to use for
-interpolating over blank mesh elements (those that had a significant
-contribution of signal), see @ref{Grid interpolation and smoothing}.
-
address@hidden --fullinterpolation
-Do interpolation irrespective of the channels in the image, see
address@hidden interpolation and smoothing}.
-
address@hidden -T INT
address@hidden --smoothwidth=INT
-The width of the average filter kernel used to smooth the interpolated
-image in units of pixels. See @ref{Grid interpolation and smoothing}. If
-this option is given a value of 1 (one), then no smoothing will be done.
-
address@hidden --fullsmoothing
-Do smoothing irrespective of the channels in the image, see @ref{Grid
-interpolation and smoothing}.
-
address@hidden --meshbasedcheck
-Store the fixed value in each mesh in one check image pixel, see
address@hidden grid values}. Note that this image has no world
-coordinate system.
-
address@hidden table
-
-
-
address@hidden Invoking astsubtractsky, , Tiling an image, SubtractSky
address@hidden Invoking SubtractSky
-
-SubtractSky will find the sky value on a grid in an image and subtract
-it. The executable name is @file{astsubtractsky} with the following
-general template:
-
address@hidden
-$ astsubtractsky [OPTION ...] Image
address@hidden example
-
address@hidden
-One line examples:
-
address@hidden
-$ astsubtractsky image.fits
-$ astsubtractsky --hdu=0 --mhdu=1 image.fits
-$ astsubtractsky --nch1=2 --nch2=2 image.fits
address@hidden example
-
address@hidden Blank pixel
-The only required input to SubtractSky is the input data file that
-currently has to be only a 2D image. But in the future it might be useful
-to use it for 1D or 3D data too. Any pixels in the image with a blank value
-will be ignored in the analysis, see @ref{Blank pixels}. The common options
-to all Gnuastro programs can be seen in @ref{Common options} and input data
-formats are explained in @ref{Arguments}.
-
-SubtractSky uses a mesh grid to tile the image. This enables it to deal
-with possible gradients, see @ref{Tiling an image}. The mesh grid options
-are common to all the programs using it, and are listed in @ref{Mesh grid
-options}. In order to ignore some pixels during the analysis, you can use
-Gnuastro's @ref{Arithmetic} program to give those pixels a blank value
-before feeding them into SubtractSky (see @ref{Blank pixels}). For a more
-accurate result, a kernel file name can be specified so the image is first
-convolved, see @ref{Quantifying signal in a mesh}, see @ref{Convolution
-kernel} for the kernel related options. The @option{--kernel} option is not
-mandatory and if not specified anywhere prior to running, the original
-image will be used. Please run SubtractSky with the @option{--help} option
-to list all the recognized options, irrespective of which part of the book
-they are fully explained in.
-
address@hidden @option
address@hidden -u FLT
address@hidden --sigclipmultip=FLT
-The multiple of the standard deviation to clip from the distribution in
address@hidden This is necessary to remove the effect of cosmic
-rays, see @ref{Sky value} and @ref{Sigma clipping}.
address@hidden -t FLT
address@hidden --sigcliptolerance=FLT
-The tolerance of sigma clipping. If the fractional change in the standard
-deviation before and after @mymath{\sigma}-clipping is less than the value
-given to this option, @mymath{\sigma}-clipping will stop, see @ref{Sigma
-clipping}.
address@hidden --checksky
-A two extension FITS image ending with @file{_smooth.fits} will be
-created showing how the interpolated sky value is smoothed.
address@hidden --checkskystd
-In the interpolation and sky checks above, include the sky standard
-deviation too. By default, only the sky value is shown in all the
-checks. However with this option, an extension will be added showing
-how the standard deviation on each mesh is finally found too.
address@hidden table
@@ -11413,6 +10883,7 @@ Statistics program is designed for such situations.
@menu
* Histogram and Cumulative Frequency Plot:: Basic definitions.
* Sigma clipping:: Definition of @mymath{\sigma}-clipping
+* Sky value::
* Invoking aststatistics:: Arguments and options to Statistics.
@end menu
@@ -11490,7 +10961,7 @@ determined from the actual dataset (none of these
options is called), the
last element in the dataset is included in the last bin's count.
address@hidden Sigma clipping, Invoking aststatistics, Histogram and
Cumulative Frequency Plot, Statistics
address@hidden Sigma clipping, Sky value, Histogram and Cumulative Frequency
Plot, Statistics
@subsection Sigma clipping
Let's assume that you have pure noise (centered on zero) with a clear
@@ -11542,8 +11013,8 @@ Go back to step 1, unless the selected exit criteria is
reached.
@noindent
The reason the median is used as a reference and not the mean is that the
mean is too significantly affected by the presence of outliers, while the
-median is less affected, see @ref{Finding the sky value}. As you can tell
-from this algorithm, besides the condition above (that the signal have
+median is less affected, see @ref{Quantifying signal in a tile}. As you can
+tell from this algorithm, besides the condition above (that the signal have
clear high signal to noise boundaries) @mymath{\sigma}-clipping is only
useful when the signal does not cover more than half of the full data
set. If they do, then the median will lie over the outliers and
@@ -11585,7 +11056,289 @@ ASCII histogram that is printed on the command-line
with
@option{--asciihist} is good enough in most cases.
@end cartouche
address@hidden Invoking aststatistics, , Sigma clipping, Statistics
+
+
+
address@hidden Sky value, Invoking aststatistics, Sigma clipping, Statistics
address@hidden Sky value
+
address@hidden Sky
+One of the most important aspects of a dataset is its reference value: the
+value of the dataset where there is no signal. Without knowing, and thus
+removing the effect of, this value it is impossible to compare the derived
+results of many high-level analyses over the dataset with other datasets
+(in the attempt to associate our results with the ``real'' world). In
+astronomy, this reference value is known as the ``Sky'' value: the value
+where there is no signal from objects (for example galaxies, stars, planets
+or comets). Depending on the dataset, the Sky value maybe a fixed value
+over the whole dataset, or it may vary based on location. For an example of
+the latter case, see Figure 11 in @url{https://arxiv.org/abs/1505.01664,
+Akhlaghi and Ichikawa (2015)}.
+
+Because of the significance of the Sky value in astronomical data analysis,
+we have devoted this subsection to it for a thorough review. We start with
+a thorough discussion on its definition (@ref{Sky value definition}). In
+the astronomical literature, researchers use a variety of methods to
+estimate the Sky value, so in @ref{Sky value misconceptions}) we review
+those and discuss their biases. From the definition of the Sky value, the
+most accurate way to estimate the Sky value is to run a detection algorithm
+(for example @ref{NoiseChisel}) over the dataset and use the un-detected
+pixels. However, there is also a more crude method that maybe useful when
+good direct detection is not initially possible (for example due to too
+many cosmic rays in a shallow image). A more crude (but simpler method)
+that is usable in such situations is discussed in @ref{Quantifying signal
+in a tile}.
+
address@hidden
+* Sky value definition:: Definition of the Sky/reference value.
+* Sky value misconceptions:: Wrong methods to estimate the Sky value.
+* Quantifying signal in a tile:: Method to estimate the presence of signal.
address@hidden menu
+
address@hidden Sky value definition, Sky value misconceptions, Sky value, Sky
value
address@hidden Sky value definition
+
address@hidden Sky value
+This analysis is taken from @url{https://arxiv.org/abs/1505.01664, Akhlaghi
+and Ichikawa (2015)}. Let's assume that all instrument defects -- bias,
+dark and flat -- have been corrected and the brightness (see @ref{Flux
+Brightness and magnitude}) of a detected object, @mymath{O}, is
+desired. The sources of flux on pixel @address@hidden this analysis
+the dimension of the data (image) is irrelevant. So if the data is an image
+(2D) with width of @mymath{w} pixels, then a pixel located on column
address@hidden and row @mymath{y} (where all counting starts from zero and (0,
+0) is located on the bottom left corner of the image), would have an index:
address@hidden of the image can be written as follows:
+
address@hidden
address@hidden
+Contribution from the target object, (@mymath{O_i}).
address@hidden
+Contribution from other detected objects, (@mymath{D_i}).
address@hidden
+Undetected objects or the fainter undetected regions of bright
+objects, (@mymath{U_i}).
address@hidden
address@hidden Cosmic rays
+A cosmic ray, (@mymath{C_i}).
address@hidden
address@hidden Background flux
+The background flux, which is defined to be the count if none of the
+others exists on that pixel, (@mymath{B_i}).
address@hidden itemize
address@hidden
+The total flux in this pixel (@mymath{T_i}) can thus be written as:
+
address@hidden
+
address@hidden Cosmic ray removal
address@hidden
+By definition, @mymath{D_i} is detected and it can be assumed that it is
+correctly estimated (deblended) and subtracted, thus @mymath{D_i=0}. There
+are also methods to detect and remove cosmic rays, for example the method
+described in van Dokkum (2001)@footnote{van Dokkum,
+P. G. (2001). Publications of the Astronomical Society of the Pacific.
+113, 1420.}, or by comparing multiple exposures. This allows us to set
address@hidden Note that in practice, @mymath{D_i} and @mymath{U_i} are
+correlated, because they both directly depend on the detection algorithm
+and its input parameters. Also note that no detection or cosmic ray removal
+algorithm is perfect. With these limitations in mind, the observed Sky
+value for this pixel (@mymath{S_i}) can be defined as
+
address@hidden Sky value
address@hidden
+
address@hidden
+Therefore, as the detection process (algorithm and input parameters)
+becomes more accurate, or @mymath{U_i\to0}, the sky value will tend to
+the background value or @mymath{S_i\to B_i}. Therefore, while
address@hidden is an inherent property of the data (pixel in an image),
address@hidden depends on the detection process. Over a group of pixels,
+for example in an image or part of an image, this equation translates
+to the average of undetected pixels. With this definition of sky, the
+object flux in the data can be calculated with
+
address@hidden T_{i}=S_{i}+O_{i} \quad\rightarrow\quad
+ O_{i}=T_{i}-S_{i}.}
+
address@hidden
address@hidden photo-electrons
+Hence, the more accurately @mymath{S_i} is measured, the more accurately
+the brightness (sum of pixel values) of the target object can be measured
+(photometry). Any under-(over-)estimation in the sky will directly
+translate to an over-(under-)estimation of the measured object's
+brightness. In the fainter outskirts of an object a very small fraction of
+the photo-electrons in the pixels actually belong to objects (see Figure 1b
+in @url{https://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa
+(2015)}). Therefore even a small over estimation of the sky value will
+result in the loss of a very large portion of most galaxies. Besides the
+lost area/brightness, this will also cause an over-estimation of the Sky
+value and thus even more under-estimation of the object's brightness. It is
+thus very important to detect the diffuse flux of a target, even if they
+are not your primary target.
+
address@hidden
address@hidden
+The @strong{Sky value} is only correctly found when all the detected
+objects (@mymath{D_i} and @mymath{C_i}) have been removed from the data.
address@hidden cartouche
+
+
+
+
address@hidden Sky value misconceptions, Quantifying signal in a tile, Sky
value definition, Sky value
address@hidden Sky value misconceptions
+
+As defined in @ref{Sky value}, the sky value is only accurately defined
+when the detection algorithm is not significantly reliant on the sky
+value. In particular its detection threshold. However, most signal-based
+detection address@hidden to Akhlaghi and Ichikawa (2015),
+signal-based detection is a detection process that relies heavily on
+assumptions about the to-be-detected objects. This method was the most
+heavily used technique prior to the introduction of NoiseChisel in that
+paper.} use the sky value as a reference to define the detection
+threshold. So these old techniques had to rely on approximations based on
+other assumptions about the data. A review of those other techniques can be
+seen in Appendix A of Akhlaghi and Ichikawa (2015)@footnote{Akhlaghi M.,
+Ichikawa. T. (2015). Astrophysical Journal Supplement Series.}. Since they
+were extensively used in astronomical data analysis for several decades,
+such approximations have given rise to a lot of misconceptions, ambiguities
+and disagreements about the sky value and how to measure it. As a summary,
+the major methods used until now were an approximation of the mode of the
+image pixel distribution and @mymath{\sigma}-clipping.
+
address@hidden
address@hidden Histogram
address@hidden Distribution mode
address@hidden Mode of a distribution
address@hidden Probability density function
address@hidden
+To find the mode of a distribution those methods would either have to
+assume (or find) a certain probability density function (PDF) or use the
+histogram. But astronomical datasets can have any distribution, making it
+almost impossible to define a generic function. Also, histogram-based
+results are very inaccurate (there is a large dispersion) and it depends on
+the histogram bin-widths.
+
address@hidden sigma-clipping
address@hidden
+Another approach was to iteratively clip the brightest pixels in the
+image (which is known as @mymath{\sigma}-clipping, since the reference
+was found from the image mean and its standard deviation or
address@hidden). See @ref{Sigma clipping} for a complete
+explanation. The problem with @mymath{\sigma}-clipping was that real
+astronomical objects have diffuse and faint wings that penetrate
+deeply into the noise. So only removing their brightest parts is
+completely useless in removing the systematic bias an object's fainter
+parts cause in the sky value.
address@hidden itemize
+
+As discussed in @ref{Sky value}, the sky value can only be correctly
+defined as the average of undetected pixels. Therefore all such
+approaches that try to approximate the sky value prior to detection
+are ultimately poor approximations.
+
+
+
address@hidden Quantifying signal in a tile, , Sky value misconceptions, Sky
value
address@hidden Quantifying signal in a tile
+
address@hidden Data
address@hidden Noise
address@hidden Signal
address@hidden Gaussian distribution
+Put simply, noise can characterized with a certain spread about a
+characteristic value. In the Gaussian distribution (most commonly used to
+model noise) the spread is defined by the standard deviation about the
+characteristic mean. Before continuing let's clarify some definitions
+first: @emph{Data} is defined as the combination of signal and noise (so a
+noisy image is one @emph{data}-set). @emph{Signal} is defined as the mean
+of the noise on each element (after sky subtraction, see @ref{Sky value
+definition}).
+
+Let's assume that the @emph{background} (see @ref{Sky value definition}) is
+subtracted and is zero. When a data set doesn't have any signal (only
+noise), the mean, median and mode of the distribution are equal within
+statistical errors and approximately equal to the background value. Signal
+always has a positive value and will never become negative, see Figure 1 in
address@hidden://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa
+(2015)}. Therefore, as more signal is added to the raw noise, the mean,
+median and mode of the dataset (which has both signal and noise) shift to
+the positive. The mean's shift is the largest. The median shifts less,
+since it is defined based on an ordered distribution and so is not affected
+by a small number of outliers. The distribution's mode shifts the least to
+the positive.
+
address@hidden Mode
address@hidden Median
+Inverting the argument above gives us a robust method to quantify the
+significance of signal in a dataset. Namely, when the mode and median of a
+distribution are approximately equal, we can argue that there is no
+significant signal. To allow for gradients (which are commonly present in
+ground-based images), we can consider the image to be made of a grid of
+tiles (see @address@hidden options to customize the
+tessellation are discussed in @ref{Processing options}.}). Hence, from the
+difference of the mode and median on each tile, we can `detect' the
+significance of signal in it. The median of a distribution is defined to be
+the value of the distribution's middle point after sorting (or 0.5
+quantile). Thus, to estimate the presence of signal, we'll compare with the
+quantile of the mode with 0.5, if the difference is larger than the value
+given to the @option{--modmedqdiff} option, this tile will be ignored. You
+can read this option as ``mode-median-quantile-diff''.
+
+This method to use the input's skewness is possible because of a new
+algorithm to find the mode of a distribution that was defined in Appendix C
+of Akhlaghi and Ichikawa (2015). However, the raw dataset's distribution is
+noisy (noise also affects the sorting), so using the argument above on the
+raw input will give a noisy result. To decrease the noise/error in
+estimating the mode, we will use convolution (see @ref{Convolution
+process}). Convolution decreases the range of the dataset and enhances its
+skewness, See Section 3.1.1 and Figure 4 in Akhlaghi and Ichikawa
+(2015). This enhanced skewness can be interpreted as an increase in the
+Signal to noise ratio of the objects buried in the noise. Therefore, to
+obtain an even better measure of the presence of signal in a mesh, the
+image can be convolved with a given kernel first.
+
+Note that through the difference of the mode and median we have actually
+`detected' data in the distribution. However this ``detection'' was only
+based on the total distribution of the data in each tile (a much lower
+resolution). This is the main limitation of this technique. The best
+approach is thus to do detection over the dataset, mask all the detected
+pixels and use the undetected regions to estimate the sky and its standard
+deviation.
+
address@hidden Cosmic rays
+The mean value of the tiles that have an approximately equal mode and
+median will be the Sky value. However there is one final hudle:
+astronomical datasets are commonly plagued with Cosmic rays. Images of
+Cosmic rays aren't smoothed by the atmosphere or telescope aperture, so
+they have sharp boundaries. Also, since they don't occupy too many pixels,
+they don't affect the mode and median calculation. But their very high
+values can greatly bias the calculation of the mean (recall how the mean
+shifts the fastest in the presence of outliers), see Figure 15 in Akhlaghi
+and Ichikawa (2015) for one example.
+
+The effect of outliers like cosmic rays on the mean and standard deviation
+can be removed through @mymath{\sigma}-clipping, see @ref{Sigma clipping}
+for a complete explanation. Therefore, after asserting that the mode and
+median are approximately equal in a tile (see @ref{Tessellation}), the
+final Sky value and its standard deviation are determined after
address@hidden with the @option{--sigmaclip} option.
+
+
+
+
+
+
+
+
+
+
+
+
+
+
address@hidden Invoking aststatistics, , Sky value, Statistics
@subsection Invoking Statistics
Statistics will print statistical measures of an input dataset (table
@@ -11603,19 +11356,21 @@ One line examples:
$ aststatistics image.fits
$ aststatistics tab.fits -c/MAG/ --histogram
$ aststatistics catalog.fits -h1 --column=MAG_F160W
+$ aststatistics image.fits --sky --kernel=kernel.fits
$ aststatistics cat.fits -cMAG -g26 -l27 --sigmaclip=3,0.2
$ aststatistics tab.txt -rPHOTO_Z -g3 -cMAG_F160W --median
@end example
@noindent
-An input image or table is necessary when processing is to be done. If a
-name is given to the @option{--output} option, it is used as the base name
-for the generated files (when applicable). Without @option{--output}, the
-input name will be used to generate an output name, see @ref{Automatic
+An input image or table is necessary when processing is to be done. If any
+output file is to be created, the value to the @option{--output} option, is
+used as the base name for the generated files. Without @option{--output},
+the input name will be used to generate an output name, see @ref{Automatic
output}. The options described below are particular to Statistics, but for
general operations, it shares a large collection of options with the other
-Gnuastro programs, see @ref{Common options}. Options can also be given in
-configuration files, for more, please see @ref{Configuration files}.
+Gnuastro programs, see @ref{Common options} for the full list. Options can
+also be given in configuration files, for more, please see
address@hidden files}.
The input dataset may have blank values (see @ref{Blank pixels}), in this
case, all blank pixels are ignored during the calculation. Initially, the
@@ -11623,56 +11378,13 @@ full dataset will be read, but it is possible to
select a specific range of
data elements to use in the analysis of each run. You can either directly
specify a minimum and maximum value for the range of data elements to use
(with @option{--greaterequal} or @option{--lessthan}), or specify the range
-using quantiles (with @option{--qrange}). If a range is specified all
+using quantiles (with @option{--qrange}). If a range is specified, all
pixels outside of it are ignored before any processing.
address@hidden ASCII plot
-When no operation is requested, Statistics will print some general basic
-properties of the input dataset on the command-line like the example below
-(ran on one of the output images of @command{make check}). It includes the
-basic information about the input dataset: filename, FITS extension (if a
-FITS file), the column name (if it was a table) or number if it the column
-doesn't have a name, along with the range and units of the data (if they
-were present in the input). Following that, the basic statistics of the
-dataset like number of points, minimum, maximum, and etc are printed. The
-outputs are finalized with an ASCII histogram to give you a general feeling
-of how the data are distributed in the given range.
-
address@hidden
-Statistics (GNU Astronomy Utilities) X.X
--------
-Input: convolve_spatial_scaled_noised.fits (hdu: 0)
-Range: from (inclusive) 9500, upto (exclusive) 11000.
-Unit: Brightness
--------
- Number of elements: 9074
- Minimum: 9622.35
- Maximum: 10999.7
- Mode: 10055.45996
- Mode quantile: 0.4001983908
- Median: 10093.7
- Mean: 10143.98257
- Standard deviation: 221.80834
--------
-Histogram:
- | **
- | ******
- | *******
- | *********
- | *************
- | **************
- | ******************
- | ********************
- | *************************** *
- | ***************************************** ***
- |* **************************************************************
- |-----------------------------------------------------------------
address@hidden example
-
The following set of options are for specifying the input/outputs of
-Statistics. Recall that besides these options that are only meaningful in
-Statistics, there are many other options that are common to all Gnuastro
-programs including Statistics, see @ref{Common options} for those.
+Statistics. There are many other input/output options that are common to
+all Gnuastro programs including Statistics, see @ref{Input output options}
+for those.
@table @option
@@ -11727,26 +11439,78 @@ it.
@end table
address@hidden ASCII plot
+When no operation is requested, Statistics will print some general basic
+properties of the input dataset on the command-line like the example below
+(ran on one of the output images of @command{make address@hidden can
+try it by running the command in the @file{tests} directory, open the image
+with a FITS viewer and have a look at it to get a sence of how these
+statistics relate to the input image/dataset.}). This default behavior is
+designed to help give you a general feeling of how the data are distributed
+and help in narrowing down your analysis.
+
address@hidden
+$ aststatistics convolve_spatial_scaled_noised.fits \
+ --greaterequal=9500 --lessthan=11000
+Statistics (GNU Astronomy Utilities) X.X
+-------
+Input: convolve_spatial_scaled_noised.fits (hdu: 0)
+Range: from (inclusive) 9500, upto (exclusive) 11000.
+Unit: Brightness
+-------
+ Number of elements: 9074
+ Minimum: 9622.35
+ Maximum: 10999.7
+ Mode: 10055.45996
+ Mode quantile: 0.4001983908
+ Median: 10093.7
+ Mean: 10143.98257
+ Standard deviation: 221.80834
+-------
+Histogram:
+ | **
+ | ******
+ | *******
+ | *********
+ | *************
+ | **************
+ | ******************
+ | ********************
+ | *************************** *
+ | ***************************************** ***
+ |* **************************************************************
+ |-----------------------------------------------------------------
address@hidden example
+
+Gnuastro's Statistics is a very general purpose program, so to be able to
+easily understand this diversity in its operations (and how to possibly run
+them together), we'll devided the operations into two types: those that
+don't respect the position of the elements and those that do (by
+tessellating the input on a tile grid, see @ref{Tessellation}). The former
+treat the whole dataset as one and can re-arrange all the elements (for
+example sort them), but the former do their processing on each tile
+independently. First, we'll review the operations that work on the whole
+dataset.
+
@cindex AWK
@cindex GNU AWK
-When you only need a single valued statistical measurement, you can use the
-group of options below. They will print only the requested value as one
-field in a line/row, like the @option{--mean}, @option{--median}
-options. These options can be called any number of times and in any
-order. The outputs of all such options will be printed on one line
-following each other (with a space character between them). This feature
-makes these options very useful in scripts, or to redirect into programs
-like GNU AWK for further immediate higher-level processing, or when you
-don't want the clutter of all the information when no option is
-specified. These are some of the most basic measures, Gnuastro is still
-under heavy development, if you want another statistical parameter, please
-contact us and we will do out best to add it to this list.
+The group of options below can used to get a single value measurement of
+the whole dataset. They will print only the requested value as one field in
+a line/row, like the @option{--mean}, @option{--median} options. These
+options can be called any number of times and in any order. The outputs of
+all such options will be printed on one line following each other (with a
+space character between them). This feature makes these options very useful
+in scripts, or to redirect into programs like GNU AWK for higher-level
+processing. These are some of the most basic measures, Gnuastro is still
+under heavy development and this list will grow. If you want another
+statistical parameter, please contact us and we will do out best to add it
+to this list, see @ref{Suggest new feature}.
@table @option
@item -n
@itemx --number
-Print the number of all used elements.
+Print the number of all used (non-blank and in range) elements.
@item --minimum
Print the minimum value of all used elements.
@@ -11838,8 +11602,8 @@ more.
The list of options below are for those statistical operations that output
more than one value. So while they can be called together in one run, their
-outputs will be distinct (each one's output will usually take in more than
-one line).
+outputs will be distinct (each one's output will usually be printed in more
+than one line).
@table @option
@@ -11934,19 +11698,13 @@ and nice histogram than what the raw command-line can
offer you (with the
Save the cumulative frequency plot of the usable values in the input
dataset into a table, similar to @option{--histogram}.
address@hidden -s FLT,FLT
address@hidden --sigmaclip=FLT,FLT
address@hidden -s
address@hidden --sigmaclip
Do @mymath{\sigma}-clipping on the usable pixels of the input dataset. See
@ref{Sigma clipping} for a full description on @mymath{\sigma}-clipping and
-also to better understand this option.
-
-This option takes two values which are separated by a comma (@key{,}). Each
-value can either be written as a single number or as a fraction of two
-numbers (for example @code{3,1/10}). The first value to this option is the
-multiple of @mymath{\sigma} that will be clipped (@mymath{\alpha} in that
-section). The second value is the exit criteria. If it is less than 1, then
-it is interpretted as tolerance and if it is larger than one it is a
-specific number. Hence, in the latter case the value must be an integer.
+also to better understand this option. The @mymath{\sigma}-clipping
+parameters can be set through the @option{--sclipparams} option (see
+below).
@item --mirror=FLT
Make a histogram and cumulative frequency plot of the mirror distribution
@@ -11975,8 +11733,8 @@ it is possible to make plots like Figure 21 of
@end table
-The final set of options in Statistics describe the histogram and
-cumulative frequency plot settings (for the @option{--histogram},
+The list of options below allow customization of the histogram and
+cumulative frequency plots (for the @option{--histogram},
@option{--cumulative}, @option{--asciihist}, and @option{--asciicfp}
options).
@@ -12010,14 +11768,14 @@ and the rest are similarly scaled. In some situations
(for example if you
want to plot the histogram and cumulative frequency plot in one plot) this
can be very useful.
address@hidden --onebinvalue=FLT
address@hidden --onebinstart=FLT
Make sure that one bin starts with the value to this option. In practice,
this will shift the bins used to find the histogram and cumulative
frequency plot such that one bin's lower interval becomes this value. For
example when the histogram range includes negative and positive values and
zero has a special significance in your analysis, then zero will be
somewhere in one bin and will mix counts of positive and negative. By
-setting @option{--onebinvalue=0}, you can make sure that the viewers of the
+setting @option{--onebinstart=0}, you can make sure that the viewers of the
histogram will not be confused without doing the math of setting a range
and number of bins.
@@ -12028,11 +11786,100 @@ you will not see the 0.000 in the first column, you
will see two symmetric
values. If the value is not within the usable input range, this option will
be ignored.
address@hidden table
+
+
+All the options described until now were from the first class of operations
+discusssed above: those that treat the whole dataset as one. However. It
+often happens that the relative position of the dataset elements over the
+dataset is significant. For example you don't want one median value for the
+whole input image, you want to know how the median changes over the
+image. For such operations, the input has to be tessellated (see
address@hidden). Thus this class of options can't currently be called
+along with the options above in one run of Statistics.
+
address@hidden @option
+
address@hidden -t
address@hidden --ontile
+Do the respective single-valued calculation over one tile of the input
+dataset, not the whole dataset. This option must be called with atleast one
+of the single valued options discussed above (for example @option{--mean}
+or @option{--quantile}). The output will be a file in the same format as
+the input. If the @option{--oneelempertile} option is called, then one
+element/pixel will be used for each tile (see @ref{Processing
+options}). Otherwise, the output will have the same size as the input, but
+each element will have the value corresponding to that tile's value. If
+multiple single valued operations are called, then for each operation there
+will be one extension in the output FITS file.
+
address@hidden -y
address@hidden --sky
+Estimate the Sky value on each tile as fully described in @ref{Quantifying
+signal in a tile}. As described in that section, several options are
+necessary to configure the Sky estimation which are listed below. The
+output file will have two extensions: the first is the Sky value and the
+second is the Sky standard deviation on each tile. Similar to
address@hidden, if the @option{--oneelempertile} option is called, then
+one element/pixel will be used for each tile (see @ref{Processing
+options}).
@end table
+The parameters for estimating the sky value can be set with the following
+options, except for the @option{--sclipparams} option (which is also used
+by the @option{--sigmaclip}), the rest are only used for the Sky value
+estimation.
+
address@hidden @option
+
address@hidden -k=STR
address@hidden --kernel=STR
+File name of kernel to help in estimating the significance of signal in a
+tile, see @ref{Quantifying signal in a tile}.
+
address@hidden --khdu=STR
+Kernel HDU to help in estimating the significance of signal in a tile, see
address@hidden signal in a tile}.
+
address@hidden --mirrordist=FLT
+Maximum distance (as a multiple of error) to estimate the difference
+between the input and mirror distributions in finding the mode, see
+Appendix C of @url{https://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa
+2015}, also see @ref{Quantifying signal in a tile}.
+
address@hidden --modmedqdiff=FLT
+The maximum acceptable distance between the mode and median, see
address@hidden signal in a tile}.
+
address@hidden --sclipparams=FLT,FLT
+The @mymath{\sigma}-clipping parameters. This option takes two values which
+are separated by a comma (@key{,}). Each value can either be written as a
+single number or as a fraction of two numbers (for example
address@hidden,1/10}). The first value to this option is the multiple of
address@hidden that will be clipped (@mymath{\alpha} in that section). The
+second value is the exit criteria. If it is less than 1, then it is
+interpretted as tolerance and if it is larger than one it is a specific
+number. Hence, in the latter case the value must be an integer.
+
address@hidden --smoothwidth=INT
+Width of a flat kernel to convolve the interpolated tile values. Tile
+interpolation is done using the median of the @option{--interpnumngb}
+neighbors of each tile (see @ref{Processing options}). If this option is
+given a value of zero or one, no smoothing will be done. Without smoothing,
+strong boundaries will probably be created between the values estimated for
+each tile. It is thus good to smooth the interpolated image so strong
+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.
address@hidden --checksky
+Create a multi-extension FITS file showing the steps that were used to
+estimate the Sky value over the input, see @ref{Quantifying signal in a
+tile}. The file will have two extensions for each step (one for the Sky and
+one for the Sky standard deviation).
address@hidden table
@node NoiseChisel, MakeCatalog, Statistics, Data analysis
@section NoiseChisel
@@ -12119,27 +11966,26 @@ Akhlaghi and Ichikawa (2015) to learn why this
particular kernel was
chosen as default. See @ref{Convolution kernel} for kernel related
options.
-NoiseChisel uses a mesh grid to tile the image. This enables it to
-deal with possible gradients, see @ref{Tiling an image}. The mesh grid
-options are common to all the programs using it, and are listed in
address@hidden grid options}. In particular, NoiseChisel uses two mesh
-grids: a large and a small one. The sizes of the meshs in each grid
-can be specified with the following two options (the
address@hidden option is not recognized by NoiseChisel). The rest
-of the options explained in @ref{Mesh grid options} apply to both
-grids.
+NoiseChisel uses a mesh grid to tile the image. This enables it to deal
+with possible gradients, see @ref{Tessellation}. The mesh grid options are
+common to all the programs using it, and are listed in @ref{Processing
+options}. In particular, NoiseChisel uses two mesh grids: a large and a
+small one. The sizes of the meshs in each grid can be specified with the
+following two options (the @option{--meshsize} option is not recognized by
+NoiseChisel). The rest of the options explained in @ref{Processing options}
+apply to both grids.
@table @option
@item -s
@itemx --smeshsize
-Similar to @option{--meshsize} in @ref{Mesh grid options}, but for the
+Similar to @option{--meshsize} in @ref{Processing options}, but for the
smaller mesh grid, which is most dependent on the gradients in the
image.
@item -l
@itemx --lmeshsize
-Similar to @option{--meshsize} in @ref{Mesh grid options}, but for the
+Similar to @option{--meshsize} in @ref{Processing options}, but for the
larger mesh grid, used for detection and segmentation Signal to noise
ratio analysis.
@@ -12217,7 +12063,7 @@ many detected regions should be excluded.
The minimum number of `psudo-detections' (in identifying false detections)
or clumps (in identifying false clumps) in each large mesh grid. If their
number is less than this value, this mesh will be left blank and filled
-during mesh interpolation, see @ref{Grid interpolation and smoothing}.
+during mesh interpolation.
The Signal to noise ratio of false detections and clumps in each mesh
is found using the quantile of the Signal to noise ratio distribution
@@ -12240,9 +12086,9 @@ can customize the detection process in NoiseChisel.
@itemx --qthresh=FLT
The quantile threshold to apply to the convolved image. The detection
process begins with applying a quantile threshold to each of the small mesh
-grid elements, see @ref{Tiling an image}. The quantile is only calculated
-for those meshs that don't have any significant signal within them, see
address@hidden signal in a mesh}.
+grid elements, see @ref{Tessellation}. The quantile is only calculated for
+those meshs that don't have any significant signal within them, see
address@hidden signal in a tile}.
@cindex Quantile
@cindex Binary image
@@ -12260,9 +12106,8 @@ are known as background (have a value of 0).
@cindex NaN
@item --checkthreshold
-Check the quantile threshold values on the mesh grid. A file suffixed
-with @file{_thresh.fits} will be created, see @ref{Checking grid
-values}.
+Check the quantile threshold values on the mesh grid. A file suffixed with
address@hidden will be created, see @ref{Tessellation}.
@item -e INT
@itemx --erode=INT
@@ -12353,13 +12198,12 @@ smaller than the value given to this option. Note that
smaller, in which case they can be considered to be approximately equal.
@item --checkdetectionsky
-Check the initial approximation of the sky value and its standard
-deviation in a FITS file ending with @file{_detsky.fits}. See
address@hidden grid values} for more information. If
address@hidden is not called, then the first extension will
-be the the binary image with initial detections labeled one and
-background labeled zero. The mesh values will be in the subsequent
-extensions.
+Check the initial approximation of the sky value and its standard deviation
+in a FITS file ending with @file{_detsky.fits}. See @ref{Tessellation} for
+more information. If @option{--meshbasedcheck} is not called, then the
+first extension will be the the binary image with initial detections
+labeled one and background labeled zero. The mesh values will be in the
+subsequent extensions.
@item -R FLT
@itemx --dthresh=FLT
@@ -12636,8 +12480,7 @@ detected region, then the whole detected region is
given a label of 1.
@item
The final sky value on each pixel. See @ref{Sky value} for a complete
-explanation. See @ref{Grid interpolation and smoothing} for an
-explanation on the boxy appearance of this image.
+explanation.
@item
Similar to the previous mesh but for the standard deviation on each
@@ -12773,10 +12616,10 @@ brightness) can be seen more clearly.
@cindex Depth
The outputs of NoiseChisel include the Sky standard deviation
(@mymath{\sigma}) on every group of pixels (a mesh) that were calculated
-from the undetected pixels in that mesh, see @ref{Tiling an image} and
+from the undetected pixels in that mesh, see @ref{Tessellation} and
@ref{NoiseChisel output}. Let's take @mymath{\sigma_m} as the median
@mymath{\sigma} over the successful meshes in the image (prior to
-interpolation or smoothing, see @ref{Grid interpolation and smoothing}).
+interpolation or smoothing).
On different instruments pixels have different physical sizes (for example
in micro-meters, or spatial angle over the sky), nevertheless, a pixel is
@@ -17282,11 +17125,11 @@ within the nodes will also be discarded.
@subsection Mesh grid for an image (@file{mesh.h})
The basic concepts behind the mesh and channel grids in Gnuastro were
-explained in @ref{Tiling an image}. Besides the improved measurement
-purposes introduced there, tiling an image can also greatly speed up some
-processing operations over the image, since each mesh can be given to a
-different thread and so the CPU can be working on multiple regions of the
-image simultaneously. The mesh structure and functions introduced here are
+explained in @ref{Tessellation}. Besides the improved measurement purposes
+introduced there, tiling an image can also greatly speed up some processing
+operations over the image, since each mesh can be given to a different
+thread and so the CPU can be working on multiple regions of the image
+simultaneously. The mesh structure and functions introduced here are
declared in @file{gnuastro/mesh.h} and allow you to also use operate on an
image mesh.
@@ -17320,16 +17163,15 @@ functions is the @code{garray}, where each mesh has a
unique
@code{gid}. @code{garray} stands for grid-array. It is used for operations
on the mesh grid, where one value is to be assigned for each mesh. It has
one element for each mesh in the image. However, due to the (possible)
-existence of channels (see @ref{Tiling an image}), each channel needs its
-own contiguous part (group of meshes) in the full @code{garray}. Each
-channel has @code{gs0*gs1} (dimensions of meshes in each channel)
-elements. There are @code{nch} parts (or channels) in total. In short, the
-meshs in each channel have to be contiguous to facilitate the neighbor
-analysis in interpolation and other channel specific jobs. So, the over-all
-structure can be viewed as a 3D array. The operations on the meshs might
-need more than one output, for example the mean and the standard
-deviation. So we have two @code{garray}s and two @code{nearest} arrays (for
-interpolation).
+existence of channels (see @ref{Tessellation}), each channel needs its own
+contiguous part (group of meshes) in the full @code{garray}. Each channel
+has @code{gs0*gs1} (dimensions of meshes in each channel) elements. There
+are @code{nch} parts (or channels) in total. In short, the meshs in each
+channel have to be contiguous to facilitate the neighbor analysis in
+interpolation and other channel specific jobs. So, the over-all structure
+can be viewed as a 3D array. The operations on the meshs might need more
+than one output, for example the mean and the standard deviation. So we
+have two @code{garray}s and two @code{nearest} arrays (for interpolation).
@end deftp
@deftypefun size_t gal_mesh_ch_based_id_from_gid (struct gal_mesh_params
@code{*mp}, size_t @code{gid})
@@ -19700,10 +19542,6 @@ shaped signal in noise (galaxies in the sky), using
thresholds that are
below the Sky value, see @url{http://arxiv.org/abs/1505.01664,
arXiv:1505.01664}.
address@hidden SubtractSky
-(@file{astsubtractsky}, @ref{SubtractSky}) Find and subtract sky value by
-comparing the mode and median on a mesh grid.
-
@item Table
(@file{asttable}, @ref{Table}) Convert FITS binary and ASCII tables into
other such tables, print them on the command-line, save them in a plain
diff --git a/lib/Makefile.am b/lib/Makefile.am
index ffa5cb4..7fd8ef9 100644
--- a/lib/Makefile.am
+++ b/lib/Makefile.am
@@ -22,12 +22,15 @@
-## Necessary flags. NOTE: $(top_srcdir)/bootstrapped/lib is only necessary
-## for internally compiled utilities and libraries. It must not be included
-## during the tests since the bootstrapped libraries are not installed.
+## Necessary flags.
##
-## The SYSCONFIG_DIR is necessary in `options.c'
-AM_CPPFLAGS = -I\$(top_srcdir)/bootstrapped/lib \
+## $(top_srcdir)/bootstrapped/lib: only necessary for the libraries since
+## the Gnulib functions will be statically linked to the Gnuastro
+## library so linking to Gnuastro is enough to access them also.
+##
+## SYSCONFIG_DIR: only necessary in `options.c' to get the system
+## installation directory.
+AM_CPPFLAGS = -I\$(top_srcdir)/bootstrapped/lib \
-DSYSCONFIG_DIR=\"$(sysconfdir)\"
@@ -89,13 +92,6 @@ EXTRA_DIST = $(headersdir)/README gnuastro.pc.in
arithmetic-binary.h \
-# Common options to all Gnuastro programs.
-dist_sysconf_DATA = gnuastro.conf
-
-
-
-
-
# Definitions for Gnuastro's the pkg-config file (inspired from GSL's
# Makefile.am)
pkgconfig_DATA = gnuastro.pc
diff --git a/lib/commonopts.h b/lib/commonopts.h
index 48eb245..4edaca9 100644
--- a/lib/commonopts.h
+++ b/lib/commonopts.h
@@ -86,7 +86,7 @@ struct argp_option gal_commonopts_options[] =
- /* Tile grid options. */
+ /* Tile grid (tessellation) options. */
{
0, 0, 0, 0,
"Tessellation (tile grid):",
@@ -172,6 +172,32 @@ struct argp_option gal_commonopts_options[] =
GAL_OPTIONS_NOT_MANDATORY,
GAL_OPTIONS_NOT_SET
},
+ {
+ "interponlyblank",
+ GAL_OPTIONS_KEY_INTERPONLYBLANK,
+ 0,
+ 0,
+ "Only interpolate over the blank tiles.",
+ GAL_OPTIONS_GROUP_TESSELLATION,
+ &cp->interponlyblank,
+ GAL_OPTIONS_NO_ARG_TYPE,
+ GAL_OPTIONS_RANGE_0_OR_1,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+ {
+ "interpnumngb",
+ GAL_OPTIONS_KEY_INTERPNUMNGB,
+ "INT",
+ 0,
+ "Number of neighbors to use for interpolation.",
+ GAL_OPTIONS_GROUP_TESSELLATION,
+ &cp->interpnumngb,
+ GAL_DATA_TYPE_SIZE_T,
+ GAL_OPTIONS_RANGE_GT_0,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
diff --git a/lib/convolve.c b/lib/convolve.c
index 880704a..3b59e16 100644
--- a/lib/convolve.c
+++ b/lib/convolve.c
@@ -33,7 +33,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gnuastro/convolve.h>
#include <gnuastro/dimension.h>
-
+#include <checkset.h>
@@ -419,12 +419,16 @@ convolve_spatial_on_thread(void *inparam)
-/* Convolve a dataset with a given kernel in the spatial domain. Since
- spatial convolution can be very series of tiles arranged as an array. */
+/* Convolve a dataset with a given kernel in the spatial domain. Spatial
+ convolution can be greatly sped up if it is done on separate tiles over
+ the image (on multiple threads). So as input, you can either give tile
+ values or one full array. Just note that if you give a single array as
+ input, the `next' element has to be `NULL'.*/
gal_data_t *
gal_convolve_spatial(gal_data_t *tiles, gal_data_t *kernel,
size_t numthreads, int edgecorrection, int convoverch)
{
+ char *name;
struct spatial_params params;
gal_data_t *out, *block=gal_tile_block(tiles);
@@ -438,8 +442,10 @@ gal_convolve_spatial(gal_data_t *tiles, gal_data_t *kernel,
"`float32' type input and kernel");
/* Allocate the output convolved dataset. */
+ name = ( block->name
+ ? gal_checkset_malloc_cat("CONVL_", block->name) : NULL );
out=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT32, block->ndim, block->dsize,
- block->wcs, 0, block->minmapsize, "CONVOLVED",
+ block->wcs, 0, block->minmapsize, name,
block->unit, NULL);
/* Set the pointers in the parameters structure. */
@@ -454,5 +460,6 @@ gal_convolve_spatial(gal_data_t *tiles, gal_data_t *kernel,
gal_data_num_in_ll(tiles), numthreads);
/* Clean up and return the output array. */
+ free(name);
return out;
}
diff --git a/lib/data.c b/lib/data.c
index e23ed40..28ab527 100644
--- a/lib/data.c
+++ b/lib/data.c
@@ -1142,6 +1142,29 @@ gal_data_free_ll(gal_data_t *list)
/*************************************************************
************** Copying ***************
*************************************************************/
+/* Wrapper for `gal_data_copy_to_new_type', but will copy to the same type
+ as the input. Recall that if the input is a tile (a part of the input,
+ which is not-contiguous if it has more than one dimension), then the
+ output will have only the elements that cover the tile.*/
+gal_data_t *
+gal_data_copy(gal_data_t *in)
+{
+ return gal_data_copy_to_new_type(in, gal_tile_block(in)->type);
+}
+
+
+
+
+/* Copy the input dataset into the already allocated space of out. */
+void
+gal_data_copy_to_allocated(gal_data_t *in, gal_data_t *out)
+{
+ gal_data_copy_to_new_type_to_allocated(in, out, out->type);
+}
+
+
+
+
/* Only to be used in `data_copy_from_string'. */
static void
@@ -1334,7 +1357,7 @@ data_copy_to_string(gal_data_t *from, gal_data_t *to)
OT ob, *o=out->array; \
size_t increment=0, num_increment=1; \
IT ib, *ist, *i=in->array, *f=i+in->size; \
- size_t mclen, contig_len=in->dsize[in->ndim-1]; \
+ size_t mclen=0, contig_len=in->dsize[in->ndim-1]; \
size_t s_e_ind[2]={0,iblock->size-1}; /* -1: this is INCLUSIVE */ \
\
/* If we are on a tile, the default values need to change. */ \
@@ -1461,7 +1484,7 @@ data_copy_to_string(gal_data_t *from, gal_data_t *to)
gal_data_t *
gal_data_copy_to_new_type(gal_data_t *in, uint8_t newtype)
{
- gal_data_t *out, *iblock=gal_tile_block(in);
+ gal_data_t *out;
/* Allocate space for the output type */
out=gal_data_alloc(NULL, newtype, in->ndim, in->dsize, in->wcs,
@@ -1480,6 +1503,36 @@ gal_data_copy_to_new_type(gal_data_t *in, uint8_t
newtype)
*/
/* Fill in the output array: */
+ gal_data_copy_to_new_type_to_allocated(in, out, newtype);
+
+ /* Return the created array */
+ return out;
+}
+
+
+
+
+
+/* Copy an array into the already allocated space in `out'. */
+void
+gal_data_copy_to_new_type_to_allocated(gal_data_t *in, gal_data_t *out,
+ uint8_t newtype)
+{
+ gal_data_t *iblock=gal_tile_block(in);
+
+ /* Two small sanity checks to avoid segmentation faults. */
+ if(out->size<in->size)
+ error(EXIT_FAILURE, 0, "the output dataset must be equal or larger than "
+ "the input in `gal_data_copy_to_new_type_allocated', the sizes "
+ "are %zu and %zu respectively", out->size, in->size);
+ if(out->type!=newtype)
+ error(EXIT_FAILURE, 0, "the output dataset must have the same type as "
+ "the requested type in `gal_data_copy_to_new_type_allocated', "
+ "the types are %s and %s respectively",
+ gal_data_type_as_string(out->type, 1),
+ gal_data_type_as_string(newtype, 1));
+
+ /* Do the copying. */
switch(out->type)
{
case GAL_DATA_TYPE_UINT8: COPY_OT_SET( uint8_t ); break;
@@ -1492,7 +1545,7 @@ gal_data_copy_to_new_type(gal_data_t *in, uint8_t newtype)
case GAL_DATA_TYPE_INT64: COPY_OT_SET( int64_t ); break;
case GAL_DATA_TYPE_FLOAT32: COPY_OT_SET( float ); break;
case GAL_DATA_TYPE_FLOAT64: COPY_OT_SET( double ); break;
- case GAL_DATA_TYPE_STRING: data_copy_to_string(in, out); break;
+ case GAL_DATA_TYPE_STRING: data_copy_to_string(in, out); break;
case GAL_DATA_TYPE_BIT:
case GAL_DATA_TYPE_STRLL:
@@ -1507,9 +1560,6 @@ gal_data_copy_to_new_type(gal_data_t *in, uint8_t newtype)
error(EXIT_FAILURE, 0, "type %d not recognized for "
"for `out->type' in gal_data_copy_to_new_type", out->type);
}
-
- /* Return the created array */
- return out;
}
@@ -1549,20 +1599,6 @@ gal_data_copy_to_new_type_free(gal_data_t *in, uint8_t
type)
-/* Wrapper for `gal_data_copy_to_new_type', but will copy to the same type
- as the input. Recall that if the input is a tile (a part of the input,
- which is not-contiguous if it has more than one dimension), then the
- output will have only the elements that cover the tile.*/
-gal_data_t *
-gal_data_copy(gal_data_t *in)
-{
- return gal_data_copy_to_new_type(in, gal_tile_block(in)->type);
-}
-
-
-
-
-
/* Copy/read the element at `index' of the array in `data' into the space
pointed to by `ptr'. */
#define COPY_ELEM(IT) { IT *o=ptr, *a=input->array; *o = a[index]; }
diff --git a/lib/fits.c b/lib/fits.c
index b1ccb59..71469c3 100644
--- a/lib/fits.c
+++ b/lib/fits.c
@@ -36,6 +36,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gnuastro/git.h>
#include <gnuastro/wcs.h>
#include <gnuastro/fits.h>
+#include <gnuastro/tile.h>
#include <gnuastro/blank.h>
#include "checkset.h"
@@ -1387,30 +1388,34 @@ gal_fits_img_read_kernel(char *filename, char *hdu,
size_t minmapsize)
WCS information) into a FITS file, but will not close it. Instead it
will pass along the FITS pointer for further modification. */
fitsfile *
-gal_fits_img_write_to_ptr(gal_data_t *data, char *filename)
+gal_fits_img_write_to_ptr(gal_data_t *input, char *filename)
{
- size_t i;
void *blank;
char *wcsstr;
fitsfile *fptr;
long fpixel=1, *naxes;
- int nkeyrec, status=0, datatype=gal_fits_type_to_datatype(data->type);
+ size_t i, ndim=input->ndim;
+ gal_data_t *towrite, *block=gal_tile_block(input);
+ int nkeyrec, status=0, datatype=gal_fits_type_to_datatype(block->type);
+
+ /* If the input is a tile (isn't a contiguous region of memory), then
+ copy it into a contiguous region. */
+ towrite = input==block ? input : gal_data_copy(input);
/* Allocate the naxis area. */
naxes=gal_data_malloc_array( ( sizeof(long)==8
? GAL_DATA_TYPE_INT64
- : GAL_DATA_TYPE_INT32 ), data->ndim);
+ : GAL_DATA_TYPE_INT32 ), ndim);
/* Open the file for writing */
fptr=gal_fits_open_to_write(filename);
/* Fill the `naxes' array (in opposite order, and `long' type): */
- for(i=0;i<data->ndim;++i)
- naxes[data->ndim-1-i]=data->dsize[i];
+ for(i=0;i<ndim;++i) naxes[ndim-1-i]=towrite->dsize[i];
/* Create the FITS file. */
- fits_create_img(fptr, gal_fits_type_to_bitpix(data->type),
- data->ndim, naxes, &status);
+ fits_create_img(fptr, gal_fits_type_to_bitpix(towrite->type),
+ ndim, naxes, &status);
gal_fits_io_error(status, NULL);
/* Remove the two comment lines put by CFITSIO. Note that in some cases,
@@ -1422,12 +1427,13 @@ gal_fits_img_write_to_ptr(gal_data_t *data, char
*filename)
status=0;
/* Write the image into the file. */
- fits_write_img(fptr, datatype, fpixel, data->size, data->array, &status);
+ fits_write_img(fptr, datatype, fpixel, towrite->size, towrite->array,
+ &status);
/* If we have blank pixels, we need to define a BLANK keyword when we are
dealing with integer types. */
- if(gal_blank_present(data))
- switch(data->type)
+ if(gal_blank_present(towrite))
+ switch(towrite->type)
{
case GAL_DATA_TYPE_FLOAT32:
case GAL_DATA_TYPE_FLOAT64:
@@ -1436,7 +1442,7 @@ gal_fits_img_write_to_ptr(gal_data_t *data, char
*filename)
break;
default:
- blank=gal_blank_alloc_write(data->type);
+ blank=gal_blank_alloc_write(towrite->type);
if(fits_write_key(fptr, datatype, "BLANK", blank,
"Pixels with no data.", &status) )
gal_fits_io_error(status, "adding the BLANK keyword");
@@ -1444,25 +1450,25 @@ gal_fits_img_write_to_ptr(gal_data_t *data, char
*filename)
}
/* Write the extension name to the header. */
- if(data->name)
- fits_write_key(fptr, TSTRING, "EXTNAME", data->name, "", &status);
+ if(towrite->name)
+ fits_write_key(fptr, TSTRING, "EXTNAME", towrite->name, "", &status);
/* Write the units to the header. */
- if(data->unit)
- fits_write_key(fptr, TSTRING, "BUNIT", data->unit, "", &status);
+ if(towrite->unit)
+ fits_write_key(fptr, TSTRING, "BUNIT", towrite->unit, "", &status);
/* Write comments if they exist. */
- if(data->comment)
- fits_write_comment(fptr, data->comment, &status);
+ if(towrite->comment)
+ fits_write_comment(fptr, towrite->comment, &status);
/* If a WCS structure is present, write it in */
- if(data->wcs)
+ if(towrite->wcs)
{
/* Decompose the `PCi_j' matrix and `CDELTi' vector. */
- gal_wcs_decompose_pc_cdelt(data->wcs);
+ gal_wcs_decompose_pc_cdelt(towrite->wcs);
/* Convert the WCS information to text. */
- status=wcshdo(WCSHDO_safe, data->wcs, &nkeyrec, &wcsstr);
+ status=wcshdo(WCSHDO_safe, towrite->wcs, &nkeyrec, &wcsstr);
if(status)
error(EXIT_FAILURE, 0, "wcshdo ERROR %d: %s", status,
wcs_errmsg[status]);
@@ -1472,6 +1478,7 @@ gal_fits_img_write_to_ptr(gal_data_t *data, char
*filename)
/* Report any errors if we had any */
free(naxes);
gal_fits_io_error(status, NULL);
+ if(towrite!=input) gal_data_free(towrite);
return fptr;
}
diff --git a/lib/gnuastro/data.h b/lib/gnuastro/data.h
index 4663c49..4381d74 100644
--- a/lib/gnuastro/data.h
+++ b/lib/gnuastro/data.h
@@ -344,13 +344,20 @@ gal_data_free_ll(gal_data_t *list);
************** Copying ***************
*************************************************************/
gal_data_t *
-gal_data_copy_to_new_type(gal_data_t *in, uint8_t newtype);
+gal_data_copy(gal_data_t *in);
+
+void
+gal_data_copy_to_allocated(gal_data_t *in, gal_data_t *out);
gal_data_t *
-gal_data_copy_to_new_type_free(gal_data_t *in, uint8_t type);
+gal_data_copy_to_new_type(gal_data_t *in, uint8_t newtype);
+
+void
+gal_data_copy_to_new_type_to_allocated(gal_data_t *in, gal_data_t *out,
+ uint8_t newtype);
gal_data_t *
-gal_data_copy(gal_data_t *in);
+gal_data_copy_to_new_type_free(gal_data_t *in, uint8_t type);
void
gal_data_copy_element_same_type(gal_data_t *input, size_t index, void *ptr);
diff --git a/lib/gnuastro/interpolate.h b/lib/gnuastro/interpolate.h
index 30aab88..2c878fb 100644
--- a/lib/gnuastro/interpolate.h
+++ b/lib/gnuastro/interpolate.h
@@ -51,7 +51,7 @@ gal_data_t *
gal_interpolate_close_neighbors(gal_data_t *input,
struct gal_tile_two_layer_params *tl,
size_t numneighbors, size_t numthreads,
- int onlyblank);
+ int onlyblank, int aslinkedlist);
diff --git a/lib/gnuastro/linkedlist.h b/lib/gnuastro/linkedlist.h
index 296c175..a67adff 100644
--- a/lib/gnuastro/linkedlist.h
+++ b/lib/gnuastro/linkedlist.h
@@ -233,6 +233,31 @@ void
gal_linkedlist_free_ill(struct gal_linkedlist_ill *list);
+
+
+/******************* void * */
+struct gal_linkedlist_vll
+{
+ void* v;
+ struct gal_linkedlist_vll *next;
+};
+
+void
+gal_linkedlist_add_to_vll(struct gal_linkedlist_vll **list, void *value);
+
+void
+gal_linkedlist_pop_from_vll(struct gal_linkedlist_vll **list, void **value);
+
+void
+gal_linkedlist_reverse_vll(struct gal_linkedlist_vll **list);
+
+void
+gal_linkedlist_free_vll(struct gal_linkedlist_vll *list, int freevalue);
+
+
+
+
+
/******************* Ordered size_t: */
struct gal_linkedlist_osll
{
diff --git a/lib/gnuastro/statistics.h b/lib/gnuastro/statistics.h
index cecab44..0701135 100644
--- a/lib/gnuastro/statistics.h
+++ b/lib/gnuastro/statistics.h
@@ -181,7 +181,7 @@ gal_statistics_cfp(gal_data_t *data, gal_data_t *bins, int
normalize);
****************************************************************/
gal_data_t *
gal_statistics_sigma_clip(gal_data_t *input, float multip, float param,
- int quiet);
+ int inplace, int quiet);
diff --git a/lib/gnuastro/tile.h b/lib/gnuastro/tile.h
index 766e08a..505fb6a 100644
--- a/lib/gnuastro/tile.h
+++ b/lib/gnuastro/tile.h
@@ -80,7 +80,8 @@ gal_tile_block_write_const_value(gal_data_t *tilevalues,
gal_data_t *tilesll,
gal_data_t *
gal_tile_block_check_tiles(gal_data_t *tiles);
-
+void *
+gal_tile_block_relative_to_other(gal_data_t *tile, gal_data_t *other);
@@ -127,9 +128,6 @@ gal_tile_full_two_layers(gal_data_t *input,
struct gal_tile_two_layer_params *tl);
void
-gal_tile_full_free_contents(struct gal_tile_two_layer_params *tl);
-
-void
gal_tile_full_permutation(struct gal_tile_two_layer_params *tl);
void
@@ -137,6 +135,13 @@ gal_tile_full_values_write(gal_data_t *tilevalues,
struct gal_tile_two_layer_params *tl,
char *filename, char *program_string);
+gal_data_t *
+gal_tile_full_values_smooth(gal_data_t *tilevalues,
+ struct gal_tile_two_layer_params *tl,
+ size_t width, size_t numthreads);
+
+void
+gal_tile_full_free_contents(struct gal_tile_two_layer_params *tl);
@@ -148,7 +153,7 @@ gal_tile_full_values_write(gal_data_t *tilevalues,
#define GAL_TILE_PO_OISET(IT, OT, OP, IN, OUT, PARSE_OUT, CHECK_BLANK) { \
gal_data_t *out_w=OUT; /* Since `OUT' may be NULL. */ \
OT *ost, *o = out_w ? out_w->array : NULL; \
- IT b, *st, *i=IN->array, *f=i+IN->size; \
+ IT b=0, *st=NULL, *i=IN->array, *f=i+IN->size; \
\
/* Write the blank value for the input type into `b'). Note that */ \
/* a tile doesn't necessarily have to have a type. */ \
diff --git a/lib/interpolate.c b/lib/interpolate.c
index 5ae15da..b963173 100644
--- a/lib/interpolate.c
+++ b/lib/interpolate.c
@@ -37,7 +37,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gnuastro/interpolate.h>
#include <gnuastro/permutation.h>
-
+#include <checkset.h>
@@ -58,12 +58,13 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
struct interpolate_params
{
gal_data_t *input;
+ size_t num;
gal_data_t *out;
gal_data_t *blanks;
size_t numneighbors;
uint8_t *thread_flags;
- void *ngb_vals;
int onlyblank;
+ struct gal_linkedlist_vll *ngb_vals;
struct gal_tile_two_layer_params *tl;
};
@@ -79,25 +80,24 @@ interpolate_close_neighbors_on_thread(void *in_prm)
struct gal_threads_params *tprm=(struct gal_threads_params *)in_prm;
struct interpolate_params *prm=(struct interpolate_params *)(tprm->params);
struct gal_tile_two_layer_params *tl=prm->tl;
- int correct_index=(tl && tl->totchannels>1 && tl->workoverch==0);
+ int correct_index=(tl && tl->totchannels>1 && !tl->workoverch);
gal_data_t *input=prm->input;
/* Rest of variables. */
+ void *nv;
float pdist;
- gal_data_t *nearest, *single;
- size_t ngb_counter, dist, *dinc;
+ uint8_t *b, *bf, *bb;
+ struct gal_linkedlist_vll *tvll;
struct gal_linkedlist_tosll *lQ, *sQ;
- uint8_t *b, *bf, *bb, type=input->type;
- void *values=input->array, *out=prm->out->array;
- size_t pind, twidth=gal_data_sizeof(input->type);
- size_t i, index, fullind, chstart, ndim=input->ndim;
+ size_t ngb_counter, dist, pind, *dinc;
+ size_t i, index, fullind, chstart=0, ndim=input->ndim;
+ gal_data_t *median, *tin, *tout, *tnear, *nearest=NULL;
size_t size = (correct_index ? tl->tottilesinch : input->size);
size_t *icoord=gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, ndim);
size_t *ncoord=gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, ndim);
size_t *dsize = (correct_index ? tl->numtilesinch : input->dsize);
uint8_t *fullflag=&prm->thread_flags[tprm->id*input->size], *flag=fullflag;
- void *nv=gal_data_ptr_increment(prm->ngb_vals, tprm->id*prm->numneighbors,
- input->type);
+
/* Initializations. */
bb=prm->blanks->array;
@@ -106,10 +106,18 @@ interpolate_close_neighbors_on_thread(void *in_prm)
do *b = *bb++ ? INTERPOLATE_FLAGS_BLANK : 0; while(++b<bf);
- /* Put the nearest neighbor values into a structure for easy processing
- later. */
- nearest=gal_data_alloc(nv, input->type, 1, &prm->numneighbors, NULL, 0,
- 0, NULL, NULL, NULL);
+ /* Put the allocated space to keep the neighbor values into a structure
+ for easy processing. */
+ tin=input;
+ for(tvll=prm->ngb_vals; tvll!=NULL; tvll=tvll->next)
+ {
+ nv=gal_data_ptr_increment(tvll->v, tprm->id*prm->numneighbors,
+ input->type);
+ gal_data_add_to_ll(&nearest, nv, tin->type, 1, &prm->numneighbors,
+ NULL, 0, -1, NULL, NULL, NULL);
+ tin=tin->next;
+ }
+ gal_data_reverse_ll(&nearest);
/* Go over all the points given to this thread. */
@@ -125,12 +133,14 @@ interpolate_close_neighbors_on_thread(void *in_prm)
next element. */
if(prm->onlyblank && !(fullflag[fullind] & INTERPOLATE_FLAGS_BLANK) )
{
- /* It should be `input->array', not `values'! Because when
- channels are treated separately, `values' is going to
- change. */
- memcpy(gal_data_ptr_increment(out, fullind, type),
- gal_data_ptr_increment(input->array, fullind, type),
- twidth);
+ tin=input;
+ for(tout=prm->out; tout!=NULL; tout=tout->next)
+ {
+ memcpy(gal_data_ptr_increment(tout->array, fullind, tin->type),
+ gal_data_ptr_increment(tin->array, fullind, tin->type),
+ gal_data_sizeof(tin->type));
+ tin=tin->next;
+ }
continue;
}
@@ -150,14 +160,15 @@ interpolate_close_neighbors_on_thread(void *in_prm)
/* Index of the first tile in this channel. */
chstart = (fullind / tl->tottilesinch) * tl->tottilesinch;
- /* Correct the values and flag pointers so we can only work in
- this channel. */
- values = gal_data_ptr_increment(input->array, chstart, type);
+ /* Set the channel's starting pointer for the flags. */
flag = gal_data_ptr_increment(fullflag, chstart,
GAL_DATA_TYPE_UINT8);
}
else
- index=fullind;
+ {
+ chstart=0;
+ index=fullind;
+ }
/* Reset all checked bits in the flags array to 0. */
@@ -180,14 +191,21 @@ interpolate_close_neighbors_on_thread(void *in_prm)
/* Pop-out (p) an index from the queue: */
gal_linkedlist_pop_from_tosll_start(&lQ, &sQ, &pind, &pdist);
- /* If this isn't a blank value (recall that `index' was the first
- element added to the list which might be blank). */
+ /* If this isn't a blank value then add its values to the list of
+ neighbor values. Note that we didn't check whether the values
+ were blank or not when adding this pixel to the queue. */
if( !(flag[pind] & INTERPOLATE_FLAGS_BLANK) )
{
- /* Copy the value into the `nv' array. */
- memcpy(gal_data_ptr_increment(nv, ngb_counter, type),
- gal_data_ptr_increment(values, pind, type),
- twidth);
+ tin=input;
+ for(tnear=nearest; tnear!=NULL; tnear=tnear->next)
+ {
+ memcpy(gal_data_ptr_increment(tnear->array, ngb_counter,
+ tin->type),
+ gal_data_ptr_increment(tin->array, chstart+pind,
+ tin->type),
+ gal_data_sizeof(tin->type));
+ tin=tin->next;
+ }
/* If we have filled all the elements, break out. */
if(++ngb_counter>=prm->numneighbors) break;
@@ -227,21 +245,30 @@ interpolate_close_neighbors_on_thread(void *in_prm)
"interpolation", ngb_counter, prm->numneighbors);
}
- /* Calculate the median of the values and write it in. */
- single=gal_statistics_median(nearest, 1);
- memcpy(gal_data_ptr_increment(out, fullind, type), single->array,
- twidth);
-
- /* Clean up. */
- gal_data_free(single);
+ /* Calculate the median of the values and write it in the output. */
+ tout=prm->out;
+ for(tnear=nearest; tnear!=NULL; tnear=tnear->next)
+ {
+ /* Find the median and copy it. */
+ median=gal_statistics_median(tnear, 1);
+ memcpy(gal_data_ptr_increment(tout->array, fullind, tout->type),
+ median->array, gal_data_sizeof(tout->type));
+
+ /* Clean up and go to next array. */
+ gal_data_free(median);
+ tout=tout->next;
+ }
}
- /* Clean up, wait for all the other threads to finish and return. */
+ /* Clean up. */
+ for(tnear=nearest; tnear!=NULL; tnear=tnear->next) tnear->array=NULL;
+ gal_data_free_ll(nearest);
free(icoord);
free(ncoord);
- nearest->array=NULL;
- gal_data_free(nearest);
+
+
+ /* Wait for all the other threads to finish and return. */
if(tprm->b) pthread_barrier_wait(tprm->b);
return NULL;
}
@@ -257,29 +284,28 @@ gal_data_t *
gal_interpolate_close_neighbors(gal_data_t *input,
struct gal_tile_two_layer_params *tl,
size_t numneighbors, size_t numthreads,
- int onlyblank)
+ int onlyblank, int aslinkedlist)
{
+ char *name;
+ gal_data_t *tin, *tout;
struct interpolate_params prm;
+ size_t ngbvnum=numthreads*numneighbors;
int permute=(tl && tl->totchannels>1 && tl->workoverch);
/* Initialize the constant parameters. */
- prm.tl=tl;
- prm.input=input;
- prm.onlyblank=onlyblank;
- prm.numneighbors=numneighbors;
+ prm.tl = tl;
+ prm.ngb_vals = NULL;
+ prm.input = input;
+ prm.onlyblank = onlyblank;
+ prm.numneighbors = numneighbors;
+ prm.num = aslinkedlist ? gal_data_num_in_ll(input) : 1;
/* Flag the blank values. */
prm.blanks=gal_blank_flag(input);
- /* Allocate space for the output. */
- prm.out=gal_data_alloc(NULL, input->type, input->ndim, input->dsize,
- input->wcs, 0, input->minmapsize, "INTERPOLATED",
- input->unit, NULL);
-
-
/* If the input is from a tile structure and the user has asked to ignore
channels, then re-order the values. */
if(permute)
@@ -290,19 +316,60 @@ gal_interpolate_close_neighbors(gal_data_t *input,
/* Re-order values to ignore channels (if necessary). */
gal_permutation_apply(input, tl->permutation);
gal_permutation_apply(prm.blanks, tl->permutation);
+
+ /* If this is a linked list, then permute remaining nodes. */
+ if(aslinkedlist)
+ for(tin=input->next; tin!=NULL; tin=tin->next)
+ gal_permutation_apply(tin, tl->permutation);
}
+ /* Allocate space for the (first) output. */
+ name = input->name ? gal_checkset_malloc_cat("INTRP_", input->name) : NULL;
+ prm.out=gal_data_alloc(NULL, input->type, input->ndim, input->dsize,
+ input->wcs, 0, input->minmapsize, name, input->unit,
+ NULL);
+ free(name);
+ gal_linkedlist_add_to_vll(&prm.ngb_vals,
+ gal_data_malloc_array(input->type, ngbvnum));
+
+
+ /* If we are given a list of datasets, make the necessary
+ allocations. The reason we are doing this after a check of
+ `aslinkedlist' is that the `input' might have a `next' element, but
+ the caller might not have called `aslinkedlist'. */
+ prm.out->next=NULL;
+ if(aslinkedlist)
+ for(tin=input->next; tin!=NULL; tin=tin->next)
+ {
+ /* A small sanity check. */
+ if( gal_data_dsize_is_different(input, tin) )
+ error(EXIT_FAILURE, 0, "The other datasets in the list must "
+ "have the same dimension and size in "
+ "`gal_interpolate_close_neighbors'");
+
+ /* Allocate the output array for this node. */
+ name = ( tin->name
+ ? gal_checkset_malloc_cat("INTRP_", tin->name) : NULL );
+ gal_data_add_to_ll(&prm.out, NULL, tin->type, tin->ndim, tin->dsize,
+ tin->wcs, 0, tin->minmapsize, name, tin->unit,
+ NULL);
+ if(name) free(name);
+
+ /* Allocate the space for the neighbor values of this input. */
+ gal_linkedlist_add_to_vll(&prm.ngb_vals,
+ gal_data_malloc_array(tin->type, ngbvnum));
+ }
+ gal_data_reverse_ll(&prm.out);
+ gal_linkedlist_reverse_vll(&prm.ngb_vals);
+
+
/* Allocate space for all the flag values of all the threads here (memory
in each thread is limited) and this is cleaner. */
prm.thread_flags=gal_data_malloc_array(GAL_DATA_TYPE_UINT8,
numthreads*input->size);
- /* Allocate space for the list of neighbor values in each thread. */
- prm.ngb_vals=gal_data_malloc_array(input->type, numthreads*numneighbors);
-
-
/* Spin off the threads. */
gal_threads_spin_off(interpolate_close_neighbors_on_thread, &prm,
input->size, numthreads);
@@ -314,7 +381,8 @@ gal_interpolate_close_neighbors(gal_data_t *input,
if(permute)
{
gal_permutation_apply_inverse(input, tl->permutation);
- gal_permutation_apply_inverse(prm.out, tl->permutation);
+ for(tout=prm.out; tout!=NULL; tout=tout->next)
+ gal_permutation_apply_inverse(tout, tl->permutation);
}
diff --git a/lib/linkedlist.c b/lib/linkedlist.c
index b02b873..30193fe 100644
--- a/lib/linkedlist.c
+++ b/lib/linkedlist.c
@@ -820,6 +820,96 @@ gal_linkedlist_free_ill(struct gal_linkedlist_ill *list)
/****************************************************************
+ ***************** void * ********************
+ ****************************************************************/
+void
+gal_linkedlist_add_to_vll(struct gal_linkedlist_vll **list, void *value)
+{
+ struct gal_linkedlist_vll *newnode;
+
+ errno=0;
+ newnode=malloc(sizeof *newnode);
+ if(newnode==NULL)
+ error(EXIT_FAILURE, errno,
+ "linkedlist: New element in gal_linkedlist_ill");
+
+ newnode->v=value;
+ newnode->next=*list;
+ *list=newnode;
+}
+
+
+
+
+
+void
+gal_linkedlist_pop_from_vll(struct gal_linkedlist_vll **list, void **value)
+{
+ struct gal_linkedlist_vll *tmp;
+ tmp=*list;
+ *value=tmp->v;
+ *list=tmp->next;
+ free(tmp);
+}
+
+
+
+
+
+void
+gal_linkedlist_reverse_vll(struct gal_linkedlist_vll **list)
+{
+ void *thisptr;
+ struct gal_linkedlist_vll *correctorder=NULL;
+
+ if( *list && (*list)->next )
+ {
+ while(*list!=NULL)
+ {
+ gal_linkedlist_pop_from_vll(list, &thisptr);
+ gal_linkedlist_add_to_vll(&correctorder, thisptr);
+ }
+ *list=correctorder;
+ }
+}
+
+
+
+
+
+void
+gal_linkedlist_free_vll(struct gal_linkedlist_vll *list, int freevalue)
+{
+ struct gal_linkedlist_vll *tmp=list, *ttmp;
+ while(tmp!=NULL)
+ {
+ if(freevalue) free(tmp->v);
+ ttmp=tmp->next;
+ free(tmp);
+ tmp=ttmp;
+ }
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/****************************************************************
****************** Ordered SLL ********************
***************** size_t ********************
****************************************************************/
diff --git a/lib/options.c b/lib/options.c
index a6e9ba8..4ceab43 100644
--- a/lib/options.c
+++ b/lib/options.c
@@ -594,6 +594,81 @@ gal_options_parse_sizes_reverse(struct argp_option
*option, char *arg,
+/* Two numbers must be provided as an argument. This function will read
+ them as the sigma-clipping multiple and parameter and store the two in a
+ 2-element array. `option->value' must point to an already allocated
+ 2-element array of double type. */
+void *
+gal_options_read_sigma_clip(struct argp_option *option, char *arg,
+ char *filename, size_t lineno, void *junk)
+{
+ char *str;
+ gal_data_t *in;
+ double *sigmaclip=option->value;
+
+ /* Caller wants to print the option values. */
+ if(lineno==-1)
+ {
+ asprintf(&str, "%g,%g", sigmaclip[0], sigmaclip[1]);
+ return str;
+ }
+
+ /* Caller wants to read the values into memory, so parse the inputs. */
+ in=gal_options_parse_list_of_numbers(arg, filename, lineno);
+
+ /* Check if there was only two numbers. */
+ if(in->size!=2)
+ error_at_line(EXIT_FAILURE, 0, filename, lineno, "the `--%s' "
+ "option takes two values (separated by a comma) for "
+ "defining the sigma-clip. However, %zu numbers were "
+ "read in the string `%s' (value to this option).\n\n"
+ "The first number is the multiple of sigma, and the "
+ "second is either the tolerance (if its is less than "
+ "1.0), or a specific number of times to clip (if it "
+ "is equal or larger than 1.0).", option->name, in->size,
+ arg);
+
+ /* Copy the sigma clip parameters into the space the caller has given (as
+ the `value' element of `option'). */
+ memcpy(option->value, in->array, 2*sizeof *sigmaclip);
+
+ /* Multiple of sigma must be positive. */
+ if( sigmaclip[0] <= 0 )
+ error_at_line(EXIT_FAILURE, 0, filename, lineno, "the first value to "
+ "the `--%s' option (multiple of sigma), must be "
+ "greater than zero. From the string `%s' (value to "
+ "this option), you have given a value of %g for the "
+ "first value", option->name, arg, sigmaclip[0]);
+
+ /* Second value must also be positive. */
+ if( sigmaclip[1] <= 0 )
+ error_at_line(EXIT_FAILURE, 0, filename, lineno, "the second value "
+ "to the `--%s' option (tolerance to stop clipping or "
+ "number of clips), must be greater than zero. From "
+ "the string `%s' (value to this option), you have "
+ "given a value of %g for the second value",
+ option->name, arg, sigmaclip[1]);
+
+ /* if the second value is larger or equal to 1.0, it must be an
+ integer. */
+ if( sigmaclip[1] >= 1.0f && ceil(sigmaclip[1]) != sigmaclip[1])
+ error_at_line(EXIT_FAILURE, 0, filename, lineno, "when the second "
+ "value to the `--%s' option is >=1, it is interpretted "
+ "as an absolute number of clips. So it must be an "
+ "integer. However, your second value is a floating "
+ "point number: %g (parsed from `%s')", option->name,
+ sigmaclip[1], arg);
+
+ /* Clean up and return. */
+ gal_data_free(in);
+ return NULL;
+}
+
+
+
+
+
+
diff --git a/lib/options.h b/lib/options.h
index 293bf6c..c6cf69f 100644
--- a/lib/options.h
+++ b/lib/options.h
@@ -106,6 +106,8 @@ enum options_common_keys
GAL_OPTIONS_KEY_WORKOVERCH,
GAL_OPTIONS_KEY_CHECKTILES,
GAL_OPTIONS_KEY_ONEELEMPERTILE,
+ GAL_OPTIONS_KEY_INTERPONLYBLANK,
+ GAL_OPTIONS_KEY_INTERPNUMNGB,
};
@@ -161,6 +163,8 @@ struct gal_options_common_params
{
/* Tessellation. */
struct gal_tile_two_layer_params tl; /* Two layer tessellation params. */
+ uint8_t interponlyblank; /* Only interpolate over blank values. */
+ size_t interpnumngb; /* Number of neighbors for interpolation. */
/* Input. */
char *hdu; /* Image extension. */
@@ -254,7 +258,9 @@ void *
gal_options_parse_sizes_reverse(struct argp_option *option, char *arg,
char *filename, size_t lineno, void *params);
-
+void *
+gal_options_read_sigma_clip(struct argp_option *option, char *arg,
+ char *filename, size_t lineno, void *junk);
/**********************************************************************/
diff --git a/lib/statistics.c b/lib/statistics.c
index 71b6843..9a236f5 100644
--- a/lib/statistics.c
+++ b/lib/statistics.c
@@ -522,9 +522,11 @@ mode_mirror_max_index_diff(struct statistics_mode_params
*p, size_t m)
/*
printf("###############\n###############\n");
- printf("### Mirror pixel: %zu\n", m);
+ printf("### Mirror pixel: %zu (mirrordist: %f, sqrt(m): %f)\n", m,
+ p->mirrordist, sqrt(m));
printf("###############\n###############\n");
*/
+
/* Go over the mirrored points. */
for(i=1; i<p->numcheck && i<=m && m+i<size ;i+=p->interval)
{
@@ -552,6 +554,7 @@ mode_mirror_max_index_diff(struct statistics_mode_params
*p, size_t m)
printf("i:%-5zu j:%-5zu diff:%-5d maxdiff: %zu\n",
i, j, (int)j-(int)i, maxdiff);
*/
+
/* The index of the actual CDF corresponding the the mirrored flux
has been found. We want the mirrored distribution to be within the
actual distribution, not beyond it, so the only acceptable results
@@ -633,7 +636,8 @@ mode_golden_section(struct statistics_mode_params *p)
-------------------------------------------------------------------*/
/*
- printf("%-5zu\t%-5zu(%d)\t%-5zu ----> dq: %-5zu di: %d\n",
+ printf("lowi:%-5zu\tmidi:%-5zu(midd: %d)\thighi:%-5zu ----> "
+ "dq: %-5zu di: %d\n",
p->lowi, p->midi, (int)p->midd, p->highi,
di, (int)dd);
*/
@@ -831,6 +835,13 @@ gal_statistics_mode(gal_data_t *input, float mirrordist,
int inplace)
gal_data_t *out=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT64, 1, &dsize,
NULL, 1, -1, NULL, NULL, NULL);
+
+ /* A small sanity check. */
+ if(mirrordist<=0)
+ error(EXIT_FAILURE, 0, "%f not acceptable as a value to `mirrordist'. "
+ "Only positive values can be given to it", mirrordist);
+
+
/* Make sure the input doesn't have blank values and is sorted. */
p.data=gal_statistics_no_blank_sorted(input, inplace);
@@ -876,9 +887,14 @@ gal_statistics_mode(gal_data_t *input, float mirrordist,
int inplace)
oa[2] = mode_symmetricity(&p, modeindex, b_val->array);
- /* Put the end of the symmetricity in. */
- b_val=gal_data_copy_to_new_type_free(mode, GAL_DATA_TYPE_FLOAT64);
- gal_data_copy_element_same_type(mode, 0, &oa[3]);
+ /* If the symmetricity is good, put it in the output, otherwise set all
+ output values to NaN. */
+ if(oa[2]>GAL_STATISTICS_MODE_GOOD_SYM)
+ {
+ b_val=gal_data_copy_to_new_type_free(mode, GAL_DATA_TYPE_FLOAT64);
+ gal_data_copy_element_same_type(b_val, 0, &oa[3]);
+ }
+ else oa[0]=oa[1]=oa[2]=oa[3]=NAN;
/* For a check:
@@ -1152,8 +1168,10 @@ gal_statistics_sort_decreasing(gal_data_t *data)
gal_data_t *
gal_statistics_no_blank_sorted(gal_data_t *input, int inplace)
{
+ int sortstatus;
gal_data_t *contig, *noblank, *sorted;
+
/* If this is a tile, then first we have to copy it into a contiguous
piece of memory. After this step, we will only be dealing with
`contig' (for a contiguous patch of memory). */
@@ -1169,6 +1187,7 @@ gal_statistics_no_blank_sorted(gal_data_t *input, int
inplace)
}
else contig=input;
+
/* Make sure there is no blanks in the array that will be used. After
this step, we won't be dealing with `input' any more, but with
`noblank'.*/
@@ -1187,9 +1206,15 @@ 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'. */
- if( gal_statistics_is_sorted(noblank) ) sorted=noblank;
+ sortstatus=gal_statistics_is_sorted(noblank);
+ if( sortstatus )
+ {
+ sorted=noblank;
+ sorted->status=sortstatus;
+ }
else
{
if(inplace) sorted=noblank;
@@ -1201,8 +1226,11 @@ gal_statistics_no_blank_sorted(gal_data_t *input, int
inplace)
sorted=gal_data_copy(noblank);
}
gal_statistics_sort_increasing(sorted);
+ sorted->status=GAL_STATISTICS_SORTED_INCREASING;
}
+
+ /* Return final array. */
return sorted;
}
@@ -1637,9 +1665,10 @@ gal_statistics_cfp(gal_data_t *data, gal_data_t *bins,
int normalize)
/****************************************************************
***************** Outliers ********************
****************************************************************/
-/* Return a data structure with an array of four values: the final number
- of points used, the median, average and standard deviation. The number
- of clips is put into the `status' element of the data structure.
+/* Return a float32 type data structure with an array of four values: the
+ final number of points used, the median, average and standard
+ deviation. The number of clips is put into the `status' element of the
+ data structure.
Inputs:
@@ -1657,10 +1686,9 @@ gal_statistics_cfp(gal_data_t *data, gal_data_t *bins,
int normalize)
point of the array and its size, calcluating the basic statistics in each
round to define the new starting point and size.
*/
-
#define SIGCLIP(IT) { \
- IT *a = sorted->array, *af = a + sorted->size; \
- IT *bf = sorted->array, *b = bf + sorted->size - 1; \
+ IT *a = nbs->array, *af = a + nbs->size; \
+ IT *bf = nbs->array, *b = bf + nbs->size - 1; \
\
/* Remove all out-of-range elements from the start of the array. */ \
if(sortstatus==GAL_STATISTICS_SORTED_INCREASING) \
@@ -1685,17 +1713,19 @@ gal_statistics_cfp(gal_data_t *data, gal_data_t *bins,
int normalize)
gal_data_t *
gal_statistics_sigma_clip(gal_data_t *input, float multip, float param,
- int quiet)
+ int inplace, int quiet)
{
- int sortstatus;
+ void *start, *nbs_array;
double *med, *mean, *std;
- void *start, *sorted_array;
- size_t num=0, dsize=4, size, oldsize;
uint8_t bytolerance = param>=1.0f ? 0 : 1;
double oldmed=NAN, oldmean=NAN, oldstd=NAN;
- gal_data_t *sorted=NULL, *median_i, *median_d, *out, *meanstd, *noblank;
+ size_t num=0, one=1, four=4, size, oldsize;
+ gal_data_t *median_i, *median_d, *out, *meanstd;
+ int sortstatus, type=gal_tile_block(input)->type;
+ gal_data_t *nbs=gal_statistics_no_blank_sorted(input, inplace);
size_t maxnum = param>=1.0f ? param : GAL_STATISTICS_SIG_CLIP_MAX_CONVERGE;
+
/* Some sanity checks. */
if( multip<=0 )
error(EXIT_FAILURE, 0, "`multip', must be greater than zero in "
@@ -1710,43 +1740,11 @@ gal_statistics_sigma_clip(gal_data_t *input, float
multip, float param,
"your given value %g", param);
- /* If there are blank elements, remove them (from a copied array). NOTE:
- we don't want to change the input. */
- if( gal_blank_present(input) )
- {
- noblank = gal_data_copy(input);
- gal_blank_remove(noblank);
- }
- else noblank=input;
-
-
- /* We want to have sorted (increasing) array. */
- sortstatus=gal_statistics_is_sorted(noblank);
- switch(sortstatus)
- {
- case GAL_STATISTICS_SORTED_NOT:
- /* Only copy (or allocate) space if it wasn't already allocated. */
- sorted = (noblank==input) ? gal_data_copy(input) : noblank;
- gal_statistics_sort_increasing(sorted);
- sortstatus=GAL_STATISTICS_SORTED_INCREASING;
- break;
-
- case GAL_STATISTICS_SORTED_DECREASING:
- case GAL_STATISTICS_SORTED_INCREASING:
- sorted=noblank;
- break;
-
- default:
- error(EXIT_FAILURE, 0, "sorted code %d not recognized in "
- "`gal_statistics_sigma_clip'", gal_statistics_is_sorted(input));
- }
-
-
/* Allocate the necessary spaces. */
- out=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT32, 1, &dsize, NULL, 0,
+ out=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT32, 1, &four, NULL, 0,
input->minmapsize, NULL, NULL, NULL);
- median_i=gal_data_alloc(NULL, sorted->type, 1, &dsize, NULL, 0,
- input->minmapsize, NULL, NULL, NULL);
+ median_i=gal_data_alloc(NULL, type, 1, &one, NULL, 0, input->minmapsize,
+ NULL, NULL, NULL);
/* Print the comments. */
@@ -1757,19 +1755,20 @@ gal_statistics_sigma_clip(gal_data_t *input, float
multip, float param,
/* 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=sorted->size;
- sorted_array=start=sorted->array;
+ size=nbs->size;
+ sortstatus=nbs->status;
+ nbs_array=start=nbs->array;
while(num<maxnum)
{
/* Find the median. */
- statistics_median_in_sorted_no_blank(sorted, median_i->array);
+ statistics_median_in_sorted_no_blank(nbs, median_i->array);
median_d=gal_data_copy_to_new_type(median_i, GAL_DATA_TYPE_FLOAT64);
/* Find the average and Standard deviation, note that both `start'
and `size' will be different in the next round. */
- sorted->array = start;
- sorted->size = oldsize = size;
- meanstd=gal_statistics_mean_std(sorted);
+ 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;
@@ -1792,7 +1791,7 @@ gal_statistics_sigma_clip(gal_data_t *input, float
multip, float param,
/* 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(sorted->type)
+ switch(type)
{
case GAL_DATA_TYPE_UINT8: SIGCLIP( uint8_t ); break;
case GAL_DATA_TYPE_INT8: SIGCLIP( int8_t ); break;
@@ -1806,7 +1805,7 @@ gal_statistics_sigma_clip(gal_data_t *input, float
multip, float param,
case GAL_DATA_TYPE_FLOAT64: SIGCLIP( double ); break;
default:
error(EXIT_FAILURE, 0, "type code %d not recognized in "
- "`gal_statistics_sigma_clip'", sorted->type);
+ "`gal_statistics_sigma_clip'", type);
}
/* Set the values from this round in the old elements, so the next
@@ -1840,8 +1839,8 @@ gal_statistics_sigma_clip(gal_data_t *input, float
multip, float param,
}
/* Clean up and return. */
- sorted->array=sorted_array;
- if(sorted!=input) gal_data_free(sorted);
+ nbs->array=nbs_array;
+ if(nbs!=input) gal_data_free(nbs);
return out;
}
diff --git a/lib/tile.c b/lib/tile.c
index 330fe02..3a6dd7c 100644
--- a/lib/tile.c
+++ b/lib/tile.c
@@ -411,6 +411,21 @@ gal_tile_block_check_tiles(gal_data_t *tilesll)
+/* Return the pointer corresponding to the tile in another data
+ structure (can have another type). */
+void *
+gal_tile_block_relative_to_other(gal_data_t *tile, gal_data_t *other)
+{
+ gal_data_t *block=gal_tile_block(tile);
+ return gal_data_ptr_increment(other->array,
+ gal_data_ptr_dist(block->array,
+ tile->array, block->type),
+ other->type);
+}
+
+
+
+
@@ -564,10 +579,10 @@ gal_tile_full(gal_data_t *input, size_t *regular,
size_t i, d, tind, numtiles, *start=NULL;
gal_data_t *tiles, *block=gal_tile_block(input);
size_t *last = gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, input->ndim);
- size_t *tsize = gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, input->ndim);
size_t *first = gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, input->ndim);
size_t *coord = gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, input->ndim);
size_t *tcoord = gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, input->ndim);
+ size_t *tsize = gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, input->ndim+1);
/* Set the first tile size and total number of tiles along each
@@ -671,7 +686,8 @@ gal_tile_full(gal_data_t *input, size_t *regular,
free(coord);
free(tcoord);
if(start) free(start);
- return tsize;
+ tsize[input->ndim]=-1; /* So it can be used for another tessellation, */
+ return tsize; /* see `gal_tile_full_sanity_check'. */
}
@@ -717,7 +733,7 @@ gal_tile_full_sanity_check(char *filename, char *hdu,
gal_data_t *input,
if(tl->numchannels[i]==0)
error(EXIT_FAILURE, 0, "the number of channels in all dimensions must "
"be larger than zero. The number for dimension %zu was zero",
- ndim);
+ i+1);
if(i!=ndim)
error(EXIT_FAILURE, 0, "%s (hdu: %s): has %zu dimensions, but only %zu "
"value(s) given for the number of channels", filename, hdu, ndim,
@@ -782,16 +798,19 @@ gal_tile_full_two_layers(gal_data_t *input,
gal_data_t *t;
size_t i, *junk, ndim=tl->ndim=input->ndim;
- /* Initialize. Note that `numchannels might have already been
- allocated. */
+ /* Initialize. */
tl->channels=tl->tiles=NULL;
- if(tl->numchannels) free(tl->numchannels);
/* Initialize necessary values and do the channels tessellation. */
- tl->numchannels = gal_tile_full(input, tl->channelsize, tl->remainderfrac,
- &tl->channels, 1);
+ junk = gal_tile_full(input, tl->channelsize, tl->remainderfrac,
+ &tl->channels, 1);
tl->totchannels = gal_dimension_total_size(ndim, tl->numchannels);
-
+ for(i=0;i<ndim;++i)
+ if(junk[i]!=tl->numchannels[i])
+ error(EXIT_FAILURE, 0, "the input and output number of channels in "
+ "`gal_tile_full_two_layers' don't match in dimension %zu, with "
+ "values of %zu and %zu respectively.", ndim-i,
+ tl->numchannels[i], junk[i]);
/* Tile each channel. While tiling the first channel, we are also going
to allocate the space for the other channels. Then pass those pointers
@@ -946,12 +965,20 @@ gal_tile_full_values_write(gal_data_t *tilevalues,
{
gal_data_t *disp;
-
/* Make the dataset to be displayed. */
if(tl->oneelempertile)
{
if(tl->ndim>1 && tl->totchannels>1)
{
+ /* A small sanity check. */
+ if(tl->permutation==NULL)
+ error(EXIT_FAILURE, 0, "no permutation defined for the input "
+ "tessellation to `gal_tile_full_values_write'");
+
+ /* Writing tile values to disk is not done for checking, not for
+ efficiency. So to be safe (allow the caller to work on
+ multiple threads), we will copy the tile values, then permute
+ those. */
disp = gal_data_copy(tilevalues);
gal_permutation_apply(disp, tl->permutation);
}
@@ -970,6 +997,74 @@ gal_tile_full_values_write(gal_data_t *tilevalues,
+/* Smooth the given values with a flat kernel of the given width. */
+gal_data_t *
+gal_tile_full_values_smooth(gal_data_t *tilevalues,
+ struct gal_tile_two_layer_params *tl,
+ size_t width, size_t numthreads)
+{
+ size_t *kdsize, knum, i;
+ gal_data_t *kernel, *smoothed;
+ struct gal_tile_two_layer_params ttl={0};
+ int permute=tl->ndim>1 && tl->totchannels>1;
+
+
+ /* Check if the width is odd. */
+ if(width%2==0)
+ error(EXIT_FAILURE, 0, "%zu not acceptable as width to "
+ "`gal_tile_full_values_smooth'. It has to be an odd number", width);
+
+
+ /* Prepare the kernel size along every dimension. */
+ kdsize=gal_data_malloc_array(GAL_DATA_TYPE_SIZE_T, tl->ndim);
+ for(i=0;i<tl->ndim;++i) kdsize[i]=width;
+
+
+ /* Make the kernel. */
+ kernel=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT32, tilevalues->ndim,
+ kdsize, NULL, 0, -1, NULL, NULL, NULL);
+ knum=gal_dimension_total_size(tl->ndim, kernel->dsize);
+ for(i=0;i<knum;++i) ((float *)(kernel->array))[i]=1/((double)knum);
+
+ /* Permute (if necessary). */
+ if(permute)
+ {
+ gal_tile_full_permutation(tl);
+ gal_permutation_apply(tilevalues, tl->permutation);
+ }
+
+ /* Do the smoothing. */
+ if(tl->workoverch)
+ smoothed=gal_convolve_spatial(tilevalues, kernel, numthreads, 1, 1);
+ else
+ {
+ /* Create the tile structure. */
+ ttl.tilesize=tl->numtilesinch;
+ ttl.numchannels=tl->numchannels;
+ gal_tile_full_sanity_check("IMPOSSIBLE", "IMP_HDU", tilevalues, &ttl);
+ gal_tile_full_two_layers(tilevalues, &ttl);
+
+ /* Do the convolution separately on each channel. */
+ smoothed=gal_convolve_spatial(ttl.tiles, kernel, numthreads, 1, 0);
+
+ /* Clean up. */
+ ttl.tilesize=ttl.numchannels=NULL;
+ gal_tile_full_free_contents(&ttl);
+ }
+
+ /* Reverse the permutation. */
+ if(permute) gal_permutation_apply_inverse(smoothed, tl->permutation);
+
+ /* Clean up and return; */
+ free(kdsize);
+ gal_data_free(kernel);
+ return smoothed;
+}
+
+
+
+
+
/* Clean up the allocated spaces in the parameters. */
void
gal_tile_full_free_contents(struct gal_tile_two_layer_params *tl)
diff --git a/tests/Makefile.am b/tests/Makefile.am
index a1837d2..167ba79 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -140,14 +140,10 @@ if COND_NOISECHISEL
noisechisel/noisechisel.sh: mknoise/addnoise.sh.log
endif
if COND_STATISTICS
- MAYBE_STATISTICS_TESTS = statistics/basicstats.sh
+ MAYBE_STATISTICS_TESTS = statistics/basicstats.sh statistics/estimate_sky.sh
- statistics/basicstats.sh: mknoise/addnoise.sh.log
-endif
-if COND_SUBTRACTSKY
- MAYBE_SUBTRACTSKY_TESTS = subtractsky/subtractsky.sh
-
- subtractsky/subtractsky.sh: mknoise/addnoise.sh.log
+ statistics/basicstats.sh: mknoise/addnoise.sh.log
+ statistics/estimate_sky.sh: mknoise/addnoise.sh.log
endif
if COND_TABLE
MAYBE_TABLE_TESTS = table/txt-to-fits-binary.sh \
diff --git a/tests/during-dev.sh b/tests/during-dev.sh
index 5918d10..eaa2399 100755
--- a/tests/during-dev.sh
+++ b/tests/during-dev.sh
@@ -85,8 +85,6 @@ options=
-
-
# RUN THE PROCEDURES
# ==================
@@ -140,7 +138,7 @@ if make -C "$builddir"; then
# the last line in the configuration file doesn't actualy end with a
# new line (in which case the appended string will be added to the end
# of the last line).
- cp "$srcdir/lib/gnuastro.conf" "$srcdir/bin/$utilname/ast$utilname.conf" \
+ cp "$srcdir/bin/gnuastro.conf" "$srcdir/bin/$utilname/ast$utilname.conf" \
.gnuastro/
echo "" >> .gnuastro/gnuastro.conf
echo " lastconfig 1" >> .gnuastro/gnuastro.conf
diff --git a/tests/prepconf.sh b/tests/prepconf.sh
index f62f105..d5a1b76 100755
--- a/tests/prepconf.sh
+++ b/tests/prepconf.sh
@@ -61,7 +61,7 @@ cat > addedoptions.txt <<EOF
lastconfig 1
log 0
EOF
-cat $topsrc/lib/gnuastro.conf addedoptions.txt > .gnuastro/gnuastro.conf
+cat $topsrc/bin/gnuastro.conf addedoptions.txt > .gnuastro/gnuastro.conf
rm addedoptions.txt
@@ -75,7 +75,7 @@ rm addedoptions.txt
# by `make check'.
for prog in arithmetic convertt convolve cosmiccal crop fits \
mkcatalog mknoise mkprof noisechisel statistics \
- subtractsky table warp
+ table warp
do
cp $topsrc/bin/$prog/ast$prog.conf .gnuastro/ast$prog.conf
done
diff --git a/tests/subtractsky/subtractsky.sh b/tests/statistics/estimate_sky.sh
similarity index 77%
rename from tests/subtractsky/subtractsky.sh
rename to tests/statistics/estimate_sky.sh
index 1d3db7b..086dc7f 100755
--- a/tests/subtractsky/subtractsky.sh
+++ b/tests/statistics/estimate_sky.sh
@@ -1,4 +1,4 @@
-# Add noise to a large image then subtract the sky.
+# Estimate the Sky value on input image.
#
# See the Tests subsection of the manual for a complete explanation
# (in the Installing gnuastro section).
@@ -22,7 +22,7 @@
# Set the variabels (The executable is in the build tree). Do the
# basic checks to see if the executable is made or if the defaults
# file exists (basicchecks.sh is in the source tree).
-prog=subtractsky
+prog=statistics
execname=../bin/$prog/ast$prog
img=convolve_spatial_noised.fits
@@ -48,4 +48,9 @@ if [ ! -f $execname ] || [ ! -f $img ]; then exit 77; fi
# Actual test script
# ==================
-$execname $img --checksky --numnearest=4
+#
+# Note that to keep things simple we are not convolving the image, so the
+# result will not be too accurate! Here we just want to see if the full
+# tessellation, estimation, interpolation and smoothing go nicely without
+# any errors.
+$execname $img --sky --checksky
diff --git a/tmpfs-config-make b/tmpfs-config-make
index 9f31ef0..1be490f 100755
--- a/tmpfs-config-make
+++ b/tmpfs-config-make
@@ -130,12 +130,12 @@ cd $build_dir
#
# ####################################
if [ ! -f Makefile ]; then
- $srcdir/configure --srcdir=$srcdir
\
- --enable-arithmetic --enable-convertt
--enable-convolve\
- --enable-cosmiccal --enable-crop --enable-fits
\
- --enable-mknoise --enable-statistics --enable-mkprof
\
- --enable-table --enable-warp
\
- --disable-shared CFLAGS="-g -O0"
+ $srcdir/configure --srcdir=$srcdir \
+ --enable-arithmetic --enable-convertt --enable-convolve \
+ --enable-cosmiccal --enable-crop --enable-fits \
+ --enable-mknoise --enable-statistics --enable-mkprof \
+ --enable-table --enable-warp \
+ --disable-shared CFLAGS="-g -O0"
fi
- [gnuastro-commits] master 60c126f 120/125: New implementation of over-segmentation to find clumps, (continued)
- [gnuastro-commits] master 60c126f 120/125: New implementation of over-segmentation to find clumps, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master c250a6d 123/125: NoiseChisel's clumps grown and objects defined, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master ffaa059 101/125: Tessellation options now common to all programs, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master dec2b22 078/125: ImageCrop works with gal_data_t features, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master f8b1820 061/125: Other configuration files read immediately when seen, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 46b4651 085/125: Image prefix removed from some program names, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master df716bc 097/125: Tesselation with the new data structure: work started, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master db39cee 088/125: Statistics returns basic results, no array lib, new blank lib, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 60cd76e 090/125: Statistics mostly complete, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 15c524d 124/125: New implementation of NoiseChisel complete, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 3b26d5c 114/125: SubtractSky program removed, Statistics estimates Sky,
Mohammad Akhlaghi <=
- [gnuastro-commits] master 99ac7dc 083/125: Gnuastro now uses fixed width (from stdint.h) types, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master eaf19a3 115/125: NoiseChisel's initial detection step implemented, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 3339b92 125/125: New implementation of MakeCatalog now complete, Mohammad Akhlaghi, 2017/04/23