\(\newcommand{L}[1]{\| #1 \|}\newcommand{VL}[1]{\L{ \vec{#1} }}\newcommand{R}[1]{\operatorname{Re}\,(#1)}\newcommand{I}[1]{\operatorname{Im}\, (#1)}\)

Making coordinate arrays with meshgrid

affine_transform works by using voxel coordinate implied by the output_shape, and transforming those. See: Resampling with images of different shapes.

numpy.meshgrid is a way of making an actual coordinate grid.

This is particularly useful when we want to use the more general form of image resampling in scipy.ndimage.map_coordinates.

If we have some shape – say output_shape – then this implies a set of coordinates. Let’s say output_shape = (5, 4) – implying a 2D array.

The implied coordinate grid will therefore have one coordinate for each pixel (2D voxel) in the (5, 4) array.

Because this array is 2D, there are two coordinate values for each pixel. For example, the coordinate of the first element in the array is (0, 0). We can make these i- and j- coordinates with meshgrid:

>>> import numpy as np
>>> i_coords, j_coords = np.meshgrid(range(5), range(4), indexing='ij')
>>> i_coords
array([[0, 0, 0, 0],
       [1, 1, 1, 1],
       [2, 2, 2, 2],
       [3, 3, 3, 3],
       [4, 4, 4, 4]])
>>> j_coords
array([[0, 1, 2, 3],
       [0, 1, 2, 3],
       [0, 1, 2, 3],
       [0, 1, 2, 3],
       [0, 1, 2, 3]])

We can make this into a shape (2, 5, 4) array where the first axis contains the (i, j) coordinate.

>>> coordinate_grid = np.array([i_coords, j_coords])
>>> coordinate_grid.shape
(2, 5, 4)

Because we have not done any transformation on the coordinate, the i, j coordinate will be the same as the index we use to get the i, j coordinate:

>>> coordinate_grid[:, 0, 0]
array([0, 0])
>>> coordinate_grid[:, 1, 0]
array([1, 0])
>>> coordinate_grid[:, 0, 1]
array([0, 1])

This is the coordinate grid implied by a shape of (5, 4).

Now imagine I wanted to do a transformation on these coordinates. Say I wanted to add 2 to the first (i) coordinate:

>>> coordinate_grid[0, :, :] += 2

Now my coordinate grid expresses a mapping between a given (\(i, j\)) coordinate, and the new coordinate (\(i', j'\). I look up the new coordinate using the \(i, j\) index into the coordinate grid:

>>> coordinate_grid[:, 0, 0]  # look up new coordinate for (0, 0)
array([2, 0])
>>> coordinate_grid[:, 1, 0]  # look up new coordinate for (1, 0)
array([3, 0])
>>> coordinate_grid[:, 0, 1]  # look up new coordinate for (0, 1)
array([2, 1])

This means we can use these coordinate grids as a mapping from an input set of coordinates to an output set of coordinates, for each pixel / voxel.

As you can imagine, meshgrid extends to three dimensions or more:

>>> output_shape = (5, 6, 7)
>>> I, J, K = output_shape
>>> i_coords, j_coords, k_coords = np.meshgrid(range(I),
...                                            range(J),
...                                            range(K),
...                                            indexing='ij')
>>> coordinate_grid = np.array([i_coords, j_coords, k_coords])
>>> coordinate_grid.shape
(3, 5, 6, 7)

>>> coordinate_grid[:, 0, 0, 0]
array([0, 0, 0])
>>> coordinate_grid[:, 1, 0, 0]
array([1, 0, 0])
>>> coordinate_grid[:, 0, 1, 0]
array([0, 1, 0])
>>> coordinate_grid[:, 0, 0, 1]
array([0, 0, 1])