Nonnegative Diffusion Orientation Distribution Function

Abstract

Because of the well-known limitations of diffusion tensor imaging (DTI) in regions of low anisotropy and multiple fiber crossing, high angular resolution diffusion imaging (HARDI) and Q-Ball Imaging (QBI) are used to estimate the probability density function (PDF) of the average spin displacement of water molecules. In particular, QBI is used to obtain the diffusion orientation distribution function (ODF) of these multiple fiber crossing. As a probability distribution function, the orientation distribution function should be nonnegative which is not guaranteed in the existing methods. This paper proposes a novel technique to guarantee the nonnegative property of ODF by solving a convex optimization problem, which has a convex quadratic objective function and a constraint involving the nonnegativity requirement on the smallest Z-eigenvalue of the diffusivity tensor. Using convex analysis and optimization techniques, we first derive the optimality conditions of this convex optimization problem. Then, we propose a gradient descent algorithm to solve this problem. We also present formulas for determining the principal directions (maxima) of the ODF. Numerical examples on synthetic data as well as MRI data are displayed to demonstrate the significance of our approach.

Appendix: The Solution Method for (7)

According to optimization theory, the optimality conditions of (7) have the form:

$$ \left \{ \begin{array}{l} \sum_{i=1}^m \sum_{j=0}^{m-i}ib_{ij}x_1^{i-1}x_2^jx_3^{m(i,j)}=m\lambda x_1,\\[5pt] \sum_{i=0}^m \sum_{j=1}^{m-i}jb_{ij}x_1^ix_2^{j-1}x_3^{m(i,j)}=m\lambda x_2,\\[5pt] \sum_{i=0}^m \sum_{j=0}^{m-i-1}m(i,j)b_{ij}x_1^ix_2^jx_3^{m(i,j)-1}=m\lambda x_3,\\[5pt] x_1^2 + x_2^2 + x_3^2=1. \end{array} \right . $$

(18)

Here m(i,j)=m−i−j. The additional “m” on the right hand sides of the first three equations make it the same as the definition of Z-eigenvalues [7, 21–24] for the symmetric tensor x. If (x,λ) is a solution of (18), then x is a stationary point of (7) and

$$ \lambda= \varPsi(\mathbf{x}) $$

(19)

is a Z-eigenvalue of u. Then, the smallest Z-eigenvalue of u is the optimal value of (7).

We may solve (18) in the following way:

Case 1

x ₃=x ₂=0. By (18), this only happens if b _m−1,1=b _m−1,0=0. In this case, x ₁=±1, λ=b _m,0.

Case 2

x ₃=x ₁=0. By (18), this only happens if b _1,m−1=b _0,m−1=0. In this case, x ₂=±1, λ=b _0,m .

Case 3

x ₃=0, $x_{1} \not= 0$ and $x_{2} \not= 0$. Then (18) becomes

$$ \left \{ \begin{array}{l} \sum_{i=1}^m ib_{i,m-i}x_1^{i-1}x_2^{m-i}=m\lambda x_1,\\ [5pt] \sum_{i=0}^{m-1} (m-i)b_{i,m-i}x_1^ix_2^{m-i-1} =m\lambda x_2,\\[5pt] \sum_{i=0}^{m-1} b_{i,m-i-1}x_1^ix_2^{m-i-1}= 0,\\ [5pt] x_1^2 + x_2^2=1. \end{array} \right . $$

(20)

We may eliminate λ in (20) and have the following equations of x ₁ and x ₂:

$$\left \{ \begin{array}{l} \sum_{i=1}^m ib_{i,m-i}x_1^{i-1}x_2^{m-i+1}\\ [5pt] \quad=\sum_{i=0}^{m-1} (m-i)b_{i,m-i}x_1^{i+1}x_2^{m-i-1},\\ [5pt] \sum_{i=0}^{m-1} b_{i,m-i-1}x_1^ix_2^{m-i-1}=0,\\ [5pt] x_1^2 + x_2^2 =1. \end{array} \right . $$

Let t=x ₁/x ₂. We have

$$ \left \{ \begin{array}{l} \sum_{i=1}^mib_{i,m-i}t^{i-1} =\sum_{i=0}^{m-1}(m-i)b_{i,m-i}t^{i+1},\\ [5pt] \sum_{i=0}^{m-1}b_{i,m-i-1}t^i = 0. \end{array} \right . $$

(21)

We may solve the two one-variable equations of (21) separately. If they have common solutions t, then (18) has solutions

Case 4

$x_{3} \not= 0$. We may eliminate λ in (18) and have the following equations of x:

$$ \left \{ \begin{array}{l} \sum_{i=1}^m \sum_{j=0}^{m-i}ib_{ij}x_1^{i-1}x_2^jx_3^{m(i,j)+1} \\[5pt] \quad=\sum_{i=0}^m \sum_{j=0}^{m-i-1}m(i,j)b_{ij}x_1^{i+1}x_2^jx_3^{m(i,j)-1},\\ [5pt] \sum_{i=0}^m \sum_{j=1}^{m-i}jb_{ij}x_1^ix_2^{j-1}x_3^{m(i,j)+1} \\ [5pt] \quad=\sum_{i=0}^m \sum_{j=0}^{m-i-1}m(i,j)b_{ij}x_1^ix_2^{j+1}x_3^{m(i,j)-1},\\ [5pt] x_1^2 + x_2^2 + x_3^2 =1. \end{array} \right . $$

(22)

Let w=x ₁/x ₃, v=x ₂/x ₃. Then we have

$$ \left \{ \begin{array}{l} \sum_{i=1}^m \sum_{j=0}^{m-i}ib_{ij}w^{i-1}v^j \\ [5pt] \quad=\sum_{i=0}^m \sum_{j=0}^{m-i-1}m(i,j)b_{ij}w^{i+1}v^j,\\ [5pt] \sum_{i=0}^m \sum_{j=1}^{m-i}jb_{ij}w^iv^{j-1} \\[5pt] \quad=\sum_{i=0}^m \sum_{j=0}^{m-i-1} m(i,j) b_{ij}w^iv^{j+1}. \end{array} \right . $$

(23)

For solving system (23), we first regard it as a system of polynomial equations of variable w and rewrite it as

$$\left \{ \begin{array}{l} \gamma_0w^m+\gamma_1w^{m-1}+\cdots+\gamma_m=0,\\ \tau_0w^{m-1}+\tau_1w^{m-2}+\cdots +\tau_{m-1}=0, \end{array} \right . $$

where γ ₀,…,γ _m , τ ₀,…,τ _m−1 are polynomials of v, which can be calculated by (19). By the Sylvester theorem, the above system of polynomial equations in w possesses solutions if and only if its resultant vanishes [9]. The resultant of this system of polynomial equations is the determinant of the following (2m−1)×(2m−1) matrix

$$ V:= \left ( \begin{array}{c@{\ }c@{\ }c@{\ }c@{\ }c@{\ }c@{\ }c@{\ }c@{\ }c} \gamma_0&\gamma_1&\cdots&\gamma_{m-2}&\gamma_{m-1}&\gamma_m&\cdots&0&0\\ 0&\gamma_0&\cdots&\gamma_{m-3}&\gamma_{m-2}&\gamma_{m-1}&\cdots &0&0\\ \cdot&\cdot&\cdots&\cdot&\cdot&\cdot&\cdots&\cdot&\cdot\\ 0&0&\cdots&\gamma_1&\gamma_2&\gamma_3&\cdots&\gamma_m&0\\ 0&0&\cdots&\gamma_0&\gamma_1&\gamma_2&\cdots&\gamma_{m-1}&\gamma_m\\ \tau_0&\tau_1&\cdots&\tau_{m-2}&\tau_{m-1}&0&\cdots&0&0\\ 0&\tau_0&\cdots&\tau_{m-3}&\tau_{m-2}&\tau_{m-1}&\cdots&0&0\\ \cdot&\cdot&\cdots&\cdot&\cdot&\cdot&\cdots&\cdot&\cdot\\ 0&0&\cdots&\tau_0&\tau_1&\tau_2&\cdots&\tau_{m-1}&0\\ 0&0&\cdots&0&\tau_0&\tau_1&\cdots&\tau_{m-2}&\tau_{m-1} \end{array} \right ), $$

which is a polynomial equation in variable v. After finding all real roots of this polynomial, we can substitute them to (23) to find all the real solutions of w. Then, using $x_{1} = \frac{w}{\sqrt{1+w^{2}+v^{2}}}$, $x_{2} = \frac{v}{\sqrt{1+w^{2}+v^{2}}}$, $x_{3} = \frac{\pm1}{\sqrt{1+w^{2}+v^{2}}}$, λ=Ψ(x), we may find all the solutions of (18) in this case.

Combine all the possible solutions of (18) in these four cases, and find λ _min(As), the smallest value of λ of these solutions. This solves (7).