```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, comment = "", fig.height = 5.5, fig.width = 5.5, cache = TRUE)
```

## Overall Pipeline

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


## Intensity normalization
- Conventional MRI intensites (T1-w, T2-w, PD, FLAIR) are acquired in arbitrary units
- Images are not comparable across scanners, subjects, and visits, even when the same protocol is used. 
    - This affects algorithm performance, prediction, inference.
    - Even simple things like thresholding an image
- Intensity normalization brings the intensities to a common scale across people.
- In this tutorial we will normalize intensities within subject using two methods:
    - Whole-brain normalization
    - White Stripe normalization [@shinohara2014statistical].


```{r t1, echo=FALSE, warning=FALSE, message=FALSE}
library(ms.lesion)
library(neurobase)
library(WhiteStripe)
fnames = get_image_filenames_list_by_subject(group = "training", 
  type = "coregistered")
t1s = lapply(fnames, function(x) readnii(x["T1"]))
tissues = lapply(fnames, function(x) readnii(x["Tissue_Classes"]))
masks = lapply(fnames, function(x) readnii(x["Brain_Mask"]))

vals = mapply(function(t1, mask){
  mask_vals(t1, mask)
}, t1s, masks, SIMPLIFY = FALSE)

dens = lapply(vals, density)

plot_densities = function(dens, xlab = "Raw Intensities", 
                          main = "Whole Brain") {
  range_x = sapply(dens, function(d) range(d$x))
  range_x = range(range_x)
  range_y = sapply(dens, function(d) range(d$y))
  range_y = range(range_y)
  plot(dens[[1]], xlim = range_x, ylim = range_y, 
       xlab = xlab, main = main)
  for (idens in 2:length(dens)) {
    lines(dens[[idens]], col = idens)
  }
}
plot_boxplots = function(vals, 
                          main = "Whole Brain") {
  boxplots <- lapply(vals, boxplot, outline = FALSE, plot = FALSE)
  boxplots = lapply(boxplots, function(x) x$stats)
  boxplots <- do.call(cbind, boxplots)
  boxplot(boxplots, main = main)
}
```

## Visualizing whole-brain intensities (each line is a person)
- We will work with the T1-w images from the training data. 
- Full brain densities are mixtures of the three tissue class distributions.

```{r t1viz, warning=FALSE, message=FALSE, echo = FALSE}
plot_densities(dens, main = "Distribution of all Voxels in Brain Mask")
```


<!-- ## Visualizing the intensities by tissue class -->

<!-- <div class="columns-2"> -->
<!-- ```{r t1viz1, warning=FALSE, message=FALSE, echo = FALSE} -->
<!-- csf_vals = mapply(function(t1, mask){ -->
<!--   mask_vals(t1, mask == 1) -->
<!-- }, t1s, tissues, SIMPLIFY = FALSE) -->
<!-- csf_dens = lapply(csf_vals, density) -->
<!-- plot_densities(csf_dens, main = "CSF") -->
<!-- ``` -->

<!-- ```{r t1viz2, warning=FALSE, message=FALSE, echo = FALSE} -->
<!-- gm_vals = mapply(function(t1, mask){ -->
<!--   mask_vals(t1, mask == 2) -->
<!-- }, t1s, tissues, SIMPLIFY = FALSE) -->
<!-- gm_dens = lapply(gm_vals, density) -->
<!-- plot_densities(gm_dens, main = "Gray Matter") -->
<!-- ``` -->
<!-- </div> -->

<!-- ## Visualizing the intensities by tissue class -->
<!-- - Notice the complete non-overlap between some subjects' white matter distributions. -->

<!-- <div class="columns-2"> -->

<!-- ```{r t1viz3, warning=FALSE, message=FALSE, echo = FALSE} -->
<!-- wm_vals = mapply(function(t1, mask){ -->
<!--   mask_vals(t1, mask == 3) -->
<!-- }, t1s, tissues, SIMPLIFY = FALSE) -->
<!-- wm_dens = lapply(wm_vals, density) -->
<!-- plot_densities(wm_dens, main = "White Matter") -->
<!-- ``` -->

<!-- ```{r wm_raw_box, echo = FALSE} -->
<!-- plot_boxplots(wm_vals, main = "White Matter") -->
<!-- ``` -->
<!-- </div> -->


## Visualizing the intensities by tissue class

```{r threeway, warning=FALSE, message=FALSE, echo = FALSE, fig.height = 4, fig.width = 10}
csf_vals = mapply(function(t1, mask){
  mask_vals(t1, mask == 1)
}, t1s, tissues, SIMPLIFY = FALSE)
csf_dens = lapply(csf_vals, density)

wm_vals = mapply(function(t1, mask){
  mask_vals(t1, mask == 3)
}, t1s, tissues, SIMPLIFY = FALSE)
wm_dens = lapply(wm_vals, density)

gm_vals = mapply(function(t1, mask){
  mask_vals(t1, mask == 2)
}, t1s, tissues, SIMPLIFY = FALSE)
gm_dens = lapply(gm_vals, density)

par(mfrow = c(1, 3))
plot_densities(csf_dens, main = "CSF")
plot_densities(gm_dens, main = "Gray Matter")
plot_densities(wm_dens, main = "White Matter")
par(mfrow = c(1, 1))
```

And these are all the same scanner/protocol!

## Whole-brain normalization

- Let's Z-score each voxel using mean $\mu_{WB}$ and standard deviation $\sigma_{WB}$ computed from all voxels in the brain mask.

$$ T1_{WB} = \frac{T1 - \mu_{WB}}{\sigma_{WB}}$$


- `zscore_img` is a function in `neurobase` that does this.
- It takes an image and a binary mask. The default is to use all voxels in the brain mask.

```{r wbViz, echo=FALSE, warning=FALSE, message=FALSE}
t1_norm = mapply(function(img, mask){
  zscore_img(img = img, mask = mask)
}, t1s, masks, SIMPLIFY = FALSE)
```

```{r wbViz_show, eval=FALSE, warning=FALSE, message=FALSE}
zscore_img(img = img, mask = mask)
```

## Whole-brain normalized intensities

<div class="columns-2">
```{r t1viz1a, warning=FALSE, message=FALSE, echo = FALSE}
csf_vals = mapply(function(t1, mask){
  mask_vals(t1, mask == 1)
}, t1s, tissues, SIMPLIFY = FALSE)
csf_dens = lapply(csf_vals, density)
plot_densities(csf_dens, main = "CSF Before")
```

```{r wbViz1, warning=FALSE, message=FALSE, echo = FALSE}
csf_norm_vals = mapply(function(t1, mask){
  mask_vals(t1, mask == 1)
}, t1_norm, tissues, SIMPLIFY = FALSE)
csf_norm_dens = lapply(csf_norm_vals, density)
plot_densities(csf_norm_dens, 
               xlab = "Whole-brain Normalized Intensities", 
               main = "CSF After")
```
</div>


## Whole-brain normalized intensities
- Gray matter distributions are more comparable.

<div class="columns-2">
```{r t1viz2a, warning=FALSE, message=FALSE, echo = FALSE}
gm_vals = mapply(function(t1, mask){
  mask_vals(t1, mask == 2)
}, t1s, tissues, SIMPLIFY = FALSE)
gm_dens = lapply(gm_vals, density)
plot_densities(gm_dens, main = "Gray Matter Before")
```

```{r wbViz2, warning=FALSE, message=FALSE, echo = FALSE}
gm_norm_vals = mapply(function(t1, mask){
  mask_vals(t1, mask == 2)
}, t1_norm, tissues, SIMPLIFY = FALSE)
gm_norm_dens = lapply(gm_norm_vals, density)
plot_densities(gm_norm_dens, 
               xlab = "Whole-brain Normalized Intensities", 
               main = "Gray Matter After")
```
</div>


## Whole-brain normalized intensities
- White matter distributions are more comparable.

<div class="columns-2">
```{r t1viz3a, warning=FALSE, message=FALSE, echo = FALSE}
wm_vals = mapply(function(t1, mask){
  mask_vals(t1, mask == 3)
}, t1s, tissues, SIMPLIFY = FALSE)
wm_dens = lapply(wm_vals, density)
plot_densities(wm_dens, main = "White Matter Before")
```

```{r wbViz3, warning=FALSE, message=FALSE, echo = FALSE}
wm_norm_vals = mapply(function(t1, mask){
  mask_vals(t1, mask == 3)
}, t1_norm, tissues, SIMPLIFY = FALSE)
wm_norm_dens = lapply(wm_norm_vals, density)
plot_densities(wm_norm_dens, 
               xlab = "Whole-brain Normalized Intensities", 
               main = "White Matter After")
```
</div>


## Other Normalizations: White Stripe

- Whole-brain normalization may be sensitive to outliers.
- Lesions in MS can have very high intensities, which lead to bad estimates of mean/variance
    - Other more robust transformations may be used, such as using the median to center, IQR to scale, etc.
- White Stripe [@shinohara2014statistical] is based on parameters obtained from a sample of normal appearing white matter (NAWM), which is robust to outliers. 
    - The idea is to make normal appearing white matter comparable across subjects and visits.


## White Stripe normalization

Procedure:

1. Find white matter area on histogram

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

## White Stripe normalization

Procedure:

1. Find white matter area on histogram

2. Estimate mean $\mu_{WS}$ and variance $\sigma_{WS}$ of voxel intensities in that area

3. Normalize with those means/variances:
$$ T1_{WS} = \frac{T1 - \mu_{WS}}{\sigma_{WS}}$$

## White Stripe normalization

- After normalization, NAWM will have a standard normal distribution and units will be in standard deviations of NAWM.
- Gray matter and CSF distributions may not be comparable after White Stripe.

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


## White Stripe normalization

Procedure:

1. Find white matter area on histogram

2. Estimate mean $\mu_{WS}$ and variance $\sigma_{WS}$ of voxel intensities in that area

3. Normalize with those means/variances:
$$ T1_{WS} = \frac{T1 - \mu_{WS}}{\sigma_{WS}}$$

- After normalization, NAWM will have a standard normal distribution and units will be in standard deviations of NAWM.
- Gray matter and CSF distributions may not be comparable after White Stripe.

## White Stripe normalization code

```{r ws_show, eval=FALSE, warning = FALSE, message = FALSE, results='hide'}
library(WhiteStripe)
ind = whitestripe(img = t1, type = "T1", stripped = TRUE)$whitestripe.ind
ws_t1 = whitestripe_norm(t1, indices = ind)
```

- The `whitestripe` function takes an image, image type (in our case T1), and a logical indicating whether the image has been skull stripped.
- The indicies of voxels in the NAWM used for estimating the normalization parameters are located in the list element `$whitestripe.ind`.
- The function `whitestripe_norm` takes an image and the indicies from a call to `whitestripe` and returns the White Stripe normalized image as a nifti.

```{r ws, echo=FALSE, warning = FALSE, message = FALSE, results='hide'}
ws_norm = function(t1) {
  ind = whitestripe(img = t1,
                    type = "T1", 
                    stripped = TRUE)$whitestripe.ind
  whitestripe_norm(t1, indices = ind)
}
t1_ws_norm = lapply(t1s, ws_norm)
```


## WhiteStripe normalized intensities

<div class="columns-2">
```{r wbViz1a, warning=FALSE, message=FALSE, echo = FALSE}
plot_densities(csf_norm_dens, 
               xlab = "Whole-brain Normalized Intensities", main = "Whole-brain: CSF")
```

```{r ws_viz_csf, warning=FALSE, message=FALSE, echo = FALSE}
csf_ws_vals = mapply(function(t1, mask){
  mask_vals(t1, mask == 1)
}, t1_ws_norm, tissues, SIMPLIFY = FALSE)
csf_ws_dens = lapply(csf_ws_vals, density)
plot_densities(csf_ws_dens, 
               xlab = "WhiteStripe Normalized Intensities", main = "White Stripe: CSF")
```
</div>


## WhiteStripe normalized intensities

<div class="columns-2">
```{r wbViz2a, warning=FALSE, message=FALSE, echo = FALSE}
plot_densities(gm_norm_dens, 
               xlab = "Whole-brain Normalized Intensities", main = "Whole-brain: Gray Matter")
```

```{r ws_viz_gm, warning=FALSE, message=FALSE, echo = FALSE}
gm_ws_vals = mapply(function(t1, mask){
  mask_vals(t1, mask == 2)
}, t1_ws_norm, tissues, SIMPLIFY = FALSE)
gm_ws_dens = lapply(gm_ws_vals, density)
plot_densities(gm_ws_dens, 
               xlab = "WhiteStripe Normalized Intensities", main = "White Stripe: Gray Matter")
```
</div>

## WhiteStripe normalized intensities

<div class="columns-2">
```{r wbViz3a, warning=FALSE, message=FALSE, echo = FALSE}
plot_densities(wm_norm_dens, 
               xlab = "Whole-brain Normalized Intensities", main = "Whole-brain: White Matter")
```

```{r ws_viz_wm, warning=FALSE, message=FALSE, echo = FALSE}
wm_ws_vals = mapply(function(t1, mask){
  mask_vals(t1, mask == 3)
}, t1_ws_norm, tissues, SIMPLIFY = FALSE)
wm_ws_dens = lapply(wm_ws_vals, density)
plot_densities(wm_ws_dens, 
               xlab = "WhiteStripe Normalized Intensities", main = "White Stripe: White Matter")
```
</div>


## Conclusions
- Intensity normalization is an important step in any image analysis with more than one subject or time point to ensure comparability across images.
- White Stripe normalization may work better and have better interpretation than whole-brain normalization for subsequent lesion segmentation algorithms and analysis.
- Other intensity normalization methods that make intensites comparable across subjects for all tissues exist.
    - RAVEL, which is an extension of WhiteStripe is one example [@fortin2016removing].
    - Located at https://github.com/Jfortin1/RAVEL
      - This was shown to have better comparability than histogram matching 


## Website

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



## References {.smaller}


<!-- ## Intensity normalization with RAVEL -->

<!-- - WhiteStripe is great, it normalized the WM intensities across samples -->
<!-- - Grey matter (GM) and cerebrospinal fluid (CSF) intensites are still not comparable -->
<!-- - RAVEL extends WhiteStripe to make intensites comparable across subjects for all tissues -->






