---
title: Coordinate Systems and Spatial Transforms
output: rmarkdown::html_vignette
vignette: |
  %\VignetteIndexEntry{Coordinate Systems and Spatial Transforms}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
params:
  family: red
  preset: homage
css: albers.css
resource_files:
- albers.css
- albers.js
includes:
  in_header: |-
    <script src="albers.js"></script>
    <script>document.addEventListener('DOMContentLoaded',function(){document.body.classList.add('palette-red');});</script>

---

```{r setup, include = FALSE}
if (requireNamespace("ggplot2", quietly = TRUE) && requireNamespace("albersdown", quietly = TRUE)) ggplot2::theme_set(albersdown::theme_albers(family = params$family, preset = params$preset))
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  message = FALSE,
  warning = FALSE
)
```

```{r load-pkg}
library(neuroim2)
```

Every neuro object in neuroim2 — a `NeuroVol`, a `NeuroVec`, an ROI — carries a
`NeuroSpace` that answers one question: *where in the scanner room does this
voxel sit?* Get that mapping right and anatomical coordinates, atlas overlays,
and multi-subject registration all fall into place. Get it wrong and your
activations end up in the wrong hemisphere.

This vignette builds the mental model from the ground up: what a `NeuroSpace`
is, how the affine transform works, and how to move fluently between voxel
indices, grid coordinates, and millimetre world coordinates.

---

## What is a NeuroSpace?

A `NeuroSpace` is the spatial reference frame attached to every neuro image. It
bundles:

- **Grid dimensions** — how many voxels along each axis (`dim`)
- **Voxel spacing** — physical size of each voxel in millimetres (`spacing`)
- **Origin** — world coordinates of voxel `(1, 1, 1)` in mm (`origin`)
- **Affine transform** — the 4 × 4 matrix that maps voxel indices to mm (`trans`)
- **Axis orientation** — which anatomical direction each image axis points (`axes`)

```{r create-space}
sp <- NeuroSpace(
  dim     = c(64L, 64L, 40L),
  spacing = c(2, 2, 2),
  origin  = c(-90, -126, -72)
)
sp
```

```{r inspect-space}
dim(sp)       # grid dimensions
spacing(sp)   # voxel sizes in mm
origin(sp)    # world coords of voxel (1,1,1)
ndim(sp)      # number of spatial dimensions
```

The `NeuroSpace` is also the bridge to any volume or time-series built on top of
it. Call `space(x)` on any neuro object to retrieve it.

---

## The Affine Transform

The affine transform is a 4 × 4 homogeneous matrix stored in `trans(sp)`. It
maps a voxel position (expressed as a zero-based column vector) to millimetre
world coordinates:

```
[x_mm]   [M  t] [i]
[y_mm] = [      [j]
[z_mm]    0  1] [k]
[ 1  ]          [1]
```

where `M` is the 3 × 3 linear part (encodes spacing, rotation, and possible
shear) and `t` is the 3-element translation (the origin in mm).

For an axis-aligned image built from `spacing` and `origin`, `M` is diagonal:

```{r show-trans}
trans(sp)
```

The translation column (`trans(sp)[1:3, 4]`) recovers the origin, and the
diagonal of the linear block recovers the voxel sizes:

```{r decompose-affine}
# Translation column = origin
trans(sp)[1:3, 4]

# Diagonal of linear block = voxel sizes
diag(trans(sp)[1:3, 1:3])
```

The inverse affine (world -> voxel) is cached and accessible via
`inverse_trans(sp)`:

```{r inverse-trans}
inverse_trans(sp)
```

### Passing an explicit affine

You can supply a full 4 × 4 affine directly to `NeuroSpace()`. neuroim2 will
derive `spacing` and `origin` from it automatically:

```{r explicit-affine}
aff <- diag(c(3, 3, 4, 1))          # 3 × 3 × 4 mm voxels
aff[1:3, 4] <- c(-90, -126, -72)    # origin

sp_aff <- NeuroSpace(dim = c(60L, 60L, 35L), trans = aff)
spacing(sp_aff)   # derived from column norms of linear block
origin(sp_aff)    # derived from translation column
```

When the affine is oblique (rotated relative to the scanner axes), `spacing()`
returns the column norms of the linear block — the true physical voxel sizes —
while the diagonal of the linear block would be smaller. See the *Oblique
affines* section below.

---

## Voxel vs World Coordinates

neuroim2 uses **1-based grid indices** everywhere, matching R's array
conventions. The affine convention is therefore:

> `grid_to_coord` subtracts 1 from each 1-based index before applying the
> affine, placing voxel `(1, 1, 1)` at the origin.

There are two flavours of voxel addressing:

| Address type | Range | Description |
|:-------------|:------|:------------|
| Linear index | `1 ... prod(dim)` | Single integer; column-major (as in R arrays) |
| Grid index | `(1...d1, 1...d2, 1...d3)` | 3-tuple of 1-based integers |

And one world-coordinate system: millimetres, defined by the affine.

```{r coord-systems-diagram, echo=FALSE, fig.cap="The three addressing schemes and the functions that convert between them.", fig.width=5.5, fig.height=2.2}
op <- par(mar = c(0, 0, 0, 0))
plot.new()
plot.window(xlim = c(0, 10), ylim = c(0, 3))

# Boxes
rect(0.2, 0.8, 2.8, 2.2, col = "#dce8f5", border = "#3a7abf", lwd = 1.5)
rect(3.7, 0.8, 6.3, 2.2, col = "#dce8f5", border = "#3a7abf", lwd = 1.5)
rect(7.2, 0.8, 9.8, 2.2, col = "#dce8f5", border = "#3a7abf", lwd = 1.5)

text(1.5, 1.7, "Linear\nindex", cex = 0.85, font = 2)
text(1.5, 1.15, "1 ... prod(dim)", cex = 0.68, col = "#555555")
text(5.0, 1.7, "Grid\nindex", cex = 0.85, font = 2)
text(5.0, 1.15, "(i, j, k)  1-based", cex = 0.68, col = "#555555")
text(8.5, 1.7, "World\ncoords", cex = 0.85, font = 2)
text(8.5, 1.15, "x, y, z  mm", cex = 0.68, col = "#555555")

# Arrows and labels
arrows(2.85, 1.5, 3.65, 1.5, length = 0.08, lwd = 1.4, col = "#3a7abf")
arrows(3.65, 1.3, 2.85, 1.3, length = 0.08, lwd = 1.4, col = "#888888")
text(3.25, 1.75, "index_to_grid", cex = 0.58, col = "#3a7abf")
text(3.25, 1.05, "grid_to_index", cex = 0.58, col = "#888888")

arrows(6.35, 1.5, 7.15, 1.5, length = 0.08, lwd = 1.4, col = "#3a7abf")
arrows(7.15, 1.3, 6.35, 1.3, length = 0.08, lwd = 1.4, col = "#888888")
text(6.75, 1.75, "grid_to_coord", cex = 0.58, col = "#3a7abf")
text(6.75, 1.05, "coord_to_grid", cex = 0.58, col = "#888888")

# Shortcut arrows below
arrows(2.85, 0.9, 7.15, 0.9, length = 0.08, lwd = 1.2, col = "#3a7abf", lty = 2)
arrows(7.15, 0.7, 2.85, 0.7, length = 0.08, lwd = 1.2, col = "#888888", lty = 2)
text(5.0, 1.0, "index_to_coord", cex = 0.55, col = "#3a7abf")
text(5.0, 0.6, "coord_to_index", cex = 0.55, col = "#888888")

par(op)
```

---

## Coordinate Conversions

All conversion functions accept a `NeuroSpace` (or any neuro object, via its
attached space) as the first argument.

### Grid <-> linear index

```{r grid-index}
# Which linear index does grid position (10, 12, 5) map to?
grid_to_index(sp, matrix(c(10, 12, 5), nrow = 1))

# And back to grid
index_to_grid(sp, 8394L)
```

### Grid -> world coordinates

`grid_to_coord` subtracts 1 from the 1-based grid before applying the affine,
so grid `(1, 1, 1)` maps exactly to the origin:

```{r grid-to-coord}
grid_to_coord(sp, matrix(c(1, 1, 1), nrow = 1))   # should equal origin(sp)
grid_to_coord(sp, matrix(c(10, 12, 5), nrow = 1))
```

Multiple points are passed as a matrix with one row per point:

```{r grid-to-coord-multi}
pts <- matrix(c(
   1,  1,  1,
  32, 32, 20,
  64, 64, 40
), ncol = 3, byrow = TRUE)

grid_to_coord(sp, pts)
```

### World -> grid coordinates

```{r coord-to-grid}
coord_to_grid(sp, c(0, 0, 0))          # near isocenter
coord_to_grid(sp, matrix(c(0, 0, 0,
                           10, -20, 8), ncol = 3, byrow = TRUE))
```

### Convenience: linear index <-> world

```{r index-coord}
# Linear index 1 should land at the origin (pass integer)
index_to_coord(sp, 1L)

# World coord back to linear index
coord_to_index(sp, matrix(c(-90, -126, -72), nrow = 1))
```

### A concrete round-trip

```{r roundtrip}
idx   <- 12345L
grid  <- index_to_grid(sp, idx)
world <- grid_to_coord(sp, grid)
back  <- coord_to_index(sp, world)

cat("index:", idx, "-> grid:", grid, "-> world:", round(world, 2),
    "-> index:", back, "\n")
```

---

## Orientation Codes

Neuroimaging images can be stored in many orientations. The **orientation code**
(also called axis codes or RAS codes) tells you which anatomical direction each
image axis points:

| Letter | Anatomical direction |
|:------:|:---------------------|
| R | increasing -> Right |
| L | increasing -> Left |
| A | increasing -> Anterior |
| P | increasing -> Posterior |
| S | increasing -> Superior |
| I | increasing -> Inferior |

A code of `"RAS"` means: axis 1 runs left-to-right, axis 2 runs
posterior-to-anterior, axis 3 runs inferior-to-superior. This is the standard
NIfTI/MNI convention. `"LPI"` (sometimes called "neurological") reverses all
three.

`affine_to_axcodes()` reads the orientation directly from the affine matrix:

```{r axcodes}
affine_to_axcodes(trans(sp))
```

The default axis-aligned `NeuroSpace` built above uses the nearest-anatomy
convention inferred from the affine. You can also inspect the `axes` slot
directly:

```{r axes-slot}
axes(sp)
```

### Reorienting a space

`reorient()` flips and permutes the affine to match a target orientation string:

```{r reorient}
sp_ras <- reorient(sp, c("R", "A", "S"))
affine_to_axcodes(trans(sp_ras))
```

---

## Creating NeuroSpaces

### From dim + spacing + origin (axis-aligned)

The simplest case: an isotropic or anisotropic grid with no rotation.

```{r create-simple}
# 2 mm isotropic, MNI-ish origin
sp_mni <- NeuroSpace(
  dim     = c(91L, 109L, 91L),
  spacing = c(2, 2, 2),
  origin  = c(-90, -126, -72)
)
sp_mni
```

### From an explicit affine

When you have a NIfTI sform/qform matrix, pass it directly:

```{r create-from-affine}
# Oblique affine: slight off-diagonal terms (scanner tilt)
aff_obl <- matrix(c(
   2.0,  0.2,  0.0,  -90,
   0.0,  2.0,  0.1, -126,
   0.0,  0.0,  2.0,  -72,
   0.0,  0.0,  0.0,    1
), nrow = 4, byrow = TRUE)

sp_obl <- NeuroSpace(dim = c(91L, 109L, 91L), trans = aff_obl)

# spacing() returns column norms — the true physical voxel sizes
spacing(sp_obl)

# diagonal is not exactly 2 mm when the image is tilted
diag(aff_obl[1:3, 1:3])
```

### When does spacing() differ from the affine diagonal?

`spacing()` always returns the column norms of the 3 × 3 linear block — the
physical length of each voxel edge. For a pure diagonal affine these equal the
diagonal entries. For a rotated or oblique affine they differ. Always use
`spacing()` for physical voxel sizes; never rely on the diagonal directly.

You can also compute voxel sizes from any affine matrix with `voxel_sizes()`:

```{r voxel-sizes}
voxel_sizes(aff_obl)
```

---

## Dimension Manipulation

### Adding a time dimension: add_dim

`add_dim()` extends a 3D `NeuroSpace` to 4D by appending a new dimension. The
spatial affine is preserved unchanged:

```{r add-dim}
sp_3d <- NeuroSpace(c(64L, 64L, 40L), spacing = c(2, 2, 2),
                    origin = c(-90, -126, -72))
sp_4d <- add_dim(sp_3d, 200)   # 200 time points
dim(sp_4d)
trans(sp_4d)                   # 4x4 spatial affine unchanged
```

### Dropping the time dimension: drop_dim

`drop_dim()` removes the last (or a named) dimension:

```{r drop-dim}
sp_back <- drop_dim(sp_4d)
dim(sp_back)
all.equal(trans(sp_back), trans(sp_3d))   # affine preserved
```

This round-trip is used internally whenever a 4D volume is sliced to a single
time point.

---

## Resampling and Reorientation

### resample — change voxel size, keep geometry

`resample()` resamples a `NeuroVol` into a target space (or another volume's
space). To demonstrate, build a small volume and resample it to a coarser grid:

```{r resample-space}
vol_mni <- NeuroVol(array(rnorm(prod(dim(sp_mni))), dim(sp_mni)), sp_mni)
sp_coarse <- NeuroSpace(c(45L, 54L, 45L), spacing = c(4, 4, 4),
                        origin = origin(sp_mni))
vol_coarse <- resample(vol_mni, sp_coarse)
dim(vol_coarse)
spacing(space(vol_coarse))
```

### deoblique — remove scanner tilt

Many scanners acquire data with a slight tilt relative to the MNI axes. The
resulting affine has non-zero off-diagonal elements (an *oblique* affine).
`deoblique()` builds an axis-aligned output space that encloses the original
field of view, using the minimum input voxel size by default (AFNI-style):

```{r deoblique}
sp_deob <- deoblique(sp_obl)
affine_to_axcodes(trans(sp_deob))   # now axis-aligned
obliquity(trans(sp_deob))           # near zero
```

When passed a `NeuroVol`, `deoblique()` also resamples the image data into the
new space.

---

## Common Gotchas

**1. Oblique affines and spacing()**
If `diag(trans(sp)[1:3, 1:3])` does not match `spacing(sp)`, the affine is
oblique. Use `obliquity(trans(sp))` to quantify the tilt angle (in radians) per
axis.

```{r gotcha-obliquity}
obliquity(aff_obl)        # non-zero: image is slightly tilted
obliquity(trans(sp_mni))  # near zero: axis-aligned
```

**2. 1-based grid indices**
R uses 1-based indexing. `grid_to_coord()` handles this by subtracting 1 before
applying the affine, so voxel `(1, 1, 1)` lands at `origin(sp)`. If you ever
use raw affine arithmetic, remember to subtract 1 from your 1-based grid
coordinates first.

**3. NIfTI sform vs qform**
NIfTI files store two possible affines (sform and qform). `read_vol()` and
`read_vec()` follow the NIfTI priority rules: sform_code > 0 -> use sform;
otherwise use qform. The resulting `trans(space(img))` reflects whichever was
chosen. You can inspect the header with `read_header()` if you need to see both.

**4. Float32 precision**
NIfTI stores affine coefficients as 32-bit floats. neuroim2 rounds the affine
to 7 significant figures on construction (`signif(trans, 7)`) to match this
precision. Round-trip coordinates through `index_to_coord` / `coord_to_index`
may therefore show sub-voxel floating-point noise at the ~0.001 mm level.

---

## Quick Reference

| Function | Input -> Output | Typical use |
|:---------|:---------------|:------------|
| `grid_to_index(sp, mat)` | grid -> linear index | looking up voxel data |
| `index_to_grid(sp, idx)` | linear index -> grid | converting R array subscripts |
| `grid_to_coord(sp, mat)` | grid -> mm | overlay on anatomy |
| `coord_to_grid(sp, mat)` | mm -> grid | atlas lookup |
| `index_to_coord(sp, idx)` | linear index -> mm | shortcut past grid |
| `coord_to_index(sp, mat)` | mm -> linear index | mask extraction |
| `affine_to_axcodes(aff)` | affine -> `"RAS"` etc. | orientation check |
| `voxel_sizes(aff)` | affine -> mm vector | physical voxel size |
| `obliquity(aff)` | affine -> radians | tilt check |
| `add_dim(sp, n)` | 3D space -> 4D space | attach time axis |
| `drop_dim(sp)` | 4D space -> 3D space | strip time axis |
| `resample(sp, spacing)` | space -> resampled space | change resolution |
| `reorient(sp, codes)` | space -> reoriented space | standardise orientation |
| `deoblique(sp)` | oblique space -> aligned space | remove scanner tilt |

## Where to go next

- `vignette("ImageVolumes")` — creating and manipulating `NeuroVol` objects
- `vignette("NeuroVector")` — working with 4D time-series (`NeuroVec`)
- `vignette("Resampling")` — image resampling in depth
- `vignette("regionOfInterest")` — ROI construction and coordinate-based queries
- `?NeuroSpace`, `?affine_to_axcodes`, `?deoblique` — function reference pages