neuroim2

CRAN

Core neuroimaging data structures for fMRI analysis in R

Installinstall.packages('neuroim2')

Sagittal, coronal and axial views of the MNI152 template with a thresholded simulated statistic map: positive values in red-yellow, negative in blue

A simulated statistic map on the package’s MNI152 template, drawn with plot_ortho(style = "dark") in make-figures.R.

Classes NeuroVol · NeuroVec · NeuroSpace · ROI

Input / output NIfTI read and write, memory-mapped 4D data

Spatial operations Smoothing, resampling, connected components

Analysis building blocks Spherical ROIs, searchlight iterators, ROI time series

What it does

Images carry their coordinate space

A NeuroVol stores the voxel array together with a NeuroSpace: grid dimensions, voxel spacing, origin and orientation. grid_to_coord() and coord_to_grid() convert between voxel indices and millimetre coordinates; the cluster table further down uses grid_to_coord() to report its peaks in millimetres.

Twelve axial slices of the MNI152 template arranged in two rows

library(neuroim2)

# A downsampled MNI152 template (about 4 mm voxels)
# ships with the package.
f <- system.file("extdata", "mni_downsampled.nii.gz",
                 package = "neuroim2")
mni <- read_vol(f)
space(mni)  # grid, spacing, origin, orientation

fig <- plot_montage(mni, n_slices = 12, ncol = 6,
                    label_size = 4.4)

Smoothing in millimetres

gaussian_blur() takes the kernel width as a FWHM in the image’s spatial units, here millimetres, so a given fwhm gives the same physical smoothing at any voxel size. plot_overlay() draws the thresholded map on the anatomy with a colour scale centred on zero.

Six axial slices of the template with a thresholded simulated z map: positive values in red-yellow, negative in blue

Simulated data: smoothed Gaussian noise, thresholded at |z| > 2.

# Simulated statistic map: Gaussian noise inside the
# brain, smoothed to 16 mm FWHM
set.seed(2)
mask  <- mni > quantile(mni[mni > 0], 0.25)
noise <- NeuroVol(rnorm(prod(dim(mni))), space(mni))
stat  <- gaussian_blur(noise, mask, fwhm = 16)
stat  <- stat / sd(stat[mask > 0])  # unit variance

fig <- plot_overlay(mni, stat, ov_thresh = 2,
                    ov_symmetric = TRUE, n_slices = 6,
                    ncol = 3, cbar_title = "z",
                    label_size = 4.4)

Clusters and their table

conn_comp() labels the connected components of voxels above threshold and, by default, returns a cluster table with the size, volume and peak voxel of each component. The table below is built from it.

Orthogonal views through the peak of the largest cluster; cluster 1 in coral and cluster 2 in teal, matching the table

Cluster	Voxels	Volume (mm³)	Peak (x, y, z mm)	Peak z-score
1	306	19,875	(-16, -56, -2)	3.45
2	126	8,184	(21, 13, -2)	3.38
3	63	4,092	(-48, 37, 27)	3.07
4	62	4,027	(21, -76, 59)	3.04
5	32	2,078	(21, 5, 67)	2.89

The five largest clusters where z > 2 in the simulated map; colours match the figure. Clusters 3–5 lie outside these three slices. Volumes use the template’s voxels, about 4 mm on a side.

# Connected components where z > 2, with a cluster table
cc   <- conn_comp(stat, threshold = 2, cluster_table = TRUE)
top5 <- head(cc$cluster_table, 5)

# Label the five largest clusters 1-5 (NA is transparent)
labels <- NeuroVol(array(NA_real_, dim(mni)), space(mni))
idx <- as.array(cc$index)
for (k in 1:5) labels[idx == top5$index[k]] <- k

# One colour per label: a stepped palette over 0.5-5.5
cols <- c("#c4452b", "#17677a", "#d99a00",
          "#7b3fa0", "#2e8b57")
pal  <- rep(cols, each = 52)[1:256]

peak <- unlist(top5[1, c("x", "y", "z")])  # voxel indices
fig  <- plot_ortho(mni, coord = peak, overlay = labels,
                   ov_cmap = pal, ov_range = c(0.5, 5.5),
                   colorbar = FALSE, label_size = 4.4)

ROIs and searchlights

spherical_roi() (radius in mm) and cuboid_roi() (half-width in voxels) return an ROIVol holding the coordinates and value of each voxel. The figure shows a sphere of 16 mm radius.

searchlight() returns a lazy iterator with one sphere of 8 mm radius for each voxel in the mask (25,277 here), each built only when it is requested. rMVPA’s volumetric searchlight analyses use these iterators.

Orthogonal views of the template with a spherical ROI of 16 mm radius highlighted in teal

# A spherical ROI and a searchlight iterator
centre <- c(24, 30, 30)                  # voxel indices
roi <- spherical_roi(mask, centroid = centre,
                     radius = 16, fill = 1)  # radius, mm
sl  <- searchlight(mask, radius = 8)   # lazy iterator

sphere <- NeuroVol(array(0, dim(mni)), space(mni))
sphere[coords(roi)] <- 1
fig <- plot_ortho(mni, coord = centre, overlay = sphere,
                  ov_thresh = 0.5, colorbar = FALSE,
                  ov_cmap = c("#17677a", "#17677a"),
                  label_size = 4.4)

4D runs read from disk on demand

read_vec(mode = "mmap") memory-maps an uncompressed .nii file instead of loading it, and series() returns the time series for any set of voxels.

The run has 120 volumes on the template grid. Loaded into R as doubles it takes 120.2 MB; written as 4-byte floats it is a 60 MB file. The memory-mapped object takes 0.01 MB, because values are read from the file only when they are accessed. The mean time course of the 251-voxel ROI follows the simulated design (r = 0.98).

Mean ROI signal over 240 seconds rising during each shaded 20-second block and returning to baseline between blocks

Simulated 4D run: 120 volumes (TR 2 s) with a block design added inside the ROI.

# Simulated 4D run: 120 volumes (TR = 2 s), with a
# 20 s on / 20 s off block design inside the ROI
n_t    <- 120
design <- rep(rep(c(0, 1), each = 10), length.out = n_t)
run    <- array(rnorm(prod(dim(mni)) * n_t),
                c(dim(mni), n_t))
in_roi <- rep(as.array(sphere) > 0, n_t)
run[in_roi] <- run[in_roi] +
  0.8 * rep(design, each = sum(as.array(sphere) > 0))
vec <- NeuroVec(run, add_dim(space(mni), n_t))

# Write to NIfTI, then memory-map it back:
# voxels are read from disk on demand
f <- tempfile(fileext = ".nii")
write_vec(vec, f, data_type = "FLOAT")
mapped <- read_vec(f, mode = "mmap")

roi_ts  <- series(mapped, coords(roi))  # time x voxel
mean_ts <- rowMeans(roi_ts)

library(ggplot2)
df  <- data.frame(t = seq_len(n_t) * 2, y = mean_ts,
                  on = design)
fig <- ggplot(df) +
  geom_rect(aes(xmin = t - 2, xmax = t, ymin = -Inf,
                ymax = Inf, alpha = on * 0.12),
            fill = "#c4452b") +
  scale_alpha_identity() +
  geom_line(aes(t, y), linewidth = 0.7, colour = "#0f1d26") +
  labs(x = "Time (s)", y = "Mean ROI signal") +
  theme_minimal(base_size = 13)

Built on neuroim2

fmrireg fMRI regression models

rMVPA Multivariate pattern analysis

neurosurf Cortical surface data

neuroatlas Atlases and parcellations

Development version

The install command above gives the CRAN release. To install the development version from GitHub:

pak::pak("bbuchsbaum/neuroim2")

Every figure on this page is drawn by software/showcase/neuroim2/make-figures.R in the site’s source, and the code shown with each figure is the code that drew it.