Tissue Class Segmentation
John Muschelli
2021-02-16
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.
1 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")
}2 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]3 Using information from the T1 image
3.1 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.
3.2 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])
}
}3.2.1 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.
3.2.2 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")
rm(list = "df")3.3 Tissue-Class Segmentation with FAST
library(fslr)
outfile = nii.stub(t1_fname, bn = TRUE)
outfile = paste0(outfile, "_FAST")
prob_files = paste0(outfile,
"_pve_", 0:2,
".nii.gz")
seg_outfile = paste0(outfile, "_Seg.nii.gz")
if (!all(file.exists(
c(seg_outfile, prob_files)
))) {
fast_hard_seg = fast(file = ss_red,
outfile = outfile,
out_type = "seg",
opts = "--nobias")
writenii(fast_hard_seg, seg_outfile)
} else {
fast_hard_seg = readnii(seg_outfile)
}
fast_seg = vector(mode = "list", length = 3)
for (i in 1:3) {
fast_seg[[i]] = readnii(prob_files[i])
} 3.3.1 FAST results
Let’s see the results of the FAST segmentation:
double_ortho(ss_red, hard_seg)
3.3.2 FAST intensity histograms
Again, we can look at the distribution of values, ad now we can compare distributions of the values from FAST to that of atropos.
df = data.frame(value = ss_red[ss_red > 0],
fast = fast_hard_seg[ss_red > 0],
atropos = hard_seg[ss_red > 0],
ind = which(ss_red > 0)
)
df = reshape2::melt(df, id.vars = c("ind", "value"),
measure.vars = c("fast", "atropos"),
value.name = "class",
variable.name = "segmentation")
df = df %>% arrange(ind)
df$ind = NULL
df$class = c("CSF", "GM", "WM")[df$class]
ggplot(df, aes(x = value, colour = factor(class))) +
geom_line(stat = "density") +
xlim(c(0, 1e6)) +
facet_wrap(~ segmentation, ncol = 1)Warning: Removed 596 rows containing non-finite values (stat_density).
rm(list = "df")3.4 Tissue-Class Segmentation with SPM
In the brain extraction tutorial we discuss the SPM segmenation procedures and show how to use them to produce hard segmentations, probability maps, and a brain extracted image. We will use the results of that tutorial to compare to that of atropos. We will exclude any tissues outside of GM, WM, and CSF (those > 3).
outfile = nii.stub(t1_fname, bn = TRUE)
outfile = file.path("..", "brain_extraction", outfile)
outfile = paste0(outfile, "_SPM_Seg.nii.gz")
spm_hard_seg = readnii(outfile)
spm_hard_seg[ spm_hard_seg > 3] = 0
dd = dropEmptyImageDimensions(
ss,
other.imgs = spm_hard_seg)
spm_hard_seg_red = dd$other.imgs3.4.1 SPM results
double_ortho(ss_red, spm_hard_seg_red)
Remember however, in the SPM segmentation, 1 is GM, 2 is WM, 3 is CSF, and in Atropos/FAST, 1 is CSF, 2 is GM, 3 is WM, .
spm_recode = niftiarr(spm_hard_seg_red, 0)
spm_recode[ spm_hard_seg_red %in% 1 ] = 2
spm_recode[ spm_hard_seg_red %in% 2 ] = 3
spm_recode[ spm_hard_seg_red %in% 3 ] = 1double_ortho(ss_red, spm_recode)
df = data.frame(spm = spm_recode[spm_recode > 0 | hard_seg > 0],
atropos = hard_seg[spm_recode > 0 | hard_seg > 0],
value = ss_red[spm_recode > 0 | hard_seg > 0])
df$spm = c("Background", "CSF", "GM", "WM")[df$spm + 1]
df$atropos = c("Background", "CSF", "GM", "WM")[df$atropos + 1]
df$spm = factor(df$spm, levels = c("Background", "CSF", "GM", "WM"))
df$atropos = factor(df$atropos, levels = c("Background", "CSF", "GM", "WM"))
tab = with(df, table(spm, atropos))
print(tab) atropos
spm Background CSF GM WM
Background 0 27972 1053 756
CSF 50255 102999 7146 0
GM 23839 66622 450565 60141
WM 455 0 946 339557We can also compare the 2 segmentations. Here, if we assume the SPM segmentation as the “gold standard” and the Atropos one as another “prediction”, we can look at the differences. Anywhere they both agree (both are a 1) it will be deemed a true positive and will be in green. Anywhere the Atropos segmentation includes a voxel but the SPM segmentation did not, it will deemed a false positive and will be in blue, vice versa in red will be a false negative.
compare = spm_recode == hard_seg
compare[ (spm_recode > 0 | hard_seg > 0) & !compare ] = 2
compare[ spm_recode == 0 & hard_seg == 0 ] = 0ortho2(ss_red, compare, col.y = alpha(c("blue", "red"), 0.5))
double_ortho(ss_red, compare, col.y = alpha(c("blue", "red"), 0.5))
x = list(ss_red,
ss_red)
y = list(spm = spm_recode,
atropos = hard_seg)
z = floor(nsli(ss_red)/2)
multi_overlay(x, y, z = z, col.y = alpha(hotmetal(), 0.25))
3.4.2 SPM intensity histograms
Although there may be places in the brain where SPM calls a class CSF, WM, or GM where the brain mask is zero, we will exclude these in the comparison to fast and atropos for a common comparison. We will make sure that if voxels within the brain mask are labeled as zero in the SPM segmentation, we will denote these as Background.
df = data.frame(value = ss_red[ss_red > 0],
fast = fast_hard_seg[ss_red > 0],
atropos = hard_seg[ss_red > 0],
spm = spm_recode[ss_red > 0],
ind = which(ss_red > 0)
)
df = reshape2::melt(df, id.vars = c("ind", "value"),
measure.vars = c("fast", "atropos", "spm"),
value.name = "class",
variable.name = "segmentation")
df = df %>% arrange(ind)
df$ind = NULL
df$class = c("Background", "CSF", "GM", "WM")[df$class + 1]
ggplot(df, aes(x = value, colour = factor(class))) +
geom_line(stat = "density") +
xlim(c(0, 1e6)) +
facet_wrap(~ segmentation, ncol = 1)Warning: Removed 894 rows containing non-finite values (stat_density).
rm(list = "df")3.5 Discussion
Note, atropos and fast generally require a skull-stripped image. Many skull-stripping algorithms remove the extra-cortical areas of the brain but inside the skull, which generally are CSF spaces with meninges. These CSF spaces are dropped after skull-stripping/brain extraction. If we are trying to consistently measure the CSF or “whole brain volume” (if that includes CSF spaces), this may cause issues. The SPM segmentation usually includes more CSF spaces, but we have shown in the brain extraction tutorial that there are areas that BET denotes as brain and SPM does not on the surface.
4 Session Info
devtools::session_info()─ Session info ───────────────────────────────────────────────────────────────
setting value
version R version 4.0.2 (2020-06-22)
os macOS Catalina 10.15.7
system x86_64, darwin17.0
ui X11
language (EN)
collate en_US.UTF-8
ctype en_US.UTF-8
tz America/New_York
date 2021-02-16
─ Packages ───────────────────────────────────────────────────────────────────
package * version date lib
abind 1.4-5 2016-07-21 [2]
animation * 2.6 2018-12-11 [2]
ANTsR * 0.5.6.1 2020-06-01 [2]
ANTsRCore * 0.7.4.6 2020-07-07 [2]
assertthat 0.2.1 2019-03-21 [2]
base64enc 0.1-3 2015-07-28 [2]
bindr 0.1.1 2018-03-13 [2]
bindrcpp * 0.2.2 2018-03-29 [2]
bitops 1.0-6 2013-08-17 [2]
cachem 1.0.4 2021-02-13 [1]
callr 3.5.1 2020-10-13 [1]
cli 2.3.0 2021-01-31 [1]
codetools 0.2-18 2020-11-04 [1]
colorout * 1.2-2 2020-06-01 [2]
colorspace 2.0-0 2020-11-11 [1]
crayon 1.4.1 2021-02-08 [1]
DBI 1.1.1 2021-01-15 [1]
dcm2niir * 0.6.9.1 2020-06-01 [2]
desc 1.2.0 2020-06-01 [2]
devtools 2.3.2 2020-09-18 [1]
digest 0.6.27 2020-10-24 [1]
dplyr * 1.0.4 2021-02-02 [1]
ellipsis 0.3.1 2020-05-15 [2]
evaluate 0.14 2019-05-28 [2]
EveTemplate * 1.0.0 2020-06-01 [2]
extrantsr * 3.9.13.1 2020-09-03 [2]
farver 2.0.3 2020-01-16 [2]
fastmap 1.1.0 2021-01-25 [1]
fs 1.5.0 2020-07-31 [2]
fslr * 2.25.0 2021-02-16 [1]
generics 0.1.0 2020-10-31 [1]
ggplot2 * 3.3.3 2020-12-30 [1]
git2r * 0.28.0 2021-01-11 [1]
glue 1.4.2 2020-08-27 [1]
gtable 0.3.0 2019-03-25 [2]
hexbin * 1.28.2 2021-01-08 [1]
highr 0.8 2019-03-20 [2]
htmltools 0.5.1.1 2021-01-22 [1]
htmlwidgets 1.5.3 2020-12-10 [1]
httr 1.4.2 2020-07-20 [2]
ichseg * 0.17.5 2020-08-13 [2]
ITKR 0.5.3.3.0 2021-02-15 [1]
jsonlite 1.7.2 2020-12-09 [1]
kirby21.base * 1.7.4 2020-10-01 [1]
kirby21.fmri * 1.7.0 2018-08-13 [2]
kirby21.t1 * 1.7.3.2 2021-01-09 [1]
knitr 1.31 2021-01-27 [1]
labeling 0.4.2 2020-10-20 [1]
lattice 0.20-41 2020-04-02 [2]
lifecycle 1.0.0 2021-02-15 [1]
magic 1.5-9 2018-09-17 [2]
magrittr 2.0.1 2020-11-17 [1]
matlabr * 1.6.0 2020-07-01 [2]
Matrix 1.3-2 2021-01-06 [1]
matrixStats * 0.58.0 2021-01-29 [1]
memoise 2.0.0 2021-01-26 [1]
mgcv 1.8-33 2020-08-27 [1]
mmand 1.6.1 2020-03-03 [2]
MNITemplate * 1.0.0 2020-06-01 [2]
munsell 0.5.0 2018-06-12 [2]
neurobase * 1.31.0 2020-10-07 [1]
neurohcp * 0.9.0 2020-10-19 [1]
nlme 3.1-152 2021-02-04 [1]
oasis * 3.0.4 2018-02-21 [2]
oro.nifti * 0.11.0 2020-09-04 [2]
papayaWidget * 0.7.1 2020-09-14 [1]
pillar 1.4.7 2020-11-20 [1]
pkgbuild 1.2.0 2020-12-15 [1]
pkgconfig 2.0.3 2019-09-22 [2]
pkgload 1.1.0 2020-05-29 [2]
plyr 1.8.6 2020-03-03 [2]
prettyunits 1.1.1 2020-01-24 [2]
processx 3.4.5 2020-11-30 [1]
ps 1.5.0 2020-12-05 [1]
purrr 0.3.4 2020-04-17 [2]
R.matlab 3.6.2 2018-09-27 [2]
R.methodsS3 * 1.8.1 2020-08-26 [1]
R.oo * 1.24.0 2020-08-26 [1]
R.utils * 2.10.1 2020-08-26 [1]
R6 2.5.0 2020-10-28 [1]
randomForest 4.6-14 2018-03-25 [2]
RColorBrewer * 1.1-2 2014-12-07 [2]
Rcpp 1.0.6 2021-01-15 [1]
RcppEigen 0.3.3.9.1 2020-12-17 [1]
remotes 2.2.0 2020-07-21 [2]
reshape2 * 1.4.4 2020-04-09 [2]
rlang 0.4.10 2020-12-30 [1]
rmarkdown * 2.6 2020-12-14 [1]
RNifti * 1.3.0 2020-12-04 [1]
rprojroot 2.0.2 2020-11-15 [1]
rstudioapi 0.13 2020-11-12 [1]
rvest * 0.3.6 2020-07-25 [2]
scales * 1.1.1 2020-05-11 [2]
sessioninfo 1.1.1 2018-11-05 [2]
spm12r * 2.8.2 2021-01-11 [1]
stapler 0.7.2 2020-07-09 [2]
stringi 1.5.3 2020-09-09 [1]
stringr * 1.4.0 2019-02-10 [2]
TCIApathfinder * 1.0.7 2021-02-04 [1]
testthat 3.0.2 2021-02-14 [1]
tibble 3.0.6 2021-01-29 [1]
tidyselect 1.1.0 2020-05-11 [2]
usethis 2.0.1 2021-02-10 [1]
vctrs 0.3.6 2020-12-17 [1]
WhiteStripe 2.3.2 2019-10-01 [2]
withr 2.4.1 2021-01-26 [1]
xfun 0.21 2021-02-10 [1]
xml2 * 1.3.2 2020-04-23 [2]
yaml * 2.2.1 2020-02-01 [2]
zoo * 1.8-8 2020-05-02 [2]
source
CRAN (R 4.0.0)
CRAN (R 4.0.0)
Github (ANTsX/ANTsR@9c7c9b7)
Github (muschellij2/ANTsRCore@61c37a1)
CRAN (R 4.0.0)
CRAN (R 4.0.0)
CRAN (R 4.0.0)
CRAN (R 4.0.0)
CRAN (R 4.0.0)
CRAN (R 4.0.2)
CRAN (R 4.0.2)
CRAN (R 4.0.2)
CRAN (R 4.0.2)
Github (jalvesaq/colorout@726d681)
CRAN (R 4.0.2)
CRAN (R 4.0.2)
CRAN (R 4.0.2)
Github (muschellij2/dcm2niir@4515762)
Github (muschellij2/desc@b0c374f)
CRAN (R 4.0.2)
CRAN (R 4.0.2)
CRAN (R 4.0.2)
CRAN (R 4.0.0)
CRAN (R 4.0.0)
Github (muschellij2/EveTemplate@ed54115)
Github (muschellij2/extrantsr@00c75ad)
CRAN (R 4.0.0)
CRAN (R 4.0.2)
CRAN (R 4.0.2)
local
CRAN (R 4.0.2)
CRAN (R 4.0.2)
Github (ropensci/git2r@4e342ca)
CRAN (R 4.0.2)
CRAN (R 4.0.0)
CRAN (R 4.0.2)
CRAN (R 4.0.0)
CRAN (R 4.0.2)
CRAN (R 4.0.2)
CRAN (R 4.0.2)
Github (muschellij2/ichseg@2741696)
Github (stnava/ITKR@ea0ac19)
CRAN (R 4.0.2)
local
CRAN (R 4.0.0)
local
CRAN (R 4.0.2)
CRAN (R 4.0.2)
CRAN (R 4.0.2)
CRAN (R 4.0.2)
CRAN (R 4.0.0)
CRAN (R 4.0.2)
local
CRAN (R 4.0.2)
CRAN (R 4.0.2)
CRAN (R 4.0.2)
CRAN (R 4.0.2)
CRAN (R 4.0.0)
Github (muschellij2/MNITemplate@46e8f81)
CRAN (R 4.0.0)
local
local
CRAN (R 4.0.2)
CRAN (R 4.0.0)
local
local
CRAN (R 4.0.2)
CRAN (R 4.0.2)
CRAN (R 4.0.0)
CRAN (R 4.0.0)
CRAN (R 4.0.0)
CRAN (R 4.0.0)
CRAN (R 4.0.2)
CRAN (R 4.0.2)
CRAN (R 4.0.0)
CRAN (R 4.0.0)
CRAN (R 4.0.2)
CRAN (R 4.0.2)
CRAN (R 4.0.2)
CRAN (R 4.0.2)
CRAN (R 4.0.0)
CRAN (R 4.0.0)
CRAN (R 4.0.2)
CRAN (R 4.0.2)
CRAN (R 4.0.2)
CRAN (R 4.0.0)
CRAN (R 4.0.2)
CRAN (R 4.0.2)
CRAN (R 4.0.2)
CRAN (R 4.0.2)
CRAN (R 4.0.2)
CRAN (R 4.0.2)
CRAN (R 4.0.0)
CRAN (R 4.0.0)
local
Github (muschellij2/stapler@79e23d2)
CRAN (R 4.0.2)
CRAN (R 4.0.0)
local
CRAN (R 4.0.2)
CRAN (R 4.0.2)
CRAN (R 4.0.0)
CRAN (R 4.0.2)
CRAN (R 4.0.2)
CRAN (R 4.0.0)
CRAN (R 4.0.2)
CRAN (R 4.0.2)
CRAN (R 4.0.0)
CRAN (R 4.0.0)
CRAN (R 4.0.0)
[1] /Users/johnmuschelli/Library/R/4.0/library
[2] /Library/Frameworks/R.framework/Versions/4.0/Resources/library