---
title: "Tutorial on anomaly detection, Part R: brain imaging"
author: "Pavlo Mozharovskyi"
output: html_notebook
---

This is an [R Markdown](http://rmarkdown.rstudio.com) Notebook on the tutorial 
on "Anomaly detection" given on Wednesday the 13th of April 2022.

# 1) Load brain data

## 1.1) Load R libraries and a few other functions to read DTI data

```{r}
library(curveDepth)
library(misc3d)
library(rgl)
library(AnalyzeFMRI)
source("DTI.R")
```

## 1.2) Load the brain imaging data

Load the anatomical brain volume data
```{r}
X <- f.read.volume("../data/101_1_dwi_fa.nii")
X.hdr <- f.read.nifti.header("../data/101_1_dwi_fa.nii")
X.hdr$dim
X <- X[, , , 1] # Since there is no time involved in DTI data
dim.X <- dim(X)
```

Load the DTI fibers data
```{r}
fibers.Left <- f.read.fibers(fibers.file = "../data/101_1_dwi.voxelcoordsL.txt")
fibers.Left.coord <- fibers.Left$fibers.coordinates
fibers.Right <- f.read.fibers(fibers.file = "../data/101_1_dwi.voxelcoordsR.txt")
fibers.Right.coord <- fibers.Right$fibers.coordinates
```

This data set contains two bundles (left and right) of 1,000 fibers each
```{r}
nrow(fibers.Left$indices.of.fibers)
nrow(fibers.Right$indices.of.fibers)
# But we will only display the first 100 fibers
n <- 100
leftFibers <- list("")
rightFibers <- list("")
```

```{r}
# Extract 100 first left fibers
for (j in 1:n) {
  f1 <- fibers.Left$indices.of.fibers[j]
  f2 <- fibers.Left$indices.of.fibers[j + 1] - 1
  x <- fibers.Left.coord[f1:f2, 1]
  y <- fibers.Left.coord[f1:f2, 2]
  z <- fibers.Left.coord[f1:f2, 3]
  leftFibers[[j]] <- list(voxels = 0, coords = cbind(x, y, z))
}
# Extract 100 first right fibers
for (j in 1:n) {
  f1 <- fibers.Right$indices.of.fibers[j]
  f2 <- fibers.Right$indices.of.fibers[j + 1] - 1
  x <- fibers.Right.coord[f1:f2, 1]
  y <- fibers.Right.coord[f1:f2, 2]
  z <- fibers.Right.coord[f1:f2, 3]
  rightFibers[[j]] <- list(voxels = 0, coords = cbind(x, y, z))
}
```

# 2) Calculate depths

```{r}
set.seed(1)
# The calculation can take quite some time (for both hemispheres, 
# depending on the computer speed, parameters 'm' and 'nDirs')
dssLeft <- depthc.Tukey(leftFibers, leftFibers, 
                        m = 50L, nDirs = 100, exactEst = FALSE)
dssRight <- depthc.Tukey(rightFibers, rightFibers, 
                         m = 50L, nDirs = 100, exactEst = FALSE)
```

# 3) Visualise 3D-interactive depth-based ordering

## 3.1) Create depth-based color palettes an determine color of each curve

```{r}
leftNColors <- length(dssLeft) + 1
my_palette <- colorRampPalette(c("yellow", "red"))(n = leftNColors)
dssLeft.colors <- my_palette[floor((dssLeft - min(dssLeft)) /
                                     (max(dssLeft) - min(dssLeft)) * 
                                     (leftNColors - 1)) + 1]
dssLeft.colors[which.max(dssLeft)] <- "blue"
rightNColors <- length(dssRight) + 1
my_palette <- colorRampPalette(c("yellow", "red"))(n = rightNColors)
dssRight.colors <- my_palette[floor((dssRight - min(dssRight)) /
                                      (max(dssRight) - min(dssRight)) * 
                                      (rightNColors - 1)) + 1]
dssRight.colors[which.max(dssRight)] <- "blue"

```

## 3.2) Display the brain volum using rgl
```{r}
contour3d(X, x = 1:dim.X[1], y = 1:dim.X[2], z = 1:dim.X[3], lev = 0.4, 
          alpha = 0.2, color = "grey", add = TRUE, fill = FALSE)
```

## 3.3) Visualise fibers
```{r}
# Plot left fibers
for (j in 1:n) {
  lines3d(leftFibers[[j]]$coords[, 1], 
          leftFibers[[j]]$coords[, 2], 
          leftFibers[[j]]$coords[, 3], 
          col = dssLeft.colors[j], add = TRUE, lwd = 3)
}
# Plot right fibers
for (j in 1:n) {
  lines3d(rightFibers[[j]]$coords[, 1], 
          rightFibers[[j]]$coords[, 2], 
          rightFibers[[j]]$coords[, 3], 
          col = dssRight.colors[j], add = TRUE, lwd = 3)
}

```

## 3.4) Change the orientation and viewpoint of the scene for suitable view
```{r}
rgl.viewpoint(theta = -25, phi = -45, fov = 30, zoom = 0.7)
```

# 4) Plot curve boxplot for the right bundle

## 4.1) Split/trim the curves of the right part of the brain w.r.t. their depth
```{r}
barplot(sort(dssRight, decreasing = TRUE))
abline(h = 0.1) # Not shown in the paper
iOrder <- order(dssRight, decreasing = TRUE)
iMostCentral <- which.max(dssRight)
iCentral <- order(dssRight, decreasing = TRUE)[1:floor(n / 2)]
iCentral <- iCentral[-1]
iOutsider <- which(dssRight < 0.1)
iOuter <- (1:length(dssRight))[-c(iMostCentral, iCentral, iOutsider)]
```

## 4.2) Plot the box-plot for the right part of the brain

```{r}
rgl.open()
rgl.bg(color = "white")
col1 <- c("royalblue", "red", "royalblue4")
```

Plot outer fibers
```{r}
for (i in iOuter){
  lines3d(rightFibers[[i]]$coords[, 1], 
          rightFibers[[i]]$coords[, 2], 
          rightFibers[[i]]$coords[, 3], 
          col = col1[1], add =TRUE, lwd = 3, alpha = 0.25)
}
```

Plot central fibers
```{r}
for (i in iCentral){
  lines3d(rightFibers[[i]]$coords[, 1], 
          rightFibers[[i]]$coords[, 2], 
          rightFibers[[i]]$coords[, 3], 
          col = col1[1], add = TRUE, lwd = 3, alpha = 0.75)
}
```

Plot outliers
```{r}
for (i in iOutsider){
  lines3d(rightFibers[[i]]$coords[, 1], 
          rightFibers[[i]]$coords[, 2], 
          rightFibers[[i]]$coords[, 3], 
          col = col1[2], add = TRUE, lwd = 3, alpha = 0.5)
}
```

Plot the deepest curve
```{r}
lines3d(rightFibers[[iMostCentral]]$coords[, 1], 
        rightFibers[[iMostCentral]]$coords[, 2], 
        rightFibers[[iMostCentral]]$coords[, 3], 
        col = col1[3], add = TRUE, lwd = 3, alpha = 1)
```

Change the orientation and viewpoint of the scene
```{r}
rgl.viewpoint(theta = -25, phi = -45, fov = 30, zoom = 0.6)
```

# ?) Try functional depths with length parametrisation

## For comparison, implement further depths, e.g. from R-package ddalpha

Add a new chunk by clicking the *Insert Chunk* button on the toolbar or by pressing *Cmd+Option+I*.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the *Preview* button or press *Cmd+Shift+K* to preview the HTML file). 

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike *Knit*, *Preview* does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.
