Beyond reading and writing, ieegio can also map data between volumes and surfaces, chain coordinate transforms, and compare regions of interest; those are covered at the end of this article.

To start, please load ieegio. This vignette uses sample data which requires extra download.

Volume files

ieegio::read_volume and ieegio::write_volume provides high-level interfaces for reading and writing volume data such as MRI, CT. fMRI, etc.

Each volume data (NIfTI, MGH, …) contains a header, a data, and a transforms list.

volume <- read_volume(nifti_file)
volume

The transforms contain transforms from volume (column, row, slice) index to other coordinate systems. The most commonly used one is vox2ras, which is a 4x4 matrix mapping the voxels to scanner (usually T1-weighted) RAS (right-anterior-superior) system.

Accessing the image values via [ operator. For example,

volume[128, , ]

Plotting the anatomical slices:

par(mfrow = c(1, 3), mar = c(0, 0, 3.1, 0))

ras_position <- c(-50, -10, 15)

ras_str <- paste(sprintf("%.0f", ras_position), collapse = ",")

for (which in c("coronal", "axial", "sagittal")) {
  plot(x = volume, position = ras_position, crosshair_gap = 10,
       crosshair_lty = 2, zoom = 3, which = which,
       main = sprintf("%s T1RAS=[%s]", which, ras_str))
}

Surface files

Reading surface file using read_surface supports multiple data types

library(ieegio)
# geometry
geometry <- read_surface(geom_file)

# measurements
measurement <- read_surface(shape_file)

# time series
time_series <- read_surface(ts_file)

You can merge them to a single object, making an object with multiple embedding data-sets:

merged <- merge(geometry, measurement, time_series)
print(merged)

Plot the surfaces in 3D viewer, colored by shape measurement

# plot the first column in measurements section
plot(merged, name = list("measurements", 1))

Plot the normalized time-series data

ts_demean <- apply(
  merged$time_series$value,
  MARGIN = 1L,
  FUN = function(x) {
    x - mean(x)
  }
)
merged$time_series$value <- t(ts_demean)
plot(
  merged, name = "time_series",
  col = c(
    "#053061", "#2166ac", "#4393c3",
    "#92c5de", "#d1e5f0", "#ffffff",
    "#fddbc7", "#f4a582", "#d6604d",
    "#b2182b", "#67001f"
  )
)

Streamline files

Reading streamlines via universal entry function read_streamlines

trk <- read_streamlines(trk_file, half_voxel_offset = TRUE)
tck <- read_streamlines(tck_file)
tt <- read_streamlines(tt_file)

To obtain the streamline data

message("Total number of streamlines: ", length(trk))

head(trk[[1]]$coords)

To preview the streamline data

pal <- colorRampPalette(c("navy", "grey", "red"))
plot(trk, col = pal(length(trk)))

To write the streamlines, for example, write tck file to trk file:

# Create a temporary file
tfile <- tempfile(fileext = ".trk")
write_streamlines(x = tck, con = tfile)

Surface annotations from `AFNI`/`SUMA`

NIML datasets (file names ending with .niml.dset) are read by the same read_surface entry point. The data type is resolved from the dataset itself: datasets carrying a label table are read as annotations, the rest as measurements.

std141_geometry <- read_surface(std141_geom_file)

annotation <- read_surface(niml_file)
annotation

The sample geometry and annotation live on the same std.141 mesh, so they merge directly:

labeled <- merge(std141_geometry, annotation)
labeled

plot(labeled, name = "annotations")

If you need the raw NIML element tree rather than a surface object, use the low-level io_read_niml together with niml_find.

Regions of interest

A region of interest is described first and computed later. as_ieegio_roi records the criteria without applying them, and resolve_roi_as carries them out, returning geometry in world (RAS) coordinates.

atlas <- read_volume(atlas_file)

# describe: which voxels count as the region
roi <- as_ieegio_roi(atlas, threshold_lb = 1, threshold_ub = 5)
roi

# compute: turn that description into geometry
resolve_roi_as(roi, "pointcloud")

Because both sides are resolved to world coordinates first, regions stored in different ways can be compared directly. Here the whole atlas is tested against the facial-nerve tracts read earlier:

overlap <- detect_roi_overlap(
  as_ieegio_roi(atlas, threshold_lb = 1),
  trk,
  radius = 2
)
overlap

The result carries the annotated streamlines back in overlap$annotated, so each tract knows whether it reached the region, and overlap$hit_ratio reports the proportion that did.

Volume to surface

volume_to_surface turns a mask or a set of atlas labels into a smoothed mesh:

volume_to_surface(atlas, threshold_lb = 1, threshold_ub = 5)

Coordinate spaces and transforms

new_space names a coordinate space, and surface_to_surface moves a surface into it, recording the target on the surface transform list.

mni <- new_space("MNI152", orientation = "RAS")
mni

surface_to_surface(
  geometry,
  space_from = "scanner",
  space_to = mni,
  transform = diag(1, 4)
)

Transforms themselves can come from other tools: io_read_ants_transform reads ANTs affine and displacement field transforms, io_read_flirt_transform reads FSL FLIRT matrices, and transform_flirt2ras converts a FLIRT transform into world (RAS) coordinates.

Read imaging data