Monday, June 4, 2018

tutorial: plotting GIfTI images in R (knitr)

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 "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.

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