PURIFY
Next-generation radio interferometric imaging
wproj_utilities.h
Go to the documentation of this file.
1 #ifndef PURIFY_wproj_utilities_H
2 #define PURIFY_wproj_utilities_H
3 
4 #include "purify/config.h"
5 #include <fstream>
6 #include <iostream>
7 #ifdef PURIFY_OPENMP
8 #include <omp.h>
9 #endif
10 #include "purify/types.h"
11 #include <random>
12 #include <set>
13 #include <stdio.h>
14 #include <string>
15 #include <time.h>
16 #include <Eigen/Sparse>
17 #include <sys/stat.h>
18 #include "purify/logging.h"
19 #include "purify/operators.h"
20 #include "purify/uvw_utilities.h"
21 
22 namespace purify {
23 namespace utilities {
25 vis_params sort_by_w(const vis_params &uv_data);
27 template <typename T>
28 t_real w_lambert(T x, const t_real &relative_difference);
29 } // namespace utilities
30 namespace wproj_utilities {
32 std::tuple<t_real, t_real> fov_cosines(t_real const &cell_x, t_real const &cell_y,
33  t_uint const &x_size, t_uint const &y_size);
35 template <typename T>
36 Sparse<t_complex> convert_sparse(const Eigen::MatrixBase<T> &M, const t_real &threshold = 0.);
38 Matrix<t_complex> generate_dde(const std::function<t_complex(t_real, t_real)> &dde,
39  const t_real &cell_x, const t_real &cell_y, const t_uint &x_size,
40  const t_uint &y_size);
42 Matrix<t_complex> generate_chirp(const t_real &w_rate, const t_real &cell_x, const t_real &cell_y,
43  const t_uint &x_size, const t_uint &y_size);
44 Matrix<t_complex> generate_chirp(const std::function<t_complex(t_real, t_real)> &dde,
45  const t_real &w_rate, const t_real &cell_x, const t_real &cell_y,
46  const t_uint &x_size, const t_uint &y_size);
48 Sparse<t_complex> create_chirp_row(const t_real &w_rate, const t_real &cell_x, const t_real &cell_y,
49  const t_uint &ftsizev, const t_uint &ftsizeu,
50  const t_real &energy_fraction,
51  const sopt::OperatorFunction<Vector<t_complex>> &fftop);
52 Sparse<t_complex> create_chirp_row(const Vector<t_complex> &chirp_image,
53  const t_real &energy_fraction,
54  const sopt::OperatorFunction<Vector<t_complex>> &fftop);
56 template <typename T>
57 t_real sparsify_row_thres(const Eigen::SparseMatrixBase<T> &row, const t_real &energy);
59 template <typename T>
60 t_real sparsify_row_dense_thres(const Eigen::MatrixBase<T> &row, const t_real &energy);
62 template <class T>
64  const t_uint &x_size, const t_uint &y_size);
66 Sparse<t_complex> wprojection_matrix(const Sparse<t_complex> &G, const t_uint &x_size,
67  const t_uint &y_size, const Vector<t_real> &w_components,
68  const t_real &cell_x, const t_real &cell_y,
69  const t_real &energy_fraction_chirp,
70  const t_real &energy_fraction_wproj,
72  const t_uint order = 1,
73  const t_real &interpolation_error = 1e-2);
75 t_real snr_metric(const Image<t_real> &model, const Image<t_real> &solution);
77 t_real mr_metric(const Image<t_real> &model, const Image<t_real> &solution);
79 t_real sparsity_sp(const Sparse<t_complex> &Gmat);
81 t_real sparsity_im(const Image<t_complex> &Cmat);
83 Sparse<t_complex> generate_vect(const t_uint &x_size, const t_uint &y_size, const t_real &sigma,
84  const t_real &mean);
86 t_real upsample_ratio_sim(const utilities::vis_params &uv_vis, const t_real &L, const t_real &M,
87  const t_int &x_size, const t_int &y_size, const t_int &multipleOf);
88 
89 template <typename T>
90 Sparse<t_complex> convert_sparse(const Eigen::MatrixBase<T> &M, const t_real &threshold) {
91  const Sparse<t_complex> M_sparse =
92  M.sparseView(std::max(M.cwiseAbs().maxCoeff() * 1e-10, threshold), 1);
93  Sparse<t_complex> out(M.rows(), M.cols());
94  out.reserve(M_sparse.nonZeros());
95  for (t_int k = 0; k < M_sparse.outerSize(); ++k)
96  for (Sparse<t_complex>::InnerIterator it(M_sparse, k); it; ++it) {
97  out.coeffRef(it.row(), it.col()) = it.value();
98  }
99  return out;
100 }
101 template <typename T>
102 t_real sparsify_row_thres(const Eigen::SparseMatrixBase<T> &row, const t_real &energy) {
103  /*
104  Takes in a row of G and returns indexes of coeff to keep in the row sparse version
105  energy:: how much energy - in l2 sens - to keep after hard-thresholding
106  */
107 
108  if (energy >= 1) return 0;
109 
110  const Sparse<t_real> row_abs = row.cwiseAbs();
111 
112  t_real abs_row_max = 0;
113  for (t_uint i = 0; i < row_abs.nonZeros(); i++) {
114  if (*(row_abs.valuePtr() + i) > abs_row_max) abs_row_max = *(row_abs.valuePtr() + i);
115  }
116  const t_real row_total_energy = row_abs.cwiseProduct(row_abs).sum();
117  const t_real target_threshold_energy_sum = row_total_energy * energy;
118  const t_int niters = 200;
119  t_real lower = 0;
120  t_real upper = abs_row_max;
121  t_real current_threshold = upper * 0.5;
122  t_real old_threshold = upper;
123  bool converged = false;
124  for (t_int i = 0; i < niters; i++) {
125  t_real threshold_energy_sum = 0;
126  for (t_int j = 0; j < row_abs.nonZeros(); j++) {
127  if (*(row_abs.valuePtr() + j) > current_threshold)
128  threshold_energy_sum += *(row_abs.valuePtr() + j) * *(row_abs.valuePtr() + j);
129  }
130  if ((std::abs(current_threshold - old_threshold) / std::abs(old_threshold) < 1e-3) and
131  ((threshold_energy_sum >= target_threshold_energy_sum) and
132  (((threshold_energy_sum / target_threshold_energy_sum) - 1) < 0.01))) {
133  converged = true;
134  break;
135  }
136  if (threshold_energy_sum < target_threshold_energy_sum)
137  upper = current_threshold;
138  else
139  lower = current_threshold;
140  old_threshold = current_threshold;
141  current_threshold = (lower + upper) * 0.5;
142  }
143  if (!converged) current_threshold = lower;
144  return current_threshold;
145 }
146 template <typename T>
147 t_real sparsify_row_dense_thres(const Eigen::MatrixBase<T> &row, const t_real &energy) {
148  /*
149  Takes in a row of G and returns indexes of coeff to keep in the row sparse version
150  energy:: how much energy - in l2 sens - to keep after hard-thresholding
151  */
152 
153  if (energy >= 1) return 0;
154 
155  const Vector<t_real> row_abs = row.cwiseAbs();
156  const t_real abs_row_max = row.cwiseAbs().maxCoeff();
157  const t_real row_total_energy = row_abs.cwiseProduct(row_abs).sum();
158  const t_real target_threshold_energy_sum = row_total_energy * energy;
159  const t_int niters = 200;
160  t_real lower = 0;
161  t_real upper = abs_row_max;
162  t_real current_threshold = upper * 0.5;
163  t_real old_threshold = upper;
164  bool converged = false;
165  for (t_int i = 0; i < niters; i++) {
166  t_real threshold_energy_sum = 0;
167  for (t_int j = 0; j < row_abs.size(); j++) {
168  if (row_abs(j) > current_threshold) threshold_energy_sum += row_abs(j) * row_abs(j);
169  }
170  if ((std::abs(current_threshold - old_threshold) / std::abs(old_threshold) < 1e-3) and
171  ((threshold_energy_sum >= target_threshold_energy_sum) and
172  (((threshold_energy_sum / target_threshold_energy_sum) - 1) < 0.01))) {
173  converged = true;
174  break;
175  }
176  if (threshold_energy_sum < target_threshold_energy_sum)
177  upper = current_threshold;
178  else
179  lower = current_threshold;
180  old_threshold = current_threshold;
181  current_threshold = (lower + upper) * 0.5;
182  }
183  if (!converged) current_threshold = lower;
184  return current_threshold;
185 }
186 template <class T>
188  const t_uint &x_size, const t_uint &y_size) {
189  assert((Grid_.nonZeros() > 0) and (chirp_.nonZeros() > 0));
190  Matrix<t_complex> output_row = Matrix<t_complex>::Zero(1, x_size * y_size);
191  const t_int x_sizec = std::floor(x_size * 0.5);
192  const t_int y_sizec = std::floor(y_size * 0.5);
193  for (t_int k = 0; k < output_row.size(); k++) {
194  // shifting indicies to workout linear convolution
195  t_int i_shift_W = 0;
196  t_int j_shift_W = 0;
197  std::tie(j_shift_W, i_shift_W) = utilities::ind2sub(k, y_size, x_size);
198  const t_int i_W = utilities::mod(i_shift_W + x_sizec, x_size);
199  const t_int j_W = utilities::mod(j_shift_W + y_sizec, y_size);
200  for (Sparse<t_complex>::InnerIterator it(Grid_, 0); it; ++it) {
201  t_int i_shift_G = 0;
202  t_int j_shift_G = 0;
203  std::tie(j_shift_G, i_shift_G) = utilities::ind2sub(it.index(), y_size, x_size);
204  const t_int i_G = utilities::mod(i_shift_G + x_sizec, x_size);
205  const t_int j_G = utilities::mod(j_shift_G + y_sizec, y_size);
206  const t_int i_C = i_W - i_G;
207  const t_int j_C = j_W - j_G;
208  // Checking that convolution result is going to be within bounds
209  if (-x_sizec <= i_C and i_C < std::ceil(x_size * 0.5) and -y_sizec <= j_C and
210  j_C < std::ceil(y_size * 0.5)) {
211  // FFTshift of result to go into gridding matrix
212  const t_int i_shift_C = utilities::mod(i_C, x_size);
213  const t_int j_shift_C = utilities::mod(j_C, y_size);
214  const t_int pos = utilities::sub2ind(j_shift_C, i_shift_C, y_size, x_size);
215  assert(0 <= pos and pos < x_size * y_size);
216  output_row(0, k) += it.value() * chirp_.coeff(0, pos);
217  }
218  }
219  }
220  return convert_sparse(output_row);
221 }
222 template <class T>
224  const Sparse<T> &chirp_, const t_uint &x_size,
225  const t_uint &y_size) {
226  assert((Grid_.nonZeros() > 0) and (chirp_.nonZeros() > 0));
227  const t_int x_sizec = std::floor(x_size * 0.5);
228  const t_int y_sizec = std::floor(y_size * 0.5);
229  std::set<t_int> non_zero_W;
230  for (Sparse<t_complex>::InnerIterator jt(chirp_, 0); jt; ++jt) {
231  // shifting indicies to workout linear convolution
232  t_int i_shift_C = 0;
233  t_int j_shift_C = 0;
234  std::tie(j_shift_C, i_shift_C) = utilities::ind2sub(jt.index(), y_size, x_size);
235  const t_int i_C = utilities::mod(i_shift_C + x_sizec, x_size) - x_sizec;
236  const t_int j_C = utilities::mod(j_shift_C + y_sizec, y_size) - y_sizec;
237  for (Sparse<t_complex>::InnerIterator it(Grid_, 0); it; ++it) {
238  t_int i_shift_G = 0;
239  t_int j_shift_G = 0;
240  std::tie(j_shift_G, i_shift_G) = utilities::ind2sub(it.index(), y_size, x_size);
241  const t_int i_G = utilities::mod(i_shift_G + x_sizec, x_size) - x_sizec;
242  const t_int j_G = utilities::mod(j_shift_G + y_sizec, y_size) - y_sizec;
243  const t_int i_W = i_G + i_C;
244  const t_int j_W = j_G + j_C;
245  if (-x_sizec <= i_W and i_W < std::ceil(x_size * 0.5) and -y_sizec <= j_W and
246  j_W < std::ceil(y_size * 0.5)) {
247  const t_int i_shift_W = utilities::mod(i_W, x_size);
248  const t_int j_shift_W = utilities::mod(j_W, y_size);
249  const t_real pos = utilities::sub2ind(j_shift_W, i_shift_W, y_size, x_size);
250  non_zero_W.insert(pos);
251  }
252  }
253  }
254  if (non_zero_W.size() < std::min(Grid_.nonZeros(), chirp_.nonZeros()))
255  throw std::runtime_error("Gridding kernel or chirp kernel is zero.");
256  Sparse<t_complex> output_row(1, x_size * y_size);
257  output_row.reserve(non_zero_W.size());
258 
259  for (t_int k = 0; k < output_row.size(); k++) {
260  // shifting indicies to workout linear convolution
261  t_int i_shift_W = 0;
262  t_int j_shift_W = 0;
263  std::tie(j_shift_W, i_shift_W) = utilities::ind2sub(k, y_size, x_size);
264  const t_int i_W = utilities::mod(i_shift_W + x_sizec, x_size);
265  const t_int j_W = utilities::mod(j_shift_W + y_sizec, y_size);
266  for (Sparse<t_complex>::InnerIterator it(Grid_, 0); it; ++it) {
267  t_int i_shift_G = 0;
268  t_int j_shift_G = 0;
269  std::tie(j_shift_G, i_shift_G) = utilities::ind2sub(it.index(), y_size, x_size);
270  const t_int i_G = utilities::mod(i_shift_G + x_sizec, x_size);
271  const t_int j_G = utilities::mod(j_shift_G + y_sizec, y_size);
272  const t_int i_C = i_W - i_G;
273  const t_int j_C = j_W - j_G;
274  // Checking that convolution result is going to be within bounds
275  if (-x_sizec <= i_C and i_C < std::ceil(x_size * 0.5) and -y_sizec <= j_C and
276  j_C < std::ceil(y_size * 0.5)) {
277  // FFTshift of result to go into gridding matrix
278  const t_int i_shift_C = utilities::mod(i_C, x_size);
279  const t_int j_shift_C = utilities::mod(j_C, y_size);
280  const t_int pos = utilities::sub2ind(j_shift_C, i_shift_C, y_size, x_size);
281  assert(0 <= pos and pos < x_size * y_size);
282  output_row.coeffRef(0, k) += it.value() * chirp_.coeff(0, pos);
283  }
284  }
285  }
286  return output_row;
287 }
288 } // namespace wproj_utilities
289 namespace utilities {
290 
291 template <typename T>
292 t_real w_lambert(const T &x, const t_real &relative_difference = 1e-8) {
293  if (std::is_integral<T>::value) return w_lambert<t_real>(x, relative_difference);
294  T z = T(x);
295  T diff = T(1e4 * relative_difference);
296  if (x < -1 / std::exp(1)) throw std::runtime_error("Out of bounds for W lambert function!");
297  while (std::abs(diff) > relative_difference) {
298  const T f = z * std::exp(z) - x;
299  const T w = z - f / (std::exp(z) * (z + 1) - f * (z + 2) / (2 * z + 2));
300  diff = w - z;
301  z = w;
302  }
303  return z;
304 }
305 
306 } // namespace utilities
307 } // namespace purify
308 
309 #endif
t_real w_lambert(T x, const t_real &relative_difference)
Calculate W Lambert function.
K::Scalar mean(const K x)
Calculate mean of vector.
Definition: utilities.h:21
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
t_real mod(const t_real &x, const t_real &y)
Mod function modified to wrap circularly for negative numbers.
Definition: utilities.cc:41
vis_params sort_by_w(const vis_params &uv_data)
Sort visibilities to be from w_max to w_min.
std::tuple< t_int, t_int > ind2sub(const t_int &sub, const t_int &cols, const t_int &rows)
Converts from index to subscript for matrix.
Definition: utilities.cc:26
series
Type of series approximation.
Definition: types.h:64
t_real snr_metric(const Image< t_real > &model, const Image< t_real > &solution)
SNR calculation.
t_real sparsity_sp(const Sparse< t_complex > &Gmat)
return fraction of non zero values from sparse matrix
Matrix< t_complex > generate_chirp(const std::function< t_complex(t_real, t_real)> &dde, 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)
t_real sparsify_row_thres(const Eigen::SparseMatrixBase< T > &row, const t_real &energy)
Returns threshold to keep a fraction of energy in the sparse row.
std::tuple< t_real, t_real > fov_cosines(t_real const &cell_x, t_real const &cell_y, t_uint const &x_size, t_uint const &y_size)
Work out max L and M directional cosines from image parameters.
t_real sparsify_row_dense_thres(const Eigen::MatrixBase< T > &row, const t_real &energy)
Returns threshold to keep a fraction of energy in the dense row.
Sparse< t_complex > create_chirp_row(const t_real &w_rate, const t_real &cell_x, const t_real &cell_y, const t_uint &ftsizev, const t_uint &ftsizeu, const t_real &energy_fraction, const sopt::OperatorFunction< Vector< t_complex >> &fftop)
Generates row of chirp matrix from image of chirp.
Sparse< t_complex > row_wise_convolution(const Sparse< t_complex > &Grid_, const Sparse< T > &chirp_, const t_uint &x_size, const t_uint &y_size)
Perform convolution with gridding matrix row and chirp matrix row.
Sparse< t_complex > convert_sparse(const Eigen::MatrixBase< T > &M, const t_real &threshold=0.)
Convert from dense to sparse.
t_real sparsity_im(const Image< t_complex > &Cmat)
return faction of non zero values from matrix
Matrix< t_complex > generate_dde(const std::function< t_complex(t_real, t_real)> &dde, const t_real &cell_x, const t_real &cell_y, const t_uint &x_size, const t_uint &y_size)
Generate image of DDE for A or W projection.
Sparse< t_complex > generate_vect(const t_uint &x_size, const t_uint &y_size, const t_real &sigma, const t_real &mean)
Genereates an image of a Gaussian as a sparse matrice.
t_real upsample_ratio_sim(const utilities::vis_params &uv_vis, const t_real &L, const t_real &M, const t_int &x_size, const t_int &y_size, const t_int &multipleOf)
Calculate upsample ratio from bandwidth (only needed for simulations)
t_real mr_metric(const Image< t_real > &model, const Image< t_real > &solution)
MR calculation.
Sparse< t_complex > row_wise_sparse_convolution(const Sparse< t_complex > &Grid_, const Sparse< T > &chirp_, const t_uint &x_size, const t_uint &y_size)
Sparse< t_complex > wprojection_matrix(const Sparse< t_complex > &G, const t_uint &x_size, const t_uint &y_size, const Vector< t_real > &w_components, const t_real &cell_x, const t_real &cell_y, const t_real &energy_fraction_chirp, const t_real &energy_fraction_wproj, const expansions::series series, const t_uint order, const t_real &interpolation_error)
Produce Gridding matrix convovled with chirp matrix for wprojection.
Eigen::SparseMatrix< T, Eigen::RowMajor, I > Sparse
A matrix of a given type.
Definition: types.h:40