R function setup

To use , a working installation of Freesurfer is required (downloads available: ). The following code was run using Freesurfer version “NA”. The Freesurfer version can be accessed using the function. The path of Freesurfer must also be set. When using R from a shell environment, after the environment variable is set (which is done when installing Freesurfer), will use this as the path to Freesurfer. If using R through a graphical user interface (GUI) such as RStudio (RStudio, Boston, MA), environmental variables and paths are not explicitly exported. Therefore, is not set and will try the default directories of Mac OSX and Linux. Freesurfer is only available on Windows via a virtual machine. If the user did not perform an standard installation of Freesurfer, the path to Freesurfer can be specified using . The function tests whether a user has a Freesurfer installation, returning a logical, which is useful for statements within examples. If function returns is , the function will return the directory of the Freesurfer installation.

Reconstruction pipeline in Freesurfer

The Freesurfer pipeline and analysis workflow for neuroanatomical images is designed to work with T1-weighted structural MRI of the brain. The full pipeline is implemented in the Freesurfer function, where the “recon” stands for reconstruction (). The recon-all function is the main workhorse of Freesurfer and is the most commonly used. Using the flag in the the function performs over 30 different steps and takes 20-40 hours to fully process a subject (https://surfer.nmr.mgh.harvard.edu/fswiki/recon-all). This process is the recommended way of fully processing an T1-weighted image in Freesurfer, and is implemented in the recon_all function.

In the recon_all function, users must specify the input file (a T1-weighted image), the output directory (if different than ), and the subject identifier. The results will be written in the individual subject directory, a sub-directory of . The syntax is:

recon_all(infile, outdir, subjid)

If there are problems with the result of this processing, there are multiple steps where users can edit certain parts of the processing, such as brain extraction, where non-brain tissues are removed from the image. The remainder of the pipeline can be run after these steps are corrected. The full pipeline is broken down into 3 separate sets of steps, referred to as , , and , which correspond to the same-named flags in used to initiate these steps. We have written wrapper functions , , and , respectively, so users can run pieces of the pipeline if desired or restart a failed process after correction to the data.

Imaging formats in and R

The package relies on the package implementation of images (referred to as objects) that are in the Neuroimaging Informatics Technology Initiative (NIfTI) format. For Freesurfer functions that require an image, the R functions that call those Freesurfer functions will take in a file name or a object. The R code will convert the to the corresponding input required for Freesurfer. From the user’s perspective, the input/output process is all within R, with one object format (). The advantage of this approach is that the user can read in an image, do manipulations of the object using standard syntax for arrays, and pass this object into the R function. Thus, users can use R functionality to manipulate objects while seamlessly passing these object to Freesurfer through .

Other Freesurfer functions require imaging formats other than NIfTI, such as the Medical Imaging NetCDF (MINC) format (). The Freesurfer installation provides functions to convert from MINC to NIfTI formats and these are implemented in functions such as and in R. Moreover, the Freesurfer function has been interfaced in (same function name), which allows for a more general conversion tool of imaging types for R users than currently implemented in native R. Thus, many formats can be converted to NIfTI and then read into R using the function from .

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):

recon_all(infile = "/path/to/T1.nii", subjid = "bert")

We see the result of this output in the “bert” directory, which includes a series of sub-directories:

if (have_fs()) {
  list.files(path  = file.path(fs_subj_dir(), "bert"))
}

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:

if (have_fs()) {
  img_mgz = readmgz(t1_mgz)
  identical(img, img_mgz)
}

Now that we have the image in R, we can plot it using the standard plotting tools for nifti objects:

if (have_fs()) {
  neurobase::ortho2(img, add.orient = FALSE, mask = img > 40)
}

The result is in figure , 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 ≥ 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:

if (have_fs()) {
  neurobase::ortho2(reoriented_img, mask = reoriented_img > 40)
}

The result is in figure , 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 :

if (have_fs()) {
  mnc = nii2mnc(reoriented_img)
  print(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.

if (have_fs()) {
  nu_from_mnc = nu_correct(file = mnc)
  class(nu_from_mnc)
}

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

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.

if (have_fs()) {
  nu_masked = nu_correct(file = reoriented_img, mask = mask)
}

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.

if (have_fs()) {
  ss = mri_watershed(img)
  ortho2(ss, mask = ss)
}

The result is in figure , 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”:

if (have_fs()) {
  mask = ss > 0
}

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:

if (have_fs()) {
  head(freesurfer::fs_lut, 3)
}

This 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 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)
}

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

if (have_fs()) {
  head(out$measures[, c("meaning", "value", "units")], n = 3)
}

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

if (have_fs()) {
  head(out$structures, n = 3)
}

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

Label files

Although we have not thoroughly tested reading in a label file from Freesurfer, we have provided the function. Here we will read a label file for the left hemisphere cortex:

if (have_fs()) {
  file = file.path(fs_subj_dir(), "bert", "label", "lh.cortex.label")
  out = read_fs_label(file)
  head(out)
}

The coordinates are mostly used in these files, not the value assigned (in this case all zeroes). These files can be used for registration as well. Overall, we have provided a reader for users, but have not used them in our research and therefore have not tested them extensively to discuss them further.