Chapter_2.lyx

#LyX 2.0 created this file. For more info see http://www.lyx.org/
\lyxformat 413
\begin_document
\begin_header
\textclass article
\begin_preamble
\numberwithin{figure}{section}
\numberwithin{table}{section}
\end_preamble
\use_default_options true
\maintain_unincluded_children false
\language english
\language_package default
\inputencoding auto
\fontencoding global
\font_roman palatino
\font_sans default
\font_typewriter default
\font_default_family default
\use_non_tex_fonts false
\font_sc false
\font_osf false
\font_sf_scale 100
\font_tt_scale 100

\graphics default
\default_output_format default
\output_sync 0
\bibtex_command default
\index_command default
\paperfontsize 12
\spacing onehalf
\use_hyperref false
\papersize default
\use_geometry true
\use_amsmath 1
\use_esint 1
\use_mhchem 1
\use_mathdots 1
\cite_engine basic
\use_bibtopic false
\use_indices false
\paperorientation portrait
\suppress_date false
\use_refstyle 0
\index Index
\shortcut idx
\color #008000
\end_index
\leftmargin 20page%
\topmargin 15page%
\rightmargin 15page%
\bottommargin 15page%
\secnumdepth 3
\tocdepth 3
\paragraph_separation indent
\paragraph_indentation default
\quotes_language english
\papercolumns 1
\papersides 1
\paperpagestyle default
\tracking_changes false
\output_changes false
\html_math_output 0
\html_css_as_file 0
\html_be_strict false
\end_header

\begin_body

\begin_layout Section
Cartesian Lattice 
\begin_inset Newline newline
\end_inset

Q-space Reconstructions
\begin_inset CommandInset label
LatexCommand label
name "sec:Cartesian-Lattice-Q-space"

\end_inset


\end_layout

\begin_layout Subsection
Overview
\end_layout

\begin_layout Standard
Between one to two thirds of imaging voxels in the human brain's white matter
 are thought to contain multiple fibre bundle crossings
\begin_inset space ~
\end_inset


\begin_inset CommandInset citation
LatexCommand cite
key "Behrens2007NeuroImage"

\end_inset

,
\begin_inset CommandInset citation
LatexCommand cite
key "Tournier2008"

\end_inset

 in which case the Diffusion Tensor model proposed by Basser et al.
\begin_inset space ~
\end_inset


\begin_inset CommandInset citation
LatexCommand cite
key "Basser1994"

\end_inset

 breaks down.
 High Angular Resolution Diffusion Imaging (HARDI) 
\begin_inset CommandInset citation
LatexCommand cite
key "Tuch2002"

\end_inset

, Diffusion Spectrum Imaging (DSI)
\begin_inset space ~
\end_inset


\begin_inset CommandInset citation
LatexCommand cite
key "callaghan1988nmr"

\end_inset

, 
\begin_inset CommandInset citation
LatexCommand cite
key "wedeen2005mapping"

\end_inset

 or Higher Order Tensors
\begin_inset space ~
\end_inset


\begin_inset CommandInset citation
LatexCommand cite
key "ozarslan2003generalized"

\end_inset

, 
\begin_inset CommandInset citation
LatexCommand cite
key "barmpoutis2009regularized"

\end_inset

 and many more reconstruction methods have been proposed to overcome the
 limitations of the Diffusion Tensor.
 These methods can be divided into those which need specific acquisition
 parametrizations, and those which can be used independently of q-space
 structure.
 For instance, for Q-ball Imaging
\begin_inset space ~
\end_inset


\begin_inset CommandInset citation
LatexCommand cite
key "Tuch2004"

\end_inset

 sampling needs to be on one or more spherical grids, and in Generalized
 Q-sampling Imaging (GQI)
\begin_inset space ~
\end_inset


\begin_inset CommandInset citation
LatexCommand cite
key "Yeh2010"

\end_inset

, requires sampling on a Cartesian grid; by contrast DTI can be used independent
ly of q-space structure.
 A further division considers the level of model assumptions for the diffusion
 process.
 Although all methods have some underlying assumptions we generally separate
 them in model-based and model-free.
 Model-based methods like the Single Tensor or Multi Tensor require a number
 of parameters to be fitted.
 By contrast, in model-free methods fitting is not necessary and the directional
ity of the underlying tissue can be approximated by some re-parametrization
 or re-transformation of the signal.
 The latter is usually more efficient than fitting models with many parameters
 which typically call for iterative methods.
\end_layout

\begin_layout Standard
This chapter presents, evaluates and compares different model-free methods
 for the reconstruction of orientation distribution functions using diffusion
 MRI data sampled on a Cartesian lattice in q-space.
 This non-parametric nature of the algorithms described here allows for
 the identification of multiple fibre crossings.
 In addition, a new method is presented named Diffusion Nabla Imaging (DNI)
 and a family of methods is defined called the Equatorial Inversion Transform
 (EIT).
 The EIT is a new way to represent and reconstruct the diffusion signal.
 Our results show that EIT can perform better or as well as the current
 state-of-the art methods i.e.
 DSI and GQI.
\end_layout

\begin_layout Subsection
Theory
\end_layout

\begin_layout Standard
We start from the classical formulation shown in Eq.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "eq:kq"

\end_inset

 of joint k-space and q-space imaging described in Callaghan
\begin_inset space ~
\end_inset


\begin_inset CommandInset citation
LatexCommand cite
key "Callaghan1991OUP"

\end_inset

, 
\begin_inset CommandInset citation
LatexCommand cite
key "callaghan1988nmr"

\end_inset

 using the narrow pulse gradient spin echo (PGSE) sequence of Tanner and
 Stejskal 
\begin_inset Formula 
\begin{eqnarray}
RF(\mathbf{k},\mathbf{q}) & = & \int\rho(\mathbf{v})\exp(i2\pi\mathbf{k}\cdot\mathbf{v})\int P_{\Delta}(\mathbf{v},\mathbf{r})\exp(i2\pi\mathbf{q}\cdot\mathbf{r})\, d\mathbf{r\,}d\mathbf{v}\label{eq:kq}
\end{eqnarray}

\end_inset


\end_layout

\begin_layout Standard
\begin_inset ERT
status open

\begin_layout Plain Layout


\backslash
noindent
\end_layout

\end_inset

 Here 
\begin_inset Formula $RF$
\end_inset

 represents the complex RF signal measured at spatial wave vector 
\begin_inset Formula $\mathbf{k}$
\end_inset

 and magnetic gradient wave vector 
\begin_inset Formula $\mathbf{q}$
\end_inset

, 
\begin_inset Formula $\rho$
\end_inset

 is the local spin density (number of protons per unit volume contributing
 to the RF signal), 
\begin_inset Formula $\Delta$
\end_inset

 is the time between diffusion gradients, 
\begin_inset Formula $P_{\Delta}$
\end_inset

 is the average diffusion propagator (transition probability distribution),
 
\series bold

\begin_inset Formula $\mathbf{v}$
\end_inset


\series default
 is the voxel coordinate, and 
\begin_inset Formula $\mathbf{r}$
\end_inset

 is the diffusion displacement.
 
\begin_inset Note Note
status collapsed

\begin_layout Plain Layout
the narrow pulse approximation allows us to easily understand the intrinsic
 Fourier relationship between the MR signal and the underlying water mobility
 that is a mirror of the tissue compartmentation.
 However, in practice, the duration of the diffusion-encoding gradient is
 not negligible compared to the diffusion time (i.e., ), so the formalism
 developed above must be reexamined.
 Interestingly, the reconstruction described in Eq.
 [4]
\end_layout

\end_inset


\end_layout

\begin_layout Standard
The k-space reconstruction with the narrow pulse approximation 
\begin_inset CommandInset citation
LatexCommand cite
key "Wedeen"

\end_inset

 gives us diffusion weighted image data 
\begin_inset Formula $S$
\end_inset

 which reveal the average propagator 
\begin_inset Formula $P_{\Delta}$
\end_inset

 of each voxel 
\begin_inset Formula 
\begin{eqnarray}
S(\mathbf{v},\mathbf{q}) & = & \int\rho(\mathbf{v})P_{\Delta}(\mathbf{v},\mathbf{r})\exp(i2\pi\mathbf{q}\cdot\mathbf{r})d\mathbf{r}\label{eq:W}
\end{eqnarray}

\end_inset


\end_layout

\begin_layout Standard
For the rest of the chapter we consider each voxel independently and assume
 intra-voxel spatial homogeneity so we can drop explicit reference to 
\begin_inset Formula $\mathbf{v}$
\end_inset

 and 
\begin_inset Formula $\Delta$
\end_inset

.
 We note in passing that the shape of 
\begin_inset Formula $P_{\Delta}$
\end_inset

 and hence of the ODF may change with different values of 
\begin_inset Formula $\Delta$
\end_inset

.
 We will not pursue this matter further here.
 We can also replace the spin density 
\family roman
\series medium
\shape up
\size normal
\emph off
\bar no
\noun off
\color none
\lang british

\begin_inset Formula $\rho(\mathbf{v})$
\end_inset


\family default
\series default
\shape default
\size default
\emph default
\bar default
\noun default
\color inherit
\lang english
 with 
\begin_inset Formula $S_{0}$
\end_inset

 i.e.
 the measured signal without diffusion weighting 
\begin_inset Formula $\mathbf{q}=\mathbf{0}$
\end_inset

.
 Therefore we can write
\begin_inset Formula 
\begin{eqnarray}
S(\mathbf{q}) & = & S_{0}\int P(\mathbf{r})\exp(i2\pi\mathbf{q}\cdot\mathbf{r})d\mathbf{r}\label{eq:S}
\end{eqnarray}

\end_inset


\end_layout

\begin_layout Standard
By applying the 3D Fourier transform in Eq.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "eq:S"

\end_inset

 we can reconstruct the average propagator also known as the diffusion spectrum
 
\begin_inset CommandInset citation
LatexCommand cite
key "Wedeen"

\end_inset

 or diffusion propagator
\begin_inset Formula 
\begin{eqnarray}
P(\mathbf{r}) & = & S_{0}^{-1}\int S(\mathbf{q})\exp(-i2\pi\mathbf{q}\cdot\mathbf{r})d\mathbf{q}\label{eq:P}
\end{eqnarray}

\end_inset


\end_layout

\begin_layout Standard
\begin_inset ERT
status open

\begin_layout Plain Layout


\backslash
noindent
\end_layout

\end_inset

 It was shown by Wedeen et al.
 
\begin_inset CommandInset citation
LatexCommand cite
key "Wedeen"

\end_inset

 that the dMRI signal is positive for any type of spin motion without net
 flux (i.e.
\begin_inset ERT
status open

\begin_layout Plain Layout

~
\end_layout

\end_inset

spin displacements due to thermal molecular agitation) or other random fluxes
 such as intravoxel incoherent motion.
 Under this assumption we can replace the complex signal 
\begin_inset Formula $S$
\end_inset

 with its modulus 
\begin_inset Formula $|S|$
\end_inset

 in Eq.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "eq:P"

\end_inset

 
\begin_inset Formula 
\begin{eqnarray}
P(\mathbf{r}) & = & S_{0}^{-1}\int|S(\mathbf{q})|\exp(-i2\pi\mathbf{q}\cdot\mathbf{r})d\mathbf{q}\label{eq:P_modulus}
\end{eqnarray}

\end_inset


\end_layout

\begin_layout Standard
The modulus of the signal coincides with the output of the standard MRI
 scanners as dMRI and that simplifies the acquisition procedure.
 It represents the density of the average relative spin displacement in
 a voxel.
 In other words, 
\begin_inset Formula $P(\mathbf{r})$
\end_inset

 is a measure of the probability that a spin in a chosen voxel , during
 the experimental mixing time 
\begin_inset Formula $\Delta$
\end_inset

, would make a vector displacement 
\begin_inset Formula $\mathbf{r}$
\end_inset

.
 We can visualize the propagator for every voxel as a 3D density volume
 (see Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:ODF"

\end_inset

).
\end_layout

\begin_layout Standard
In the classical DSI acquisition, at each location, diffusion-weighted images
 are acquired for 
\begin_inset Formula $N=515$
\end_inset

 or fewer values of q-encoding, comprising in q-space the points of a cubic
 lattice within the sphere of five lattice units in radius.
 Therefore,
\begin_inset Formula 
\begin{eqnarray}
\mathbf{q} & = & \alpha\mathbf{q}_{x}+\beta\mathbf{q}_{y}+\gamma\mathbf{q}_{z}\label{eq:q_lattice}
\end{eqnarray}

\end_inset


\end_layout

\begin_layout Standard
\noindent
\align block
with 
\begin_inset Formula $\alpha,\beta,\gamma\in\mathbb{Z}^{+}$
\end_inset

 and 
\begin_inset Formula $(\alpha^{2}+\beta^{2}+\gamma^{2})^{1/2}\leq5$
\end_inset

.
 The signal is premultiplied by a Hanning window before Fourier transform
 in order to ensure a smooth attenuation of the signal at high 
\begin_inset Formula $q$
\end_inset

 values.
 Often, to obtain data for the complete grid of 
\begin_inset Formula $515$
\end_inset

 q-vectors (which also means that we need to collect 
\begin_inset Formula $515$
\end_inset

 diffusion weighted volumes), the overall acquisition time would be too
 long and a smaller number of unique q-vectors are employed for just a single
 hemisphere usually between 
\begin_inset Formula $101$
\end_inset

 to 
\begin_inset Formula $257$
\end_inset

 points
\begin_inset space ~
\end_inset


\begin_inset CommandInset citation
LatexCommand cite
key "Kuo"

\end_inset

.
 This is valid because the underlying self-diffusion process is symmetric
 and so the signal is symmetric, therefore the vectors can be mapped on
 the other hemisphere to create the full q-space.
 
\end_layout

\begin_layout Standard
Since we are mainly interested in the angular structure of the underlying
 tissue, we further simplify the data by taking the weighted radial summation
 of 
\begin_inset Formula $P(\mathbf{r})$
\end_inset


\lang british

\begin_inset Formula 
\begin{equation}
\psi_{DSI}(\hat{\mathbf{u}})=\int_{0}^{\infty}P(r\hat{\mathbf{u}})r^{2}dr\label{eq:ODF_DSI}
\end{equation}

\end_inset


\end_layout

\begin_layout Standard
\begin_inset ERT
status open

\begin_layout Plain Layout


\backslash
noindent
\end_layout

\end_inset

 This defines the 
\lang british
orientation density function (ODF) for DSI which measures the quantity of
 diffusion in the direction of the unit vector 
\begin_inset Formula $\mathbf{\hat{u}}$
\end_inset

 where 
\begin_inset Formula $\mathbf{r=}r\hat{\mathbf{u}}$
\end_inset

.
\end_layout

\begin_layout Standard

\lang british
Note at this point that in order to find the ODF we have to first create
 the diffusion propagator by 
\lang english
applying
\lang british
 the Fourier transform on the lattice.
 
\lang english
Yeh et al.
\begin_inset space ~
\end_inset


\begin_inset CommandInset citation
LatexCommand cite
key "Yeh2010"

\end_inset

 proposed a direct way to calculate a slightly different ODF using the Cosine
 transform.
\end_layout

\begin_layout Standard
\align block
In order to represent the average propagator in the scale of spin quantity
 Yeh et al.
\begin_inset space ~
\end_inset


\begin_inset CommandInset citation
LatexCommand cite
key "Yeh2010"

\end_inset

 introduced the 
\emph on
spin density function 
\emph default

\begin_inset Formula $Q$
\end_inset

 which is estimated by scaling the average propagator 
\begin_inset Formula $P_{\Delta}$
\end_inset

 with the spin density 
\begin_inset Formula $\rho$
\end_inset

, i.e.
 
\begin_inset Formula $Q(\mathbf{r})=\rho P(\mathbf{r})=S_{0}P(\mathbf{r})$
\end_inset

.
 From Eq.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "eq:S"

\end_inset

 we obtain
\family roman
\series medium
\shape up
\size normal
\emph off
\bar no
\noun off
\color none
\lang british

\begin_inset Formula 
\begin{eqnarray}
S(\mathbf{q}) & = & \int Q(\mathbf{r})\exp(i2\pi\mathbf{q}\cdot\mathbf{r})d\mathbf{r}\label{eq:W2Q}
\end{eqnarray}

\end_inset


\family default
\series default
\shape default
\size default
\emph default
\bar default
\noun default
\color inherit
\lang english
We can apply the Fourier transform again to Eq.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "eq:W2Q"

\end_inset

 and obtain
\begin_inset Formula 
\begin{eqnarray}
Q(\mathbf{r}) & = & \int S(\mathbf{q})exp(-i2\pi\mathbf{q}\cdot\mathbf{r})d\mathbf{q}\label{eq:Q2S_complex}
\end{eqnarray}

\end_inset

Because
\family roman
\series medium
\shape up
\size normal
\emph off
\bar no
\noun off
\color none
\lang british
 
\begin_inset Formula $Q(\mathbf{r})$
\end_inset


\family default
\series default
\shape default
\size default
\emph default
\bar default
\noun default
\color inherit
\lang english
 is real and 
\begin_inset Formula $S(\mathbf{q})$
\end_inset

 is symmetric (even), i.e.
 
\begin_inset Formula $S(\mathbf{q})=S(-\mathbf{q})$
\end_inset

, we can use directly the Fourier Cosine transform (see section 
\begin_inset CommandInset ref
LatexCommand ref
reference "sub:The-cosine-transform"

\end_inset

) to calculate
\begin_inset Formula 
\begin{eqnarray}
Q(\mathbf{r}) & = & \int S(\mathbf{q})cos(2\pi\mathbf{q}\cdot\mathbf{r})d\mathbf{q}\label{eq:cosine_transform}
\end{eqnarray}

\end_inset


\end_layout

\begin_layout Standard
\noindent
and obtain the 
\begin_inset Quotes eld
\end_inset

spin
\begin_inset Quotes erd
\end_inset

 orientation distribution function (SDF) 
\begin_inset Formula $\psi_{GQI}$
\end_inset

 from an unweighted truncated radial projection
\end_layout

\begin_layout Standard
\begin_inset Formula 
\begin{eqnarray}
\psi_{GQI}(\mathbf{\hat{u}}) & = & \intop_{0}^{\lambda}Q(r\mathbf{\hat{u}})dr\label{eq:SDF}\\
 & = & \intop_{0}^{\lambda}\int S(\mathbf{q})\cos(2\pi r\mathbf{q}\cdot\mathbf{\hat{u}})d\mathbf{q}dr\\
 & = & \lambda\int S(\mathbf{q})\mathtt{sinc}(2\pi\lambda\mathbf{q}\cdot\mathbf{\hat{u}})d\mathbf{q}
\end{eqnarray}

\end_inset


\end_layout

\begin_layout Standard
\noindent
\align block
\begin_inset ERT
status open

\begin_layout Plain Layout


\backslash
noindent
\end_layout

\end_inset

 where 
\begin_inset Formula $\lambda$
\end_inset

 is a constant called the diffusion sampling length.
 This parameter acts as a smoothing factor.
 The higher 
\begin_inset Formula $\lambda$
\end_inset

 the more detailed the SDF will be but also more noisy.
 This ODF was used as the basis of the analysis of the GQI method.
 It provides a simple direct analytical solution of the ODF which can be
 written in a simple matrix form
\begin_inset Formula 
\begin{eqnarray}
\bm{\psi}_{GQI}= & \mathbf{s}\cdot\mathtt{sinc}((6D\cdot G\circ\mathbf{b}\circ\mathbb{1})\cdot U^{T})\lambda/\pi\label{eq:GQI_analytical}
\end{eqnarray}

\end_inset


\end_layout

\begin_layout Standard
\noindent
\align block
where 
\begin_inset Formula $\cdot$
\end_inset

 denotes standard matrix or vector dot product, 
\family roman
\series medium
\shape up
\size normal
\emph off
\bar no
\noun off
\color none
\lang british

\begin_inset Formula $\circ$
\end_inset

 
\family default
\series default
\shape default
\size default
\emph default
\bar default
\noun default
\color inherit
\lang english
denotes the Hadamard product, 
\begin_inset Formula $\bm{\psi}{}_{GQI}$
\end_inset

 as a 
\begin_inset Formula $M$
\end_inset

-dimensional vector with components corresponding to the selected directions
 
\begin_inset Formula $\hat{\mathbf{u}}$
\end_inset

 on the ODF sphere, 
\begin_inset Formula $\mathbf{s}$
\end_inset

 is a vector with all the signal values, D=0.00251
\begin_inset space ~
\end_inset


\begin_inset CommandInset citation
LatexCommand cite
key "canalesrodriguez2009mdq"

\end_inset

 
\family roman
\series medium
\shape up
\size normal
\emph off
\bar no
\noun off
\color none
\lang british
where 
\begin_inset Formula $D$
\end_inset

 is a constant known as the free water diffusion coefficient
\family default
\series default
\shape default
\size default
\emph default
\bar default
\noun default
\color inherit
\lang english
, 
\begin_inset Formula $G$
\end_inset

 is the 
\begin_inset Formula $N\times3$
\end_inset

 matrix with the gradient vectors, 
\begin_inset Formula $\mathbf{b}$
\end_inset

 is the 
\begin_inset Formula $N\times1$
\end_inset

 matrix with the b-values and 
\begin_inset Formula $\mathbb{1}$
\end_inset

 is the 
\begin_inset Formula $N\times3$
\end_inset

 incidence matrix where all values are equal to 
\begin_inset Formula $1$
\end_inset

.
\end_layout

\begin_layout Standard
For a similar ODF like the one produced using DSI we need to take the weighted
 truncated radial projection.
 This will give us a different 
\begin_inset Quotes eld
\end_inset

spin
\begin_inset Quotes erd
\end_inset

 ODF which we symbolize with 
\begin_inset Formula $\psi_{GQI_{2}}$
\end_inset


\begin_inset Formula 
\begin{eqnarray}
\psi_{GQI2}(\mathbf{\hat{u}}) & = & \intop_{0}^{\lambda}Q(r\mathbf{\hat{u}})r^{2}dr\label{eq:SDF2}\\
 & = & \lambda^{3}\int S(\mathbf{q})H(2\pi\lambda\mathbf{q}\cdot\mathbf{\hat{u}})d\mathbf{q}
\end{eqnarray}

\end_inset


\end_layout

\begin_layout Standard
\begin_inset ERT
status open

\begin_layout Plain Layout


\backslash
noindent
\end_layout

\end_inset

 where 
\begin_inset Formula $H(x)=\begin{cases}
\frac{2\cos(x)}{x^{2}} & +\frac{(x^{2}-2)\sin(x)}{x^{3}},x\neq0\\
 & 1/3\qquad\qquad,x=0
\end{cases}$
\end_inset

.
\end_layout

\begin_layout Standard
\align left
This equation can be similarly implemented with a simple matrix transform
\begin_inset Formula 
\begin{eqnarray*}
\bm{\psi}_{GQI2} & = & \mathbf{s}\cdot\mathtt{H}((6D\cdot G\circ\mathbf{b}\circ\mathbb{1})\cdot U^{T})\lambda^{3}/\pi
\end{eqnarray*}

\end_inset


\end_layout

\begin_layout Standard
\align left
and has not to date been published with real or simulated data sets.
\end_layout

\begin_layout Standard
The addition of the spin density plays a very important role on normalizing
 the ODF and providing scalar or vector metrics for the analysis of dMRI
 data sets.
 GQI, similarly to DSI, expects the q-vectors to sit on a cubic lattice
 within a sphere.
 However, because of the direct analytical formulation of the GQI ODFs;
 the creation of the volumetric grid with the signal values is not necessary.
 This makes GQI advantageous on memory and CPU efficiency.
 Furthermore, no Hanning filter is necessary.
\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename last_figures/triple_crossing_gqi_gqi2_dsi_dni.png
	lyxscale 30
	scale 30

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption

\begin_layout Plain Layout
Showing the ODFs from two randomly oriented simulated 
\begin_inset Formula $3$
\end_inset

-fibre crossings at 
\begin_inset Formula $50^{\circ}$
\end_inset

(top) and 
\begin_inset Formula $90^{\circ}$
\end_inset

 angles between each pair of fibres using different Cartesian lattice q-space
 reconstruction methods.
 
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset CommandInset label
LatexCommand label
name "Flo:beautiful-triple-crossing"

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Subsection
Other methods
\end_layout

\begin_layout Standard
Pickalov et al.
\begin_inset space ~
\end_inset


\begin_inset CommandInset citation
LatexCommand cite
key "pickalov2006tra"

\end_inset

 proposed a new method for reconstructing the diffusion propagator by applying
 an iterative inverse Radon transform on measurements along many radial
 lines; computing 1D tomographic projections to reconstruct the 3D EAP.
 This technique measures DW images along a few radial lines of q-space but
 still requires hundreds of samples to reliably recover the EAP.
 Currently, to reconstruct the EAP, the state-of-the-art model-free techniques
 apart from diffusion spectrum imaging are hybrid diffusion imaging (HYDI)
\begin_inset space ~
\end_inset


\begin_inset CommandInset citation
LatexCommand cite
key "wu2007hybrid"

\end_inset

 and multiple q-shell diffusion propagator imaging (mq-DPI)
\begin_inset space ~
\end_inset


\begin_inset CommandInset citation
LatexCommand cite
key "descoteaux2010multiple"

\end_inset

.
 HYDI acquires the signal values on five concentric spherical q-space shells,
 then interpolates onto a cubic grid and applies the standard Fourier transform
 in the same way as DSI.
 In mq-DPI the EAP is calculated by solving Laplace's equation for the diffusion
 signal using a real and symmetric modified spherical harmonic basis.
 The EAP can be found analytically by the inversion of a linear system using
 Laplace-Beltrami regularization.
 In addition, Exact Q-ball imaging (EQBI)
\begin_inset space ~
\end_inset


\begin_inset CommandInset citation
LatexCommand cite
key "canalesrodriguez2009mdq"

\end_inset

 provides a different method to calculate the ODF analytically using multiple
 spherical q-space shells.
 Similarly, Aganj et al.
\begin_inset space ~
\end_inset


\begin_inset CommandInset citation
LatexCommand cite
key "Aganj2010"

\end_inset

, proposed an analytical solution for the multi-shell case by incorporating
 a mono-exponential or bi-exponential model (CSA-ODF).
 Another method for finding a directional distribution on the sphere was
 proposed by Özarslan et al.
\begin_inset space ~
\end_inset


\begin_inset CommandInset citation
LatexCommand cite
key "ozarslan2006resolution"

\end_inset

 called the diffusion orientation transform (DOT).
 This method calculates a different statistic 
\begin_inset Formula $P(r_{0}\mathbf{\hat{u}})$
\end_inset

, the probability of finding a particle initially at the origin at the point
 
\begin_inset Formula $r_{0}\mathbf{\hat{u}}$
\end_inset

, using spherical harmonics.
 Not surprisingly there is a relationship connecting CSA with DOT which
 is 
\lang british

\begin_inset Formula 
\begin{equation}
\psi_{CSA}(\hat{\mathbf{u}})=\int_{0}^{\infty}DOT(r\hat{\mathbf{u}})r^{2}dr\label{eq:CSAandDOT}
\end{equation}

\end_inset


\end_layout

\begin_layout Standard

\lang british
J
\lang english
ansons et al.
\begin_inset space ~
\end_inset


\begin_inset CommandInset citation
LatexCommand cite
key "JansonsPAS2003"

\end_inset

 proposed a different function on the sphere than the ODFs described above,
 to be used on data sets acquired on a single spherical q-space shell.
 They called this spherical function persistent angular structure (PAS).
 This method has very good angular resolution.
 It uses the principle of maximum entropy however, it is rather slow as
 it nonlinear fitting is used in order to identify many parameters.
 PAS is a statistic on the sphere defined as 
\begin_inset Formula $PAS(\hat{\mathbf{u}})=\exp(\lambda_{0}+\sum_{j=1}^{N}\lambda_{j}\cos(\mathbf{q}_{j}\cdot k\hat{\mathbf{u}}))$
\end_inset

 where 
\begin_inset Formula $\lambda$
\end_inset

 are the unknown parameters, 
\begin_inset Formula $k$
\end_inset

 is constant and 
\begin_inset Formula $N$
\end_inset

 is the number of DWIs.
 The relationship 
\begin_inset Formula $\int PAS(\hat{\mathbf{u}})\exp(i\mathbf{q}_{j}\cdot k\hat{\mathbf{u}})d\hat{\mathbf{u}}=E(\mathbf{q}_{j})$
\end_inset

 provides the bridge between PAS and the diffusion signal (
\begin_inset Formula $E(\mathbf{q})$
\end_inset

).
\end_layout

\begin_layout Standard
The first publication of using spherical harmonic expansions with diffusivity
 profiles, which are now quite common in the literature, was by Alexander
 et al.
\begin_inset space ~
\end_inset


\begin_inset CommandInset citation
LatexCommand cite
key "alexander2002detection"

\end_inset

.
 Q-ball imaging was introduced by Tuch
\begin_inset space ~
\end_inset


\begin_inset CommandInset citation
LatexCommand cite
key "Tuch2004"

\end_inset

 and a new ODF defined as 
\begin_inset Formula $\psi(\mathbf{\hat{u}})=\frac{1}{Z}\int_{0}^{\infty}P(r\mathbf{\hat{u}})dr$
\end_inset

 where 
\begin_inset Formula $Z$
\end_inset

 is a normalization constant.
 It was later provided for Q-Ball imaging a fast and analytical solution
 using spherical harmonics (SH) and Laplace-Beltrami regularization 
\begin_inset CommandInset citation
LatexCommand cite
key "Descoteaux2007"

\end_inset

.
 Tournier et al.
 
\begin_inset CommandInset citation
LatexCommand cite
key "tournier2004direct"

\end_inset

, 
\begin_inset CommandInset citation
LatexCommand cite
key "Tournier2007"

\end_inset

 introduced a spherical deconvolution method where first the SH coefficients
 were estimated, then single fiber ODFs were used as a deconvolution kernel
 estimated from the real data.
 Then, the sharper fODF (fiber orientation distribution function) was obtained
 by a simple linear transformation 
\begin_inset CommandInset citation
LatexCommand cite
key "descoteaux2007deterministic"

\end_inset

.
 Other deconvolution approaches were proposed in 
\begin_inset CommandInset citation
LatexCommand cite
key "sakaie2007objective"

\end_inset

 and 
\begin_inset CommandInset citation
LatexCommand cite
key "yeh2011estimation"

\end_inset

.
\end_layout

\begin_layout Standard
On Tensor related methods we have the classical Single Tensor 
\begin_inset CommandInset citation
LatexCommand cite
key "Basser1994"

\end_inset

, Sticks and Ball
\begin_inset CommandInset citation
LatexCommand cite
key "Behrens2007NeuroImage"

\end_inset

, Multi-Tensor 
\begin_inset CommandInset citation
LatexCommand cite
key "pasternak2008variational"

\end_inset

, 
\begin_inset CommandInset citation
LatexCommand cite
key "liu2004characterizing"

\end_inset

 and Higher Rank Tensors 
\begin_inset CommandInset citation
LatexCommand cite
key "ozarslan2003generalized"

\end_inset

, 
\begin_inset CommandInset citation
LatexCommand cite
key "barmpoutis2009regularized"

\end_inset

.
 In addition there are also model based methods which try to calculate non-Gauss
ian properties, for example the Kurtosis Tensor
\begin_inset space ~
\end_inset


\begin_inset CommandInset citation
LatexCommand cite
key "jensen2005diffusional"

\end_inset

,
\begin_inset space ~
\end_inset


\begin_inset CommandInset citation
LatexCommand cite
key "lu2006three"

\end_inset

 which is used in Diffusion Kurtosis Imaging (DKI).
\end_layout

\begin_layout Standard
Finally, new model-based methods are emerging which are trying to calculate
 statistics like the axonal thickness distribution from dMRI data sets.
 These are usually based on model free and restricted components; CHARMED
\begin_inset space ~
\end_inset


\begin_inset CommandInset citation
LatexCommand cite
key "assaf2005composite"

\end_inset

, 
\begin_inset CommandInset citation
LatexCommand cite
key "assaf2004new"

\end_inset

, AxCaliber
\begin_inset space ~
\end_inset


\begin_inset CommandInset citation
LatexCommand cite
key "assaf2008axcaliber"

\end_inset

 and the orientation invariant ActiveAx 
\begin_inset CommandInset citation
LatexCommand cite
key "alexander2010orientationally"

\end_inset

 are some well known methods of this type.
 Q-space Imaging (QSI) can be used to identify distributions of axon-diameter
 too 
\begin_inset CommandInset citation
LatexCommand cite
key "ong2010quantifying"

\end_inset

.
\end_layout

\begin_layout Subsection
Diffusion Nabla Imaging
\end_layout

\begin_layout Standard
A new method for the calculation of the real ODF is proposed here.
 This is based on the theoretical work done by Aganj et al.
\begin_inset space ~
\end_inset


\begin_inset CommandInset citation
LatexCommand cite
key "aganj2010reconstruction"

\end_inset

 and Canales-Rodriguez et al.
\begin_inset space ~
\end_inset


\begin_inset CommandInset citation
LatexCommand cite
key "Canales-Rodriguez2009"

\end_inset

 using two important theorems from Fourier Analysis
\end_layout

\begin_layout Enumerate
The Fourier transform of 
\begin_inset Formula $P(\mathbf{r})r^{2}=-\nabla^{2}E(\mathbf{q})$
\end_inset

 where 
\begin_inset Formula $\nabla^{2}$
\end_inset

 is the Laplacian operator (for proof see section 
\begin_inset CommandInset ref
LatexCommand ref
reference "sub:Fourier-transform-ofPr2"

\end_inset

).
\end_layout

\begin_layout Enumerate
For a symmetric function 
\begin_inset Formula $E:\mathbb{R}^{3}\rightarrow\mathbb{R}$
\end_inset

 and for the arbitrary unit vector 
\begin_inset Formula $\hat{\mathbf{u}}$
\end_inset

 we have 
\begin_inset Formula $\intop_{0}^{\infty}E(r\hat{u})dr=\frac{1}{8\pi^{2}}\int\int_{\hat{\mathbf{u}}^{\perp}}E(q)\, q\, dq\, d\phi$
\end_inset

 where 
\begin_inset Formula $\hat{\mathbf{u}}^{\perp}$
\end_inset

 is the plane perpendicular to 
\begin_inset Formula $\hat{\mathbf{u}}$
\end_inset

 (for proof see section 
\begin_inset CommandInset ref
LatexCommand ref
reference "sub:Radial-projection-of-symmetric"

\end_inset

).
\end_layout

\begin_layout Standard
From Eq.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "eq:ODF_DSI"

\end_inset

 we see that the integration is over 
\begin_inset Formula $P(r\hat{\mathbf{u}})r^{2}$
\end_inset

 , therefore we can write
\end_layout

\begin_layout Standard

\lang british
\begin_inset Formula 
\begin{equation}
\psi_{DNI}(\hat{\mathbf{u}})=-\frac{1}{8\pi^{2}}\int_{\hat{\mathbf{u}}^{\perp}}\int_{0}^{\infty}\nabla^{2}E(q)\, q\, dq\, d\phi\label{eq:ODF_DNI}
\end{equation}

\end_inset


\end_layout

\begin_layout Standard
\noindent
\align block

\lang british
where 
\begin_inset Formula $\phi$
\end_inset

 is the angular rotation component operating on the plane perpendicular
 to 
\begin_inset Formula $\hat{\mathbf{u}}$
\end_inset

, 
\family roman
\series medium
\shape up
\size normal
\emph off
\bar no
\noun off
\color none
\lang english

\begin_inset Formula $\nabla^{2}$
\end_inset

 is the Laplacian operator and 
\begin_inset Formula $E(q)=S(q)/S_{0}$
\end_inset

 is the normalized diffusion signal.
 Eq.
 
\begin_inset CommandInset ref
LatexCommand ref
reference "eq:ODF_DNI"

\end_inset

 has the advantage that no Fourier transform is necessary.
 We need however, a way to calculate the Laplacian of the signal.
 This can be analytically derived for a spherical grid 
\family default
\series default
\shape default
\size default
\emph default
\bar default
\noun default
\color inherit

\begin_inset CommandInset citation
LatexCommand cite
key "aganj2010reconstruction"

\end_inset

 and we propose here that it can be directly calculated in a cubic grid
 using the standard 3D discrete Laplacian filter which is given by the 3D
 kernel defined by the following 
\begin_inset Formula $3\times3\times3$
\end_inset

 array: 
\end_layout

\begin_layout Standard
\begin_inset Formula 
\[
\left[\left(\begin{array}{ccc}
0 & 0 & 0\\
0 & 1 & 0\\
0 & 0 & 0
\end{array}\right),\left(\begin{array}{ccc}
0 & \:\,1 & \:\,0\\
1 & -6 & \:\,1\\
0 & \:\,1 & \:\,0
\end{array}\right),\left(\begin{array}{ccc}
0 & 0 & 0\\
0 & 1 & 0\\
0 & 0 & 0
\end{array}\right)\right].
\]

\end_inset


\end_layout

\begin_layout Standard
\noindent
\align block
This is a filter commonly used for image processing.
 From now on when we use the Laplacian operator in order to measure the
 directionality of the diffusion signal we will call this reconstruction
 method Diffusion Nabla Imaging as nabla-squared 
\family roman
\series medium
\shape up
\size normal
\emph off
\bar no
\noun off
\color none
(
\begin_inset Formula $\nabla^{2}$
\end_inset

) 
\family default
\series default
\shape default
\size default
\emph default
\bar default
\noun default
\color inherit
is the symbol for the Laplacian operator.
 In Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:beautiful-triple-crossing"

\end_inset

 we present the ODFs from two randomly oriented simulated 
\begin_inset Formula $3$
\end_inset

-fibre crossings at 
\begin_inset Formula $50^{\circ}$
\end_inset

 and 
\begin_inset Formula $90^{\circ}$
\end_inset

 angles between each other using different grid based reconstruction methods.
 The parameters used here are for DSI: radial sampling range 
\begin_inset Formula $2.1-6$
\end_inset

 with 
\begin_inset Formula $0.2$
\end_inset

 and Hanning filter width 
\begin_inset Formula $36$
\end_inset

, GQI: 
\begin_inset Formula $\text{\lambda=1.2}$
\end_inset

, GQI2: 
\begin_inset Formula $\lambda=3$
\end_inset

 and DNI: radial sampling 
\begin_inset Formula $0-5$
\end_inset

 with 
\begin_inset Formula $0.2$
\end_inset

 steps.
 All methods used the same reconstruction sphere with 
\begin_inset Formula $642$
\end_inset

 vertices and 
\begin_inset Formula $1,280$
\end_inset

 faces.
 
\end_layout

\begin_layout Subsection
Equatorial Inversion Transform
\end_layout

\begin_layout Standard
We propose an important theoretical construction called the Equatorial Inversion
 Transform (EIT) which creates a general formulation for the interpretation
 of the directionality of the diffusion signal.
 This idea is founded on two general properties of the diffusion signal:
 (a) If we visualize the diffusion signal for a single fibre for all gradient
 directions we see the generated shape to be lowest towards the direction
 of the fibre and highest on the plane perpendicular to that direction (see
 Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:single_fibre_spherical_grid"

\end_inset

).
 (b) The diffusion signal is additive i.e.
 
\begin_inset Formula $S(\hat{\mathbf{f}}_{1})+S(\hat{\mathbf{f}}_{2})=S(\hat{\mathbf{f}}_{1}+\hat{\mathbf{f}}_{2})$
\end_inset

, where 
\begin_inset Formula $\hat{\mathbf{f}}_{1}$
\end_inset

, 
\begin_inset Formula $\hat{\mathbf{f}}_{2}$
\end_inset

 are the unit directions of the fibres.
 In simple terms the signal of 2-fibre crossing can be decomposed linearly
 to the signals of the two fibres that create the crossing.
 The same holds for any number of fibres in a crossing.
\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename last_figures/single_fiber_signal_sphere_grid2.png
	lyxscale 30
	scale 50

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption

\begin_layout Plain Layout
The diffusion signal has the beautiful property that it is minimum along
 the direction of a fibre with unit direction 
\begin_inset Formula $\mathbf{\hat{\bm{f}}}$
\end_inset

 and maximum along the equator defined by the plane perpendicular to that
 fibre direction.
 This property is the basic inspiration behind the EIT.
 In this picture the 3D surface plot of a simulated signal for a spherical
 grid acquisition with b-value 2,000 is shown using a yellow-red colourmap.
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset CommandInset label
LatexCommand label
name "Flo:single_fibre_spherical_grid"

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
These are two very important geometric properties of the signal that we
 can try to exploit to its limit by calculating equatorial integrals in
 order to identify the directionality of the signal.
 
\end_layout

\begin_layout Standard
Apart from the visual confirmation, further supporting evidence that equatorial
 integration is crucial for derivation of directionality can be seen in
 Eq.
\begin_inset space ~
\end_inset


\family roman
\series medium
\shape up
\size normal
\emph off
\bar no
\noun off
\color none

\begin_inset CommandInset ref
LatexCommand ref
reference "eq:ODF_DNI"

\end_inset

 where an equatorial integral creates a connection between the real ODF
 and the signal.
 The Funk-Radon Transform (FRT) used by 
\family default
\series default
\shape default
\size default
\emph default
\bar default
\noun default
\color inherit

\begin_inset CommandInset citation
LatexCommand cite
key "Tuch2004"

\end_inset

 is another example where equatorial integration is employed using the reconstru
ction sphere.
 We will see next that DNI and FRT are just a subset of the EIT.
\end_layout

\begin_layout Standard
With EIT the most important goal is to try to identify the orientational
 variation in the signal in the most accurate way by generating a spherical
 density.
 However, it is possible to calculate additionally the classical ODF as
 defined by Wedeen et al.
\begin_inset space ~
\end_inset


\begin_inset CommandInset citation
LatexCommand cite
key "Wedeen"

\end_inset

.
\end_layout

\begin_layout Standard
The EIT shown in Eq.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "eq:ODF_EIT"

\end_inset

 consists of an integration along the equator and along radial lines.
 A function 
\begin_inset Formula $F$
\end_inset

 of the signal is multiplied by a radial weighting function 
\begin_inset Formula $O$
\end_inset

.
 This construction is a generalization of the previous ODFs and it can support
 successfully many different function families for 
\begin_inset Formula $F$
\end_inset

 and 
\begin_inset Formula $O$
\end_inset

 which can all more or less accurately identify the directional distribution
 of the signal.
 More precisely the EIT is defined as
\end_layout

\begin_layout Standard

\lang british
\begin_inset Formula 
\begin{equation}
\psi_{EIT}(\hat{\mathbf{u}})=\int_{\hat{\mathbf{u}}^{\perp}}\int_{0}^{\infty}F(E(q))O(q)dqd\phi\label{eq:ODF_EIT}
\end{equation}

\end_inset


\end_layout

\begin_layout Standard
\align left

\lang british
where 
\begin_inset Formula $F$
\end_inset

 could be for example any of the following functions
\lang english

\begin_inset Formula 
\begin{equation}
F(E(q))=\begin{cases}
\,\,\,\,\,\,\,\,\, E(q) & (I)\\
-\nabla^{2}(E(q)) & (II)\\
\,\,\,\,\,\nabla^{4}(E(q)) & (III)\\
\,\,\,\,\,\ldots
\end{cases}\label{eq:eit_operator}
\end{equation}

\end_inset

and 
\begin_inset Formula $O$
\end_inset

 could be for example any of the following functions
\begin_inset Formula 
\begin{equation}
O(q)=\begin{cases}
1 & (0)\\
q & (1)\\
q^{2} & (2)\\
\ldots
\end{cases}\label{eq:eit_radial_weighting}
\end{equation}

\end_inset


\begin_inset Float table
wide false
sideways false
status open

\begin_layout Plain Layout
\align center
\begin_inset Tabular
<lyxtabular version="3" rows="5" columns="4">
<features tabularvalignment="middle">
<column alignment="center" valignment="top" width="0">
<column alignment="center" valignment="top" width="0">
<column alignment="center" valignment="top" width="0">
<column alignment="center" valignment="top" width="0">
<row>
<cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
\begin_inset Formula $F$
\end_inset


\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
\begin_inset Formula $O$
\end_inset


\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
Name
\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" rightline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
Comment
\end_layout

\end_inset
</cell>
</row>
<row>
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout

\family roman
\series medium
\shape up
\size normal
\emph off
\bar no
\noun off
\color none
\begin_inset Formula $-\nabla^{2}(E(q))$
\end_inset


\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
\begin_inset Formula $q$
\end_inset


\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
DNI
\begin_inset Formula $\equiv$
\end_inset

EITL
\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" rightline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
\begin_inset Tabular
<lyxtabular version="3" rows="2" columns="1">
<features tabularvalignment="middle">
<column alignment="center" valignment="top" width="0">
<row>
<cell alignment="center" valignment="top" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
calculates the real ODF without
\end_layout

\end_inset
</cell>
</row>
<row>
<cell alignment="center" valignment="top" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
 the complications of the FFT
\end_layout

\end_inset
</cell>
</row>
</lyxtabular>

\end_inset


\end_layout

\end_inset
</cell>
</row>
<row>
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
\begin_inset Formula $\nabla^{4}(E(q))$
\end_inset


\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
\begin_inset Formula $q$
\end_inset


\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
EITL2
\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" rightline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
high resolution at low angles
\end_layout

\end_inset
</cell>
</row>
<row>
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout

\family roman
\series medium
\shape up
\size normal
\emph off
\bar no
\noun off
\color none
\begin_inset Formula $E(q)$
\end_inset


\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
\begin_inset Formula $q$
\end_inset


\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
EITS
\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" rightline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
\begin_inset Tabular
<lyxtabular version="3" rows="2" columns="1">
<features tabularvalignment="middle">
<column alignment="center" valignment="top" width="0">
<row>
<cell alignment="center" valignment="top" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
impressive resolution without
\end_layout

\end_inset
</cell>
</row>
<row>
<cell alignment="center" valignment="top" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
any preprocessing of the signal
\end_layout

\end_inset
</cell>
</row>
</lyxtabular>

\end_inset

 
\end_layout

\end_inset
</cell>
</row>
<row>
<cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout

\family roman
\series medium
\shape up
\size normal
\emph off
\bar no
\noun off
\color none
\begin_inset Formula $E(q)$
\end_inset


\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
\begin_inset Formula $1$
\end_inset


\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
'QBI'-like
\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" rightline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
similar to the Funk Radon Transform
\end_layout

\end_inset
</cell>
</row>
</lyxtabular>

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption

\begin_layout Plain Layout
The Equatorial inversion transform (EIT) can be used to explain many other
 reconstructions algorithms.
 
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset CommandInset label
LatexCommand label
name "Tab:EIT_table"

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
In Tab.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Tab:EIT_table"

\end_inset

 we see that by choosing different functions for 
\begin_inset Formula $F$
\end_inset

 or 
\begin_inset Formula $O$
\end_inset

 we can generate both old and new distribution functions on the sphere.
 With 
\begin_inset Formula $F(E(q))=-\nabla^{2}(E(q))$
\end_inset

 and 
\begin_inset Formula $O(q)=q$
\end_inset

 we can generate 
\begin_inset Formula $\psi_{DNI}$
\end_inset

 which is theoretically identical to the DSI real ODF(
\begin_inset Formula $\psi_{DSI}$
\end_inset

).
 If 
\begin_inset Formula $F(E(q))=E(q)$
\end_inset

 and 
\begin_inset Formula $O(q)=1$
\end_inset

 then this is similar to the Funk Radon Transform (used in Q-ball imaging)
 but applied to multiple spherical shells.
 However, we can also try to use different functions like 
\begin_inset Formula $F(E(q))=-\nabla^{4}(E(q))$
\end_inset

 and 
\begin_inset Formula $O(q)=q$
\end_inset

 which can potentially increase the amount of directional information beyond
 that of the standard ODFs.
 Before starting investigating the realms of EIT we will first give a short
 overview of other methods commonly found in the literature.
\end_layout

\begin_layout Subsection
Implementation 
\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename last_figures/fast_DNI.png
	scale 50

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption

\begin_layout Plain Layout
A: Standard EIT.
 B: Fast EIT.
 Fast EIT is an order of magnitude faster than standard EIT.
 The key idea here is that we can reduce computations by storing the sum
 of the radial integrals for every vertex in the reconstruction sphere and
 then we can also precompute the indices of the vertices that are near the
 equator of every vertex (inside an equatorial zone).
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset CommandInset label
LatexCommand label
name "Flo:sEITvsEIT"

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Subsubsection
Standard EIT
\begin_inset CommandInset label
LatexCommand label
name "sub:Standard-EIT"

\end_inset


\end_layout

\begin_layout Standard
Eq.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "eq:ODF_DNI"

\end_inset

 and 
\begin_inset CommandInset ref
LatexCommand ref
reference "eq:ODF_EIT"

\end_inset

 can be implemented in a standard way by evaluating the 3D signal on the
 grid multiple times for every direction 
\begin_inset Formula $\hat{\bm{u}}$
\end_inset

 as shown in Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:sEITvsEIT"

\end_inset

A.
 This suggests that if for example we use a reconstruction sphere of 
\begin_inset Formula $642$
\end_inset

 vertices and the signal is centered inside a cubic grid of size 
\begin_inset Formula $17\times17\times17$
\end_inset

 where the radial integration (
\begin_inset Formula $q$
\end_inset

) takes place in 
\begin_inset Formula $30$
\end_inset

 steps and the equatorial (
\begin_inset Formula $\phi$
\end_inset

) in 
\begin_inset Formula $63$
\end_inset

 steps, then we need to interpolate 
\begin_inset Formula $642\times30\times30\simeq1.2$
\end_inset

 million times on the cubic grid.
 For this reason we invented Fast EIT, a new method that needs an order
 of magnitude less evaluations.
\end_layout

\begin_layout Standard
In this document whenever we use the prefix 's' in front of a method it
 means it was calculated with the standard EIT algorithm.
 For example if standard EIT is used for DNI we will write sDNI or sEITL.
 ('L' stands for 'Laplacian' or 'Nabla').
 Of course, sDNI and sEITL are equivalent.
\end_layout

\begin_layout Subsubsection
Fast EIT
\end_layout

\begin_layout Standard
A much faster algorithm than the standard EIT is described here.
 The main idea is that we can store the sum of the radial integrals for
 every vertex in the reconstruction sphere so we can then precompute the
 indices of the vertices that are near the equator of every vertex (inside
 an equatorial zone).
 This becomes clear in Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:sEITvsEIT"

\end_inset

B.
 Following these calculations, the spherical distribution function can be
 approximated with much less operations.
 The full algorithm is given in Alg.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Alg:Fast_EIT"

\end_inset

.
 The input is the vertices 
\begin_inset Formula $\hat{\bm{u}}_{i}$
\end_inset

 of the reconstruction sphere and the normalized signal 
\begin_inset Formula $E.$
\end_inset

 Then, for every point of the reconstruction sphere 
\begin_inset Formula $\hat{\bm{u}}_{i}$
\end_inset

, we save the indices of the vertices 
\begin_inset Formula $j$
\end_inset

 of 
\begin_inset Formula $\hat{\bm{u}}_{j}$
\end_inset

, which are inside an equatorial zone, in list 
\begin_inset Formula $J_{i}$
\end_inset

.
 The width of the equatorial zone 
\begin_inset Formula $z$
\end_inset

 is a constant set empirically to 
\begin_inset Formula $5^{\circ}$
\end_inset

 .
 If a very highly dense reconstruction sphere is used with more than 
\begin_inset Formula $642$
\end_inset

 vertices, which is the one we used, then the zone can be smaller.
 That can potentially increase the angular resolution of the method.
\end_layout

\begin_layout Standard
\begin_inset Float algorithm
wide false
sideways false
status open

\begin_layout Plain Layout
\begin_inset ERT
status open

\begin_layout Plain Layout


\backslash
textbf{Input} $U=
\backslash
{
\backslash
hat{
\backslash
bm{u}}_{1}
\backslash
ldots
\backslash
hat{
\backslash
bm{u}}_{m}
\backslash
}$, E
\backslash

\backslash

\end_layout

\begin_layout Plain Layout


\backslash
textbf{Output} $
\backslash
bm{
\backslash
psi}_{EIT}$
\backslash

\backslash

\end_layout

\begin_layout Plain Layout


\backslash

\backslash

\end_layout

\begin_layout Plain Layout


\backslash
textbf{For} $
\backslash
hat{
\backslash
bm{u}}_{i}$ in $U$ 
\backslash
textbf{Do}
\backslash

\backslash

\end_layout

\begin_layout Plain Layout


\backslash
hspace*{2em} $J_{i}=
\backslash
{j : |
\backslash
arccos(
\backslash
hat{
\backslash
bm{u}}_{i}
\backslash
cdot
\backslash
hat{
\backslash
bm{u}}_{j})|
\backslash
leq z 
\backslash
}$
\backslash

\backslash

\end_layout

\begin_layout Plain Layout


\backslash
hspace*{2em} $
\backslash
mathbb{B}(
\backslash
hat{
\backslash
bm{u}}_{i})=
\backslash
sum_{k=0}^{n}F(E(q_{k}
\backslash
hat{
\backslash
bm{u}}_{i}))O(q_{k}
\backslash
hat{
\backslash
bm{u}}_{i})$
\backslash

\backslash

\end_layout

\begin_layout Plain Layout


\backslash
hspace*{4em} where $F(E(q_{k}
\backslash
hat{
\backslash
bm{u}}_{i}))$ is interpolated on the lattice.
\backslash

\backslash

\end_layout

\begin_layout Plain Layout


\backslash
textbf{EndFor} 
\backslash

\backslash

\end_layout

\begin_layout Plain Layout


\backslash
textbf{For} $
\backslash
hat{
\backslash
bm{u}}_{i}$ in $U$ 
\backslash
textbf{Do}
\backslash

\backslash

\end_layout

\begin_layout Plain Layout


\backslash
hspace*{2em} $
\backslash
psi_{EIT}(
\backslash
hat{
\backslash
bm{u}}_{i})=
\backslash
frac{1}{N_i}
\backslash
sum_{j
\backslash
in J_i}
\backslash
mathbb{B}(
\backslash
hat{
\backslash
bm{u}}_{j})$
\backslash

\backslash
 
\end_layout

\begin_layout Plain Layout


\backslash
hspace*{4em} where $N_i$ is the number of indices in $J_i$.
\backslash

\backslash

\end_layout

\begin_layout Plain Layout


\backslash
textbf{EndFor}
\end_layout

\begin_layout Plain Layout

\end_layout

\begin_layout Plain Layout


\backslash
caption{Fast Equatorial Inversion Transform}
\end_layout

\begin_layout Plain Layout

\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset CommandInset label
LatexCommand label
name "Alg:Fast_EIT"

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
At the next stage we calculate sums along every radius on the direction
 of 
\begin_inset Formula $\hat{\bm{u}}_{i}$
\end_inset

 in the following way: 
\begin_inset ERT
status open

\begin_layout Plain Layout

$
\backslash
mathbb{B}(
\backslash
hat{
\backslash
bm{u}}_{i})=
\backslash
sum_{k=0}^{n}F(E(q_{k}
\backslash
hat{
\backslash
bm{u}}_{i}))O(q_{k}
\backslash
hat{
\backslash
bm{u}}_{i})$
\end_layout

\end_inset

 and obtain the final EIT ODF as the average of the sums in the equator
 
\begin_inset ERT
status open

\begin_layout Plain Layout

$
\backslash
psi_{EIT}(
\backslash
hat{
\backslash
bm{u}}_{i})=
\backslash
frac{1}{N_i}
\backslash
sum_{j
\backslash
in J_i}
\backslash
mathbb{B}(
\backslash
hat{
\backslash
bm{u}}_{j})$
\end_layout

\end_inset

 where 
\begin_inset Formula $F$
\end_inset

 is evaluated with trilinear for example interpolation on the lattice and
 
\begin_inset Formula $N_{i}$
\end_inset

 is the number of indices in 
\begin_inset Formula $J_{i}.$
\end_inset


\end_layout

\begin_layout Standard
In section 
\begin_inset CommandInset ref
LatexCommand ref
reference "sub:Multi-fiber-Simulations"

\end_inset

 the standard EITL (sEITL) is compared with fast EITL.
 In Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:FastvsStandardEITL"

\end_inset

 it is shown that the fast EIT has very similar results with the standard
 EIT.
 From now on whenever we write EIT we refer to the fast version.
 
\end_layout

\begin_layout Subsection
Peak Finding
\begin_inset CommandInset label
LatexCommand label
name "sub:Peak-Finding"

\end_inset


\end_layout

\begin_layout Standard
After we have generated the ODFs we need to find the peaks (local maxima)
 from which we can easily approximate the direction of the fibres.
 Peak finding is not trivial if there are many local maxima in the ODFs
 or the ODFs are noisy.
 Here we present an algorithm (see Alg.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Alg:PeakFinding"

\end_inset

) which reduces the amount of small local variations and returns a number
 of sorted peaks and their indices in the reconstruction sphere.
 The input of this algorithm is 
\begin_inset Formula $\bm{\psi}$
\end_inset

 (ODF) and the faces of a symmetric on the z-axis evenly distributed sphere
 (see Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:sphere_642"

\end_inset

C).
\end_layout

\begin_layout Standard
\begin_inset Float algorithm
wide false
sideways false
status open

\begin_layout Plain Layout
\begin_inset ERT
status open

\begin_layout Plain Layout


\backslash
textbf{Input} ODF $
\backslash
bm
\backslash
psi$, faces $
\backslash
Phi$ 
\backslash

\backslash

\end_layout

\begin_layout Plain Layout


\backslash
textbf{Output} peaks $P$ and indices $I$ 
\backslash

\backslash

\end_layout

\begin_layout Plain Layout


\backslash

\backslash

\end_layout

\begin_layout Plain Layout


\backslash
textbf{For} face $
\backslash
Phi_i$ in $
\backslash
Phi$ 
\backslash
textbf{Do}
\backslash

\backslash

\end_layout

\begin_layout Plain Layout


\backslash
hspace*{2em} $f_0,f_1,f_2 = 
\backslash
Phi_i$ 
\backslash

\backslash
	
\end_layout

\begin_layout Plain Layout


\backslash
hspace*{2em} $d_0,d_1,d_2 =
\backslash
bm{
\backslash
psi}[f_0],
\backslash
bm{
\backslash
psi}[f_1],
\backslash
bm{
\backslash
psi}[f_2]$ 
\backslash

\backslash
	
\end_layout

\begin_layout Plain Layout


\backslash
hspace*{2em} 
\backslash
textbf{If} $d_0 
\backslash
geq d_1$ 
\backslash
textbf{and} $d_2$ 
\backslash
textbf{Do}
\backslash

\backslash

\end_layout

\begin_layout Plain Layout


\backslash
hspace*{4em} $P[f_1]=P[f_2]=0$ 
\backslash

\backslash

\end_layout

\begin_layout Plain Layout


\backslash
hspace*{4em} 
\backslash
textbf{continue} 
\backslash

\backslash

\end_layout

\begin_layout Plain Layout


\backslash
hspace*{2em} 
\backslash
textbf{If} $d_1 
\backslash
geq d_0$ 
\backslash
textbf{and} $d_2$ 
\backslash
textbf{Do}
\backslash

\backslash

\end_layout

\begin_layout Plain Layout


\backslash
hspace*{4em} $P[f_0]=P[f_2]=0$ 
\backslash

\backslash

\end_layout

\begin_layout Plain Layout


\backslash
hspace*{4em} 
\backslash
textbf{continue} 
\backslash

\backslash

\end_layout

\begin_layout Plain Layout


\backslash
hspace*{2em} 
\backslash
textbf{If} $d_2 
\backslash
geq d_0$ 
\backslash
textbf{and} $d_1$ 
\backslash
textbf{Do}
\backslash

\backslash

\end_layout

\begin_layout Plain Layout


\backslash
hspace*{4em} $P[f_0]=P[f_1]=0$ 
\backslash

\backslash

\end_layout

\begin_layout Plain Layout


\backslash
hspace*{4em} 
\backslash
textbf{continue} 
\backslash

\backslash

\end_layout

\begin_layout Plain Layout


\backslash
hspace*{2em} 
\backslash
textbf{EndIf} 
\backslash

\backslash

\end_layout

\begin_layout Plain Layout


\backslash
textbf{EndFor}
\end_layout

\begin_layout Plain Layout

\end_layout

\begin_layout Plain Layout


\backslash
caption{Peak Finding with a Symmetric Ordered Sphere}
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset CommandInset label
LatexCommand label
name "Alg:PeakFinding"

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
We have used a triangulation of the unit sphere (which we refer to simply
 as 'sphere') obtained by triangular subdivision of a regular icosahedron.
 It is symmetric over the z-axis, i.e.
\begin_inset space ~
\end_inset

for each vertex 
\begin_inset Formula $(x,\, y,\, z)$
\end_inset

 there is a corresponding vertex 
\begin_inset Formula $(x,\, y,\,-z)$
\end_inset

 in the opposite hemisphere.
 The same sphere was used in 
\begin_inset CommandInset citation
LatexCommand cite
key "Yeh2010"

\end_inset

 for GQI reconstructions.
 Every face (triangle) corresponds to a list of the 
\begin_inset Formula $3$
\end_inset

 indices of the 
\begin_inset Formula $3$
\end_inset

 vertices on the sphere.
 The idea here is that we can travel from face to face and nullify all points
 on a face which are lower that the higher value of the face.
 At the end only local maxima will survive the procedure.
 The algorithm is presented in detail in Alg.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Alg:PeakFinding"

\end_inset

.
\end_layout

\begin_layout Standard
The sphere we use is of course discrete and therefore it adds some constraints
 on the angular resolution (worst case
\begin_inset Formula $\pm4.96^{\circ}$
\end_inset


\begin_inset Note Note
status collapsed

\begin_layout Plain Layout
\begin_inset Formula $f+v-e=2$
\end_inset

, 
\begin_inset Formula $e=3f/2$
\end_inset

, assume equilateral triangle with area 
\begin_inset Formula $A=\frac{4\pi}{1280}=\frac{\sqrt{3}}{4}l^{2}$
\end_inset

 then 
\begin_inset Formula $l=0.15$
\end_inset

 now the circumradius is given by 
\begin_inset Formula $R=l/\sqrt{3}=.008$
\end_inset

 in radians and in degrees is 
\begin_inset Formula $4.96^{\circ}.$
\end_inset


\end_layout

\end_inset

) of the peaks found from the ODF.
 In addition, the proposed Peak Finding algorithm can reduce slightly more
 the angular resolution.
 For example, in Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:sphere_642"

\end_inset

A, B we show that if point 
\series bold
a
\series default
 was a local maxima then only points 
\series bold
b
\series default
 could be alternative local maxima for 
\series bold
a
\series default
 but none of the unlabeled points could be a second peak.
 Nevertheless, we found Alg.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Alg:PeakFinding"

\end_inset

 to be extremely useful and fast.
 The same algorithm was used also in 
\begin_inset CommandInset citation
LatexCommand cite
key "Yeh2010"

\end_inset

 but it was not documented as such.
\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename last_figures/peak_finding_sphere.png
	lyxscale 30
	scale 60
	rotateOrigin center

\end_inset


\end_layout

\begin_layout Plain Layout
\align center
\begin_inset Caption

\begin_layout Plain Layout
A and B: Each point 
\series bold
a
\series default
 is a local maximum for all its neighboring faces, then only at 
\series bold
b
\series default
 are other possible local maxima.
 This simple illustration shows that the triangulation of the sphere is
 important for the determination of closed peaks and that peaks which belong
 to the same triangle cannot be determined.
 C: the sphere used for ODF reconstructions consisting of 
\begin_inset Formula $642$
\end_inset

 vertices and 
\begin_inset Formula $1,280$
\end_inset

 faces produced by subdivisions of the icosahedron.
\end_layout

\end_inset


\begin_inset CommandInset label
LatexCommand label
name "Flo:sphere_642"

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Subsection
Spherical Angular Smoothing
\begin_inset CommandInset label
LatexCommand label
name "sub:Spherical-Angular-Smoothing"

\end_inset


\end_layout

\begin_layout Standard
All current non-parametric dMRI reconstruction algorithms use some type
 of 
\begin_inset Quotes eld
\end_inset

smoothing
\begin_inset Quotes erd
\end_inset

 to reduce the effect of noise in the real data sets.
 DSI uses a Hanning filter and then avoids sampling from low values in 
\series bold
r
\series default
-space.
 In GQI, smoothing is controlled from a scalar parameter; the diffusion
 sampling length and in spherical harmonic inversion methods 
\begin_inset CommandInset citation
LatexCommand cite
key "Descoteaux2007MagResMed"

\end_inset

, 
\begin_inset CommandInset citation
LatexCommand cite
key "aganj2010reconstruction"

\end_inset

 the amount of smoothing is controlled by using only a number of the first
 components of the SH series.
 
\end_layout

\begin_layout Standard
All these approaches smooth and calculate the ODFs simultaneously.
 Our appoach differs in that we propose that the ODF is first calculated
 and then smoothed.
 For example, using the operator shown below in matrix form
\begin_inset Formula 
\[
W=\exp(\frac{U\cdot U^{T}}{\sigma})
\]

\end_inset


\end_layout

\begin_layout Standard
\noindent
\align block
where 
\begin_inset Formula $U$
\end_inset

 is the an 
\begin_inset Formula $N\times3$
\end_inset

 matrix holding the 
\begin_inset Formula $N$
\end_inset

 points of the ODF reconstruction sphere and 
\begin_inset Formula $\sigma$
\end_inset

 is a smoothing parameter acting like the variance.
 At the next step we can smooth any ODF(
\begin_inset Formula $\bm{\psi}$
\end_inset

) creating a new ODF(
\begin_inset Formula $\bm{\psi}'$
\end_inset

) in the following way
\family roman
\series medium
\shape up
\size normal
\emph off
\bar no
\noun off
\color none

\begin_inset Formula 
\begin{equation}
\bm{\psi}'=\bm{\psi}\cdot\frac{W}{\sum_{j}W_{j}}\label{eq:spherical_gaussian_angular_smoothing}
\end{equation}

\end_inset


\end_layout

\begin_layout Standard
\noindent
\align block

\family roman
\series medium
\shape up
\size normal
\emph off
\bar no
\noun off
\color none
where 
\begin_inset Formula $j$
\end_inset

 denotes row indexing, 
\begin_inset Formula $\sum_{j}W_{j}$
\end_inset

 acts as a normalization for the angular weighting 
\begin_inset Formula $W$
\end_inset

, 
\begin_inset Formula $\bm{\psi}$
\end_inset

 is the initial ODF and
\family default
\series default
\shape default
\size default
\emph default
\bar default
\noun default
\color inherit
 
\family roman
\series medium
\shape up
\size normal
\emph off
\bar no
\noun off
\color none

\begin_inset Formula $\bm{\psi}'$
\end_inset

 is the smoothed ODF 
\family default
\series default
\shape default
\size default
\emph default
\bar default
\noun default
\color inherit
.
 The advantage of this method is that it is more comprehensive and direct.
 It also uses information from all directions simultaneously.
 Similar operators can be constructed that weight more lower or higher peaks.
 The operator shown here weighs more peaks that are closer in angular distance.
 In Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:Smoothing_Example"

\end_inset

 we see the effect of this equation on a simulated triple-fibre crossing;
 distorted with Gaussian noise with SNR 
\begin_inset Formula $20$
\end_inset

 and reconstructed as an EITL density function.
 The simulation was used using a Sticks and Ball model with diffusivity
 value 
\begin_inset Formula $0.0015$
\end_inset

 and S0=
\begin_inset Formula $100$
\end_inset

.
\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename last_figures/triple_crossing_laplacian_snr20.png
	lyxscale 50
	scale 20

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption

\begin_layout Plain Layout
An example of spherical angular Gaussian smoothing applied with different
 smoothing factors on the ODF of a triple-fibre crossing on the left.
 
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\align center
\begin_inset CommandInset label
LatexCommand label
name "Flo:Smoothing_Example"

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
We can easily observe that when we increase the smoothing factor 
\begin_inset Formula $\sigma$
\end_inset

, small noisy peaks, as seen in the center of the unsmoothed spherical function,
 can easily be removed.
 However, with too much smoothing even the strongest peaks can lose their
 definition.
 This spherical operator can help to set the trade-off between noise and
 signal and it can also simplify the peak finding process, i.e.
\begin_inset space ~
\end_inset

finding the underlying primary fibre directions as this problem is much
 easier on smooth surfaces.
\end_layout

\begin_layout Standard
Finally, decoupling the smoothing from the reconstruction step gives an
 important advantage: reducing the effect of the noise to our data more
 strongly and independently.
 Many spherical operators can be added as plugins independently of the reconstru
ction phase, and these can work with any function on the sphere (see Eq.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "eq:spherical_gaussian_angular_smoothing"

\end_inset

).
\end_layout

\begin_layout Subsection
Comparisons and Results
\end_layout

\begin_layout Standard
Validation of reconstruction and tractography algorithms is not straightforward
 due to the lack of relevant gold standards.
 Simulated voxels and software phantoms are a useful way to overcome this
 difficulty and test new methods.
 Following the simulation results, we show results with real human data
 sets.
\end_layout

\begin_layout Subsubsection
Multi-fibre Simulations
\begin_inset CommandInset label
LatexCommand label
name "sub:Multi-fiber-Simulations"

\end_inset


\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\begin_inset ERT
status open

\begin_layout Plain Layout

[th!]
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename last_figures/low_angles_with_labels.png
	lyxscale 20
	scale 60

\end_inset


\begin_inset Caption

\begin_layout Plain Layout
Visualizing ODFs created from different reconstruction methods sDNI (sEITL),
 EITL, EITL2, DSI, GQI and GQI2.
 We can see that standard DNI (sDNI), EITL and EITL2 can resolve the correct
 angular fibre directions at lower angles than the other methods.
 For example see column at angle of 
\begin_inset Formula $25^{\circ}$
\end_inset

.
 
\end_layout

\end_inset


\begin_inset CommandInset label
LatexCommand label
name "Flo:multi-fibre-sims"

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
\align block
For single voxel simulations we used the model proposed in Behrens et al.
\begin_inset space ~
\end_inset


\begin_inset CommandInset citation
LatexCommand cite
key "Behrens2007NeuroImage"

\end_inset

; the multi-compartment model also known as Sticks and Ball which simulates
 the diffusion signal as
\begin_inset Formula 
\begin{equation}
S_{i}=S_{0}((1-\sum_{j=1}^{N}f_{j})\exp(-b_{i}d)+\sum_{j=1}^{N}f_{j}\exp(-b_{i}d\cos(\theta_{ij})^{2}))\label{eq:sticks_ball_eq}
\end{equation}

\end_inset


\begin_inset ERT
status open

\begin_layout Plain Layout


\backslash
noindent
\end_layout

\end_inset

 where 
\begin_inset Formula $\theta_{ij}$
\end_inset

 is the angle between gradient direction 
\begin_inset Formula $\mathbf{\hat{g}}_{i}$
\end_inset

 and fibre (stick) unit direction 
\begin_inset Formula $\mathbf{\hat{u}}_{j}$
\end_inset

.
 The amount of representation for every fibre is given by 
\begin_inset Formula $f$
\end_inset

 and 
\begin_inset Formula $d$
\end_inset

 is the diffusivity value for the entire model.
 A Multi Tensor 
\begin_inset CommandInset citation
LatexCommand cite
key "liu2004characterizing"

\end_inset

 approach was also created for software phantoms using the formula
\family roman
\series medium
\shape up
\size normal
\emph off
\bar no
\noun off
\color none
\lang british

\begin_inset Formula 
\begin{equation}
S_{i}=S_{0}\sum_{j=1}^{N}\exp(-b\mathbf{\hat{g}}^{T}D_{j}\mathbf{\hat{g}})\label{eq:multitensor}
\end{equation}

\end_inset


\family default
\series default
\shape default
\size default
\emph default
\bar default
\noun default
\color inherit
where 
\begin_inset Formula $D_{j}$
\end_inset

 is the diffusion tensor for every fibre 
\begin_inset Formula $j$
\end_inset

.
\end_layout

\begin_layout Standard
In Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:multi-fibre-sims"

\end_inset

 we present the outcome of an experiment of two crossing fibres using different
 reconstruction methods: sDNI (sEITL), EITL, EITL2, DSI, GQI and GQI2.
 These are based on simulations of 2-fibre crossings from 
\begin_inset Formula $0^{\circ}$
\end_inset

 to 
\begin_inset Formula $90^{\circ}$
\end_inset

 using Eq.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "eq:sticks_ball_eq"

\end_inset

 with diffusivity value of 
\begin_inset Formula $1.5\times10^{-3}$
\end_inset

 
\lang british

\begin_inset Formula $\textrm{mm}^{2}/\textrm{sec}$
\end_inset


\lang english
 and 257 b-values with maximum b-value 
\begin_inset Formula $11,000$
\end_inset

.
 However, all these methods will perform accurately beyond 
\begin_inset Formula $50^{\circ}$
\end_inset

 therefore in this figure we present only the lower angles.
 We can see that sDNI, EITL and EITL2 performed better than the other methods.
 Especially EITL2 was able to resolve crossing at 
\begin_inset Formula $25^{\circ}$
\end_inset

 which is lower than the accuracy resolved in the current state of the art
 methods.
 In order to confirm this fascinating result we created a more general experimen
t with many iterations (of 2-fibre and 3-fibre crossings) which is presented
 below.
\end_layout

\begin_layout Standard
A comparison metric is needed in order to evaluate the new and old reconstructio
n methods discussed in this chapter.
 The standard procedure is to calculate the similarity between the measured
 and simulated ground truth data sets.
 We want to calculate the angular precision of the ODFs from simulations
 derived from Eq.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "eq:sticks_ball_eq"

\end_inset

.
 We define a new similarity metric called Angular Similarity (AS) which
 computes the cosine distance of the best match between the set of measured
 fibre orientations and the known set of simulated fibres.
 This metric will be used to compare 2-fibre and 3-fibre crossings.
 AS is 
\begin_inset Formula $0$
\end_inset

 when there is no match i.e.
 angular distance is 
\begin_inset Formula $0$
\end_inset

, 
\begin_inset Formula $1$
\end_inset

 when one fibre is matched (
\begin_inset Formula $0^{\circ}$
\end_inset

), 
\begin_inset Formula $2$
\end_inset

 when two fibres are matched and 
\begin_inset Formula $3$
\end_inset

 when three fibres are matched.
 In table 
\begin_inset CommandInset ref
LatexCommand ref
reference "Tab:AS_behaviour"

\end_inset

 we show a few examples of AS behaviour with simple unit vector sets.
 
\end_layout

\begin_layout Standard
\begin_inset Float table
wide false
sideways false
status open

\begin_layout Plain Layout
\align center
\begin_inset Tabular
<lyxtabular version="3" rows="5" columns="3">
<features tabularvalignment="middle">
<column alignment="center" valignment="top" width="0">
<column alignment="center" valignment="top" width="0">
<column alignment="center" valignment="top" width="0">
<row>
<cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
Known
\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
Measured
\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" rightline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
AS
\end_layout

\end_inset
</cell>
</row>
<row>
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
\begin_inset Formula $(1,0,0),(0,1,0)$
\end_inset


\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
\begin_inset Formula $(0,0,1)$
\end_inset


\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" rightline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
\begin_inset Formula $0$
\end_inset


\end_layout

\end_inset
</cell>
</row>
<row>
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
\begin_inset Formula $(1,0,0),(0,1,0)$
\end_inset


\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
\begin_inset Formula $(0,1,0)$
\end_inset


\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" rightline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
\begin_inset Formula $1$
\end_inset


\end_layout

\end_inset
</cell>
</row>
<row>
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
\begin_inset Formula $(1,0,0),(0,1,0)$
\end_inset


\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout

\family roman
\series medium
\shape up
\size normal
\emph off
\bar no
\noun off
\color none
\lang british
\begin_inset Formula $(0,\nicefrac{\sqrt{2}}{2},\nicefrac{\sqrt{2}}{2})$
\end_inset


\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" leftline="true" rightline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout

\family roman
\series medium
\shape up
\size normal
\emph off
\bar no
\noun off
\color none
\lang british
\begin_inset Formula $\nicefrac{\sqrt{2}}{2}$
\end_inset


\end_layout

\end_inset
</cell>
</row>
<row>
<cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
\begin_inset Formula $(1,0,0),(0,1,0),(0,0,1)$
\end_inset


\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
\begin_inset Formula $(1,0,0),(0,0,1)$
\end_inset


\end_layout

\end_inset
</cell>
<cell alignment="center" valignment="top" topline="true" bottomline="true" leftline="true" rightline="true" usebox="none">
\begin_inset Text

\begin_layout Plain Layout
\begin_inset Formula $2$
\end_inset


\end_layout

\end_inset
</cell>
</row>
</lyxtabular>

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption

\begin_layout Plain Layout
Examples of angular similarity (AS) behaviour with simple unit vector sets.
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset CommandInset label
LatexCommand label
name "Tab:AS_behaviour"

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
\align block
If our ground truth set consists of 
\begin_inset Formula $g=[(1,0,0),(0,1,0)]=[g_{0},g_{1}]$
\end_inset

 and the measured set consists of 
\begin_inset Formula $m=[(0,0,1)]$
\end_inset

 then AS=
\begin_inset Formula $0$
\end_inset

.
 If the measured set was 
\begin_inset Formula $m=[(0,\nicefrac{\sqrt{2}}{2},\nicefrac{\sqrt{2}}{2})]=[m_{0}]$
\end_inset

 then AS
\family roman
\series medium
\shape up
\size normal
\emph off
\bar no
\noun off
\color none
\lang british
 is
\begin_inset Formula $\nicefrac{\sqrt{2}}{2}$
\end_inset

.
 This is because according to the AS definition we have 
\begin_inset Formula $AS(g,m)=\max(|g_{0}\cdot m_{0}|,|g_{1},m_{0}|)$
\end_inset

.
 Which is equal to 
\begin_inset Formula $\nicefrac{\sqrt{2}}{2}$
\end_inset

.
\end_layout

\begin_layout Standard
\align block

\family roman
\series medium
\shape up
\size normal
\emph off
\bar no
\noun off
\color none
\lang british
\begin_inset ERT
status open

\begin_layout Plain Layout


\backslash
noindent
\end_layout

\end_inset


\begin_inset space ~
\end_inset

If 
\begin_inset Formula $g=[(1,0,0),(0,1,0)]=[g_{0},g_{1}]$
\end_inset

 and 
\begin_inset Formula $m=g$
\end_inset

 then 
\begin_inset Formula $AS(g,m)=\max(|g_{0}\cdot m_{0}|+|g_{1}\cdot m_{1}|,|g_{0}\cdot m_{1}|+|g_{1}\cdot m_{0})=2.$
\end_inset

 We created an experiment where we set two fibres at an increasing angle
 of 
\family default
\series default
\shape default
\size default
\emph default
\bar default
\noun default
\color inherit
\lang english

\begin_inset Formula $2.5^{\circ}$
\end_inset


\family roman
\series medium
\shape up
\size normal
\emph off
\bar no
\noun off
\color none
\lang british
from 
\begin_inset Formula $0^{\circ}$
\end_inset

 to 
\begin_inset Formula $90^{\circ}$
\end_inset

 and then rotate them uniformly around 
\begin_inset Formula $200$
\end_inset

 random axes.
 This operation produces 
\begin_inset Formula $7,400$
\end_inset

 simulated ODFs and the results are shown in Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:2_100"

\end_inset

 and Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:2_20"

\end_inset

 with different signal to noise ratio.
 For these simulations noise was normally distributed.
 What we see in the figures is the average angular similarity where the
 average is calculated from the 
\begin_inset Formula $200$
\end_inset

 random orientations for the same angle.
\end_layout

\begin_layout Standard
\align block
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename last_figures/2_100.png
	lyxscale 30
	scale 50

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption

\begin_layout Plain Layout
Average angular similarity of 2-fibre crossings with SNR 100
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset CommandInset label
LatexCommand label
name "Flo:2_100"

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
\align block
We can easily observe in Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:2_100"

\end_inset

 that EITL2 can resolve more accurately fibre crossings at low angles and
 continues to perform well even at higher angles 
\begin_inset Formula $>50^{\circ}$
\end_inset

 .
 EITL performs better than DSI, GQI, GQI2 and EITS at low angles and very
 well at high angles as well.
 GQI2 performs better than DSI, GQI, and ETS.
 It is also impressive that EITS can have such a good performance although
 it is such a simple operation.
 In summary we see from the graphs that EITL2
\begin_inset Formula $>$
\end_inset

EITL
\begin_inset Formula $>$
\end_inset

GQI2
\begin_inset Formula $>$
\end_inset

DSI
\begin_inset Formula $>$
\end_inset

GQI
\begin_inset Formula $>$
\end_inset

EITS where 
\begin_inset Formula $>$
\end_inset

 means higher average angular similarity.
 The same pattern takes place even when we increase the noise level see
 for example Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:2_20"

\end_inset

.
 We will see next that the same pattern appears even with 3-fibre crossings
 and high levels of noise.
\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename last_figures/2_20.png
	lyxscale 30
	scale 50

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption

\begin_layout Plain Layout
Average angular similarity of 
\begin_inset Formula $2$
\end_inset

-fibre crossings with SNR 20.
 
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset CommandInset label
LatexCommand label
name "Flo:2_20"

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
\align block
We also measured the accuracy in 
\begin_inset Formula $3$
\end_inset

-fibre crossings.
 In this experiment the 
\begin_inset Formula $3$
\end_inset

-fibres will always have the same angular distance between each other.
 That distance will increase from 
\family roman
\series medium
\shape up
\size normal
\emph off
\bar no
\noun off
\color none
\lang british

\begin_inset Formula $0^{\circ}$
\end_inset

 to 
\begin_inset Formula $90^{\circ}$
\end_inset

 
\family default
\series default
\shape default
\size default
\emph default
\bar default
\noun default
\color inherit
\lang english
with steps of 
\begin_inset Formula $2.3^{\circ}$
\end_inset

 on average and all 
\begin_inset Formula $3$
\end_inset

 fibres will be reoriented 
\begin_inset Formula $200$
\end_inset

 times.
 That gave 
\begin_inset Formula $8,000$
\end_inset

 simulated crossings.
 
\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename last_figures/3_100.png
	lyxscale 30
	scale 50

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption

\begin_layout Plain Layout
Average angular similarity of 
\begin_inset Formula $3$
\end_inset

-fibre crossings with SNR 
\begin_inset Formula $100$
\end_inset

.
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset CommandInset label
LatexCommand label
name "Flo:3_100"

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
\align block
The results of the 3-fibre crossings shown in Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:3_100"

\end_inset

 and Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:3_20"

\end_inset

 were very similar to those of the 2-fibre crossings; EITL2 performed better
 at low angles, showing slightly reduced performance at high angles.
 EITL performed better with low angles than the rest of the methods, having
 also high accuracy on larger angles.
\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename last_figures/3_20.png
	lyxscale 30
	scale 50

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption

\begin_layout Plain Layout
Average angular similarity of 
\begin_inset Formula $3$
\end_inset

-fibre crossings with SNR 
\begin_inset Formula $20$
\end_inset

.
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset CommandInset label
LatexCommand label
name "Flo:3_20"

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
\align block
These summarization plots give strong evidence that both DNI (EITL) and
 in general EIT can be used to accurately generate spherical distribution
 functions for the determination of the directional information of the diffusion
 signal.
 They performed equivalently or better than the current state-of-the-art
 grid-based reconstruction methods i.e DSI and GQI.
 The determination of the fibre directions was not affected considerably
 by noise.
 
\end_layout

\begin_layout Standard
Furthermore, we can also see that GQI2 can do better than DSI, GQI and that
 EITS gives results that are very similar to GQI.
 The parameters used for these simulations were DSI: radial sampling 
\begin_inset Formula $2.1-6$
\end_inset

, hanning filter width: 
\begin_inset Formula $36$
\end_inset

, GQI: 
\begin_inset Formula $\text{\lambda=1.2}$
\end_inset

, GQI2: 
\begin_inset Formula $\lambda=3$
\end_inset

, and EITS, EITL, EITL2 were all calculated with the standard options zonal
 width (
\begin_inset Formula $z=5^{\circ}$
\end_inset

), grid size 
\begin_inset Formula $17\times17\times17$
\end_inset

, radial sampling 
\begin_inset Formula $0-5$
\end_inset

 with 
\begin_inset Formula $0.1$
\end_inset

 steps and no further post-processing or smoothing was used.
 All methods were using the same reconstruction sphere with 
\begin_inset Formula $642$
\end_inset

 vertices and 
\begin_inset Formula $1,280$
\end_inset

 faces.
\end_layout

\begin_layout Standard
In these tests, EIT and fast EIT produced very similar results.
 For example, a simple test for the 
\begin_inset Formula $3$
\end_inset

-fibre case as seen in Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:FastvsStandardEITL"

\end_inset

 shows that there is close agreement between the two methods i.e.
\begin_inset space ~
\end_inset

their results are nearly equivalent.
 We can therefore conclude that the fast EIT is an acceptable approximation
 of the standard EIT.
\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename last_figures/EITL_sEITL.png
	lyxscale 30
	scale 50

\end_inset


\begin_inset Caption

\begin_layout Plain Layout
This diagram shows that when we compute EITL with the fast or standard method
 the results are nearly equivalent.
 The mean angular similarity for the case of 3-fibres crossings is very
 similar when using standard DNI or fast DNI (EITL).
 
\end_layout

\end_inset


\begin_inset CommandInset label
LatexCommand label
name "Flo:FastvsStandardEITL"

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Subsubsection
Software Phantoms
\begin_inset CommandInset label
LatexCommand label
name "sub:Digital-Phantoms"

\end_inset


\end_layout

\begin_layout Standard
A software phantom generation tool was developed which can simulate the
 diffusion weighted signal for one or more fibres represented by different
 discrete 3D orbital functions.
 This work is an extension of the phantom developed by Correia et al.
\begin_inset space ~
\end_inset


\begin_inset CommandInset citation
LatexCommand cite
key "Correia2009"

\end_inset

 which supported only paths with analytically calculated derivatives.
\end_layout

\begin_layout Standard
The idea here is that we first create any orbital function 
\begin_inset Formula $f(t):\mathbb{R}\rightarrow\mathbb{R}^{3}$
\end_inset

 and calculate numerically its derivatives at small steps 
\begin_inset Formula $\Delta t$
\end_inset

.
 We can then scale it and centre it so that it fits in an image volume of
 the desired size.
 We expect that many segments of the discrete function 
\begin_inset Formula $f$
\end_inset

 will fall into every voxel in the volume and that more curved parts of
 
\begin_inset Formula $f$
\end_inset

 will have higher representation in the voxel than less curved parts.
 For every segment, we can find the main direction of the orbit 
\begin_inset Formula $\mathbf{v}=\frac{f(t+1)-f(t)}{\Delta t}$
\end_inset

 and calculate the rotation matrix 
\begin_inset Formula $\mathbf{R}$
\end_inset

 that rotates 
\begin_inset Formula $\hat{\mathbf{x}}=(1,0,0)$
\end_inset

 to 
\begin_inset Formula $\mathbf{v}$
\end_inset

.
 Then, the signal for each element of the fibre for a given b-value 
\begin_inset Formula $b$
\end_inset

 and a given gradient sampling direction 
\begin_inset Formula $\hat{\mathbf{g}}$
\end_inset

, is given by the following Single Tensor formula
\end_layout

\begin_layout Standard
\begin_inset Formula 
\begin{equation}
\Delta S=S_{0}\exp(-b\hat{\mathbf{g}}^{T}\mathbf{R}\bm{\Lambda}\mathbf{R}^{T}\hat{\mathbf{g}})\label{eq:step_signal}
\end{equation}

\end_inset


\end_layout

\begin_layout Standard
\align left
where 
\family roman
\series medium
\shape up
\size normal
\emph off
\bar no
\noun off
\color none
\lang british

\begin_inset Formula $S_{0}$
\end_inset

 is the unattenuated signal of the fibre, and the diffusion tensor is given
 by
\end_layout

\begin_layout Standard
\align center
\begin_inset Formula 
\begin{equation}
\bm{\Lambda}=\left(\begin{array}{ccc}
\lambda_{\parallel} & 0 & 0\\
0 & \lambda_{\perp} & 0\\
0 & 0 & \lambda_{\perp}
\end{array}\right)\label{eq:step_tensor}
\end{equation}

\end_inset


\end_layout

\begin_layout Standard
\noindent
\align block
Therefore, the total signal of the voxel for one gradient direction is given
 by the summations of all the contributions of the 
\begin_inset Formula $K$
\end_inset

 elements in the voxel
\end_layout

\begin_layout Standard
\begin_inset Formula 
\begin{equation}
S_{vox}=\sum_{i=1}^{K}\Delta S_{i}\label{eq:digital_phantom_signal}
\end{equation}

\end_inset


\end_layout

\begin_layout Standard
In addition, we can generate simulations of more than one fibre by generating
 a single volume for every orbit and then add them all together to create
 complex configurations in the final volume.
 This is acceptable, under the assumption that the diffusion is Gaussian
 in all compartments, because the diffusion signal is additive i.e.
 the signal of a crossing of two fibres is equal to the sum of the the signals
 of the individual fibres.
 In this way, we can simulate phantoms with Multi Tensor based diffusion
 signals as that described in Eq.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "eq:multitensor"

\end_inset

.
 We can increase the thickness of the fibres using a typical smoothing kernel
 or duplicate the fibres radially.
 At the end we can add different levels of noise e.g.
 Rician or Gaussian noise with a prespecified SNR.
 
\end_layout

\begin_layout Standard
The method we use to create these software phantoms offers the opportunity
 to simulate partial volume effects.
 If partial volume effects are not desired, we need to normalize by dividing
 by the number of fibre elements for each voxel.
 In Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:cool_phantoms"

\end_inset

 we can see the volume renderings of two different phantoms created with
 the method described here.
 This function is implemented in module 
\begin_inset Formula $\texttt{dipy.sims.phantom}$
\end_inset

.
\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename last_figures/digital_phantom.png
	lyxscale 40
	scale 60

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption

\begin_layout Plain Layout
Volume renderings of the unattenuated signals of two digital phantoms.
 On the left 3 fibres intersect on regular angles with Rician noise of SNR=
\begin_inset Formula $20$
\end_inset

.
 On the right a helicoidal fibre is shown clear of noise.
 For both phantoms 
\begin_inset Formula $S_{0}=100$
\end_inset

 and prolate tensors with eigenvalues 
\begin_inset Formula $\lambda_{\parallel}=1.4\cdot10^{-3}$
\end_inset

 
\lang british

\begin_inset Formula $m^{2}/sec$
\end_inset


\lang english
 and 
\begin_inset Formula $\lambda_{\perp}=.35\cdot10^{-3}m^{2}/sec$
\end_inset

 were used.
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\align center
\begin_inset CommandInset label
LatexCommand label
name "Flo:cool_phantoms"

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Subsubsection
Results with software phantoms
\end_layout

\begin_layout Standard
With the purpose of comparing and visualizing the differences between the
 reconstruction methods described in this chapter a phantom of two crossing
 bundles was created.
 The bundles are crossing at an angle of 
\begin_inset Formula $90^{\circ}$
\end_inset

.
 The phantom was generated using the method described in the previous section.
 Here we describe the basic steps: (a) We first represented the first bundle
 as a discrete straight path starting from point 
\begin_inset Formula $(-1,-1,0)$
\end_inset

 and ending at point 
\begin_inset Formula $(1,1,0)$
\end_inset

 with using 
\begin_inset Formula $1,000$
\end_inset

 time steps.
 (b) We scaled, centred and radially expanded this path so that it fits
 a volume of size 
\begin_inset Formula $64\times64\times64$
\end_inset

.
 This volume corresponds to the diffusion volume without any weighting.
 (c) We then applied the weightings for all the following volumes corresponding
 to non-zero b-values.
 (d) We replicated the same procedure for the other bundle which initially
 started as an orbit from position 
\begin_inset Formula $(-1,1,0)$
\end_inset

 and ended at position 
\begin_inset Formula $(1,-1,0)$
\end_inset

.
 (e) We added the two volumes together to create an 'x' shape (see Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:x-shape-thin-tensor"

\end_inset

,
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:x-shape-fat-tensor"

\end_inset

).
 (f) We added Rician noise with SNR=
\begin_inset Formula $5$
\end_inset

.
 As in this chapter we concentrate on Cartesian Lattice Q-space acquisitions
 we generated b-vectors and b-values using a keyhole Cartesian sampling
 grid 
\begin_inset CommandInset citation
LatexCommand cite
key "Tuch2002ThesisMIT"

\end_inset

 with 
\begin_inset Formula $515$
\end_inset

 q-vectors.
 The maximum b-value was 
\begin_inset Formula $11,538$
\end_inset

 and the minimum was 
\begin_inset Formula $0$
\end_inset

.
 Two sets of simulation experiments were performed each using a Tensor of
 different shapes.
\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\begin_inset ERT
status open

\begin_layout Plain Layout

[th!]
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename last_figures/software_phantom_comparisons_rician_05.png
	lyxscale 20
	scale 14

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption

\begin_layout Plain Layout
Results with an 'x' shape software phantom.
 Every single tensor compartment had the following eigenvalues 
\begin_inset Formula $\lambda_{\parallel}=1.4\times10^{-3}$
\end_inset

 
\lang british

\begin_inset Formula $\textrm{mm}^{2}/\textrm{sec}$
\end_inset


\lang english
 and 
\begin_inset Formula $\lambda_{\perp}=0.1\times10^{-3}\textrm{mm}^{2}/\textrm{sec}$
\end_inset

.
 Rician noise was added with SNR = 
\begin_inset Formula $5$
\end_inset

.
 GQI is very similar to EITS, GQI2 is very similar to EITL and DSI is very
 similar to EITL.
 In Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:x-shape-thin-tensor-zoomed"

\end_inset

 the regions at the centers of the phantoms are depicted in higher resolution.
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\align center
\begin_inset CommandInset label
LatexCommand label
name "Flo:x-shape-thin-tensor"

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
In the first experiment shown in Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:x-shape-thin-tensor"

\end_inset

,
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:x-shape-thin-tensor-zoomed"

\end_inset

 we used a more aniso- tropic prolate tensor for the simulation with eigenvalues
 
\begin_inset Formula $\lambda_{\parallel}=1.4\times10^{-3}$
\end_inset

 
\lang british

\begin_inset Formula $\textrm{mm}^{2}/\textrm{sec}$
\end_inset


\lang english
 and 
\begin_inset Formula $\lambda_{\perp}=0.1\times10^{-3}\textrm{mm}^{2}/\textrm{sec}$
\end_inset

.
 In the second experiment shown in Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:x-shape-fat-tensor"

\end_inset

 and Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:x-shape-fat-tensor-zoomed"

\end_inset

 we used a much less anisotropic prolate tensor with
\begin_inset Formula $\lambda_{\parallel}=1.7\times10^{-3}$
\end_inset

 
\lang british

\begin_inset Formula $\textrm{mm}^{2}/\textrm{sec}$
\end_inset


\lang english
 and 
\begin_inset Formula $\lambda_{\perp}=0.3\times10^{-3}\textrm{mm}^{2}/\textrm{sec}$
\end_inset

.
 The values of 
\begin_inset Formula $\lambda$
\end_inset

 are based in 
\begin_inset CommandInset citation
LatexCommand cite
key "Canales-Rodriguez2009"

\end_inset

.
 It is well known that noise has higher effect on less anisotropic areas.
 We can see this effect by comparing the overlapped FAs of these two two
 figures (
\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:x-shape-thin-tensor"

\end_inset

, 
\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:x-shape-fat-tensor"

\end_inset

).
 We can also see that all six methods (DSI, GQI, GQI2, EITL, EITL2, EITS)
 can resolve correctly the fibre directions by looking at their spherical
 distribution functions using a standard colour map.
 For visualization purposes all ODFs are shown in relative size as they
 have been scaled so that their maximum values correspond to 
\begin_inset Formula $1$
\end_inset

.
 
\end_layout

\begin_layout Standard
Furthermore, we can easily observe that GQI is mostly similar to EITS, GQI2
 is very similar to EITL and DSI is mostly similar to EITL.
 The fact that DSI ODFs are very similar to those of EITL ODFs is to be
 expected as the two methods create theoretically the same real ODFs.
 Remarkably, EITL can create these ODFs without using the Fourier Transform
 neither using any filter or thresholds in r-space which are necessary in
 DSI.
\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\begin_inset ERT
status open

\begin_layout Plain Layout

[ht!]
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename last_figures/software_phantom_comparisons_rician_05_zoomed.png
	lyxscale 50
	scale 150

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption

\begin_layout Plain Layout
Same as Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:x-shape-thin-tensor"

\end_inset

 showing in higher resolution the spherical distributions in the centers
 of the phantoms.
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\align center
\begin_inset CommandInset label
LatexCommand label
name "Flo:x-shape-thin-tensor-zoomed"

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\begin_inset ERT
status open

\begin_layout Plain Layout

[ht!]
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename last_figures/software_phantom_comparisons_rician_05_fat.png
	lyxscale 20
	scale 14

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption

\begin_layout Plain Layout
Showing the spherical distribution functions (DSI, GQI, GQI2, EITL, EITL2,
 EITS) of a software phantom generated by two bundles where each bundle
 contains single tensors along the direction of the phantom.
 On the crossing area, a dual Tensor effect in every voxel is observed.
 Every single Tensor compartment had the following eigenvalues 
\begin_inset Formula $\lambda_{\parallel}=1.7\times10^{-3}$
\end_inset

 
\lang british

\begin_inset Formula $\textrm{mm}^{2}/\textrm{sec}$
\end_inset


\lang english
 and 
\begin_inset Formula $\lambda_{\perp}=0.3\times10^{-3}\textrm{mm}^{2}/\textrm{sec}$
\end_inset

.
 Rician noise was added with SNR=
\begin_inset Formula $5$
\end_inset

.
 We also visualize simultaneously the FA for this slice.
 We can see that in the crossing area (gray background) the FA values drop
 considerably.
 However, the ODFs represent precisely the crossing.
 
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\align center
\begin_inset CommandInset label
LatexCommand label
name "Flo:x-shape-fat-tensor"

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:x-shape-thin-tensor"

\end_inset

, 
\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:x-shape-fat-tensor"

\end_inset

 show that all the different grid-based reconstruction methods can reconstruct
 correctly the underlying fibre directions even when noise is present.
 However, we can see that when tensors are less anisotropic the noise has
 a stronger effect in the resulting spherical distributions.
 We can also see that GQI & EITS are less sharp than DSI & EITL and these
 are less sharp than GQI2 & EITL2.
 Also DSI, GQI2, EITL, EITL2 have much lower minima than GQI and EITS.
 
\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\begin_inset ERT
status open

\begin_layout Plain Layout

[th]
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename last_figures/software_phantom_comparisons_rician_05_fat_zoomed.png
	lyxscale 50
	scale 150

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption

\begin_layout Plain Layout
A zoomed version of Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:x-shape-fat-tensor"

\end_inset

 showing the spherical distributions in the centers of the phantoms at higher
 resolution.
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\align center
\begin_inset CommandInset label
LatexCommand label
name "Flo:x-shape-fat-tensor-zoomed"

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
\align block
In the EIT-based reconstruction results shown in Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:x-shape-thin-tensor-zoomed"

\end_inset

 and Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:x-shape-fat-tensor-zoomed"

\end_inset

 we do not use any amount of smoothing as used in DSI (through hanning filter),
 GQI, GQI2 (through sampling length) and it is extraordinary that we still
 obtain such well defined distributions.
 If we want to apply some weighting/smoothing/denoising in EIT-based methods
 that is simply possible through the spherical angular smoothing approach
 described in section 
\begin_inset CommandInset ref
LatexCommand ref
reference "sub:Spherical-Angular-Smoothing"

\end_inset

.
 
\end_layout

\begin_layout Standard
The parameters used for these simulations were for DSI: radial sampling
 
\begin_inset Formula $2.1-6$
\end_inset

, hanning filter width: 
\begin_inset Formula $36$
\end_inset

 , GQI: 
\begin_inset Formula $\text{\lambda=1.2}$
\end_inset

, GQI2: 
\begin_inset Formula $\lambda=3$
\end_inset

, and EITS, EITL, EITL2 were all calculated with the standard options (
\begin_inset Formula $z=\pm5$
\end_inset

) and no further post-processing or smoothing was used.
 All methods were using the same reconstruction sphere with 
\begin_inset Formula $642$
\end_inset

 vertices and 
\begin_inset Formula $1,280$
\end_inset

 faces.
\end_layout

\begin_layout Subsubsection
Results with humans
\end_layout

\begin_layout Standard
We want to compare reconstruction methods on Cartesian grid-based acquisitions
 first with data sets which are rich on directions and commonly used for
 DSI processing.
 For this purpose we used a data set which was available online at 
\begin_inset ERT
status open

\begin_layout Plain Layout

$
\backslash
texttt{cmtk.org}$
\end_layout

\end_inset

 from the Diffusion Group at Ecole Polytechnique Fédérale de Lausanne (EPFL),
 Switzerland.
 So, this data set was obtained from a 3T scanner (TIM Trio, Siemens) with
 a 32 channels head coil.
 The field of view was 
\begin_inset Formula $210\times210\,\textrm{mm}^{2}$
\end_inset

, matrix size 
\begin_inset Formula $96\times96$
\end_inset

, and slice thickness 
\begin_inset Formula $3\textrm{\,\textrm{mm}}$
\end_inset

.
 
\begin_inset Formula $44$
\end_inset

 slices were acquired and the voxel resolution was 
\begin_inset Formula $2.2\times2.2\times3.0\,\textrm{mm}^{3}$
\end_inset

.
 A 
\begin_inset Formula $258$
\end_inset

-point half grid acquisition scheme with a maximum b-value of 
\begin_inset Formula $8011\,\textrm{s}/\textrm{mm}^{2}$
\end_inset

 also known as DSI515
\begin_inset space ~
\end_inset


\begin_inset CommandInset citation
LatexCommand cite
key "Wedeen2008"

\end_inset

 was used.
 The total acquisition time was 
\begin_inset Formula $34\,\textrm{min}$
\end_inset

 with TR=
\begin_inset Formula $8200\,\textrm{ms}$
\end_inset

 and TE=
\begin_inset Formula $165\,\textrm{ms}$
\end_inset

.
\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\begin_inset ERT
status open

\begin_layout Plain Layout

[th!]
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename last_figures/real_data_257_all_grid_methods_purple.png
	lyxscale 20
	scale 60

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption

\begin_layout Plain Layout
Showing the same slice of a human brain reconstructed with 6 different Cartesian
 grid q-space based methods.
 The ODFs are visualized on top of the FA slice.
 A clearer presentation of a region near the left upper corner (with purple
 shading) is given in Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:real_257_zoomed"

\end_inset

 for all the 
\begin_inset Formula $6$
\end_inset

 methods.
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\align center
\begin_inset CommandInset label
LatexCommand label
name "Flo:real_257"

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
The parameters used for these simulations were for DSI: radial sampling
 
\begin_inset Formula $2.1-6$
\end_inset

, Hanning filter width: 
\begin_inset Formula $36$
\end_inset

, GQI: 
\begin_inset Formula $\text{\lambda=1.2}$
\end_inset

, GQI2: 
\begin_inset Formula $\lambda=3$
\end_inset

, and for EITS, EITL, EITL2 were all calculated using the standard options
 for zonal width (
\begin_inset Formula $z=5^{o}$
\end_inset

) and spherical angular smoothing (
\begin_inset Formula $s=0.05)$
\end_inset

.
 All methods were using the same reconstruction sphere 
\begin_inset Formula $642$
\end_inset

 vertices and 
\begin_inset Formula $1,280$
\end_inset

 faces.
 The results of this experiment are shown on top of an FA slice of a healthy
 human in Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:real_257"

\end_inset

 and in higher resolution in Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:real_257_zoomed"

\end_inset

.
 It is observed that EITL, EITL2 and EITS can be used for reconstructing
 these data sets as their results appear very similar to the results given
 by DSI, GQI and GQI2.
 We can also easily see that EITL and EITL2 are relatively sharp which can
 be of an advantage for the purpose of recovering correctly the underlying
 real fibre directions.
\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\begin_inset ERT
status open

\begin_layout Plain Layout

[th!]
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename last_figures/real_data_257_all_grid_methods_zoomed.png
	lyxscale 50
	scale 150

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption

\begin_layout Plain Layout
The upper-left corners (purple shading region) of the panels of Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:real_257"

\end_inset

 are shown here at higher resolution.
 These data sets belong to a real human.
 In contrast with the results shown in simulations (see Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:x-shape-thin-tensor-zoomed"

\end_inset

) we applied spherical angular smoothing with 
\begin_inset Formula $s=0.05$
\end_inset

 for EITL, EITL2 and EITS in order to remove small noisy spikes in the distribut
ions.
 In agreement with the results of Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:x-shape-thin-tensor-zoomed"

\end_inset

, EITS is very similar to GQI.
 The difference between DSI, GQI2 and EITL, EITL2 is smaller as a result
 of the application of angular weighting.
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\align center
\begin_inset CommandInset label
LatexCommand label
name "Flo:real_257_zoomed"

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
We also tested our results with another human brain data set generated at
 a 3T scanner (TIM Trio, Siemens) at the Medical Research Council Cognition
 and Brain Sciences Unit, Cambridge, UK.
 We used Siemens advanced diffusion work-in-progress sequence, and STEAM
 
\begin_inset CommandInset citation
LatexCommand cite
key "merboldt1992diffusion,MAB04"

\end_inset

 as the diffusion preparation method.
 The field of view was 
\begin_inset Formula $240\times240\,\textrm{mm}^{2}$
\end_inset

, matrix size 
\begin_inset Formula $96\times96$
\end_inset

, and slice thickness 
\begin_inset Formula $2.5\,\textrm{mm}$
\end_inset

 (no gap).
 
\begin_inset Formula $55$
\end_inset

 slices were acquired to achieve full brain coverage, and the voxel resolution
 was 
\begin_inset Formula $2.5\times2.5\times2.5\,\textrm{mm{}^{3}}$
\end_inset

.
 In this experiment a smaller number of gradient vectors were used.
 A 
\begin_inset Formula $102$
\end_inset

-point half grid acquisition with a maximum b-value of 
\begin_inset Formula $4,000\,\textrm{s}/\textrm{mm^{2}}$
\end_inset

 was used.
 The total acquisition time was only 
\begin_inset Formula $14\,\textrm{min}\,21\,\textrm{s}$
\end_inset

 with TR=
\begin_inset Formula $8,200\,\textrm{ms}$
\end_inset

 and TE=
\begin_inset Formula $69\,\textrm{ms}$
\end_inset

.
\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\begin_inset ERT
status open

\begin_layout Plain Layout

[th!]
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename last_figures/real_101_eitl_0.05.png
	lyxscale 30
	scale 30

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption

\begin_layout Plain Layout
EITL ODFs rendered on top of FA of a human brain data set.
 A small 
\begin_inset Formula $102$
\end_inset

-point half grid acquisition with a maximum b-value of 
\begin_inset Formula $4,000\,\textrm{s}/\textrm{mm^{2}}$
\end_inset

 was used.
 Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:real_101_CBU_left_up"

\end_inset

 and 
\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:real_101_CBU_left_up_zoomed"

\end_inset

 are zoomed versions of the same figure.
 We can see clearly single fibres on the CC and CST areas but also crossing
 fibres at the Centrum Semiovale and at the areas where big bundles cross.
 Also the non-white matter areas are evidently more isotropic.
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\align center
\begin_inset CommandInset label
LatexCommand label
name "Flo:real_101_CBU"

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\begin_inset ERT
status open

\begin_layout Plain Layout

[th!]
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename last_figures/real_101_eitl_0.05_left_up_purple.png
	lyxscale 50
	scale 80

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption

\begin_layout Plain Layout
The upper part of Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:real_101_CBU"

\end_inset

 is shown here at higher resolution.
 The purple shaded part is given in higher resolution in Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:real_101_CBU_left_up_zoomed"

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\align center
\begin_inset CommandInset label
LatexCommand label
name "Flo:real_101_CBU_left_up"

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\begin_inset ERT
status open

\begin_layout Plain Layout

[th!]
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename last_figures/real_101_eitl_0.05_left_up_zoomed.png
	lyxscale 50
	scale 150

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption

\begin_layout Plain Layout
EITL ODFs of 1-fibre, 2-fibre and 3-fibre crossings from a real human data
 set of 
\begin_inset Formula $101$
\end_inset

 applied weighted diffusion volume and 
\begin_inset Formula $1$
\end_inset

 without weighting 
\begin_inset Formula $(b0)$
\end_inset

 .
 This picture is a zoomed version of the purple shaded area shown in Fig.
 
\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:real_101_CBU_left_up"

\end_inset

.
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\align center
\begin_inset CommandInset label
LatexCommand label
name "Flo:real_101_CBU_left_up_zoomed"

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
In Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:real_101_CBU"

\end_inset

 a slice is shown where different parts of white matter are visible with
 the FA background image.
 We can clearly see structures like the Corpus Callosum (CC) and Cortical-Spinal
 Tract (CST) and Centrum Semiovale areas.
 The ODFs of EITL are shown superimposed on the FA.
 The parameters used for EITL were: a standard zonal width 
\begin_inset Formula $z=5^{\circ}$
\end_inset

 and spherical angular smoothing 
\begin_inset Formula $s=0.05$
\end_inset

 with the same reconstruction sphere (
\begin_inset Formula $642$
\end_inset

 vertices, 
\begin_inset Formula $1,280$
\end_inset

 faces) as before.
\end_layout

\begin_layout Standard
For illustration purposes the upper part of Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:real_101_CBU"

\end_inset

 is depicted again in Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:real_101_CBU_left_up"

\end_inset

, and the region with purple shading from Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:real_101_CBU_left_up"

\end_inset

 is given at an even higher resolution in Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:real_101_CBU_left_up_zoomed"

\end_inset

.
 We used Mayavi 
\begin_inset CommandInset citation
LatexCommand cite
key "ramachandran2011mayavi"

\end_inset

, a Python visualization library based on VTK to make the visualizations
 shown in the figures of this section.
\end_layout

\begin_layout Standard
Although less directions were used in this acquisition scheme we obtain
 a similarly accurate depiction of the underlying white matter structure
 in comparison with that of 
\begin_inset Formula $258$
\end_inset

 directions.
 This gives great hope that we can use grid-based reconstruction methods
 with half-grid sequences with 
\begin_inset Formula $~100$
\end_inset

 gradient directions.
 This was also shown by 
\begin_inset CommandInset citation
LatexCommand cite
key "Kuo"

\end_inset

 and 
\begin_inset CommandInset citation
LatexCommand cite
key "Yeh2010"

\end_inset

 who used similar number of directions.
 
\end_layout

\begin_layout Standard
In all the figures with real data sets we can see single fibres as those
 usually found at the center of CC, and 
\begin_inset Formula $2$
\end_inset

 or 
\begin_inset Formula $3$
\end_inset

-fibre crossings in the intersection areas of CC with the CST and other
 bundles.
 
\end_layout

\begin_layout Standard
We can see for example in Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:real_101_CBU_left_up_zoomed"

\end_inset

 that the effect of spherical angular smoothing can help alleviate the noise
 effects and focus our concentration on depicting the major directions which
 are also of highest concern.
 
\end_layout

\begin_layout Subsection
Anisotropy metrics
\end_layout

\begin_layout Standard
Until this moment we discussed about density functions on the sphere as
 a way to represent complex fibre directionality in the voxel.
 These density functions are represented as multidimensional vectors containing
 
\begin_inset Formula $200$
\end_inset

 or more dimensions in each voxel and it can be cumbersome to use them directly
 for subject comparisons or visualization purposes.
 For this purpose most people use simple scalar summarizing metrics e.g.
 Tensor-based FA, MD or ODF-based such as the Generalized FA (GFA) 
\begin_inset CommandInset citation
LatexCommand cite
key "Tuch2002"

\end_inset

.
 In this section we will show that a similar scalar function like FA can
 be constructed non-parametrically.
 We call this NPA which stands for non-parametric anisotropy.
 We will also start experimenting with metrics that have more than one scalar
 value and can represent more accurately the directionality in each voxel
 that is lost with FA, GFA and MD.
 We will investigate and explain here the realms and robustness of Quantitative
 Anisotropy which was first introduced by Yeh et al.
 
\begin_inset CommandInset citation
LatexCommand cite
key "Yeh2010"

\end_inset

.
\end_layout

\begin_layout Subsubsection
Non-parametric Anisotropy
\end_layout

\begin_layout Standard
Local voxelwise measures such as fractional anisotropy (FA), apparent diffusivit
y coefficient (ADC), or mean diffusivity (MD) have been extensively adopted
 in clinical and applied research practice based on diffusion weighted MR
 imaging (dMRI).
 This underlines the need for valid and reliable measures which can indicate
 the degree of local organisation of white matter in the brain.
 The measures listed above are based on the parametric simple diffusion
 tensor (SDT or DTI) model 
\begin_inset CommandInset citation
LatexCommand cite
key "Basser1994"

\end_inset

 which works well when there is a single dominant fibre direction.
 When the local organization is more complex however, the information it
 provides is not so valid 
\begin_inset CommandInset citation
LatexCommand cite
key "Yeh2010,Tuch2002ThesisMIT"

\end_inset

.
 We show how model-free, alternatives can yield non-parametric anisotropy
 (NPA).
 These are constructed from the GQI ODF.
 We apply exact analytical results which show the form of the GQI-ODF when
 the single tensor model is correct, and further indicate how the tensor's
 parameters may be estimated from this model-free approach.
 We compare the performance of these parametric and non-parametric measures
 for simulated data.
\end_layout

\begin_layout Standard
Simulations were computed for a 
\begin_inset Formula $102$
\end_inset

-point grid sampling scheme, with a maximum b-value of 
\begin_inset Formula $4,000\,\textrm{s}/\textrm{mm}^{2}$
\end_inset

.
 The simulated fibre was aligned with the gradient frame of reference, and
 the diagonal elements of the diffusion tensor, D, where chosen to match
 typical values for white matter: 
\begin_inset Formula $\lambda_{1}=1.4\times10^{-3}$
\end_inset

 
\begin_inset Formula $\textrm{mm}^{2}/\textrm{s}$
\end_inset

, and 
\begin_inset Formula $\lambda_{2}=\lambda_{3}=0.35\times10^{-3}$
\end_inset

 
\begin_inset Formula $\textrm{mm}^{2}/\textrm{s}$
\end_inset

.
 Variable fibre orientation was realised by spatially rotating the simulated
 fibres at discrete orientations.
 
\begin_inset Formula $100$
\end_inset

 orientations were used, which spanned uniformly the space of 
\begin_inset Formula $(\theta,\phi)$
\end_inset

.
 
\end_layout

\begin_layout Standard
In addition to the SDT a two compartment model with an isotropic component
 was added with volume fraction 
\begin_inset Formula $0.5$
\end_inset

 and diffusivity 
\begin_inset Formula $0.7\times10^{-3}$
\end_inset

 
\begin_inset Formula $\textrm{mm}^{2}/\textrm{s}$
\end_inset

.
 For each acquisition scheme and fibre type, the “ideal” (noise-free) diffusion
 weighted signals were calculated according to the SDT mo- del, assuming
 a constant ideal value of the baseline signal 
\begin_inset Formula $S_{0}=100$
\end_inset

.
 Complex Gaussian noise was then superimposed upon the ideal signals to
 provide the complex noise-contaminated signals and their magnitude was
 then obtained.
 This results in noisy values with a Rician distribution, which can be scaled
 in order to set the signal to noise ratio to any desired level.
 In this study the SNRs were 
\begin_inset Formula $20$
\end_inset

, 
\begin_inset Formula $40$
\end_inset

, 
\begin_inset Formula $60$
\end_inset

, 
\begin_inset Formula $80$
\end_inset

 and 
\begin_inset Formula $100$
\end_inset

.
 The GQI ODF and SDT were fitted using 
\begin_inset ERT
status open

\begin_layout Plain Layout

$
\backslash
texttt{DIPY}$
\end_layout

\end_inset

 (
\begin_inset ERT
status open

\begin_layout Plain Layout

$
\backslash
texttt{dipy.org}$
\end_layout

\end_inset

).
 
\end_layout

\begin_layout Standard
The GQI ODF was calculated for a tessellated spherical icosahedron with
 
\begin_inset Formula $362$
\end_inset

 vertices and 
\begin_inset Formula $720$
\end_inset

 faces.
 Two values (
\begin_inset Formula $1.2$
\end_inset

 and 
\begin_inset Formula $3.5$
\end_inset

) were used for 
\begin_inset Formula $\lambda$
\end_inset

, the diffusion sampling length.
 Non-parametric FA, NPA, was calculated from the ODF by: 
\end_layout

\begin_layout Enumerate
Locating the vertex 
\begin_inset Formula $V_{1}$
\end_inset

 with maximum GQI ODF value 
\begin_inset Formula $\mathrm{max}_{1}$
\end_inset

.
\end_layout

\begin_layout Enumerate
With 
\begin_inset Formula $V_{1}$
\end_inset

 as pole, locating the vertex 
\begin_inset Formula $V_{2}$
\end_inset

 on the corresponding equatorial band of width ±5 degrees with maximum GQI
 ODF value 
\begin_inset Formula $\textrm{max}_{2}$
\end_inset

.
\end_layout

\begin_layout Enumerate
Locating a vertex 
\begin_inset Formula $V_{3}$
\end_inset

 in the equatorial band at approximately 
\begin_inset Formula $90$
\end_inset

 degrees away from 
\begin_inset Formula $V_{2}$
\end_inset

, denoting the GQI ODF value of 
\begin_inset Formula $\textrm{max}_{3}$
\end_inset

 at 
\begin_inset Formula $V_{3}$
\end_inset

.
 
\end_layout

\begin_layout Enumerate
With 
\begin_inset Formula $\textrm{npd}_{1}=\textrm{max}_{1}^{2}$
\end_inset

, 
\begin_inset Formula $\textrm{npd}_{2}=\textrm{max}_{2}^{2}$
\end_inset

, and 
\begin_inset Formula $\textrm{npd}_{3}=\textrm{max}_{3}^{2}$
\end_inset

, non-parametric anisotropy (NPA) was calculated by applying the classical
 FA formula 
\begin_inset CommandInset ref
LatexCommand ref
reference "eq:FA_classical"

\end_inset

 to the 3 values (
\begin_inset Formula $\textrm{npd}_{1},\textrm{\,\ npd}_{2},\,\textrm{npd}_{3}$
\end_inset

).
 
\end_layout

\begin_layout Standard

\series medium
The rationale for the squared ODF values is based on Tuch's formula (Eq.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "eq:tuchs_gaussian_odf"

\end_inset

) for the ODF in the SDT case which implies that the ODF in the 
\begin_inset Formula $3$
\end_inset

 principal axes directions of the tensor is proportional to the square root
 of the corresponding eigenvalue of the tensor.
 We have further derived an exact formula (see section
\begin_inset space ~
\end_inset


\begin_inset ERT
status open

\begin_layout Plain Layout


\backslash
ref{sub:Tensor-in-GQI}
\end_layout

\end_inset

): 
\begin_inset Formula $\mathrm{max_{j}}\propto\sqrt{\lambda_{j}}\,[\Phi(cL_{\Delta}/\sqrt{\lambda_{j}})-.5]$
\end_inset

 where 
\begin_inset Formula $c$
\end_inset

 is a constant that depends on the acquisition parameters, and 
\begin_inset Formula $\Phi$
\end_inset

 is the cumulative distribution function of the standard Gaussian distribution.
 
\end_layout

\begin_layout Standard

\series medium
The average NPA and FA are presented in Fig.
\begin_inset space ~
\end_inset


\series default

\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:NPA_1"

\end_inset

 
\series medium
and 
\series default

\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:NPA_2"

\end_inset


\series medium
 for 
\begin_inset Formula $200$
\end_inset

 simulations for each noise level, and single fibres with or without an
 isotro- pic component and with different diffusion sampling length.
 We can see that NPA gives very similar results with FA and as expected
 it is modulated by the degree of smoothing controlled by the value of the
 diffusion sampling length.
 
\end_layout

\begin_layout Standard

\series medium
We plan to extend this approach with voxels containing multiple peaks where
 FA would be unable to give an informative result and also extend it to
 other types of ODFs.
 In summary, we have shown that an informative new scalar anisotropy function
 (NPA) can be calculated without fitting just from the GQI ODF which promises
 to be a model-free proxy for FA.
 NPA differs from GFA 
\begin_inset CommandInset citation
LatexCommand cite
key "Tuch2004"

\end_inset

 in that it uses just 
\begin_inset Formula $3$
\end_inset

 values of the GQI ODF with a geometric relationship instead of the entire
 ODF.
 
\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\begin_inset ERT
status open

\begin_layout Plain Layout

[th!]
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename figures/NPA_L_Delta1.2.png
	scale 70

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption

\begin_layout Plain Layout
Comparison of NPA with FA for single fiber with and without an isotropic
 compartment at a range of signal to noise ratios.
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset CommandInset label
LatexCommand label
name "Flo:NPA_1"

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\begin_inset ERT
status open

\begin_layout Plain Layout

[th!]
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename figures/NPA_L_Delta3.5.png
	scale 70

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption

\begin_layout Plain Layout
As in Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:NPA_1"

\end_inset

 but with higher diffusion sampling length - less smoothing.
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset CommandInset label
LatexCommand label
name "Flo:NPA_2"

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Subsubsection
Quantitative Anisotropy
\begin_inset CommandInset label
LatexCommand label
name "sub:Quantitative-Anisotropy"

\end_inset


\end_layout

\begin_layout Standard
Quantitative anisotropy (QA) was first used by Yeh et al.
\begin_inset space ~
\end_inset


\begin_inset CommandInset citation
LatexCommand cite
key "Yeh2010"

\end_inset

 as a way of representing the peaks of the ODF with as few values as possible.
 This works in the following way: a) we create the ODF, b) we find the peaks
 using Alg.
 
\begin_inset CommandInset ref
LatexCommand ref
reference "Alg:PeakFinding"

\end_inset

, c) then 
\begin_inset Formula $QA_{i}$
\end_inset

 is equal to the peak 
\begin_inset Formula $i$
\end_inset

 minus the minimum value for the entire ODF.
 This is illustrated in Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:QA_sketch"

\end_inset

 where we can see a star-shaped ODF with three peaks (symmetric) (
\begin_inset Formula $PK$
\end_inset

).
 This ODF can be represented just with 
\begin_inset Formula $3$
\end_inset

 QA values where for example the highest value will be: 
\begin_inset Formula $QA_{0}=\max(\bm{\psi}_{GQI})-\min(\bm{\psi}_{GQI})$
\end_inset

, where 
\begin_inset Formula $\max(\bm{\psi}_{GQI})$
\end_inset

 is the value of the first peak 
\begin_inset Formula $PK_{0}$
\end_inset

, and 
\begin_inset Formula $QA_{2}=PK_{2}-\min(\bm{\psi}_{GQI})$
\end_inset

 with 
\begin_inset Formula $PK_{0}\geq PK_{1}\geq PK_{2}$
\end_inset

.
\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\begin_inset ERT
status open

\begin_layout Plain Layout

[th!]
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename QA.png
	lyxscale 50
	scale 60

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption

\begin_layout Plain Layout
QA is calculated from an ODF.
 The sphere represents the 
\begin_inset Quotes eld
\end_inset

isotropic
\begin_inset Quotes erd
\end_inset

 component (minimum value) of a GQI ODF (star) which will be removed from
 the calculation of QA.
 QA acts like a differential component with higher values in anisotropic
 areas and lower in isotropic.
 The big advantage over FA is that it can represent crossings.
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset CommandInset label
LatexCommand label
name "Flo:QA_sketch"

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
\align block
QA acts like a differential operator which is higher on anisotropic ODFs
 and lower on more isotropic.
 Actually, for a purely isotropic ODF, 
\begin_inset Formula $QA=\bm{0}$
\end_inset

.
 QA can be also easily normalized by the maximum ODF value of all voxels
 which is usually at the CSF where there is a great amount of water.
 If this normalization is in effect then we can very easily remove the backgroun
d noise i.e.
 non-white matter areas, scalp, skin, muscles etc.
 just because these will have very low QA values.
 We can see this interesting property of QA in Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:QA_directions"

\end_inset

.
 Of course the most important property of QA is that it can represent crossings
 and assign a weight for every peak.
 We will make great use of these weightings in Chapter 3.
 for the creation of tractographies.
 Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:QA_directions"

\end_inset

 was created using DSI Studio 
\begin_inset Foot
status open

\begin_layout Plain Layout
\align block
\begin_inset Formula $\texttt{dsi-studio.labsolver.org}$
\end_inset


\end_layout

\end_inset

 and the sequence parametrization is the same with the one provided in the
 experiments of the next section.
\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\begin_inset ERT
status open

\begin_layout Plain Layout

[th!]
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename figures/dsi_studio_result.png
	lyxscale 40
	scale 40

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption

\begin_layout Plain Layout
Multiple crossings of a real human data set using Quantitative Anisotropy.
 The first component of QA (
\begin_inset Formula $QA_{0}$
\end_inset

) is also shown in the background.
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset CommandInset label
LatexCommand label
name "Flo:QA_directions"

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Subsubsection
Robustness of QA
\end_layout

\begin_layout Standard
GQI was shown to have comparable accuracy to other well established q-space
 methods when it comes to resolving crossing fibres.
 In addition, this is achievable with as few as 
\begin_inset Formula $102$
\end_inset

 points on a grid sampling scheme, bringing the total acquisition time down
 to a clinically acceptable level.
 Another advantage of GQI is that it is also applicable to a shell sampling
 scheme.
 Despite their successes in tractography applications, q-space techniques
 have until now failed to produce scalar metrics that could replace the
 ones derived from the diffusion tensor model (e.g.
 mean diffusivity, MD, and fractional anisotropy, FA) in terms of their
 multi-subject comparability and specificity to pathology.
 The data acquired with a grid sampling scheme can still be used to estimate
 a diffusion tensor and respective scalar parameters, but the effects of
 the high b-values required for q-space imaging (
\begin_inset Formula $>2,000s/mm^{2}$
\end_inset

) in the accuracy of the resulting DTI parameters has not been well characterize
d.
 The authors of GQI have also proposed a new scalar metric called quantitative
 anisotropy (QA) which was described in the previous sections, but its propertie
s have not been compared to those of FA.
 In this section we will compare the estimated values of MD, FA and 
\begin_inset Formula $QA_{0}$
\end_inset

 (first component of QA) obtained with grid and shell sampling schemes,
 in terms of their precision and ability to differentiate between different
 brain fibre populations.
 Therefore, we will try to test next the hypothesis that 
\begin_inset Formula $QA_{0}$
\end_inset

 can be as used for subject comparisons.
\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\begin_inset ERT
status open

\begin_layout Plain Layout

[th!]
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename figures/MartaISMRM2.png
	lyxscale 50
	scale 130

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption

\begin_layout Plain Layout
Sample results of the paired t-tests comparing CSR (FA) and CSR (
\begin_inset Formula $QA_{0}$
\end_inset

)
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset CommandInset label
LatexCommand label
name "Flo:CSR_t-test"

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
Twelve healthy volunteers aged between 
\begin_inset Formula $18$
\end_inset

 and 
\begin_inset Formula $40$
\end_inset

 were scanned on a 3T scanner (TIM Trio, Siemens), using Siemens advanced
 diffusion work-in-progress sequence, and STEAM 
\begin_inset CommandInset citation
LatexCommand cite
key "merboldt1992diffusion,MAB04"

\end_inset

 as the diffusion preparation method.
 The field of view was 
\begin_inset Formula $240\times240\,\textrm{mm}^{2}$
\end_inset

, matrix size 
\begin_inset Formula $96\times96$
\end_inset

, and slice thickness 
\begin_inset Formula $2.5\,\textrm{mm}$
\end_inset

 (no gap).
 
\begin_inset Formula $55$
\end_inset

 slices were acquired to achieve full brain coverage, and the voxel resolution
 was 
\begin_inset Formula $2.5\times2.5\times2.5\,\textrm{mm}^{3}$
\end_inset

.
 Two sampling schemes were considered: a 102-point grid acquisition with
 a maximum b-value of 
\begin_inset Formula $4,000\,\textrm{s}/\textrm{mm}^{2}$
\end_inset

, and a single shell acquisition using 
\begin_inset Formula $118$
\end_inset

 non-collinear gradient directions and a b-value of 
\begin_inset Formula $1,000\,\textrm{s}/\textrm{mm}^{2}$
\end_inset

 
\begin_inset CommandInset citation
LatexCommand cite
key "correia_garyfallidis2011"

\end_inset

.
 The two acquisition schemes were matched for total acquisition time in
 (14 min 37
\begin_inset space ~
\end_inset

s), voxel resolution, and bandwidth.
 FA, MD and 
\begin_inset Formula $QA_{0}$
\end_inset

 maps were then generated for each acquisition scheme and for the 
\begin_inset Formula $12$
\end_inset

 volunteers using 
\begin_inset Formula $\texttt{DIPY}$
\end_inset

 
\begin_inset CommandInset citation
LatexCommand cite
key "garyfallidis2011dipy"

\end_inset

.
 All the FA datasets were non-linearly registered into MNI space using FSL
 tools, and the same transformation parameters were applied to MD and 
\begin_inset Formula $QA_{0}$
\end_inset

 maps .
 Fourteen ROIs of different brain regions were drawn in MNI space: Putamen
 (left and right), Caudate (left and right), Thalamus (left and right),
 Para-sagittal white matter (left and right), Pons, Internal Capsule (left
 and right), and Genu, Body and Splenium of the Corpus Callosum.
 Small cubic ROIs were also constructed by finding the centroid of each
 anatomical ROI and using it as the centre for a 
\begin_inset Formula $3\times3\times3\,\mathrm{mm}^{3}$
\end_inset

 ROI.
 For each ROI we calculated the mean value for each metric, and the spatial
 coefficient of variation (CV) within the ROI (see Eq.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "eq:CV_roi"

\end_inset

).
 
\begin_inset Formula 
\begin{eqnarray}
CV_{ROI} & = & \frac{\sigma_{x}}{<x>}=\frac{N_{voxels}\sqrt{\sum_{x_{i}\in ROI}(x_{i}-<x>)^{2}}}{\sqrt{N_{voxels}-1}\sum_{x_{i}\in ROI}x_{i}}\label{eq:CV_roi}
\end{eqnarray}

\end_inset


\end_layout

\begin_layout Standard
The coefficient of variation of each ROI mean across subjects was also calculate
d, as a measure of each metric’s comparability between subjects.
 The contrast-to-scatter ratio (CSR) (calculated for FA in Eq.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "eq:CSR_FA"

\end_inset

) is a good measure of a metric’s ability to differentiate between different
 brain fibre populations 
\begin_inset CommandInset citation
LatexCommand cite
key "correia_garyfallidis2011"

\end_inset

.
 
\begin_inset Formula 
\begin{eqnarray}
CSR(FA) & = & \frac{mean(FA)_{ROI_{1}}-mean(FA)_{ROI_{2}}}{\sqrt{var(FA)_{ROI_{1}}+var(FA)_{ROI_{2}}}}\label{eq:CSR_FA}
\end{eqnarray}

\end_inset


\end_layout

\begin_layout Standard
Combining the left and right versions of each ROI, we have 
\begin_inset Formula $9$
\end_inset

 ROIs of different brain populations, which can be used to define 
\begin_inset Formula $36$
\end_inset

 pairs of ROIs, and the CSR of all metrics was calculated for each of these
 pairs.
 Paired t-tests were then conducted to compare the performance of each metric
 with the two acquisition schemes, and also to compare FA and 
\begin_inset Formula $QA_{0}$
\end_inset

 directly for each acquisition scheme.
 
\end_layout

\begin_layout Standard
\begin_inset Float figure
wide false
sideways false
status open

\begin_layout Plain Layout
\begin_inset ERT
status open

\begin_layout Plain Layout

[th!]
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\align center
\begin_inset Graphics
	filename figures/MartaISMRM1.png

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset Caption

\begin_layout Plain Layout
Results of the paired t-test comparing the CVs across subjects for MD, FA
 and QA0.
\end_layout

\end_inset


\end_layout

\begin_layout Plain Layout
\begin_inset CommandInset label
LatexCommand label
name "Flo:t-test_QA0"

\end_inset


\end_layout

\end_inset


\end_layout

\begin_layout Standard
The 
\begin_inset Formula $102$
\end_inset

 grid sampling scheme produces significantly higher mean FA and 
\begin_inset Formula $QA_{0}$
\end_inset

 values than the ones obtained with the 
\begin_inset Formula $118$
\end_inset

 shell scheme, while the opposite was observed for MD.
 The CSR results for FA and QA0 were not significantly different between
 the two acquisition schemes, but the 
\begin_inset Formula $102$
\end_inset

 grid scheme produces significantly higher CSRs for MD for 26/36 ROI pairs
 (Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:CSR_t-test"

\end_inset

).
 For MD, no significant different was found for the CV across subjects,
 but for FA and QA0 the 102 scheme produced results more comparable across
 the different volunteers (Fig.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "Flo:t-test_QA0"

\end_inset

).
 For FA and MD the 
\begin_inset Formula $102$
\end_inset

 scheme showed lower CV within ROIs, especially for white matter, but no
 difference was found for 
\begin_inset Formula $QA_{0}$
\end_inset

.
 When comparing FA and 
\begin_inset Formula $QA_{0}$
\end_inset

 directly, our results show that FA produces higher CSRs than QA0 for 
\begin_inset Formula $23/36$
\end_inset

 ROI pairs for the 
\begin_inset Formula $102$
\end_inset

 grid sampling, and for 
\begin_inset Formula $19/36$
\end_inset

 ROI pairs for the 
\begin_inset Formula $118$
\end_inset

 scheme.
 FA also shows lower variation across subjects for both acquisition schemes.
 Finally, FA showes lower CVs within white matter ROIs, while 
\begin_inset Formula $QA_{0}$
\end_inset

 shows less variability for grey matter.
 The results described and shown above were obtained with the cubic ROIs,
 but do not differ significantly when the same analysis was applied to larger
 anatomical ROIs.
\end_layout

\begin_layout Standard
Our results indicate that the MD and FA maps generated from a grid sampling
 scheme designed for GQI are still suitable for analysis, since they do
 not show poorer performance when compared to a single shell and low b-value
 acquisition.
 In fact, the overall results suggest that the 
\begin_inset Formula $102$
\end_inset

 grid sampling produces slightly more robust results than the 
\begin_inset Formula $118$
\end_inset

 shell acquisition.
 A previous study 
\begin_inset CommandInset citation
LatexCommand cite
key "correia2009looking"

\end_inset

 has shown that metrics such as MD and FA benefit from the use of multiple
 b-values, which could explain the better performance of the 
\begin_inset Formula $102$
\end_inset

 grid scheme.
\end_layout

\begin_layout Subsection
Discussion and Conclusion
\end_layout

\begin_layout Standard
Non-parametric methods have the advantage of representing the signal with
 minimum number of assumptions and without needing any fitting.
 For many years there has been a trend in science to prefer model-based
 rather than model free (non-parametric) methods.
 This is perhaps because model-based can be easier to describe, and allow
 the use of popular Bayesian approaches more readily.
 However, there are some crucial issues with fitting: (a) Usually the interestin
g models have many parameters and that makes fitting very slow.
 (b) Commonly non-linear fitting is needed and accurate fitting is not trivial.
 (c) Often the model does not represent precisely the complexity of the
 real problem.
 (d) The more complex the model, the more difficult to fit 
\begin_inset CommandInset citation
LatexCommand cite
key "rice2006mathematical"

\end_inset

, 
\begin_inset CommandInset citation
LatexCommand cite
key "lee1997bayesian"

\end_inset

,
\begin_inset CommandInset citation
LatexCommand cite
key "montgomery2001introduction"

\end_inset

.
 
\end_layout

\begin_layout Standard
Non-parametric methods avoid fitting model parameters and that gives them
 a big advantage.
 The focus of this chapter was on introducing and developing new non-parametric
 methods (EIT) or comparing and extending existing ones (GQI2).
 We showed that a simple, fast and comprehensive transform exists that we
 call the Equatorial Inversion Transform (EIT).
 With this transform we showed that we can represent accurately the directional
 information of the diffusion signal.
 Furthermore, we showed that there are many different functions (
\begin_inset Formula $F$
\end_inset

 and 
\begin_inset Formula $O$
\end_inset

 see Eq.
\begin_inset space ~
\end_inset


\begin_inset CommandInset ref
LatexCommand ref
reference "eq:ODF_EIT"

\end_inset

, 
\begin_inset CommandInset ref
LatexCommand ref
reference "eq:eit_operator"

\end_inset

, 
\begin_inset CommandInset ref
LatexCommand ref
reference "eq:eit_radial_weighting"

\end_inset

) which can be used in order to create spherical density functions and use
 these to find the primary fibre directions.
 With a correct choose of 
\begin_inset Formula $F$
\end_inset

 and 
\begin_inset Formula $O$
\end_inset

 we can create theoretically the same ODF as the real ODF (DSI ODF).
 This can be done using EITL which is a type of EIT.
 Nonetheless, other density functions can be created that can identify the
 leading fibre directions without being real ODFs but they are still different
 types of spherical densities.
 EITL2 and EITS are examples of this last case.
 
\end_layout

\begin_layout Standard
The EIT concept opens new doors for the investigation of dMRI where many
 new functionals can be invented in the future that emphasise different
 properties of the signal.
 We have already illustrated and measured that EIT has the best performance
 with simulations against the state-of-the-art methods like DSI and GQI
 and that empirically, EIT gives as good results with real data sets.
\end_layout

\begin_layout Standard
The EIT finds the ODF directly without creating the diffusion propagator.
 If for some purpose the diffusion propagator is still required, then DSI
 or DPI 
\begin_inset CommandInset citation
LatexCommand cite
key "descoteaux2010multiple"

\end_inset

 are favourable.
 It could be interesting in the future to try and recover the propagator
 using ideas from the EIT.
 However, nearly always the propagator is not needed for the analysis.
 Furthermore, comparing 4D densities like the propagator is a non-trivial
 problem, also storing the propagator for every voxel is very inefficient.
\end_layout

\begin_layout Standard
We discussed that GQI can be used for creating Quantitative Anisotropy (QA).
 We observed that QA acts like a differential operator.
 It bears some similarities with FA as it is maximum on anisotropic and
 
\begin_inset Formula $0$
\end_inset

 on isotropic voxels.
 QA assumes that a substantial isotropic part can always be removed from
 the ODF and that makes it more favourable for spherical functions like
 those of GQI and EITS.
 This is in contrast to sharper densities like those of DSI and GQI2 where
 QA is not as useful because the minimum value of these densities will be
 usually near 
\begin_inset Formula $0$
\end_inset

.
 
\end_layout

\begin_layout Standard
GQI needs a manually set parameter; the diffusion sampling length and in
 contrast the EIT is fully automatic i.e.
 we always just used the few default parameters for all experiments.
 The diffusion sampling length can be slightly different from experiment
 to experiment.
 The asset of GQI and GQI2 (which was presented together with GQI but not
 investigated until today) is that they are fast to compute and have simple
 analytical solutions.
 GQI2 seems robust and smooth and it has good performance both with simulations
 and real data.
 
\end_layout

\begin_layout Standard
It is important to stress that there are similarities between all these
 methods; DSI is similar to EITL, GQI to EITS, GQI2 to EITL2.
 In addition, we showed that we can denoise the signal using a Gaussian
 Spherical Angular method which operates on spherical densities and has
 a single parameter which is similar to the variance.
\end_layout

\begin_layout Standard
Finally, we showed that the first component of QA (highest QA value) can
 be used for subject comparisons in a similar way to FA.
 We also showed that NPA could replace FA if we want to calculate anisotropy
 in a completely geometric way.
 
\end_layout

\begin_layout Standard
The source code for all the methods analyzed in this chapter is available
 at 
\begin_inset Formula $\texttt{dipy.org}$
\end_inset

 under module 
\begin_inset Formula $\texttt{dipy.reconst}$
\end_inset

.
\end_layout

\begin_layout Standard

\lang british
\begin_inset Note Note
status open

\begin_layout Plain Layout

\lang british
\begin_inset CommandInset bibtex
LatexCommand bibtex
bibfiles "diffusion"
options "ieeetr"

\end_inset


\end_layout

\end_inset


\end_layout

\end_body
\end_document