---
title: 'neuroim2: a tour'
output:
  rmarkdown::html_vignette:
    css: albers.css
    includes:
      in_header:
      - albers-header.inc
      - albers-header.html
    toc: yes
    toc_depth: 2.0
resource_files:
- albers.css
- albers-fonts.css
- albers.js
- albers-header.inc
- albers-header.html
- fonts
params:
  family: red
  preset: interaction

vignette: |
  %\VignetteIndexEntry{neuroim2: a tour}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
if (requireNamespace("ragg", quietly = TRUE)) knitr::opts_chunk$set(dev = "ragg_png")
if (requireNamespace("systemfonts", quietly = TRUE)) albersdown::albers_register_fonts()
if (requireNamespace("ggplot2", quietly = TRUE) && requireNamespace("albersdown", quietly = TRUE)) ggplot2::theme_set(albersdown::theme_albers(family = params$family, preset = params$preset))
source("_common.R")
```

```{r albers-classes, echo=FALSE, results='asis'}
cat(sprintf(
  paste0(
    '<script>document.addEventListener("DOMContentLoaded",function(){',
    'document.body.classList.remove("palette-red","palette-lapis","palette-ochre","palette-teal","palette-green","palette-violet","preset-homage","preset-interaction","preset-study","preset-structural","preset-adobe","preset-midnight");',
    'document.body.classList.add("palette-%s","preset-%s");',
    '});</script>'
  ),
  params$family,
  params$preset
))
```

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

## Look at an image

```{r anatomy, fig.cap = "plot() on a NeuroVol: a montage of evenly spaced axial slices.", fig.alt = "Nine axial slices of an anatomical brain image.", fig.height = 6.5}
anat <- read_vol(system.file("extdata", "mni_downsampled.nii.gz", package = "neuroim2"))
plot(anat, zlevels = round(seq(6, 40, length.out = 9)))
```

That is a brain, read from a NIfTI file in one line and drawn in another. It is
also an ordinary R array:

```{r array-like}
dim(anat)
anat[24, 28, 24]
max(anat / 2)
```

Indexing, arithmetic and comparison all work. The part that matters is what
comes back: a comparison gives an image, not a bare logical array.

```{r still-an-image}
class(anat > 100)[1]
```

Everything in this package is built on that — objects that behave like arrays
but never forget their geometry.

## From millimetres to voxels

An image is not just numbers; it is numbers *somewhere*. That somewhere is a
`NeuroSpace`, attached to every object, and it is what separates neuroimaging
from array programming.

Here we switch images. Anatomical and functional scans are separate
acquisitions on separate grids, and the rest of this tour works in the
functional one — a brain mask on the grid an EPI run was acquired on.

```{r mask}
mask <- read_vol(system.file("extdata", "global_mask2.nii.gz", package = "neuroim2"))

spacing(mask)
affine_to_axcodes(trans(mask))
```

`spacing()` says these voxels are 3.5 x 3.5 x 3.7 mm. The axis codes say the
first array axis increases towards the left (`L`), the second towards the
anterior (`A`), and the third towards the superior (`S`). `space(mask)` prints all of it
at once, and `origin()` and `trans()` pull out the remaining pieces.

Those codes describe the voxel grid. The millimetre coordinates themselves are
always RAS — x increasing to the right, y forward, z up — so a negative x is on
the left whichever way the grid happens to run.

With that, the package can answer the question you actually care about: *which
voxel sits at this location?*

```{r coord}
g <- coord_to_grid(mask, c(-34, -28, 10))
g
```

A point 34 mm left, 28 mm behind and 10 mm above this image's world origin
becomes a position on its grid. The answer is fractional because a millimetre
coordinate rarely lands on a voxel centre. Going back is exact:

```{r coord-back}
grid_to_coord(mask, matrix(g, nrow = 1))
```

Keep the fractional position if you are going to interpolate, and round it only
when you need an integer index — rounding moves you up to half a voxel, which
here is about 1.8 mm on any one axis:

```{r coord-round}
grid_to_coord(mask, matrix(round(g), nrow = 1))
```

This image is in native scanner space, not a template space, so its origin is
wherever the scanner put it. Reading anatomical labels off these numbers would
require a registered image; `vignette("spaces-and-coordinates")` covers how to
tell what space you are actually in.

## Add a time axis

Functional data has a fourth dimension. `read_vec()` reads 4D files; here we
simulate a run instead, so the time series in this article have the temporal
structure real BOLD data has.

```{r bold}
bold <- simulate_fmri(mask, n_time = 60, seed = 1)

dim(bold)
```

A `NeuroVec` is 60 volumes sharing one spatial frame. Pull the time course of a
single voxel with `series()`:

```{r series, fig.cap = "One voxel's simulated BOLD time course.", fig.alt = "Line plot of a single voxel time series over 60 scans."}
vox <- round(g)
ts <- series(bold, vox[1], vox[2], vox[3])

plot(ts, type = "l", xlab = "scan", ylab = "signal", main = "Single voxel")
```

## Summarise a region

Single voxels are noisy, so analyses usually work over a region. Put a 10 mm
sphere at that voxel and extract every time course inside it.

```{r roi}
roi <- spherical_roi(mask, vox, radius = 10, nonzero = TRUE)
length(roi)
```

Note the mixed units: the centre is a voxel index, the radius is millimetres.
`nonzero = TRUE` drops voxels outside the mask — this sphere sits well inside
the brain so it removes none, but near the edge it does the work.

```{r roi-values}
roi_ts <- series_roi(bold, roi)
dim(values(roi_ts))
```

`values()` returns **scans by voxels**, so averaging across columns gives the
region's mean time course:

```{r roi-mean, fig.cap = "Single voxel and region mean on a shared axis. Averaging halves the amplitude.", fig.alt = "Two time series on the same axes; the region mean has visibly smaller swings.", fig.height = 3.6}
roi_mean <- rowMeans(values(roi_ts))

plot(ts, type = "l", col = "grey60", xlab = "scan", ylab = "signal",
     main = paste("1 voxel vs mean of", length(roi)))
lines(roi_mean, lwd = 2)
legend("topright", c("voxel", "ROI mean"), col = c("grey60", "black"),
       lwd = c(1, 2), bty = "n")
```

```{r roi-sd}
c(voxel = sd(ts), roi_mean = sd(roi_mean))
```

Averaging cuts the amplitude roughly in half rather than by the square root of
85, because the simulation gives nearby voxels shared spatial structure — the
same reason real BOLD voxels are not independent samples.

That is the shape of most work in this package: **define a spatial support,
extract values from it, reduce them to something smaller.**

## Make a result and write it out

Reductions run the other way too — from a 4D series down to one volume. Here is
each voxel's temporal standard deviation, a standard quality-control map:

```{r sdmap, fig.cap = "Temporal standard deviation per voxel, on a robust intensity range.", fig.alt = "Nine axial slices of a temporal standard deviation map.", fig.height = 6.5}
mat <- as.matrix(bold)
sd_map <- NeuroVol(apply(mat, 1, sd), drop_dim(space(bold)))

plot(sd_map,
  zlevels = round(seq(4, 22, length.out = 9)),
  irange = c(0, quantile(sd_map[sd_map > 0], 0.99))
)
```

Three things are worth pulling out of those four lines. `as.matrix()` flattens a
`NeuroVec` to **voxels by time** — the transpose of what `values()` gave us
above, so check the orientation whenever you switch between them. `drop_dim()`
takes the 4D space down to the matching 3D one, which is what keeps the result
aligned with its source. And QC maps need a robust intensity range: a handful of
edge voxels here run to five times the 99th percentile, and on the default full
range they would flatten everything else to black.

The result is a proper `NeuroVol`, so writing it produces a NIfTI any other tool
can read — including this one:

```{r write}
out <- tempfile(fileext = ".nii.gz")
write_vol(sd_map, out)

back <- read_vol(out)
all.equal(spacing(back), spacing(sd_map))
max(abs(back - sd_map))
```

The values differ in the seventh decimal place because NIfTI stored them as
32-bit floats. Compare images with a tolerance, not `identical()`.

```{r cleanup, include = FALSE}
unlink(out)
```

## Where to go next

You have now read images, inspected their geometry, moved between millimetres
and voxels, extracted a region's time course, reduced a series to a map, and
written it back to disk. Three articles finish the foundations, in this order:

- `vignette("spaces-and-coordinates")` — the affine, orientation codes and
  conversions in full. Read this next; everything else assumes it.
- `vignette("volumes-and-vectors")` — building, slicing and combining the
  containers.
- `vignette("reading-and-writing")` — file formats, headers, and the header
  problem that silently produces wrong coordinates.

Then pick a task: `vignette("regions-and-searchlights")`,
`vignette("resampling-and-orientation")`, `vignette("smoothing-and-filtering")`,
`vignette("visualization")`, or `vignette("large-data")`.