numpy_meshgrid.Rmd

---
jupyter:
  jupytext:
    text_representation:
      extension: .Rmd
      format_name: rmarkdown
      format_version: '1.2'
      jupytext_version: 1.11.5
---

$\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`:

```{python}
import numpy as np
i_coords, j_coords = np.meshgrid(range(5), range(4), indexing='ij')
i_coords
j_coords
```

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

```{python}
coordinate_grid = np.array([i_coords, j_coords])
coordinate_grid.shape
```

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:

```{python}
coordinate_grid[:, 0, 0]
coordinate_grid[:, 1, 0]
coordinate_grid[:, 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:

```{python}
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:

```{python}
coordinate_grid[:, 0, 0]  # look up new coordinate for (0, 0)
coordinate_grid[:, 1, 0]  # look up new coordinate for (1, 0)
coordinate_grid[:, 0, 1]  # look up new coordinate for (0, 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:

```{python}
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
```

```{python}
coordinate_grid[:, 0, 0, 0]
coordinate_grid[:, 1, 0, 0]
coordinate_grid[:, 0, 1, 0]
coordinate_grid[:, 0, 0, 1]
```

<!-- vim:ft=rst -->
<!-- Course -->
<!-- BIC -->
<!-- Python distributions -->
<!-- Version control -->
<!-- Editors -->
<!-- Python and common libraries -->
<!-- IPython -->
<!-- Virtualenv and helpers -->
<!-- Pypi and packaging -->
<!-- Mac development -->
<!-- Windows development -->
<!-- Nipy and friends -->
<!-- FMRI datasets -->
<!-- Languages -->
<!-- Imaging software -->
<!-- Installation -->
<!-- Tutorials -->
<!-- MB tutorials -->
<!-- Ideas -->
<!-- Psych-214 -->
<!-- People -->
<!-- Licenses -->
<!-- Neuroimaging stuff -->
<!-- OpenFMRI projects -->
<!-- Unix -->
<!-- Substitutions -->