PURIFY
Next-generation radio interferometric imaging
measurement_operator.cc
Go to the documentation of this file.
1 #include <iomanip>
2 #include "catch2/catch_all.hpp"
3 
4 #include "purify/directories.h"
5 #include "purify/kernels.h"
6 #include "purify/operators.h"
7 #include "purify/pfitsio.h"
8 #include "purify/test_data.h"
9 #include "purify/utilities.h"
10 #include "purify/wproj_operators.h"
11 #include <sopt/power_method.h>
12 
13 using namespace purify;
14 using Catch::Approx;
15 
16 TEST_CASE("regression_degrid") {
17  const t_int imsizex = 256;
18  const t_int imsizey = 256;
19  Vector<t_real> u = Vector<t_real>::Random(10) * imsizex / 2;
20  Vector<t_real> v = Vector<t_real>::Random(10) * imsizey / 2;
21  Vector<t_complex> input = Vector<t_complex>::Zero(imsizex * imsizey);
22  input(static_cast<t_int>(imsizex * 0.5 + imsizey * 0.5 * imsizex)) = 1.;
23  const t_uint M = u.size();
24  const t_real oversample_ratio = 2;
25  const t_int ftsizev = std::floor(imsizey * oversample_ratio);
26  const t_int ftsizeu = std::floor(imsizex * oversample_ratio);
27  const t_uint Ju = 8;
28  const t_uint Jv = 8;
29 
30  const Vector<t_complex> y = Vector<t_complex>::Ones(u.size());
31  CAPTURE(u);
32  CAPTURE(v);
33 
34  SECTION("kb") {
35  const kernels::kernel kernel = kernels::kernel::kb;
36  const std::shared_ptr<sopt::LinearTransform<Vector<t_complex>>> measure_op =
37  measurementoperator::init_degrid_operator_2d<Vector<t_complex>>(
38  u, v, Vector<t_real>::Zero(M), Vector<t_complex>::Ones(M), imsizey, imsizex,
39  oversample_ratio, kernel, Ju, Jv);
40  const Vector<t_complex> y_test = *measure_op * input;
41  CAPTURE(y_test.array() / y.array());
42  const t_real max_test = y_test.cwiseAbs().mean();
43  CAPTURE(y_test / max_test);
44  CAPTURE(y);
45  CHECK((y_test / max_test).isApprox((y), 1e-6));
46  }
47  SECTION("pswf") {
48  const kernels::kernel kernel = kernels::kernel::pswf;
49  const auto measure_op = measurementoperator::init_degrid_operator_2d<Vector<t_complex>>(
50  u, v, Vector<t_real>::Zero(M), Vector<t_complex>::Ones(M), imsizey, imsizex,
51  oversample_ratio, kernel, 6, 6);
52  const Vector<t_complex> y_test = *measure_op * input;
53  CAPTURE(y_test.array() / y.array());
54  const t_real max_test = y_test.cwiseAbs().mean();
55  CAPTURE(y_test / max_test);
56  CAPTURE(y);
57  CHECK((y_test / max_test).isApprox((y), 1e-3));
58  }
59  SECTION("gauss") {
60  const kernels::kernel kernel = kernels::kernel::gauss;
61  const auto measure_op = measurementoperator::init_degrid_operator_2d<Vector<t_complex>>(
62  u, v, Vector<t_real>::Zero(M), Vector<t_complex>::Ones(M), imsizey, imsizex,
63  oversample_ratio, kernel, 10, 10);
64  const Vector<t_complex> y_test = *measure_op * input;
65  CAPTURE(y_test.array() / y.array());
66  const t_real max_test = y_test.cwiseAbs().mean();
67  CAPTURE(y_test / max_test);
68  CAPTURE(y);
69  CHECK((y_test / max_test).isApprox((y), 1e-4));
70  }
71 }
72 
73 TEST_CASE("flux units") {
74  const t_real oversample_ratio = 2;
75  Vector<t_real> u = Vector<t_real>::Random(10);
76  Vector<t_real> v = Vector<t_real>::Random(10);
77  const t_uint M = u.size();
78  const Vector<t_complex> y = Vector<t_complex>::Ones(u.size());
79  SECTION("kb") {
80  const kernels::kernel kernel = kernels::kernel::kb;
81  for (auto& J : {4, 5, 6, 7, 8}) {
82  for (auto& imsize : {128, 256, 512}) {
83  const std::shared_ptr<sopt::LinearTransform<Vector<t_complex>>> measure_op =
84  measurementoperator::init_degrid_operator_2d<Vector<t_complex>>(
85  u * imsize / 2, v * imsize / 2, Vector<t_real>::Zero(M), Vector<t_complex>::Ones(M),
86  imsize, imsize, oversample_ratio, kernel, J, J);
87  Vector<t_complex> input = Vector<t_complex>::Zero(imsize * imsize);
88  input(static_cast<t_int>(imsize * 0.5 + imsize * 0.5 * imsize)) = 1.;
89  const Vector<t_complex> y_test = (*measure_op * input).eval();
90  CAPTURE(y_test.cwiseAbs().mean());
91  CAPTURE(y);
92  CAPTURE(y_test);
93  CAPTURE(J);
94  CAPTURE(imsize);
95  CHECK(y_test.isApprox(y, 1e-3));
96  const Vector<t_complex> psf = (measure_op->adjoint() * y) * 1. / M;
97  CHECK(std::real(psf(static_cast<t_int>(imsize * 0.5 + imsize * 0.5 * imsize))) ==
98  Approx(1.).margin(0.001));
99  }
100  }
101  }
102  SECTION("pswf") {
103  const kernels::kernel kernel = kernels::kernel::pswf;
104  for (auto& J : {4, 5, 6, 7, 8}) {
105  for (auto& imsize : {128, 256, 512}) {
106  if (J != 6)
107  CHECK_THROWS(measurementoperator::init_degrid_operator_2d<Vector<t_complex>>(
108  u * imsize / 2, v * imsize / 2, Vector<t_real>::Zero(M), Vector<t_complex>::Ones(M),
109  imsize, imsize, oversample_ratio, kernel, J, J));
110  else {
111  const std::shared_ptr<sopt::LinearTransform<Vector<t_complex>>> measure_op =
112  measurementoperator::init_degrid_operator_2d<Vector<t_complex>>(
113  u * imsize / 2, v * imsize / 2, Vector<t_real>::Zero(M),
114  Vector<t_complex>::Ones(M), imsize, imsize, oversample_ratio, kernel, J, J);
115  Vector<t_complex> input = Vector<t_complex>::Zero(imsize * imsize);
116  input(static_cast<t_int>(imsize * 0.5 + imsize * 0.5 * imsize)) = 1.;
117  const Vector<t_complex> y_test = (*measure_op * input).eval();
118  CAPTURE(y_test.cwiseAbs().mean());
119  CAPTURE(y);
120  CAPTURE(y_test);
121  CAPTURE(J);
122  CAPTURE(imsize);
123  CHECK(y_test.isApprox(y, 1e-3));
124  const Vector<t_complex> psf = (measure_op->adjoint() * y) * 1. / M;
125  CHECK(std::real(psf(static_cast<t_int>(imsize * 0.5 + imsize * 0.5 * imsize))) ==
126  Approx(1.).margin(0.001));
127  }
128  }
129  }
130  }
131  SECTION("gauss") {
132  const kernels::kernel kernel = kernels::kernel::gauss;
133  for (auto& J : {4, 5, 6, 7, 8}) {
134  for (auto& imsize : {128, 256, 512}) {
135  const std::shared_ptr<sopt::LinearTransform<Vector<t_complex>>> measure_op =
136  measurementoperator::init_degrid_operator_2d<Vector<t_complex>>(
137  u * imsize / 2, v * imsize / 2, Vector<t_real>::Zero(M), Vector<t_complex>::Ones(M),
138  imsize, imsize, oversample_ratio, kernel, J, J);
139  Vector<t_complex> input = Vector<t_complex>::Zero(imsize * imsize);
140  input(static_cast<t_int>(imsize * 0.5 + imsize * 0.5 * imsize)) = 1.;
141  const Vector<t_complex> y_test = (*measure_op * input).eval();
142  CAPTURE(y_test.cwiseAbs().mean());
143  CAPTURE(y);
144  CAPTURE(y_test);
145  CAPTURE(J);
146  CAPTURE(imsize);
147  CHECK(y_test.isApprox(y, 1e-2));
148  const Vector<t_complex> psf = (measure_op->adjoint() * y) * 1. / M;
149  CHECK(std::real(psf(static_cast<t_int>(imsize * 0.5 + imsize * 0.5 * imsize))) ==
150  Approx(1.).margin(0.01));
151  }
152  }
153  }
154  SECTION("box") {
155  const kernels::kernel kernel = kernels::kernel::box;
156  for (auto& J : {4, 5, 6, 7, 8}) {
157  for (auto& imsize : {128, 256, 512}) {
158  const std::shared_ptr<sopt::LinearTransform<Vector<t_complex>>> measure_op =
159  measurementoperator::init_degrid_operator_2d<Vector<t_complex>>(
160  u * imsize / 2, v * imsize / 2, Vector<t_real>::Zero(M), Vector<t_complex>::Ones(M),
161  imsize, imsize, oversample_ratio, kernel, J, J);
162  Vector<t_complex> input = Vector<t_complex>::Zero(imsize * imsize);
163  input(static_cast<t_int>(imsize * 0.5 + imsize * 0.5 * imsize)) = 1.;
164  const Vector<t_complex> y_test = (*measure_op * input).eval();
165  CAPTURE(y_test.cwiseAbs().mean());
166  CAPTURE(y);
167  CAPTURE(y_test);
168  CAPTURE(J);
169  CAPTURE(imsize);
170  CHECK(y_test.isApprox(y, 1e-3));
171  const Vector<t_complex> psf = (measure_op->adjoint() * y) * 1. / M;
172  CHECK(std::real(psf(static_cast<t_int>(imsize * 0.5 + imsize * 0.5 * imsize))) ==
173  Approx(1.).margin(0.001));
174  }
175  }
176  }
177  SECTION("wproj kb") {
178  const kernels::kernel kernel = kernels::kernel::kb;
179  for (auto& J : {4, 5, 6, 7, 8}) {
180  for (auto& imsize : {128, 256, 512}) {
181  const std::shared_ptr<sopt::LinearTransform<Vector<t_complex>>> measure_op =
182  measurementoperator::init_degrid_operator_2d<Vector<t_complex>>(
183  u * imsize / 2, v * imsize / 2, Vector<t_real>::Zero(M), Vector<t_complex>::Ones(M),
184  imsize, imsize, oversample_ratio, kernel, J, J, false, 1e-6, 1e-6);
185  Vector<t_complex> input = Vector<t_complex>::Zero(imsize * imsize);
186  input(static_cast<t_int>(imsize * 0.5 + imsize * 0.5 * imsize)) = 1.;
187  const Vector<t_complex> y_test = (*measure_op * input).eval();
188  CAPTURE(y_test.cwiseAbs().mean());
189  CAPTURE(y);
190  CAPTURE(y_test);
191  CAPTURE(J);
192  CAPTURE(imsize);
193  CHECK(y_test.isApprox(y, 1e-3));
194  const Vector<t_complex> psf = (measure_op->adjoint() * y) * 1. / M;
195  CHECK(std::real(psf(static_cast<t_int>(imsize * 0.5 + imsize * 0.5 * imsize))) ==
196  Approx(1.).margin(0.001));
197  }
198  }
199  }
200  SECTION("wproj pswf") {
201  const kernels::kernel kernel = kernels::kernel::pswf;
202  for (auto& J : {4, 5, 6, 7, 8}) {
203  for (auto& imsize : {128, 256, 512}) {
204  if (J != 6)
205  CHECK_THROWS(measurementoperator::init_degrid_operator_2d<Vector<t_complex>>(
206  u * imsize / 2, v * imsize / 2, Vector<t_real>::Zero(M), Vector<t_complex>::Ones(M),
207  imsize, imsize, oversample_ratio, kernel, J, J, false, 1e-6, 1e-6));
208  else {
209  const std::shared_ptr<sopt::LinearTransform<Vector<t_complex>>> measure_op =
210  measurementoperator::init_degrid_operator_2d<Vector<t_complex>>(
211  u * imsize / 2, v * imsize / 2, Vector<t_real>::Zero(M),
212  Vector<t_complex>::Ones(M), imsize, imsize, oversample_ratio, kernel, J, J, false,
213  1e-6, 1e-6);
214  Vector<t_complex> input = Vector<t_complex>::Zero(imsize * imsize);
215  input(static_cast<t_int>(imsize * 0.5 + imsize * 0.5 * imsize)) = 1.;
216  const Vector<t_complex> y_test = (*measure_op * input).eval();
217  CAPTURE(y_test.cwiseAbs().mean());
218  CAPTURE(y);
219  CAPTURE(y_test);
220  CAPTURE(J);
221  CAPTURE(imsize);
222  CHECK(y_test.isApprox(y, 1e-3));
223  const Vector<t_complex> psf = (measure_op->adjoint() * y) * 1. / M;
224  CHECK(std::real(psf(static_cast<t_int>(imsize * 0.5 + imsize * 0.5 * imsize))) ==
225  Approx(1.).margin(0.001));
226  }
227  }
228  }
229  }
230  SECTION("wproj gauss") {
231  const kernels::kernel kernel = kernels::kernel::gauss;
232  for (auto& J : {4, 5, 6, 7, 8}) {
233  for (auto& imsize : {128, 256, 512}) {
234  const std::shared_ptr<sopt::LinearTransform<Vector<t_complex>>> measure_op =
235  measurementoperator::init_degrid_operator_2d<Vector<t_complex>>(
236  u * imsize / 2, v * imsize / 2, Vector<t_real>::Zero(M), Vector<t_complex>::Ones(M),
237  imsize, imsize, oversample_ratio, kernel, J, J, false, 1e-6, 1e-6);
238  Vector<t_complex> input = Vector<t_complex>::Zero(imsize * imsize);
239  input(static_cast<t_int>(imsize * 0.5 + imsize * 0.5 * imsize)) = 1.;
240  const Vector<t_complex> y_test = (*measure_op * input).eval();
241  CAPTURE(y_test.cwiseAbs().mean());
242  CAPTURE(y);
243  CAPTURE(y_test);
244  CAPTURE(J);
245  CAPTURE(imsize);
246  CHECK(y_test.isApprox(y, 1e-2));
247  const Vector<t_complex> psf = (measure_op->adjoint() * y) * 1. / M;
248  CHECK(std::real(psf(static_cast<t_int>(imsize * 0.5 + imsize * 0.5 * imsize))) ==
249  Approx(1.).margin(0.01));
250  }
251  }
252  }
253  SECTION("wproj box") {
254  const kernels::kernel kernel = kernels::kernel::box;
255  for (auto& J : {4, 5, 6, 7, 8}) {
256  for (auto& imsize : {128, 256, 512}) {
257  const std::shared_ptr<sopt::LinearTransform<Vector<t_complex>>> measure_op =
258  measurementoperator::init_degrid_operator_2d<Vector<t_complex>>(
259  u * imsize / 2, v * imsize / 2, Vector<t_real>::Zero(M), Vector<t_complex>::Ones(M),
260  imsize, imsize, oversample_ratio, kernel, J, J, false, 1e-6, 1e-6);
261  Vector<t_complex> input = Vector<t_complex>::Zero(imsize * imsize);
262  input(static_cast<t_int>(imsize * 0.5 + imsize * 0.5 * imsize)) = 1.;
263  const Vector<t_complex> y_test = (*measure_op * input).eval();
264  CAPTURE(y_test.cwiseAbs().mean());
265  CAPTURE(y);
266  CAPTURE(y_test);
267  CAPTURE(J);
268  CAPTURE(imsize);
269  CHECK(y_test.isApprox(y, 1e-2));
270  const Vector<t_complex> psf = (measure_op->adjoint() * y) * 1. / M;
271  CHECK(std::real(psf(static_cast<t_int>(imsize * 0.5 + imsize * 0.5 * imsize))) ==
272  Approx(1.).margin(0.001));
273  }
274  }
275  }
276 }
277 TEST_CASE("normed operator") {
278  const t_real oversample_ratio = 2;
279  const t_int power_iters = 1000;
280  const t_real power_tol = 1e-4;
281  Vector<t_real> u = Vector<t_real>::Random(10);
282  Vector<t_real> v = Vector<t_real>::Random(10);
283  const t_uint M = u.size();
284  const Vector<t_complex> y = Vector<t_complex>::Ones(u.size());
285  SECTION("kb") {
286  const kernels::kernel kernel = kernels::kernel::kb;
287  for (auto& J : {4, 5, 6, 7, 8}) {
288  for (auto& imsize : {128, 256, 512}) {
289  const Vector<t_complex> init = Vector<t_complex>::Ones(imsize * imsize);
290  auto power_method_result = sopt::algorithm::normalise_operator<Vector<t_complex>>(
291  measurementoperator::init_degrid_operator_2d<Vector<t_complex>>(
292  u * imsize / 2, v * imsize / 2, Vector<t_real>::Zero(M), Vector<t_complex>::Ones(M),
293  imsize, imsize, oversample_ratio, kernel, J, J),
294  power_iters, power_tol, init);
295  const t_real norm = std::get<0>(power_method_result);
296  const std::shared_ptr<sopt::LinearTransform<Vector<t_complex>>> measure_op =
297  std::get<2>(power_method_result);
298 
299  Vector<t_complex> input = Vector<t_complex>::Zero(imsize * imsize);
300  input(static_cast<t_int>(imsize * 0.5 + imsize * 0.5 * imsize)) = 1.;
301  const Vector<t_complex> y_test = (*measure_op * input).eval();
302  CAPTURE(y_test.cwiseAbs().mean());
303  CAPTURE(y);
304  CAPTURE(y_test);
305  CAPTURE(J);
306  CAPTURE(imsize);
307  CHECK((y_test * norm).isApprox(y, 1e-3));
308  const Vector<t_complex> psf = (measure_op->adjoint() * y) * 1. / M * norm;
309  CHECK(std::real(psf(static_cast<t_int>(imsize * 0.5 + imsize * 0.5 * imsize))) ==
310  Approx(1.).margin(0.001));
311  }
312  }
313  }
314 }
CHECK
#define CHECK(CONDITION, ERROR)
Definition: casa.cc:6
CHECK_THROWS
#define CHECK_THROWS(STATEMENT, ERROR)
Definition: casa.cc:8
operators_test::u
const std::vector< t_real > u
data for u coordinate
Definition: operators.cc:18
operators_test::v
const std::vector< t_real > v
data for v coordinate
Definition: operators.cc:20
purify::measurementoperator::init_degrid_operator_2d
std::shared_ptr< sopt::LinearTransform< T > > init_degrid_operator_2d(const Vector< t_real > &u, const Vector< t_real > &v, const Vector< t_real > &w, const Vector< t_complex > &weights, const t_uint &imsizey, const t_uint &imsizex, const t_real &oversample_ratio=2, const kernels::kernel kernel=kernels::kernel::kb, const t_uint Ju=4, const t_uint Jv=4, const bool w_stacking=false, const t_real &cellx=1, const t_real &celly=1)
Returns linear transform that is the standard degridding operator.
Definition: operators.h:608
TEST_CASE
TEST_CASE("regression_degrid")
Definition: measurement_operator.cc:16
wproj_operators.h
Generated by doxygen 1.9.1