---
title: "Analyzing and Visualizing Freesurfer Outputs"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Analyzing and Visualizing Freesurfer Outputs}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
bibliography: ref.bib
---



Welcome to this tutorial on analyzing and visualizing Freesurfer outputs using the `freesurfer` R package.
This guide will walk you through accessing and interpreting the rich data generated by Freesurfer, including anatomical segmentations, statistical measurements, and brain surfaces.

# Exploring Freesurfer Output: A Practical Guide

After Freesurfer completes its `recon-all` processing, as introduced in our "Setting up `freesurfer`" vignette, it organizes all the generated output into a neat, structured directory for each subject.
For this tutorial, we will focus on a subject named "bert," which is typically included with Freesurfer installations, making it ideal for reproducible examples.

If you were to run a full analysis using `recon_all`, your output would look something like this:


``` r
library(freesurfer)
```

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

Let's start by inspecting the main sub-directories within the "bert" subject's output folder.
This will give us an overview of where different types of data are stored.


``` r
list.files(path = file.path(fs_subj_dir(), "bert"))
```

```
## [1] "label"   "mri"     "scripts" "stats"  
## [5] "surf"    "tmp"     "touch"   "trash"
```

In this tutorial, we will specifically explore data found in three key directories:

* "mri": This folder contains various imaging data, such as T1-weighted images and anatomical segmentations.
* "stats": Here, you'll find statistical measures derived from different brain structures.
* "surf": This directory stores surface and curvature outputs generated by Freesurfer's processing.

## Visualizing Segmented Brain Structures

One of Freesurfer's most powerful features is its ability to segment cortical and subcortical brain structures.
You can visualize these segmentations, which are typically found in the "mri" folder (e.g., `aseg.mgz`).

To ensure your visualizations are consistent with Freesurfer's own output, we'll use the Freesurfer look-up table (LUT).
This LUT is a crucial tool that maps numerical indices to specific anatomical labels and their corresponding RGB color values.
You can explore the full LUT on the [fswiki](https://surfer.nmr.mgh.harvard.edu/fswiki/FsTutorial/AnatomicalROI/FreeSurferColorLUT).
Conveniently, the `freesurfer` package includes this `fs_lut` object directly for easy access:


``` r
head(fs_lut, 3)
```

```
##   index                      label   R   G
## 1     0                    Unknown   0   0
## 2     1     Left-Cerebral-Exterior  70 130
## 3     2 Left-Cerebral-White-Matter 245 245
##     B A
## 1   0 0
## 2 180 0
## 3 245 0
```

As you can see, the `fs_lut` object provides indices, anatomical labels, and color representations for a wide range of brain structures.
Note that the alpha channel in `fs_lut` is set to 0 for all regions.
We will manually adjust the transparency when converting the RGB values to colors, ensuring that your R-generated visualizations match Freesurfer's native output.

Now, let's load the `aseg.mgz` file, which contains the anatomical segmentation, and overlay it onto the brain-extracted image.
As covered in the [Image processing vignette](image_processing.html), you will likely want to reorient images one read in.
Since the `neurobase::orient_rpi` returns a list with an image and the orientation of the image, we suggest making a custom function to reorient on need.


``` r
#' Helper function to reorient images to RPI
#' @param ... Arguments passed to neurobase::orient_rpi
#' @return Reoriented image
reorient <- function(...) {
  neurobase::orient_rpi(...)$img
}

# Read in t1 for background image
# and reorient it
t1_file = file.path(fs_subj_dir(), "bert", "mri", "T1.mgz")
bert_t1 = read_mgz(t1_file) |>
  reorient()

seg_file = file.path(fs_subj_dir(), "bert", "mri", "aseg.mgz")
seg = read_mgz(seg_file) |>
  reorient()

breaks = c(-1, fs_lut$index)

colors = rgb(
  fs_lut$R,
  fs_lut$G,
  fs_lut$B,
  alpha = 255 / 2, # Apply transparency
  maxColorValue = 255
)

neurobase::ortho2(
  bert_t1,
  seg,
  col.y = colors,
  ybreaks = breaks,
  mask = bert_t1
)
```

![Overlay of anatomical segmentation from Freesurfer recon-all command. Different colors represent distinct brain structures according to the Freesurfer LUT.](figure/seg_file-1.png)

As shown in the figure above, the overlay clearly displays various segmented brain structures, each distinguished by a unique color from the Freesurfer LUT.
A small tip: remember that the number of `breaks` must always be one greater than the number of `colors` you provide.
Also, the indices for `fs_lut` start at zero, so we add an extra element to the `breaks` to account for this.

## Analyzing Anatomical Statistics for Brain Structures

Beyond just visualizing segmented images, Freesurfer also provides a wealth of quantitative information about whole-brain and sub-structures in its statistical files.

The "aseg.stats" file, located in the "stats" folder for a subject like "bert," contains various measures and statistics derived from the anatomical segmentation.
The `freesurfer` package offers the `read_aseg_stats` function to easily parse this file.
This function returns a list containing two useful `data.frame`s:


``` r
file = file.path(fs_subj_dir(), "bert", "stats", "aseg.stats")
out = read_aseg_stats(file)
names(out)
```

```
## [1] "measures"   "structures"
```

The `measures` element within the output list provides global measurements of the brain, such as total brain volume, as well as measures for broad anatomical categories like gray matter volume.


``` r
head(out$measures[, c("meaning", "value", "units")], n = 3)
```

```
##                                        meaning
## 1                    brain segmentation volume
## 2 brain segmentation volume without ventricles
## 3      volume of ventricles and choroid plexus
##            value units
## 1 1218061.000000  mm^3
## 2 1199378.000000  mm^3
## 3   15041.000000  mm^3
```

These large-scale brain volume measures are frequently used in neuroimaging analyses, especially for comparing brain characteristics across different groups or tracking changes over time.
On the other hand, the `structures` element offers a more detailed set of measures and statistics for a predefined set of specific anatomical structures:


``` r
head(out$structures, n = 3)
```

```
##   Index SegId NVoxels Volume_mm3
## 1     1     4    6943     7316.9
## 2     2     5     187      237.7
## 3     3     7   15043    15559.6
##                     StructName normMean
## 1       Left-Lateral-Ventricle  36.0349
## 2            Left-Inf-Lat-Vent  51.8984
## 3 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        61
```

Similar to the global measures, these structure-specific metrics are crucial for conducting more detailed analyses.
For example, researchers might use these to investigate differences in hippocampus volumes between patients with Alzheimer's disease and healthy controls.
Furthermore, a significant deviation in volume (either globally or for a specific structure) in an individual subject could indicate atrophy or a potential segmentation error that warrants further inspection.

## Converting and Visualizing Brain Surfaces with `mris_convert`

Freesurfer excels not only at generating volumetric segmentations but also at providing detailed segmentations of different brain surfaces.
While `mri_convert` is your go-to for handling image volumes, `mris_convert` (notice the "s"!) is the specialized tool for converting between various image surface formats.
These surface files typically store sets of vertices and faces, which are essential for generating compelling 3D plots of the brain.

The `freesurfer` package includes an implementation of `mris_convert` (with a function of the same name) and also provides helpful functions like `surface_to_triangles`.
This latter function is particularly useful for converting Freesurfer surfaces into a format that's perfectly suited for 3D plotting in R.
Let's demonstrate this by reading in the left and right pial surfaces of the brain and visualizing them using the powerful `rgl` package [@rgl].




``` r
rgl.useNULL()
library(rgl)
right_file = file.path(
  fs_subj_dir(),
  "bert",
  "surf",
  "rh.pial"
)
right_triangles = surface_to_triangles(right_file)

left_file = file.path(
  fs_subj_dir(),
  "bert",
  "surf",
  "lh.pial"
)
left_triangles = surface_to_triangles(left_file)
```


``` r
# Open an rgl window
open3d()

triangles3d(
  right_triangles,
  color = rainbow(nrow(right_triangles))
)
triangles3d(
  left_triangles,
  color = rainbow(nrow(left_triangles))
)

# Adjust viewpoint for better visualization
view3d(0, -70)
```


```
## [1] FALSE
```

```
## glX 
##   3
```

![3D rendering of the pial surfaces of the left and right hemispheres from Freesurfer output, displayed using rgl.](rgl.png)


The 3D plot beautifully illustrates how the `freesurfer` package empowers you to load and visualize complex surface data directly within R.
By leveraging `rgl`'s interactive capabilities, you gain a powerful way to explore and analyze Freesurfer's diverse outputs, which include not only images and surfaces but also tabular metrics.

## Working with Label Files

While extensive testing of Freesurfer label file reading has not been performed, the `freesurfer` package does provide the `read_fs_label` function for this purpose.
Below is a practical example of how to read a label file for the left hemisphere cortex:


``` r
file = file.path(fs_subj_dir(), "bert", "label", "lh.cortex.label")
out = read_fs_label(file)
head(out)
```

```
##   vertex_num r_coord  a_coord s_coord
## 1          0 -11.127 -101.103  -1.062
## 2          1 -10.776 -100.985  -1.818
## 3          2 -10.967 -101.124  -1.850
## 4          3 -11.642 -101.242  -1.580
## 5          4 -12.188 -101.283  -2.012
## 6          5 -10.848 -101.056  -2.769
##          value
## 1 0.0000000000
## 2 0.0000000000
## 3 0.0000000000
## 4 0.0000000000
## 5 0.0000000000
## 6 0.0000000000
```

In these label files, the coordinates are primarily used, rather than the assigned value (which is typically zero in this specific example).
These files can also be quite useful for various registration purposes.
While a reader for these files is provided, their specific applications in research have not been thoroughly explored or tested for further discussion in this vignette.

