From c2478a598e2346f44cc4b4dc7a8fbf008c0885dd Mon Sep 17 00:00:00 2001
From: Andre Offringa <offringa@astron.nl>
Date: Thu, 17 Aug 2023 10:49:42 +0000
Subject: [PATCH] Add absolute auto-masking threshold
---
cpp/radler.cc | 24 +++++++++++++++---------
cpp/settings.h | 8 +++++++-
cpp/test/test_python_deconvolution.cc | 2 +-
cpp/test/test_radler.cc | 15 ++++++++-------
python/docstrings/settings_docstrings.h | 18 +++++++++++-------
python/pysettings.cc | 8 ++++++--
python/test/test_radler.py | 4 ++--
python/test/test_settings.py | 4 ++--
python/test/test_work_table.py | 4 ++--
9 files changed, 54 insertions(+), 33 deletions(-)
diff --git a/cpp/radler.cc b/cpp/radler.cc
index 78ee214a..7510eb3f 100644
--- a/cpp/radler.cc
+++ b/cpp/radler.cc
@@ -154,7 +154,9 @@ void Radler::Perform(bool& reached_major_threshold,
double stddev = integrated.StdDevFromMAD();
Logger::Info << "Estimated standard deviation of background noise: "
<< FluxDensity::ToNiceString(stddev) << '\n';
- if (settings_.auto_mask_sigma && auto_mask_is_finished_) {
+ const bool auto_mask_is_enabled =
+ settings_.auto_mask_sigma || settings_.absolute_auto_mask_threshold;
+ if (auto_mask_is_enabled && auto_mask_is_finished_) {
// When we are in the second phase of automasking, don't use
// the RMS background anymore
parallel_deconvolution_->SetRmsFactorImage(Image());
@@ -208,12 +210,16 @@ void Radler::Perform(bool& reached_major_threshold,
parallel_deconvolution_->SetRmsFactorImage(std::move(rms_image));
}
}
- if (settings_.auto_mask_sigma && !auto_mask_is_finished_) {
+ if (auto_mask_is_enabled && !auto_mask_is_finished_) {
+ const double combined_auto_mask_threshold =
+ std::max(stddev * settings_.auto_mask_sigma.value_or(0.0),
+ settings_.absolute_auto_mask_threshold.value_or(0.0));
parallel_deconvolution_->SetThreshold(
- std::max(stddev * (*settings_.auto_mask_sigma), settings_.threshold));
+ std::max(combined_auto_mask_threshold, settings_.absolute_threshold));
} else if (settings_.auto_threshold_sigma) {
- parallel_deconvolution_->SetThreshold(std::max(
- stddev * (*settings_.auto_threshold_sigma), settings_.threshold));
+ parallel_deconvolution_->SetThreshold(
+ std::max(stddev * (*settings_.auto_threshold_sigma),
+ settings_.absolute_threshold));
}
integrated.Reset();
@@ -222,7 +228,7 @@ void Radler::Perform(bool& reached_major_threshold,
residual_set.LoadAndAveragePsfs();
if (settings_.algorithm_type == AlgorithmType::kMultiscale) {
- if (settings_.auto_mask_sigma) {
+ if (auto_mask_is_enabled) {
if (auto_mask_is_finished_) {
parallel_deconvolution_->SetAutoMaskMode(false, true);
} else {
@@ -230,7 +236,7 @@ void Radler::Perform(bool& reached_major_threshold,
}
}
} else {
- if (settings_.auto_mask_sigma && auto_mask_is_finished_) {
+ if (auto_mask_is_enabled && auto_mask_is_finished_) {
if (auto_mask_.empty()) {
// Generate the auto-mask from the model image(s)
auto_mask_.assign(image_width_ * image_height_, false);
@@ -251,7 +257,7 @@ void Radler::Perform(bool& reached_major_threshold,
residual_set, model_set, psf_images, table_->PsfOffsets(),
reached_major_threshold);
- if (!reached_major_threshold && settings_.auto_mask_sigma &&
+ if (!reached_major_threshold && auto_mask_is_enabled &&
!auto_mask_is_finished_) {
Logger::Info << "Auto-masking threshold reached; continuing next major "
"iteration with deeper threshold and mask.\n";
@@ -339,7 +345,7 @@ void Radler::InitializeDeconvolutionAlgorithm(
}
algorithm->SetMaxIterations(settings_.minor_iteration_count);
- algorithm->SetThreshold(settings_.threshold);
+ algorithm->SetThreshold(settings_.absolute_threshold);
algorithm->SetMinorLoopGain(settings_.minor_loop_gain);
algorithm->SetMajorLoopGain(settings_.major_loop_gain);
algorithm->SetCleanBorderRatio(settings_.border_ratio);
diff --git a/cpp/settings.h b/cpp/settings.h
index 7ddf08a0..1cc6e136 100644
--- a/cpp/settings.h
+++ b/cpp/settings.h
@@ -150,7 +150,7 @@ struct Settings {
* which means that Radler will keep continuing until another criterion (e.g.
* nr. of iterations) is reached.
*/
- double threshold = 0.0;
+ double absolute_threshold = 0.0;
/**
* Gain value for minor loop iterations.
@@ -198,6 +198,12 @@ struct Settings {
*/
std::optional<double> auto_mask_sigma = std::nullopt;
+ /**
+ * @brief Like @ref auto_mask_sigma, but instead specifies an absolute
+ * level where to stop generation of the auto-mask.
+ */
+ std::optional<double> absolute_auto_mask_threshold;
+
/**
* If @c true, maintain a list of components while performing deconvolution.
* This works with the @ref AlgorithmType::kGenericClean and @ref
diff --git a/cpp/test/test_python_deconvolution.cc b/cpp/test/test_python_deconvolution.cc
index 013a91b1..b2e5ad68 100644
--- a/cpp/test/test_python_deconvolution.cc
+++ b/cpp/test/test_python_deconvolution.cc
@@ -32,7 +32,7 @@ class PythonFileFixture {
settings.pixel_scale.x = kPixelScale;
settings.pixel_scale.y = kPixelScale;
settings.minor_iteration_count = 1000;
- settings.threshold = 1.0e-8;
+ settings.absolute_threshold = 1.0e-8;
settings.algorithm_type = AlgorithmType::kPython;
settings.python.filename = kPythonFilename;
}
diff --git a/cpp/test/test_radler.cc b/cpp/test/test_radler.cc
index 9ddd8319..e51bd8d5 100644
--- a/cpp/test/test_radler.cc
+++ b/cpp/test/test_radler.cc
@@ -86,7 +86,7 @@ struct SettingsFixture {
settings.pixel_scale.x = kPixelScale;
settings.pixel_scale.y = kPixelScale;
settings.minor_iteration_count = 1000;
- settings.threshold = 1.0e-8;
+ settings.absolute_threshold = 1.0e-8;
}
Settings settings;
@@ -191,7 +191,7 @@ BOOST_AUTO_TEST_CASE(diffuse_source) {
Settings settings;
settings.algorithm_type = AlgorithmType::kMultiscale;
- settings.threshold = 1.0e-8;
+ settings.absolute_threshold = 1.0e-8;
settings.major_iteration_count = 30;
settings.pixel_scale.x = imgReader.PixelSizeX();
settings.pixel_scale.y = imgReader.PixelSizeY();
@@ -199,7 +199,7 @@ BOOST_AUTO_TEST_CASE(diffuse_source) {
settings.trimmed_image_height = imgReader.ImageHeight();
settings.minor_iteration_count = 300;
settings.minor_loop_gain = 0.8;
- settings.auto_mask_sigma = 4;
+ settings.auto_mask_sigma = 4.0;
const double beamScale = imgReader.BeamMajorAxisRad();
@@ -209,18 +209,19 @@ BOOST_AUTO_TEST_CASE(diffuse_source) {
const int major_iteration_count = 0;
radler.Perform(reached_threshold, major_iteration_count);
- BOOST_REQUIRE(radler.IterationNumber() <= settings.minor_iteration_count);
+ BOOST_CHECK_LE(radler.IterationNumber(), settings.minor_iteration_count);
+ BOOST_CHECK_GE(radler.IterationNumber(), 100);
const double max_residual = residual_image.Max();
const double rms_residual = residual_image.RMS();
// Checks that RMS value in residual image went down at least by 25%. A
// reduction is expected as we remove the peaks from the dirty image.
- BOOST_REQUIRE(rms_residual < 0.75 * rms_dirty_image);
+ BOOST_CHECK_LT(rms_residual, 0.75 * rms_dirty_image);
// Checks that the components in the dirty image are reduced in the first
- // iteration. We expect thehighest peak to be reduced by 90% in this case.
- BOOST_REQUIRE(max_residual < 0.1 * max_dirty_image);
+ // iteration. We expect the highest peak to be reduced by 90% in this case.
+ BOOST_CHECK_LT(max_residual, 0.1 * max_dirty_image);
}
BOOST_AUTO_TEST_SUITE_END()
diff --git a/python/docstrings/settings_docstrings.h b/python/docstrings/settings_docstrings.h
index 0881342d..f76c0c88 100644
--- a/python/docstrings/settings_docstrings.h
+++ b/python/docstrings/settings_docstrings.h
@@ -205,6 +205,17 @@ static const char *__doc_radler_Settings_SpectralFitting_terms =
R"doc(Number of spectral terms to constrain the channels to, or zero to
disable.)doc";
+static const char *__doc_radler_Settings_absolute_auto_mask_threshold =
+R"doc(Like auto_mask_sigma, but instead specifies an absolute level where to
+stop generation of the auto-mask.)doc";
+
+static const char *__doc_radler_Settings_absolute_threshold =
+R"doc(Value in Jy that defines when to stop cleaning. Radler::Perform() will
+stop its major iteration and set ``reached_major_threshold``=false
+when the peak residual flux is below the given threshold. The default
+value is 0.0, which means that Radler will keep continuing until
+another criterion (e.g. nr. of iterations) is reached.)doc";
+
static const char *__doc_radler_Settings_algorithm_type =
R"doc(@{ These deconvolution settings are algorithm-specific. For each
algorithm type, a single struct holds all algorithm-specific settings
@@ -357,13 +368,6 @@ set ``reached_major_threshold``=false.)doc";
static const char *__doc_radler_Settings_thread_count = R"doc()doc";
-static const char *__doc_radler_Settings_threshold =
-R"doc(Value in Jy that defines when to stop cleaning. Radler::Perform() will
-stop its major iteration and set ``reached_major_threshold``=false
-when the peak residual flux is below the given threshold. The default
-value is 0.0, which means that Radler will keep continuing until
-another criterion (e.g. nr. of iterations) is reached.)doc";
-
static const char *__doc_radler_Settings_trimmed_image_height = R"doc(Trimmed image height)doc";
static const char *__doc_radler_Settings_trimmed_image_width =
diff --git a/python/pysettings.cc b/python/pysettings.cc
index 35b33be7..56854b9f 100644
--- a/python/pysettings.cc
+++ b/python/pysettings.cc
@@ -71,8 +71,9 @@ void init_settings(py::module& m) {
DOC(radler_Settings_linked_polarizations))
.def_readwrite("parallel", &radler::Settings::parallel,
DOC(radler_Settings_parallel))
- .def_readwrite("threshold", &radler::Settings::threshold,
- DOC(radler_Settings_threshold))
+ .def_readwrite("absolute_threshold",
+ &radler::Settings::absolute_threshold,
+ DOC(radler_Settings_absolute_threshold))
.def_readwrite("minor_loop_gain", &radler::Settings::minor_loop_gain,
DOC(radler_Settings_minor_loop_gain))
.def_readwrite("major_loop_gain", &radler::Settings::major_loop_gain,
@@ -82,6 +83,9 @@ void init_settings(py::module& m) {
DOC(radler_Settings_auto_threshold_sigma))
.def_readwrite("auto_mask_sigma", &radler::Settings::auto_mask_sigma,
DOC(radler_Settings_auto_mask_sigma))
+ .def_readwrite("absolute_auto_mask_threshold",
+ &radler::Settings::absolute_auto_mask_threshold,
+ DOC(radler_Settings_absolute_auto_mask_threshold))
.def_readwrite("save_source_list", &radler::Settings::save_source_list,
DOC(radler_Settings_save_source_list))
.def_readwrite("minor_iteration_count",
diff --git a/python/test/test_radler.py b/python/test/test_radler.py
index 414b95ed..61d27431 100644
--- a/python/test/test_radler.py
+++ b/python/test/test_radler.py
@@ -45,7 +45,7 @@ def get_settings():
settings.pixel_scale.x = PIXEL_SCALE
settings.pixel_scale.y = PIXEL_SCALE
settings.minor_iteration_count = MINOR_ITERATION_COUNT
- settings.threshold = 1e-8
+ settings.absolute_threshold = 1e-8
return settings
@@ -126,7 +126,7 @@ def test_input_dtype(settings):
residual = residual.astype(np.float32)
rd.Radler(settings, psf, residual, model, BEAM_SIZE)
- model = model.astype(np.int)
+ model = model.astype(int)
with pytest.raises(TypeError):
rd.Radler(settings, psf, residual, model, BEAM_SIZE)
model = model.astype(np.float32)
diff --git a/python/test/test_settings.py b/python/test/test_settings.py
index 3eb4e61b..3d54814e 100644
--- a/python/test/test_settings.py
+++ b/python/test/test_settings.py
@@ -25,7 +25,7 @@ def test_layout():
)
assert nested_structs == nested_structs_ref
- n_properties_ref = 33
+ n_properties_ref = 34
properties = set(filter(lambda x: re.match("^[a-z]+", x), dir(settings)))
assert len(properties) == n_properties_ref
@@ -45,7 +45,7 @@ def test_default():
assert settings.parallel.grid_width == 1
assert settings.parallel.grid_height == 1
assert settings.parallel.max_threads > 0
- assert settings.threshold == 0.0
+ assert settings.absolute_threshold == 0.0
assert settings.minor_loop_gain == 0.1
assert settings.major_loop_gain == 1.0
assert settings.auto_threshold_sigma == None
diff --git a/python/test/test_work_table.py b/python/test/test_work_table.py
index 5a114dfe..a8ca36fd 100644
--- a/python/test/test_work_table.py
+++ b/python/test/test_work_table.py
@@ -166,8 +166,8 @@ def test_add_entries_wrong_type():
entries = [rd.WorkTableEntry(), rd.WorkTableEntry(), rd.WorkTableEntry()]
psf = np.ones((4, 4), dtype=np.float64)
- residual = np.ones((4, 4), dtype=np.int)
- model = np.ones((4, 4), np.complex)
+ residual = np.ones((4, 4), dtype=int)
+ model = np.ones((4, 4), complex)
with pytest.raises(TypeError):
entries[0].psfs.append(psf)
--
GitLab