PURIFY
Next-generation radio interferometric imaging
Functions
purify::details Namespace Reference

Functions

Sparse< t_complex > init_gridding_matrix_2d (const Vector< t_real > &u, const Vector< t_real > &v, const Vector< t_complex > &weights, const t_uint &imsizey_, const t_uint &imsizex_, const t_real &oversample_ratio, const std::function< t_real(t_real)> kernelu, const std::function< t_real(t_real)> kernelv, const t_uint Ju, const t_uint Jv)
 Construct gridding matrix. More...
 
Image< t_complex > init_correction2d (const t_real &oversample_ratio, const t_uint &imsizey_, const t_uint &imsizex_, const std::function< t_real(t_real)> ftkernelu, const std::function< t_real(t_real)> ftkernelv, const t_real &w_mean, const t_real &cellx, const t_real &celly)
 
Sparse< t_complex > init_gridding_matrix_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, const std::function< t_real(t_real)> &ftkerneluv, const std::function< t_real(t_real)> &kerneluv, const t_uint Ju, const t_uint Jw, const t_real cellx, const t_real celly, const t_real abs_error, const t_real rel_error, const dde_type dde)
 Construct gridding matrix with wprojection. More...
 
template<class STORAGE_INDEX_TYPE = t_int>
Sparse< t_complex, STORAGE_INDEX_TYPE > init_gridding_matrix_2d (const t_uint number_of_images, const std::vector< t_int > &image_index, const Vector< t_real > &u, const Vector< t_real > &v, const Vector< t_complex > &weights, const t_uint &imsizey_, const t_uint &imsizex_, const t_real &oversample_ratio, const std::function< t_real(t_real)> kernelu, const std::function< t_real(t_real)> kernelv, const t_uint Ju, const t_uint Jv)
 Construct all to all gridding matrix. More...
 
template<class STORAGE_INDEX_TYPE = t_int>
Sparse< t_complex, STORAGE_INDEX_TYPE > init_gridding_matrix_2d (const t_uint number_of_images, const std::vector< t_int > &image_index, const std::vector< t_real > &w_stacks, 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, const std::function< t_real(t_real)> &ftkerneluv, const std::function< t_real(t_real)> &kerneluv, const t_uint Ju, const t_uint Jw, const t_real cellx, const t_real celly, const t_real abs_error, const t_real rel_error, const dde_type dde)
 Construct all to all gridding matrix with wprojection. More...
 
template<class T , class... ARGS>
Sparse< t_complex > init_gridding_matrix_2d (const Sparse< T > &mixing_matrix, ARGS &&...args)
 Construct gridding matrix with mixing. More...
 
Image< t_complex > init_correction_radial_2d (const t_real oversample_ratio, const t_uint imsizey_, const t_uint imsizex_, const std::function< t_real(t_real)> &ftkerneluv, const t_real w_mean, const t_real cellx, const t_real celly)
 

Function Documentation

◆ init_correction2d()

Image< t_complex > purify::details::init_correction2d ( const t_real &  oversample_ratio,
const t_uint &  imsizey_,
const t_uint &  imsizex_,
const std::function< t_real(t_real)>  ftkernelu,
const std::function< t_real(t_real)>  ftkernelv,
const t_real &  w_mean,
const t_real &  cellx,
const t_real &  celly 
)

Given the Fourier transform of a gridding kernel, creates the scaling image for gridding correction.

Definition at line 47 of file operators.cc.

51  {
52  const t_uint ftsizeu_ = std::floor(imsizex_ * oversample_ratio);
53  const t_uint ftsizev_ = std::floor(imsizey_ * oversample_ratio);
54  const t_uint x_start = std::floor(ftsizeu_ * 0.5 - imsizex_ * 0.5);
55  const t_uint y_start = std::floor(ftsizev_ * 0.5 - imsizey_ * 0.5);
56 
57  Array<t_real> range;
58  range.setLinSpaced(std::max(ftsizeu_, ftsizev_), 0.5, std::max(ftsizeu_, ftsizev_) - 0.5);
59  return ((1e0 / range.segment(y_start, imsizey_).unaryExpr(ftkernelv)).matrix() *
60  (1e0 / range.segment(x_start, imsizex_).unaryExpr(ftkernelu)).matrix().transpose())
61  .array() *
62  t_complex(1., 0.) *
63  widefield::generate_chirp(w_mean, cellx, celly, imsizex_, imsizey_).array() * imsizex_ *
64  imsizey_;
65 }
purify::widefield::generate_chirp
Matrix< t_complex > generate_chirp(const t_real w_rate, const t_real cell_x, const t_real cell_y, const t_uint x_size, const t_uint y_size)
Generates image of chirp.
Definition: wide_field_utilities.cc:34

References purify::widefield::generate_chirp().

Referenced by purify::operators::base_padding_and_FFT_2d(), main(), and TEST_CASE().

◆ init_correction_radial_2d()

Image< t_complex > purify::details::init_correction_radial_2d ( const t_real  oversample_ratio,
const t_uint  imsizey_,
const t_uint  imsizex_,
const std::function< t_real(t_real)> &  ftkerneluv,
const t_real  w_mean,
const t_real  cellx,
const t_real  celly 
)

Given the Fourier transform of a radially symmetric gridding kernel, creates the scaling image for gridding correction.

Definition at line 116 of file wproj_operators.cc.

120  {
121  const t_uint ftsizeu_ = std::floor(imsizex_ * oversample_ratio);
122  const t_uint ftsizev_ = std::floor(imsizey_ * oversample_ratio);
123  const t_uint x_start = std::floor(ftsizeu_ * 0.5 - imsizex_ * 0.5);
124  const t_uint y_start = std::floor(ftsizev_ * 0.5 - imsizey_ * 0.5);
125 
126  Image<t_complex> gridding_correction = Image<t_complex>::Zero(imsizey_, imsizex_);
127  for (int i = 0; i < gridding_correction.rows(); i++)
128  for (int j = 0; j < gridding_correction.cols(); j++)
129  gridding_correction(i, j) =
130  1. / ftkerneluv(std::sqrt(std::pow((y_start + i + 0.5) / ftsizev_ - 0.5, 2) +
131  std::pow((j + 0.5 + x_start) / ftsizeu_ - 0.5, 2)));
132 
133  return gridding_correction.array() *
134  widefield::generate_chirp(w_mean, cellx, celly, imsizex_, imsizey_).array() * imsizex_ *
135  imsizey_;
136 }

References purify::widefield::generate_chirp().

Referenced by purify::operators::base_padding_and_FFT_2d().

◆ init_gridding_matrix_2d() [1/5]

template<class T , class... ARGS>
Sparse<t_complex> purify::details::init_gridding_matrix_2d ( const Sparse< T > &  mixing_matrix,
ARGS &&...  args 
)

Construct gridding matrix with mixing.

Definition at line 224 of file operators.h.

224  {
225  if (mixing_matrix.rows() * mixing_matrix.cols() < 2)
226  return init_gridding_matrix_2d(std::forward<ARGS>(args)...);
227  const Sparse<t_complex> G = init_gridding_matrix_2d(std::forward<ARGS>(args)...);
228  if (mixing_matrix.cols() != G.rows())
229  throw std::runtime_error(
230  "The columns of the mixing matrix do not match the number of visibilities");
231  return mixing_matrix * init_gridding_matrix_2d(std::forward<ARGS>(args)...);
232 };
purify::operators::init_gridding_matrix_2d
std::tuple< sopt::OperatorFunction< T >, sopt::OperatorFunction< T > > init_gridding_matrix_2d(ARGS &&...args)
constructs lambdas that apply degridding matrix with adjoint
Definition: operators.h:312

References init_gridding_matrix_2d().

◆ init_gridding_matrix_2d() [2/5]

template<class STORAGE_INDEX_TYPE = t_int>
Sparse<t_complex, STORAGE_INDEX_TYPE> purify::details::init_gridding_matrix_2d ( const t_uint  number_of_images,
const std::vector< t_int > &  image_index,
const std::vector< t_real > &  w_stacks,
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,
const std::function< t_real(t_real)> &  ftkerneluv,
const std::function< t_real(t_real)> &  kerneluv,
const t_uint  Ju,
const t_uint  Jw,
const t_real  cellx,
const t_real  celly,
const t_real  abs_error,
const t_real  rel_error,
const dde_type  dde 
)

Construct all to all gridding matrix with wprojection.

Definition at line 105 of file operators.h.

112  {
113  const t_uint ftsizev_ = std::floor(imsizey_ * oversample_ratio);
114  const t_uint ftsizeu_ = std::floor(imsizex_ * oversample_ratio);
115  const t_real du = widefield::pixel_to_lambda(cellx, imsizex_, oversample_ratio);
116  const t_real dv = widefield::pixel_to_lambda(celly, imsizey_, oversample_ratio);
117  PURIFY_HIGH_LOG("du x du: {}, {}", du, dv);
118  if (std::abs(du - dv) > 1e-6)
119  throw std::runtime_error(
120  "Field of view along l and m is not the same, this assumption is required for "
121  "w-projection.");
122  const t_uint rows = u.size();
123  const STORAGE_INDEX_TYPE cols = ftsizeu_ * ftsizev_ * number_of_images;
124  if (u.size() != v.size())
125  throw std::runtime_error(
126  "Size of u and v vectors are not the same for creating gridding matrix.");
127  if (u.size() != w.size())
128  throw std::runtime_error(
129  "Size of u and w vectors are not the same for creating gridding matrix.");
130  if (u.size() != weights.size())
131  throw std::runtime_error(
132  "Size of u and w vectors are not the same for creating gridding matrix.");
133  if (Ju > Jw)
134  throw std::runtime_error(
135  "w kernel size must be at least the size of w=0 kernel, must have Ju <= Jw.");
136  // count gridding coefficients with variable support size
137  Vector<t_int> total_coeffs = Vector<t_int>::Zero(w.size());
138  for (int i = 0; i < w.size(); i++) {
139  const t_int Ju_max = widefield::w_support(w(i) - w_stacks.at(image_index.at(i)), du,
140  static_cast<t_int>(Ju), static_cast<t_int>(Jw));
141  total_coeffs(i) = Ju_max * Ju_max;
142  }
143  const t_real num_of_coeffs = total_coeffs.array().cast<t_real>().sum();
144  PURIFY_HIGH_LOG("Using {} rows (coefficients per a row {}), and memory of {} MegaBytes",
145  total_coeffs.size(), total_coeffs.array().cast<t_real>().mean(),
146  16. * num_of_coeffs / std::pow(10., 6));
147  Sparse<t_complex, STORAGE_INDEX_TYPE> interpolation_matrix(static_cast<STORAGE_INDEX_TYPE>(rows),
148  cols);
149  try {
150  interpolation_matrix.reserve(total_coeffs);
151  } catch (std::bad_alloc e) {
152  throw std::runtime_error(
153  "Not enough memory for coefficients, choose upper limit on support size Jw.");
154  }
155 
156  const t_complex I(0., 1.);
157 
158  const t_uint max_evaluations = 1e8;
159  const t_real relative_error = rel_error;
160  const t_real absolute_error = abs_error;
161 
162  auto const ftkernel_radial = [&](const t_real l) -> t_real { return ftkerneluv(l); };
163 
164  std::int64_t coeffs_done = 0;
165  std::int64_t total = 0;
166 
167 #pragma omp parallel for
168  for (t_int m = 0; m < rows; ++m) {
169  // w_projection convolution setup
170  const t_real w_val = w(m) - w_stacks.at(image_index.at(m));
171  const t_int Ju_max = widefield::w_support(w_val, du, Ju, Jw);
172  t_uint evaluations = 0;
173  const t_int kwu = std::floor(u(m) - Ju_max * 0.5);
174  const t_int kwv = std::floor(v(m) - Ju_max * 0.5);
175 
176  for (t_int ju = 1; ju < Ju_max + 1; ++ju) {
177  for (t_int jv = 1; jv < Ju_max + 1; ++jv) {
178  const t_uint q = utilities::mod(kwu + ju, ftsizeu_);
179  const t_uint p = utilities::mod(kwv + jv, ftsizev_);
180  const STORAGE_INDEX_TYPE index =
181  static_cast<STORAGE_INDEX_TYPE>(utilities::sub2ind(p, q, ftsizev_, ftsizeu_)) +
182  static_cast<STORAGE_INDEX_TYPE>(image_index.at(m)) *
183  static_cast<STORAGE_INDEX_TYPE>(ftsizev_ * ftsizeu_);
184  interpolation_matrix.insert(static_cast<STORAGE_INDEX_TYPE>(m), index) =
185  std::exp(-2 * constant::pi * I * ((kwu + ju) * 0.5 + (kwv + jv) * 0.5)) * weights(m) *
186  ((dde == dde_type::wkernel_radial)
187  ? projection_kernels::exact_w_projection_integration_1d(
188  (u(m) - (kwu + ju)), (v(m) - (kwv + jv)), w_val, du, oversample_ratio,
189  ftkernel_radial, max_evaluations, absolute_error, relative_error,
190  (du > 1.) ? integration::method::p : integration::method::h, evaluations)
191  : projection_kernels::exact_w_projection_integration(
192  (u(m) - (kwu + ju)), (v(m) - (kwv + jv)), w_val, du, dv, oversample_ratio,
193  ftkernel_radial, ftkernel_radial, max_evaluations, absolute_error,
194  relative_error, integration::method::h, evaluations));
195 #pragma omp critical(add_eval)
196  {
197  coeffs_done++;
198  if (num_of_coeffs > 100) {
199  if ((coeffs_done % (static_cast<t_int>(num_of_coeffs / 100)) == 0)) {
200  PURIFY_LOW_LOG(
201  "w = {}, support = {}x{}, coeffs: {} of {}, {}%", w_val, Ju_max, Ju_max,
202  coeffs_done, num_of_coeffs,
203  static_cast<t_real>(coeffs_done) / static_cast<t_real>(num_of_coeffs) * 100.);
204  }
205  }
206  }
207  }
208  }
209  }
210  interpolation_matrix.makeCompressed();
211  return interpolation_matrix;
212 }
PURIFY_LOW_LOG
#define PURIFY_LOW_LOG(...)
Low priority message.
Definition: logging.h:207
PURIFY_HIGH_LOG
#define PURIFY_HIGH_LOG(...)
High priority message.
Definition: logging.h:203
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::constant::pi
const t_real pi
mathematical constant
Definition: types.h:70
purify::projection_kernels::exact_w_projection_integration_1d
t_complex exact_w_projection_integration_1d(const t_real u, const t_real v, const t_real w, const t_real du, const t_real oversample_ratio, const std::function< t_complex(t_real)> &ftkerneluv, const t_uint &max_evaluations, const t_real &absolute_error, const t_real &relative_error, const integration::method method, t_uint &evaluations)
Definition: wkernel_integration.cc:25
purify::projection_kernels::exact_w_projection_integration
t_complex exact_w_projection_integration(const t_real u, const t_real v, const t_real w, const t_real du, const t_real dv, const t_real oversample_ratio, const std::function< t_complex(t_real)> &ftkernelu, const std::function< t_complex(t_real)> &ftkernelv, const t_uint &max_evaluations, const t_real &absolute_error, const t_real &relative_error, const integration::method method, t_uint &evaluations)
numerical integration of chirp and kernel in image domain
Definition: wkernel_integration.cc:46
purify::utilities::sub2ind
t_int sub2ind(const t_int &row, const t_int &col, const t_int &rows, const t_int &cols)
Converts from subscript to index for matrix.
Definition: utilities.cc:12
purify::utilities::mod
t_real mod(const t_real &x, const t_real &y)
Mod function modified to wrap circularly for negative numbers.
Definition: utilities.cc:41
purify::widefield::w_support
t_int w_support(const t_real w, const t_real du, const t_int min, const t_int max)
estimate support size of w given u resolution du
Definition: wide_field_utilities.cc:7
purify::widefield::pixel_to_lambda
t_real pixel_to_lambda(const t_real cell, const t_uint imsize, const t_real oversample_ratio)
return factors to convert between arcsecond pixel size image space and lambda for uv space
Definition: wide_field_utilities.cc:11

References purify::projection_kernels::exact_w_projection_integration(), purify::projection_kernels::exact_w_projection_integration_1d(), purify::integration::h, purify::I, purify::utilities::mod(), purify::integration::p, purify::constant::pi, purify::widefield::pixel_to_lambda(), PURIFY_HIGH_LOG, PURIFY_LOW_LOG, purify::utilities::sub2ind(), operators_test::u, operators_test::v, purify::widefield::w_support(), and purify::wkernel_radial.

◆ init_gridding_matrix_2d() [3/5]

template<class STORAGE_INDEX_TYPE = t_int>
Sparse<t_complex, STORAGE_INDEX_TYPE> purify::details::init_gridding_matrix_2d ( const t_uint  number_of_images,
const std::vector< t_int > &  image_index,
const Vector< t_real > &  u,
const Vector< t_real > &  v,
const Vector< t_complex > &  weights,
const t_uint &  imsizey_,
const t_uint &  imsizex_,
const t_real &  oversample_ratio,
const std::function< t_real(t_real)>  kernelu,
const std::function< t_real(t_real)>  kernelv,
const t_uint  Ju,
const t_uint  Jv 
)

Construct all to all gridding matrix.

Definition at line 56 of file operators.h.

61  {
62  if (std::any_of(image_index.begin(), image_index.end(), [&number_of_images](int index) {
63  return index < 0 or index > (number_of_images - 1);
64  }))
65  throw std::runtime_error("Image index is out of bounds");
66  const t_uint ftsizev_ = std::floor(imsizey_ * oversample_ratio);
67  const t_uint ftsizeu_ = std::floor(imsizex_ * oversample_ratio);
68  const t_uint rows = u.size();
69  const STORAGE_INDEX_TYPE cols = ftsizeu_ * ftsizev_ * number_of_images;
70  if (u.size() != v.size())
71  throw std::runtime_error(
72  "Size of u and v vectors are not the same for creating gridding matrix.");
73 
74  Sparse<t_complex, STORAGE_INDEX_TYPE> interpolation_matrix(rows, cols);
75  interpolation_matrix.reserve(Vector<t_int>::Constant(rows, Ju * Jv));
76 
77  const t_complex I(0, 1);
78  const t_int ju_max = std::min(Ju, ftsizeu_);
79  const t_int jv_max = std::min(Jv, ftsizev_);
80  // If I collapse this for loop there is a crash when using MPI... Sparse<>::insert() doesn't work
81  // right
82 #pragma omp parallel for
83  for (t_int m = 0; m < rows; ++m) {
84  for (t_int ju = 1; ju < ju_max + 1; ++ju) {
85  for (t_int jv = 1; jv < jv_max + 1; ++jv) {
86  const t_real k_u = std::floor(u(m) - ju_max * 0.5);
87  const t_real k_v = std::floor(v(m) - jv_max * 0.5);
88  const t_uint q = utilities::mod(k_u + ju, ftsizeu_);
89  const t_uint p = utilities::mod(k_v + jv, ftsizev_);
90  assert(image_index.at(m) < number_of_images);
91  const STORAGE_INDEX_TYPE index =
92  static_cast<STORAGE_INDEX_TYPE>(utilities::sub2ind(p, q, ftsizev_, ftsizeu_)) +
93  static_cast<STORAGE_INDEX_TYPE>(image_index.at(m)) *
94  static_cast<STORAGE_INDEX_TYPE>(ftsizev_ * ftsizeu_);
95  interpolation_matrix.insert(static_cast<STORAGE_INDEX_TYPE>(m), index) =
96  std::exp(-2 * constant::pi * I * ((k_u + ju) * 0.5 + (k_v + jv) * 0.5)) *
97  kernelu(u(m) - (k_u + ju)) * kernelv(v(m) - (k_v + jv)) * weights(m);
98  }
99  }
100  }
101  return interpolation_matrix;
102 }

References purify::I, purify::utilities::mod(), purify::constant::pi, purify::utilities::sub2ind(), operators_test::u, and operators_test::v.

◆ init_gridding_matrix_2d() [4/5]

Sparse< t_complex > purify::details::init_gridding_matrix_2d ( const Vector< t_real > &  u,
const Vector< t_real > &  v,
const Vector< t_complex > &  weights,
const t_uint &  imsizey_,
const t_uint &  imsizex_,
const t_real &  oversample_ratio,
const std::function< t_real(t_real)>  kernelu,
const std::function< t_real(t_real)>  kernelv,
const t_uint  Ju,
const t_uint  Jv 
)

Construct gridding matrix.

Definition at line 7 of file operators.cc.

12  {
13  const t_uint ftsizev_ = std::floor(imsizey_ * oversample_ratio);
14  const t_uint ftsizeu_ = std::floor(imsizex_ * oversample_ratio);
15  const t_uint rows = u.size();
16  const t_uint cols = ftsizeu_ * ftsizev_;
17  if (u.size() != v.size())
18  throw std::runtime_error(
19  "Size of u and v vectors are not the same for creating gridding matrix.");
20 
21  Sparse<t_complex> interpolation_matrix(rows, cols);
22  interpolation_matrix.reserve(Vector<t_int>::Constant(rows, Ju * Jv));
23 
24  const t_complex I(0, 1);
25  const t_int ju_max = std::min(Ju, ftsizeu_);
26  const t_int jv_max = std::min(Jv, ftsizev_);
27  // If I collapse this for loop there is a crash when using MPI... Sparse<>::insert() doesn't work
28  // right
29 #pragma omp parallel for
30  for (t_int m = 0; m < rows; ++m) {
31  for (t_int ju = 1; ju < ju_max + 1; ++ju) {
32  for (t_int jv = 1; jv < jv_max + 1; ++jv) {
33  const t_real k_u = std::floor(u(m) - ju_max * 0.5);
34  const t_real k_v = std::floor(v(m) - jv_max * 0.5);
35  const t_uint q = utilities::mod(k_u + ju, ftsizeu_);
36  const t_uint p = utilities::mod(k_v + jv, ftsizev_);
37  const t_uint index = utilities::sub2ind(p, q, ftsizev_, ftsizeu_);
38  interpolation_matrix.insert(m, index) =
39  std::exp(-2 * constant::pi * I * ((k_u + ju) * 0.5 + (k_v + jv) * 0.5)) *
40  kernelu(u(m) - (k_u + ju)) * kernelv(v(m) - (k_v + jv)) * weights(m);
41  }
42  }
43  }
44  return interpolation_matrix;
45 }

References purify::I, purify::utilities::mod(), purify::constant::pi, purify::utilities::sub2ind(), operators_test::u, and operators_test::v.

Referenced by purify::operators::init_gridding_matrix_2d(), and init_gridding_matrix_2d().

◆ init_gridding_matrix_2d() [5/5]

Sparse< t_complex > purify::details::init_gridding_matrix_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,
const std::function< t_real(t_real)> &  ftkerneluv,
const std::function< t_real(t_real)> &  kerneluv,
const t_uint  Ju,
const t_uint  Jw,
const t_real  cellx,
const t_real  celly,
const t_real  abs_error,
const t_real  rel_error,
const dde_type  dde 
)

Construct gridding matrix with wprojection.

Definition at line 10 of file wproj_operators.cc.

18  {
19  const t_uint ftsizev_ = std::floor(imsizey_ * oversample_ratio);
20  const t_uint ftsizeu_ = std::floor(imsizex_ * oversample_ratio);
21  const t_real du = widefield::pixel_to_lambda(cellx, imsizex_, oversample_ratio);
22  const t_real dv = widefield::pixel_to_lambda(celly, imsizey_, oversample_ratio);
23  PURIFY_HIGH_LOG("du x du: {}, {}", du, dv);
24  if (std::abs(du - dv) > 1e-6)
25  throw std::runtime_error(
26  "Field of view along l and m is not the same, this assumption is required for "
27  "w-projection.");
28  const t_uint rows = u.size();
29  const t_uint cols = ftsizeu_ * ftsizev_;
30  if (u.size() != v.size())
31  throw std::runtime_error(
32  "Size of u and v vectors are not the same for creating gridding matrix.");
33  if (u.size() != w.size())
34  throw std::runtime_error(
35  "Size of u and w vectors are not the same for creating gridding matrix.");
36  if (u.size() != weights.size())
37  throw std::runtime_error(
38  "Size of u and w vectors are not the same for creating gridding matrix.");
39  if (Ju > Jw)
40  throw std::runtime_error(
41  "w kernel size must be at least the size of w=0 kernel, must have Ju <= Jw.");
42  // count gridding coefficients with variable support size
43  Vector<t_int> total_coeffs = Vector<t_int>::Zero(w.size());
44  for (int i = 0; i < w.size(); i++) {
45  const t_int Ju_max =
46  widefield::w_support(w(i), du, static_cast<t_int>(Ju), static_cast<t_int>(Jw));
47  total_coeffs(i) = Ju_max * Ju_max;
48  }
49  const t_real num_of_coeffs = total_coeffs.array().cast<t_real>().sum();
50  PURIFY_HIGH_LOG("Using {} rows (coefficients per a row {}), and memory of {} MegaBytes",
51  total_coeffs.size(), total_coeffs.array().cast<t_real>().mean(),
52  16. * num_of_coeffs / std::pow(10., 6));
53  Sparse<t_complex> interpolation_matrix(rows, cols);
54  try {
55  interpolation_matrix.reserve(total_coeffs);
56  } catch (std::bad_alloc e) {
57  throw std::runtime_error(
58  "Not enough memory for coefficients, choose upper limit on support size Jw.");
59  }
60 
61  const t_complex I(0., 1.);
62 
63  const t_uint max_evaluations = 1e8;
64  const t_real relative_error = rel_error;
65  const t_real absolute_error = abs_error;
66 
67  auto const ftkernel_radial = [&](const t_real l) -> t_real { return ftkerneluv(l); };
68 
69  t_int coeffs_done = 0;
70  t_uint total = 0;
71 
72 #pragma omp parallel for
73  for (t_int m = 0; m < rows; ++m) {
74  // w_projection convolution setup
75  const t_int Ju_max = widefield::w_support(w(m), du, Ju, Jw);
76  t_uint evaluations = 0;
77  const t_int kwu = std::floor(u(m) - Ju_max * 0.5);
78  const t_int kwv = std::floor(v(m) - Ju_max * 0.5);
79  const t_real w_val = w(m); //((0. < w(m)) ? 1: -1) * std::min(Ju_max * du, std::abs(w(m)));
80 
81  for (t_int ju = 1; ju < Ju_max + 1; ++ju) {
82  for (t_int jv = 1; jv < Ju_max + 1; ++jv) {
83  const t_uint q = utilities::mod(kwu + ju, ftsizeu_);
84  const t_uint p = utilities::mod(kwv + jv, ftsizev_);
85  const t_uint index = utilities::sub2ind(p, q, ftsizev_, ftsizeu_);
86  interpolation_matrix.insert(m, index) =
87  std::exp(-2 * constant::pi * I * ((kwu + ju) * 0.5 + (kwv + jv) * 0.5)) * weights(m) *
88  ((dde == dde_type::wkernel_radial)
89  ? projection_kernels::exact_w_projection_integration_1d(
90  (u(m) - (kwu + ju)), (v(m) - (kwv + jv)), w_val, du, oversample_ratio,
91  ftkernel_radial, max_evaluations, absolute_error, relative_error,
92  (du > 1.) ? integration::method::p : integration::method::h, evaluations)
93  : projection_kernels::exact_w_projection_integration(
94  (u(m) - (kwu + ju)), (v(m) - (kwv + jv)), w_val, du, dv, oversample_ratio,
95  ftkernel_radial, ftkernel_radial, max_evaluations, absolute_error,
96  relative_error, integration::method::h, evaluations));
97 #pragma omp critical(add_eval)
98  {
99  coeffs_done++;
100  if (num_of_coeffs > 100) {
101  if ((coeffs_done % (static_cast<t_int>(num_of_coeffs / 100)) == 0)) {
102  PURIFY_LOW_LOG(
103  "w = {}, support = {}x{}, coeffs: {} of {}, {}%", w_val, Ju_max, Ju_max,
104  coeffs_done, num_of_coeffs,
105  static_cast<t_real>(coeffs_done) / static_cast<t_real>(num_of_coeffs) * 100.);
106  }
107  }
108  }
109  }
110  }
111  }
112  interpolation_matrix.makeCompressed();
113  return interpolation_matrix;
114 }

References purify::projection_kernels::exact_w_projection_integration(), purify::projection_kernels::exact_w_projection_integration_1d(), purify::integration::h, purify::I, purify::utilities::mod(), purify::integration::p, purify::constant::pi, purify::widefield::pixel_to_lambda(), PURIFY_HIGH_LOG, PURIFY_LOW_LOG, purify::utilities::sub2ind(), operators_test::u, operators_test::v, purify::widefield::w_support(), and purify::wkernel_radial.

Generated by doxygen 1.9.1