All code for this document is located at here.

In this tutorial we will discuss within a visit registration (co-registration) of multi-sequence MRI images.

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.smri" %in% packages) {
  devtools::install_github("muschellij2/kirby21.smri")
}
if (!"EveTemplate" %in% packages) {
  devtools::install_github("jfortin1/EveTemplate")
}

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.smri)
library(kirby21.base)
run_mods = c("T1", "T2", "FLAIR")
fnames = get_image_filenames_list_by_visit(
  ids = 113, 
  modalities = run_mods, 
  visits = c(1,2 ))
visit_1 = fnames$`1`$`113`
visit_2 = fnames$`2`$`113`
mods = visit_1 %>% nii.stub(bn = TRUE) %>% 
  strsplit("-") %>% 
  sapply(dplyr::last)
names(visit_1) = names(visit_2) = mods
visit_1 = visit_1[run_mods]
visit_2 = visit_2[run_mods]

Processing images within a visit

The function preprocess_mri_within from extrantsr wraps a series of steps. The function below will perform:

  1. N4 inhomogeneity (Tustison et al. 2010) correction to each image.
  2. Estimate the transformation to the first file (T1 image).
  3. Perform this transformation, registering the images, and interpolating the images using a Lanczos windowed sinc interpolator.

preprocess_mri_within can also perform skull stripping using BET, but we have shown in the brain extraction tutorial that running BET without running neck removal.

library(extrantsr)
outfiles = nii.stub(visit_1, bn = TRUE)
proc_files = paste0(outfiles, "_proc.nii.gz")
names(proc_files) = names(outfiles)
if (!all(file.exists(proc_files))) {
  extrantsr::preprocess_mri_within(
    files = visit_1,
    outfiles = proc_files,
    correct = TRUE,
    retimg = FALSE,
    correction = "N4")
}
proc_imgs = lapply(proc_files, readnii)
lapply(proc_imgs, ortho2)

$T1
NULL

$T2
NULL

$FLAIR
NULL

Brain Extraction of the T1 image

Similar to the brain extraction tutorial, we can run fslbet_robust to get a good brain mask for the T1 image.

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

Applying the brain mask

As each processed image is now registered to the T1 image, we can apply this mask to all image sequences. Here we will mask the images using mask_img and then we will drop empty dimensions based on the mask.

mask = ss > 0
proc_imgs = lapply(proc_imgs, mask_img, mask = mask)
dd = dropEmptyImageDimensions(mask, other.imgs = proc_imgs)
mask = dd$outimg
proc_imgs = dd$other.imgs

Again we will plot the masked images:

lapply(proc_imgs, ortho2)

$T1
NULL

$T2
NULL

$FLAIR
NULL

Inhomogeneity correction on skull-stripped images

Some believe that inhomogeneity correction should be done after skull stripping (even if done before as well). Here we will run extrantsr::bias_correct which calls ANTsR. We will pass in the mask:

n4_proc_imgs = plyr::llply(
  proc_imgs, 
  bias_correct, 
  correction = "N4", 
  mask = mask,
  retimg = TRUE,
  .progress = "text")

  |                                                                       
  |                                                                 |   0%
  |                                                                       
  |======================                                           |  33%
  |                                                                       
  |===========================================                      |  67%
  |                                                                       
  |=================================================================| 100%

Here we will write out the processed images for later use (in other tutorials):

outfiles = nii.stub(visit_1, bn = TRUE)
outfiles = paste0(outfiles, "_proc_N4_SS.nii.gz")
if (!all(file.exists(outfiles))) {
  mapply(function(img, outfile) {
    writenii(img, filename = outfile)
  }, n4_proc_imgs, outfiles)
}

Intensity Normalization

Here we will do intensity normalization of the brain. We will do z-score normalization, where the estimates of the mean and standard deviation are based on all voxels within the mask. This is referred to as whole-brain normalization.

norm_imgs = plyr::llply(
  n4_proc_imgs, 
  zscore_img,
  margin = NULL,
  centrality = "mean",
  variability = "sd",
  mask = mask,
  .progress = "text")

  |                                                                       
  |                                                                 |   0%
  |                                                                       
  |======================                                           |  33%
  |                                                                       
  |===========================================                      |  67%
  |                                                                       
  |=================================================================| 100%

Visualizing the marginal intensities

Here we can make a data.frame of the normalized images. We also create a long data.frame (called long) for plotting in ggplot2.

df = sapply(norm_imgs, function(x){
  x[ mask == 1 ]
})
long = reshape2::melt(df)
colnames(long) = c("ind", "sequence", "value")
long$ind = NULL
df = data.frame(df)

Marginal distributions

Here we plot the distributions of each imaging sequence separately.

ggplot(long, aes(x = value, colour = factor(sequence))) + 
  geom_line(stat = "density")

Bi-variate distributions

Here we make binned hexagrams to represent the 2-dimensional distributions of each imaging sequence against the other.

g = ggplot(df) + stat_binhex()
g + aes(x = T1, y = T2)

g + aes(x = T1, y = FLAIR)

g + aes(x = T2, y = FLAIR)