We have presented a hemispherical mesh-based formulation for spherical deconvolution to recover fiber orientations directly from DWI signal attenuations. This method has several unique characteristics, which we have highlighted in this report. It ensures that the estimated FOD will strictly satisfy the physical constraints of realness, symmetry, and non-negativity. Moreover, the final solution will have unit mass and be globally optimal, subject to constraints, with respect to the objective function embodied in Eq. (1)
. Our method also provides flexibility in the form of the regularization parameters τ
, as well as in the choice of the regularization matrix D
. This flexibility does not imply the lack of robustness, however, as we have shown that an appropriately-chosen set of regularization parameters can provide reasonable results for multiple data sets, both real and simulated.
The primary advantages of our mesh-based formulation are in its guarantee of physically-valid FOD reconstruction and its straightforward, well-characterized approach to obtaining these results. As noted, previous attempts have been made to minimize negative values in FOD reconstruction obtained using the spherical harmonics framework, but such efforts require an ad hoc approach which has a more obscure interpretation than the principled gradient projection method we propose here. We also note that enforcing a non-negativity criterion through the frequency domain may have an undesirable behavior because the spherical harmonic basis functions are periodic and globally supported, and thus manipulating the coefficients to diminish negative values in one region will necessarily change the FOD values over all of S2. By applying the non-negativity constraint directly in the spatial domain, the mesh-based approach enables local adjustments to FOD shape without distant consequences. Moreover, even the modified spherical harmonics-based deconvolution approaches cannot completely eliminate negative values in the final solution, an important consideration for probabilistic or front-evolution-based tractography methods that sample the FOD at each step, which may be forced to invoke ad hoc logic themselves to handle such events. Moreover, as we have shown that such ad hoc methods may reduce the accuracy of the FOD, especially in regions of high anisotropy ( and , right), it is conceivable that the mesh-based spherical deconvolution approach has the potential to improve downstream tractography results and ultimately our understanding of brain connectivity.
Secondary advantages include beneficial side effects of the mesh-based representation itself. Determining the most likely fiber direction, for use in deterministic tractography, becomes trivial. Similarly, the mesh-based representation permits straightforward comparison of orientation differences between FODs in neighboring voxels, a metric which can be used to enable broader-scale spatial regularization of fiber tracts. For these reasons, and because we have not observed significant differences in angular resolution between our mesh-based approach and CSD, for example, we suggest that the mesh-based representation may expedite future analyses without sacrificing accuracy, albeit at the cost of a less compact representation.
The mesh-based approach also presents some unique challenges, particularly with regards to discretization. For the i
th vertex, the value xiwi
represents the likelihood of fibers oriented in directions which lie within some solid angle surrounding that vertex; the size of this solid angle is determined by the mesh geometry. This is similar to the interpretation of the discrete dODF in Tuch (2004)
. For some applications, however, it may be desirable to know the value of the FOD in any arbitrary direction. This can be obtained from piecewise linear interpolation of x
on the triangulated mesh. Because x
as weighted by w
has unit integral over the mesh, the continuous function resulting from interpolating all arbitrary directions on S2
in this manner will also have unit mass and retain the properties of a probability distribution, including non-negativity.
The maxima of the FOD, however, will always lie on the original mesh points following piecewise linear interpolation. For downstream applications such as deterministic tractography methods which rely solely on the maxima of the distribution, this raises a significant concern about possible discretization error. Here we have attempted to minimize this error by using a very dense mesh. For a 1281 point mesh, the angle between neighboring vertices is slightly greater than 4°, resulting in a maximum discretization error of approximately 2°. This represents less than one standard deviation of the error we observe in simulations, even at high SNR (). However, computation of the solution for a dense mesh is necessarily expensive as the convolution matrix A contains as many columns as there are mesh vertices (for the 60-directional data sets in this work, a single Pipeline CPU solves 100 voxels in just under 30 s). If possible, disregarding non-brain (or perhaps non-white-matter) voxels can reduce the computation time. We propose two additional approaches to further mitigate this burden. First, as the deconvolution operation treats each voxel independently, once the convolution matrix A is determined, the operation can be parallelized in a straightforward manner to compute multiple voxels simultaneously. This is the approach we have utilized in our Pipeline implementation. Second, one might use an adaptive re-meshing strategy where the deconvolution is carried out initially on a sparse mesh, and then again on a mesh which contains additional vertices near the previously detected maxima. Indeed, as we have observed little difference in angular accuracy and resolution between the “projected” and “clipped” methods, the mesh-based approach may not be the most economical choice for applications which make use of only the maxima of the FOD, rather than the full distribution.
The process of choosing appropriate regularization parameters also warrants further discussion. Throughout this report, we have defined D
as the weighted spatial difference matrix over neighboring mesh vertices. As noted previously, however, the regularization matrix is a free parameter, and the user may select it in accordance with prior knowledge about the form of the solution. Accordingly, a seemingly obvious choice for D
might be the identity, suggesting that the desired FOD is sparse. We note that this choice, along with letting p
=1, converts Eq. (1)
into the form commonly encountered in “compressed-sensing” methods (Candes and Tao, 2006
; Donoho, 2006
). While it is true that the FOD is zero over most of its support, this choice of regularization parameters, in our experience, produces appropriately-oriented, but highly noisy FODs due to the ill-posed nature of the spherical deconvolution. Furthermore, we observe that for p
=2, Eq. (1)
takes the classical Tikhonov-regularized form, and an explicit solution is available. Although this closed-form solution does not guarantee non-negativity, it may be useful as an initial guess for beginning projected gradient descent. Also, in this paper, we have relied on the multi-fiber simulation of Eq. (6)
to select optimal regularization parameters for non-synthetic data sets. We have managed to obtain reasonable results with this approach, although this model is certainly an oversimplification of the true diffusion process (Inglis et al., 2001
; Maier et al., 2001
; Assaf et al., 2004
). In addition, we note that, while we have tuned these parameters at the higher b
-value (3000 s/mm2
) used for simulation in previous works (Sakaie and Lowe, 2007
; Tournier et al., 2007
), we generate satisfactory results from data sets acquired at more conventional b
-values (1159 s/mm2
). This is likely largely a consequence of the fact that we obtain the single fiber response function by sampling each data set individually. Despite this robustness, it may be worthwhile, in some applications, to optimize the tuning for a particular data set by examining results for several values of τ
for a few selected voxels from regions with known anatomy before computing the solution for the full volume.
We have noted a relative scarcity of clinical reports utilizing HARDI reconstructions and associated metrics, relative to those incorporating more traditional DTI analyses. No doubt this imbalance results partially from the increased technical requirements of HARDI acquisition, but also in part due to a sparsity of widely available tools for conducting such analyses. To this end, we hope that the software we have made available via the Pipeline will be of use to others in testing our spherical deconvolution approach on a variety of data sets for methodological comparisons and validations, as well as for quantitative investigations of clinical populations.
Future directions for investigation using the mesh-based formulation include analysis of more complex convex optimization methods in order to minimize computational load. Newton–Raphson methods, quasi-Newton methods, and conjugate gradient methods all have theoretical speed advantages for solving Eq. (1)
; it will be important to examine whether these can be realized in practice while still imposing the necessary constraints. In addition, the incorporation of more sophisticated prior knowledge into the regularization matrix D
may lead to improved FOD reconstructions. Finally, we note that the mesh-based representation lends itself naturally to neighborhood-guided spatial regularization, and further exploration along may lead to improved or more robust tractography results.