Skip to content
Snippets Groups Projects
Commit c2478a59 authored by Andre Offringa's avatar Andre Offringa
Browse files

Add absolute auto-masking threshold

parent fd49e089
No related branches found
No related tags found
1 merge request!119Add absolute auto-masking threshold
Pipeline #60230 passed
......@@ -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);
......
......@@ -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
......
......@@ -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;
}
......
......@@ -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 the highest peak to be reduced by 90% in this case.
BOOST_REQUIRE(max_residual < 0.1 * max_dirty_image);
BOOST_CHECK_LT(max_residual, 0.1 * max_dirty_image);
}
BOOST_AUTO_TEST_SUITE_END()
......
......@@ -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 =
......
......@@ -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",
......
......@@ -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)
......
......@@ -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
......
......@@ -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)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment