convolve_weights#

regridding.convolve_weights(weights, kernel, axis_input=None, axis_output=None)[source]#

Convolve the output of a set of weights with a kernel, such as a point-spread function.

If the weights computed by regridding.weights() are the sparse matrix \(W\), mapping each input cell \(j\) onto the output cells \(i\), this computes the sparse matrix \(P W\), where

\[P_{i' i} = K_i(i' - i)\]

spreads whatever lands in output cell \(i\) over the cells \(i'\) around it. Applying the result with regridding.regrid_from_weights() is the same as resampling and then convolving with \(K\), but as a single set of weights it can be reused, and transposed by regridding.transpose_weights() or regridding.transpose_weights_conservative() to give the transpose of the whole operation.

Parameters:

  • weights (tuple[ndarray, tuple[int, ...], tuple[int, ...]]) – Weights computed by regridding.weights().

  • kernel (ndarray) –

    The kernel, \(K\). Its last len(axis_output) axes are the kernel itself, one for each of axis_output, in the same order. Along an axis of length \(n\), the element at index \(\lfloor n / 2 \rfloor\) is the center, which is where scipy.ndimage.convolve() places it.

    Any axes before those are broadcast against the output grid, shape_output, so the kernel may be different for each element of the orthogonal axes, such as for each wavelength, or it may vary across the output grid along the resampled axes. A kernel which varies across the output grid is indexed by the cell the light lands in before it is spread, \(i\) above.

    The kernel must be dimensionless.

  • axis_input (None | int | Sequence[int]) – The resampled axes of the input grid, as given to regridding.weights() and regridding.regrid_from_weights(). If None, all the axes of the input grid are resampled. This is only needed when the kernel varies along an orthogonal axis which the weights are broadcast along, since the shape of the input grid then has to be broadcast too.

  • axis_output (None | int | Sequence[int]) – The resampled axes of the output grid, as given to regridding.weights() and regridding.regrid_from_weights(). If None, all the axes of the output grid are resampled.

Returns:

  • The convolved weights, with the same shapes as weights, unless the

  • kernel varies along an orthogonal axis the weights are broadcast along,

  • in which case the weights and both shapes are broadcast along it.

Raises:

ValueError – If kernel has too few axes, if its leading axes cannot be broadcast to the output grid, or if it is not dimensionless; or if axis_output does not describe the grid the weights were built for, which shows as the weights having more orthogonal axes than it leaves or as an output index outside the grid it selects.

Return type:

tuple[ndarray, tuple[int, …], tuple[int, …]]

Notes

The kernel acts on the output cells, so it describes how the contents of a cell are redistributed among its neighbors, rather than how a point is blurred. For a point-spread function \(h\) sampled on cells of width \(\Delta\), that is \(h\) convolved with a cell twice, once for the cell the light lands in and once for the cell it is collected in, and then sampled at the cell centers.

Light which the kernel spreads beyond the edge of the output grid is lost, as it would be off the edge of a sensor, so the total of the result is reduced near the edges.

Each input cell is convolved independently of every other, in a scratch box covering its footprint: the bounding box of the output cells it reaches, grown by the kernel. Each of its weights is spread through the kernel into the box, and the box is then read out in order, so the cost is one multiply-add for each weight and element of the kernel, plus one visit to each cell of the box. Weights ordered by input cell, as regridding.weights() returns them, come out ordered and with each (input, output) pair once.

A kernel which varies along some of the resampled axes is stored with one row for each cell along those axes only, so a kernel which varies along one axis of a large grid costs no more than that axis.

The result has more weights than weights does, by roughly the factor by which the kernel grows the footprint of an input cell: for a cell which lands on two output cells across, a kernel three cells across makes the weights about four times as long, and seven cells across about eighteen times.

Weights built with device="cuda" are convolved on the device, and the result is left there. Their empty slots, which carry an index of -1, are dropped.

See also

regridding.weights(), regridding.regrid_from_weights(), regridding.transpose_weights_conservative()

Examples

Rotate an array onto a new grid, and blur it with a Gaussian kernel in the same set of weights.

import numpy as np
import matplotlib.pyplot as plt
import regridding

# Define the input grid
x_input = np.linspace(-4, 4, num=17)
y_input = np.linspace(-4, 4, num=17)
x_input, y_input = np.meshgrid(x_input, y_input, indexing="ij")

# Rotate the input grid
angle = 0.3
x_rotated = x_input * np.cos(angle) - y_input * np.sin(angle)
y_rotated = x_input * np.sin(angle) + y_input * np.cos(angle)

# Define a uniform output grid
x_output = np.linspace(-6, 6, num=25)
y_output = np.linspace(-6, 6, num=25)
x_output, y_output = np.meshgrid(x_output, y_output, indexing="ij")

# Define an array of values with two bright cells
values_input = np.zeros((16, 16))
values_input[4, 4] = 1
values_input[10, 8] = 1

# Save the weights which rotate the input array onto the output grid
weights = regridding.weights(
    coordinates_input=(x_rotated, y_rotated),
    coordinates_output=(x_output, y_output),
    method="conservative",
)

# Define a Gaussian kernel, five cells across
offset = np.arange(-2, 3)
kernel = np.exp(-np.square(offset) / 2)
kernel = kernel[:, np.newaxis] * kernel[np.newaxis, :]
kernel = kernel / kernel.sum()

# Blur the output of the weights with the kernel
weights_blurred = regridding.convolve_weights(weights, kernel)

# Apply both sets of weights
values_rotated = regridding.regrid_from_weights(
    *weights,
    values_input=values_input,
)
values_blurred = regridding.regrid_from_weights(
    *weights_blurred,
    values_input=values_input,
)

# Plot the original, rotated, and blurred arrays
fig, axs = plt.subplots(
    ncols=3,
    sharex=True,
    sharey=True,
    figsize=(9, 3.4),
    constrained_layout=True,
)
axs[0].pcolormesh(x_rotated, y_rotated, values_input);
axs[0].set_title(f"original, total {values_input.sum():.2f}");
axs[1].pcolormesh(x_output, y_output, values_rotated);
axs[1].set_title(f"rotated, total {values_rotated.sum():.2f}");
axs[2].pcolormesh(x_output, y_output, values_blurred);
axs[2].set_title(f"rotated and blurred, total {values_blurred.sum():.2f}");
for ax in axs:
    ax.set_aspect("equal");
../_images/regridding.convolve_weights_0_1.png