intro to working with volume and surface brain data

When preparing to update my surface (gifti) tutorial knitr, I realized that some of the confusion I was trying to address wasn't due to the code so much as the general differences between volume and surface datasets. This post tries to fill in some of those gaps for people using my R scripts. Much of the information in this post is general, but there are many topics I am not going to cover at all - what I describe here is not the only way to work with fMRI brain data; there are many other file formats and software packages.

first, volumes.

Volumetric brain data files are first, both historically (all fMRI analyses were originally volumetric) and fundamentally: we (including our brains) are physically in three dimensions, as is the data that comes off the scanner (often as DICOM files, which neuroimagers usually promptly convert to NIfTI, such as via dcm2niix when making a BIDS dataset).

Consider a standard anatomical image, such as the HCP 1200 Group Average, which makes a good default underlay for adult MNI datasets (S1200_AverageT1w_81x96x81.nii.gz is a downsampled version used in my tutorials). Opening S1200_AverageT1w_81x96x81.nii.gz with MRIcron shows us slices  for each of the three planes (below left), and the schematic at right is how the data comes out of a NIfTI image:

This image's data is a big cube subdivided into a grid of little cubes (voxels; volumetric pixels), each 2.4 mm on a side. It has 81 voxels along the i axis, 96 along the j, and 81 along the k, for 81*96*81 = 629,856 voxels total. There is a number in every one of these 629,856 voxels: the voxels fill the entire cube, not just the brain part in the middle. Thus, when you plot slices through a NIfTI it hopefully makes a brain-shaped picture, but that's not guaranteed - you can put any shape you want into a NIfTI file. (And hence why you should always include a "does it look like a brain" QC step.)

The code below opens the NIfTI in R with the splendid RNifti package and shows us the same information as the MRIcron screen capture above: the voxel under the cursor (at 25, 56, 37; red scribble) has value 605.887 (green scribble). (Warning: these match since both R and MRIcron are 1-based; if using something 0-based like afni or python you will need to shift indices.)

 library(RNifti); # package for NIfTI image reading https://github.com/jonclayden/RNifti  
   
 fname <- paste0(in.path, "S1200_AverageT1w_81x96x81.nii.gz");  # https://osf.io/6hdxv/  
 if (file.exists(fname)) { under.img <- readNifti(fname); } else { stop(paste("missing:", fname)); }  
   
 str(under.img)  
 # 'niftiImage' num [1:81, 1:96, 1:81] 0 0 0 0 0 0 0 0 0 0 ...  
 # - attr(*, "pixdim")= num [1:3] 2.4 2.4 2.4  
 # - attr(*, "pixunits")= chr [1:2] "mm" "Unknown"  
 # - attr(*, ".nifti_image_ptr")=<externalptr>   
 # - attr(*, ".nifti_image_ver")= int 2  
   
 under.img[25,56,37]  
 # [1] 605.887  

Note that under.img is a 3d array of numbers with various attributes; its pixdim and pixunits match those in the MRIcron "NIfTI Header Information" dialog in the screen capture (open with Window -> Information). There are many more NIfTI header fields; the header information is what makes it possible for software like MRIcron to match the 3d array of numbers in the NIfTI with actual space. 

We don't need to implement every NIfTI header field or work out warping and realigning on our own, thankfully! The key is to think of the volumetric images as a 3d (or 4d, if functional timeseries - lots of 3d images stuck together along the 4th dimension) array of numbers with important properties attached. If  you have multiple NIfTI files you can treat them as normal arrays of numbers, so long as all the files have matching properties (same voxel sizes, orientation, etc. - compatible headers). Thankfully, all the main fMRI software programs have functions for making image properties match; here's my afni example.

why surfaces?

I wrote above that neuroimaging data (and our brains!) start in 3d volume space; why would we want to think about it any other way? Well, when doing fMRI the goal is often to analyze the cortical grey matter BOLD signal, and we know the human cortex is physically a folded sheet (cortical "ribbon"); a convoluted surface.

Figure 2 from doi:10.1007/s00429-012-0411-8, PMID: 22488096 (below) shows a cortical surface (lots more info in the paper and elsewhere); the blue scribble I added to pane c follows the grey matter ribbon.

In principle, if the preprocessing algorithms find this grey matter ribbon we can concentrate our analyses there and have better precision and stronger results. There are a lot of complications in finding and transforming the grey matter ribbon of any particular person into a surface (a few of my musings; Andrew Jahn's), but we won't explore any of those here. The explanations below use standard images from FreeSurfer; to get a surface version of a new dataset I recommend preprocessing with software (e.g., fmriprep) which will output gifti files (or whatever format you want) directly (don't, for example, preprocess completely as volumes and then convert to surfaces for analysis).

introduction to gifti files

A single NIfTI file is enough for an entire brain, but not a single gifti file: gifti files come in paired sets (see also). The "paired" is because each gifti file has info for one hemisphere only: for the whole cortex you need a pair of files, one left and one right. And the "sets" is because gifti files come in multiple types, labeled by the bit of the filename before the .gii. Here we'll only use two types: .surf.gii surface "underlay" images, and .label.gii parcel labels.

For my gifti knitr tutorial I shared several sets of paired gifti files; here we'll look at the Schaefer2018 parcellation fsaverage5 ones in R, using the gifti library's readGIfTI function. (I converted the original files to gifti using FreeSurfer's mris_convert function.)

First, read four files: a pair each of .surf.gii and .label.gii:

 library(gifti);  # for readGIfTI  https://github.com/muschellij2/gifti  
   
 # surface "underlay anatomy"  
 fname <- paste0(in.path, "fsaverage5_lh.inflated_avg.surf.gii");  # https://osf.io/nxpa7/  
 if (file.exists(fname)) { surf.L <- readGIfTI(fname); } else { stop("can't find fsaverage5_lh.inflated_avg.surf.gii"); }  
 fname <- paste0(in.path, "fsaverage5_rh.inflated_avg.surf.gii");  # https://osf.io/24y9q/  
 if (file.exists(fname)) { surf.R <- readGIfTI(fname); } else { stop("can't find fsaverage5_rh.inflated_avg.surf.gii"); }  
    
 # Schaefer parcellation "overlay" surfaces for each hemisphere   
 fname <- paste0(in.path, "Schaefer2018_400Parcels_7Networks_order_10K_L.label.gii");  # https://osf.io/7wk8v/  
 if (file.exists(fname)) { parcels.L <- readGIfTI(fname); } else { stop(paste("missing:", fname)); }  
 fname <- paste0(in.path, "Schaefer2018_400Parcels_7Networks_order_10K_R.label.gii");  # https://osf.io/7tms8/  
 if (file.exists(fname)) { parcels.R <- readGIfTI(fname); } else { stop(paste("missing:", fname)); }   
   

The next code looks a bit at the list objects made with readGIfTI(). For plotting we only need the $data fields, but I suggest investigating the other entries as well to gain a sense of the type of information in them.

 str(surf.L$data);  
 # List of 2  
 # $ pointset: num [1:10242, 1:3] -0.685 22.359 41.099 9.658 -36.773 ...  
 # $ triangle: int [1:20480, 1:3] 0 0 0 0 0 3 3 3 4 4 ...  
   
 str(surf.R$data)  
 # List of 2  
 # $ pointset: num [1:10242, 1:3] -5.74 8.43 43.5 12.6 -40.33 ...  
 # $ triangle: int [1:20480, 1:3] 0 0 0 0 0 3 3 3 4 4 ..  
   
 str(parcels.L$data)  
 # List of 1  
 # $ labels: int [1:10242, 1] 57 81 49 180 40 132 14 22 0 143 ...  
   
 str(parcels.R$data)  
 # List of 1  
 # $ labels: int [1:10242, 1] 61 81 46 110 111 200 21 95 104 119 ...  

Notice that the two surf objects each have two subpart arrays ($pointset and $triangle), while each parcel object only has a vector of integers. The parcel vectors' lengths are 10,242, which matches the number of rows in the surface pointset arrays; for plotting, the gifti data sizes must match like this. Why 10,242? That's the number of vertices per hemisphere in fsaverage5; other surface spaces have a different number of vertices.

Triangles and vertices? These terms come from how surface meshes can be represented in computer science: points (vertices, orange arrow) connected by edges of specific lengths, making lots of little triangles, as in this part of the FreeSurfer website brain I enlarged:

using and plotting giftis

The key idea when working with gifti images using my scripts is that we manipulate the vectors of vertex-wise data, not the .surf.gii files (nor any other spatial information). In a given surface space (here, freesurfer5) the location of each vertex as ordered in the 10,242-element vectors is fixed, so we can manipulate its value as needed, as shown in the next examples.

This code produces the image below. The left hemisphere is blank: all its vertices (vals.L) are 0. For the right I again started with a vector of 10,242 0s, but then changed its 1,800th entry to 1, producing a brain picture with a red splotch at the location of vertex 1800. 

 vals.L <- rep(0, 10242); # start with all vertices 0  
 map.L <- gifti.map(surf.L.infl$data$pointset, surf.L.infl$data$triangle, vals.L); # map vals.L onto left surface  
   
 vals.R <- rep(0, 10242); # start with all vertices 0  
 vals.R[1800] <- 1; # change vertex 1800 to 1  
 map.R <- gifti.map(surf.R.infl$data$pointset, surf.R.infl$data$triangle, vals.R);  # map vals.R onto right surface  
   
 # plot lateral and medial of both hemispheres  
 plot.surface(map.L, ttl="left hemisphere, blank (all vertices 0)", pos.lims=c(0.5,5.5), which.surface='fsaverage5');   
 plot.surface(map.R, ttl="right hemisphere, vertex 1800 in red", pos.lims=c(0.5,5.5), which.surface='fsaverage5');   

And a final example (more are in the knitr tutorial and other demo code): the next code plots (Schaefer 400x7) parcel 60 on the left hemisphere, and all 200 parcels on the right hemisphere. How to know which of the 10,242 vertices correspond to left hemisphere parcel 60? The key is parcels.L (created above from Schaefer2018_400Parcels_7Networks_order_10K_L.label.gii), which has a 10,242-element vector of parcel labels: we just need to find the 60s.

 vals.L <- rep(0, 10242); # start with all vertices 0  
 inds <- which(parcels.L$data[[1]] == 60); # find indices; where 60s are in the 10,242-element vector parcels5.L$data  
 vals.L[inds] <- 4;  # change indices in vals.L for parcel 60 to 4  
 map.L <- gifti.map(surf.L.infl$data$pointset, surf.L.infl$data$triangle, vals.L);  # map vals.L onto left surface  
 map.R <- gifti.map(surf.R.infl$data$pointset, surf.R.infl$data$triangle, parcels.R$data[[1]]); # map all parcels.R vertices onto right surface  
   
 # plot lateral and medial of both hemispheres  
 plot.surface(map.L, ttl="left hemisphere, parcel 60 in orange", pos.lims=c(0.5,5.5), which.surface='fsaverage5'); # plot lateral and medial, left  
 plot.surface(map.R, ttl="right hemisphere, parcels 1 to 200", pos.lims=c(0.5,200.5), which.surface='fsaverage5'); # plot lateral and medial, right  

volume and surface brain plotting knitr tutorial

Here is the second of my pair of posts introducing my updated knitr brain plotting tutorials. The first tutorial describes setting up RStudio for knitr compilation, with base R graphics examples - start there. This post adds a pair of brain plotting tutorials, one for volumetric images and the other for surfaces. The source (knitrIntro_NIfTI.rnw, knitrIntro_gifti.rnw) and image files needed for compilation can be downloaded from the blog osf site, "knitr tutorials" section. The compiled pdfs are included as well, at knitrIntro_NIfTI.pdf and knitrIntro_gifti.pdf.

The surface and volume examples are roughly parallel, covering plotting both continuous statistical overlays and parcellations. The example code includes assigning values to parcels (e.g., to show the results of an analysis), adjusting the plot appearance, and doing math with the images in R. Please read both the text in the pdf and the code comments.

The brain images in the knitrs (a few of which are above) are added with a pair of functions I wrote: plot.volume() and plot.surface(). These are updated versions of the functions in previous posts on this blog, with changes to improve the appearance of the surface images and make the color scaling parameters for volumes more similar to those for surfaces.

My intention is that the plotting functions and usage examples in these knitrs will enable beginners to quickly start making useful (and attractive!) knitr documents to summarize, explore, and describe their own analyses. Please let me know if you encounter any bugs, have a suggestion for an improvement, or have a new feature suggestion.

UPDATE 13 January 2025: I did a major update of knitrIntro_gifti.rnw and replaced the gifti screen capture in this post to match. See also a new blog post introducing volume and surface datasets.

UPDATE 31 March 2023: I uploaded corrected versions of niftiPlottingFunctions.R and knitrIntro_NIfTI.pdf to osf and the volumetric screenshot above. The correction was needed because of a bug in my plot.volume() function causing the negative colors to be reversed: instead of values near zero being dark blue, they were cyan. Major thanks to Tan Nguyen for spotting the unexpected coloring! If you downloaded a copy of niftiPlottingFunctions.R before 30 March 2023, please replace it with this version.

UPDATE 29 June 2021: Changed the volumetric screenshot to show the new (diverging) cool colors (and a legend); niftiPlottingFunctions.R and giftiPlottingFunctions.R now use the same color scales.

working with surfaces: musings and a suggestion

First, some musings about surface datasets. I'm not convinced that surface analysis is the best strategy for human fMRI, particularly task-based fMRI. I don't dispute that the cortex is a sheet/ribbon/surface; what I'm skeptical about is whether fMRI resolution is (or could be) anywhere close to good enough to accurately keep track of the surface over hours of scanning. (More grey matter musings here.)

I also have concerns about how exactly to do the surface interpolation: MRI collects volumes, from which surfaces must be interpolated, and there are multiple algorithms for finding the surface (in the anatomical images) and performing the transformation (including of the corresponding fMRI images). None of these are remotely trivial procedures, and I have a general bias against adding complexity (and experimenter degrees of freedom) to preprocessing pipelines/analyses.

For a practical issue, quality control and interpreting results on the surface is more difficult than in the volume. I prefer omitting masking during volumetric analysis, which allows checking if the "blobs" fall in an anatomically-sensible way (e.g., are the results stronger in the grey matter or the ventricles?). This type of judgment is not possible on the surface, because it only includes grey matter. For a more extreme example, in the image below the top row is from buggy code that randomized the vertex assignments, while the second row are results with a weak effect, but correct assignments (ignore the different underlay surfaces). The two images don't look all that different on the surface, but a parallel error in the volume would make salt-and-pepper all over the 3d frame: something clearly "not brain-shaped" and so easier to spot as a mistake.

Even with those reservations, I often need to use surface images (and the benchmark GLM results look reasonable, so it does seem to work ok); here are some suggestions if you need to as well.

Key: surfaces are surfaces (made up of vertices) and volumes are volumes (made up of voxels), do not transform back and forth. If you want surface GLM results, conduct GLMs on vertex timeseries; for volume GLM results, conduct GLMs on voxel timeseries.

It is possible to quickly make a surface from a volume (e.g., with a couple of wb_command -volume-to-surface-mapping commands), but that is a very rough process, and should be used only for making an illustration (e.g., for a poster, though even that's not great - see below), NOT for data that will be analyzed further. Instead, generating the surfaces should be done during preprocessing (e.g., with fmriprep), so that you have a pair of gifti timeseries files (one for each hemisphere) in your chosen output space (e.g., subject space, fsaverage5) along with the preprocessed 4d NIfTI file (e.g., subject space, MNI). (Note: I have found it easiest to work with surfaces as gifti files; cifti files can be separated into a pair of giftis and a NIfTI.)

Why this emphasis on keeping surfaces as surfaces and volumes as volumes? The main reason is related to the fuzziness and complexity in creating the surfaces. In both the fmriprep and HCP surface preprocessing pipelines vertices and voxels do not correspond at the end: you can't find a voxel timeseries to match each vertex timeseries (or the reverse). Instead, the vertex timeseries are a weighted mix of voxels determined to correspond to each bit of cortical ribbon, perhaps followed by parcel-constrained (or other) smoothing. It may be theoretically possible to undo all that weighting and interpolation, but no documented method or program currently exists, as far as I know. More generally, information is lost when transforming volumes to surfaces (and the reverse), so it should be done as rarely as possible, and the "rarest" possible is once (during preprocessing).

What can you do with giftis? Pretty much everything you can do with NIfTIs (plotting, extracting stats, running GLMs, etc.), though a bit of tweaking may be needed. For example, we use afni for GLMs, and many of its functions can now take gifti images. I use R for nearly all gifti plotting and manipulation, as in this tutorial, which includes my plotting function, and turning vertex vectors into surface images.

This snippet of code illustrates working with giftis in R (calculating vertex-wise mean, standard deviation, and tsnr for QC): in.fname is the path to a gifti timeseries (such as produced by fmriprep). It is read with readGIfTI from the gifti library, than rearranged into a vertex x frame matrix. Once in that form it's trivial to calculate vertex-wise statistics, which this code writes out in plain text; they could also be plotted directly.

      in.img <- readGIfTI(in.fname)$data;  # just save the vertex timeseries part  
      in.tbl <- matrix(unlist(in.img), nrow=length(in.img), byrow=TRUE); # now vertices in rows, timepoints in columns  
   
      out.tbl <- array(NA, c(ncol(in.tbl), 3));  # one row for each vertex  
      colnames(out.tbl) <- c("sd", "mean", "tsnr");  
   
      out.tbl[,1] <- apply(in.tbl, 2, "sd");  
      out.tbl[,2] <- apply(in.tbl, 2, "mean");  
      out.tbl[,3] <- out.tbl[,2]/out.tbl[,1];  # tsnr: mean/sd  
      write.table(out.tbl, out.fname);  

I try to carry this "surfaces as surfaces and volumes as volumes" strategy throughout the process, to keep the distinction obvious. For example, if a parcel-average timecourse was derived from surface data, show it next to a surface illustration of the parcel. If a GLM was run in the volume, show the beta maps on volumetric slices. Many people (myself included!) have made surface versions of volume data for a final presentation because they look "prettier". This was perhaps acceptable when only volumetric analyses were possible - we knew the surface image was just an illustration. But now that entire analyses can be carried out in the surface or volume I think we should stop this habit and show the data as it actually is - volumes can be "pretty", too!

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.

Connectome Workbench: making a surface version of a volumetric image

This is a tutorial for making a surface version of a volumetric (NIfTI) image (e.g., a ROI) for visualizing with the Connectome Workbench. This post replaces a tutorial I wrote back in 2013, and assumes a bit of familiarity with Workbench. If you're using Workbench for the first time I suggest that you complete the tutorial in this post before trying this one.

First, a warning: as advertised, the steps in this post will make a surface version of a volumetric image. However, the surface version will be an approximation, and likely only suitable for visualization purposes (e.g., making an illustration of a set of ROIs for a talk). If you have an entire dataset that you want to prepare for surface analysis (e.g., running GLMs), you need different procedures (e.g., SUMA, FreeSurfer). Again, I suggest the directions in this post (wb_command -volume-to-surface-mapping) be used cautiously, for quick visualizations, and accompanied by careful confirmation that the mapping produced a reasonable result.

needed files

Before we can make a surface version of a volumetric image, we need to know what it's aligned to, so that we can pick the proper corresponding surface template. Recall that gifti surface files (*.surf.gii) are sort of the underlay anatomy for surface images (e.g., in my Getting Started post on Workbench we load surf.gii files to get a blank brain), so we'll need gifti surface files that will work with our volumetric image (and to serve as the underlay when we're ready to plot the converted volumetric ROI).

For this demo, we'll use fakeBrain.nii.gz (it should let you download without signing in; this is the same NIfTI as shown in other posts), which is aligned to MNI. One MNI dataset with the necessary surface files is the HCP 1200 Subjects Group Average Data; this post describes the files and gives the direct ConnectomeDB download link.

The HCP 1200 Subject Group Average download contains multiple *.surf.gii for each hemisphere, including midthickness, pial, and very_inflated. We can use any of these for visualization in Workbench, but which we pick for the volume to surface conversion does make a difference in what the resulting surface image will look like. It seems best to start with the midthickness surface for the conversion, then try others if the projection seems off.

using wb_command

The wb_command -volume-to-surface-mapping function does the conversion. wb_command.exe (on my Windows machine; file extension may vary) should be in the same directory as wb_view.exe, which you use to start the Workbench GUI. Don't double-click wb_command.exe - it's a command line program. Instead, open up a command prompt and navigate to the directory containing wb_command.exe (on my machine, /bin_windows64/). If you  type wb_command at the prompt it should print out some version and help information; if not, check if you're in the correct directory, and try ./wb_command if you're on a linux-type system.

Now we're ready: we give the function our input NIfTI (fakeBrain.nii.gz), our surface gifti (S1200.L.midthickness_MSMAll.32k_fs_LR.surf.gii), the output file we want it to make (demoL.shape.gii), and the options for it to use (-trilinear). Since surface gifti files are just for one hemisphere, we have to do the command twice, once for each. (I included the full path to each file below; update for your machine.)

 wb_command -volume-to-surface-mapping d:/temp/fakeBrain.nii.gz d:/Workbench/HCP_S1200_GroupAvg_v1/S1200.L.midthickness_MSMAll.32k_fs_LR.surf.gii d:/temp/demoL.shape.gii -trilinear  
 wb_command -volume-to-surface-mapping d:/temp/fakeBrain.nii.gz d:/Workbench/HCP_S1200_GroupAvg_v1/S1200.R.midthickness_MSMAll.32k_fs_LR.surf.gii d:/temp/demoR.shape.gii -trilinear  

We now have demoL.shape.gii and demoR.shape.gii, surface versions of fakeBrain.nii.gz, which can be viewed in Workbench (or other gifti-aware programs). Check the surface projection carefully: does the ROI align properly with the anatomy? If not, try a different .surf.gii or fitting option (e.g., -enclosing) in the wb_command call; these can make a big difference.

Below the jump are the surface images from the above commands, plotted on the S1200 midthickness in Workbench.

assigning arbitrary values to surface parcels: method 1

This tutorial describes a method for plotting arbitrary per-parcel values (such as from an analysis) on the surface. For example, let's say I want to display MMP parcels 1, 10, and 15 (only) in red, or (more usefully) to assign continuous numbers to each parcel, and then display parcels with larger numbers in hotter colors.

This post describes a method using matlab and creating GIfTI files; see the next post for a method using R, wb_command functions, and creating a CIFTI file. Both methods work, but one or the other might be more convenient in particular situations.

I assume that you have a surface version of the parcellation to use as a template. For example, the MMP parcellation was released in CIFTI format as part of the HCP 1200 release, and the Gordon (2014) parcellation can be downloaded here.

I'll be using the MMP in this example; if you want to follow along, download a copy of the S1200 Group Average Data Release; I put mine at d:/Workbench/HCP_S1200_GroupAvg_v1/. The MMP template is named Q1-Q6_RelatedValidation210.CorticalAreas_dil_Final_Final_Areas_Group_Colors.32k_fs_LR.dlabel.nii. (If you're not sure how to use the files in the S1200 release, try this tutorial to get started.)

I have the values to assign to the parcels in a text file with 180 lines (one line for each MMP parcel). For this tutorial, let's do the simple example of assigning a color to parcels 1, 10, and 15 only. An easy way to do this is to make a text file with 1s on these rows and 0s on all the other rows. I prefer R, but since the GIfTI library is in matlab, here's matlab code for making the little text file:

 out = repelem(0, 180);  
 out(1) = 1;  
 out(10) = 1;  
 out(15) = 1;  
 csvwrite('D:\temp\parcelNums.csv', out')  

The MMP template is in CIFTI format, but we can extract GIfTI files for each hemisphere using wb_command cifti-separate:

 wb_command -cifti-separate D:\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:\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  

This matlab code reads in the text file with the new parcel values and the GIfTI templates, then writes GIfTI files with the new parcel values substituted for the parcel label numbers:

 addpath 'C:/Program Files/MATLAB/gifti-1.6';  % so matlab can find the library  
 mmpL = gifti('D:\temp\mmpL.func.gii');  % load the left side gifti MMP atlas  
 mmpR = gifti('D:\temp\mmpR.func.gii');  % and the right side MMP atlas  
 newvals = csvread('D:\temp\parcelNums.csv');  % 180 integers; new value for each parcel  
 Lout = mmpL;  % output gifti  
 Lout.cdata(:,1) = repelem(0, size(mmpL.cdata,1));  % replace the values with zeros  
 Rout = mmpR;  
 Rout.cdata(:,1) = repelem(0, size(mmpR.cdata,1));   
 for i=1:180  % i = 1;  
   inds = find(mmpR.cdata == i);  % find vertices for parcel i  
   Rout.cdata(inds,1) = newvals(i);  % put the new value in for this parcel's vertices  
   inds = find(mmpL.cdata == (i+180)); % in MMP, left hemisphere vertices are 181:360  
   Lout.cdata(inds,1) = newvals(i);    
 end  
 save(Lout,'D:\temp\plotL.func.gii','Base64Binary');  % save the gifti  
 save(Rout,'D:\temp\plotR.func.gii','Base64Binary');  

You can now to plot these GIfTIs in Workbench (see this post if you're not sure how); I plotted them on the S1200 Group Average (MNI) anatomy:

I clicked to put a marker in the left visual parcel. The value at this vertex is 1, as assigned (green lines). I loaded in the MMP atlas as well (blue lines), so it tells me (correctly!) that the marker is in the L_V1_ROI.

UPDATE 9 November 2018: I added making this same gifti on the fly and plotting in a knitr to my gifti plotting demo.

NIfTI, CIFTI, GIFTI in the HCP and Workbench: a primer

The HCP is releasing preprocessed data in both volumetric NIfTI and surface/volumetric CIFTI formats. Working with the HCP files, or doing much of anything with the Workbench, requires navigating through a plethora of .*.nii and .*.gii files. In this post I'll explain why we need all these files, and how they relate to each other. Disclaimer: I'm writing this as a primer from the viewpoint of someone familiar with volumetric fMRI data analysis; it is not at all a full description of everything the files can be used for. Also, though I'm referring to the HCP and Workbench, these file formats are used by other projects and software.

For a starting point, consider how we work with volumetric NIfTI files. Neuroimagers often think about volumetric NIfTI files as storing functional data in a 4d matrix (x,y,z, and time). Libraries such as  RNifti make reading NIfTI files fairly easy: they create a 3d or 4d matrix of voxel values, plus a object with the header information.

While you can get an idea of the anatomy by looking at slices of the 4d functional data matrix, analyses generally rely on having a 3d matrix of anatomical data (binary mask of regions, anatomic scan, etc) perfectly aligned to the 4d functional data. So, the 4d NIfTI file doesn't contain everything we need: we get some alignment information out of the header (qfactor, etc), but also need the registered 3d anatomical data. For a concrete example, I had to provide two files for the little ROI-based analysis demo: the dataset (4d NIfTI with preprocessed BOLD) and the ROI mask (binary 3d NIfTI showing the voxels corresponding to the anatomical region of interest), plus stating that the dataset was normalized to the MNI anatomical atlas (so that we can overlay the data on the correct anatomical template).

Now, on to CIFTI. CIFTI-2 files follow the NIfTI-2 file format specification (CIFTI-2 is a "flavor" of NIfTI-2, so both use the *.nii file extension), and both consist of a data matrix and headers. In the case of the HCP data, the functional timecourses are in the data matrix part of *.dtseries.nii CIFTI files. Like NIfTI volume files, the CIFTI file contains information about where voxels are, though this information is stored in a different place (in the extension containing the CIFTI XML). But, paralleling how you need an anatomic file to figure out exactly where the voxels in a volumetric NIfTI lie, you need other files (not just the CIFTI) to tell you where the surface vertices lie, and how they're connected (the "triangles", etc). Aside: While I wrote "surface vertices" in this paragraph, note that the HCP CIFTIs store both surface vertices (for the cortical sheet) and volumetric voxels (for sub-cortical structures).

These "other files" are not a single file but multiple; as many as necessary. Having all of these files is akin to having multiple ROI files available for an analysis: you won't use each ROI in each analysis, just the ones corresponding to the anatomical area (or whatever) you need for a particular test. The "other files" for the HCP are not just ROIs, but can also be underlying anatomy at different inflation levels, maps of tissue types, etc.

For example, at left is a screenshot showing some of the "other files" provided for each HCP person in the released datasets. These files are from /100307_Q3/MNINonLinear/Native/: the maps are in subject space. Many files with similar names are in /100307_Q3/MNINonLinear/fsaverage_LR32k/: maps of the same structures/types, but aligned to the MNI template anatomy (specifically, the 32k Conte69 mesh, see page 112 of Glasser, et. al 2013).

And now we're encountering GIFTI files: many of the "other files" are in GIFTI format, with the extension .*.gii. The naming of the "other files" (the last bit before the .gii) in the HCP tends to follow the CARET conventions, and gives a hint as to what sort of information they contain:

*.surf.gii, "gifti surface files", contain only vertex coordinates and triangles (which vertices are connected). The HCP *.surf.gii files are mostly structures that you might want to overlay data onto, such as 100307.L.inflated.native.surf.gii (left hemisphere, inflated) and 100307.L.midthickness.native.surf.gii.(left hemisphere, not inflated at all, but rather halfway through the thickness of the cortical ribbon).

*.func.gii and *.shape.gii, "metric files", contain data values for every vertex. Essentially, these are data arrays whose indices correspond to a surface file - you need a matching surface file to know where in the brain to put the data stored in a metric file. For example, a metric file from the HCP release is 100307.L.corrThickness.native.shape.gii: the cortical thickness at each vertex.

For an example of how these files work together, my tutorial on plotting a NIfTI image with the Workbench uses the wb_command -volume-to-surface-mapping program to create .shape.gii files aligned to Conte69.*.midthickness.32k_fs_LR.surf.gii. The data from the volumetric NIfTI (e.g. searchlight accuracies at each voxel) is stored (by vertex) in .shape.gii files, but a shape.gii file by itself isn't enough to plot the data properly on a surface: you need an aligned .surf.gii file as well. Paralleling how you need an aligned anatomy to properly overlay a volumetric NIfTI ROI, you need an aligned surf.gii to know how to properly locate the data from a metric file.

Whew! Hopefully this primer helps explain why so many files are released with the HCP data, and a bit about how they work together. For additional information see the Workbench Glossary, as well as Glasser, et. al 2013. If you've found any references particularly useful that I haven't already linked to, please send them along and I'll add links.

I want to end this post with a BIG thank you to Tim Coalson, who patiently (and repeatedly) walked me through these file types and how they relate to each other.

UPDATE 17 November 2016: Emma Robinson describes the various surface files, including the MSM related ones.

UPDATE 20 March 2026: Reviewing this old post, I think it still serves as a general introduction to the three file types, though many of its examples are outdated. My brain plotting knitr tutorials have much more current examples of reading and plotting both surface (gifti) and volume (nifti) files, and should be useful even if you aren't interested in the knitr aspect. There are now R packages like ciftiTools which can directly work with cifti files, and many people use ciftis for non-HCP datasets. Personally, I find it more straightforward to use giftis for surfaces and niftis for volumes than the combined ciftis, but converting is easy, for example with wb_command -cifti-separate.