Creating 3D visualizations of MRI data: A brief guide

Procedure #1: Visualizing cluster maps in a glass brain

As a first exercise in visualizing MRI data in 3D, we will start with a statistical map. Depending on where your maps are coming from, you may need to apply a height threshold (i.e., t- or Z-critical) and/or a minimum cluster extend threshold (k). As a starting point and to make this guide more general and more reproducible, I will start with a statistical map obtained from NeuroSynth (http://www.neurosynth.org; Yarkoni et al., 2011), which will be in NIfTI format (Neuroimaging Informatics Technology Initiative; http://nifti.nimh.nih.gov). Briefly, NeuroSynth conducts automated meta-analyses across thousands of fMRI studies by calculating a frequency metric for how often specific terms are mentioned in the paper (e.g., “memory”, “emotion”) in relation to voxels reported in the results tables. See Yarkoni et al. (2011) for further details. As an example of how to obtain thresholded statistical maps from SPM analyses, see Appendix B.

For this example I searched the online version of NeuroSynth for the term “memory” (http://www.neurosynth.org/analyses/terms/memory/; see Figure 3A). I used the “forward inference” map as the example statistical map, along with the anatomical volume provided (obtained by clicking the download buttons displayed to the right of the layer names). For these examples, un-gzip the NIfTI volumes from NeuroSynth. Rename the forward inference map file from memory_pAgF_z_FDR_0̷.0̷1.nii to statmap.nii.

Figure 3. Images of the statistical map used in Procedure #1.

(a) Obtaining the statistical map from NeuroSynth.org. (b) Coronal slices of the thresholded activation clusters. (c) 3D renderings of the clusters from two different perspectives. (d) Stereoscopic anaglyph 3D rendering of the first perspective shown in panel C, to be viewed using red-blue 3D glasses.

Obtaining the anatomical ‘glass brain’ image. Since the anatomical 3D surface meshes are generally usable, in addition to outlining the steps for creating this glass brain volume, the resulting surface mesh file is also provided as Supplementary material (see Appendix A). While there is an abundance of anatomical volumes in normalized template space, here we will use the one provided on NeuroSynth (click “anatomical” where shown in Figure 2A).

Obtaining the thresholded cluster image. Before the map can be rendered in 3D, both the height and cluster-extent thresholds should be applied. In some fMRI analysis packages this can be output directly (e.g., see Appendix B). If this is the case, export the threshold cluster image and skip to section 1.3; if this is not the case, we will manually apply these thresholds ourselves. Here I will use examples of how to manually apply these thresholds using MATLAB (R2013a; The MathWorks Inc., Natick, MA) and SPM8 (http://www.fil.ion.ucl.ac.uk/spm/; Wellcome Department of Cognitive Neurology, UCL, London, UK), though other packages are able to do this as well. (See Madan, 2014, for an introductory guide to MATLAB.)

Manually applying the height and cluster-extent threshold is a bit cumbersome. Using the imcalc function in SPM, we can easily apply a height threshold to our volume by outputting a binary volume, where the voxel intensity statistic (i.e., from a t-, F-, or Z-statistic map) is above the threshold. Since our NeuroSynth memory image has a large number of highly significant clusters, we will threshold our statmap.nii to isolate the voxels where the statistic (Z-value) is above 12 (MATLAB code shown below). In the current case, lower thresholds yielded large clusters, which made the figure less interpretable (i.e., many regions comprised a single cluster, making it difficult to view the topology of the regions from the 3D view). When plotting results from your own fMRI study, you would likely use a threshold around 3 for the t- or Z-statistical map (corresponding to approximately p<.001).

>> spm_imcalc_ui('statmap.nii','statmapH.nii','i1>12');

To apply the cluster-extent threshold, we will use the nii_threshreslicecluster function (freely available from http://www.mccauslandcenter.sc.edu/CRNL/tools/spm8-scripts) to isolate clusters of voxels of at least a minimum volume of 400 mm³. Again, this value can be adjusted, and usually would be set higher than you would use for your statistical analyses, as the 3D rendering is intended more to provide a global view of the significant clusters, and is encumbered by the inclusion of many small clusters. A volume of 400 mm³ corresponds to 50 voxels where the voxel size is 2 mm-isotropic. The function can also apply height thresholds, but it thresholds rather than binarizes the image (i.e., converting it to a mask), which is not as useful for our current purposes.

>> nii_threshreslicecluster ('statmapH.nii','statmapH.nii',.5,40̷0̷);

The output file from this command will automatically be named rstatmapH.nii, rename it to statmapThresh.nii.

The height threshold should also be applied to the anatomical volume, anatomical.nii, but the cluster threshold is unnecessary. The resulting output file will be named anatomicalH.nii, rename it to glassbrain.nii.

Convert to VTK. To visualize the NIfTI volumes in 3D, we need to convert the voxel data into a 3D surface mesh in the VTK (Visualization ToolKit) format. Designed for anatomical tracing, ITK-SNAP includes this functionality. The simplest way to do so is to load each volume as both the main volume and as the segmentation volume. If you use the structural volume as the main volume and the statistical map as the segmentation, you may have issues with the bounding boxes not matching. Since we will move to another program with our 3D surfaces, it does not matter if the bounding boxes match or not.

Make sure that the volume is loaded correctly, as shown in Figure 2A. While ITK-SNAP can render 3D volumes, as shown in the bottom left portion of the screenshot shown in Figure 2A, its rendering options are limited. For instance, if you want to render several volumes in 3D from a consistent perspective/camera angle, ITK-SNAP is unable to accommodate; while ITK-SNAP can temporarily store camera information, this perspective information is lost if the program is closed or crashes, and it cannot be saved for later use nor can it be manually specified. As a result, it will be impossible to obtain the exact same camera angle. To rectify this shortcoming, we will make our 3D renderings in ParaView, which also has additional useful features. To export the meshes in VTK format from ITK-SNAP, use the menus to navigate to Segmentation, then “Export as Surface Mesh…”. Next, choose “Export meshes for all labels as a single scene” and save the file as a “VTK PolyData File”. In the current example there is only one surface mesh in each volume, but this is not always the case, such as in the ROI example discussed later. Note, it is possible to export volume data, rather than surface mesh data, as a VTK file in ITK-SNAP, but these files will not work with ParaView in the next step. If your VTK file does not work, double check that it was correctly exported as a surface mesh.

Repeat these steps for both the statistical map and anatomical volume.

Render in 3D. Start ParaView and open your two new VTK files within the same scene. ParaView can be a bit overwhelming at first, but it has many useful features for rendering and setting up the camera. With some adjustment of the colors and opacity for the two surfaces, it should be fairly easy to produce a set up in ParaView similar to Figure 2B. You can rotate the camera manually using the mouse, and can reset the camera position with the buttons labelled “+X” through “−Z”. When the scene state is saved, the camera position is preserved in the scene file, allowing you to easily load another statistical map at a later time. The scene as a whole can be saved by selecting “File” then “Save State…” (PVSM format). The final renderings produced here are shown in Figure 3C, corresponding to the series of coronal slices shown in Figure 3B. Renderings can be saved using either “File” then “Save Screenshot…” or “Export Scene…”. Screenshots will always be exported as raster (i.e., pixel) images, while ‘exported scenes’ are vector/polygon based. Note that exported PDFs can also be based on “rasterize 3D geometry” (there is a checkbox). If you are unsure what you require, a screenshot is likely sufficient, but do try and experiment to find out what settings best meet your needs, as this overview of ParaView’s functionality is far from comprehensive.

ParaView can also render stereoscopic 3D figures (e.g., anaglyph [red-blue], side-by-side) with a variety of 3D-compatible glasses options. An example of a red-blue stereoscopic render is shown in Figure 3D.

Procedure #2: Visualizing anatomical ROIs

Obtain ROI volume. For this example, I extracted several regions of the medial temporal cortex (hippocampus, amygdala, parahippocampal gyrus, fusiform gyrus) from the right hemisphere of the Hammers et al. (2003) maximum probability atlas (n30r83; http://biomedic.doc.ic.ac.uk/brain-development/index.php?n=Main.AdultMaxProb). Regions were extracted using the imcalc tool included in SPM8, such that each ROI corresponded to a unique intensity value (1=hippocampus, 2=amygdala, 3=parahippocampal gyrus, 4=fusiform gyrus):

>> spm_imcalc_ui({'Hammers_mith_atlas_n30r83_SPM5.nii'},'roiMTL.nii','(i1==1)*1 + (i1==3)*2 + (i1==9)*3 + (i1==15)*4');

The ROI is shown plotted over a structural volume in Figure 4A.

Figure 4. Images of the anatomical ROIs used in Procedure #2.

(a) Coronal slices of the anatomical ROIs. Panels B-E depict different 3D rendering settings of the ROIs, (b) points along the surface of the ROIs, (c) wireframe, (d) surface, and (e) surface rendering with specular. Note that the structural image used in panel A is from a different source than the traced ROIs, and thus they do not perfectly align.

Convert to VTK. As before, use ITK-SNAP to load the NIfTI volume and convert it to a VTK surface mesh. If you have multiple surfaces in the same volume, as we do here, be sure to select “Export meshes for all labels as a single scene” when exporting the surfaces.

Render in 3D. Start ParaView and load the VTK file, as done previously. As shown in Figure 4B–E, the volumes can be rendered as points, wireframes, and surfaces. Furthermore, many settings can be customized to adjust the rendering properties, such as the lighting/reflectance properties shown in Figure 4D and 4E.

ParaView can also create cameras that move over time, allowing for the generation of animations of the structures rotating. This can be done using the “Animation View” panel in the bottom-center of ParaView: select “Camera”, “Orbit”, and then “+”. The default settings for the camera positions are usually sufficient. If desired, the camera path can also be edited afterwards by inputting specific coordinates (the best way to preview the path is to simply press ‘play’ at the top and see how it looks). Even without rendering the animation itself, having a camera path allows for later reproduction of 3D renderings from the same camera positions.

Using a camera path, an animation can be rendered by going “File”, “Save Animation”. An example rendered video is shown in Movie 1. (Note, videos here were re-compressed with Handbrake [https://www.handbrake.fr; freely available for Windows and Mac] to reduce their file size).