Working with Volumes and Vectors
In neuroim2: Data Structures for Brain Imaging Data

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
)
suppressPackageStartupMessages(library(neuroim2))

neuroim2 has two basic data containers you reach for constantly: NeuroVol for one 3D image and NeuroVec for a 4D stack of volumes. This article shows how to move between them, inspect geometry, and extract the small pieces you actually analyze.

Start with a single volume

Use read_vol() when the file is a 3D image or when you want one volume at a time.

vol_file <- system.file("extdata", "global_mask2.nii.gz", package = "neuroim2")
vol <- read_vol(vol_file)

dim(vol)
spacing(vol)
origin(vol)

stopifnot(length(dim(vol)) == 3L)
stopifnot(all(spacing(vol) > 0))

The returned object behaves like an array, but it also carries a NeuroSpace so you can keep voxel values tied to real spatial coordinates.

space(vol)
vol[1, 1, 1]

Read a 4D time series

Use read_vec() when the file contains a time series or another stack of aligned volumes.

vec_file <- system.file("extdata", "global_mask_v4.nii", package = "neuroim2")
vec <- read_vec(vec_file)

dim(vec)
vec

stopifnot(length(dim(vec)) == 4L)
stopifnot(dim(vec)[4] > 1L)

A NeuroVec still has one shared 3D spatial frame. The fourth dimension is the axis you usually split, subset, and summarize.

Pull out the volumes you need

sub_vector() is the direct way to keep only selected timepoints.

first_two <- sub_vector(vec, 1:2)

dim(first_two)

stopifnot(dim(first_two)[4] == 2L)

Convert to a matrix when you need table-like work

as.matrix() is useful when you want one row per voxel and one column per timepoint.

mat <- as.matrix(vec)

dim(mat)
mat[1:4, 1:2]

stopifnot(nrow(mat) == prod(dim(vec)[1:3]))
stopifnot(ncol(mat) == dim(vec)[4])

That representation is convenient for temporal summaries, scaling, or passing the data into general R modeling code.

Extract one voxel or one ROI time series

For a single voxel, use series().

voxel_ts <- series(vec, 12, 12, 12)
voxel_ts

For a spatial region, create an ROI and use series_roi().

roi <- spherical_roi(drop(first_two[[1]]), c(12, 12, 12), radius = 6)
roi_ts <- series_roi(vec, roi)

dim(values(roi_ts))

stopifnot(length(roi) > 0L)
stopifnot(nrow(values(roi_ts)) == dim(vec)[4])

Switch to sparse storage when the support matters

If only a subset of voxels should be considered present, read through a mask and keep the result sparse.

mask_vol <- read_vol(vec_file) > 0
sparse_vec <- read_vec(vec_file, mask = mask_vol)

class(sparse_vec)
dim(sparse_vec)

stopifnot(inherits(sparse_vec, "SparseNeuroVec"))
stopifnot(identical(dim(sparse_vec), dim(vec)))

This is the right move when absent voxels should be treated as missing support, not as literal zeros.

Where to go next

vignette("ChoosingBackends") for storage tradeoffs
vignette("coordinate-systems") for affine and orientation details
vignette("Resampling") for grid changes and reorientation
vignette("AnalysisWorkflows") for ROIs, searchlights, and split-map-reduce patterns