[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 4a08b3f: NoiseChisel & Statistics: improved ou
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 4a08b3f: NoiseChisel & Statistics: improved outlier rejection in Sky |
|
Date: |
Fri, 4 Dec 2020 01:17:44 -0500 (EST) |
branch: master
commit 4a08b3fc205f966a00b42607c33156638fb20707
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>
NoiseChisel & Statistics: improved outlier rejection in Sky
Until now, once the detection was done, the tiles that weren't
significantly covered by detections (based on '--minskyfrac') would be
directly used for the interpolation phase: there was no outlier rejection
beyond the median-based interpolation. Also, when removing outliers during
the finding of the quantile threshold, we were using global measures to
find and reject outliers. So it wouldn't opeate too well when there were
strong gradients for example.
With this commit, to find outliers in the quantile threshold, we now
calculate the distance between the second-smallest and second-largest
values in the '--interpnumngb' tiles around it. This allows for a more
local rejection of some tiles. This is done in the newly added
'gal_tileinternal_no_outlier_local' funciton. We then look at the
distribution of those values to find outliers. Also, outlier rejection is
now also implemented on the final Sky image.
In the process, some other issues have been fixed:
- While looking at previous code to remember the low-level implementation
details, I noticed that the lines are longer than 75 characters, making
the functions hard to read.
- The internal 'INTERPOLATE_FLAGS_CHECKED' macro in 'lib/interpolate.c'
was renamed to 'INTERPOLATE_FLAGS_NGB_CHECKED' to be more clear. When
reviewing the code, initially I mistakenly thought it was the fact that
a pixel is checked or not. But it was actually a check to see if the
pixel has been checked during the search to find the nearest neighbor.
- In the 'gal_statistics_quantile_function_index' function, I noticed that
we weren't accounting for cases where the lowest two elements have the
same value (to see if the array is sorted in increasing or decreasing
order). As a result, some results were being mistakenly returned as
'inf'. With this commit, we are now using the relevant flag that
explicitly keeps this information.
- Added sanity checks on the 'value' argument in
'gal_statistics_quantile_function' and its index finding function.
---
NEWS | 13 +
bin/gnuastro.conf | 2 +-
bin/noisechisel/astnoisechisel.conf | 2 +-
bin/noisechisel/sky.c | 20 +-
bin/noisechisel/threshold.c | 24 +-
bin/statistics/aststatistics.conf | 2 +-
bin/statistics/sky.c | 8 +-
doc/gnuastro.texi | 93 ++++---
lib/gnuastro-internal/tile-internal.h | 8 +
lib/gnuastro/interpolate.h | 8 +-
lib/interpolate.c | 50 ++--
lib/statistics.c | 33 ++-
lib/tile-internal.c | 481 +++++++++++++++++++++++++++++++++-
tests/during-dev.sh | 2 +
tests/noisechisel/noisechisel.sh | 2 +-
15 files changed, 646 insertions(+), 102 deletions(-)
diff --git a/NEWS b/NEWS
index 32268ab..da28615 100644
--- a/NEWS
+++ b/NEWS
@@ -131,12 +131,25 @@ See the end of the file for license conditions.
(which can dramatically slow down the program). Please see the newly
added "Memory management" section of the book for a complete
explanation of Gnuastro's new memory management strategy.
+ --interpnumngb: the default value has been increased to 15 (from 9). The
+ reason for this is that we now have a more robust outlier removal
+ algorithm (see description under "NoiseChisel & Statistics").
Fits:
- The '--pixelscale' option also prints the pixel area (for 2D inputs,
or 2D slices of 3D inputs) and the voxel volume (for 3D inputs). Until
now, it would only print the pixel scale along each dimension.
+ NoiseChisel & Statistics:
+ - New algorithm used to reject outlying tiles. In NoiseChisel this is
+ done when estimating the quantile threshold, the pseudo-detection
+ threshold and the final Sky value. In Statistics, its just the Sky
+ value. Unlike the previous method that used the global distribution of
+ tile values to identify outliers, the new algorithm uses a relative
+ measure. See the book for more. Since outliers are now rejected more
+ robustly, the default value of 'outliersigma' has been decreased to
+ '5' (previously it was 10).
+
Table:
- Column arithmetic operators 'degree-to-ra' and 'degree-to-dec' will
return the sexagesimal format of '_h_m_s' and '_d_m_s'
diff --git a/bin/gnuastro.conf b/bin/gnuastro.conf
index 40d5f8f..d7b86be 100644
--- a/bin/gnuastro.conf
+++ b/bin/gnuastro.conf
@@ -31,7 +31,7 @@
remainderfrac 0.1
workoverch 0
interpmetric radial
- interpnumngb 9
+ interpnumngb 15
interponlyblank 0
# Output:
diff --git a/bin/noisechisel/astnoisechisel.conf
b/bin/noisechisel/astnoisechisel.conf
index add04ee..74dc854 100644
--- a/bin/noisechisel/astnoisechisel.conf
+++ b/bin/noisechisel/astnoisechisel.conf
@@ -30,7 +30,7 @@
# Detection:
meanmedqdiff 0.005
qthresh 0.3
- outliersigma 10
+ outliersigma 5
outliersclip 3,0.2
smoothwidth 3
erode 2
diff --git a/bin/noisechisel/sky.c b/bin/noisechisel/sky.c
index f200969..b2ca9d2 100644
--- a/bin/noisechisel/sky.c
+++ b/bin/noisechisel/sky.c
@@ -197,12 +197,12 @@ sky_and_std(struct noisechiselparams *p, char *checkname)
/* Allocate space for the mean and standard deviation. */
- p->sky=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, p->input->ndim, tl->numtiles,
- NULL, 0, cp->minmapsize, p->cp.quietmmap, NULL,
- p->input->unit, NULL);
- p->std=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, p->input->ndim, tl->numtiles,
- NULL, 0, cp->minmapsize, p->cp.quietmmap, NULL,
- p->input->unit, NULL);
+ p->sky=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, p->input->ndim,
+ tl->numtiles, NULL, 0, cp->minmapsize,
+ p->cp.quietmmap, NULL, p->input->unit, NULL);
+ p->std=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, p->input->ndim,
+ tl->numtiles, NULL, 0, cp->minmapsize,
+ p->cp.quietmmap, NULL, p->input->unit, NULL);
/* Find the Sky and its STD on proper tiles. */
@@ -251,6 +251,14 @@ sky_and_std(struct noisechiselparams *p, char *checkname)
p->cpscorr = p->minstd>1 ? 1.0f : p->minstd;
+ /* Remove outlier tiles */
+ gal_tileinternal_no_outlier_local(p->sky, p->std, NULL, &p->cp.tl,
+ p->cp.interpmetric, p->cp.interpnumngb,
+ p->cp.numthreads, p->outliersclip,
+ p->outliersigma, checkname);
+
+
+
/* Interpolate and smooth the derived values. */
threshold_interp_smooth(p, &p->sky, &p->std, NULL, checkname);
diff --git a/bin/noisechisel/threshold.c b/bin/noisechisel/threshold.c
index 56159a4..4ef5137 100644
--- a/bin/noisechisel/threshold.c
+++ b/bin/noisechisel/threshold.c
@@ -33,6 +33,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gnuastro/threads.h>
#include <gnuastro/pointer.h>
#include <gnuastro/statistics.h>
+#include <gnuastro/permutation.h>
#include <gnuastro/interpolate.h>
#include <gnuastro-internal/timing.h>
@@ -460,9 +461,9 @@ qthresh_on_tile(void *in_prm)
tarray=tile->array; tblock=tile->block;
tile->array=gal_tile_block_relative_to_other(tile, p->conv);
tile->block=p->conv;
- usage->ndim=ndim; /* Since usage was modified in */
- usage->size=p->maxtcontig; /* place, it needs to be */
- gal_data_copy_to_allocated(tile, usage); /* re-initialized. */
+ usage->ndim=ndim; /* Since usage was modified in */
+ usage->size=p->maxtcontig; /* place, it needs to be */
+ gal_data_copy_to_allocated(tile, usage);/* re-initialized. */
tile->array=tarray; tile->block=tblock;
}
@@ -609,10 +610,12 @@ threshold_quantile_find_apply(struct noisechiselparams *p)
/* Allocate space for the quantile threshold values. */
qprm.erode_th=gal_data_alloc(NULL, p->input->type, p->input->ndim,
tl->numtiles, NULL, 0, cp->minmapsize,
- p->cp.quietmmap, NULL, p->input->unit, NULL);
+ p->cp.quietmmap, NULL, p->input->unit,
+ NULL);
qprm.noerode_th=gal_data_alloc(NULL, p->input->type, p->input->ndim,
tl->numtiles, NULL, 0, cp->minmapsize,
- p->cp.quietmmap, NULL, p->input->unit, NULL);
+ p->cp.quietmmap, NULL, p->input->unit,
+ NULL);
qprm.expand_th = ( p->detgrowquant!=1.0f
? gal_data_alloc(NULL, p->input->type, p->input->ndim,
tl->numtiles, NULL, 0, cp->minmapsize,
@@ -666,11 +669,12 @@ threshold_quantile_find_apply(struct noisechiselparams *p)
}
- /* Remove outliers if requested. */
- if(p->outliersigma!=0.0)
- gal_tileinternal_no_outlier(qprm.erode_th, qprm.noerode_th,
- qprm.expand_th, &p->cp.tl, p->outliersclip,
- p->outliersigma, p->qthreshname);
+ /* Remove the outliers. */
+ gal_tileinternal_no_outlier_local(qprm.erode_th, qprm.noerode_th,
+ qprm.expand_th, &cp->tl,
+ cp->interpmetric, cp->interpnumngb,
+ cp->numthreads, p->outliersclip,
+ p->outliersigma, p->qthreshname);
/* Check if the number of acceptable tiles is more than the minimum
diff --git a/bin/statistics/aststatistics.conf
b/bin/statistics/aststatistics.conf
index ee26c09..19807c2 100644
--- a/bin/statistics/aststatistics.conf
+++ b/bin/statistics/aststatistics.conf
@@ -24,7 +24,7 @@
# Sky and its STD settings
khdu 1
meanmedqdiff 0.005
- outliersigma 10
+ outliersigma 5
outliersclip 3,0.2
smoothwidth 3
sclipparams 3,0.1
diff --git a/bin/statistics/sky.c b/bin/statistics/sky.c
index df1c437..9f34fa7 100644
--- a/bin/statistics/sky.c
+++ b/bin/statistics/sky.c
@@ -212,10 +212,10 @@ sky(struct statisticsparams *p)
/* Remove outliers if requested. */
- if(p->outliersigma!=0.0)
- gal_tileinternal_no_outlier(p->sky_t, p->std_t, NULL, &p->cp.tl,
- p->outliersclip, p->outliersigma,
- p->checkskyname);
+ gal_tileinternal_no_outlier_local(p->sky_t, p->std_t, NULL, &cp->tl,
+ cp->interpmetric, cp->interpnumngb,
+ cp->numthreads, p->outliersclip,
+ p->outliersigma, p->checkskyname);
/* Interpolate the Sky and its standard deviation. */
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 6105e9f..9bae32b 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -493,6 +493,7 @@ Sky value
NoiseChisel
+* NoiseChisel changes after publication:: Updates since published papers.
* Invoking astnoisechisel:: Options and arguments for NoiseChisel.
Invoking NoiseChisel
@@ -13928,27 +13929,37 @@ However, prior to interpolating over the failed
tiles, another point should be c
In some cases, the gradient over these wings can be on scales that is larger
than the tiles.
The mean-median distance test will pass on such tiles and will cause a strong
peak in the interpolated tile grid, see @ref{Detecting large extended targets}.
-The tiles that exist over the wings of large galaxies or bright stars are
outliers in the distribution of tiles that passed the mean-median quantile
distance test.
+These tiles actually have signal in them and are outliers in the distribution
of tiles that passed the mean-median quantile distance test.
Therefore, the final step of ``quantifying signal in a tile'' is to look at
this distribution and remove the outliers.
@mymath{\sigma}-clipping is a good solution for removing a few outliers, but
the problem with outliers of this kind is that there may be many such tiles
(depending on the large/bright stars/galaxies in the image).
-Therefore a novel outlier rejection algorithm will be used.
+
+For each tile, we find the nearest @mymath{N_{ngb}} tiles that had a usable
value (@mymath{N_{ngb}} is the value given to @option{--interpnumngb}).
+We then sort them and find the difference between the largest and
second-to-smallest elements (The minimum isn't used because the scatter can be
large).
+Let's call this the tile's @emph{slope}.
+In effect it is a very crude estimate of the slope of the @mymath{N_{ngb}}
tiles in the vicinity of that tile.
+All the tiles that are on a region of flat noise will have similar slope
values, but if a few tiles fall on the wings of a bright star or galaxy, their
slope will be significantly larger than the flat region tiles.
+We just have to find the smallest tile slope value that is an outlier compared
to the rest and reject all tiles with a slope larger than that.
@cindex Outliers
@cindex Identifying outliers
-To identify the first outlier, we'll use the distribution of distances between
sorted elements.
-If there are @mymath{N} successful tiles, for every tile, the distance between
the adjacent @mymath{N/2} previous elements is found, giving a distribution
containing @mymath{N/2-1} points.
-The @mymath{\sigma}-clipped median and standard deviation of this distribution
is then found (@mymath{\sigma}-clipping is configured with
@option{--outliersclip}).
-Finally, if the distance between the element and its previous element is more
than @option{--outliersigma} multiples of the @mymath{\sigma}-clipped standard
deviation added with the @mymath{\sigma}-clipped median, that element is
considered an outlier and all tiles larger than that value are ignored.
+To identify the smallest outlier, we'll use the distribution of distances
between sorted elements.
+Let's assume there are @mymath{N} successful tiles in the image.
+We sort them in increasing order based on their @emph{slope} (defined above).
+So the array of tile slopes can be written as @mymath{@{s[0], s[1], ...,
s[N]@}}, where @mymath{s[0]<s[1]} and so on.
+We then start from element @mymath{M} and calculate the distances between all
two adjacent values before it: @mymath{@{s[1]-s[0], s[2]-s[1], s[M-1]-s[M-2]@}}.
+
+The @mymath{\sigma}-clipped median and standard deviation of this distribution
is then found (@mymath{\sigma}-clipping is configured with
@option{--outliersclip} which takes two values, see @ref{Sigma clipping}).
+If the distance between the element and its previous element is more than
@option{--outliersigma} multiples of the @mymath{\sigma}-clipped standard
deviation added with the @mymath{\sigma}-clipped median, that element is
considered an outlier and all tiles larger than that value are ignored.
Formally, if we assume there are @mymath{N} elements.
They are first sorted.
-Searching for the outlier starts on element @mymath{N/2} (integer division).
-Let's take @mymath{v_i} to be the @mymath{i}-th element of the sorted input
(with no blank values) and @mymath{m} and @mymath{\sigma} as the
@mymath{\sigma}-clipped median and standard deviation from the distances of the
previous @mymath{N/2-1} elements (not including @mymath{v_i}).
+Searching for the outlier starts on element @mymath{N/3} (integer division).
+Let's take @mymath{v_i} to be the @mymath{i}-th element of the sorted input
(with no blank values) and @mymath{m} and @mymath{\sigma} as the
@mymath{\sigma}-clipped median and standard deviation from the distances of the
previous @mymath{N/3-1} elements (not including @mymath{v_i}).
If the value given to @option{--outliersigma} is displayed with @mymath{s},
the @mymath{i}-th element is considered as an outlier when the condition below
is true.
@dispmath{{(v_i-v_{i-1})-m\over \sigma}>s}
-Since @mymath{i} begins from the median, the outlier has to be larger than the
median.
+Since @mymath{i} begins from the @mymath{N/3}-th element in the sorted array
(a quantile of @mymath{1/3=0.33}), the outlier has to be larger than the
@mymath{0.33} quantile value of the dataset.
You can use the check images (for example @option{--checksky} in the
Statistics program or @option{--checkqthresh}, @option{--checkdetsky} and
@option{--checksky} options in NoiseChisel for any of its steps that uses this
outlier rejection) to inspect the steps and see which tiles have been discarded
as outliers prior to interpolation.
@@ -14523,9 +14534,9 @@ The file will have two extensions for each step (one
for the Sky and one for the
@cindex Segmentation
Once instrumental signatures are removed from the raw data (image) in the
initial reduction process (see @ref{Data manipulation}).
You are naturally eager to start answering the scientific questions that
motivated the data collection in the first place.
-However, the raw dataset/image is just an array of values/pixels, that is all!
These raw values cannot directly be used to answer your scientific questions:
for example ``how many galaxies are there in the image?''.
+However, the raw dataset/image is just an array of values/pixels, that is all!
These raw values cannot directly be used to answer your scientific questions:
for example ``how many galaxies are there in the image?'', ``What is their
brightness?'' and etc.
-The first high-level step in the analysis of your dataset will thus be to
classify, or label, the dataset elements (pixels) into two classes:
+The first high-level step your analysis will therefore be to classify, or
label, the dataset elements (pixels) into two classes:
1) Noise, where random effects are the major contributor to the value, and
2) Signal, where non-random factors (for example light from a distant galaxy)
are present.
This classification of the elements in a dataset is formally known as
@emph{detection}.
@@ -14548,7 +14559,7 @@ $ astarithmetic in.fits 100 gt 2 connected-components
Since almost no astronomical target has such sharp edges, we need a more
advanced detection methodology.
NoiseChisel uses a new noise-based paradigm for detection of very extended and
diffuse targets that are drowned deeply in the ocean of noise.
It was initially introduced in @url{https://arxiv.org/abs/1505.01664, Akhlaghi
and Ichikawa [2015]} and improvements after the first four were published in
@url{https://arxiv.org/abs/1909.11230, Akhlaghi [2019]}.
-Please take the time to go through these papers to most effectively use
NoiseChisel.
+Please take the time to go through these papers to most effectively understand
the need of NoiseChisel and how best to use it.
The name of NoiseChisel is derived from the first thing it does after
thresholding the dataset: to erode it.
In mathematical morphology, erosion on pixels can be pictured as carving-off
boundary pixels.
@@ -14577,8 +14588,8 @@ In this case, the output won't be a binary image any
more, the signal will have
You can then directly feed NoiseChisel's output into MakeCatalog for
measurements over the detections and the production of a catalog (see
@ref{MakeCatalog}).
Thanks to the published papers mentioned above, there is no need to provide a
more complete introduction to NoiseChisel in this book.
-@c However, published papers cannot be updated any more, but the software has
evolved/changed.
-@c The changes since publication are documented in @ref{NoiseChisel changes
after publication}.
+However, published papers cannot be updated any more, but the software has
evolved/changed.
+The changes since publication are documented in @ref{NoiseChisel changes after
publication}.
In @ref{Invoking astnoisechisel}, the details of running NoiseChisel and its
options are discussed.
As discussed above, detection is one of the most important steps for your
scientific result.
@@ -14597,24 +14608,41 @@ Below, in @ref{Invoking astnoisechisel}, we will
review NoiseChisel's input, det
If you have used NoiseChisel within your research, please run it with
@option{--cite} to list the papers you should cite and how to acknowledge its
funding sources.
@menu
+* NoiseChisel changes after publication:: Updates since published papers.
* Invoking astnoisechisel:: Options and arguments for NoiseChisel.
@end menu
-@c @node NoiseChisel changes after publication, Invoking astnoisechisel,
NoiseChisel, NoiseChisel
-@c @subsection NoiseChisel changes after publication
+@node NoiseChisel changes after publication, Invoking astnoisechisel,
NoiseChisel, NoiseChisel
+@subsection NoiseChisel changes after publication
-@c NoiseChisel was initially introduced in
@url{https://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa [2015]} and
updates after the first four years were published in
@url{https://arxiv.org/abs/1909.11230, Akhlaghi [2019]}.
-@c It is thus strongly recommended to read this paper for a good understanding
of what it does and how each parameter influences the output.
-@c To help in understanding how it works, those papers have many figures
showing every step on multiple mock and real examples.
+NoiseChisel was initially introduced in @url{https://arxiv.org/abs/1505.01664,
Akhlaghi and Ichikawa [2015]} and updates after the first four years were
published in @url{https://arxiv.org/abs/1909.11230, Akhlaghi [2019]}.
+It is thus strongly recommended to read this paper for a good understanding of
what it does and how each parameter influences the output.
+To help in understanding how it works, those papers have many figures showing
every step on multiple mock and real examples.
-@c However, the papers cannot be updated anymore, but NoiseChisel has evolved
(and will continue to do so): better algorithms or steps have been found, thus
options will be added or removed.
-@c This book is thus the final and definitive guide to NoiseChisel.
-@c The aim of this section is to make the transition from the paper to the
installed version on your system, as smooth as possible with the list below.
-@c For a more detailed list of changes in previous Gnuastro releases/versions,
please see the @file{NEWS} file@footnote{The @file{NEWS} file is present in the
released Gnuastro tarball, see @ref{Release tarball}.}.
+However, the papers cannot be updated anymore, but NoiseChisel has evolved
(and will continue to do so): better algorithms or steps have been found, thus
options will be added or removed.
+This book is thus the final and definitive guide to NoiseChisel.
+The aim of this section is to make the transition from the papers above to the
installed version on your system, as smooth as possible with the list below.
+For a more detailed list of changes in previous Gnuastro releases/versions,
please see the @file{NEWS} file@footnote{The @file{NEWS} file is present in the
released Gnuastro tarball, see @ref{Release tarball}.}.
+@itemize
+@item
+Outlier rejection is now also applied on the tiles with good Sky value (see
below for more on the new outlier rejection algorithm).
+
+@item
+Improved outlier rejection for identifying tiles without any signal (improving
thresholds and Sky estimation):
+Until version 0.14, outliers were defined globally by finding the first
outlier-by-distance using the raw tile values in all the tiles within the image.
+But given its global nature, this method would fail to find the outliers in
the vicinity of bright sources easily.
+Therefore the faint wings of galaxies/stars could be mistakenly identified as
Sky and in the sky extension the footprint of larger objects became visible.
+
+It was possible to play with the parameters to correct this, but a better way
was found and implemented in version 0.14: instead of finding outliers in the
raw tile values, we now measure the @emph{slope} of the tile's nearby tiles and
find outlier in the slopes.
+For more on the outlier-by-distance algorithm and the definition of
@emph{slope}, see @ref{Quantifying signal in a tile}.
+In our tests, this gave a much improved estimate of the internal thresholds
and final Sky values.
+@end itemize
-@node Invoking astnoisechisel, , NoiseChisel, NoiseChisel
+
+
+@node Invoking astnoisechisel, , NoiseChisel changes after publication,
NoiseChisel
@subsection Invoking NoiseChisel
NoiseChisel will detect signal in noise producing a multi-extension dataset
containing a binary detection map which is the same size as the input.
@@ -25763,13 +25791,16 @@ function will return @code{GAL_BLANK_SIZE_T}.
@end deftypefun
@deftypefun {gal_data_t *} gal_statistics_quantile_function (gal_data_t
@code{*input}, gal_data_t @code{*value}, int @code{inplace})
-Return a single-element (@code{double} or @code{float64}) dataset
-containing the quantile function of the non-blank values in @code{input} at
-@code{value}. In other words, this function will return the quantile of
-@code{value} in @code{input}. If the value is smaller than the input's
-smallest element, the returned value will be zero. If the value is larger
-than the input's largest element, then the returned value is 1. See
-@code{gal_statistics_median} for a description of @code{inplace}.
+
+Return a single-element dataset containing the quantile function of the
non-blank values in @code{input} at @code{value} (a single-element dataset).
+The numerical data type is of the returned dataset is @code{float64} (or
@code{double}).
+In other words, this function will return the quantile of @code{value} in
@code{input}.
+@code{value} has to have the same type as @code{input}.
+See @code{gal_statistics_median} for a description of @code{inplace}.
+
+When all elements are blank, the returned value will be NaN.
+If the value is smaller than the input's smallest element, the returned value
will be negative infinity.
+If the value is larger than the input's largest element, then the returned
value will be positive infinity
@end deftypefun
@deftypefun {gal_data_t *} gal_statistics_unique (gal_data_t @code{*input},
int @code{inplace})
diff --git a/lib/gnuastro-internal/tile-internal.h
b/lib/gnuastro-internal/tile-internal.h
index 0465fc5..30dd6a3 100644
--- a/lib/gnuastro-internal/tile-internal.h
+++ b/lib/gnuastro-internal/tile-internal.h
@@ -66,6 +66,14 @@ gal_tileinternal_no_outlier(gal_data_t *first, gal_data_t
*second,
double *outliersclip, float outliersigma,
char *filename);
+void
+gal_tileinternal_no_outlier_local(gal_data_t *input, gal_data_t *second,
+ gal_data_t *third,
+ struct gal_tile_two_layer_params *tl,
+ uint8_t metric, size_t numneighbors,
+ size_t numthreads, double *outliersclip,
+ double outliersigma, char *filename);
+
__END_C_DECLS /* From C++ preparations */
#endif /* __GAL_TILE_INTERNAL_H__ */
diff --git a/lib/gnuastro/interpolate.h b/lib/gnuastro/interpolate.h
index 2b84d0d..21c7e13 100644
--- a/lib/gnuastro/interpolate.h
+++ b/lib/gnuastro/interpolate.h
@@ -88,10 +88,10 @@ enum gal_interpolate_1D_types
gal_data_t *
gal_interpolate_neighbors(gal_data_t *input,
- struct gal_tile_two_layer_params *tl,
- uint8_t metric, size_t numneighbors,
- size_t numthreads, int onlyblank,
- int aslinkedlist, int function);
+ struct gal_tile_two_layer_params *tl,
+ uint8_t metric, size_t numneighbors,
+ size_t numthreads, int onlyblank,
+ int aslinkedlist, int function);
gsl_spline *
gal_interpolate_1d_make_gsl_spline(gal_data_t *X, gal_data_t *Y, int type_1d);
diff --git a/lib/interpolate.c b/lib/interpolate.c
index 545f817..100a753 100644
--- a/lib/interpolate.c
+++ b/lib/interpolate.c
@@ -48,9 +48,9 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
/*************** (Dimension agnostic) ****************/
/*********************************************************************/
/* These are bit-flags, so we're using hexadecimal notation. */
-#define INTERPOLATE_FLAGS_NO 0
-#define INTERPOLATE_FLAGS_CHECKED 0x1
-#define INTERPOLATE_FLAGS_BLANK 0x2
+#define INTERPOLATE_FLAGS_NO 0
+#define INTERPOLATE_FLAGS_NGB_CHECKED 0x1
+#define INTERPOLATE_FLAGS_BLANK 0x2
@@ -96,8 +96,8 @@ interpolate_neighbors_on_thread(void *in_prm)
float dist, pdist;
uint8_t *b, *bf, *bb;
gal_list_void_t *tvll;
+ size_t ngb_counter, pind;
gal_list_dosizet_t *lQ, *sQ;
- size_t ngb_counter, pind, *dinc;
size_t i, index, fullind, chstart=0, ndim=input->ndim;
gal_data_t *tin, *tout, *tnear, *value=NULL, *nearest=NULL;
size_t size = (correct_index ? tl->tottilesinch : input->size);
@@ -106,13 +106,21 @@ interpolate_neighbors_on_thread(void *in_prm)
"icoord");
size_t *ncoord=gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__,
"ncoord");
- uint8_t *fullflag=&prm->thread_flags[tprm->id*input->size], *flag=fullflag;
+ uint8_t *flag, *fullflag=&prm->thread_flags[tprm->id*input->size];
+ /* Based on the above. */
+ size_t *dinc=gal_dimension_increment(ndim, dsize);
- /* Initializations. */
+
+ /* Initialize the flags array. We need two flags during this processing:
+ 1) to see if there are blanks. 2) to see if a neighbor has been
+ checked. These are both binary (0 or 1). So to avoid wasting space, we
+ will use bits to store them. We start with only setting the blank flag
+ once for the whole thread. Then for each interpolated pixel, we reset
+ the neighbor-check flag. */
+ flag=fullflag;
bb=prm->blanks->array;
bf=(b=fullflag)+input->size;
- dinc=gal_dimension_increment(ndim, dsize);
do *b = *bb++ ? INTERPOLATE_FLAGS_BLANK : 0; while(++b<bf);
@@ -123,8 +131,9 @@ interpolate_neighbors_on_thread(void *in_prm)
{
nv=gal_pointer_increment(tvll->v, tprm->id*prm->numneighbors,
input->type);
- gal_list_data_add_alloc(&nearest, nv, tin->type, 1, &prm->numneighbors,
- NULL, 0, -1, 1, NULL, NULL, NULL);
+ gal_list_data_add_alloc(&nearest, nv, tin->type, 1,
+ &prm->numneighbors, NULL, 0, -1, 1,
+ NULL, NULL, NULL);
tin=tin->next;
}
gal_list_data_reverse(&nearest);
@@ -183,7 +192,7 @@ interpolate_neighbors_on_thread(void *in_prm)
/* Reset all checked bits in the flags array to 0. */
ngb_counter=0;
bf=(b=flag)+size;
- do *b &= ~(INTERPOLATE_FLAGS_CHECKED); while(++b<bf);
+ do *b &= ~(INTERPOLATE_FLAGS_NGB_CHECKED); while(++b<bf);
/* Get the coordinates of this pixel (to be interpolated). */
@@ -233,7 +242,7 @@ interpolate_neighbors_on_thread(void *in_prm)
IMPORTANT: we must not check for blank values here,
otherwise we won't be able to parse over extended blank
regions. */
- if( !(flag[nind] & INTERPOLATE_FLAGS_CHECKED) )
+ if( !(flag[nind] & INTERPOLATE_FLAGS_NGB_CHECKED) )
{
/* Get the coordinates of this neighbor. */
gal_dimension_index_to_coord(nind, ndim, dsize, ncoord);
@@ -245,7 +254,7 @@ interpolate_neighbors_on_thread(void *in_prm)
gal_list_dosizet_add(&lQ, &sQ, nind, dist);
/* Flag this neighbor as checked. */
- flag[nind] |= INTERPOLATE_FLAGS_CHECKED;
+ flag[nind] |= INTERPOLATE_FLAGS_NGB_CHECKED;
}
} );
@@ -256,7 +265,8 @@ interpolate_neighbors_on_thread(void *in_prm)
if(sQ==NULL)
error(EXIT_FAILURE, 0, "%s: only %zu neighbors found while "
"you had asked to use %zu neighbors for close neighbor "
- "interpolation", __func__, ngb_counter, prm->numneighbors);
+ "interpolation", __func__, ngb_counter,
+ prm->numneighbors);
}
/* Calculate the desired statistic, and write it in the output. */
@@ -348,10 +358,10 @@ interpolate_copy_input(gal_data_t *input, int
aslinkedlist)
'input[i]' corresponds to 'tiles[i]' in the tessellation. */
gal_data_t *
gal_interpolate_neighbors(gal_data_t *input,
- struct gal_tile_two_layer_params *tl,
- uint8_t metric, size_t numneighbors,
- size_t numthreads, int onlyblank,
- int aslinkedlist, int function)
+ struct gal_tile_two_layer_params *tl,
+ uint8_t metric, size_t numneighbors,
+ size_t numthreads, int onlyblank,
+ int aslinkedlist, int function)
{
gal_data_t *tin, *tout;
struct interpolate_ngb_params prm;
@@ -361,9 +371,9 @@ gal_interpolate_neighbors(gal_data_t *input,
/* If there are no blank values in the array, AND we should only fill
blank values, then simply copy the input and abort. */
- if( (input->flag | GAL_DATA_FLAG_BLANK_CH) /* Zero bit is meaningful.*/
- && !(input->flag | GAL_DATA_FLAG_HASBLANK) /* There are no blanks. */
- && onlyblank ) /* Only interpolate blank.*/
+ if( (input->flag | GAL_DATA_FLAG_BLANK_CH) /* Zero bit is meaningful.*/
+ && !(input->flag | GAL_DATA_FLAG_HASBLANK)/* There are no blanks. */
+ && onlyblank ) /* Only interpolate blank.*/
return interpolate_copy_input(input, aslinkedlist);
diff --git a/lib/statistics.c b/lib/statistics.c
index bafab53..2f55605 100644
--- a/lib/statistics.c
+++ b/lib/statistics.c
@@ -430,7 +430,7 @@ gal_statistics_quantile(gal_data_t *input, double quantile,
int inplace)
r=a++; \
\
/* Increasing array: */ \
- if( *a < *(a+1) ) \
+ if( nbs->flag & GAL_DATA_FLAG_SORTED_I ) \
{ \
if( v>=*r ) \
{ \
@@ -455,17 +455,21 @@ gal_statistics_quantile(gal_data_t *input, double
quantile, int inplace)
if(parsed && a<af) index = a-r; \
}
size_t
-gal_statistics_quantile_function_index(gal_data_t *input, gal_data_t *value,
- int inplace)
+gal_statistics_quantile_function_index(gal_data_t *input,
+ gal_data_t *invalue, int inplace)
{
int parsed=0;
+ gal_data_t *value;
size_t index=GAL_BLANK_SIZE_T;
gal_data_t *nbs=gal_statistics_no_blank_sorted(input, inplace);
- /* A sanity check. */
- if(nbs->type!=value->type)
- error(EXIT_FAILURE, 0, "%s: the types of the input dataset and requested "
- "value have to be the same", __func__);
+ /* Make sure the value has the same type. */
+ if(invalue->size>1)
+ error(EXIT_FAILURE, 0, "%s: the 'value' argument must only have "
+ "one element", __func__);
+ value = ( (nbs->type==invalue->type)
+ ? invalue
+ : gal_data_copy_to_new_type(invalue, nbs->type) );
/* Only continue processing if we have non-blank elements. */
if(nbs->size)
@@ -494,6 +498,7 @@ gal_statistics_quantile_function_index(gal_data_t *input,
gal_data_t *value,
}
/* Clean up and return. */
+ if(value!=invalue) gal_data_free(value);
if(nbs!=input) gal_data_free(nbs);
return index;
}
@@ -519,11 +524,19 @@ gal_statistics_quantile_function(gal_data_t *input,
gal_data_t *value,
int inplace)
{
double *d;
- size_t dsize=1;
+ size_t ind, dsize=1;
gal_data_t *nbs=gal_statistics_no_blank_sorted(input, inplace);
gal_data_t *out=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &dsize,
NULL, 1, -1, 1, NULL, NULL, NULL);
- size_t ind=gal_statistics_quantile_function_index(input, value, inplace);
+
+ /* Sanity checks. */
+ if(value->size>1)
+ error(EXIT_FAILURE, 0, "%s: the 'value' argument must only have "
+ "one element", __func__);
+
+ /* Calculate the index of the value. */
+ ind=gal_statistics_quantile_function_index(input, value, inplace);
+ //printf("ind: %zu (%zu)\n", ind, input->size);
/* Only continue processing if there are non-blank values. */
if(nbs->size)
@@ -2396,7 +2409,7 @@ gal_statistics_sigma_clip(gal_data_t *input, float
multip, float param,
(float)(arr[i]-arr[i-1]), sarr[1], sarr[3], \
(((double)(arr[i]-arr[i-1])) - sarr[1])/sarr[3]); \
\
- /* Terminate the loop if the dist. is larger than requested. */ \
+ /* Terminate the loop if the dist is larger than requested. */ \
/* This shows we have reached the first outlier's position. */ \
if( (((double)(arr[i]-arr[i-1])) - sarr[1]) > sigma*sarr[3] ) \
{ \
diff --git a/lib/tile-internal.c b/lib/tile-internal.c
index 7965aa8..82a9161 100644
--- a/lib/tile-internal.c
+++ b/lib/tile-internal.c
@@ -26,11 +26,15 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <stdio.h>
#include <errno.h>
#include <error.h>
+#include <string.h>
#include <stdlib.h>
#include <gnuastro/tile.h>
+#include <gnuastro/threads.h>
#include <gnuastro/pointer.h>
#include <gnuastro/statistics.h>
+#include <gnuastro/interpolate.h>
+#include <gnuastro/permutation.h>
#include <gnuastro-internal/tile-internal.h>
@@ -71,7 +75,8 @@ tileinternal_no_outlier_work(gal_data_t *first, gal_data_t
*second,
second->array=gal_pointer_increment(second->array, start,
second->type);
if(third)
- third->array=gal_pointer_increment(third->array, start, third->type);
+ third->array=gal_pointer_increment(third->array, start,
+ third->type);
/* Correct their sizes. */
first->size=tottilesinch;
@@ -83,12 +88,12 @@ tileinternal_no_outlier_work(gal_data_t *first, gal_data_t
*second,
first array. */
arr1=first->array;
nbs=gal_statistics_no_blank_sorted(first, 0);
- outlier_p=gal_statistics_outlier_bydistance(1, nbs, nbs->size/2,
outliersigma,
- outliersclip[0], outliersclip[1],
- 1, 1);
- outlier_n=gal_statistics_outlier_bydistance(0, nbs, nbs->size/2,
outliersigma,
- outliersclip[0], outliersclip[1],
- 1, 1);
+ outlier_p=gal_statistics_outlier_bydistance(1, nbs, nbs->size/2,
+ outliersigma, outliersclip[0],
+ outliersclip[1], 1, 1);
+ outlier_n=gal_statistics_outlier_bydistance(0, nbs, nbs->size/2,
+ outliersigma, outliersclip[0],
+ outliersclip[1], 1, 1);
/* For a check.
{
float *med;
@@ -142,12 +147,12 @@ tileinternal_no_outlier_work(gal_data_t *first,
gal_data_t *second,
of them. */
arr2=second->array;
nbs=gal_statistics_no_blank_sorted(second, 0);
- outlier_p=gal_statistics_outlier_bydistance(1, nbs, nbs->size, outliersigma,
- outliersclip[0], outliersclip[1],
- 1, 1);
- outlier_n=gal_statistics_outlier_bydistance(0, nbs, nbs->size, outliersigma,
- outliersclip[0], outliersclip[1],
- 1, 1);
+ outlier_p=gal_statistics_outlier_bydistance(1, nbs, nbs->size,
+ outliersigma, outliersclip[0],
+ outliersclip[1], 1, 1);
+ outlier_n=gal_statistics_outlier_bydistance(0, nbs, nbs->size,
+ outliersigma, outliersclip[0],
+ outliersclip[1], 1, 1);
gal_data_free(nbs);
if(outlier_p)
{
@@ -272,3 +277,453 @@ gal_tileinternal_no_outlier(gal_data_t *first, gal_data_t
*second,
}
}
}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/*************************************************************/
+/************ Local outlier removal ************/
+/*************************************************************/
+#define TILEINTERNAL_OUTLIER_FLAGS_NO 0
+#define TILEINTERNAL_OUTLIER_FLAGS_NGB_CHECKED 0x1
+#define TILEINTERNAL_OUTLIER_FLAGS_BLANK 0x2
+
+struct tileinternal_outlier_local
+{
+ gal_data_t *input;
+ gal_data_t *measure;
+ gal_data_t *blanks;
+ size_t numneighbors;
+ uint8_t *thread_flags;
+ gal_list_void_t *ngb_vals;
+ float (*metric)(size_t *, size_t *, size_t );
+
+ struct gal_tile_two_layer_params *tl;
+};
+
+
+
+
+
+/* Run the interpolation on many threads. */
+static void *
+gal_tileinternal_no_outlier_local_on_thread(void *in_prm)
+{
+ /* Low-level variables that others depend on. */
+ struct gal_threads_params *tprm=(struct gal_threads_params *)in_prm;
+ struct tileinternal_outlier_local *prm=
+ (struct tileinternal_outlier_local *)(tprm->params);
+
+ /* Higher-level variables. */
+ struct gal_tile_two_layer_params *tl=prm->tl;
+ int correct_index=(tl && tl->totchannels>1 && !tl->workoverch);
+ gal_data_t *input=prm->input;
+
+ /* Rest of variables. */
+ void *nv;
+ uint8_t *b, *bf, *bb;
+ gal_list_void_t *tvll;
+ size_t ngb_counter, pind;
+ gal_list_dosizet_t *lQ, *sQ;
+ gal_data_t *tin, *tnear, *nearest=NULL;
+ float dist, pdist, *tnarr, *marr=prm->measure->array;
+ size_t i, index, fullind, chstart=0, ndim=input->ndim;
+ size_t size = (correct_index ? tl->tottilesinch : input->size);
+ size_t *dsize = (correct_index ? tl->numtilesinch : input->dsize);
+ size_t *icoord=gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__,
+ "icoord");
+ size_t *ncoord=gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__,
+ "ncoord");
+ uint8_t *flag, *fullflag=&prm->thread_flags[tprm->id*input->size];
+
+ /* Based on the above. */
+ size_t *dinc=gal_dimension_increment(ndim, dsize);
+
+ /* Initialize the flags array. We need two flags during this processing:
+ 1) to see if there are blanks. 2) to see if a neighbor has been
+ checked. These are both binary (0 or 1). So to avoid wasting space, we
+ will use bits to store them. We start with only setting the blank flag
+ once for the whole thread. Then for each interpolated pixel, we reset
+ the neighbor-check flag. */
+ flag=fullflag;
+ bb=prm->blanks->array;
+ bf=(b=fullflag)+input->size;
+ do *b = *bb++ ? TILEINTERNAL_OUTLIER_FLAGS_BLANK : 0; while(++b<bf);
+
+
+ /* 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_pointer_increment(tvll->v, tprm->id*prm->numneighbors,
+ input->type);
+ gal_list_data_add_alloc(&nearest, nv, tin->type, 1,
+ &prm->numneighbors, NULL, 0, -1, 1,
+ NULL, NULL, NULL);
+ tin=tin->next;
+ }
+ gal_list_data_reverse(&nearest);
+
+
+ /* Go over all the points given to this thread. */
+ for(i=0; tprm->indexs[i] != GAL_BLANK_SIZE_T; ++i)
+ {
+ /* For easy reading. */
+ fullind=tprm->indexs[i];
+
+
+ /* If we are on a blank element, then ignore this pixel. */
+ if( (fullflag[fullind] & TILEINTERNAL_OUTLIER_FLAGS_BLANK) )
+ { marr[fullind]=NAN; continue; }
+
+
+ /* Correct the index (if necessary). When the values come from a
+ tiled dataset, the caller might want to interpolate the values of
+ each channel separately (not mix values from different
+ channels). In such a case, the tiles of each channel (and their
+ values in 'input' are contiguous. So we need to correct
+ 'tprm->indexs[i]' (which is the index over the whole tessellation,
+ including all channels). */
+ if(correct_index)
+ {
+ /* Index of this tile in its channel. */
+ index = fullind % tl->tottilesinch;
+
+ /* Index of the first tile in this channel. */
+ chstart = (fullind / tl->tottilesinch) * tl->tottilesinch;
+
+ /* Set the channel's starting pointer for the flags. */
+ flag = gal_pointer_increment(fullflag, chstart, GAL_TYPE_UINT8);
+ }
+ else
+ {
+ chstart=0;
+ index=fullind;
+ }
+
+
+ /* Reset all checked bits in the flags array to 0. */
+ ngb_counter=0;
+ bf=(b=flag)+size;
+ do *b &= ~(TILEINTERNAL_OUTLIER_FLAGS_NGB_CHECKED); while(++b<bf);
+
+
+ /* Get the coordinates of this pixel (to be interpolated). */
+ gal_dimension_index_to_coord(index, ndim, dsize, icoord);
+
+
+ /* Start parsing the neighbors. We will use a two-way ordered linked
+ list structure. To start from the nearest and go out to the
+ farthest. */
+ lQ=sQ=NULL;
+ gal_list_dosizet_add(&lQ, &sQ, index, 0.0f);
+ while(sQ)
+ {
+ /* Pop-out (p) an index from the queue: */
+ pind=gal_list_dosizet_pop_smallest(&lQ, &sQ, &pdist);
+
+ /* 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] & TILEINTERNAL_OUTLIER_FLAGS_BLANK) )
+ {
+ tin=input;
+ for(tnear=nearest; tnear!=NULL; tnear=tnear->next)
+ {
+ memcpy(gal_pointer_increment(tnear->array, ngb_counter,
+ tin->type),
+ gal_pointer_increment(tin->array, chstart+pind,
+ tin->type),
+ gal_type_sizeof(tin->type));
+ tin=tin->next;
+ }
+
+ /* If we have filled all the elements clean up the linked
+ list and break out. */
+ if(++ngb_counter>=prm->numneighbors)
+ {
+ if(lQ) gal_list_dosizet_free(lQ);
+ break;
+ }
+ }
+
+ /* Go over all the neighbors of this popped pixel and add them to
+ the list of neighbors to be checked. */
+ GAL_DIMENSION_NEIGHBOR_OP(pind, ndim, dsize, 1, dinc,
+ {
+ /* Only look at neighbors that have not been checked. VERY
+ IMPORTANT: we must not check for blank values here,
+ otherwise we won't be able to parse over extended blank
+ regions. */
+ if( !(flag[nind] & TILEINTERNAL_OUTLIER_FLAGS_NGB_CHECKED) )
+ {
+ /* Get the coordinates of this neighbor. */
+ gal_dimension_index_to_coord(nind, ndim, dsize, ncoord);
+
+ /* Distance of this neighbor to the one to be filled. */
+ dist=prm->metric(icoord, ncoord, ndim);
+
+ /* Add this neighbor to the list. */
+ gal_list_dosizet_add(&lQ, &sQ, nind, dist);
+
+ /* Flag this neighbor as checked. */
+ flag[nind] |= TILEINTERNAL_OUTLIER_FLAGS_NGB_CHECKED;
+ }
+ } );
+
+ /* If there are no more meshes to add to the queue, then this
+ shows, there were not enough points for
+ interpolation. Normally, this loop should only be exited
+ through the 'currentnum>=numnearest' check above. */
+ if(sQ==NULL)
+ error(EXIT_FAILURE, 0, "%s: only %zu neighbors found while "
+ "you had asked to use %zu neighbors for close neighbor "
+ "interpolation", __func__, ngb_counter,
+ prm->numneighbors);
+ }
+
+ /* Calculate the desired statistic, and write it in the output. */
+ for(tnear=nearest; tnear!=NULL; tnear=tnear->next)
+ {
+ /* First, reset the sorting flags (which remain from the last
+ time). */
+ tnear->flag &= ~(GAL_DATA_FLAG_SORT_CH | GAL_DATA_FLAG_BLANK_CH);
+
+ /* For a check on the values.
+ { size_t i; float *f=tnear->array;
+ for(i=0;i<tnear->size;++i) printf("%f\n", f[i]); } */
+
+ /* Sort the elements, then find the difference between the
+ maximium and the value that is just after the minimum. We are
+ doing this because the scatter in the minimum can be large. */
+ tnarr=tnear->array;
+ gal_statistics_sort_increasing(tnear);
+ marr[fullind] = tnarr[tnear->size-1]-tnarr[1];
+ }
+ }
+
+
+ /* Clean up. */
+ for(tnear=nearest; tnear!=NULL; tnear=tnear->next) tnear->array=NULL;
+ gal_list_data_free(nearest);
+ free(icoord);
+ free(ncoord);
+ free(dinc);
+
+
+ /* Wait for all the other threads to finish and return. */
+ if(tprm->b) pthread_barrier_wait(tprm->b);
+ return NULL;
+}
+
+
+
+
+
+void
+gal_tileinternal_no_outlier_local(gal_data_t *input, gal_data_t *second,
+ gal_data_t *third,
+ struct gal_tile_two_layer_params *tl,
+ uint8_t metric, size_t numneighbors,
+ size_t numthreads, double *outliersclip,
+ double outliersigma, char *filename)
+{
+ gal_data_t *othresh;
+ float *base, *f, *ff, thresh;
+ struct tileinternal_outlier_local prm;
+ size_t owindow, ngbvnum=numthreads*numneighbors;
+ int permute=(tl && tl->totchannels>1 && tl->workoverch);
+
+
+ /* Sanity checks. */
+ if(numneighbors<=3)
+ error(EXIT_FAILURE, 0, "interpnumngb has to be larger than 3, but "
+ "is currently %zu", numneighbors);
+ if(input->type!=GAL_TYPE_FLOAT32)
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix "
+ "the problem. The input to this function (not NoiseChisel) "
+ "should be in 32-bit floating point, but it is %s", __func__,
+ PACKAGE_BUGREPORT, gal_type_name(input->type, 1));
+ if(second && second->type!=GAL_TYPE_FLOAT32)
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix "
+ "the problem. The 'second' argument to this function (not "
+ "NoiseChisel) should be in 32-bit floating point, but it is "
+ "%s", __func__, PACKAGE_BUGREPORT, gal_type_name(input->type, 1));
+ if(third && third->type!=GAL_TYPE_FLOAT32)
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix "
+ "the problem. The 'third' argument to this function (not "
+ "NoiseChisel) should be in 32-bit floating point, but it is "
+ "%s", __func__, PACKAGE_BUGREPORT, gal_type_name(input->type, 1));
+ if(second && gal_dimension_is_different(input, second) )
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix "
+ "the problem. The 'second' argument to this function (not "
+ "NoiseChisel) doesn't have the same size as the input",
+ __func__, PACKAGE_BUGREPORT);
+ if(third && gal_dimension_is_different(input, third) )
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix "
+ "the problem. The 'third' argument to this function (not "
+ "NoiseChisel) doesn't have the same size as the input",
+ __func__, PACKAGE_BUGREPORT);
+
+
+ /* Initialize the constant parameters. */
+ prm.tl = tl;
+ prm.ngb_vals = NULL;
+ prm.input = input;
+ prm.numneighbors = numneighbors;
+
+
+ /* Set the distance metric. */
+ switch(metric)
+ {
+ case GAL_INTERPOLATE_NEIGHBORS_METRIC_RADIAL:
+ prm.metric=gal_dimension_dist_radial;
+ break;
+ case GAL_INTERPOLATE_NEIGHBORS_METRIC_MANHATTAN:
+ prm.metric=gal_dimension_dist_manhattan;
+ break;
+ default:
+ error(EXIT_FAILURE, 0, "%s: %d is not a valid metric identifier",
+ __func__, metric);
+ }
+
+
+ /* Flag the blank values. */
+ prm.blanks=gal_blank_flag(input);
+
+
+ /* If the input is from a tile structure and the user has asked to ignore
+ channels, then re-order the values. */
+ if(permute)
+ {
+ /* Prepare the permutation (if necessary/not already defined). */
+ gal_tile_full_permutation(tl);
+
+ /* Re-order values to ignore channels (if necessary). */
+ gal_permutation_apply(input, tl->permutation);
+ gal_permutation_apply(prm.blanks, tl->permutation);
+ }
+
+
+ /* Necessary allocations and basic checks. 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.measure=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, input->ndim,
+ input->dsize, input->wcs, 0, input->minmapsize,
+ input->quietmmap, NULL, input->unit, NULL);
+ gal_list_void_add(&prm.ngb_vals,
+ gal_pointer_allocate(input->type, ngbvnum, 0,
+ __func__, "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_pointer_allocate(GAL_TYPE_UINT8,
+ numthreads*input->size, 0,
+ __func__, "prm.thread_flags");
+
+
+ /* Spin off the threads. */
+ gal_threads_spin_off(gal_tileinternal_no_outlier_local_on_thread,
+ &prm, input->size, numthreads, input->minmapsize,
+ input->quietmmap);
+
+
+ /* Find the outliers in the distribution, we will start from the first
+ third of the cases to find the first outlier. Note that this should
+ not be done in-place because we need the 'measure' arrray
+ afterwards. */
+ owindow=(prm.measure->size - gal_blank_number(prm.measure, 1))/3;
+ othresh=gal_statistics_outlier_bydistance(1, prm.measure, owindow,
+ outliersigma, outliersclip[0],
+ outliersclip[1], 0, 1);
+
+
+ /* If an outlier threshold was actually found, then mask all the tiles
+ larger than that value. */
+ if(othresh)
+ {
+ base=prm.measure->array;
+ ff=(f=input->array)+input->size;
+ thresh=((float *)(othresh->array))[0];
+ do { *f = isnan(*f) ? *f : (*base>thresh ? NAN : *f); ++base; }
+ while(++f<ff);
+ }
+
+
+ /* For a check.
+ printf("measure-threshold: %f\n", thresh);
+ if(permute)
+ gal_permutation_apply_inverse(prm.measure, tl->permutation);
+ gal_tile_full_values_write(prm.measure, tl, 1, "measure.fits",
+ NULL, NULL);
+ */
+
+
+ /* If the values were permuted for the interpolation, then re-order the
+ values back to their original location (so they correspond to their
+ tile indexs. */
+ if(permute)
+ gal_permutation_apply_inverse(input, tl->permutation);
+
+
+ /* If the 'other' arrays are given, set all the blank elements here to
+ blank there too. */
+ if(second)
+ {
+ base=input->array; ff=(f=second->array)+second->size;
+ do { *f = isnan(*base++) ? NAN : *f ;} while(++f<ff);
+ }
+ if(third)
+ {
+ base=input->array; ff=(f=third->array)+third->size;
+ do { *f = isnan(*base++) ? NAN : *f ;} while(++f<ff);
+ }
+
+
+ /* Write the check images if necessary. */
+ if(filename)
+ {
+ input->name="VALUE1_NO_OUTLIER";
+ gal_tile_full_values_write(input, tl, 1, filename, NULL, NULL);
+ input->name=NULL;
+ if(second)
+ {
+ second->name="VALUE2_NO_OUTLIER";
+ gal_tile_full_values_write(second, tl, 1, filename, NULL, NULL);
+ second->name=NULL;
+ }
+ if(third)
+ {
+ third->name="VALUE3_NO_OUTLIER";
+ gal_tile_full_values_write(third, tl, 1, filename, NULL, NULL);
+ third->name=NULL;
+ }
+ }
+
+
+ /* Clean up and return. */
+ gal_data_free(othresh);
+ free(prm.thread_flags);
+ gal_data_free(prm.blanks);
+ gal_data_free(prm.measure);
+ gal_list_void_free(prm.ngb_vals, 1);
+}
diff --git a/tests/during-dev.sh b/tests/during-dev.sh
index 65b9de9..abfc22d 100755
--- a/tests/during-dev.sh
+++ b/tests/during-dev.sh
@@ -88,6 +88,8 @@ options=
+
+
# RUN THE PROCEDURES
# ==================
diff --git a/tests/noisechisel/noisechisel.sh b/tests/noisechisel/noisechisel.sh
index 383983b..13c5db4 100755
--- a/tests/noisechisel/noisechisel.sh
+++ b/tests/noisechisel/noisechisel.sh
@@ -54,6 +54,6 @@ if [ ! -f $img ]; then echo "$img does not exist.";
exit 77; fi
# 'check_with_program' can be something like Valgrind or an empty
# string. Such programs will execute the command if present and help in
# debugging when the developer doesn't have access to the user's system.
-$check_with_program $execname $img --tilesize=100,100 --snquant=0.999 \
+$check_with_program $execname $img --detgrowquant=0.7 \
--cleangrowndet --checkdetection \
--continueaftercheck
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [gnuastro-commits] master 4a08b3f: NoiseChisel & Statistics: improved outlier rejection in Sky,
Mohammad Akhlaghi <=