Institute of Physics and Engineering in Medicine Phys. Med. Biol. 60 (2015) 2339–2354
Physics in Medicine & Biology doi:10.1088/0031-9155/60/6/2339
Resolving intravoxel fiber architecture using nonconvex regularized blind compressed sensing C Y Chu1,2, J P Huang1,2, C Y Sun2, W Y Liu1,2 and Y M Zhu1,2 1
HIT-INSA Sino French Research Centre fssor Biomedical Imaging, Harbin Institute of Technology, Harbin, People’s Republic of China 2 CREATIS, CNRS UMR 5220, Inserm U630, INSA of Lyon, University of Lyon, Villeurbanne, France E-mail: [email protected] Received 21 August 2014, revised 10 November 2014 Accepted for publication 20 January 2015 Published 26 February 2015 Abstract
In diffusion magnetic resonance imaging, accurate and reliable estimation of intravoxel fiber architectures is a major prerequisite for tractography algorithms or any other derived statistical analysis. Several methods have been proposed that estimate intravoxel fiber architectures using low angular resolution acquisitions owing to their shorter acquisition time and relatively low b-values. But these methods are highly sensitive to noise. In this work, we propose a nonconvex regularized blind compressed sensing approach to estimate intravoxel fiber architectures in low angular resolution acquisitions. The method models diffusion-weighted (DW) signals as a sparse linear combination of unfixed reconstruction basis functions and introduces a nonconvex regularizer to enhance the noise immunity. We present a general solving framework to simultaneously estimate the sparse coefficients and the reconstruction basis. Experiments on synthetic, phantom, and real human brain DW images demonstrate the superiority of the proposed approach. Keywords: blind compressed sensing, diffusion-weighted imaging, intravoxel fiber orientations, nonconvex regularization, low angular resolution diffusion (Some figures may appear in colour only in the online journal) 1. Introduction Diffusion tensor imaging (DTI) is a recent noninvasive technique based on measuring the diffusion of water molecules in different directions in a given voxel (Basser et al 1994). It has 0031-9155/15/062339+16$33.00 © 2015 Institute of Physics and Engineering in Medicine Printed in the UK
2339
C Y Chu et al
Phys. Med. Biol. 60 (2015) 2339
been recently intensively studied to investigate the microstructures of tissues or organs such as brain (Alexander et al 2007, Assaf and Pasternak 2008) and heart (Wu et al 2006, 2007, Frindel et al 2009, 2010, Yang et al 2011, 2012, Wei et al 2013). The tensor model represents only one independent and dominant direction per voxel. Despite its simplicity and robustness, it has been shown to be incorrect in regions containing intravoxel orientational heterogeneity such as crossing and branching of fiber bundles. A number of methods have been proposed to overcome the single fiber orientation limitation of the DTI approach. These include the high angular resolution diffusion imaging (Frank 2002, Tuch et al 2002), generalized diffusion tensors imaging (Ozarslan and Mareci 2003, Liu et al 2004, 2010, Barmpoutis et al 2009), diffusion kurtosis imaging (Jensen et al 2005), Q-ball imaging (Tuch 2004, Hess et al 2006, Descoteaux et al 2007, Barnett 2009), diffusion spectrum imaging (Wedeen et al 2005, 2008, Gilbert et al 2006, Gaige et al 2010), hybrid diffusion imaging (Wu and Alexander 2007), etc. More details can be found in (Assemlal et al 2011, Hasan et al 2011). However, most of these approaches require relatively high b-values and a large number of gradient directions to acquire more detailed information. Such requirements are generally hard to meet in common clinical applications due to longer scan time and hardware constraints. There have been several methods, with which it is possible to resolve crossing fibers from conventional DTI acquisitions (with low b-values and about 30 gradient directions). A geometrically constrained two tensor model was proposed to fit diffusion-weighted (DW) signals at voxels susceptible to containing crossing fiber populations (Peled et al 2006). Independent component analysis was also introduced to fit a prolate tensor mixture (Singh 2005). In (Tournier et al 2007), nonnegativity constrained superresolved spherical deconvolution (CSD) was used to estimate the fiber orientation distribution. In (Jing et al 2012), a single channel blind source separation technique was employed to enhance the performance of the spherical deconvolution approach in low angular resolution DW data. Basis function decomposition (BFD) method was developed in (Ramirez-Manzanares and Rivera 2006) and (Ramirez-Manzanares et al 2007), in which the intravoxel information was recovered at voxels containing fiber crossings via the use of a linear combination of a predefined discrete tensor basis set, each of which represents a possible single fiber population. In (Landman et al 2012), a similar approach was introduced, where a constrained compressed sensing (CCS) method was used to infer complex fiber architecture through estimation of a set of sparse weights. However, these methods are highly sensitive to noise and risk over- or under-fitting. In general, the BFD (Ramirez-Manzanares et al 2007) and CCS (Landman et al 2012) methods have advantages due to not having to specify the number of distinct fiber populations, to their capacity of producing directly intravoxel fiber orientations and to their computational efficiency. However, in these methods, the fixed eigenvalues were used to generate the diffusion basis functions in the BFD method or to construct the reconstruction basis in the CCS method. The use of fixed eigenvalues introduces inaccuracies and unnatural limitations in them, since different tissue regions may present different diffusion properties. In addition, to address the noise sensitivity problem, regularization techniques have been used to reduce noise artifacts. For example, in (Ramirez-Manzanares et al 2007), a regularizer related to fiber’s spatial smoothness was defined as a convex function that makes the model resolving simpler but at the same time limits the effect of regularization. In this paper, we propose a nonconvex regularized blind compressed sensing (BCS) approach to resolve intravoxel fiber architectures from conventional low angular resolution DTI acquisitions. The idea is to use a multicompartment model similar to that in the 2340
C Y Chu et al
Phys. Med. Biol. 60 (2015) 2339
CCS and BFD methods, but in which the eigenvalues are not fixed, but estimated as for the other unknown parameters by employing the BCS approach. Meanwhile, we define a novel nonconvex regularizer allowing improving the quality of the results for noise corrupted data. 2. Methodology Our aim is to develop an approach, in which the nonconvex regularized BCS is used to better resolve intravoxel fiber orientations from low angular resolution DW imaging. To this end, two problems have to be resolved. First, low signal-to-noise ratio (SNR) of DW data makes the BCS model hardly convergent, therefore an efficient regularization must be used to improve the convergence. Second, the nonconvex regularization method requires solving high dimension global optimization problem which still remains a major challenge to the optimization research community. To tackle these problems, we decrease the dimension of nonconvex regularized BCS model by taking into account the fiber architectures of neighboring voxels, which allows the nonconvex regularized BCS model to be easily solved. 2.1. BCS-based fiber architecture resolving
Using multicompartment models (Frank 2002, Tuch et al 2002), an observed DW signal Si can be expressed as a finite mixture of Stejskal–Tanner tensor formulation M
S i = S0 ∑ f j exp (−bgTi Djgi) + ηi, i = 1, 2, … , N ,
(1)
j=1
where S0 is the nondiffusion-weighted signal of the voxel, gi the ith diffusion gradient, b the diffusion sensitization factor, N the number of gradient directions, M the number of possible compartments within each voxel, fj the apparent volume fraction of the voxel with diffusion tensor Dj, and ηi the Rician noise. Using eigen decomposition, the diffusion tensor can be expressed as D = RΛRT, where R is a rotation matrix whose column is the eigenvector of D, Λ is a diagonal matrix whose diagonal elements, λ1, λ2 and λ3, are the corresponding eigenvalues. In the present study, we assumed λ1 > λ2 = λ3, i.e. a cylindrical symmetry tensor. A diffusion tensor can then be expressed as ⎡ λ2 ⎤ ⎡ Δλ ⎤ ⎢ ⎥, T ⎢ ⎥ D R +⎢ λ2 0 =R ⎥ ⎢⎣ ⎥⎦ 0 λ2 ⎦ ⎣
(2)
⎡ Δλ ⎤ where Δλ = λ1 − λ2. Introducing (2) in (1) and let D Δλ = R ⎢ 0 ⎥ RT, the multicompart⎢⎣ ⎥ ment model can be written as 0⎦ M
(
)
S i = S0 ∑ f ′ j exp −bgiTD Δλ j g i + ηi, i = 1, 2, … , N , j=1
2341
(3)
C Y Chu et al
Phys. Med. Biol. 60 (2015) 2339
(
)
where f ′ j = f j exp(−bλ2 ). Set, BijΔλ = exp −bgiTD Δλ j g i , the model (3) can be rewritten in matrix form as y = BΔλ f ′ + η,
(4)
where each column of N × M matrix BΔλ is an atom of the reconstruction basis, the ith element of basis is sampling from a given tensor with eigenvalues (Δλ, 0, 0) along the ith gradient direction, M × 1 vector f′ is mixing coefficients, and N × 1 vector η is a noise term. Both Δλ and f′ can be estimated simultaneously by solving the BCS problem (Gleichman and Eldar 2011) as argmin ′ Δλ, f
{
BΔλ f ′ − y
2
+ β f′
1
}
s . t . f ′ j ≥ 0, ∀ j ,
(5)
where β > 0 is the sparsity regularization parameter which provides a tradeoff between the precision of model fitting and the sparsity requirement. Similar to classical CCS schemes (Landman et al 2012), the DW signal in a voxel was modeled as a sparse linear combination of basis tensors in a parametric dictionary. However, instead of assuming a dictionary with fixed parameters, the BCS scheme determines the parameters in a data driven way. In the present study, we use the set Ψ =
{ (f , v ) j = 1, 2, … , M} with v j
j
j
designating
the main eigenvector of diffusion tensor Dj to represent the fiber architecture in a voxel. Generally, most elements of f are zero. By removing these zero elements, Ψ can be simplified as Φ = αj, d j j = 1, 2, … , Ms , where αj is the jth nonzero element in f, dj the jth fiber orientation in a voxel, and Ms the number of nonzero elements in f (i.e. the number of fibers or fiber bundles in a voxel).
{(
}
)
2.2. Nonconvex regularization
We introduce the regularization technique in the process of fiber architecture estimation to reduce noise effect by assuming smoothness: neighboring voxels are assumed to be similar. To do that, we define a nonconvex regularizer which penalizes the dissimilarity of the coefficient distribution between neighboring voxels. We assume that Φ p = αi p, di p , i = 1, 2, … , Mp and Φq = α j q, d j q , j = 1, 2, … , Mq
{(
}
)
{(
)
}
represent the fiber architectures of voxels p and q, respectively. Their corresponding total Mp
Mq
i=1
j=1
amount of fiber volume fraction is Ap = ∑ αi p and Aq = ∑ α j q. To compute the dissimilarity between Φ p and Φq we need to match the closest fiber orientations because there are multiple fibers in a voxel. Generally, Ap is not equal to Aq. Therefore, we normalize αip ∼ = α / A and α ∼ = α / A , respectively. The fiber architecture dissimilarand αjq as α p q ip ip jq jq
(
)
ity Df Φ p, Φq is calculated as the solution of the following linear programming problem, in which the flow that minimizes an overall cost function subject to a set of constraints is computed
2342
C Y Chu et al
Phys. Med. Biol. 60 (2015) 2339 Mp Mq ⎧ ⎫ ⎪ ⎪ Df Φ p, Φq = min ⎨∑ ∑ xijDv (di p, dj q) ⎬ , ⎪ ⎪ x ⎩ i=1 j=1 ⎭
(
)
xij ≥ 0
s.t.
Mp
∑ xij = α∼jq
(6)
i=1 Mq
∑ xij = α∼ip j=1
Mp Mq
∑ ∑ xij = 1 i=1 j=1
(
)
where Dv (di p, dj q) = arccos diTpdj q is the angle between two fiber bundles with orientations dip and djq, respectively. The value range of Dv(dip, djq) is [0, π/2]. The equation system (6) is known as the Earth mover’s distance (EMD) (Rubner et al 2000) in computer science or Wasserstein metric in mathematics, which allows evaluating the dissimilarity between two multidimensional distributions in some feature space, i.e. a distance measure between single features. Intuitively, if each distribution is viewed as a unit amount of ‘dirt’ piled on a region, the EMD is the minimum ‘cost’ of turning one pile into the other, which is assumed to be the amount of dirt that needs to be moved times the distance it has to be moved. In the present study, we solve the optimization problem (6) using the fast algorithm in (Pele and Werman 2009). Based on the EMD distance, we define the regularizer as
( )
R Φ = nc p
∑
q ∈ Up
(
)
ω′pqDf Φ p, Φq / Z ′,
(7)
where the weighted factors ω′pq = p − q 2 is the Euclidean distance between p and q, Up is the neighborhood of p with radius h, and Z ′ = ∑ ω′pq is the sum of weighted factors. q ∈ Up
Introducing regularizer (7) in (5), we obtain the nonconvex regularized BCS problem as follows. argmin ′ Δλ, f
{
BΔλ f ′ − y
2
+ β1 f ′ 1 + β2Rnc (Φ )
}
s . t . f ′ j ≥ 0, ∀ j ,
(8)
where the nonnegative control parameters β1 and β2 allow us to tune the amount of regularization. 2.3. Numerical solution framework
Generally, to reduce the reconstruction error, the dimension of f′ (which equals the number of atoms of reconstruction basis) in (8) should be greater than 200 (Landman et al 2012), which makes it difficult to solve the optimization problem (8). To overcome such difficulty, we propose a novel numerical solution framework, the main idea of which consists of decreasing the dimension of the optimization problem by taking into account the fiber architecture of neighboring voxels, thus making the problem to be more easily solved. 2343
C Y Chu et al
Phys. Med. Biol. 60 (2015) 2339
For the initialization step of the procedure, we select a fixed value Δλp(0) for each voxel p, and the fiber architecture Φ p(0) of voxel p is estimated by solving the L1 optimization problem argmin (0) fp
{ Bf p − y (0)
2
+ β f (0) p
1
} s.t.
f j(0) p ≥ 0, ∀ j .
(9)
Then, at each iteration k, to each voxel p, considering every neighboring voxel q ∈ Up with fiber architecture Φq(k − 1) = α (j kq− 1), d(jkq− 1) , j = 1, 2, … , Ms(qk − 1) , we define a fiber orienta-
{(
}
)
(k )
tion set Θ p (voter) which contains all fiber orientations of the neighboring voxel set Up. So Θ p(k ) = d(jkq−i 1), j = 1, 2, … , Ms(qk −i 1), i = 1, 2, … , NUp , where NUp is the number of voxels
{
}
in Up. After that, we vote for reconstruction matrix B (a N × M matrix, where N designates the number of gradient directions and M the number of possible compartments within each voxel): assuming the rth atom Br (rth column of B) with main eigenvector vr, if vrTd(jkq−i 1) is smaller than a voting threshold ta, we increase the vote for Br by one. So the total vote of Br represents how many fiber orientations in Θ p(k ) are close to vr. By choosing Mv columns having the largest (number of total votes from B, we construct a dimension-reduced reconstruction Δλ k − 1) matrix Bp p . We can thus solve the following optimization problem using global optimization methods such as an evolutionary search strategy named SP-UCI (Chu et al 2011) we used in the present study ⎧ Δλ (k−1) argmin ⎨ Bp p f (pk ) − y (k ) ⎩ fp
( )
where R′nc Φ p(k ) =
∑
q ∈ Up
2
⎫ + β1 f (pk ) + β2R′nc Φ p(k ) ⎬ s . t . f j(pk ) ≥ 0, ∀ j , 1 ⎭
( )
(
(10)
)
ω′pqDf Φ p(k ), Φq(k − 1) / Z ′.
Finally, at kth iteration, to each voxel p, we solve the following problem to estimate Δλp(k ) ⎧ ⎫ ⎪ ⎪ Δλp(k ) (k ) ⎨ ∑ Bx f x − y ⎬ , argmin ⎪ 2⎪ Δλp(k ) ⎩ x ∈ U ′p ⎭
(11)
where U ′p is the neighborhood of voxel p with radius H. The problem (11) can be solved efficiently since only one variable is concerned. The iteration stopping criterion used is when k equals the predefined number of iterations kN. 2.4. Performance optimization
To improve the convergence speed, we propose a group estimation strategy (GES) that is based on two ideas. First, whenever possible, we use more reliable fiber architecture estimation Φq(k ) instead of Φq(k − 1) of neighboring voxels to estimate the fiber architecture of voxel p at kth iteration. Second, we preferentially estimate the voxels with simple architectures (single orientation, estimated at (k − 1)th iteration). Based on the first idea, we rewrite the ‘voter’ Θ p(k ) and the regularizer R′nc Φ p(k ) respectively as
( )
2344
C Y Chu et al
Phys. Med. Biol. 60 (2015) 2339
⎧ k, Φq(ki ) has been estimated , Θ′p(k ) = d(jxq)i, j = 1, 2, … , Ms(qx )i , i = 1, 2, … , NUp , x = ⎨ otherwise ⎩ k − 1,
{
and
( )
R? nc Φ p(k ) =
}
⎧ k, Φq(k ) has been estimated . ω? pqDf Φ p(k ), Φq(x ) / Z ′, x = ⎨ otherwise ⎩ k − 1, q ∈ Up
(
∑
)
(12)
⎪
⎪
(13)
Based on the second idea, we define the group priority of a voxel p by PG = NS + DC, where NS is the number of voxels with only one fiber orientation in cubic 3 × 3 × 3 neighborhood of p and DC is the distance between p and the nearest voxel with multiple fiber orientations. We preferentially estimate the voxels with high group priority, while the voxels with the same priority estimated in order from left-top-anterior to right-bottom-posterior. The optimized numerical solution framework is summarized as follows: (a) Initialization: select the parameters Δλp(0), β1, β2, ta, Mv, h, H, kN, and estimate fiber structure Φ p(0) for every voxel p using (9). Set k = 1. (b) Iteration: Group voxels using the proposed grouping method. For each voxel p ordered by the grouping priority, vote to construct dimension-reduced Δλ (k − 1)
reconstruction basis Bp p ⎧ Δλ (k−1) argmin ⎨ Bp p f (pk ) − y (k ) ⎩ fp
2
using Θ′p(k ), and then estimate Φ p(k ) by
⎫ + β1 f (pk ) + β2R″nc Φ p(k ) ⎬ s.t. f j(pk ) ≥ 0, ∀ j , 1 ⎭
( )
(14)
For each voxel p, estimate Δλp(k ) via (11). (c) Stopping: stop if k = kN, otherwise increase k by 1 and continue with step (b). 2.5. The choice of parameters
The proposed approach involves several parameters, so we briefly discuss in the following their use in practice. The initialization parameter Δλp(0) does not affect the results owing to the adaptivity of the proposed framework. However, a reasonable initialization value can increase the convergence speed. Our default choice is Δλp(0) = 0.8. The parameter Mv represents the number of atoms in the reduced reconstruction matrix. A too small Mv will lead to abnormal smoothing of the fibers, whereas a too large Mv will increase computational workload since larger Mv means higher dimension in the optimization problem (14). In the present study, we choose Mv = 10–15 empirically, knowing that a slightly higher Mv can be a good choice for higher noise level or complex fiber architecture situations. The voting threshold ta should be chosen such that the number of atoms having votes should be slightly greater than Mv. As described in section 2.3, we compare vrTd(jkq−i 1) with the voting threshold ta. This is equivalent to comparing the angle between vr and d(jkq−i 1) with arccos(ta). If vrTd(jkq−i 1)
Physics in Medicine & Biology doi:10.1088/0031-9155/60/6/2339
Resolving intravoxel fiber architecture using nonconvex regularized blind compressed sensing C Y Chu1,2, J P Huang1,2, C Y Sun2, W Y Liu1,2 and Y M Zhu1,2 1
HIT-INSA Sino French Research Centre fssor Biomedical Imaging, Harbin Institute of Technology, Harbin, People’s Republic of China 2 CREATIS, CNRS UMR 5220, Inserm U630, INSA of Lyon, University of Lyon, Villeurbanne, France E-mail: [email protected] Received 21 August 2014, revised 10 November 2014 Accepted for publication 20 January 2015 Published 26 February 2015 Abstract
In diffusion magnetic resonance imaging, accurate and reliable estimation of intravoxel fiber architectures is a major prerequisite for tractography algorithms or any other derived statistical analysis. Several methods have been proposed that estimate intravoxel fiber architectures using low angular resolution acquisitions owing to their shorter acquisition time and relatively low b-values. But these methods are highly sensitive to noise. In this work, we propose a nonconvex regularized blind compressed sensing approach to estimate intravoxel fiber architectures in low angular resolution acquisitions. The method models diffusion-weighted (DW) signals as a sparse linear combination of unfixed reconstruction basis functions and introduces a nonconvex regularizer to enhance the noise immunity. We present a general solving framework to simultaneously estimate the sparse coefficients and the reconstruction basis. Experiments on synthetic, phantom, and real human brain DW images demonstrate the superiority of the proposed approach. Keywords: blind compressed sensing, diffusion-weighted imaging, intravoxel fiber orientations, nonconvex regularization, low angular resolution diffusion (Some figures may appear in colour only in the online journal) 1. Introduction Diffusion tensor imaging (DTI) is a recent noninvasive technique based on measuring the diffusion of water molecules in different directions in a given voxel (Basser et al 1994). It has 0031-9155/15/062339+16$33.00 © 2015 Institute of Physics and Engineering in Medicine Printed in the UK
2339
C Y Chu et al
Phys. Med. Biol. 60 (2015) 2339
been recently intensively studied to investigate the microstructures of tissues or organs such as brain (Alexander et al 2007, Assaf and Pasternak 2008) and heart (Wu et al 2006, 2007, Frindel et al 2009, 2010, Yang et al 2011, 2012, Wei et al 2013). The tensor model represents only one independent and dominant direction per voxel. Despite its simplicity and robustness, it has been shown to be incorrect in regions containing intravoxel orientational heterogeneity such as crossing and branching of fiber bundles. A number of methods have been proposed to overcome the single fiber orientation limitation of the DTI approach. These include the high angular resolution diffusion imaging (Frank 2002, Tuch et al 2002), generalized diffusion tensors imaging (Ozarslan and Mareci 2003, Liu et al 2004, 2010, Barmpoutis et al 2009), diffusion kurtosis imaging (Jensen et al 2005), Q-ball imaging (Tuch 2004, Hess et al 2006, Descoteaux et al 2007, Barnett 2009), diffusion spectrum imaging (Wedeen et al 2005, 2008, Gilbert et al 2006, Gaige et al 2010), hybrid diffusion imaging (Wu and Alexander 2007), etc. More details can be found in (Assemlal et al 2011, Hasan et al 2011). However, most of these approaches require relatively high b-values and a large number of gradient directions to acquire more detailed information. Such requirements are generally hard to meet in common clinical applications due to longer scan time and hardware constraints. There have been several methods, with which it is possible to resolve crossing fibers from conventional DTI acquisitions (with low b-values and about 30 gradient directions). A geometrically constrained two tensor model was proposed to fit diffusion-weighted (DW) signals at voxels susceptible to containing crossing fiber populations (Peled et al 2006). Independent component analysis was also introduced to fit a prolate tensor mixture (Singh 2005). In (Tournier et al 2007), nonnegativity constrained superresolved spherical deconvolution (CSD) was used to estimate the fiber orientation distribution. In (Jing et al 2012), a single channel blind source separation technique was employed to enhance the performance of the spherical deconvolution approach in low angular resolution DW data. Basis function decomposition (BFD) method was developed in (Ramirez-Manzanares and Rivera 2006) and (Ramirez-Manzanares et al 2007), in which the intravoxel information was recovered at voxels containing fiber crossings via the use of a linear combination of a predefined discrete tensor basis set, each of which represents a possible single fiber population. In (Landman et al 2012), a similar approach was introduced, where a constrained compressed sensing (CCS) method was used to infer complex fiber architecture through estimation of a set of sparse weights. However, these methods are highly sensitive to noise and risk over- or under-fitting. In general, the BFD (Ramirez-Manzanares et al 2007) and CCS (Landman et al 2012) methods have advantages due to not having to specify the number of distinct fiber populations, to their capacity of producing directly intravoxel fiber orientations and to their computational efficiency. However, in these methods, the fixed eigenvalues were used to generate the diffusion basis functions in the BFD method or to construct the reconstruction basis in the CCS method. The use of fixed eigenvalues introduces inaccuracies and unnatural limitations in them, since different tissue regions may present different diffusion properties. In addition, to address the noise sensitivity problem, regularization techniques have been used to reduce noise artifacts. For example, in (Ramirez-Manzanares et al 2007), a regularizer related to fiber’s spatial smoothness was defined as a convex function that makes the model resolving simpler but at the same time limits the effect of regularization. In this paper, we propose a nonconvex regularized blind compressed sensing (BCS) approach to resolve intravoxel fiber architectures from conventional low angular resolution DTI acquisitions. The idea is to use a multicompartment model similar to that in the 2340
C Y Chu et al
Phys. Med. Biol. 60 (2015) 2339
CCS and BFD methods, but in which the eigenvalues are not fixed, but estimated as for the other unknown parameters by employing the BCS approach. Meanwhile, we define a novel nonconvex regularizer allowing improving the quality of the results for noise corrupted data. 2. Methodology Our aim is to develop an approach, in which the nonconvex regularized BCS is used to better resolve intravoxel fiber orientations from low angular resolution DW imaging. To this end, two problems have to be resolved. First, low signal-to-noise ratio (SNR) of DW data makes the BCS model hardly convergent, therefore an efficient regularization must be used to improve the convergence. Second, the nonconvex regularization method requires solving high dimension global optimization problem which still remains a major challenge to the optimization research community. To tackle these problems, we decrease the dimension of nonconvex regularized BCS model by taking into account the fiber architectures of neighboring voxels, which allows the nonconvex regularized BCS model to be easily solved. 2.1. BCS-based fiber architecture resolving
Using multicompartment models (Frank 2002, Tuch et al 2002), an observed DW signal Si can be expressed as a finite mixture of Stejskal–Tanner tensor formulation M
S i = S0 ∑ f j exp (−bgTi Djgi) + ηi, i = 1, 2, … , N ,
(1)
j=1
where S0 is the nondiffusion-weighted signal of the voxel, gi the ith diffusion gradient, b the diffusion sensitization factor, N the number of gradient directions, M the number of possible compartments within each voxel, fj the apparent volume fraction of the voxel with diffusion tensor Dj, and ηi the Rician noise. Using eigen decomposition, the diffusion tensor can be expressed as D = RΛRT, where R is a rotation matrix whose column is the eigenvector of D, Λ is a diagonal matrix whose diagonal elements, λ1, λ2 and λ3, are the corresponding eigenvalues. In the present study, we assumed λ1 > λ2 = λ3, i.e. a cylindrical symmetry tensor. A diffusion tensor can then be expressed as ⎡ λ2 ⎤ ⎡ Δλ ⎤ ⎢ ⎥, T ⎢ ⎥ D R +⎢ λ2 0 =R ⎥ ⎢⎣ ⎥⎦ 0 λ2 ⎦ ⎣
(2)
⎡ Δλ ⎤ where Δλ = λ1 − λ2. Introducing (2) in (1) and let D Δλ = R ⎢ 0 ⎥ RT, the multicompart⎢⎣ ⎥ ment model can be written as 0⎦ M
(
)
S i = S0 ∑ f ′ j exp −bgiTD Δλ j g i + ηi, i = 1, 2, … , N , j=1
2341
(3)
C Y Chu et al
Phys. Med. Biol. 60 (2015) 2339
(
)
where f ′ j = f j exp(−bλ2 ). Set, BijΔλ = exp −bgiTD Δλ j g i , the model (3) can be rewritten in matrix form as y = BΔλ f ′ + η,
(4)
where each column of N × M matrix BΔλ is an atom of the reconstruction basis, the ith element of basis is sampling from a given tensor with eigenvalues (Δλ, 0, 0) along the ith gradient direction, M × 1 vector f′ is mixing coefficients, and N × 1 vector η is a noise term. Both Δλ and f′ can be estimated simultaneously by solving the BCS problem (Gleichman and Eldar 2011) as argmin ′ Δλ, f
{
BΔλ f ′ − y
2
+ β f′
1
}
s . t . f ′ j ≥ 0, ∀ j ,
(5)
where β > 0 is the sparsity regularization parameter which provides a tradeoff between the precision of model fitting and the sparsity requirement. Similar to classical CCS schemes (Landman et al 2012), the DW signal in a voxel was modeled as a sparse linear combination of basis tensors in a parametric dictionary. However, instead of assuming a dictionary with fixed parameters, the BCS scheme determines the parameters in a data driven way. In the present study, we use the set Ψ =
{ (f , v ) j = 1, 2, … , M} with v j
j
j
designating
the main eigenvector of diffusion tensor Dj to represent the fiber architecture in a voxel. Generally, most elements of f are zero. By removing these zero elements, Ψ can be simplified as Φ = αj, d j j = 1, 2, … , Ms , where αj is the jth nonzero element in f, dj the jth fiber orientation in a voxel, and Ms the number of nonzero elements in f (i.e. the number of fibers or fiber bundles in a voxel).
{(
}
)
2.2. Nonconvex regularization
We introduce the regularization technique in the process of fiber architecture estimation to reduce noise effect by assuming smoothness: neighboring voxels are assumed to be similar. To do that, we define a nonconvex regularizer which penalizes the dissimilarity of the coefficient distribution between neighboring voxels. We assume that Φ p = αi p, di p , i = 1, 2, … , Mp and Φq = α j q, d j q , j = 1, 2, … , Mq
{(
}
)
{(
)
}
represent the fiber architectures of voxels p and q, respectively. Their corresponding total Mp
Mq
i=1
j=1
amount of fiber volume fraction is Ap = ∑ αi p and Aq = ∑ α j q. To compute the dissimilarity between Φ p and Φq we need to match the closest fiber orientations because there are multiple fibers in a voxel. Generally, Ap is not equal to Aq. Therefore, we normalize αip ∼ = α / A and α ∼ = α / A , respectively. The fiber architecture dissimilarand αjq as α p q ip ip jq jq
(
)
ity Df Φ p, Φq is calculated as the solution of the following linear programming problem, in which the flow that minimizes an overall cost function subject to a set of constraints is computed
2342
C Y Chu et al
Phys. Med. Biol. 60 (2015) 2339 Mp Mq ⎧ ⎫ ⎪ ⎪ Df Φ p, Φq = min ⎨∑ ∑ xijDv (di p, dj q) ⎬ , ⎪ ⎪ x ⎩ i=1 j=1 ⎭
(
)
xij ≥ 0
s.t.
Mp
∑ xij = α∼jq
(6)
i=1 Mq
∑ xij = α∼ip j=1
Mp Mq
∑ ∑ xij = 1 i=1 j=1
(
)
where Dv (di p, dj q) = arccos diTpdj q is the angle between two fiber bundles with orientations dip and djq, respectively. The value range of Dv(dip, djq) is [0, π/2]. The equation system (6) is known as the Earth mover’s distance (EMD) (Rubner et al 2000) in computer science or Wasserstein metric in mathematics, which allows evaluating the dissimilarity between two multidimensional distributions in some feature space, i.e. a distance measure between single features. Intuitively, if each distribution is viewed as a unit amount of ‘dirt’ piled on a region, the EMD is the minimum ‘cost’ of turning one pile into the other, which is assumed to be the amount of dirt that needs to be moved times the distance it has to be moved. In the present study, we solve the optimization problem (6) using the fast algorithm in (Pele and Werman 2009). Based on the EMD distance, we define the regularizer as
( )
R Φ = nc p
∑
q ∈ Up
(
)
ω′pqDf Φ p, Φq / Z ′,
(7)
where the weighted factors ω′pq = p − q 2 is the Euclidean distance between p and q, Up is the neighborhood of p with radius h, and Z ′ = ∑ ω′pq is the sum of weighted factors. q ∈ Up
Introducing regularizer (7) in (5), we obtain the nonconvex regularized BCS problem as follows. argmin ′ Δλ, f
{
BΔλ f ′ − y
2
+ β1 f ′ 1 + β2Rnc (Φ )
}
s . t . f ′ j ≥ 0, ∀ j ,
(8)
where the nonnegative control parameters β1 and β2 allow us to tune the amount of regularization. 2.3. Numerical solution framework
Generally, to reduce the reconstruction error, the dimension of f′ (which equals the number of atoms of reconstruction basis) in (8) should be greater than 200 (Landman et al 2012), which makes it difficult to solve the optimization problem (8). To overcome such difficulty, we propose a novel numerical solution framework, the main idea of which consists of decreasing the dimension of the optimization problem by taking into account the fiber architecture of neighboring voxels, thus making the problem to be more easily solved. 2343
C Y Chu et al
Phys. Med. Biol. 60 (2015) 2339
For the initialization step of the procedure, we select a fixed value Δλp(0) for each voxel p, and the fiber architecture Φ p(0) of voxel p is estimated by solving the L1 optimization problem argmin (0) fp
{ Bf p − y (0)
2
+ β f (0) p
1
} s.t.
f j(0) p ≥ 0, ∀ j .
(9)
Then, at each iteration k, to each voxel p, considering every neighboring voxel q ∈ Up with fiber architecture Φq(k − 1) = α (j kq− 1), d(jkq− 1) , j = 1, 2, … , Ms(qk − 1) , we define a fiber orienta-
{(
}
)
(k )
tion set Θ p (voter) which contains all fiber orientations of the neighboring voxel set Up. So Θ p(k ) = d(jkq−i 1), j = 1, 2, … , Ms(qk −i 1), i = 1, 2, … , NUp , where NUp is the number of voxels
{
}
in Up. After that, we vote for reconstruction matrix B (a N × M matrix, where N designates the number of gradient directions and M the number of possible compartments within each voxel): assuming the rth atom Br (rth column of B) with main eigenvector vr, if vrTd(jkq−i 1) is smaller than a voting threshold ta, we increase the vote for Br by one. So the total vote of Br represents how many fiber orientations in Θ p(k ) are close to vr. By choosing Mv columns having the largest (number of total votes from B, we construct a dimension-reduced reconstruction Δλ k − 1) matrix Bp p . We can thus solve the following optimization problem using global optimization methods such as an evolutionary search strategy named SP-UCI (Chu et al 2011) we used in the present study ⎧ Δλ (k−1) argmin ⎨ Bp p f (pk ) − y (k ) ⎩ fp
( )
where R′nc Φ p(k ) =
∑
q ∈ Up
2
⎫ + β1 f (pk ) + β2R′nc Φ p(k ) ⎬ s . t . f j(pk ) ≥ 0, ∀ j , 1 ⎭
( )
(
(10)
)
ω′pqDf Φ p(k ), Φq(k − 1) / Z ′.
Finally, at kth iteration, to each voxel p, we solve the following problem to estimate Δλp(k ) ⎧ ⎫ ⎪ ⎪ Δλp(k ) (k ) ⎨ ∑ Bx f x − y ⎬ , argmin ⎪ 2⎪ Δλp(k ) ⎩ x ∈ U ′p ⎭
(11)
where U ′p is the neighborhood of voxel p with radius H. The problem (11) can be solved efficiently since only one variable is concerned. The iteration stopping criterion used is when k equals the predefined number of iterations kN. 2.4. Performance optimization
To improve the convergence speed, we propose a group estimation strategy (GES) that is based on two ideas. First, whenever possible, we use more reliable fiber architecture estimation Φq(k ) instead of Φq(k − 1) of neighboring voxels to estimate the fiber architecture of voxel p at kth iteration. Second, we preferentially estimate the voxels with simple architectures (single orientation, estimated at (k − 1)th iteration). Based on the first idea, we rewrite the ‘voter’ Θ p(k ) and the regularizer R′nc Φ p(k ) respectively as
( )
2344
C Y Chu et al
Phys. Med. Biol. 60 (2015) 2339
⎧ k, Φq(ki ) has been estimated , Θ′p(k ) = d(jxq)i, j = 1, 2, … , Ms(qx )i , i = 1, 2, … , NUp , x = ⎨ otherwise ⎩ k − 1,
{
and
( )
R? nc Φ p(k ) =
}
⎧ k, Φq(k ) has been estimated . ω? pqDf Φ p(k ), Φq(x ) / Z ′, x = ⎨ otherwise ⎩ k − 1, q ∈ Up
(
∑
)
(12)
⎪
⎪
(13)
Based on the second idea, we define the group priority of a voxel p by PG = NS + DC, where NS is the number of voxels with only one fiber orientation in cubic 3 × 3 × 3 neighborhood of p and DC is the distance between p and the nearest voxel with multiple fiber orientations. We preferentially estimate the voxels with high group priority, while the voxels with the same priority estimated in order from left-top-anterior to right-bottom-posterior. The optimized numerical solution framework is summarized as follows: (a) Initialization: select the parameters Δλp(0), β1, β2, ta, Mv, h, H, kN, and estimate fiber structure Φ p(0) for every voxel p using (9). Set k = 1. (b) Iteration: Group voxels using the proposed grouping method. For each voxel p ordered by the grouping priority, vote to construct dimension-reduced Δλ (k − 1)
reconstruction basis Bp p ⎧ Δλ (k−1) argmin ⎨ Bp p f (pk ) − y (k ) ⎩ fp
2
using Θ′p(k ), and then estimate Φ p(k ) by
⎫ + β1 f (pk ) + β2R″nc Φ p(k ) ⎬ s.t. f j(pk ) ≥ 0, ∀ j , 1 ⎭
( )
(14)
For each voxel p, estimate Δλp(k ) via (11). (c) Stopping: stop if k = kN, otherwise increase k by 1 and continue with step (b). 2.5. The choice of parameters
The proposed approach involves several parameters, so we briefly discuss in the following their use in practice. The initialization parameter Δλp(0) does not affect the results owing to the adaptivity of the proposed framework. However, a reasonable initialization value can increase the convergence speed. Our default choice is Δλp(0) = 0.8. The parameter Mv represents the number of atoms in the reduced reconstruction matrix. A too small Mv will lead to abnormal smoothing of the fibers, whereas a too large Mv will increase computational workload since larger Mv means higher dimension in the optimization problem (14). In the present study, we choose Mv = 10–15 empirically, knowing that a slightly higher Mv can be a good choice for higher noise level or complex fiber architecture situations. The voting threshold ta should be chosen such that the number of atoms having votes should be slightly greater than Mv. As described in section 2.3, we compare vrTd(jkq−i 1) with the voting threshold ta. This is equivalent to comparing the angle between vr and d(jkq−i 1) with arccos(ta). If vrTd(jkq−i 1)
