report.tex

\documentclass{article}
\usepackage[ampersand]{easylist}
\usepackage[margin=3.5cm]{geometry}
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage{color}
\usepackage{listings}
\usepackage{parskip}
\usepackage{xcolor}
\usepackage{wrapfig}

% We will all remember this as the point where the document prelude went to shit
\definecolor{darkgray}{gray}{0.3}

\PassOptionsToPackage{hyphens}{url}\usepackage[pdftex,
    colorlinks=true,
    linkcolor=darkgray,
    urlcolor=black,
    ]{hyperref}

% Set font to something a bit neater.
\usepackage{times}
\renewcommand{\familydefault}{\sfdefault}

\title{%
Brainstorm\\
\large Visualizing brain DTI data using particle systems}

\author{Vegard Itland --- Robin Grundvåg --- Stian Soltvedt}
\date{2018--12--14}

% Helps with overfull hboxes
\emergencystretch3em%
\hfuzz.5pt

% blatantly stolen from the INF226 thing
\definecolor{codebd}{RGB}{226,228,230}
\definecolor{codebg}{RGB}{246,248,250}
\definecolor{codefg}{RGB}{36,41,46}
\newcommand{\code}[1]{\fcolorbox{codebd}{codebg}{\lstinline[basicstyle=\ttfamily\color{codefg}]{#1}}}

\newcommand{\reference}[1]{[\hyperref[ref:#1]{\textbf{#1}}]}
\newcommand{\secref}[2]{\hyperref[sec:#1]{\textbf{#2}}}

\begin{document}

\maketitle
\pagenumbering{arabic}
\includegraphics[width=\textwidth]{brainstorm.png}

\newpage
\section*{Introduction}

Brainstorm is an application designed to visualize brain diffusion tensor imaging (DTI) data using particle fields. While particle visualization has been used in other contexts such as blood streams, streamlines are the prevalent visualization method for visualization of brain DTI data, and as far as we can tell nobody has yet attempted to visualize such datasets using particles.

The application is designed to be data-agnostic, meaning that it can read and visualize any dataset (in a supported file format) using particles, and is not limited to visualizing brain data.

\subsubsection*{DTI data}

Diffusion tensor imaging (DTI) is a medical imaging technique used to measure the interconnectivity of different parts of the brain. The actual measurements are of the direction of possible movement of water molecules in the brain, which in the brain's neural pathways are aligned with the direction of the path. Based on at least six measurements at a given point, we can obtain a tensor representing the vector field's directionality at that point by performing eigenanalysis on the measurement matrix. This enables us to construct a vector field using the measurement data.

\subsection*{Technology}

Brainstorm is written using Rust. Rust was chosen due to the need for a low-level language to easily work with OpenGL directly. While C/C++ is the obvious first candidate for such a project, we wished to avoid the error-prone manual memory management and pointer arithmetic that would require. Rust also easily compiles to WebAssembly, making it possible to also deploy Brainstorm to a website. This makes it more accessible to potential users without compromising on performance for the desktop version.

\section*{Implementation}

\subsection*{Making a cross-platform application}

\subsubsection*{Compiling to all platforms}

Rust compiles to all targets that LLVM can compile to, leading to broad platform support on a language level. This includes WebAssembly, although WebAssembly currently has many limitations. These include no threading support, no DOM access, and only being able to pass primitive types and pointers to JavaScript. There's also no API for basic functionality like accessing the date and time.

To get around this, we used the \code{stdweb} \reference{std} library. \code{stdweb} sets up the webpage and handles all calls between JavaScript and Rust automatically so there's no need to manually create external interfaces and callbacks. We also used the \code{cargo-web} \reference{car} tool to automate building the project as a website as well as running the project on a local web-server for testing.

The next obstacle for cross-platform support was rendering and GUI\@. Rust GUI libraries are few and immature, and our needs were very specific as Brainstorm requires direct access to OpenGL shaders. We also did not find any GUI libraries which could run on both OpenGL and WebGL\@. That left us to implement our own basic GUI toolkit from scratch.

To avoid having to maintain separate versions of code for both OpenGL and WebGL, we created a layer of abstraction over the two APIs. Fortunately, we were not the first to do this. \code{Kiss3D} \reference{Kis} is a simple 3D graphics engine written in Rust that can render to both OpenGL and WebGL, which our OpenGL/WebGL abstraction layer heavily borrows and copies from. This made OpenGL and WebGL calls identical for the rest of the code base with the correct bindings being determined at compile-time, making our codebase much simpler.

\subsubsection*{Writing a GUI toolkit from scratch}

\begin{wrapfigure}{r}{6cm}
  \centering
    \includegraphics[width=6cm]{ui.png}
    \raggedright%
    \textit{Fig. 1: The Brainstorm GUI}
    \vspace{1cm}
\end{wrapfigure}

Brainstorm has a set of data structures which handle drawing and texturing primitive triangles to the screen, which are used as building blocks to represent buttons and sliders. Text labels were challenging; in order to render text, a font must be loaded and parsed, the correct glyphs must be chosen, laid out and rasterized with sub-pixel positioning, and then cached on the GPU for performance. We used the \code{RustType} \reference{Rus} library to do most of this for us.


All GUI code is centralized in a separate GUI subsystem that handles  initialization, rendering and event handling for all UI elements. In order to make the UI elements interact with the rest of the application, a \code{State} struct is created that holds all the settings that can change at runtime, and pass around a reference to it to the various subsystems, such as the GUI and particle engines. Each UI element can then be loaded with a function that takes a mutable reference to the \code{State} struct and changes it according to the state and interactions with the UI element.


\subsubsection*{Handling files}

File picking is a difficult problem on desktop platform, let alone on the web. On the desktop side, Brainstorm uses the \code{nfd} \reference{nfd} crate that provides bindings to the C library \code{Native File Dialog} \reference{Nat}. This library abstracts away the differences between the native file picker API's on the three main desktop platforms --- Windows, Linux and Mac OS X.

The web required a completely different approach. Using \code{stdweb} for bindings, we created a small API in JavaScript that provided several functions:

\begin{easylist}
    \ListProperties(Space*=-3pt)
    & Open the file picker and continue.
    & Polling to see if a file has been picked.
    & Get the path of the last picked file.
\end{easylist}

This API lets Brainstorm retrieve the Base64-encoded file chosen by the user as soon as the browser is finished loading it. After that it is decoded to a regular byte stream and can be used it as if it was loaded from disk.

\subsection*{Data format}

This project uses the dataset of Gordon Kindlmann's brain \reference{Dif}. The dataset is in the NRRD data format, consisting of a header which provides information about the way the data is structured, and a body which provides the actual data. NRRD is a flexible data format, but our program currently only supports a subset of it which conforms closely to the format described in \reference{Dif}.

The choice of this data format was three-fold: First, it is simple to work with from a programming perspective, due to the ability to separate the header from the data body, the human and machine readable fields, and because the data is arranged in a very simple structure --- just a 3D grid of voxels. Secondly, one of the most easily accessible, open, well-documented brain DTI datasets was in this format. Lastly, the same page offered a generated helix sample dataset in the same format, allowing us to test that the methods we used would work with other datasets than just the brain.

To translate this data into a format directly understandable by our program, we made a preprocessor accepting NRRD data of the given format, which spits out the data as a serialized byte array which can be deserialized and used directly by our program. This enables us to store preprocessed data in a file for rapid loading, while also allowing the program to load NRRD data on the fly. Furthermore, it enables an easy way to extend the data in the preprocessing stage through further analysis or even manual inspection by experts before loading it into the program, which may be made capable of interpreting the extensions of this data (see \secref{prepro}{Pre-processing data analysis} for an example).

\begin{wrapfigure}{r}{6cm}
  \centering
    \includegraphics[width=6cm]{STgeneric}
    \raggedright%
    \textit{Fig. 2: Elliptical (anisotropic) tensor representation. CC-BY-SA 2.5 via \url{https://commons.wikimedia.org/wiki/File:STgeneric.png}}
    %\vspace{-2cm}
\end{wrapfigure}

The data body consists of an array of 32 bit floating point values in raw data, which are uniformly spaced measurements arrayed in memory in a 2D-slice-by-slice manner. The Z axis orders the slices such that higher Z memory indices indicate a higher spatial position of the slice. Thus, the coordinate (X, Y, Z) refers to the measurement in the Yth row and Xth column in the Zth slice of the data from the bottom. Our particle engine (in the initial rotation) uses the programming convention that the X and Y axes correspond to the computer screen's coordinate system, and Z indicates depth. This discrepancy was solved by rearranging the structure of the data in the preprocessing stage.

Each voxel in the dataset consists of a 7-dimensional vector of this format: (Confidence, Dxx, Dxy, Dxz, Dyy, Dyz, Dzz), where:

\begin{easylist}[itemize]
\ListProperties(Space*=-3pt)
& Confidence represents the likelihood of the measurement being signal (as opposed to noise). Usually this value is zero or one, but could in principle be anything between. For now, the preprocessor ignores any voxel with confidence lower than 1, making it essentially a mask of the brain.
& The other values form a matrix in this manner:
\end{easylist}

% Must be outside easylist for easylist not to go crazy
% hspace for indent
\hspace{1cm}
\(\begin{bmatrix}
    Dxx & Dxy & Dxz \\
    Dxy & Dyy & Dyz \\
    Dxz & Dyz & Dzz \\
\end{bmatrix}\)

Performing eigenanalysis on this matrix using the \code{nalgebra} \reference{Lin} library gives us the eigenvectors and the corresponding eigenvalues. Finally, we select the most significant eigenvector \((e_x,e_y,e_z)\), calculate the FA (fractional anisotropy) value using the eigenvalues, and push the \((e_x, e_y, e_z, fa)\) vector into our own 3D array in the order our program expects.


\subsubsection*{Pre-processing data analysis}\label{sec:prepro}

As an enhancement of the data, we attempted to analyze it in order to determine interesting areas from which to perform particle seeding. We decided to maximize particle travel length while also ensuring some spread in the selected points, to avoid that all points were clustered in a particular region.

The algorithm we chose uniformly samples points in the dataset and simulates the particle's path starting from that point using a fixed number of iterations (\(n = 1000\)). Starting points were calculated greedily based on the Euclidian distance of the end point to the starting point, and any starting point resulting in a path that collided with previously selected paths was discarded. This was done in order to encourage spread among the points we selected, so that we wouldn't select several close points which merged into one stream. We also attempted to give spread some additional weight by adding in a term considering the sum of the distance of the starting point to all previously selected points.

The final function we attempted to maximize was \(\textrm{dist}(p_e^i, p_s^i) + \textrm{dist}(p_e^i, p_s^i)\cdot{(\sum_{j=1}^{i-1}{\textrm{dist}(p_s^i, p_s^j)})}^{\sqrt{i}}\), where \(p_s^i\) and \(p_e^i\) indicates the start and end point of the \(i\)th selected seeding point, respectively. This formula was obtained experimentally, and may not be appropriate for all datasets. It ensures these properties:

\begin{easylist}[itemize]
\ListProperties(Space*=-3pt)
& The initial point only considers distance traveled
& The sum of the distance between points becomes increasingly important for later seeding points, increasing the spread of the points
\end{easylist}


Finally, in order to reduce the time it would take to process datasets, we parameterized some options allowing us to sacrifice resolution for speed:

\begin{easylist}[itemize]
\ListProperties(Space*=-3pt)
& The number of seeding points to calculate
& Step size of uniform sampling (assumes close values are likely to be similar)
& Threshold of product of FA values surrounding the point (assumes only strongly directional seeding points are worth considering)
\end{easylist}

Furthermore, the following optimizations were added, possibly sacrificing accuracy for speed:

\begin{easylist}[itemize]
\ListProperties(Space*=-3pt)
& We used a nearest neighbor interpolation scheme for computational simplicity\footnote{A recently discovered bug reveals that in fact, truncation was used instead of rounding (as nearest neighbor would do), but the result is presumed to be close enough for datasets of reasonably high resolution}
& If the particle's path reached a point where the FA value was zero, it would immediately disqualify the starting point path and move on
\end{easylist}

\begin{figure}[h]
    \makebox[5cm][l] {
        \includegraphics[width=5cm]{pts1}
    }
    \makebox[0.5\textwidth][l] {
        \includegraphics[height=5.92cm]{pts2}
    }\\
    \textit{Fig. 3: The generated seeding points seen from behind the brain (left) and on the side (right).}
\end{figure}

\subsubsection*{Output: Custom data format}

The data structure produced by the parser library is as follows:

\begin{verbatim}
#[derive(Serialize, Deserialize)]
struct VectorField {
    width: usize,
    depth: usize,
    height: usize,
    field: Vec<Vec<Vec<(f32,f32,f32,f32)>>>,
    seeding_pts: Vec<(f32,f32,f32)>
}
\end{verbatim}

Using the Rust crate \code{serde} \reference{Ser}, we can serialize this structure as a byte array (using the bincode format \reference{Bin}), and deserialize it for use in the application. The actual internal representation of the vector field in the application itself is flattened due to perceived performance benefits.

\subsection*{OpenGL and WebGL}

Making a renderer work on both the desktop and the web turned out to be trickier than we thought. WebGL's feature set is essentially a subset of OpenGL, with some subtle differences on top. These differences often took significant time and effort to track down and handle properly, but once we got the framework figured out, writing the actual algorithms was pretty simple.

\textbf{Note:} From here on, `OpenGL' will be used to refer to the OpenGL/WebGL binding framework.

\subsubsection*{CPU and GPU based particle rendering}

% NOTE: this is a pretty simplistic summary and doesn't account for lifetimes, high-pass filtering, low-pass filtering, etc.
The CPU based particle renderer was our first attempt at drawing particles. It simply keeps a list of particles in memory, which is updated and drawn once per frame. For each particle, its next position is found by performing a trilinear interpolation on the vector field based on the current particle position. After all particles are updated and moved, they're sent to OpenGL to be drawn as points. Then each point is expanded to a quad by changing the point size to reflect the actual particle size, and the particle is rounded to a circle with a fragment shader which calculates whether the fragment is within a certain radius from the center of the quad.

While the CPU based approach works well, it could not support as many particles as we wanted. The approach would also make streamlets (fading trails behind the particles) difficult to implement, and would further reduce the number of particles on-screen. One approach to increase it would be to add multi-threading, but this would not work on the web due to WebAssembly limitations. We therefore decided to try a GPU based approach.

\begin{figure}[h]
    \makebox[5cm][l] {
        \includegraphics[width=5cm]{cpu.png}
    }
    \makebox[0.5\textwidth][l] {
        \includegraphics[height=5.45cm]{gpu.png}
    }\\
    \textit{Fig. 4: CPU particles (left) and GPU particles (right)}
\end{figure}

% NOTE: Again, this is a simplified summary and does not account for any implementation details of the actual shaders used...
The GPU based particle rendering keeps track of an array of 3D textures, where the length of the array determines the length of the streamlets. Each texture represents the state of the particle system at a given frame, and each pixel represents the position of a single particle. In order to update the positions, we do the following:

\begin{easylist}
\ListProperties(Space*=-3pt)
& Pick a texture containing the current particle state and a 3D texture representing the vector field as input.
& For each particle, perform a trilinear interpolation in the 3D vector field texture to find the vector at the point it currently is.
& Add the vector to the current position and save the new particle position as the color on the output texture.
\end{easylist}

In order to display the streamlets in world space we have to do another render pass. We create another array of vertices that represent the vertices that should actually be drawn on the screen, as well as an array of indices based on the streamlets. These are passed to OpenGL along with the texture array that represents the particle state. We take the positions of the vertices from the colors in the output texture we save earlier and use those as the world space positions. Then we repeat that on the streamlets using the previous world states, based on the gl\_VertexID, and simply place each segment of the streamlet on the earlier calculated position.

\textbf{Note:} This is a simplified version of how it is implemented, as WebGL does not support read and write to the same texture, even if it is an array. Thus we are actually rendering back and forth between two texture arrays, but the logic is mostly the same (just a bit more tedious).

\subsubsection*{Marching Cubes mesh generation}

In order to draw the mesh of the vector field, we implemented the Marching Cubes algorithm. It creates a mesh based on the strength of the vectors in the vector field. This mesh only makes sense when the vector field has edges (like a brain or spiral, rather than noise fields), and we found it very helpful to easily see what part of the field we were looking at. The implementation is not all that special, and is very inspired by an article by Paul Bourke \reference{Mar}. We simply pass our vector field data and calculate the length of the vectors in order to create a scalar field and use it to construct a mesh.

\section*{Conclusion}

Looking at Brainstorm as it stands right now, it is potentially a useful tool for researchers wanting to visualize DTI scans in a new way. It creates an overview of the overall structure of the field, as well as letting the user explore the field, guided by the map sliders and filters. The data agnosticism also opens up for other applications of the program, although this is limited by the narrow set of data formats it currently supports. However, support for more formats can be added fairly trivially.

\subsection*{Future work}

Here are some ideas for enhancements and improvements to the Brainstorm application:

\begin{easylist}
\ListProperties(Space*=-3pt,Space=-3pt)
& Selecting between different interpolation schemes (e.g. Euler, Runge-Kutta).
& Additional sliders for particle counts, transparency, etc.
& Data generation to play with different random vector fields.
& Improving the high- and low-pass filters to ``obstruct'' particles physically rather than destroying and respawning them.
    && Likely a difficult task, but would likely significantly improve the visibility of nerve structures in the brain DTI scans.
& Visualizing probabilistic paths by using the less significant tensor vectors as an uncertainty measure.
    && Also a time-consuming task, but would leverage a strength of visualizing using particle fields rather than streamlines.
& Manual labeling of seeding points by experts.
\end{easylist}

\section*{References}

\textbf{[Dif]}\label{ref:Dif}
Diffusion tensor MRI datasets. \url{http://www.sci.utah.edu/~gk/DTI-data/}. (Accessed on 2018--12--14). Brain dataset courtesy of Gordon Kindlmann at the Scientific Computing and Imaging Institute, University of Utah, and Andrew Alexander, W. M. Keck Laboratory for Functional Brain Imaging and Behavior, University of Wisconsin-Madison.

\textbf{[std]}\label{ref:std}
A standard library for the client-side Web. \url{https://github.com/koute/stdweb}. (Accessed on 2018--12--17).

\textbf{[car]}\label{ref:car}
A cargo subcommand for the client-side Web. \url{https://github.com/koute/cargo-web}. (Accessed on 2018--12--17).

\textbf{[Kis]}\label{ref:Kis}
Kiss3D. \url{http://kiss3d.org/}. (Accessed on 2018--12--17).

\textbf{[Rus]}\label{ref:Rus}
RustType. \url{https://github.com/redox-os/rusttype}. (Accessed on 2018--12--17).

\textbf{[nfd]}\label{ref:nfd}
nfd-rs. \url{https://github.com/saurvs/nfd-rs}. (Accessed on 2018--12--17).

\textbf{[Nat]}\label{ref:Nat}
Native File Dialog. \url{https://github.com/mlabbe/nativefiledialog} (Accessed on 2018--12--17).

\textbf{[Lin]}\label{ref:Lin}
Linear algebra library. \url{https://www.nalgebra.org/}. (Accessed on 2018--12--17).

\textbf{[Ser]}\label{ref:Ser}
Serde. \url{https://serde.rs/}. (Accessed on 2018--12--17).

\textbf{[Bin]}\label{ref:Bin}
Bincode. \url{https://github.com/TyOverby/bincode}. (Accessed on 2018--12--17).

\textbf{[Mar]}\label{ref:Mar}
Polygonising a scalar field (Marching Cubes). \url{http://paulbourke.net/geometry/polygonise/}. (Accessed on 2018--12--16)

\end{document}