```{r setup, include=FALSE}
library(methods)
library(scales)
library(fslr)
library(extrantsr)
knitr::opts_chunk$set(
  message = FALSE, echo = TRUE, 
  comment = "", 
  warning = FALSE)
```

## Overall Pipeline

<img src="flow.png" style="width: 75%; display: block; margin: auto;"> 


## Image segmentation 

- We are often interested in subdividing or segmenting the brain into meaningful biological regions of interest (ROIs) for an analysis.
- Examples: tissue segmentation, segmentation of deep gray matter structures, segmentation of pathology (MS lesions, tumors, ...)
- We will perform 3-class tissue segmentation in R using `fslr` and `ANTsR`:
    * Cerebrospinal fluid (CSF)
    * Gray Matter (GM)
    * White Matter (WM)

## Loading Data

- Let's read in the training T1 and brain mask for training subject 02.

```{r loading}
library(ms.lesion)
library(neurobase)
all_files = get_image_filenames_list_by_subject(
  group = "training", 
  type = "coregistered")
files = all_files$training02 
t1 = readnii(files["T1"])
mask = readnii(files["Brain_Mask"])
```

```{r window, echo = FALSE}
# t1 = robust_window(t1, probs = c(0, 0.9999))
# t1 = window_img(t1, window = c(0, 300))
```

```{r reduce, echo = FALSE}
run_mask = t1 > 100
dd_orig = drop_empty_dim(run_mask, keep_ind = TRUE)
rt1 = robust_window(t1)
```

## Tissue Segmentation: Large Outliers
- Many tissue class segmentations are based on k-means clustering.
- These methods can be skewed by large outliers.

```{r hist_vals}
hist(t1, mask = mask, breaks = 2000); text(x = 600, y = 40000, '"outliers"')
```

## Where are the outliers?

```{r which_big, eval = FALSE}
ortho2(rt1, t1 > 400, col.y = alpha("red", 0.5)) # xyz - cog of a region
```

```{r which_big_show, echo = FALSE}
xrt1 = apply_empty_dim(rt1, inds = dd_orig$inds)
xt1 = apply_empty_dim(t1, inds = dd_orig$inds)
ortho2(xrt1, xt1 > 400, xyz = xyz(xt1 > 400), col.y = alpha("red", 0.5)) 
```

```{r run_window, echo=FALSE}
t1[ t1 < 0 ] = 0
t1 = mask_img(t1, mask)
rt1 = robust_window(t1)
```

## What does the histogram look like after truncating outliers?

```{r which_big_after_robust}
rt1 = robust_window(t1)
hist(rt1, mask = mask, breaks = 2000); 
```

## Tissue Segmentation using FSL FAST

- FAST is based on a hidden Markov random field model and an Expectation-Maximization algorithm [@fast].
- It jointly produces a bias field corrected image and a probabilistic tissue segmentation.
- More robust to noise and outliers than finite mixture model-based methods that do not incorporate spatial information.


The `fslr` function `fast` calls fast from FSL.  The `--nobias` option tells FSL to not perform inhomogeneity correction (N4 already performed in `ANTsR`). 

```{r fast_show, eval = FALSE}
t1fast = fast(t1, 
              outfile = paste0(nii.stub(t1file), "_FAST"), 
              opts = "--nobias")
```

```{r fast_show_run, echo = FALSE}
t1file = files["T1"]
outfile = paste0(nii.stub(t1file, bn = TRUE), "_FAST.nii.gz")
if (!file.exists(outfile)) {
  t1fast = fast(t1, 
                opts = "--nobias")
  writenii(t1fast, filename = outfile)
} else {
  t1fast = readnii(outfile)
}
```

## FAST Results

FAST assumes three tissue classes and produces an image with the three labels, ordered by increasing within-class mean intensities. In a T1 image, this results in:

* Level 1: CSF
* Level 2: Gray Matter
* Level 3: White Matter

## FAST: White Matter

```{r fast_wm_nonrobust, eval = FALSE}
ortho2(rt1, t1fast == 3, col.y = alpha("red", 0.5), text = "White Matter")
```

```{r fast_wm_nonrobust_show, echo = FALSE}
xrt1 = apply_empty_dim(rt1, inds = dd_orig$inds)
xt1fast = apply_empty_dim(t1fast, inds = dd_orig$inds)
ortho2(xrt1, xt1fast == 3, xyz=xyz(xt1 > 400), col.y = alpha("red", 0.5), text = "White Matter")
```

## FAST: Gray Matter

```{r fast_gm_nonrobust, eval = FALSE}
ortho2(rt1, t1fast == 2, col.y = alpha("red", 0.5), text = "Gray Matter")
```

```{r fast_gm_nonrobust_show, echo = FALSE}
ortho2(xrt1, xt1fast == 2, xyz=xyz(xt1 > 400), col.y = alpha("red", 0.5), text = "Gray Matter")
```

## FAST: CSF

```{r fast_csf_nonrobust, eval = FALSE}
ortho2(rt1, t1fast == 1, col.y = alpha("red", 0.5), text = "CSF")
```

```{r fast_csf_nonrobust_show, echo = FALSE}
ortho2(xrt1, xt1fast == 1, xyz=xyz(xt1 > 400), col.y = alpha("red", 0.5), text = "CSF")
```


## Removing large values: Is there an effect?

```{r fast_better_show, eval = FALSE}
robust_fast = fast(rt1, # the robust_window(t1)
              outfile = paste0(nii.stub(t1file), "_FAST"), 
              opts = "--nobias")
```

```{r fast_better_run, echo = FALSE}
robust_fast = readnii(files["FAST"])
```

```{r multi_overlay, eval = FALSE, echo = FALSE}
separate = separate_img(robust_fast, levels = 1:3)
names(separate) = c("CSF", "GM", "WM")
L = c(T1 = list(t1), separate)
L$T1 = robust_window(L$T1)
r = range(L$T1)
L$T1 = (L$T1 - r[1])/(r[2] - r[1])
dd = dropEmptyImageDimensions(t1fast > 0, other.imgs = L)
L = dd$other.imgs
multi_overlay(L, z = 58, text = names(L), text.x = 0.5, text.y = 1.4,
              text.cex = 2.5)
```


## FAST with Window: White Matter

```{r prep_wm, echo = FALSE}
rt1_list = list(rt1, rt1)
m_list = list((t1fast==3), (robust_fast==3))
```
```{r fast_wm, echo = FALSE}
multi_overlay(x=rt1_list, y=m_list, col.y=alpha("red", 0.5), text=c("Raw", "Robust"), text.x=c(.3, .7), text.y=c(.1, .1))
```

## FAST with Window: Gray Matter

```{r prep_gm, echo = FALSE}
m_list = list((t1fast==2), (robust_fast==2))
```
```{r fast_gm, echo = FALSE}
multi_overlay(x=rt1_list, y=m_list, col.y=alpha("red", 0.5), text=c("Raw", "Robust"), text.x=c(.3, .7), text.y=c(.1, .1))
```

## FAST with Window: CSF

```{r prep_csf, echo = FALSE}
m_list = list((t1fast==1), (robust_fast==1))
```
```{r fast_csf, echo = FALSE}
multi_overlay(x=rt1_list, y=m_list, col.y=alpha("red", 0.5), text=c("Raw", "Robust"), text.x=c(.3, .7), text.y=c(.1, .1))
```

## FAST Results

- Overall the results look good
    - Not much difference after dampening outliers using `robust_window`

## Tissue Segmentation using ANTsR, extrantsr

- Uses Atropos [@avants2011open] 
    - 3D K-means clustering + a Markov random field
- The `extrantsr::otropos` function works with `nifti` objects
    - calls `ANTsR::atropos` function

```{r otropos_show, eval = FALSE}
t1_otropos = otropos(a = t1, x = mask) # using original data
t1seg = t1_otropos$segmentation
```

```{r otropos_run, echo = FALSE}
t1file = files["T1"]
outfile = paste0(nii.stub(t1file, bn = TRUE), "_Atropos.nii.gz")
if (!file.exists(outfile)) {
  t1_otropos = otropos(a = t1, x = mask) # using robust
  t1seg = t1_otropos$segmentation
  writenii(t1seg, filename = outfile)
} else {
  t1seg = readnii(outfile)
}
```

## Atropos: White Matter

```{r otropos_wm, eval = FALSE}
ortho2(rt1, t1seg == 3, col.y = alpha("red", 0.5), text = "White Matter")
```

```{r otropos_wm_show, echo = FALSE}
xt1seg = apply_empty_dim(t1seg, inds = dd_orig$inds)
ortho2(xrt1, xt1seg == 3, xyz=xyz(xt1 > 400), col.y = alpha("red", 0.5), text = "White Matter")
```

## Atropos: Gray Matter

```{r otropos_gm, eval = FALSE}
ortho2(rt1, t1seg == 2, col.y = alpha("red", 0.5), text = "Gray Matter")
```

```{r otropos_gm_show, echo = FALSE}
ortho2(xrt1, xt1seg == 2, xyz=xyz(xt1 > 400), col.y = alpha("red", 0.5), text = "Gray Matter")
```

## Atropos: CSF


```{r otropos_csf, eval = FALSE}
ortho2(rt1, t1seg == 1, col.y = alpha("red", 0.5), text = "CSF")
```

```{r otropos_csf_show, echo = FALSE}
ortho2(xrt1, xt1seg == 1, xyz=xyz(xt1 > 400), col.y = alpha("red", 0.5), text = "CSF")
```

## Default Atropos Results

- Overall the results look good
	- The k-means clustering is not affected by outliers in this case
- We will still try using `robust_window`


## Atropos using Windowing


```{r robust_otropos_show, eval = FALSE}
robust_t1_otropos = otropos(a = rt1, x = mask) # using robust
robust_t1seg = robust_t1_otropos$segmentation
```

```{r robust_otropos_run, echo = FALSE}
robust_t1seg = readnii(files["Tissue_Classes"])
```

```{r atropos_robust_overall_show, eval = FALSE}
double_ortho(rt1, robust_t1seg)
```

```{r atropos_robust_overall, echo = FALSE}
xrobust_t1seg = apply_empty_dim(
  robust_t1seg, 
  inds = dd_orig$inds)
double_ortho(xrt1, xrobust_t1seg)
```


## Atropos with Window: White Matter

```{r otropos_wm_robust_show, eval = FALSE}
ortho2(rt1, robust_t1seg == 3, col.y = alpha("red", 0.5), text = "White Matter")
```

```{r otropos_wm_robust, echo = FALSE}
ortho2(xrt1, xrobust_t1seg == 3, xyz=xyz(xt1 > 400), col.y = alpha("red", 0.5), text = "White Matter")
```

## Atropos with Window: Gray Matter

```{r otropos_gm_robust_show, eval = FALSE}
ortho2(rt1, robust_t1seg == 2, col.y = alpha("red", 0.5), text = "Gray Matter")
```

```{r otropos_gm_robust, echo = FALSE}
ortho2(xrt1, xrobust_t1seg == 2, xyz=xyz(xt1 > 400), col.y = alpha("red", 0.5), text = "Gray Matter")
```

## Atropos with Window: CSF

```{r otropos_csf_robust_show, eval = FALSE}
ortho2(rt1, robust_t1seg == 1, col.y = alpha("red", 0.5), text = "CSF")
```

```{r otropos_csf_robust, echo = FALSE}
ortho2(xrt1, xrobust_t1seg == 1, xyz=xyz(xt1 > 400), col.y = alpha("red", 0.5), text = "CSF")
```


## Atropos with Window Results

- Overall the results look like they reasonably separate the classes
    - No ground truth
- Winsorizing large outliers did not make a difference for this image

```{r prep_wm_atropos, echo = FALSE}
m_list = list(t1seg == 3, robust_t1seg==3)
```

```{r prep_wm_atropos_show, echo = FALSE}
multi_overlay(x=rt1_list, y=m_list, col.y=alpha("red", 0.5), text=c("Raw", "Robust"), text.x=c(.3, .7), text.y=c(.1, .1))
```


## Atropos WM vs. FAST WM


```{r prep_wm_atropos_fast, echo = FALSE}
m_list = list(t1fast == 3, robust_t1seg==3)
```

```{r prep_wm_atropos_fast_show, echo = FALSE}
multi_overlay(x=rt1_list, y=m_list, col.y=alpha("red", 0.5), text=c("FAST", "Atropos"), text.x=c(.3, .7), text.y=c(.1, .1))
```

## Estimating the Volume of Each Class

We can create a table which will count the number of voxels in each category:

```{r tabs}
tab_fsl = table(robust_fast[ robust_fast != 0])
tab_ants = table(robust_t1seg[ robust_t1seg != 0])
names(tab_fsl) = names(tab_ants) = c("CSF", "GM", "WM")
tab_fsl
tab_ants
out_mask = robust_fast != 0 | robust_t1seg != 0 
table(robust_fast[ out_mask ], robust_t1seg[ out_mask ])
```

## Estimating the Volume of Each Class

By multiplying by the voxel resolution (in cubic centimeters) using the `voxres` function, we can get volumes

```{r volumes}
vres = oro.nifti::voxres(t1, units = "cm")
print(vres)
vol_fsl = tab_fsl * vres
vol_fsl
vol_ants = tab_ants * vres
vol_ants
```


## Website

<a href="../index.html" target="_blank">http://johnmuschelli.com/imaging_in_r</a>



## References {.smaller}
