Notes on measurement operator

Measurement Operator

The measurement operator maps between the model of the sky and what the telescope would measure. For a radio interferometric telescope, this is a model of the radio sky and its Fourier coefficients. Ideally, the model of the sky could be Fourier transformed directly to get the coefficients, but this is not computationally possible for large images and many coefficients. As an approximation, an FFT is used, with interpolation of coefficients off the FFT grid as an approximation.

The code

template <class T>
    sopt::LinearTransform<T>
    init_degrid_operator_2d(const utilities::vis_params &uv_vis_input, const t_uint &imsizey,
                            const t_uint &imsizex, const t_real &cell_x, const t_real &cell_y,
                            const t_real &oversample_ratio = 2, const t_uint &power_iters = 100,
                            const t_real &power_tol = 1e-4, const std::string &kernel = "kb",
                            const t_uint Ju = 4, const t_uint Jv = 4,
                            const std::string &ft_plan = "measure", const bool w_term = false,
                            const t_real &energy_chirp_fraction = 1,
                            const t_real &energy_kernel_fraction = 1);

The returned object is a sopt::LinearTransform<T> operator. This represents a matrix of dimensions (M, N), where N = imsizey * imsizex is the size of the image vector and M = uv_vis_input.u.size() is the size of the measurements vector. Typically, T = Vector<t_complex>, and the operator can be applied to Vector<t_complex>

const sopt::LinearTransform<Vector<t_complex>> measure_op = measurementoperator::init_degrid_operator_2d(uv_vis_input, imsizey, imsizex, cell_y, cell_x);

const auto x = Vector<t_complex>::Random(imsizey * imsizex);
const auto y = measure_op * x;

const auto y_indirect = Vector<t_complex>::Random(uv_vis_input.u.size());
const auto x_indirect = measure_op.adjoint() * y_indirect;

where we have applied the linear transform measure_op and its adjoint measure_op.adjoint().

Input

The vectors uv_vis_input.u, uv_vis_input.v, uv_vis_input.w represent the coordinates in the Fourier domain. For small field of view, we assume that uv_vis_input.w=0, and the coordinates are on a 2D (u, v) plane.

The uv_vis_input.weights are the weights applied to each measurement. This would be the inverse of the covariance matrix of the measurements, which we assume to be diagonal.

imsizey and imsizex are the height and width of the model image in pixels.

cell_y and cell_x are the height and width of a pixel in arcseconds. This is only important if w_term = true or if uv_vis_input.units = lambda. If they are 0, then an automatic value can be calculated from the maximum (u,v) vector.

oversample_ratio is the ratio to zero pad the image before the FFT, to increase interpolation accuracy.

power_iters is the number of iterations needed for the power method, used to estimate the normalisation of the operator.

power_tol is the relative difference convergence criteria for the power method.

kernel is the type of interpolation kernel to use in degridding.

Ju, Jv is the support size of the interpolation kernel.

ft_plan is the type of planning done by FFTW, and can be either measure or estimate.

w_term can be true or false. It will tell the measurement operator to perform corrections to the interpolation kernel for when uv_vis_input.w != 0, using what is known as the w-projection. This will slow down the measurement operator construction greatly.

energy_chirp_fraction only applies when w_term = true, and should speed up computation. It says how much energy the w_term correction should keep.

energy_kernel_fraction only applies when w_term = true, and should speed up computation. It says how much energy the final kernel should keep after applying the w_term correction.

MPI variations of the Measurement Operator

template <class T>
    sopt::LinearTransform<T>
    init_degrid_operator_2d(const sopt::mpi::Communicator &comm,
                            const utilities::vis_params &uv_vis_input, const t_uint &imsizey,
                            const t_uint &imsizex, const t_real &cell_x, const t_real &cell_y,
                            const t_real &oversample_ratio = 2, const t_uint &power_iters = 100,
                            const t_real &power_tol = 1e-4, const std::string &kernel = "kb",
                            const t_uint Ju = 4, const t_uint Jv = 4,
                            const std::string &ft_plan = "measure", const bool w_term = false,
                            const t_real &energy_chirp_fraction = 1,
                            const t_real &energy_kernel_fraction = 1);

This is the same as the non-MPI measurement operator, but the adjoint will perform an comm.all_sum_all of the output image. This means the image is distributed. This is used in Algorithm 3.

template <class T>
    sopt::LinearTransform<T>
    init_degrid_operator_2d_mpi(const sopt::mpi::Communicator &comm,
                            const utilities::vis_params &uv_vis_input, const t_uint &imsizey,
                            const t_uint &imsizex, const t_real &cell_x, const t_real &cell_y,
                            const t_real &oversample_ratio = 2, const t_uint &power_iters = 100,
                            const t_real &power_tol = 1e-4, const std::string &kernel = "kb",
                            const t_uint Ju = 4, const t_uint Jv = 4,
                            const std::string &ft_plan = "measure", const bool w_term = false,
                            const t_real &energy_chirp_fraction = 1,
                            const t_real &energy_kernel_fraction = 1);

This has a different structure from the non-MPI measurement operator. The FFT is performed on the root node, then the FFT grid is scattered to nodes for interpolation to their measurements. The adjoint will interpolate onto the FT grid on each node, then the grids are gathered to the root node. The grid is then IFFT'd, and broadcast to all other nodes. This means the image is distributed, as well as parts of the FT grid. This is used in Algorithm 1.