Monday, June 4, 2018

tutorial: plotting GIfTI images in R (knitr)

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 "blobs" are F-statistics from a GLM targeting right-hand button pushes, so the peaks in left motor areas are sensible. The underlay surfaces are the "inflated" images from the HCP 1200 Subjects Group Average Data release; download here. The code requires you to have underlay surfaces and statistics to overlay in perfect agreement (the same number of nodes, same triangle connections, etc.), but this is a general requirement for surface plotting or analysis.


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.


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, cache=TRUE>>=  
 # 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)  
   
 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"));   
   
   
 # adapted from gifti package gifti_map_val function.  
 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. I find facets=TRUE produce better images on windows 10, facets=FALSE better on linux.  
 plot.surface <- function(gii, ttl, pos.lims=c(NA,NA), neg.lims=c(NA,NA), set.facets=TRUE) {   
  # gii <- tmp; ttl <- "title"; pos.lims <- c(5,20); neg.lims <- c(NA,NA); set.facets <- TRUE;  
  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);
  # 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(cm.colors(64), heat.colors(64)[1:60]);  # both hot and cool colors  
   sc.text <- paste("\n 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 <- heat.colors(64)[1:60];   # heat.colors() truncated so top value not white - white doesn't show well on the underlay
   sc.text <- paste("\n 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 <- cm.colors(64);  
   sc.text <- paste("\n colors:", neg.lims[1], "to", neg.lims[2]);  
  }  
    
  # plot the medial and lateral views  
  triangle3D(tri=gii$coords, theta=90, phi=0, ltheta=90, lphi=0, bty='n', colkey=FALSE, zlim=c(-60,90), d=6, colvar=gii$values, col=cols,   
        clim=c.lims, facets=set.facets, resfac=0.01);    
  mtext(text=ttl, side=3, line=-2.5, cex=0.6);  
  triangle3D(tri=gii$coords, theta=270, phi=0, ltheta=270, lphi=0, bty='n', colkey=FALSE, zlim=c(-60,90), d=6, colvar=gii$values, col=cols,  
        clim=c.lims, facets=set.facets, resfac=0.01);   
  mtext(text=paste0("range: ", round(full.min,3), " to ", round(full.max,3), sc.text), side=1, line=-2.5, cex=0.5, adj=1);  
 }  
   
 @  
   
 \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  
 \vspace{0.2 cm}   
 <<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  
  }  
 }   
   
 @  
   
   
 \newpage  
 \noindent t-statistics: both positive (hot) and negative (cool) values. \par  
 \vspace{0.2 cm}   
 <<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);  # positive 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  
  }  
 }   
   
 @  
   
 \end{document}  
   

No comments:

Post a Comment