-
Notifications
You must be signed in to change notification settings - Fork 13
Operators in the measurement operator
This page explains how the measurement operators designed, in terms of each operation. The base operator and its adjoint are constructed in terms of its linear operations. This involves a scaling and zero padding, an FFT, and a sparse matrix that interpolates FFT coefficients off the grid to the measurement locations.
The code can be found in cpp/purify/operators.h.
Below is an example of how each operator is constructed. The functions are defined in
314 sopt::OperatorFunction<T> directZ, indirectZ; //Zero padding operator
315 sopt::OperatorFunction<T> directFFT, indirectFFT; // FFT operator
316 sopt::OperatorFunction<T> directG, indirectG; // Degridding matrixThese functions are lambda functions, that have input and return type T. Typically, T = Vector<t_complex>.
They are constructed in the function
300 template <class T>
301 std::tuple<sopt::OperatorFunction<T>, sopt::OperatorFunction<T>>
302 base_degrid_operator_2d(const Vector<t_real> &u, const Vector<t_real> &v, const Vector<t_real> &w,
303 const Vector<t_complex> &weights, const t_uint &imsizey,
304 const t_uint &imsizex, const t_real &oversample_ratio = 2,
305 const std::string &kernel = "kb", const t_uint Ju = 4, const t_uint Jv = 4,
306 const std::string &ft_plan = "measure", const bool w_term = false,
307 const t_real &cellx = 1, const t_real &celly = 1,
308 const t_real &energy_chirp_fraction = 1,
309 const t_real &energy_kernel_fraction = 1) {
310
311 std::function<t_real(t_real)> kernelu, kernelv, ftkernelu, ftkernelv;
312 std::tie(kernelu, kernelv, ftkernelu, ftkernelv)
313 = purify::create_kernels(kernel, Ju, Jv, imsizey, imsizex, oversample_ratio);
314 sopt::OperatorFunction<T> directZ, indirectZ;
315 sopt::OperatorFunction<T> directFFT, indirectFFT;
316 sopt::OperatorFunction<T> directG, indirectG;
317 const Image<t_real> S = purify::details::init_correction2d(oversample_ratio, imsizey, imsizex,
318 ftkernelu, ftkernelv);
319 PURIFY_LOW_LOG("Building Measurement Operator: WGFZDB");
320 PURIFY_LOW_LOG("Constructing Zero Padding and Correction Operator: ZDB");
321 PURIFY_MEDIUM_LOG("Image size (width, height): {} x {}", imsizex, imsizey);
322 PURIFY_MEDIUM_LOG("Oversampling Factor: {}", oversample_ratio);
323 PURIFY_MEDIUM_LOG("FoV (width, height): {} deg x {} deg", imsizex * cellx / (60. * 60.),
324 imsizey * celly / (60. * 60.));
325 std::tie(directZ, indirectZ) = purify::operators::init_zero_padding_2d<T>(S, oversample_ratio);
326 PURIFY_LOW_LOG("Constructing FFT operator: F");
327 if(ft_plan == "measure") {
328 PURIFY_LOW_LOG("Measuring...");
329 } else {
330 PURIFY_LOW_LOG("Using an estimate");
331 }
332 std::tie(directFFT, indirectFFT)
333 = purify::operators::init_FFT_2d<T>(imsizey, imsizex, oversample_ratio, ft_plan);
334 PURIFY_LOW_LOG("Constructing Weighting and Gridding Operators: WG");
335 PURIFY_MEDIUM_LOG("Number of visibilities: {}", u.size());
336 PURIFY_MEDIUM_LOG("Kernel Support: {} x {}", Ju, Jv);
337 std::tie(directG, indirectG) = purify::operators::init_gridding_matrix_2d<T>(
338 u, v, w, weights, imsizey, imsizex, oversample_ratio, kernelv, kernelu, Ju, Jv, w_term, cellx,
339 celly, energy_chirp_fraction, energy_kernel_fraction);
340 auto direct = sopt::chained_operators<T>(directG, directFFT, directZ);
341 auto indirect = sopt::chained_operators<T>(indirectZ, indirectFFT, indirectG);
342 return std::make_tuple(direct, indirect);
343 }The inputs for this function are the same as for the measurement operator. This function returns the direct and indirect linear transforms of the measurement operator.
After each function is constructed, the direct and indirect functions are chained and returned
340 auto direct = sopt::chained_operators<T>(directG, directFFT, directZ);
341 auto indirect = sopt::chained_operators<T>(indirectZ, indirectFFT, indirectG);
342 return std::make_tuple(direct, indirect);The function sopt::chained_operators<T> is a function that will take lambda functions of type input and return the composition of these functions. For example, direct = sopt::chained_operators<T>(A, B, C, D) will return the function such that direct(x) == A(B(C(D(x)))). It is expected that the input type of D is T and the output type of A is T.
The zero padding operator needs to pad an image vector with zeros, and scale the image to correct for the interpolation off the grid (done later in the measurement operator).
First we need to generate the interpolation kernels and their Fourier transforms, these are lambda functions.
311 std::function<t_real(t_real)> kernelu, kernelv, ftkernelu, ftkernelv;
312 std::tie(kernelu, kernelv, ftkernelu, ftkernelv)
313 = purify::create_kernels(kernel, Ju, Jv, imsizey, imsizex, oversample_ratio);An S image is then generated, which is known as the gridding correction.
317 const Image<t_real> S = purify::details::init_correction2d(oversample_ratio, imsizey, imsizex,
318 ftkernelu, ftkernelv);The zero padding function is then created with
325 std::tie(directZ, indirectZ) = purify::operators::init_zero_padding_2d<T>(S, oversample_ratio);When you look inside this function below, you can see that it constructs two lambda functions that will scale an image by S, and zero pad (direct) or cutout (indirect) the image.
211 //! Construsts zero padding operator
212 template <class T>
213 std::tuple<sopt::OperatorFunction<T>, sopt::OperatorFunction<T>>
214 init_zero_padding_2d(const Image<t_real> &S, const t_real &oversample_ratio) {
215 const t_uint imsizex_ = S.cols();
216 const t_uint imsizey_ = S.rows();
217 const t_uint ftsizeu_ = std::floor(imsizex_ * oversample_ratio);
218 const t_uint ftsizev_ = std::floor(imsizey_ * oversample_ratio);
219 const t_uint x_start = std::floor(ftsizeu_ * 0.5 - imsizex_ * 0.5);
220 const t_uint y_start = std::floor(ftsizev_ * 0.5 - imsizey_ * 0.5);
221 auto direct = [=](T &output, const T &x) {
222 assert(x.size() == imsizex_ * imsizey_);
223 output = Vector<t_complex>::Zero(ftsizeu_ * ftsizev_);
224 #pragma omp simd collapse(2)
225 for(t_uint j = 0; j < imsizey_; j++) {
226 for(t_uint i = 0; i < imsizex_; i++) {
227 const t_uint input_index = utilities::sub2ind(j, i, imsizey_, imsizex_);
228 const t_uint output_index
229 = utilities::sub2ind(y_start + j, x_start + i, ftsizev_, ftsizeu_);
230 output(output_index) = S(j, i) * x(input_index);
231 }
232 }
233 };
234 auto indirect = [=](T &output, const T &x) {
235 assert(x.size() == ftsizeu_ * ftsizev_);
236 output = T::Zero(imsizey_ * imsizex_);
237 #pragma omp simd collapse(2)
238 for(t_uint j = 0; j < imsizey_; j++) {
239 for(t_uint i = 0; i < imsizex_; i++) {
240 const t_uint output_index = utilities::sub2ind(j, i, imsizey_, imsizex_);
241 const t_uint input_index = utilities::sub2ind(y_start + j, x_start + i, ftsizev_, ftsizeu_);
242 output(output_index) = S(j, i) * x(input_index);
243 }
244 }
245 };
246 return std::make_tuple(direct, indirect);
247 }