All code for this document is located at here.

In this tutorial we will discuss performing brain segmentation using the brain extraction tool (BET) in fsl and a robust version using a wrapper function in extrantsr, fslbet_robust.

Data Packages

For this analysis, I will use one subject from the Kirby 21 data set. The kirby21.base and kirby21.fmri packages are necessary for this analysis and have the data we will be working on. You need devtools to install these. Please refer to installing devtools for additional instructions or troubleshooting.

packages = installed.packages()
packages = packages[, "Package"]
if (!"kirby21.base" %in% packages) {
  devtools::install_github("muschellij2/kirby21.base")
}
if (!"kirby21.t1" %in% packages) {
  devtools::install_github("muschellij2/kirby21.t1")
}

Loading Data

We will use the get_image_filenames_df function to extract the filenames on our hard disk for the T1 image.

library(kirby21.t1)
library(kirby21.base)
fnames = get_image_filenames_df(ids = 113, 
                    modalities = c("T1"), 
                    visits = c(1),
                    long = FALSE)
t1_fname = fnames$T1[1]

T1 image

Let’s take a look at the T1-weighted image.

t1 = readnii(t1_fname)
ortho2(t1)

rm(list = "t1")

Here we see the brain and other parts of the image are present. Most notably, the neck of the subject was imaged. Sometimes this can cause problems with segmentation and image registration.

Attempt 1: Brain Extraction of T1 image using BET

Here we will use FSL’s Brain Extraction Tool (BET) to extract the brain tissue from the rest of the image.

library(fslr)
outfile = nii.stub(t1_fname, bn = TRUE)
outfile = paste0(outfile, "_SS_Naive.nii.gz")
if (!file.exists(outfile)) {
  ss_naive = fslbet(infile = t1_fname, outfile = outfile)
} else {
  ss_naive = readnii(outfile)
}
ortho2(ss_naive)

We see that naively, BET does not perform well for this image.

Brain Extraction of T1 image using BET

Here we will use FSL’s Brain Extraction Tool (BET) to extract the brain tissue from the rest of the image. We use the modification of BET in extrantsr, which is called through fslbet_robust. In fslbet_robust, the image is corrected using the N4 inhomogeneity correction. The neck of the T1 image is then removed and then BET is run, the center of gravity (COG) is estimated, and BET is run with this new COG. We used a procedure where the neck is removed in 2 registration steps, which is more robust than just the one (which is the default).

outfile = nii.stub(t1_fname, bn = TRUE)
outfile = paste0(outfile, "_SS.nii.gz")
if (!file.exists(outfile)) {
  ss = extrantsr::fslbet_robust(t1_fname, 
    remover = "double_remove_neck",
    outfile = outfile)
} else {
  ss = readnii(outfile)
}

Let’s look at the skull-stripped image.

ortho2(ss)

Here we see the skull-stripped image. But did we drop out “brain areas”?

alpha = function(col, alpha = 1) {
  cols = t(col2rgb(col, alpha = FALSE)/255)
  rgb(cols, alpha = alpha)
}      
ortho2(t1_fname, ss > 0, col.y = alpha("red", 0.5))

We can again use dropEmptyImageDimensions to remove extraneous slices, which helps with reducing storage of the image on disk, zooming in while plotting, and may aid registration.

ss_red = dropEmptyImageDimensions(ss)
ortho2(ss_red)

Again, we can see the zoomed-in view of the image now.

Brain Extraction of T1 image using SPM

Note, to use SPM12, you must have MATLAB. We will use the spm12r package, which calls the matlabr package to call SPM functions within MATLAB.

library(spm12r)

The spm12_segment function takes in the original image and will perform segmentation on the entire image (background, skull, etc). The result of spm12_segment is a list of length 6 with a probability for each tissue type. The order of them is gray matter (GM), white matter (WM), and cerebrospinal fluid (CSF). We will create a brain mask from these first 3 classes. We can convert the probabilities to a hard segmentation using spm_probs_to_seg, which takes the maximum class probability to assign the class for each voxel. There are additional options for ties, but the default is to use the first class (GM > WM > CSF).

outfile = nii.stub(t1_fname, bn = TRUE)
spm_prob_files = paste0(outfile,
                        "_prob_", 1:6,
                        ".nii.gz")
ss_outfile = paste0(outfile, "_SPM_SS.nii.gz")
outfile = paste0(outfile, "_SPM_Seg.nii.gz")
outfiles = c(outfile, ss_outfile, spm_prob_files)
if (!all(file.exists(outfiles))) {
  spm_seg = spm12_segment(t1_fname)$outfiles
  spm_hard_seg = spm_probs_to_seg(img = spm_seg)
  writenii(spm_hard_seg, filename = outfile)
  
  spm_ss = spm_hard_seg >= 1 & spm_hard_seg <= 3
  writenii(spm_ss, filename = ss_outfile)
  
  for (i in seq_along(spm_seg)) {
    writenii(spm_seg[[i]], spm_prob_files[i]) 
  }  
} else {
  spm_seg = vector(mode = "list", 
                   length = length(spm_prob_files))
  for (i in seq_along(spm_seg)) {
    spm_seg[[i]] = readnii(spm_prob_files[i]) 
  }
  spm_hard_seg = readnii(outfile)
  spm_ss = readnii(ss_outfile)
}

Results of SPM Tissue Segmentation

double_ortho(t1_fname, spm_hard_seg)

Results of SPM Brain Segmentation

Here we will show again the T1 image with the corresponding skull stripped mask in red.

alpha = function(col, alpha = 1) {
  cols = t(col2rgb(col, alpha = FALSE)/255)
  rgb(cols, alpha = alpha)
}      
ortho2(t1_fname, spm_ss > 0, col.y = alpha("red", 0.5))

Additional Preprocessing to do

We can also remove the neck of the image and rerun the segmentation. We will run the double_remove_neck function to perform this.

outfile = nii.stub(t1_fname, bn = TRUE)
outfile = paste0(outfile,
                 "_noneck.nii.gz")
if (!file.exists(outfile)) {
  noneck = extrantsr::double_remove_neck(
    t1_fname,
    template.file = file.path(fslr::fsldir(), "data/standard",
                              "MNI152_T1_1mm_brain.nii.gz"), 
    template.mask = file.path(fslr::fsldir(),
                              "data/standard", 
                              "MNI152_T1_1mm_brain_mask.nii.gz"))
  writenii(noneck, filename = outfile)
} else {
  noneck = readnii(outfile)
}

Here we see that most of the neck was truly removed from the original image.

double_ortho(t1_fname, noneck)