Tissue Class Segmentation

All code for this document is located at here.

In this tutorial we will discuss performing tissue class segmentation using atropos in ANTsR and it’s wrapper function in extrantsr, otropos.

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) {
  source("https://neuroconductor.org/neurocLite.R")
  neuroc_install("kirby21.base")    
}
if (!"kirby21.fmri" %in% packages) {
  source("https://neuroconductor.org/neurocLite.R")
  neuroc_install("kirby21.fmri")   
}

Loading Data

We will use the get_image_filenames_df function to extract the filenames on our hard disk for the T1 image and the fMRI images (4D).

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]

Using information from the T1 image

Brain extracted image

Please visit the brain extraction tutorial on how to extract a brain from this image. We will use the output from fslbet_robust from that tutorial.

outfile = nii.stub(t1_fname, bn = TRUE)
outfile = file.path("..", "brain_extraction", outfile)
outfile = paste0(outfile, "_SS.nii.gz")
ss = readnii(outfile)
ss_red = dropEmptyImageDimensions(ss)
ortho2(ss_red)

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

Tissue-Class Segmentation with Atropos

outfile = nii.stub(t1_fname, bn = TRUE)
prob_files = paste0(outfile,
                    "_prob_", 1:3,
                    ".nii.gz")
seg_outfile = paste0(outfile, "_Seg.nii.gz")

if (!all(file.exists(
  c(seg_outfile, prob_files)
  ))) {
  seg = extrantsr::otropos(
    ss_red, 
    x = ss_red > 0,
    v = 1)
  hard_seg = seg$segmentation
  writenii(hard_seg, seg_outfile)
  for (i in seq_along(seg$probabilityimages)) {
    writenii(seg$probabilityimages[[i]], prob_files[i]) 
  }
  # writenii(seg, )
} else {
  hard_seg = readnii(seg_outfile)
  seg = vector(mode = "list", length = 2)
  names(seg) = c("segmentation", "probabilityimages")
  seg$segmentation = hard_seg
  seg$probabilityimages = vector(mode = "list", length = 3)
  for (i in 1:3) {
    seg$probabilityimages[[i]] = readnii(prob_files[i]) 
  }  
}

Atropos results

Now we have a hard segmentation, which assigns a class with the maximum probability to that voxel. We also have a separate probability image for each tissue class.

double_ortho(ss_red, hard_seg)

We see that much of the structures have been segmented well, but there may be errors.

Atropos intensity histograms

We can also look at the distribution of intensities (marginally) for each tissue class. In atropos, the classes are ordered by mean intensity, so we can re-assign them to the corresponding tissue class

df = data.frame(value = ss_red[ss_red > 0],
                class = hard_seg[ss_red > 0])
df$class = c("CSF", "GM", "WM")[df$class]
ggplot(df, aes(x = value, colour = factor(class))) + geom_line(stat = "density")