Pixelated convergence interpolation

[AlexandreAdam]

Background

Deflection angles are computed from the convergence using the following integral

$$\boldsymbol{\alpha}(\boldsymbol{\theta}) = \frac{1}{\pi}\int_{\mathbb{R}^2}\, \kappa(\boldsymbol{\theta}') G(\boldsymbol{\theta} - \boldsymbol{\theta}') d\boldsymbol{\theta}'$$

where the kernel for the convolution is

$$G(\boldsymbol{\theta} - \boldsymbol{\theta}') = \frac{\boldsymbol{\theta} - \boldsymbol{\theta}'}{\lVert \boldsymbol{\theta} - \boldsymbol{\theta}' \rVert^2}$$

In general, we would need to solve this integral for every point $\boldsymbol{\theta}$ requested by a user. In order to speed up this computation for the pixelated convergence $\kappa_{ij}$, defined on a grid $\{\boldsymbol{\theta}_{ij}'\}$, with $i,j \in \{0,\dots, 2N-1 \}$ and $\boldsymbol{\theta}_{00}' = (0, 0)$, we can instead precompute a convolution kernel matrix as follows

$$G_{ij} = \frac{- \boldsymbol{\theta}_{ij}'}{\lVert \boldsymbol{\theta}_{ij}' \rVert^2}$$

Note that the elements $\kappa_{ij}$ is only defined on a grid with $i = \{0, \dots, N-1\}$. We apply circular boundary condition to extend it to $i = \{0, \dots, 2N-1\}$.

This kernel is particularly useful when $\boldsymbol{\theta}$ lies exactly on our coordinate grid, i.e. $\boldsymbol{\theta}_{\ell k} \in \{\boldsymbol{\theta}_{ij}\}$. In that case, we can approximate the convolution by a discrete circular convolution as follows

$$\boldsymbol{\alpha}(\boldsymbol{\theta}_{\ell k}) = \frac{(\Delta \theta)^2}{\pi}\sum_{i = 0}^{2N-1} \sum_{j = 0}^{2N-1} \kappa_{ij} G_{\ell - i, k - j}$$

Using the convolution theorem, we can evaluate this in linear time with the FFT algorithm. Denoting the Fourier transform as $\mathcal{F}$, and considering only one component of the deflection field (e.g. the $x$-component), we have

$$\alpha_x(\boldsymbol{\theta}_{\ell k}) = \frac{(\Delta \theta)^2}{\pi} \mathcal{F}^{-1} \big[\mathcal{F}[\kappa] \odot \mathcal{F}[G_x] \big]_{\ell k}$$

where $\odot$ denotes the Hadamard product (element-wise product).

Issue 1

Suppose that $\boldsymbol{\theta}_{\ell k}$ does not lie on our convolution grid. How should we evaluate the deflection angle for the pixelated convergence? In the code, we use bilinear interpolation which does not seem right and relies on strong assumptions, namely

1. The point $\boldsymbol{\theta}$ lies within the convolution grid or close enough to its edge
2. The convolution grid is spaced well enough to accurately reconstruct the underlying deflection field (which might have a Nyquist frequency smaller than the pixelated convergence)

The PixelatedConvergence needs more thorough testing to validate these assumptions, and perhaps needs an additional method to handle more exotic use cases.

Issue 2

We need to include an interpolation of the convergence map to decouple the convolution grid from the coordinates of the pixels in the pixelated convergence model.