Example analyses and use of functions
Reconstruction
For this paper, we will not run the analysis on a subject, but rather explore the output results for a subject included in the Freesurfer installation for reproducibility for the user. In particular, in the default Freesurfer subjects directory, there is a subject named “bert”. If we were to run all the analyses, we would use the code (described below):
We see the result of this output in the “bert” directory, which includes a series of sub-directories:
[1] "label" "mri" "scripts" "stats" "surf" "tmp" "touch"
[8] "trash" We will explore the results in “mri”, which contain imaging data, “stats”, which containing statistics based on structures of the brain, and “surf”, which contain the surface and curvature output from the Freesurfer processing.
MRI conversion: the function
The typical output format of brain volumes from Freesurfer is MGH/MGZ format, which is explained here: . As NIfTI formats are one of the most common formats and has been the common format for analysis in the and packages, it is useful to convert these files to a NIfTI format to read into R. The Freesurfer function will be used for that conversion. Here we will use the T1-weighted image from the “bert” subject and convert it to NIfTI, and read it into R:
library(freesurfer)
if (have_fs()) {
bert_dir = file.path(fs_subj_dir(), "bert") # subject directory
t1_mgz = file.path(bert_dir, "mri", "T1.mgz") # mgz file
t1_nii_fname = tempfile(fileext = ".nii.gz") # temporary NIfTI file
freesurfer::mri_convert(t1_mgz, t1_nii_fname) # conversion
img = neurobase::readnii(t1_nii_fname) # read in outputs
}As this is a commonly-used process, we have wrapped these two steps
into the readmgz and readmgh functions, which
combine the mri_convert and readnii functions.
Here we show that these steps are equivalent to the readmgz
function:
[1] TRUENow that we have the image in R, we can plot it using the standard
plotting tools for nifti objects:
The result is in figure \(\ref{fig:mri_plot}\), which contains 3 slices of the head: axially, viewing the brain from the top of the head (top left), sagittally, viewing the brain from the right side (top right) and coronally, viewing the brain from the back of the head (bottom left).
Note, the image is not stored in the right/posterior/inferior (RPI)
orientation which is assumed when displaying using the function http://www.grahamwideman.com/gw/brain/orientation/orientterms.htm.
We can use the rpi_orient function in (version \(\geq\) 2.4.0) or fslswapdim to
reorient the image to the assumed orientation.
if (requireNamespace("fslr", quietly = TRUE)) {
fsl_check = fslr::have_fsl()
} else {
fsl_check = FALSE
}
if (have_fs()) {
if (fsl_check) {
L = fslr::rpi_orient(img)
reoriented_img = L[["img"]]
} else {
reoriented_img = img
}
}We see that this function puts this image in the RPI orientation,
which matches the assumed orientation for ortho2:
The result is in figure \(\ref{fig:mri_plot2}\), which changes the views in reference to which panel they are located and matches the orientation markers assumed by .
Bias-field correction: the function
MRI images typically exhibit good contrast between soft tissue classes, but intensity inhomogeneities in the radio frequency field can cause differences in the ranges of tissue types at different spatial locations (e.g. top versus bottom of the brain). These inhomogeneities/non-uniformities can cause problems with algorithms based on histograms, quantiles, or raw intensities . Therefore, correction for image inhomogeneities is a crucial step in many analyses. The Freesurfer function performs the non-uniformity correction by , and the function of the same name will run the correction and return an image.
The Freesurfer function requires a MINC format (). For this to work, you can convert the object to a MINC file using :
[1] "/var/folders/1s/wrtqcpxn685_zk570bnx9_rr0000gr/T//RtmpSsiADQ/filed7173a83a9e3.mnc"We can pass this MINC file into the function, which will run the correction and then convert the output MNC to a NIfTI object.
[1] "nifti"
attr(,"package")
[1] "oro.nifti"We see that the results are indeed objects. We can plot the estimated bias field (log-transformed for display purposes) side-by-side with the image to view which areas had been differentially corrected (Figure \(\ref{fig:nu_correct_plot}\)).
In addition to the and functions above, we have a wrapper function
for reading in MINC files, after conversion to NIfTI files. If you pass
in a object in directly into nu_correct, the function will
automatically convert any NIfTI input files, and then run the correction
(shown below). We can also pass in a mask of the brain (see next
section) to run the correction only the areas of the brain.
Overall, this correction is a way to make the intensities of the brain more homogeneous spatially. This method is different from that implemented in FSL (and therefore ), so it provides an alternative method to the R user than currently available.
Brain extraction: the function
The function will segment the brain from the remainder of the image, such as extra-cranial tissues. Other imaging software in R have implemented the watershed algorithm, such as . These methods have not been directly adapted for MRI nor specifically for brain extraction. In , we can pass in the object and the output is a brain-extracted object.
The result is in figure \(\ref{fig:watershed_plot}\), where we see areas of the skull, eyes, face, and other areas of the image are removed. We do see some areas remain that may be part of some of the membranes between the brain and the skull, but this looks like an adequate brain extraction for most analyses.
As the result is a object, we can create a mask by standard logical operations for arrays. As MRI scans are typically positive-valued, the positive areas of the image are the “brain”:
We can then use this mask to perform operations on the image, such as subsetting.
Segmentations of brain structures
Freesurfer is commonly used to segment cortical and subcortical structures of the brain. We can visualize images of these segmentations, which are located in the “mri” folder. We will choose the colors based on the Freesurfer look up table (LUT), which values can be explored at https://surfer.nmr.mgh.harvard.edu/fswiki/FsTutorial/AnatomicalROI/FreeSurferColorLUT. This look up table provides a label for each structure and the color associated with it:
index label R G B A
1 0 Unknown 0 0 0 0
2 1 Left-Cerebral-Exterior 70 130 180 0
3 2 Left-Cerebral-White-Matter 245 245 245 0This object is included in and denotes the indices, labels, and color representation of the structure. We note that the alpha channel is set to \(0\) for all regions of interest, so we will not use it in the calculation of the colors from RGB space. This LUT allows visualizations produced in R to be consistent with those from Freesurfer.
if (have_fs()) {
seg_file = file = file.path(fs_subj_dir(), "bert", "mri", "aseg.mgz")
seg = readmgz(seg_file)
breaks = c(-1, fs_lut$index)
colors = rgb(fs_lut$R, fs_lut$G, fs_lut$B, alpha = 255/2, maxColorValue = 255)
ortho2(ss, seg, col.y = colors, ybreaks = breaks)
}
Note above that the number of breaks must be one larger than the number of colors and the indices start at zero, so we add an additional element to the indices. The result in figure \(\ref{fig:seg_file}\) shows the image with colors overlaid.
Reading in anatomical statistics for brain structures
We have explored the spatial results in the brain images, but not the quantitative information about the brain and sub-structures that are available from Freesurfer output. The “aseg.stats” in the “stats” folder for subject bert corresponds to measures and statistics from the anatomical segmentation. The function reads this corresponding file and creates a list of 2 different s:
if (have_fs()) {
file = file.path(fs_subj_dir(), "bert", "stats", "aseg.stats")
out = read_aseg_stats(file)
names(out)
}[1] "measures" "structures"The element corresponds to global measurements of the brain (e.g. volume of the brain) as well as measures of gross anatomical structures (e.g. gray matter).
meaning value units
1 brain segmentation volume 1218061.000000 mm^3
2 brain segmentation volume without ventricles 1199378.000000 mm^3
3 volume of ventricles and choroid plexus 15041.000000 mm^3In some imaging analyses, comparing at these large measures of brain volume over time or across groups are of interest. Alternatively, the element corresponds to a set of measures and statistics for a set of fixed anatomical structures.
Index SegId NVoxels Volume_mm3 StructName normMean
1 1 4 6943 7316.9 Left-Lateral-Ventricle 36.0349
2 2 5 187 237.7 Left-Inf-Lat-Vent 51.8984
3 3 7 15043 15559.6 Left-Cerebellum-White-Matter 86.5151
normStdDev normMin normMax normRange
1 12.2361 14 89 75
2 9.4303 21 78 57
3 5.0521 38 99 61Similarly with global measures, these structure-specific measures are used in analysis. For example, measuring differences in hippocampus volumes across patients with Alzheimer’s disease and those without. Moreover, a large deviation in volume, globally or locally, for a specific subject may indicate atrophy of a structure or an indication of a segmentation error.
Converting surfaces using
Freesurfer includes segmentations of different surfaces of the brain
alongside the volumetric segmentations above. As the function provides a
tool to convert image volumes to a series of output formats, the (note
the “s”) allows users to convert between image surface formats. These
surfaces usually store sets of vertices and faces to be plotted in 3
dimensions. has implemented (with a function of the same name) as well
as functions to convert surfaces from Freesurfer to a set of triangles
in R, such as surface_to_triangles. We will read in the
left and right side of the pial surface of the brain and display the
surface using .
if (have_fs()) {
right_file = file.path(fs_subj_dir(),
"bert", "surf", "rh.pial")
right_triangles = surface_to_triangles(infile = right_file)
left_file = file.path(fs_subj_dir(),
"bert", "surf", "lh.pial")
left_triangles = surface_to_triangles(infile = left_file)
rgl::rgl.open()
rgl::rgl.triangles(right_triangles,
color = rainbow(nrow(right_triangles)))
rgl::rgl.triangles(left_triangles,
color = rainbow(nrow(left_triangles)))
}Thus, we can read in the output images, surfaces, and the tables of output metrics from Freesurfer.
Additional features
For the initial release, we did not implement a method to read the annotation files and other surface-based files that Freesurfer uses. Reading in these files are planned for a future release and may work with the functions described above. Freesurfer can also analyze diffusion tensor imaging data and some of the functions have been adapted for but have not been thoroughly tested.