In bbuchsbaum/neurosurf: Data structures and IO for surface-based neuroimaging data

Load example low desnity brain surface meshes

The neurosurf package includes some low density surfaces that are useful for testing and demonstrating package features. We first load left and right "inflated" and "smoothwm" surfaces that are included in the extdata package folder.

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(rgl)
library(neuroim2)
library(neurosurf)

inflated_lh_asc <- system.file("extdata", "std.8_lh.inflated.asc", package="neurosurf")
inflated_rh_asc <- system.file("extdata", "std.8_rh.inflated.asc", package="neurosurf")

white_lh_asc <- system.file("extdata", "std.8_lh.smoothwm.asc", package="neurosurf")
white_rh_asc <- system.file("extdata", "std.8_rh.smoothwm.asc", package="neurosurf")

Smoothing surface geometry

It is is often useful to smooth a surface mesh for visualization and other purposes. neurosurf uses the Rvcg package to smooth surface meshes using various smoothing approaches. Below we show smoothed surfaces using the Taubin and HCLaplace methods.

white_surf <- read_surf(white_lh_asc)
white_rh_surf <- read_surf(white_rh_asc)

open3d()
layout3d(matrix(1:4, 2,2), heights=c(1,3))
t1 <- text3d(0,0,0,"Taubin smooth, lambda=.5"); next3d()
white_surf1 <- smooth(white_surf, type="taubin", lambda=.8)
p1 <- plot(white_surf1); next3d()

t2 <- text3d(0,0,0,"HCLaplace smooth, delta=.2, iteration=5"); next3d()
white_surf2 <- smooth(white_surf, type="HCLaplace", delta=.2, iteration=5)
p2 <- plot(white_surf2)
rglwidget()

Adding a color layer to a surface mesh

Here we add a colors to the surface that are a function of each surface node's x coordinate. When then map these values to a rainbow color map.

mfrow3d(1, 2, byrow = TRUE)
xvals <- white_surf1@mesh$vb[1,]
p <- plot(white_surf2, vals=xvals, cmap=rainbow(255))
next3d()
p2 <- plot(white_surf2, vals=xvals, cmap=heat.colors(255))
rglwidget()

Showing an activation map overlaid on a surface mesh

We will plot surface in a row of 3. We generate a set of random values and then smooth those values along the surface to approximate a realistic activation pattern.

In the first column we display all the values in the map. Next we threshold all values between (-2,2). In the last panel we additionally add a cluster size threshold of 30 nodes.

open3d()
mfrow3d(1, 3, byrow = TRUE)
vals <- rnorm(length(nodes(white_surf2)))
surf <- NeuroSurface(white_surf2, indices=1:length(vals), data=vals)
ssurf <- smooth(surf)
p <- plot(geometry(ssurf), vals=values(ssurf), cmap=rainbow(100), irange=c(-2,2))

next3d()
comp <- conn_comp(ssurf, threshold=c(-.2,.2))
p2 <- plot(geometry(ssurf), vals=values(ssurf), cmap=rainbow(100), irange=c(-2,2), thresh=c(-.2, .2))

next3d()
csurf <- cluster_threshold(ssurf, size=30, threshold=c(-.2,.2))
p2 <- plot(csurf, cmap=rainbow(100), irange=c(-2,2), thresh=c(-.2, .2))
rglwidget()

Showing two hemisperes in same scene

open3d()

curv_lh <- curvature(white_surf2)
white_rh_surf2 <- smooth(white_rh_surf, type="HCLaplace", delta=.2, iteration=5)
curv_rh <- curvature(white_rh_surf2)

p <- plot(white_surf2, bgcol=curv_cols(curv_lh), viewpoint="posterior")
p <- plot(white_rh_surf2,bgcol=curv_cols(curv_rh), viewpoint="posterior")
rglwidget()