1 #ifndef PURIFY_wproj_utilities_H
2 #define PURIFY_wproj_utilities_H
4 #include "purify/config.h"
16 #include <Eigen/Sparse>
25 vis_params
sort_by_w(
const vis_params &uv_data);
28 t_real
w_lambert(T x,
const t_real &relative_difference);
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);
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);
49 const t_uint &ftsizev,
const t_uint &ftsizeu,
50 const t_real &energy_fraction,
51 const sopt::OperatorFunction<Vector<t_complex>> &fftop);
53 const t_real &energy_fraction,
54 const sopt::OperatorFunction<Vector<t_complex>> &fftop);
57 t_real
sparsify_row_thres(
const Eigen::SparseMatrixBase<T> &row,
const t_real &energy);
64 const t_uint &x_size,
const t_uint &y_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);
87 const t_int &x_size,
const t_int &y_size,
const t_int &multipleOf);
92 M.sparseView(std::max(M.cwiseAbs().maxCoeff() * 1e-10, threshold), 1);
94 out.reserve(M_sparse.nonZeros());
95 for (t_int k = 0; k < M_sparse.outerSize(); ++k)
97 out.coeffRef(it.row(), it.col()) = it.value();
101 template <
typename T>
108 if (energy >= 1)
return 0;
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);
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;
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);
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))) {
136 if (threshold_energy_sum < target_threshold_energy_sum)
137 upper = current_threshold;
139 lower = current_threshold;
140 old_threshold = current_threshold;
141 current_threshold = (lower + upper) * 0.5;
143 if (!converged) current_threshold = lower;
144 return current_threshold;
146 template <
typename T>
153 if (energy >= 1)
return 0;
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;
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);
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))) {
176 if (threshold_energy_sum < target_threshold_energy_sum)
177 upper = current_threshold;
179 lower = current_threshold;
180 old_threshold = current_threshold;
181 current_threshold = (lower + upper) * 0.5;
183 if (!converged) current_threshold = lower;
184 return current_threshold;
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++) {
206 const t_int i_C = i_W - i_G;
207 const t_int j_C = j_W - j_G;
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)) {
215 assert(0 <= pos and pos < x_size * y_size);
216 output_row(0, k) += it.value() * chirp_.coeff(0, pos);
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;
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;
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)) {
250 non_zero_W.insert(pos);
254 if (non_zero_W.size() < std::min(Grid_.nonZeros(), chirp_.nonZeros()))
255 throw std::runtime_error(
"Gridding kernel or chirp kernel is zero.");
257 output_row.reserve(non_zero_W.size());
259 for (t_int k = 0; k < output_row.size(); k++) {
272 const t_int i_C = i_W - i_G;
273 const t_int j_C = j_W - j_G;
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)) {
281 assert(0 <= pos and pos < x_size * y_size);
282 output_row.coeffRef(0, k) += it.value() * chirp_.coeff(0, pos);
289 namespace utilities {
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);
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));
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.
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.
t_real mod(const t_real &x, const t_real &y)
Mod function modified to wrap circularly for negative numbers.
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.
series
Type of series approximation.
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.