--- 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.