```{r setup, include=FALSE, message = FALSE}
library(methods)
knitr::opts_chunk$set(echo = TRUE, comment = "")
library(methods)
library(ggplot2)
library(ms.lesion)
library(neurobase)
library(extrantsr)
library(scales)
```
## MS Lesion
Let's reset and read in the T1 image from a MS lesion data set:
```{r reading_in_image, eval = FALSE}
t1 = neurobase::readnii("training01_01_t1.nii.gz")
t1[ t1 < 0 ] = 0
```
```{r reading_in_image_run, echo = FALSE}
t1 = neurobase::readnii("../training01_01_t1.nii.gz")
t1[ t1 < 0 ] = 0
```
## Overall Pipeline
## Inhomogeneity correction
- Scans can have nonuniform intensities throughout the brain
- Usually low frequency - smooth over the brain (assumed)
- Referred to as bias, bias field, or inhomogeneity
Image From https://www.slicer.org/w/images/7/77/MRI_Bias_Field_Correction_Slicer3_close_up.png
## Image Data
It's hard to see subtler bias fields, but sometimes they can be seen visually.
```{r ortho2_show}
ortho2(robust_window(t1))
```
## Image Data
```{r ortho2_show_flair, eval = FALSE}
flair = neurobase::readnii("training01_01_flair.nii.gz")
ortho2(robust_window(flair))
```
```{r ortho2_run_flair, echo = FALSE}
flair = neurobase::readnii("../training01_01_flair.nii.gz")
flair[ flair < 0 ] = 0
flair = drop_empty_dim(flair > 50, other.imgs = flair)
flair = flair$other.imgs
ortho2(robust_window(flair))
```
## Image Data: Lightbox
```{r lightbox}
image(robust_window(t1), useRaster = TRUE)
```
## N4 Inhomogeneity Correction
We will use N4: Improved N3 Bias Correction [@tustison_n4itk_2010].
The model assumed in the N4 is:
$$
v(x) = u(x)f(x) + n(x)
$$
- $v$ is the given image
- $u$ is the uncorrupted image
- $f$ is the bias field
- $n$ is the noise (assumed to be independent and Gaussian)
- $x$ is a location in the image
## N4 Inhomogeneity Correction
The data is log-transformed and assuming a noise-free scenario, we have:
$$
\log(v(x)) = \log(u(x)) + \log( f(x) )
$$
- N4 uses a B-spline approximation of the bias field
- It iterates until a convergence criteria is met
- when the updated bias field is the same as the last iteration
- It outputs the data back in the original units (not log-transformed)
## Bias Field Correction
Here we will use the `bias_correct` function in `extrantsr`, which calls `n4BiasFieldCorrection` from `ANTsR`.
You can pass in the image:
```{r bc_show, message = FALSE, eval = FALSE}
library(extrantsr)
bc_t1 = bias_correct(file = t1, correction = "N4")
```
```{r bc_run, echo = FALSE}
out_fname = "../output/training01_01_t1_n4.nii.gz"
if (!file.exists(out_fname)) {
bc_t1 = bias_correct(file = t1, correction = "N4")
} else {
bc_t1 = readnii(out_fname)
}
```
or the filename (but negatives are in there):
```{r, eval = FALSE}
bc_t1 = bias_correct(file = "training01_01_t1.nii.gz", correction = "N4")
```
## Visualizing Bias Field Correction
Here we take the ratio of the images and overlay it on the original image:
```{r ratio_plot}
ratio = t1 / bc_t1; ortho2(t1, ratio)
```
```{r making_scales, echo = FALSE}
library(scales)
in_mask = (ratio < 0.999 | ratio > 1.0001) & ratio != 0
# get the quantiles
quantiles = quantile(ratio[ in_mask ], na.rm = TRUE,
probs = seq(0, 1, by = 0.1) )
quantiles = unique(quantiles)
# get a diverging gradient palette
fcol = scales::brewer_pal(type = "div", palette = "Spectral")(10)
# need one fewer color than breaks/quantiles
colors = gradient_n_pal(fcol)(seq(0,1, length = length(quantiles) - 1))
colors = scales::alpha(colors, 0.5) # colors are created
```
## Visualizing Bias Field Correction
We are breaking the ratio into quantiles:
```{r better_ratio_plot}
ortho2(t1, ratio, col.y = colors, ybreaks = quantiles, ycolorbar = TRUE)
```
## Conclusions
- Inhomogeneity correction is one of the first steps of most structural MRI pipelines
- Inhomogeneity can cause problems for other methods/segmentation
- Corrections try to make tissues of the same class to have similar intensities
- Use the `extrantsr` `bias_correct` function
- There is also `fsl_biascorrect` from `fslr` (not as effective in our experience)
- You may also want to run corrections after skull stripping on the brain only
- this is possible with the result after the brain extraction lecture
- correction before skull-stripping may be necessary and can improve after correction
## Website
http://johnmuschelli.com/imaging_in_r
## References