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
192
193
194
195
196
197
198
|
// SPDX-License-Identifier: LGPL-3.0-only
#ifndef RADLER_ALGORITHMS_IUWT_DECONVOLUTION_ALGORITHM_H_
#define RADLER_ALGORITHMS_IUWT_DECONVOLUTION_ALGORITHM_H_
#include <vector>
#include <aocommon/image.h>
#include <aocommon/staticfor.h>
#include <aocommon/uvector.h>
#include "image_set.h"
#include "algorithms/iuwt/image_analysis.h"
#include "algorithms/iuwt/iuwt_decomposition.h"
namespace radler::algorithms {
class IuwtDeconvolutionAlgorithm {
public:
IuwtDeconvolutionAlgorithm(size_t width, size_t height, float minor_loop_gain,
float major_loop_gain, float clean_border,
bool allow_negative_components, const bool* mask,
float absolute_threshold,
float threshold_sigma_level = 4.0,
float tolerance = 0.75);
IuwtDeconvolutionAlgorithm(const IuwtDeconvolutionAlgorithm&) = default;
IuwtDeconvolutionAlgorithm(IuwtDeconvolutionAlgorithm&&) = default;
IuwtDeconvolutionAlgorithm& operator=(const IuwtDeconvolutionAlgorithm&) =
default;
IuwtDeconvolutionAlgorithm& operator=(IuwtDeconvolutionAlgorithm&&) = default;
float PerformMajorIteration(size_t& iter_counter, size_t n_iter,
ImageSet& model_set, ImageSet& dirty_set,
const std::vector<aocommon::Image>& psfs,
bool& reached_major_threshold);
void Subtract(float* dest, const aocommon::Image& rhs);
void Subtract(aocommon::Image& dest, const aocommon::Image& rhs) {
Subtract(dest.Data(), rhs);
}
private:
struct ValComponent {
ValComponent() {}
ValComponent(size_t _x, size_t _y, int _scale, float _val = 0.0)
: x(_x), y(_y), scale(_scale), val(_val) {}
std::string ToString() const {
std::ostringstream str;
str << x << ',' << y << ", scale " << scale;
return str.str();
}
size_t x;
size_t y;
int scale;
float val;
};
struct ScaleResponse {
float rms;
float peak_response;
float peak_response_to_next_scale;
float convolvedPeakResponse;
double b_major;
double b_minor;
double b_pa;
size_t convolved_area;
};
float GetMaxAbsWithoutMask(const aocommon::Image& data, size_t& x, size_t& y,
size_t width);
float GetMaxAbsWithMask(const aocommon::Image& data, size_t& x, size_t& y,
size_t width);
float GetMaxAbs(const aocommon::Image& data, size_t& x, size_t& y,
size_t width) {
if (mask_ == nullptr)
return GetMaxAbsWithoutMask(data, x, y, width);
else
return GetMaxAbsWithMask(data, x, y, width);
}
void MeasureRMSPerScale(const float* image, const float* convolved_image,
float* scratch, size_t end_scale,
std::vector<ScaleResponse>& psf_response);
float Mad(const float* dest);
float DotProduct(const aocommon::Image& lhs, const aocommon::Image& rhs);
void BoundingBox(size_t& x1, size_t& y1, size_t& x2, size_t& y2,
const aocommon::Image& image, size_t width, size_t height);
void AdjustBox(size_t& x1, size_t& y1, size_t& x2, size_t& y2, size_t width,
size_t height, int end_scale);
void Trim(aocommon::Image& dest, const float* source, size_t old_width,
size_t x1, size_t y1, size_t x2, size_t y2);
void Trim(aocommon::Image& dest, const aocommon::Image& source, size_t x1,
size_t y1, size_t x2, size_t y2) {
Trim(dest, source.Data(), source.Width(), x1, y1, x2, y2);
}
void TrimPsf(aocommon::Image& dest, const aocommon::Image& source,
size_t new_width, size_t new_height) {
const size_t oldWidth = source.Width();
const size_t oldHeight = source.Height();
Trim(dest, source, (oldWidth - new_width) / 2, (oldHeight - new_height) / 2,
(oldWidth + new_width) / 2, (oldHeight + new_height) / 2);
}
void Untrim(aocommon::Image& image, size_t width, size_t height, size_t x1,
size_t y1, size_t x2, size_t y2);
float Snr(const iuwt::IuwtDecomposition& noisy_image,
const iuwt::IuwtDecomposition& model) const;
bool RunConjugateGradient(iuwt::IuwtDecomposition& iuwt,
const iuwt::IuwtMask& mask,
aocommon::Image& masked_dirty,
aocommon::Image& structure_model,
aocommon::Image& scratch,
const aocommon::Image& psf_kernel, size_t width,
size_t height);
bool FillAndDeconvolveStructure(
iuwt::IuwtDecomposition& iuwt, aocommon::Image& dirty,
ImageSet& structure_model_full, aocommon::Image& scratch,
const aocommon::Image& psf, const aocommon::Image& psf_kernel,
const std::vector<aocommon::Image>& psfs, size_t cur_end_scale,
size_t cur_min_scale, size_t width, size_t height,
const aocommon::UVector<float>& thresholds,
const iuwt::image_analysis::Component& max_comp, bool allow_trimming,
const bool* prior_mask);
bool FindAndDeconvolveStructure(
iuwt::IuwtDecomposition& iuwt, aocommon::Image& dirty,
const aocommon::Image& psf, const aocommon::Image& psf_kernel,
const std::vector<aocommon::Image>& psfs, aocommon::Image& scratch,
ImageSet& structure_model_full, size_t cur_end_scale,
size_t cur_min_scale, std::vector<ValComponent>& max_components);
void PerformSubImageFitAll(iuwt::IuwtDecomposition& iuwt,
const iuwt::IuwtMask& mask,
const aocommon::Image& structure_model,
aocommon::Image& scratch_a,
aocommon::Image& scratch_b,
const iuwt::image_analysis::Component& max_comp,
ImageSet& fitted_model, const aocommon::Image& psf,
const std::vector<aocommon::Image>& psfs,
const aocommon::Image& dirty);
void PerformSubImageFitSingle(
iuwt::IuwtDecomposition& iuwt, const iuwt::IuwtMask& mask,
const aocommon::Image& structure_model, aocommon::Image& scratch_b,
const iuwt::image_analysis::Component& max_comp,
const aocommon::Image& psf, aocommon::Image& sub_dirty,
float* fitted_sub_model, aocommon::UVector<float>& correction_factor);
float PerformSubImageComponentFitBoxed(
iuwt::IuwtDecomposition& iuwt, const iuwt::IuwtMask& mask,
const std::vector<iuwt::image_analysis::Component2D>& area,
aocommon::Image& scratch, aocommon::Image& masked_dirty,
const aocommon::Image& psf, const aocommon::Image& psf_kernel, size_t x1,
size_t y1, size_t x2, size_t y2);
float PerformSubImageComponentFit(
iuwt::IuwtDecomposition& iuwt, const iuwt::IuwtMask& mask,
const std::vector<iuwt::image_analysis::Component2D>& area,
aocommon::Image& scratch, aocommon::Image& masked_dirty,
const aocommon::Image& psf_kernel, size_t x_offset, size_t y_offset);
float CentralPeak(const aocommon::Image& data) {
return data[width_ / 2 + (height_ / 2) * width_];
}
size_t width_;
size_t height_;
size_t cur_box_x_start_;
size_t cur_box_x_end_;
size_t cur_box_y_start_;
size_t cur_box_y_end_;
float minor_loop_gain_;
float major_loop_gain_;
float clean_border_;
const bool* mask_;
float absolute_threshold_;
float threshold_sigma_level_;
float tolerance_;
aocommon::UVector<float> rmses_;
std::vector<ScaleResponse> psf_response_;
bool allow_negative_components_;
ImageSet* dirty_set_;
};
} // namespace radler::algorithms
#endif // RADLER_ALGORITHMS_IUWT_DECONVOLUTION_ALGORITHM_H_
|
