1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
|
// SPDX-License-Identifier: LGPL-3.0-only
#ifndef RADLER_ALGORITHMS_DECONVOLUTION_ALGORITHM_H_
#define RADLER_ALGORITHMS_DECONVOLUTION_ALGORITHM_H_
#include <cmath>
#include <limits>
#include <string>
#include <aocommon/image.h>
#include <aocommon/logger.h>
#include <aocommon/polarization.h>
#include <aocommon/uvector.h>
#include <schaapcommon/fitters/spectralfitter.h>
#include "image_set.h"
namespace radler::algorithms {
/**
* Class to capture information returned by
* @ref DeconvolutionAlgorithm::ExecuteMajorIteration().
*/
struct DeconvolutionResult {
/**
* The peak (in Jy) of the highest residual value, or zero if unknown
* or irrelevant.
*/
float final_peak_value = 0.0;
/**
* A value of @c true indicates that the function should be called again
* after a predict-inversion round. This is e.g. the case when the major
* iteration threshold was reached of a clean algorithm.
*/
bool another_iteration_required = false;
/**
* If @c true, the results of this iteration are worse than at the start.
* With clean algorithms, this happens when the peak value is (significantly)
* higher than at the start. If @c true, @ref another_iteration_required
* should normally be @c false indicating no progress is made. When using
* parallel deconvolution, a value of @c true will cause the results of the
* diverging sub-image to be reset. See also:
* @ref Settings::divergence_limit .
*/
bool is_diverging = false;
};
class DeconvolutionAlgorithm {
public:
virtual ~DeconvolutionAlgorithm() = default;
// Alhough deleting the copy-assignment violates the rule of three, it is not
// used. Defining it would only result in unused and untested code.
DeconvolutionAlgorithm& operator=(const DeconvolutionAlgorithm&) = delete;
DeconvolutionAlgorithm(DeconvolutionAlgorithm&&) = delete;
DeconvolutionAlgorithm& operator=(DeconvolutionAlgorithm&&) = delete;
virtual DeconvolutionResult ExecuteMajorIteration(
ImageSet& data_image, ImageSet& model_image,
const std::vector<aocommon::Image>& psf_images) = 0;
virtual std::unique_ptr<DeconvolutionAlgorithm> Clone() const = 0;
void SetMaxIterations(size_t max_iterations) {
settings_.max_iterations = max_iterations;
}
void SetThreshold(float threshold) { settings_.threshold = threshold; }
void SetMajorIterationThreshold(float major_iteration_threshold) {
settings_.major_iteration_threshold = major_iteration_threshold;
}
void SetMinorLoopGain(float minor_loop_gain) {
settings_.minor_loop_gain = minor_loop_gain;
}
void SetMajorLoopGain(float major_loop_gain) {
settings_.major_loop_gain = major_loop_gain;
}
void SetAllowNegativeComponents(bool allow_negative_components) {
settings_.allow_negative_components = allow_negative_components;
}
void SetStopOnNegativeComponents(bool stop_on_negative_component) {
settings_.stop_on_negative_component = stop_on_negative_component;
}
void SetCleanBorderRatio(float clean_border_ratio) {
settings_.clean_border_ratio = clean_border_ratio;
}
void SetCleanMask(const bool* clean_mask) {
settings_.clean_mask = clean_mask;
}
void SetDivergenceLimit(float divergence_limit) {
settings_.divergence_limit = divergence_limit;
}
void SetLogReceiver(aocommon::LogReceiver& log_receiver) {
log_receiver_ = &log_receiver;
}
size_t MaxIterations() const { return settings_.max_iterations; }
float Threshold() const { return settings_.threshold; }
float MajorIterationThreshold() const {
return settings_.major_iteration_threshold;
}
float MinorLoopGain() const { return settings_.minor_loop_gain; }
float MajorLoopGain() const { return settings_.major_loop_gain; }
float CleanBorderRatio() const { return settings_.clean_border_ratio; }
bool AllowNegativeComponents() const {
return settings_.allow_negative_components;
}
bool StopOnNegativeComponents() const {
return settings_.stop_on_negative_component;
}
const bool* CleanMask() const { return settings_.clean_mask; }
size_t IterationNumber() const { return iteration_number_; }
float DivergenceLimit() const { return settings_.divergence_limit; }
void SetIterationNumber(size_t iteration_number) {
iteration_number_ = iteration_number;
}
void SetSpectralFitter(
std::unique_ptr<schaapcommon::fitters::SpectralFitter> fitter,
const size_t n_polarizations) {
spectral_fitter_ = std::move(fitter);
n_polarizations_ = n_polarizations;
}
void SetSpectrallyForcedImages(std::vector<aocommon::Image>&& images) {
spectral_fitter_->SetForcedTerms(std::move(images));
}
const schaapcommon::fitters::SpectralFitter& Fitter() const {
return *spectral_fitter_;
}
void SetRmsFactorImage(aocommon::Image&& image) {
rms_factor_image_ = std::move(image);
}
const aocommon::Image& RmsFactorImage() const { return rms_factor_image_; }
/**
* Fit an array of values to a curve, and replace those values
* with the curve values. The position parameters are used when
* constraint fitting is used. Different polarizations are fitted
* independently.
* @param values is an array the size of the ImageSet (so npolarizations x
* nchannels).
*/
void PerformSpectralFit(float* values, size_t x, size_t y);
protected:
DeconvolutionAlgorithm();
DeconvolutionAlgorithm(const DeconvolutionAlgorithm&);
aocommon::LogReceiver& LogReceiver() { return *log_receiver_; };
private:
// Using a settings struct simplifies the constructors.
struct {
float threshold = 0.0;
float major_iteration_threshold = 0.0;
float minor_loop_gain = 0.1;
float major_loop_gain = 1.0;
float clean_border_ratio = 0.05;
size_t max_iterations = 500;
float divergence_limit = 4.0;
bool allow_negative_components = true;
bool stop_on_negative_component = false;
const bool* clean_mask = nullptr;
} settings_;
aocommon::LogReceiver* log_receiver_ = nullptr;
std::vector<float> fitting_scratch_;
std::unique_ptr<schaapcommon::fitters::SpectralFitter> spectral_fitter_;
aocommon::Image rms_factor_image_;
size_t iteration_number_ = 0;
size_t n_polarizations_ = 1;
};
} // namespace radler::algorithms
#endif // RADLER_ALGORITHMS_DECONVOLUTION_ALGORITHM_H_
|
