The idea of convolutional gridding is relatively simple. The raw data samples (which are located at irregular coordinates) are convolved with a kernel function, and compute the result on the center coordinates of the desired regular grid cells. If a Gaussian is used as kernel, one can safe a lot of computing time, because each raw data sample will only influence the output grid cells in immediate vicinity.
Mathematically, this approach is describe with the following formula:
$$R_{i,j}[s]=\frac{1}{W_{i,j}}\sum_n R_n[s](\alpha_n,\delta_n)w(\alpha_{i,j},\delta_{i,j};\alpha_n,\delta_n)$$where:
$$W_{i,j}\equiv\sum_n w(\alpha_{i,j},\delta_{i,j};\alpha_n,\delta_n),$$is called the weight map.
Here, $R_n [s]$ and $R_{i,j}[s]$ are two different representations of the true signal s. The index n runs over the list of all input samples, with respective coordinates $(\alpha_n,\delta_n)$, while the regular output grid can be parametrized via pixel coordinates $(i,j)$ with associated world coordinates $(\alpha_{i,j}, \delta_{i,j})$. The value of the weighting kernel w depends only on the input and output coordinates. In most cases a radially symmetric kernel is applied, such that:
$$w(\alpha_{i,j},\delta_{i,j};\alpha_n,\delta_n)=w\left(\mathcal{D}(\alpha_{i,j},\delta_{i,j};\alpha_n,\delta_n)\right)$$with the distance $\mathcal{D}$ between a pixel center world coordinate, $(\alpha_{i,j}, \delta_{i,j})$, and the $n$-th input coordinate, $(\alpha_n,\delta_n)$.
Since once usually doesn't want to keep all raw data samples in memory, we use two helper maps. While we iterate over the samples (sum over $n$), we compute the products $R_n[s](\alpha_n,\delta_n)w(\alpha_{i,j},\delta_{i,j};\alpha_n,\delta_n)$ for all pixels $(i, j)$ within a sphere around the $(\alpha_n,\delta_n)$ position, and add them to an empty array A[i,j], which can be identified with $R_{i,j}[s]\cdot W_{i,j}$. At the same time, we do the same for the weighting factors, only, and add these to another array B[i,j] aka the weight map. When all raw data samples have been processed, the two resulting arrays are divided (A[i,j] / B[i,j]), which gives $R_{i,j}$.
(We note that cygrid in fact follows a slightly different approach, caching some intermediary values, to allow multi-threaded processing. For details, please see the Paper.)
We demonstrate the approach described above with a little animation.
In [1]:
%load_ext autoreload
%autoreload 2
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
In [2]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation, rc
from astropy.wcs import WCS
import cygrid
np.seterr(divide='ignore', invalid='ignore')
# may need to install ffmpeg for video rendering:
# conda install -c conda-forge ffmpeg
rc('animation', html='html5')
For our first demonstration, we will sample a pre-defined function at random position (for longitude and latitude between -1 and 1) and grid to a regular map. While we do this, we display intermediary results.
In [3]:
def func(l, b):
return l * (1 - l) * np.cos(4 * np.pi * l) * np.sin(4 * np.pi * b ** 2) ** 2
In [4]:
header = {
'NAXIS': 2,
'NAXIS1': 101,
'NAXIS2': 101,
'CTYPE1': 'GLON-CAR',
'CTYPE2': 'GLAT-CAR',
'CUNIT1': 'deg',
'CUNIT2': 'deg',
'CDELT1': -2 / 101,
'CDELT2': 2 / 101,
'CRPIX1': 51,
'CRPIX2': 51,
'CRVAL1': 0,
'CRVAL2': 0,
}
In [5]:
def sample(nsize):
l, b = np.random.uniform(-1, 1, (2, nsize))
signal = func(l, b)
return l, b, signal
In [6]:
nsize = 51000
all_l, all_b, all_signal = sample(nsize)
In [7]:
l, b, signal = all_l[0:1], all_b[0:1], all_signal[0:1]
In [8]:
gridder = cygrid.WcsGrid(header)
kernelsize_fwhm = 3 / 60 # in degrees
kernelsize_sigma = kernelsize_fwhm / np.sqrt(8 * np.log(2))
kernel_support = 4 * kernelsize_sigma
gridder.set_kernel(
'gauss1d',
(kernelsize_sigma,),
kernel_support,
kernel_support / 2.
)
In [9]:
gridder.grid(l, b, signal)
In [10]:
slice_list = (
list(range(0, 20, 1)) +
list(range(20, 200, 10)) +
list(range(200, 2000, 100)) +
list(range(2000, 31001, 1000))
)
start_i = slice_list[:-1]
end_i = slice_list[1:]
In [11]:
my_wcs = WCS(header).celestial
fig = plt.figure(figsize=(12, 12))
axes = [
fig.add_subplot(2, 2, idx, projection=my_wcs)
for idx in range(1, 5)
]
imkw = dict(origin='lower', interpolation='nearest', cmap='viridis')
all_x, all_y = my_wcs.all_world2pix(all_l, all_b, 0)
points = axes[0].scatter(
all_x, all_y, c=all_signal, alpha=0,
cmap='viridis', vmin=-0.5, vmax=0.5
)
title_p = axes[0].set_title('Raw samples (# samples: {:5d})'.format(0), loc='center', fontsize=12)
im1 = axes[1].imshow(gridder.get_datacube(), **imkw, vmin=-0.5, vmax=0.5)
title1 = axes[1].set_title('Gridded data, A[i,j] / B[i,j] (# samples: {:5d})'.format(0), loc='center', fontsize=12)
im2 = axes[2].imshow(gridder.get_unweighted_datacube(), **imkw, vmin=-10, vmax=10)
title2 = axes[2].set_title('Gridded data * weight map, A[i,j] (# samples: {:5d})'.format(0), loc='center', fontsize=12)
im3 = axes[3].imshow(gridder.get_weights(), **imkw, vmin=0, vmax=10)
title3 = axes[3].set_title('Weight map, B[i,j] (# samples: {:5d})'.format(0), loc='center', fontsize=12)
for ax in axes:
ax.set_xlim((0, 100))
ax.set_ylim((0, 100))
lon, lat = ax.coords
lon.set_axislabel('R.A. [deg]')
lat.set_axislabel('Dec [deg]')
def init():
face_colors = points.get_facecolors()
face_colors[0, 3] = 1
im1.set_array(gridder.get_datacube())
im2.set_array(gridder.get_unweighted_datacube())
im3.set_array(gridder.get_weights())
return (
points, im1, im2, im3,
title_p, title1, title2, title3
)
def animate(i):
_sl = slice(start_i[i], end_i[i])
l, b, signal = all_l[_sl], all_b[_sl], all_signal[_sl]
gridder.grid(l, b, signal)
face_colors = points.get_facecolors()
face_colors[_sl, 3] = 1
im1.set_array(gridder.get_datacube())
im2.set_array(gridder.get_unweighted_datacube())
im3.set_array(gridder.get_weights())
title_p.set_text('Raw samples (# samples: {:5d})'.format(start_i[i]))
title1.set_text('Gridded data, A[i,j] / B[i,j] (# samples: {:5d})'.format(start_i[i]))
title2.set_text('Gridded data * weight map, A[i,j] (# samples: {:5d})'.format(start_i[i]))
title3.set_text('Weight map, B[i,j] (# samples: {:5d})'.format(start_i[i]))
return (
points, im1, im2, im3,
title_p, title1, title2, title3
)
anim = animation.FuncAnimation(
fig, animate, init_func=init, frames=len(start_i),
interval=200, blit=True
)
# this takes a while!
plt.close(anim._fig)
anim
Out[11]:
In [ ]: