NOTICE: this post was replaced by a newer version 13 March 2020. I will leave this post here, but strongly suggest you follow the new version instead.
In a previous tutorial I showed how to plot volumetric NIfTI images in R and knitr. Here, I give functions for plotting surface GIfTI brain images in R. If you're not familiar with surface images, this and this post may be useful. Please cite if you find this useful.
This uses the R gifti package (thank you, John Muschelli!) to read the GIfTI images, and the plot3D package for the plotting. It is possible to plot interactive 3d images in R (see for example, the gifti package vignette), but static images are more useful for the application I have in mind: batch-produced knitr files of surface (afni) GLM results.
Here is the first part of the knitr file produced in this tutorial; the full pdf can be downloaded here:
The gifti-plotting function I wrote displays the medial and lateral views, with hot colors for positive values and cool colors for negative. The function takes the positive and/or negative color scaling limits, with values closer to zero than the minimum value not shown, but those above the maximum plotted in the maximum color. This behavior is conventional for neuroimaging; for example, this is the same image, plotted in the Connectome Workbench with the same color scaling.
Full code for the demo is below the jump (or downloaded here). To run the demo you will need the two gifti overlay images, which can be downloaded here (left) and here (right), as well as the HCP 1200 underlay anatomies (see above). If you don't want to run the code within knitr the relevant functions can be separated from the latex sections, but you will need to set the image sizing and resolution.
Please let me know if you find this code useful, particularly if you develop any improvements!
UPDATE 20 July 2018: Fixed the plot.surface function to display the entire range (not truncated) in the lower right note.
UPDATE 7 September 2018: Another fix to plot.surface() to truncate the heat.colors() range so that the top values are plotted in yellow, rather than white, which wasn't very visible.
Also, Chris_Rorden posted a python and Surfice version on Neurostars.
UPDATE 9 November 2018: Added a section assigning arbitrary values to surface parcels and then plotting, matching this demo.
UPDATE 9 January 2019: Added a which.surface term to the plot.surface function, specifying whether plotting HCP or fsaverage5 surfaces. The default is HCP, for backward compatibility.
UPDATE 10 April 2019: Improved color ranges.
NOTICE: this post was replaced by a newer version 13 March 2020. I will leave this post here, but strongly suggest you follow the new version instead.
NOTICE: this post was replaced by a newer version 13 March 2020. I will leave this post here, but strongly suggest you follow the new version instead.
knitr (R) code for plotting gifti images; the same file can be downloaded here. \documentclass{article}
\addtolength{\oddsidemargin}{-1.25in}
\addtolength{\evensidemargin}{-1.25in}
\addtolength{\textwidth}{2.5in}
\addtolength{\topmargin}{-.875in}
\addtolength{\textheight}{1.75in}
\begin{document}
<<startup, echo=FALSE, message=FALSE, warning=FALSE>>=
# code written by Joset A. Etzel (jetzel@wustl.edu) and may be adapted, provided this source is cited.
# started 1 June 2018. Updated 9 November 2018 with the parcel-based plotting.
# Cognitive Control & Psychopathology Lab, Psychological & Brain Sciences, Washington University in St. Louis (USA)
library(gifti); # for readGIfTI
library(plot3D); # for triangle3D
rm(list=ls()); # clear R's memory
under.path <- "d:/localFiles/workbench/HCP_S1200_GroupAvg_v1/"; # path to the underlay surface files (surf.gii)
over.path <- "d:/temp/"; # path to the overlay files (func.gii, with statistics)
# the HCP 1200 subjects release surfaces for underlay image; see https://www.mail-archive.com/hcp-users@humanconnectome.org/msg05032.html
# https://db.humanconnectome.org/app/action/ChooseDownloadResources?project=HCP_Resources&resource=Workbench&filePath=HCP_S1200_GroupAvg_v1.zip
surf.L <- readGIfTI(file.path(under.path, "S1200.L.inflated_MSMAll.32k_fs_LR.surf.gii")); # L underlay surface
surf.R <- readGIfTI(file.path(under.path, "S1200.R.inflated_MSMAll.32k_fs_LR.surf.gii"));
# read in the gifti files with unique integers for each parcel; must be aligned with the underlay surface giftis (surf.L and surf.R)
# these gifti files extracted from the MMP cifti included with the HCP 1200 subjects release (same link as for surf.L and surf.R); see
# http://mvpa.blogspot.com/2017/11/tutorial-assigning-arbitrary-values-to.html.
# wb_command -cifti-separate D:\localFiles\Workbench\HCP_S1200_GroupAvg_v1\Q1-Q6_RelatedValidation210.CorticalAreas_dil_Final_Final_Areas_Group_Colors.32k_fs_LR.dlabel.nii COLUMN -metric CORTEX_LEFT D:\temp\mmpL.func.gii
# wb_command -cifti-separate D:\localFiles\Workbench\HCP_S1200_GroupAvg_v1\Q1-Q6_RelatedValidation210.CorticalAreas_dil_Final_Final_Areas_Group_Colors.32k_fs_LR.dlabel.nii COLUMN -metric CORTEX_RIGHT D:\temp\mmpR.func.gii
parcels.L <- readGIfTI(paste0(over.path, "mmpL.func.gii"));
parcels.R <- readGIfTI(paste0(over.path, "mmpR.func.gii"));
# so that this is a stand-alone demo, the two functions for plotting giftis are below. However, I suggest you save the code below as a separate
# file called giftiPlottingFunctions.R, then source this file in each knitr, rather than copying the functions into each knitr. The advantage of
# having a separate giftiPlottingFunctions.R file is that if you update the plotting function you'll only need to update it in that one file, not
# copy the change into multiple knitrs.
# source("d:/temp/giftiPlottingFunctions.R"); # preferred method; delete the lines below after making giftiPlottingFunctions.R.
####################### start giftiPlottingFunctions.R
# standard gifti image plotting functions, intended for calling when plotting in knitr files (or similar).
# see http://mvpa.blogspot.com/2018/06/tutorial-plotting-gifti-images-in-r.html for more explanation
# code written by Joset A. Etzel (jetzel@wustl.edu) and may be adapted, provided this source is cited. 1 June 2018.
# Cognitive Control & Psychopathology Lab, Psychological & Brain Sciences, Washington University in St. Louis (USA)
# last updated 10 April 2019 by Jo Etzel.
# colors to use for positive (warm) and negative (cool) values, generated using http://www.perbang.dk/rgbgradient/
cols.warm <- c("#E50800", "#E50C00", "#E51001", "#E61402", "#E61803", "#E61C03", "#E72004", "#E72405", "#E72806", "#E82D06",
"#E83107", "#E83508", "#E93909", "#E93D09", "#E9410A", "#EA450B", "#EA490C", "#EB4E0C", "#EB520D", "#EB560E",
"#EC5A0F", "#EC5E10", "#EC6210", "#ED6611", "#ED6A12", "#ED6E13", "#EE7313", "#EE7714", "#EE7B15", "#EF7F16",
"#EF8316", "#F08717", "#F08B18", "#F08F19", "#F19419", "#F1981A", "#F19C1B", "#F2A01C", "#F2A41C", "#F2A81D",
"#F3AC1E", "#F3B01F", "#F3B420", "#F4B920", "#F4BD21", "#F5C122", "#F5C523", "#F5C923", "#F6CD24", "#F6D125",
"#F6D526", "#F7DA26", "#F7DE27", "#F7E228", "#F8E629", "#F8EA29", "#F8EE2A", "#F9F22B", "#F9F62C", "#FAFA2D");
cols.cool <- c("#001CE5", "#001FE4", "#0123E4", "#0227E4", "#032BE3", "#032EE3", "#0432E3", "#0536E3", "#063AE2", "#063EE2",
"#0741E2", "#0845E2", "#0949E1", "#094DE1", "#0A50E1", "#0B54E1", "#0C58E0", "#0C5CE0", "#0D60E0", "#0E63E0",
"#0F67DF", "#106BDF", "#106FDF", "#1172DF", "#1276DE", "#137ADE", "#137EDE", "#1482DE", "#1585DD", "#1689DD",
"#168DDD", "#1791DD", "#1894DC", "#1998DC", "#199CDC", "#1AA0DC", "#1BA4DB", "#1CA7DB", "#1CABDB", "#1DAFDB",
"#1EB3DA", "#1FB6DA", "#20BADA", "#20BEDA", "#21C2D9", "#22C6D9", "#23C9D9", "#23CDD9", "#24D1D8", "#25D5D8",
"#26D8D8", "#26DCD8", "#27E0D7", "#28E4D7", "#29E8D7", "#29EBD7", "#2AEFD6", "#2BF3D6", "#2CF7D6", "#2DFBD6");
# adapted from gifti package gifti_map_val function. See notes in "D:\gitFiles_ccplabwustl\R01\Jo\for800msecTR\gifti_test.R"
gifti.map <- function (pointset, triangle, values) { # pointset <- surf.L$data$pointset; triangle <- surf.L$data$triangle; values <- over.L$data[[1]][,1];
tri <- as.vector(t(triangle) + 1); # transpose, unwrap + 1; +1 for 0-based (input) vs. 1-based (R)
out.vals <- values[tri[seq(from=1, to=length(tri), by=3)]]; # color triangle by FIRST point in each:
pointset <- pointset[tri, ];
return(list(coords=pointset, values=out.vals))
}
# surface plotting function
# gii is the object returned by gifti.map
# ttl is a string to display at the plot title (use "" for no title)
# pos.lims and neg.lims are the limits for the positive and negative (respectively) color scaling. Values from 0 to pos.lims[1] are not plotted,
# nor are values from neg.lims[2] to 0. Values larger than pos.lims[2] or smaller than neg.lims[1] are shown in the hottest or coolest color.
# print.scaling.label controls whether the color range label is written on the finished plot.
# which.surface indicates which underlay surface anatomy is being plotted, HCP or fsaverage5.
plot.surface <- function(gii, ttl, pos.lims=c(NA,NA), neg.lims=c(NA,NA), print.scaling.label=TRUE, which.surface="HCP") {
# gii <- tmp; ttl <- paste(lbls[i], "L", sess.ids[ssid]); pos.lims <- c(2,6); neg.lims <- c(-6,-2);
full.min <- min(gii$values, na.rm=TRUE); # store original minimum and maximum for later text label
full.max <- max(gii$values, na.rm=TRUE);
sc.text <- ""; # start scaling label with a blank string
# first fix the values for proper plotting
if (!is.na(pos.lims[1])) { # showing positive values, so fix too-big values (so don't disappear when plotting)
inds <- which(gii$values > pos.lims[2]);
if (length(inds) > 0) { gii$values[inds] <- pos.lims[2]; }
}
if (!is.na(neg.lims[1])) { # showing negative values, so fix too-small values
inds <- which(gii$values < neg.lims[1]);
if (length(inds) > 0) { gii$values[inds] <- neg.lims[1]; }
}
if (!is.na(pos.lims[1]) & !is.na(neg.lims[1])) { # both positive and negtive, so fix middle (zero-ish) values
inds <- which(gii$values > neg.lims[2] & gii$values < pos.lims[1]); # NA out middle values so don't get plotted
if (length(inds) > 0) { gii$values[inds] <- NA; }
c.lims <- c(neg.lims[1], pos.lims[2]); # set color scaling to smallest negative to biggest positive
cols <- c(rev(cols.cool), cols.warm); # both hot and cool colors
sc.text <- paste("colors:", neg.lims[1], "to", neg.lims[2], "&", pos.lims[1], "to", pos.lims[2]);
}
if (!is.na(pos.lims[1]) & is.na(neg.lims[1])) { # only positive values
c.lims <- pos.lims;
cols <- cols.warm;
sc.text <- paste("colors:", pos.lims[1], "to", pos.lims[2]);
}
if (is.na(pos.lims[1]) & !is.na(neg.lims[1])) { # only negative values
c.lims <- neg.lims;
cols <- rev(cols.cool);
sc.text <- paste("colors:", neg.lims[1], "to", neg.lims[2]);
}
# slightly different z-axis scaling and title positioning seems to work better for HCP and fsaverage5 (fMRIPrep) surface anatomies
if (which.surface == "HCP") { z.lim <- c(-60,90); ttl.line <- -2.5; }
if (which.surface == "fsaverage5") { z.lim <- c(-130,90); ttl.line <- -1; }
# actually plot
triangle3D(tri=gii$coords, theta=90, phi=0, ltheta=90, lphi=0, bty='n', colkey=FALSE, zlim=z.lim, d=6, colvar=gii$values, col=cols, clim=c.lims, facets=FALSE, resfac=0.01);
mtext(text=ttl, side=3, line=ttl.line, cex=0.6, adj=0);
triangle3D(tri=gii$coords, theta=270, phi=0, ltheta=270, lphi=0, bty='n', colkey=FALSE, zlim=z.lim, d=6, colvar=gii$values, col=cols, clim=c.lims, facets=FALSE, resfac=0.01);
if (print.scaling.label == TRUE) {
mtext(text=paste0("range: ", round(full.min,3), " to ", round(full.max,3), "\n", sc.text), side=1, line=-2.5, cex=0.5, adj=1);
}
if (print.scaling.label == "flipped") {
mtext(text=paste0(sc.text, "\nrange: ", round(full.min,3), " to ", round(full.max,3)), side=1, line=-2.5, cex=0.5, adj=1);
}
}
####################### end giftiPlottingFunctions.R
@
\section*{giftiPlottingDemo.rnw}
\noindent code written by Joset A. Etzel (\texttt{jetzel@wustl.edu}) on 1 June 2018 and released on \texttt{mvpa.blogspot.com}. \par
\noindent compiled \today\ \par
\vspace{0.3 cm}
\noindent This code demonstrates how to read a GIfTI brain image in to R and display in a standard format. The overlay images are t and F statistics from a single-subject afni GLM fitting button-pushes. Two buttons, each pressed with the right hand. \par
\vspace{0.3 cm}
\noindent F-statistics: only positive (hot) values. \par
<<code1, cache=TRUE, dev='png', dpi=150, echo=FALSE, fig.height=2, fig.width=7, fig.align='center'>>=
# NOTE: I usually do dev='pdf', but that makes these files HUGE. png is not vector, but resolution seems adequate.
layout(matrix(1:4, 1:4)); # one row with four images
par(oma=rep(0,4), mar=rep(0,4), bg='white');
# which indices of the gifti correspond to which statistics; see afni command 3dinfo
lbl.ids <- c("button1_Fstat", "button2_Fstat");
lbl.inds <- c(24, 31);
lims <- c(5,20); # color scaling for the F stats: don't show any F < 5; set maximum color to 20
fname.L <- file.path(over.path, "STATS_REML_L.func.gii");
fname.R <- file.path(over.path, "STATS_REML_R.func.gii");
if (file.exists(fname.L) & file.exists(fname.R)) {
over.L <- readGIfTI(fname.L); # overlay (STATS) for Left
over.R <- readGIfTI(fname.R); # and Right
for (i in 1:length(lbl.ids)) { # i <- 1;
tmp <- gifti.map(surf.L$data$pointset, surf.L$data$triangle, over.L$data[[lbl.inds[i]]][,1]); # get locations and colors in plotting format
plot.surface(tmp, paste(lbl.ids[i], "L"), lims); # left hemisphere
tmp <- gifti.map(surf.R$data$pointset, surf.R$data$triangle, over.R$data[[lbl.inds[i]]][,1]); # get locations and colors in plotting format
plot.surface(tmp, paste(lbl.ids[i], "R"), lims); # right hemisphere
}
}
@
<<code1b, cache=TRUE, echo=FALSE, dev='pdf', fig.height=0.5, fig.width=7.5, fig.align="center">>=
par(oma=rep(0,4), mar=rep(0,4), bg='white');
# positive colors only
plot(x=0, y=0, col='white', ylab="", xlab="", main="", bty='n', xaxt='n', yaxt='n'); # blank plot to hold legend
full.seq <- round(seq(from=lims[1], to=lims[2], length.out=length(cols.warm)),1);
show.pts <- round(seq(from=1, to=length(cols.warm), length.out=10));
legend('top', fill=cols.warm[show.pts], legend=full.seq[show.pts], bty='n', cex=0.8, title="", horiz=TRUE);
@
\newpage
\noindent t-statistics: both positive (hot) and negative (cool) values. \par
<<code2, cache=TRUE, dev='png', dpi=150, echo=FALSE, fig.height=2, fig.width=7, fig.align='center'>>=
# NOTE: I usually do dev='pdf', but that makes these files HUGE. png is not vector, but resolution seems adequate.
layout(matrix(1:4, 1:4)); # one row with four images
par(oma=rep(0,4), mar=rep(0,4), bg='white');
# which indices of the gifti correspond to which statistics; see afni command 3dinfo
lbl.ids <- c("button1#1_Tstat", "button2#1_Tstat");
lbl.inds <- c(21, 28); # which list index in the STATS gifti correspond to the entries in lbl.ids
plims <- c(2, 6); # positive color scaling: don't show any t < 2 but > 0; set maximum color to 6
nlims <- c(-6, -2); # negative color scaling: don't show any t > -2 but < 0; set maximum color to -6
fname.L <- file.path(over.path, "STATS_REML_L.func.gii");
fname.R <- file.path(over.path, "STATS_REML_R.func.gii");
if (file.exists(fname.L) & file.exists(fname.R)) {
over.L <- readGIfTI(fname.L); # overlay (STATS) for Left
over.R <- readGIfTI(fname.R); # and Right
for (i in 1:length(lbl.ids)) { # i <- 2;
tmp <- gifti.map(surf.L$data$pointset, surf.L$data$triangle, over.L$data[[lbl.inds[i]]][,1]); # get locations and colors in plotting format
plot.surface(tmp, paste(lbl.ids[i], "L"), plims, nlims); # left hemisphere
tmp <- gifti.map(surf.R$data$pointset, surf.R$data$triangle, over.R$data[[lbl.inds[i]]][,1]); # get locations and colors in plotting format
plot.surface(tmp, paste(lbl.ids[i], "R"), plims, nlims); # right hemisphere
}
}
@
<<code2b, cache=TRUE, echo=FALSE, dev='pdf', fig.height=0.75, fig.width=7.5, fig.align="center">>=
par(oma=rep(0,4), mar=rep(0,4), bg='white');
# positive and negative colors
plot(x=0, y=0, col='white', ylab="", xlab="", main="", bty='n', xaxt='n', yaxt='n'); # blank plot to hold legend
full.seq <- round(seq(from=plims[1], to=plims[2], length.out=length(cols.warm)),1);
show.pts <- round(seq(from=1, to=length(cols.warm), length.out=10));
legend('top', fill=cols.warm[show.pts], legend=full.seq[show.pts], bty='n', cex=0.8, title="", horiz=TRUE);
full.seq <- round(seq(from=nlims[1], to=nlims[2], length.out=length(cols.cool)),1);
show.pts <- round(seq(from=1, to=length(cols.cool), length.out=10));
legend('bottom', fill=rev(cols.cool)[show.pts], legend=full.seq[show.pts], bty='n', cex=0.8, title="", horiz=TRUE);
@
\vspace{2.0 cm}
\noindent Coloring MMP \textbf{parcels}, rather than vertices. This colors MMP parcels 1, 10, and 15 red (bilaterally), to match the demo at \texttt{http://mvpa.blogspot.com/2017/11/tutorial-assigning-arbitrary-values-to.html}. \par
<<code3, cache=TRUE, dev='png', dpi=150, echo=FALSE, fig.height=2, fig.width=7, fig.align='center'>>=
# NOTE: I usually do dev='pdf', but that makes these files HUGE. png is not vector, but resolution seems adequate.
layout(matrix(1:4, 1:4)); # one row with four images
par(oma=rep(0,4), mar=rep(0,4), bg='white');
# read in the values to plot. The file has one row for each parcel.
# code to create parcelNums.csv is given at http://mvpa.blogspot.com/2017/11/tutorial-assigning-arbitrary-values-to.html, and is
# out = repelem(0, 180); out(1) = 1; out(10) = 1; out(15) = 1; csvwrite('D:\temp\parcelNums.csv', out');
stats.L <- read.csv(paste0(over.path, "parcelNums.csv"), header=FALSE);
stats.R <- read.csv(paste0(over.path, "parcelNums.csv"), header=FALSE); # same in both hemispheres for this demo
# confirm that the parcel template and stats files have the same number of parcels (ROIs)
if (length(unique(parcels.L$data[[1]])) != nrow(stats.L)+1) { stop("length(unique(parcels.L$data[[1]])) != nrow(stats.L)+1"); } # +1 for 0 vertices
if (length(unique(parcels.R$data[[1]])) != nrow(stats.R)+1) { stop("length(unique(parcels.R$data[[1]])) != nrow(stats.R)+1"); }
# we want to assign a number to every parcel. The plotting functions require a number for every vertex. This uses
# the gifti parcels.L and parcels.R objects (assigned in the startup code block) to indicate which vertices make up each parcel.
# first, make vectors of 0s with entries for each vertex in each hemisphere
img.L <- rep(0, length(parcels.L$data[[1]]));
img.R <- rep(0, length(parcels.R$data[[1]]));
# iterate over the parcels, pulling the number for each parcel and assigning it to the vertices corresponding to that parcel
for (pid in 1:nrow(stats.L)) { # pid <- 1;
img.L[which(parcels.L$data[[1]] == (pid+180))] <- stats.L[pid,1]; # in MMP, left hemisphere vertices are 181:360
}
for (pid in 1:nrow(stats.R)) { img.R[which(parcels.R$data[[1]] == pid)] <- stats.R[pid,1]; }
# finally, plot
plims <- c(0.75, 2); # positive color scaling: don't show any 0 < t < 0.75; maximum hot color is 2
tmpL <- gifti.map(surf.L$data$pointset, surf.L$data$triangle, img.L); # turn the vector of vertices to a plotting structure
plot.surface(tmpL, "Left hemisphere, MMP", pos.lims=plims, neg.lims=c(NA,NA), print.scaling.label=FALSE);
tmpR <- gifti.map(surf.R$data$pointset, surf.R$data$triangle, img.R);
plot.surface(tmpR, "", pos.lims=plims, neg.lims=c(NA,NA), print.scaling.label=TRUE); # neg.lims=c(NA,NA) for no negative colors
@
\end{document}
No comments:
Post a Comment