```{r setup, include=FALSE, message = FALSE}
library(methods)
knitr::opts_chunk$set(echo = TRUE, comment = "")
library(ggplot2)
library(ms.lesion)
library(neurobase)
library(extrantsr)
library(scales)
```


## Again we read in the data

```{r colon_twice, eval = FALSE}
t1 = neurobase::readnii("training01_01_t1.nii.gz")
```

```{r colon_twice_run, echo = FALSE}
t1 = neurobase::readnii("../training01_01_t1.nii.gz")
```

## Density of an Image

Let's do a marginal density of the values:

```{r dens}
plot(density(t1)) # large spike at 0
```


## Density with a mask

You can also pass in a mask to most standard functions:

```{r dens_with_mask}
plot(density(t1, mask = t1 > 0))
```

## Similarly: a Histogram

Note the high intensities:

```{r histog, echo = TRUE}
hist(t1)
```

## Orthographic view with additions

The `neurobase::ortho2` function displays nifti objects in 3 different planes.

```{r ortho2}
neurobase::ortho2(t1)
```

## Brightening up the image

We see a dark image; a this single large value affects how values are mapped.  The function `robust_window` calculates Winsorizes an image, by default the 0 (min) and 99.9th quantile, and sets values outside of this range to that quantile (97.5th below).  

```{r ortho2_rob}
ortho2(robust_window(t1, probs = c(0, 0.975)))
```

## Robust Density

Note the x-axis:

```{r dens_robust, echo = FALSE, fig.height = 5, fig.width = 10}
par(mfrow = c(1,2))
plot(density(t1), main = "Density of T1") 
plot(density(robust_window(t1, probs = c(0, 0.975))), main = "Density of Windowed T1")
par(mfrow = c(1,1))
```



## Overlaying images in `ortho2`

For the rest of the slides we will use the robust `t1` for plotting
```{r robust_final, echo = FALSE}
t1 = robust_window(t1, probs = c(0, 0.975))
```

Here we plot the T1 and a mask of values over $300$:

```{r ortho_nona}
ortho2(t1, y = t1 > 300)
```

## Double orthographic view

Sometimes you would like to represent 2 images side by side, of the same dimensions and orientation of course (useful for checking registration), use `double_ortho`

```{r double_ortho}
double_ortho(t1, y = t1 > 300, col.y = "white")
```

## Lightbox: view all slices

The `oro.nifti::image` function shows a lightbox view, all slices of an image:

```{r all_slices}
image(t1, useRaster = TRUE) # look at average brightness over each slice
```

## Viewing specific slices

The `slice` function can plot individual slices:

```{r two_slicewslice}
oro.nifti::slice(t1, z = c(60, 80))
```

## Different Planes

We can specify `z` the same way but change the `plane` to be different to get a different slice of the brain (could also do coronal):

```{r one_slice_sag}
oro.nifti::slice(t1, z = 125, plane = "sagittal")
```

## Overlaying slices

We can also overlay one slice of an image upon another using the `oro.nifti::slice_overlay` function.  
```{r one_slice_overlay}
slice_overlay(t1, y = t1 > 300, z = 80)
```

## Smoothing an Image (not extensively covered)

If you want to do 3D Gaussian smoothing, the `extrantsr::smooth_image` is helpful:

```{r smoothed}
library(extrantsr)
sm_img = smooth_image(t1, sigma = 2)
double_ortho(t1, sm_img)
```


## Conclusions

- `ortho2` - show orthographic images (and with overlays)
- `image` - shows multiple slices of an image
- `slice` - shows only specified slices
- `slice_overlay` - similar to `image` but with an overlay
- `double_ortho` - similar to `ortho2` but side-by-side
- `robust_window` - good for setting high values to not so high

## Website

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




<!-- ## `ms.lesion` Package -->

<!-- The `ms.lesion` package was made for this course.  It has all the data we will work with and the outputs from the analyses.   -->

<!-- The main function we will use is `get_image_filenames_list_by_subject`.  It returns a `list` of filenames of the images for the 5 training and 3 test subjects from and open MS data set [@lesjak2017novel] -->

<!-- - The full data is available at https://smart-stats-tools.org.  -->


<!-- ## `ms.lesion` Package -->

<!-- ```{r} -->
<!-- library(ms.lesion) -->
<!-- files = get_image_filenames_list_by_subject() -->
<!-- length(files); names(files);  -->
<!-- ``` -->
<!-- ```{r, eval = FALSE} -->
<!-- head(files$training01) -->
<!-- ``` -->
<!-- ```{r, echo = FALSE} -->
<!-- f = files$training01 -->
<!-- f = strsplit(f, "/") -->
<!-- f = sapply(f, function(x){ -->
<!--   ind = which(x == "library") -->
<!--   x = x[ind:length(x)] -->
<!--   x = paste(x, collapse = "/") -->
<!-- }) -->
<!-- head(f) -->
<!-- rm(list = "f") -->
<!-- ``` -->


<!-- ## Reading in a T1 image -->

<!-- Let's read in the T1 image from a MS lesion data set to visualize: -->

<!-- ```{r} -->
<!-- files = files$training01 -->
<!-- t1_fname = files["T1"] -->
<!-- t1 = readnii(t1_fname) -->
<!-- t1 -->
<!-- ``` -->


<!-- ## Changing the Windowing -->

<!-- The `zlim` option can also map intensities that can be plotted: -->

<!-- ```{r ortho2_zlim} -->
<!-- ortho2(t1, zlim = quantile(t1, probs = c(0, 0.999))) -->
<!-- ``` -->




<!-- ## Dropping empty dimensions -->

<!-- In some instances, there are extraneous slices to an image.  For example, in the Eve template image we read in, it is just the brain.  Areas of the skull and extracranial tissue are removed, but the slices remain so that the brain image and the original image are in the same space with the same dimensions.  For plotting or further analyses, we can drop these empty dimensions using the `neurobase::dropEmptyImageDimensions` function.  -->


<!-- By default, if one `nifti` is passed to the function and `keep_ind = FALSE`, then the return is a `nifti` object.  -->

<!-- ```{r dd, cache=FALSE} -->
<!-- reduced_mask = dropEmptyImageDimensions(qmask) -->
<!-- dim(qmask) -->
<!-- dim(reduced_mask) -->
<!-- ``` -->

<!-- ## Plotting the Reduced Image -->

<!-- ```{r plot_red_mask} -->
<!-- ortho2(reduced_mask) -->
<!-- ``` -->

<!-- ## Dropping empty dimensions -->

<!-- Many times, you want to get empty dimensions from some sort of mask/region of interest, then apply that dropping procedure to an image you want to visualize -->
<!-- - the `*EmptyImageDimensions` functions are good for removes "extraneous" slices for zooming in -->
<!-- ```{r dd_get, cache=FALSE} -->
<!-- keep_inds = getEmptyImageDimensions(qmask) -->
<!-- reduced_img = applyEmptyImageDimensions(t1, inds = keep_inds) -->
<!-- ortho2(reduced_mask) -->
<!-- ``` -->
