- Email: [email protected]
A highly accurate symmetric optical flow based high-dimensional nonlinear spatial normalization of brain images
Magnetic Resonance Imaging xxx (2015) xxx–xxx
Contents lists available at ScienceDirect
Magnetic Resonance Imaging journal homepage: www.mrijournal.com
A highly accurate symmetric optical flow based high-dimensional nonlinear spatial normalization of brain images☆ Ying Wen a,⁎, Lili Hou a, Lianghua He b, Bradley S. Peterson c, Dongrong Xu d,⁎ a
Shanghai Key Laboratory of Multidimensional Information Processing & Department of Computer Science and Technology, East China Normal University, Shanghai, 200241, China Key Laboratory of Embedded System and Service Computing, Ministry of Education, Tongji University, Shanghai 200092, China Institute for the Developing Mind, Children's Hospital Los Angeles, The Keck School of Medicine at the University of Southern California, Los Angeles, CA 90027, U.S.A. d MRI Unit & Epidemiology Division, Department of Psychiatry, Columbia University & New York State Psychiatric Institute, New York, 10032, U.S.A. b c
a r t i c l e
i n f o
Article history: Received 9 May 2014 Revised 21 October 2014 Accepted 18 January 2015 Available online xxxx Keywords: Brain image normalization Symmetric optical flow Hierarchical strategy Optical flow constraints
a b s t r a c t Spatial normalization plays a key role in voxel-based analyses of brain images. We propose a highly accurate algorithm for high-dimensional spatial normalization of brain images based on the technique of symmetric optical flow. We first construct a three dimension optical model with the consistency assumption of intensity and consistency of the gradient of intensity under a constraint of discontinuity-preserving spatio-temporal smoothness. Then, an efficient inverse consistency optical flow is proposed with aims of higher registration accuracy, where the flow is naturally symmetric. By employing a hierarchical strategy ranging from coarse to fine scales of resolution and a method of Euler-Lagrange numerical analysis, our algorithm is capable of registering brain images data. Experiments using both simulated and real datasets demonstrated that the accuracy of our algorithm is not only better than that of those traditional optical flow algorithms, but also comparable to other registration methods used extensively in the medical imaging community. Moreover, our registration algorithm is fully automated, requiring a very limited number of parameters and no manual intervention. © 2015 Elsevier Inc. All rights reserved.
1. Introduction Many methods in image analysis have been developed to tackle diffusion tensor images (DTI) registration [1–3]. DTI registration generally falls into two categories based on setting up the correspondence between two subjects. The first involves feature-based matching, i.e. transformations that are calculated based on a number of anatomical correspondences established manually, semi-automatically, or fully automatically on a number of distinct anatomical features, including distinct landmarks [4,5] or a combination of curves and surfaces, such as sulci and gyri [6,7]. In addition, finding anatomical correspondences is a key to the success of the spatial normalization. To find anatomical correspondences, conventional methods generally first extract scalar features from each tensor individually; and then by constructing scalar maps, regional integration and other operations such as edge detection, a final set of features for correspondence matching can be constituted. For example, Timer [8] and F-Timer [9] were developed ☆ This work was supported by The National Natural Science Foundation of China (No. 61273261, 61272267, 81471734), Program for New Century Excellent Talents in University (NCET-11-0381), Shanghai Collaborative Innovation Center of Trustworthy Software for Internet of Things (ZF1213), NIDA 1R01DA027100 and NIMH 5R01MH089582. ⁎ Corresponding authors. E-mail addresses: [email protected] (Y. Wen), [email protected] (B.S. Peterson), [email protected] (D. Xu).
based on Hammer [10], which adopts so-called driving voxels in addition to geometric moments as the features upon which registration is based. The boundary-based registration uses the measure of the thickness of cortex in T1-weighted image [11]. Yang et al. guided their registration using information on a tensor’s structural geometry and orientation [12]. Xue et al. proposed a local fast marching method for DTI registration [13]. These landmark-based methods numerically seek optimal transformations (under respective model-specific assumptions) to maximize similarity across brains. The second category of registration methods is based on volumetric transformations between one brain image and a template image. These methods assume that the images are acquired using the same pulse sequence [14–16] and these methods therefore are generally more efficient, because they do not require the construction of a specific anatomical model using anatomic landmarks. Christensen et al. [17] proposed a viscous fluid flow model to implement registration. Unlike fluid flow, optical flow is another volumetric transformation method. Optical flow is usually posed as an energy minimization problem, with the energy function containing a data term and a smoothness term. The variational framework, together with coarse-to-fine refinement, is widely used in optical flow estimation. The technique of optical flow developed in computer vision for object tracking [18,19] is one such method. The two approaches differ with the constraints in their implementations. The former adopts viscous constraints and body
http://dx.doi.org/10.1016/j.mri.2015.01.013 0730-725X/© 2015 Elsevier Inc. All rights reserved.
Please cite this article as: Wen Y, et al, A highly accurate symmetric optical flow based high-dimensional nonlinear spatial normalization of brain images, Magn Reson Imaging (2015), http://dx.doi.org/10.1016/j.mri.2015.01.013
2
Y. Wen et al. / Magnetic Resonance Imaging xxx (2015) xxx–xxx
Fig. 1. Fractional anisotropy maps of the 20 randomly selected subjects from our database. The brains vary significantly in morphometry.
force while the latter adopts gradient and smooth constants. For the solution of optical flow, Xu et al. proposed an extended coarse-to-fine refinement framework which reduces the reliance of flow estimation on coarse level and corrects the optical flow in each scale [20]. Stoll et al. proposed an adaptive integration strategy for descriptor matching by selecting the positions where descriptor matching may improve the optical flow estimation [21]. Besides optical flow method, diffeomorphism methods are often used for image registration. The registration methods used extensively in the medical imaging community include symmetric diffeomorphic image registration (SyN) [22], SPM5 DARTEL [23], as well as with Thirion’s Demons algorithm [24] etc. The optimization of diffeomorphism method is implemented on a space of diffeomorphic space instead of the entire space of displacement field. SyN maximizes the cross-correlation in the space of diffeomorphic maps. DARTEL is a diffeomorphic image registration for estimating inverse consistent deformations. Demons algorithm uses an approximate elastic regularizer to solve an optical flow problem. Traditional techniques of optical flow estimation do not generally yield symmetrical solutions: the results will differ if they are applied between images f1 and f2 or between images f2 and f1. Ideally, the forward and inverse flows should be uniquely determined and should be inverses of one another. Estimating the forward and inverse flows independently very rarely results in a consistent set of transformations due to a large number of local minima. To overcome this deficiency in current optical flow, symmetrical formulations of the optical flow has been proposed in Refs. [29–31], where the solution is a constraint to being symmetric using a combination of the flow in both directions. In this work, we present a simple and
effective method to recover a dense optical flow field map from two images. The idea is that it exists as a symmetric displacement between displacement vectors from f1 to f2 and f2 to f1 and to minimize an energy function of optical flow. This variational problem is then solved using the gradient flow defined by the Euler–Lagrange equations associated to the energy. This idea is introduced for DTI brain registration in this paper. We focus on the optical flow based high-dimensional nonlinear spatial normalization of brain images. We first extend the traditional 2D optic flow to 3D form for high-dimensional nonlinear spatial normalization. Then, we propose a simple and effective symmetric variational formulation of the optical flow, where the flow is naturally symmetric for the problem of the forward and inverse flows. Finally, we incorporate the symmetric optical model into a hierarchical framework that is applied to imaging data in multiple resolution scales from coarse to fine, thereby building up a registration approach with high accuracy. 2. Materials and methods 2.1. Data acquisition We acquired all the human Diffusion-weighted imaging (DWI) datasets from a GE 3.0 Tesla MR scanner system (General Electric, Milwaukee, Wisconsin) at New York State Psychiatry Institute, Columbia University. This unit was equipped with a high strength (50 mT/m) and high-speed gradient system (slew rate = 200 T/m/s) capable of conducting DWI. A single-shot spin echo, echo planar imaging (SE–EPI) sequence was used to acquire the DWI data, with
Please cite this article as: Wen Y, et al, A highly accurate symmetric optical flow based high-dimensional nonlinear spatial normalization of brain images, Magn Reson Imaging (2015), http://dx.doi.org/10.1016/j.mri.2015.01.013
Y. Wen et al. / Magnetic Resonance Imaging xxx (2015) xxx–xxx
diffusion gradients applied at b = 1000 s/mm2 along 25 non-collinear spatial directions plus 3 baseline images without applying diffusion gradient. The thickness of each slice was 2.5 mm without gap, and about 68 axial slices were assigned to cover the entire brain. The other parameters of the DWI sequence were minimal echo time (TE) at about 77.3 ms; repetition time (TR) = 17000 ms; field of view (FOV) = 240 mm × 240 mm and acquisition data matrix =132 × 128. The final images were machine-interpolated to an in-plane resolution of 256 × 256 voxels, resulting in a nominal in-plane resolution of 0.9375 mm × 0.9375 mm. The total DWI scanning time was about 16 minutes 26 seconds with the number of excitation (NEX) = 2. Twenty data sets of healthy volunteers were randomly selected from our database. These brains vary significantly in morphology (Fig. 1). 2.2. Methods Our algorithm generally follows the conventional framework for image registration, i.e. rigid global transformation plus nonlinear registration. The global transformation roughly aligns an individual image to a selected template image in terms of its size and spatial orientation. By extending the standard 2D optic flow method to a 3D form, we propose in this paper a registration method based on the 3D optical flow method for local nonlinear registration of the DTI data which is high-dimensional and incorporates an energy function based on the three mentioned assumptions. In this process, Euler-Lagrange equation [25] is employed to optimize the energy function, and the successive over-relaxation method (SOR) [26] is used to search for the optimal solution. This framework calculates the deformation filed for normalizing DTI data using diffusion anisotropy index (DAI) images that reflect the underlying diffusion property. The same framework then normalizes the DTI data using this calculated deformation field that correlates the spatial correspondence between the two spaces. The DAI used in our case is the fractional anisotropy (FA) map.
minimizing the energy based on this L1 norm requires that the function is a convex. Because in this form ε is a small positive constant, ψ(e 2)remains a convex which meets the mentioned requirement in the process of searching minimum. Moreover, this form of ψ(e 2) does not introduce any additional parameters, because ε is only for numerical calculation and can be fixed to a constant value (0.001 in our case). Finally, we incorporate a smoothness term into the energy function to satisfy the assumption that the flow field must be piecewise smooth. The smoothness term is necessary to cope with the function’s ill-posedness and noise which introduce destructive effects in the optical flow approach. In the optical flow equation, the variation of the pixels’ intensity is interpreted as their motion. Without the motion smoothness term, the estimated motion of adjacent pixels can be completely disparate or significantly deviate from the truth. The smoothness constraint is achieved by penalizing the total variation of the flow field, which can be expressed as: Z Esmooth ðu; v; wÞ ¼
We now present our nonlinear registration method by deducting the equations using the 3D form of the optic flow. Suppose that s := (x, y, z, t) T is the voxel in the image at location (x, y, z) at time t and p := (u, v, w, 1) T is the deformation vector field, the data term is: 2
2
j f ðs þ pÞ− f ðsÞj þ βj∇f ðs þ pÞ−∇f ðsÞj
ds
ð1Þ
Because quadratic penalisers may generate outliers that may excessively bias the estimation of optic flow, we introduce into this energy function an increasing concave function ψ(e2), revising it into a more robust energy function [28], although the function ψ can also be applied separately to each of these terms: Z Eðu; v; wÞ ¼
S
2
Z Etotal ðu; v; wÞ ¼
S
2
ð3Þ
2
2
ψðj f ðs þ pÞ− f ðsÞj þ βj∇f ðs þ pÞ−∇ f ðsÞj Þ 2
2.2.2. Three dimensional optical flow algorithm
S
2
ψðj∇uj þ j∇vj þ j∇wj Þds
2
ð4Þ
2
þλψðj∇uj þ j∇vj þ j∇wj Þds
The global transformation is a rigid affine transformation. In general, an affine transformation is composed of linear transformations (rotation, scaling or shear) and a translation (or “shift”). The brain in its native image space is transformed to the space of the template using the global affine transformation [27], as a proper initialization for the subsequent coregistration process.
Z
S
using the same function for Ψ as above. The spatio-temporal gradient ∇ := (∂x, ∂y, ∂z, ∂t)′ incorporates the needed spatio-temporal smoothness assumption. For coregistration that involves only two images, this gradient term should be replaced by one that involves only the spatial gradient ∇ := (∂x, ∂y, ∂z)′. Thus, the total energy is the weighted sum of (2) and (3):
2.2.1. Global transformation
Eðu; v; wÞ ¼
3
2 2 ψ j f ðs þ pÞ− f ðsÞj þ βj∇f ðs þ pÞ−∇f ðsÞj ds
ð2Þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi We adopt the form ψ e2 ¼ e2 þ ε2 in this energy calculation, which results in (modified) L1 minimization [29]. However,
where, λ is a positive regularization parameter. Now we need to solve this equation for searching the best solution. 2.2.3. The proposed optical flow model In the standard formulation of variational optical flow, the estimated motion vector fields depend on the template image and are asymmetric. However, in most application the solution should be independent of the template image. Motivated by [30–32], we use a symmetric optical flow to settle the problem of consistency. For image fm between first image f1 and second image f2, there exist an optical flow p that satisfies:
f m ðsÞ ¼ f ðs−p=2Þ ¼ f ðs þ p=2Þ
ð5Þ
Fig. 2 shows a symmetric flow description. The images in the first row of Fig. 2 display the flow from f1 to f2 in the traditional way, and the images in the second row display the flow in a symmetric way. Then, (4) can be written as 1 2 ψðj f ðs þ p=2Þ−f ðs−p=2Þj Þþ B C Eðu; v; wÞ ¼ @ βψðj∇f ðs þ p=2Þ−∇ f ðs−p=2Þj2 Þþ Ads Ω 2 2 2 λψðj∇uj þ j∇vj þ j∇wj Þ Z
0
ð6Þ
Please cite this article as: Wen Y, et al, A highly accurate symmetric optical flow based high-dimensional nonlinear spatial normalization of brain images, Magn Reson Imaging (2015), http://dx.doi.org/10.1016/j.mri.2015.01.013
4
Y. Wen et al. / Magnetic Resonance Imaging xxx (2015) xxx–xxx
We name ft the difference image between the two consecutive image frames; fi (i = x, y, z) is the first order derivative of the movement image f(s + p); fij(i,j = x, y, z) is the second order derivative of the movement image and fit (i = x, y, z) is the first order of the difference image. To solve (7) for the variables u, v, w, we follow the construction of the numerical scheme performed in [12]. The first step is to introduce the fixed point iteration on p. Let p k = (u k, v k, w k,1) T the solution of the system, i.e. the deformation field, at the kth iteration step and let f⁎⁎k be the abbreviations defined in (8), an implicit scheme for the smoothness term and a semi-implicit scheme for the data term have been applied. This gives the following system of equations at the (k + 1)th iteration step: Fig. 2. An explanation of symmetric optical flow.
2.2.4. Numerical solution to the energy function Using the Euler–Lagrange equations, the energy function (6)can be transformed into the following form:
8 0 2 2 2 2 > > ψ f t þ β f xt þ f yt þ f zt > > > > > > f x f t þ β f xx f xt þ f xy f yt þ f xz f zt > > > > > 0 2 2 2 > > −λdiv ψ j∇uj þ j∇vj þ j∇wj ∇u ¼ 0; > > > > 0 2 2 2 2 > > ψ f t þ β f xt þ f yt þ f zt > > < f y f t þ β f xy f xt þ f yy f yt þ f yz f zt > > > 0 2 2 2 > > −λdiv ψ j∇uj þ j∇vj þ j∇wj ∇v ¼ 0; > > > > 0 2 2 2 2 > > ψ f t þ β f xt þ f yt þ f zt > > > > > > > ⋅ f z f t þ β f xz f xt þ f yz f yt þ f zz f zt > > > > > : −λdiv ψ0 j∇uj2 þ j∇vj2 þ j∇wj2 ∇w ¼ 0
8 > 0 kþ1 2 kþ1 2 kþ1 2 kþ1 2 > Þ þ β f xt þ f yt þ f zt ψ ft > > > > > > k kþ1 k kþ1 k kþ1 k kþ1 > > ⋅ f x f t þ β f xx f xt þ f xy f yt þ f xz f zt > > > > > kþ1 2 > 0 kþ1 2 kþ1 2 kþ1 > ¼ 0; −λdiv ψ ∇u þ ∇v þ ∇w ∇u > > > > > > > 0 kþ1 2 kþ1 2 kþ1 2 kþ1 2 > > ft þ β f xt þ f yt þ f zt > ψ > > < k kþ1 k kþ1 k kþ1 k kþ1 ⋅ f y f t þ β f xy f xt þ f yy f yt þ f yz f zt > > > 2 2 > kþ1 0 kþ1 kþ1 2 kþ1 > > ¼ 0; −λdiv ψ ∇u þ ∇v þ ∇w ∇v > > > > 2 2 2 2 > > 0 kþ1 kþ1 kþ1 kþ1 > > ft þ β f xt þ f yt þ f zt > > ψ > > > > > ⋅ f k f kþ1 þ β f k f kþ1 þ f k f kþ1 þ f k f kþ1 > > z t xz xt yz yt zz zt > > > > kþ1 2 0 kþ1 2 kþ1 2 kþ1 > > ¼0 : −λdiv ψ ∇u þ ∇v þ ∇w ∇w
ð9Þ
ð7Þ This system is still nonlinear because of the nonlinear function ψ′ and the symbols f⁎k + 1. To facilitate the subsequent calculation, we apply the first order Taylor expansions to remove the nonlinearity from the components f⁎k + 1: kþ1 k k k ≈f s þ p þ fx s þ p du f s þ p k k k k þfy s þ p dv þ fz s þ p dw
where we use the following abbreviation:
∂ f ðs þ p=2Þ þ ∂x f ðs−p=2Þ fx ¼ x 2 ∂y f ðs þ p=2Þ þ ∂y f ðs−p=2Þ fy ¼ 2 ∂ f ðs þ p=2Þ þ ∂z f ðs−p=2Þ fz ¼ z 2 f t ¼ f ðs þ p=2Þ−f ðs−p=2Þ ∂ f ðs þ p=2Þ þ ∂xx f ðs−p=2Þ f xx ¼ xx 2 ∂yy f ðs þ p=2Þ þ ∂yy f ðs−p=2Þ f yy ¼ 2 ∂ f ðs þ p=2Þ þ ∂zz f ðs−p=2Þ f zz ¼ zz 2 ∂xy f ðs þ p=2Þ þ ∂xy f ðs−p=2Þ f xy ¼ 2 ∂ f ðs þ p=2Þ þ ∂xz f ðs−p=2Þ f xz ¼ xz 2 ∂yz f ðs þ p=2Þ þ ∂yz f ðs−p=2Þ f yz ¼ 2 f xt ¼ ∂x f ðs þ p=2Þ−∂x f ðs−p=2Þ f yt ¼ ∂y f ðs þ p=2Þ−∂y f ðs−p=2Þ f zt ¼ ∂z f ðs þ p=2Þ−∂z f ðs−p=2Þ
The unknown iteration variable p k + 1 is thus substituted by the known variable p k and an unknown term dp k = (du k, dv k, dw k, 1) T. Therefore, the derivative of the data term and the derivative of the smoothness term in (7) become:
ð8Þ
0 k 0 k k k k k k k 2 f t þ f xt du þ f yt dv þ f zt dw ψ D :¼ ψ k k k k k k k 2 þβ f xt þ f xx du þ f xy dv þ f xz dw k k k k k k k 2 þ f yt þ f xy du þ f yy dv þ f yz dw k k k k k k k 2 ; þ f zt þ f xz du þ f yz dv þ f zz dw 2 0 k 0 k k ψ S :¼ ψ ∇ u þ du 2 k k k k 2 þ ∇ v þ dv þ ∇ w þ dw
ð10Þ
where (ψ ')Dk can be interpreted as a robustness factor in the data term, and (ψ ')Sk as diffusivity in the smoothness term.
Please cite this article as: Wen Y, et al, A highly accurate symmetric optical flow based high-dimensional nonlinear spatial normalization of brain images, Magn Reson Imaging (2015), http://dx.doi.org/10.1016/j.mri.2015.01.013
Y. Wen et al. / Magnetic Resonance Imaging xxx (2015) xxx–xxx
Subsequently, (9) can be written as 8 0 k k k k k k k k k > ψ D f x f t þ f xt du þ f yt dv þ f zt dw > > > > > k k k k k k k k > > > þβ f xx f xt þ f xx du þ f xy dv þ f xz dw > > > > > þf kxy f kyt þ f kxy duk þ f kyy dvk þ f kyz dwk > > > > k k > k k k k k k > > þf xz f zt þ f xz du þ f yz dv þ f zz dw > > > > > −λdiv ψ0 k ∇ uk þ duk ¼ 0; > > S > > > 0 k k k > k k k k k k > ψ f > y f t þ f xt du þ f yt dv þ f zt dw > D > > > > > þβ f kxy f kxt þ f kxx duk þ f kxy dvk þ f kxz dwk > > < k k k k k k k k þf yy f yt þ f xy du þ f yy dv þ f yz dw > > > k k k k k k k k > > þf yz f zt þ f xz du þ f yz dv þ f zz dw > > > > 0 k k k > > ∇ v þ dv ¼ 0; > S > −λdiv ψ > > 0 k > k k k k k k k k > > ψ D f z f t þ f xt du þ f yt dv þ f zt dw > > > > k k k k k k k k > > þβ f xz f xt þ f xx du þ f xy dv þ f xz dw > > > > k k > k k k k k k > þf f þ f xy du þ f yy dv þ f yz dw > > > yz yt > > k k k k k k k k > > þf zz f zt þ f xz du þ f yz dv þ f zz dw > > > > > : −λdiv ψ0 k ∇ wk þ dwk ¼ 0 S
multiresolution approach captures large-scale deformation in the lower-resolution steps and small deformation in the finer resolution steps. After a designated number of iterations are carried out at each resolution level, or when convergence is observed, the estimated motion field is passed to the next level of a super sampled resolution using cubic interpolation, which serves as the initialization for the new estimation at the finer level of image resolution. 3. Experiments
ð11Þ
However, (11) still contains nonlinear components ψ' D and ψ' S. In order to remove the remaining nonlinearity in ψ ', a second, inner, fixed point iteration loop is applied. Let du k,n and dv k,n denote the iteration variables at step n; and let (ψ ')Dk,n and (ψ ')Sk,n denote the robustness factor and the diffusivity defined in (11) at iteration k and n (the iteration regarding k is for removing the linearity in f⁎k + 1 and the iteration regarding n is for removing the linearity in ψ' D and ψ' S). Then, the final linear system of equations in du k,n + 1 and dv k,n + 1 reads 8 0 k;n k k k k;nþ1 k k;nþ1 k k;nþ1 > ψ D f x f t þ f xt du þ f yt dv þ f zt dw > > > > > k k k k;nþ1 k k;nþ1 k k;nþ1 > > þ f xy dv þ f xz dw > þβ f xx f xt þ f xx du > > > > > þf kxy f kyt þ f kxy duk;nþ1 þ f kyy dvk;nþ1 þ f kyz dwk;nþ1 > > > > k k > k k;nþ1 k k;nþ1 k k;nþ1 > > þf xz f zt þ f xz du þ f yz dv þ f zz dw > > > > > −λdiv ψ0 k;n ∇ uk þ duk;nþ1 ¼ 0; > > S > > > 0 k;n k k > k k;nþ1 k k;nþ1 k k;nþ1 > f þ f yt dv þ f zt dw ψ > y f t þ f xt du > D > > > > > þβ f kxy f kxt þ f kxx duk;nþ1 þ f kxy dvk;nþ1 þ f kxz dwk;nþ1 > > < k k k k;nþ1 k k;nþ1 k k;nþ1 þf yy f yt þ f xy du þ f yy dv þ f yz dw > > > k k k k;nþ1 k k;nþ1 k k;nþ1 > > þf yz f zt þ f xz du þ f yz dv þ f zz dw > > > > 0 k;n k k;nþ1 > > ¼ 0; −λdiv ψ S ∇ v þ dv > > > > 0 k;n k k > k k;nþ1 k k;nþ1 k k;nþ1 > > ψ D f z f t þ f xt du þ f yt dv þ f zt dw > > > > k k k k;nþ1 k k;nþ1 k k;nþ1 > > þβ f xz f xt þ f xx du þ f xy dv þ f xz dw > > > > > k k k k;nþ1 k k;nþ1 k k;nþ1 > þf yz f yt þ f xy du þ f yy dv þ f yz dw > > > > > k k k k;nþ1 k k;nþ1 k k;nþ1 > > þf zz f zt þ f xz du þ f yz dv þ f zz dw > > > > > : −λdiv ψ0 k;n ∇ wk þ dwk;nþ1 ¼ 0 S
5
ð12Þ
Eq. (12) is actually in the form of A(du, dv, dw) T = B. Using standard discretization for the derivatives, the resulting sparse linear system of equations can now be solved using common numerical methods, such as Gauss-Seidel or SOR iterations. We employed a hierarchical multiresolution scheme to accommodate different scales of motion for improving computational efficiency. The multiresolution approach speeds up the convergence of the iteration represented by (12). We also expect that the
We conducted three different experiments to evaluate from various aspects the performance of our algorithm. The first experiment examined the accuracy of the deformation fields generated by our algorithm using simulated datasets. The second experiment demonstrated the excellent performance of our algorithm using real-world data acquired directly from human participants. The third experiment evaluated the performance of our algorithm by conducting a comparison with other algorithms using fiber tractography. All our experiments were conducted using the following empirical parameters λ = 40; β = 60; number of SOR iterations = 500, image dimension: 256 × 256 × Z (Z varied from 50 to 65 slices), resolution = 0.9375 mm × 0.9375 mm × 2.5 mm. Datasets with motion artifacts and outliers were not included. 3.1. Simulated experiments In the experiment, we aimed to evaluate quantitatively the performance of our algorithm using simulated data in terms of the morphological accuracy of spatial normalization. We inspected the registration error of our algorithm for coregistering two brain images for their structural details, and therefore focused on the morphological information. The distinct features relevant particularly to tensor data were thus temporarily not considered, which would be the focus of the subsequent experiments. To estimate the proposed registration model accuracy, we adopted a strategy similar to that was reported in [19]. We first randomly selected from our DTI database of brain imaging data six morphological representative datasets that were acquired from 6 healthy human participants. These 6 datasets differed from each other in size, shape and structure. We then randomly chose one of the 6 as the template brain, and calculated the FA maps for all the 6 datasets. Next, we warped the FA of the template brain to each FA of the other 5 datasets individually, using SPM5 DARTEL [22]. We thus obtained 5 versions of the template brain in 5 different individual spaces (Fig. 3), along with 5 deformation fields DFT2I that maps the correspondence between the template space and the individual brains. Although the 5 configurations of the template brain differed with each other in size and shape in different spaces, they shared exactly the same structure as these 5 were basically the same brain. A good algorithm for spatial normalization should be able to accurately bring them back to the original template space with perfect registration to the original template brain. We applied our algorithm thus to test the morphological accuracy that our algorithm could offer. We compared the performance of our optical flow based algorithm with several other algorithms including the affine transform, the SPM5 DARTEL, the Horn-Schunck’s algorithm [15] and the Lucas-kanade algorithm [16]. In the template brain, we manually marked 200 landmarks on those visually distinct structural features. The deformation fields DFT2I mapped them to 1000 locations (200 each for 5 brains) in the individual spaces; and the deformation fields DFI2T generated by our algorithm took all the 5 deformed brains back to the template space, thus mapping these 1000 locations to 5 new sets of 200 locations in the original template space. By measuring the average distance between the 1000 new
Please cite this article as: Wen Y, et al, A highly accurate symmetric optical flow based high-dimensional nonlinear spatial normalization of brain images, Magn Reson Imaging (2015), http://dx.doi.org/10.1016/j.mri.2015.01.013
6
Y. Wen et al. / Magnetic Resonance Imaging xxx (2015) xxx–xxx
(a)
(b1)
(b2)
(b3)
(b4)
(b5)
Fig. 3. Images of the synthetic brains for inspecting the registration performance of our method. (a) The template brain. (b1–b5) Five deformed versions of the template brain in 5 individual spaces.
locations and their corresponding landmarks we originally plotted, we would be able to evaluate the registration error of using our algorithm. Suppose P(x, y, z) are the predefined landmarks in the template brain, P ' (x ', y ', z ') the corresponding locations in the subject spaces mapped by DFT2I, and P 00 (x 00 , y 00 , z 00 ) the final locations back to the template space generated by DFI2T, the morphological accuracy of the normalization algorithm was measured by the mean and the standard deviation of distance between P(x, y, z) and P 00 (x 00 , y 00 , z 00 ) as follows: md ¼
1 X
P ðx; y; zÞ−P 00 x00 ; y00 ; z00
jV j x;y;z∈V
ð13Þ
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
1 X
P ðx; y; zÞ−P 00 x00 ; y00 ; z00 −m 2 d jV j x;y;z∈V
ð14Þ
and σd ¼
where ‖⋅‖ is the Euclidean distance, |V| the number of landmarks we defined in the template brain. To provide a more general evaluation of the registration error, we also calculated the mean and standard deviation of the distance for every brain voxel in the whole image and compared the results generated by all the methods under comparison. The quantitative report indicated that the registration error of our method was only 0.85 mm (mean) and 0.88 mm (std. var) for the 1000 anatomical landmarks, and 0.98 mm and 0.94 mm for all the brain voxels (Table 1), both smaller than the size of one voxel. In contrast, among all the other four methods except ours (DARTEL, affine registration, the Horn-Schunck method and the Lucas-kanade method) under comparison, DARTEL performed the best but the errors were significantly larger than ours (Table 1). The results demonstrated that our method performed superior to the other methods in coregistering the image for the structural details. 3.2. Experiments on DTI data
is to inspect the sharpness of an average image obtained from a set of coregistered images. We randomly selected from our DTI database 20 healthy brains for this experiment. These brains vary significantly in morphology. One of the 20 brains was then randomly selected as the template for all the rest 19 to be spatially normalized. We applied also the affine registration and DARTEL registration to compare with the performance of our algorithm. After all datasets were normalized into the template space, we averaged them. In addition to visual inspection of the sharpness of the averaged data, we also performed quantitatively evaluation of the sharpness for the whole 3D image by measuring the standard deviation of the FA value from the reference in the template at each voxel and generated a standard deviation map for each method that was under comparison. The mean image of the registered brains showed that compared with the results generated by DARTEL and the affine registration, our method generated the clearest structural boundaries in the results (Fig. 4), indicating that our method performed the best coregisteration of the individual images. Our method not only registered nearly perfectly the general morphology of the brains, it also registered very well the small structures inside of the brain. In the following experiment, the tracking results appeared that the tracked fibers based on the data normalized using our method resampled the ground truth much better than did the results based on the data normalized using the DARTEL method (Fig. 5). We obtained 693 fibers passing through the splenium ROI and 512 fibers passing through the genu ROI of a subject shown in Fig. 5(b). After the data normalized using DARTEL and our method, we obtained 699 and 689 fibers passing through the splenium ROI and 348 and 490 fibers passing through the genu ROI shown in Fig. 5(c)
Table 1 A comparison of the registration errors using different methods (mm). Measurement
In this subsection, we wanted to evaluate the morphological accuracy of our algorithm in practice directly using real-world DTI data acquired from human participants. If a registration method is not accurate, the structural boundary of the average image of a number of registered images will definitely be fuzzy. Therefore, a direct and intuitive way of gauging registration error of an algorithm
Horn-Schunck Lucas-kanade Affine registration DARTEL Our method
Landmarks
All voxels
Mean
Std.
Mean
Std.
2.47 3.14 4.38 1.32 0.85
2.58 2.66 3.02 0.98 0.88
3.12 3.48 4.65 1.45 0.98
3.37 3.54 3.78 1.24 0.94
Please cite this article as: Wen Y, et al, A highly accurate symmetric optical flow based high-dimensional nonlinear spatial normalization of brain images, Magn Reson Imaging (2015), http://dx.doi.org/10.1016/j.mri.2015.01.013
Y. Wen et al. / Magnetic Resonance Imaging xxx (2015) xxx–xxx
(a)
(b)
(c)
7
(d)
Fig. 4. A comparison of the performance between our method and other approaches, using 20 human brain images. (a) The template brain. (b) The mean brain of the registration results generated by the affine registration. (c) The mean brain of the registration results generated by DARTEL registration. (d) The mean brain of the registration results generated by our method.
and (d), respectively. Compared the result based on DARTEL shown in Fig. 5(c), the topology of the result based on our approach as shown in Fig. 5(d) better resembled the ground truth in Fig. 5(b). For example, the y-shaped fiber branches (pointed by the white arrows) were well reserved using our method but disappeared in the data generated by DARTEL; moreover, DARTEL generated lots of originally non-existing fiber pathways (pointed by red arrows). In particular, the data processed using our method almost duplicated the fiber structures existing in the original data. In
(a)
(b)
contrast, the data using DARTEL generated lots of non-existing fiber pathways and destroyed lots that originally existed. Similarly, the quantitative analysis again demonstrated that the results based on our method approximated the ground truth quite well, while the results based on SPM5 DARTEL had significant difference from the ground truth (Table 2). In general, the statistics showed that the results by DARTEL had a much larger standard deviation than those using our method. To validate the registration performance of our method, we compared the proposed method with greyscale image based
(c)
(d)
Fig. 5. The results of fiber tracking, based on the two regions-of-interest using the diffusion tensor imaging data normalized using either DARTEL or our approach. (a) A triplanar view (viewed from top left of the brain) of the tracked fibers in the original template brain that passed through the two ROIs, viewed from above the right ear. The fibers were superimposed on the corresponding FA image. (b)The same as that shown in (b), but only the tracked fibers are displayed here, and the view was generated from top of the brain. (c)The tracked fibers in the data normalized using DARTEL. (d) The tracked fibers in the data normalized using our method.
Please cite this article as: Wen Y, et al, A highly accurate symmetric optical flow based high-dimensional nonlinear spatial normalization of brain images, Magn Reson Imaging (2015), http://dx.doi.org/10.1016/j.mri.2015.01.013
8
Y. Wen et al. / Magnetic Resonance Imaging xxx (2015) xxx–xxx
Table 2 A comparison of the statistics of the fiber tracking results based on the two seeding volumes defined on the genu and the splenium of the corpus callosum (The statistics in parentheses denotes the difference between the testing item and the original template). Object
Template
Fibers_bundle
Splenium
Genu
DARTEL Splenium
Genu
Splenium
Genu
Mean fiber length Std of length Mean fiber angle Std of angle Mean fiber count Std of count Mean volume Std of volume
49.84 / 3.09 / 693 / 14604 /
32.05 / 4.88 / 512 / 5384 /
44.69 1.76 3.06 0.025 699 18.19 14 760 262.21
31.86 1.54 4.98 0.21 348 27.65 4182 182.62
49.81 1.37 3.09 0.02 689 12.91 14 512 202.63
31.38 1.57 4.84 0.06 490 11.88 5039 152.36
methods, such as symmetric diffeomorphic image registration (SyN) [21], SPM5 DARTEL [22], as well as with Thirion’s Demons algorithm [23] used extensively in the medical imaging community. SyN algorithm was implemented with the extended version of the ITK deformable image registration frame, described in Yoo [33]. Thirion’s Demons algorithm was known to perform well in inter-subject deformable image registration and used an approximate elastic regularizer to solve an optical flow problem. Demons algorithm was freely available in the Insight ToolKit and has been optimized by the ITK community (www.itk.org). We considered the above three registration methods for three reasons: (1) we wanted to evaluate brain registration algorithms, not brain labeling algorithms or particular atlas-based approaches, (2) these methods are based on greyscale image but not feature points, and (3) their performances are good according to Klein’s report [34]. In this experiment, each of the 20 brain images was registered to the 19 others, resulting in 380 pairwise registrations. Each of the brains is represented 19 times as the registration source and 19 times as the target. We segmented the brain image into three tissues: white matter (WM), gray matter (GM), cerebrospinal fluid (CSF) using the software FSL [35–37]. Then, we adopted two measurements to estimate the effectiveness of the registration: Target overlap [38] and Dice coefficient [39] and their definitions are: T S T ρ ¼ O \ O O
=
ð15Þ
and T T S S ρ ¼ 2 O \ O = O þ O
ð16Þ
where O T and O S denote the ROI in the template and subject image, respectively. Target overlap and Dice coefficient measure varies in the range [0, 1]. Table 3 presented a registration performance comparison of four registration methods. The average target overlay and dice coefficient and their standard deviation of four methods were listed in Table 3.
Our method
The performance of Demons algorithm was weaker than that of others. It can be seen that the performances of DARTEL, SyN and our methods were very close to each other. The target overlap from our method was roughly 0.835 in case of WM, 0.722 for GM and 0.807 for CSF. The dice coefficient from our method was roughly 0.820 in case of WM, 0.683 for GM and 0.811 for CSF. These values were slightly higher than or close to the results using other algorithms. As a whole, our method was comparable to other methods used extensively in the medical imaging community.
4. Conclusion We have developed a new approach based on three dimensional symmetric optic flow models for deformable registration of brain images. The proposed algorithm takes FA maps to generate the point correspondence between the subject and template spaces. We designed it so as to integrate the DTI-specific characteristics into the registration framework. This framework may not be suitable for neuroimaging DTI studies that actually study FA variances in diseases or across populations, but will be fine for those studying other DAIs (e.g., D?, DO, MD, etc.), because this framework is supposed to already remove all possible FA and morphometric variances during the normalization procedure. However, the framework is ready to use other replacements instead of the FA images for generating the point correspondence. One alterative is the baseline data acquired along with the DWI data, which are actually heavily T2-weighted images. Our algorithm is a fully automated process of registration for brain image data and the software implements the approach presented in http://www.columbia.edu/~dx2103/brainimagescope. html. Compared to other advanced algorithms currently available such as SPM5 DARTEL [22], HAMMER [7] and Timer [5], our algorithm is relatively more convenient and user-friendly with comparably high accuracy and efficiency. We plan to further evaluate and validate the method in larger datasets, especially pediatric databases in the future.
Table 3 The registration results comparison of demons, SPM5 dartel, SYN and our method. Method
Demons
DARTEL
SyN
Our method
0.831 ± 0.038 0.705 ± 0.052 0.785 ± 0.127
0.837 ± 0.024 0.711 ± 0.067 0.788 ± 0.138
0.835 ± 0.025 0.722 ± 0.059 0.807 ± 0.135
Target overlap WM GM CSF
0752 ± 0.120 0666 ± 0.121 0.741 ± 0.166
Demons
DARTEL
SyN
Our Method
0.817 ± 0.058 0.661 ± 0.067 0.791 ± 0.106
0.821 ± 0.047 0.672 ± 0.087 0.790 ± 0.091
0.820 ± 0.049 0.683 ± 0.069 0.811 ± 0.097
Dice coefficient 0767 ± 0.164 0.620 ± 0.110 0.753 ± 0.131
Please cite this article as: Wen Y, et al, A highly accurate symmetric optical flow based high-dimensional nonlinear spatial normalization of brain images, Magn Reson Imaging (2015), http://dx.doi.org/10.1016/j.mri.2015.01.013
Y. Wen et al. / Magnetic Resonance Imaging xxx (2015) xxx–xxx
References [1] Du J, Goh A, Qiu A. Large deformation diffeomorphic metric mapping of orientation distribution functions. Int. Conf. Information Processing in Medical Imaging, 22; 2011. p. 448–62. [2] Geng X, Ross TJ, Gu H, Shin W, Zhan W, Chao YP, et al. Diffeomorphic image registration of diffusion MRI using spherical harmonics. IEEE Trans Med Imaging 2011;30(3):747–58. [3] Duarte-Carvajalino JM, Sapiro G, Harel N, Lenglet C. A Framework for Linear and Non-Linear Registration of Diffusion-Weighted MRIs Using Angular Interpolation. Front Neurosci 2013:7–41. [4] Bookstein FL. Principal warps: Thin-plate splines and the decomposition of deformations. IEEE Trans Pattern Anal Machine Intell 1989;11:567–85. [5] Rohr K. “Image registration based on thin plate splines and local estimates of anisotropic landmark localization uncertainties”. In Proc. Conf. Medical Image Computing and Computer-Assisted Intervention (MICCAI’98); 1998. p. 1174–83. [6] Vaillant M, Davatzikos C. Hierarchical matching of cortical features for deformable brain image registration. Proc. Int. Conf. Information Processing in Medical Imaging; 1999. p. 182–95. [7] Chui H, Win L, Schultz R, Duncan J, Rangarajan A. A unified feature registration method for brain mapping. Proc. Int. Conf. Information Processing in Medical Imaging; 2001. p. 300–14. [8] Yap PT, Wu GR. TIMER: Tensor Image Morphing for Elastic Registration. NeuroImage 2009;47:549–63. [9] Yap PT, Wu GR, Zhu HT, Lin W, Shen DG. F-TIMER: Fast Tensor Image Morphing for Elastic Registration. IEEE Trans Med Imaging 2010;29(5):1192–203. [10] Shen D, Davatzikos C. HAMMER: heirarchical attribute matching mechanism for elastic registration. IEEE Trans Med Imaging 2002;21(11):1421–39. [11] Greve DN, Fishl B. Accurate and robust brain image alignment using boundarybased registration. NeuroImage 2009;48:63–72. [12] Yang JZ, Shen DG, Misra C, Xu XY, Resnick S, Davatzikos C, et al. “Spatial normalization of diffusion tensor images based on anisotropic segmentation”. Medical Imaging 2008: Image Processing 2008;6914:69140 L1–10. [13] Xue Z, Hai L, Lei G, Stephen TC. A local fast marching-based diffusion tensor image registration algorithm by simultaneously considering spatial deformation and tensor orientation. NeuroImage 2010;52:119–30. [14] Memin E, Perez P. Hierarchical estimation and segmentation of dense motion fields. Int J Comput Vis 2002;46(2):129–55. [15] Black MJ, Anandan P. Robust dynamic motion estimation over time. Proc. 1991 IEEE Computer Society Conference on Computer Vision and Pattern Recognition. IEEE Computer Society Press; 1991. p. 292–302. [16] Black MJ, Anandan P. The robust estimation of multiple motions: parametric and piecewise smooth flow fields. Comput Vis Image Underst 1996;63(1): 75–104. [17] Christensen GE, Rabbitt RD, Miller MI. Deformable templates using large deformation kinematics. IEEE Trans Image Process 1996;5(10):1435–47. [18] Horn B, Schunck B. Determining optical flow. Artif Intell 1981;17:185–203.
9
[19] Lucas B, Kanade T. An iterative image registration technique with an application to stereo vision. Proc. Seventh International Joint Conference on Artificial Intelligence; 1981. p. 674–9 [Vancouver, Canada]. [20] Xu L, Jia J, Matsushita Y. Motion detail preserving optical flow estimation. IEEE Trans Pattern Anal Mach Intell 2012;34:1744–57. [21] Stoll M, Volz S, Bruhn A. Adaptive integration of feature matches into variational optical flow methods. Computer Vision–ACCV 2012. Springer; 2013. p. 1–14. [22] Avants BB, Epstein CL, Grossman M, Gee JC. Symmetric diffeomorphic image registration with cross-correlation: evaluating automated labeling of elderly and neurodegenerative brain. Med Image Anal 2008;12:26–41. [23] Ashburner J. A fast diffeomorphic image registration algorithm. NeuroImage 2007;38:95–113. [24] Thirion JP. Non-rigid matching using demons. IEEE Computer Vision and Pattern Recognition; 1996. p. 245–51. [25] Fox Charles. An introduction to the calculus of variations. Courier Dover Publications; 1987. [26] David M. Iterative methods for solving partial difference equations of elliptical type. (PhD thesis) Harvard University; 1950. [27] Hill S Francis. Computer Graphics. Prentice Hall; 2001. [28] Black MJ, Anandan P. “Robust dynamic motion estimation over time”. In Proc 1991 IEEE Computer Society Conference on Computer Vision and Patter Recognition; 1991. p. 292–302. [29] Yin W, Osher S, Goldfarb D, Darbon J. Bregman Iterative Algorithms for L1Minimization with Applications to Compressed Sensing, SIAM J. Imaging Sciences 2008;1(1):143–68. [30] Sun D, Roth S, Black M. A Quantitative Analysis of Current Practices in Optical Flow Estimation and the Principles Behind Them. Int J Comput Vis 2013:1–23. [31] Álvarez L, Castano CA, García M, Krissian K, Mazorra L, Salgado A, et al. Symmetric optical flow. Computer Aided Systems Theory–EUROCAST 2007. Springer; 2007. p. 676–83. [32] Alvarez L, Deriche R, Papadopoulo T, Sánchez J. Symmetrical dense optical flow estimation with occlusions detection. Int J Comput Vis 2007;75:371–85. [33] Yoo T. Insight Into Images: Principles and Practice for Segmentation, Registration and Image Analysis. Natick, MA: A.K. Peters Ltd; 2003. [34] Klein A, Andersson J, Ardekani BA, Ashburner J, et al. Evaluation of 14 nonlinear deformation algorithms applied to human brain MRI registration. NeuroImage 2009;46:786–802. [35] Smith SM, Jenkinson M, Woolrich MW, Beckmann CF, Behrens TEJ, JohansenBerg H, et al. Advances in functional and structural MR image analysis and implementation as FSL. NeuroImage 2004;23(S1):208–19. [36] Jenkinson M, Beckmann CF, Behrens TE, Woolrich MW, Smith SM. FSL. NeuroImage 2012;62:782–90. [37] Jenkinson M, Smith SM. Bayesian analysis of neuroimaging data in FSL. NeuroImage 2009;45:S173–86. [38] Wu G, Kim M, Wang Q, Shen D. S-HAMMER: Hierarchical attribute‐guided, symmetric diffeomorphic registration for MR brain images. Hum Brain Mapp 2014;35(3):1044–60. [39] Chang HH, Zhuang AH, Valentino DJ, Chu WC. Performance measure characterization for evaluating neuroimage segmentation algorithms. NeuroImage 2009; 47:122–35.
Please cite this article as: Wen Y, et al, A highly accurate symmetric optical flow based high-dimensional nonlinear spatial normalization of brain images, Magn Reson Imaging (2015), http://dx.doi.org/10.1016/j.mri.2015.01.013
Contents lists available at ScienceDirect
Magnetic Resonance Imaging journal homepage: www.mrijournal.com
A highly accurate symmetric optical flow based high-dimensional nonlinear spatial normalization of brain images☆ Ying Wen a,⁎, Lili Hou a, Lianghua He b, Bradley S. Peterson c, Dongrong Xu d,⁎ a
Shanghai Key Laboratory of Multidimensional Information Processing & Department of Computer Science and Technology, East China Normal University, Shanghai, 200241, China Key Laboratory of Embedded System and Service Computing, Ministry of Education, Tongji University, Shanghai 200092, China Institute for the Developing Mind, Children's Hospital Los Angeles, The Keck School of Medicine at the University of Southern California, Los Angeles, CA 90027, U.S.A. d MRI Unit & Epidemiology Division, Department of Psychiatry, Columbia University & New York State Psychiatric Institute, New York, 10032, U.S.A. b c
a r t i c l e
i n f o
Article history: Received 9 May 2014 Revised 21 October 2014 Accepted 18 January 2015 Available online xxxx Keywords: Brain image normalization Symmetric optical flow Hierarchical strategy Optical flow constraints
a b s t r a c t Spatial normalization plays a key role in voxel-based analyses of brain images. We propose a highly accurate algorithm for high-dimensional spatial normalization of brain images based on the technique of symmetric optical flow. We first construct a three dimension optical model with the consistency assumption of intensity and consistency of the gradient of intensity under a constraint of discontinuity-preserving spatio-temporal smoothness. Then, an efficient inverse consistency optical flow is proposed with aims of higher registration accuracy, where the flow is naturally symmetric. By employing a hierarchical strategy ranging from coarse to fine scales of resolution and a method of Euler-Lagrange numerical analysis, our algorithm is capable of registering brain images data. Experiments using both simulated and real datasets demonstrated that the accuracy of our algorithm is not only better than that of those traditional optical flow algorithms, but also comparable to other registration methods used extensively in the medical imaging community. Moreover, our registration algorithm is fully automated, requiring a very limited number of parameters and no manual intervention. © 2015 Elsevier Inc. All rights reserved.
1. Introduction Many methods in image analysis have been developed to tackle diffusion tensor images (DTI) registration [1–3]. DTI registration generally falls into two categories based on setting up the correspondence between two subjects. The first involves feature-based matching, i.e. transformations that are calculated based on a number of anatomical correspondences established manually, semi-automatically, or fully automatically on a number of distinct anatomical features, including distinct landmarks [4,5] or a combination of curves and surfaces, such as sulci and gyri [6,7]. In addition, finding anatomical correspondences is a key to the success of the spatial normalization. To find anatomical correspondences, conventional methods generally first extract scalar features from each tensor individually; and then by constructing scalar maps, regional integration and other operations such as edge detection, a final set of features for correspondence matching can be constituted. For example, Timer [8] and F-Timer [9] were developed ☆ This work was supported by The National Natural Science Foundation of China (No. 61273261, 61272267, 81471734), Program for New Century Excellent Talents in University (NCET-11-0381), Shanghai Collaborative Innovation Center of Trustworthy Software for Internet of Things (ZF1213), NIDA 1R01DA027100 and NIMH 5R01MH089582. ⁎ Corresponding authors. E-mail addresses: [email protected] (Y. Wen), [email protected] (B.S. Peterson), [email protected] (D. Xu).
based on Hammer [10], which adopts so-called driving voxels in addition to geometric moments as the features upon which registration is based. The boundary-based registration uses the measure of the thickness of cortex in T1-weighted image [11]. Yang et al. guided their registration using information on a tensor’s structural geometry and orientation [12]. Xue et al. proposed a local fast marching method for DTI registration [13]. These landmark-based methods numerically seek optimal transformations (under respective model-specific assumptions) to maximize similarity across brains. The second category of registration methods is based on volumetric transformations between one brain image and a template image. These methods assume that the images are acquired using the same pulse sequence [14–16] and these methods therefore are generally more efficient, because they do not require the construction of a specific anatomical model using anatomic landmarks. Christensen et al. [17] proposed a viscous fluid flow model to implement registration. Unlike fluid flow, optical flow is another volumetric transformation method. Optical flow is usually posed as an energy minimization problem, with the energy function containing a data term and a smoothness term. The variational framework, together with coarse-to-fine refinement, is widely used in optical flow estimation. The technique of optical flow developed in computer vision for object tracking [18,19] is one such method. The two approaches differ with the constraints in their implementations. The former adopts viscous constraints and body
http://dx.doi.org/10.1016/j.mri.2015.01.013 0730-725X/© 2015 Elsevier Inc. All rights reserved.
Please cite this article as: Wen Y, et al, A highly accurate symmetric optical flow based high-dimensional nonlinear spatial normalization of brain images, Magn Reson Imaging (2015), http://dx.doi.org/10.1016/j.mri.2015.01.013
2
Y. Wen et al. / Magnetic Resonance Imaging xxx (2015) xxx–xxx
Fig. 1. Fractional anisotropy maps of the 20 randomly selected subjects from our database. The brains vary significantly in morphometry.
force while the latter adopts gradient and smooth constants. For the solution of optical flow, Xu et al. proposed an extended coarse-to-fine refinement framework which reduces the reliance of flow estimation on coarse level and corrects the optical flow in each scale [20]. Stoll et al. proposed an adaptive integration strategy for descriptor matching by selecting the positions where descriptor matching may improve the optical flow estimation [21]. Besides optical flow method, diffeomorphism methods are often used for image registration. The registration methods used extensively in the medical imaging community include symmetric diffeomorphic image registration (SyN) [22], SPM5 DARTEL [23], as well as with Thirion’s Demons algorithm [24] etc. The optimization of diffeomorphism method is implemented on a space of diffeomorphic space instead of the entire space of displacement field. SyN maximizes the cross-correlation in the space of diffeomorphic maps. DARTEL is a diffeomorphic image registration for estimating inverse consistent deformations. Demons algorithm uses an approximate elastic regularizer to solve an optical flow problem. Traditional techniques of optical flow estimation do not generally yield symmetrical solutions: the results will differ if they are applied between images f1 and f2 or between images f2 and f1. Ideally, the forward and inverse flows should be uniquely determined and should be inverses of one another. Estimating the forward and inverse flows independently very rarely results in a consistent set of transformations due to a large number of local minima. To overcome this deficiency in current optical flow, symmetrical formulations of the optical flow has been proposed in Refs. [29–31], where the solution is a constraint to being symmetric using a combination of the flow in both directions. In this work, we present a simple and
effective method to recover a dense optical flow field map from two images. The idea is that it exists as a symmetric displacement between displacement vectors from f1 to f2 and f2 to f1 and to minimize an energy function of optical flow. This variational problem is then solved using the gradient flow defined by the Euler–Lagrange equations associated to the energy. This idea is introduced for DTI brain registration in this paper. We focus on the optical flow based high-dimensional nonlinear spatial normalization of brain images. We first extend the traditional 2D optic flow to 3D form for high-dimensional nonlinear spatial normalization. Then, we propose a simple and effective symmetric variational formulation of the optical flow, where the flow is naturally symmetric for the problem of the forward and inverse flows. Finally, we incorporate the symmetric optical model into a hierarchical framework that is applied to imaging data in multiple resolution scales from coarse to fine, thereby building up a registration approach with high accuracy. 2. Materials and methods 2.1. Data acquisition We acquired all the human Diffusion-weighted imaging (DWI) datasets from a GE 3.0 Tesla MR scanner system (General Electric, Milwaukee, Wisconsin) at New York State Psychiatry Institute, Columbia University. This unit was equipped with a high strength (50 mT/m) and high-speed gradient system (slew rate = 200 T/m/s) capable of conducting DWI. A single-shot spin echo, echo planar imaging (SE–EPI) sequence was used to acquire the DWI data, with
Please cite this article as: Wen Y, et al, A highly accurate symmetric optical flow based high-dimensional nonlinear spatial normalization of brain images, Magn Reson Imaging (2015), http://dx.doi.org/10.1016/j.mri.2015.01.013
Y. Wen et al. / Magnetic Resonance Imaging xxx (2015) xxx–xxx
diffusion gradients applied at b = 1000 s/mm2 along 25 non-collinear spatial directions plus 3 baseline images without applying diffusion gradient. The thickness of each slice was 2.5 mm without gap, and about 68 axial slices were assigned to cover the entire brain. The other parameters of the DWI sequence were minimal echo time (TE) at about 77.3 ms; repetition time (TR) = 17000 ms; field of view (FOV) = 240 mm × 240 mm and acquisition data matrix =132 × 128. The final images were machine-interpolated to an in-plane resolution of 256 × 256 voxels, resulting in a nominal in-plane resolution of 0.9375 mm × 0.9375 mm. The total DWI scanning time was about 16 minutes 26 seconds with the number of excitation (NEX) = 2. Twenty data sets of healthy volunteers were randomly selected from our database. These brains vary significantly in morphology (Fig. 1). 2.2. Methods Our algorithm generally follows the conventional framework for image registration, i.e. rigid global transformation plus nonlinear registration. The global transformation roughly aligns an individual image to a selected template image in terms of its size and spatial orientation. By extending the standard 2D optic flow method to a 3D form, we propose in this paper a registration method based on the 3D optical flow method for local nonlinear registration of the DTI data which is high-dimensional and incorporates an energy function based on the three mentioned assumptions. In this process, Euler-Lagrange equation [25] is employed to optimize the energy function, and the successive over-relaxation method (SOR) [26] is used to search for the optimal solution. This framework calculates the deformation filed for normalizing DTI data using diffusion anisotropy index (DAI) images that reflect the underlying diffusion property. The same framework then normalizes the DTI data using this calculated deformation field that correlates the spatial correspondence between the two spaces. The DAI used in our case is the fractional anisotropy (FA) map.
minimizing the energy based on this L1 norm requires that the function is a convex. Because in this form ε is a small positive constant, ψ(e 2)remains a convex which meets the mentioned requirement in the process of searching minimum. Moreover, this form of ψ(e 2) does not introduce any additional parameters, because ε is only for numerical calculation and can be fixed to a constant value (0.001 in our case). Finally, we incorporate a smoothness term into the energy function to satisfy the assumption that the flow field must be piecewise smooth. The smoothness term is necessary to cope with the function’s ill-posedness and noise which introduce destructive effects in the optical flow approach. In the optical flow equation, the variation of the pixels’ intensity is interpreted as their motion. Without the motion smoothness term, the estimated motion of adjacent pixels can be completely disparate or significantly deviate from the truth. The smoothness constraint is achieved by penalizing the total variation of the flow field, which can be expressed as: Z Esmooth ðu; v; wÞ ¼
We now present our nonlinear registration method by deducting the equations using the 3D form of the optic flow. Suppose that s := (x, y, z, t) T is the voxel in the image at location (x, y, z) at time t and p := (u, v, w, 1) T is the deformation vector field, the data term is: 2
2
j f ðs þ pÞ− f ðsÞj þ βj∇f ðs þ pÞ−∇f ðsÞj
ds
ð1Þ
Because quadratic penalisers may generate outliers that may excessively bias the estimation of optic flow, we introduce into this energy function an increasing concave function ψ(e2), revising it into a more robust energy function [28], although the function ψ can also be applied separately to each of these terms: Z Eðu; v; wÞ ¼
S
2
Z Etotal ðu; v; wÞ ¼
S
2
ð3Þ
2
2
ψðj f ðs þ pÞ− f ðsÞj þ βj∇f ðs þ pÞ−∇ f ðsÞj Þ 2
2.2.2. Three dimensional optical flow algorithm
S
2
ψðj∇uj þ j∇vj þ j∇wj Þds
2
ð4Þ
2
þλψðj∇uj þ j∇vj þ j∇wj Þds
The global transformation is a rigid affine transformation. In general, an affine transformation is composed of linear transformations (rotation, scaling or shear) and a translation (or “shift”). The brain in its native image space is transformed to the space of the template using the global affine transformation [27], as a proper initialization for the subsequent coregistration process.
Z
S
using the same function for Ψ as above. The spatio-temporal gradient ∇ := (∂x, ∂y, ∂z, ∂t)′ incorporates the needed spatio-temporal smoothness assumption. For coregistration that involves only two images, this gradient term should be replaced by one that involves only the spatial gradient ∇ := (∂x, ∂y, ∂z)′. Thus, the total energy is the weighted sum of (2) and (3):
2.2.1. Global transformation
Eðu; v; wÞ ¼
3
2 2 ψ j f ðs þ pÞ− f ðsÞj þ βj∇f ðs þ pÞ−∇f ðsÞj ds
ð2Þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi We adopt the form ψ e2 ¼ e2 þ ε2 in this energy calculation, which results in (modified) L1 minimization [29]. However,
where, λ is a positive regularization parameter. Now we need to solve this equation for searching the best solution. 2.2.3. The proposed optical flow model In the standard formulation of variational optical flow, the estimated motion vector fields depend on the template image and are asymmetric. However, in most application the solution should be independent of the template image. Motivated by [30–32], we use a symmetric optical flow to settle the problem of consistency. For image fm between first image f1 and second image f2, there exist an optical flow p that satisfies:
f m ðsÞ ¼ f ðs−p=2Þ ¼ f ðs þ p=2Þ
ð5Þ
Fig. 2 shows a symmetric flow description. The images in the first row of Fig. 2 display the flow from f1 to f2 in the traditional way, and the images in the second row display the flow in a symmetric way. Then, (4) can be written as 1 2 ψðj f ðs þ p=2Þ−f ðs−p=2Þj Þþ B C Eðu; v; wÞ ¼ @ βψðj∇f ðs þ p=2Þ−∇ f ðs−p=2Þj2 Þþ Ads Ω 2 2 2 λψðj∇uj þ j∇vj þ j∇wj Þ Z
0
ð6Þ
Please cite this article as: Wen Y, et al, A highly accurate symmetric optical flow based high-dimensional nonlinear spatial normalization of brain images, Magn Reson Imaging (2015), http://dx.doi.org/10.1016/j.mri.2015.01.013
4
Y. Wen et al. / Magnetic Resonance Imaging xxx (2015) xxx–xxx
We name ft the difference image between the two consecutive image frames; fi (i = x, y, z) is the first order derivative of the movement image f(s + p); fij(i,j = x, y, z) is the second order derivative of the movement image and fit (i = x, y, z) is the first order of the difference image. To solve (7) for the variables u, v, w, we follow the construction of the numerical scheme performed in [12]. The first step is to introduce the fixed point iteration on p. Let p k = (u k, v k, w k,1) T the solution of the system, i.e. the deformation field, at the kth iteration step and let f⁎⁎k be the abbreviations defined in (8), an implicit scheme for the smoothness term and a semi-implicit scheme for the data term have been applied. This gives the following system of equations at the (k + 1)th iteration step: Fig. 2. An explanation of symmetric optical flow.
2.2.4. Numerical solution to the energy function Using the Euler–Lagrange equations, the energy function (6)can be transformed into the following form:
8 0 2 2 2 2 > > ψ f t þ β f xt þ f yt þ f zt > > > > > > f x f t þ β f xx f xt þ f xy f yt þ f xz f zt > > > > > 0 2 2 2 > > −λdiv ψ j∇uj þ j∇vj þ j∇wj ∇u ¼ 0; > > > > 0 2 2 2 2 > > ψ f t þ β f xt þ f yt þ f zt > > < f y f t þ β f xy f xt þ f yy f yt þ f yz f zt > > > 0 2 2 2 > > −λdiv ψ j∇uj þ j∇vj þ j∇wj ∇v ¼ 0; > > > > 0 2 2 2 2 > > ψ f t þ β f xt þ f yt þ f zt > > > > > > > ⋅ f z f t þ β f xz f xt þ f yz f yt þ f zz f zt > > > > > : −λdiv ψ0 j∇uj2 þ j∇vj2 þ j∇wj2 ∇w ¼ 0
8 > 0 kþ1 2 kþ1 2 kþ1 2 kþ1 2 > Þ þ β f xt þ f yt þ f zt ψ ft > > > > > > k kþ1 k kþ1 k kþ1 k kþ1 > > ⋅ f x f t þ β f xx f xt þ f xy f yt þ f xz f zt > > > > > kþ1 2 > 0 kþ1 2 kþ1 2 kþ1 > ¼ 0; −λdiv ψ ∇u þ ∇v þ ∇w ∇u > > > > > > > 0 kþ1 2 kþ1 2 kþ1 2 kþ1 2 > > ft þ β f xt þ f yt þ f zt > ψ > > < k kþ1 k kþ1 k kþ1 k kþ1 ⋅ f y f t þ β f xy f xt þ f yy f yt þ f yz f zt > > > 2 2 > kþ1 0 kþ1 kþ1 2 kþ1 > > ¼ 0; −λdiv ψ ∇u þ ∇v þ ∇w ∇v > > > > 2 2 2 2 > > 0 kþ1 kþ1 kþ1 kþ1 > > ft þ β f xt þ f yt þ f zt > > ψ > > > > > ⋅ f k f kþ1 þ β f k f kþ1 þ f k f kþ1 þ f k f kþ1 > > z t xz xt yz yt zz zt > > > > kþ1 2 0 kþ1 2 kþ1 2 kþ1 > > ¼0 : −λdiv ψ ∇u þ ∇v þ ∇w ∇w
ð9Þ
ð7Þ This system is still nonlinear because of the nonlinear function ψ′ and the symbols f⁎k + 1. To facilitate the subsequent calculation, we apply the first order Taylor expansions to remove the nonlinearity from the components f⁎k + 1: kþ1 k k k ≈f s þ p þ fx s þ p du f s þ p k k k k þfy s þ p dv þ fz s þ p dw
where we use the following abbreviation:
∂ f ðs þ p=2Þ þ ∂x f ðs−p=2Þ fx ¼ x 2 ∂y f ðs þ p=2Þ þ ∂y f ðs−p=2Þ fy ¼ 2 ∂ f ðs þ p=2Þ þ ∂z f ðs−p=2Þ fz ¼ z 2 f t ¼ f ðs þ p=2Þ−f ðs−p=2Þ ∂ f ðs þ p=2Þ þ ∂xx f ðs−p=2Þ f xx ¼ xx 2 ∂yy f ðs þ p=2Þ þ ∂yy f ðs−p=2Þ f yy ¼ 2 ∂ f ðs þ p=2Þ þ ∂zz f ðs−p=2Þ f zz ¼ zz 2 ∂xy f ðs þ p=2Þ þ ∂xy f ðs−p=2Þ f xy ¼ 2 ∂ f ðs þ p=2Þ þ ∂xz f ðs−p=2Þ f xz ¼ xz 2 ∂yz f ðs þ p=2Þ þ ∂yz f ðs−p=2Þ f yz ¼ 2 f xt ¼ ∂x f ðs þ p=2Þ−∂x f ðs−p=2Þ f yt ¼ ∂y f ðs þ p=2Þ−∂y f ðs−p=2Þ f zt ¼ ∂z f ðs þ p=2Þ−∂z f ðs−p=2Þ
The unknown iteration variable p k + 1 is thus substituted by the known variable p k and an unknown term dp k = (du k, dv k, dw k, 1) T. Therefore, the derivative of the data term and the derivative of the smoothness term in (7) become:
ð8Þ
0 k 0 k k k k k k k 2 f t þ f xt du þ f yt dv þ f zt dw ψ D :¼ ψ k k k k k k k 2 þβ f xt þ f xx du þ f xy dv þ f xz dw k k k k k k k 2 þ f yt þ f xy du þ f yy dv þ f yz dw k k k k k k k 2 ; þ f zt þ f xz du þ f yz dv þ f zz dw 2 0 k 0 k k ψ S :¼ ψ ∇ u þ du 2 k k k k 2 þ ∇ v þ dv þ ∇ w þ dw
ð10Þ
where (ψ ')Dk can be interpreted as a robustness factor in the data term, and (ψ ')Sk as diffusivity in the smoothness term.
Please cite this article as: Wen Y, et al, A highly accurate symmetric optical flow based high-dimensional nonlinear spatial normalization of brain images, Magn Reson Imaging (2015), http://dx.doi.org/10.1016/j.mri.2015.01.013
Y. Wen et al. / Magnetic Resonance Imaging xxx (2015) xxx–xxx
Subsequently, (9) can be written as 8 0 k k k k k k k k k > ψ D f x f t þ f xt du þ f yt dv þ f zt dw > > > > > k k k k k k k k > > > þβ f xx f xt þ f xx du þ f xy dv þ f xz dw > > > > > þf kxy f kyt þ f kxy duk þ f kyy dvk þ f kyz dwk > > > > k k > k k k k k k > > þf xz f zt þ f xz du þ f yz dv þ f zz dw > > > > > −λdiv ψ0 k ∇ uk þ duk ¼ 0; > > S > > > 0 k k k > k k k k k k > ψ f > y f t þ f xt du þ f yt dv þ f zt dw > D > > > > > þβ f kxy f kxt þ f kxx duk þ f kxy dvk þ f kxz dwk > > < k k k k k k k k þf yy f yt þ f xy du þ f yy dv þ f yz dw > > > k k k k k k k k > > þf yz f zt þ f xz du þ f yz dv þ f zz dw > > > > 0 k k k > > ∇ v þ dv ¼ 0; > S > −λdiv ψ > > 0 k > k k k k k k k k > > ψ D f z f t þ f xt du þ f yt dv þ f zt dw > > > > k k k k k k k k > > þβ f xz f xt þ f xx du þ f xy dv þ f xz dw > > > > k k > k k k k k k > þf f þ f xy du þ f yy dv þ f yz dw > > > yz yt > > k k k k k k k k > > þf zz f zt þ f xz du þ f yz dv þ f zz dw > > > > > : −λdiv ψ0 k ∇ wk þ dwk ¼ 0 S
multiresolution approach captures large-scale deformation in the lower-resolution steps and small deformation in the finer resolution steps. After a designated number of iterations are carried out at each resolution level, or when convergence is observed, the estimated motion field is passed to the next level of a super sampled resolution using cubic interpolation, which serves as the initialization for the new estimation at the finer level of image resolution. 3. Experiments
ð11Þ
However, (11) still contains nonlinear components ψ' D and ψ' S. In order to remove the remaining nonlinearity in ψ ', a second, inner, fixed point iteration loop is applied. Let du k,n and dv k,n denote the iteration variables at step n; and let (ψ ')Dk,n and (ψ ')Sk,n denote the robustness factor and the diffusivity defined in (11) at iteration k and n (the iteration regarding k is for removing the linearity in f⁎k + 1 and the iteration regarding n is for removing the linearity in ψ' D and ψ' S). Then, the final linear system of equations in du k,n + 1 and dv k,n + 1 reads 8 0 k;n k k k k;nþ1 k k;nþ1 k k;nþ1 > ψ D f x f t þ f xt du þ f yt dv þ f zt dw > > > > > k k k k;nþ1 k k;nþ1 k k;nþ1 > > þ f xy dv þ f xz dw > þβ f xx f xt þ f xx du > > > > > þf kxy f kyt þ f kxy duk;nþ1 þ f kyy dvk;nþ1 þ f kyz dwk;nþ1 > > > > k k > k k;nþ1 k k;nþ1 k k;nþ1 > > þf xz f zt þ f xz du þ f yz dv þ f zz dw > > > > > −λdiv ψ0 k;n ∇ uk þ duk;nþ1 ¼ 0; > > S > > > 0 k;n k k > k k;nþ1 k k;nþ1 k k;nþ1 > f þ f yt dv þ f zt dw ψ > y f t þ f xt du > D > > > > > þβ f kxy f kxt þ f kxx duk;nþ1 þ f kxy dvk;nþ1 þ f kxz dwk;nþ1 > > < k k k k;nþ1 k k;nþ1 k k;nþ1 þf yy f yt þ f xy du þ f yy dv þ f yz dw > > > k k k k;nþ1 k k;nþ1 k k;nþ1 > > þf yz f zt þ f xz du þ f yz dv þ f zz dw > > > > 0 k;n k k;nþ1 > > ¼ 0; −λdiv ψ S ∇ v þ dv > > > > 0 k;n k k > k k;nþ1 k k;nþ1 k k;nþ1 > > ψ D f z f t þ f xt du þ f yt dv þ f zt dw > > > > k k k k;nþ1 k k;nþ1 k k;nþ1 > > þβ f xz f xt þ f xx du þ f xy dv þ f xz dw > > > > > k k k k;nþ1 k k;nþ1 k k;nþ1 > þf yz f yt þ f xy du þ f yy dv þ f yz dw > > > > > k k k k;nþ1 k k;nþ1 k k;nþ1 > > þf zz f zt þ f xz du þ f yz dv þ f zz dw > > > > > : −λdiv ψ0 k;n ∇ wk þ dwk;nþ1 ¼ 0 S
5
ð12Þ
Eq. (12) is actually in the form of A(du, dv, dw) T = B. Using standard discretization for the derivatives, the resulting sparse linear system of equations can now be solved using common numerical methods, such as Gauss-Seidel or SOR iterations. We employed a hierarchical multiresolution scheme to accommodate different scales of motion for improving computational efficiency. The multiresolution approach speeds up the convergence of the iteration represented by (12). We also expect that the
We conducted three different experiments to evaluate from various aspects the performance of our algorithm. The first experiment examined the accuracy of the deformation fields generated by our algorithm using simulated datasets. The second experiment demonstrated the excellent performance of our algorithm using real-world data acquired directly from human participants. The third experiment evaluated the performance of our algorithm by conducting a comparison with other algorithms using fiber tractography. All our experiments were conducted using the following empirical parameters λ = 40; β = 60; number of SOR iterations = 500, image dimension: 256 × 256 × Z (Z varied from 50 to 65 slices), resolution = 0.9375 mm × 0.9375 mm × 2.5 mm. Datasets with motion artifacts and outliers were not included. 3.1. Simulated experiments In the experiment, we aimed to evaluate quantitatively the performance of our algorithm using simulated data in terms of the morphological accuracy of spatial normalization. We inspected the registration error of our algorithm for coregistering two brain images for their structural details, and therefore focused on the morphological information. The distinct features relevant particularly to tensor data were thus temporarily not considered, which would be the focus of the subsequent experiments. To estimate the proposed registration model accuracy, we adopted a strategy similar to that was reported in [19]. We first randomly selected from our DTI database of brain imaging data six morphological representative datasets that were acquired from 6 healthy human participants. These 6 datasets differed from each other in size, shape and structure. We then randomly chose one of the 6 as the template brain, and calculated the FA maps for all the 6 datasets. Next, we warped the FA of the template brain to each FA of the other 5 datasets individually, using SPM5 DARTEL [22]. We thus obtained 5 versions of the template brain in 5 different individual spaces (Fig. 3), along with 5 deformation fields DFT2I that maps the correspondence between the template space and the individual brains. Although the 5 configurations of the template brain differed with each other in size and shape in different spaces, they shared exactly the same structure as these 5 were basically the same brain. A good algorithm for spatial normalization should be able to accurately bring them back to the original template space with perfect registration to the original template brain. We applied our algorithm thus to test the morphological accuracy that our algorithm could offer. We compared the performance of our optical flow based algorithm with several other algorithms including the affine transform, the SPM5 DARTEL, the Horn-Schunck’s algorithm [15] and the Lucas-kanade algorithm [16]. In the template brain, we manually marked 200 landmarks on those visually distinct structural features. The deformation fields DFT2I mapped them to 1000 locations (200 each for 5 brains) in the individual spaces; and the deformation fields DFI2T generated by our algorithm took all the 5 deformed brains back to the template space, thus mapping these 1000 locations to 5 new sets of 200 locations in the original template space. By measuring the average distance between the 1000 new
Please cite this article as: Wen Y, et al, A highly accurate symmetric optical flow based high-dimensional nonlinear spatial normalization of brain images, Magn Reson Imaging (2015), http://dx.doi.org/10.1016/j.mri.2015.01.013
6
Y. Wen et al. / Magnetic Resonance Imaging xxx (2015) xxx–xxx
(a)
(b1)
(b2)
(b3)
(b4)
(b5)
Fig. 3. Images of the synthetic brains for inspecting the registration performance of our method. (a) The template brain. (b1–b5) Five deformed versions of the template brain in 5 individual spaces.
locations and their corresponding landmarks we originally plotted, we would be able to evaluate the registration error of using our algorithm. Suppose P(x, y, z) are the predefined landmarks in the template brain, P ' (x ', y ', z ') the corresponding locations in the subject spaces mapped by DFT2I, and P 00 (x 00 , y 00 , z 00 ) the final locations back to the template space generated by DFI2T, the morphological accuracy of the normalization algorithm was measured by the mean and the standard deviation of distance between P(x, y, z) and P 00 (x 00 , y 00 , z 00 ) as follows: md ¼
1 X
P ðx; y; zÞ−P 00 x00 ; y00 ; z00
jV j x;y;z∈V
ð13Þ
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
1 X
P ðx; y; zÞ−P 00 x00 ; y00 ; z00 −m 2 d jV j x;y;z∈V
ð14Þ
and σd ¼
where ‖⋅‖ is the Euclidean distance, |V| the number of landmarks we defined in the template brain. To provide a more general evaluation of the registration error, we also calculated the mean and standard deviation of the distance for every brain voxel in the whole image and compared the results generated by all the methods under comparison. The quantitative report indicated that the registration error of our method was only 0.85 mm (mean) and 0.88 mm (std. var) for the 1000 anatomical landmarks, and 0.98 mm and 0.94 mm for all the brain voxels (Table 1), both smaller than the size of one voxel. In contrast, among all the other four methods except ours (DARTEL, affine registration, the Horn-Schunck method and the Lucas-kanade method) under comparison, DARTEL performed the best but the errors were significantly larger than ours (Table 1). The results demonstrated that our method performed superior to the other methods in coregistering the image for the structural details. 3.2. Experiments on DTI data
is to inspect the sharpness of an average image obtained from a set of coregistered images. We randomly selected from our DTI database 20 healthy brains for this experiment. These brains vary significantly in morphology. One of the 20 brains was then randomly selected as the template for all the rest 19 to be spatially normalized. We applied also the affine registration and DARTEL registration to compare with the performance of our algorithm. After all datasets were normalized into the template space, we averaged them. In addition to visual inspection of the sharpness of the averaged data, we also performed quantitatively evaluation of the sharpness for the whole 3D image by measuring the standard deviation of the FA value from the reference in the template at each voxel and generated a standard deviation map for each method that was under comparison. The mean image of the registered brains showed that compared with the results generated by DARTEL and the affine registration, our method generated the clearest structural boundaries in the results (Fig. 4), indicating that our method performed the best coregisteration of the individual images. Our method not only registered nearly perfectly the general morphology of the brains, it also registered very well the small structures inside of the brain. In the following experiment, the tracking results appeared that the tracked fibers based on the data normalized using our method resampled the ground truth much better than did the results based on the data normalized using the DARTEL method (Fig. 5). We obtained 693 fibers passing through the splenium ROI and 512 fibers passing through the genu ROI of a subject shown in Fig. 5(b). After the data normalized using DARTEL and our method, we obtained 699 and 689 fibers passing through the splenium ROI and 348 and 490 fibers passing through the genu ROI shown in Fig. 5(c)
Table 1 A comparison of the registration errors using different methods (mm). Measurement
In this subsection, we wanted to evaluate the morphological accuracy of our algorithm in practice directly using real-world DTI data acquired from human participants. If a registration method is not accurate, the structural boundary of the average image of a number of registered images will definitely be fuzzy. Therefore, a direct and intuitive way of gauging registration error of an algorithm
Horn-Schunck Lucas-kanade Affine registration DARTEL Our method
Landmarks
All voxels
Mean
Std.
Mean
Std.
2.47 3.14 4.38 1.32 0.85
2.58 2.66 3.02 0.98 0.88
3.12 3.48 4.65 1.45 0.98
3.37 3.54 3.78 1.24 0.94
Please cite this article as: Wen Y, et al, A highly accurate symmetric optical flow based high-dimensional nonlinear spatial normalization of brain images, Magn Reson Imaging (2015), http://dx.doi.org/10.1016/j.mri.2015.01.013
Y. Wen et al. / Magnetic Resonance Imaging xxx (2015) xxx–xxx
(a)
(b)
(c)
7
(d)
Fig. 4. A comparison of the performance between our method and other approaches, using 20 human brain images. (a) The template brain. (b) The mean brain of the registration results generated by the affine registration. (c) The mean brain of the registration results generated by DARTEL registration. (d) The mean brain of the registration results generated by our method.
and (d), respectively. Compared the result based on DARTEL shown in Fig. 5(c), the topology of the result based on our approach as shown in Fig. 5(d) better resembled the ground truth in Fig. 5(b). For example, the y-shaped fiber branches (pointed by the white arrows) were well reserved using our method but disappeared in the data generated by DARTEL; moreover, DARTEL generated lots of originally non-existing fiber pathways (pointed by red arrows). In particular, the data processed using our method almost duplicated the fiber structures existing in the original data. In
(a)
(b)
contrast, the data using DARTEL generated lots of non-existing fiber pathways and destroyed lots that originally existed. Similarly, the quantitative analysis again demonstrated that the results based on our method approximated the ground truth quite well, while the results based on SPM5 DARTEL had significant difference from the ground truth (Table 2). In general, the statistics showed that the results by DARTEL had a much larger standard deviation than those using our method. To validate the registration performance of our method, we compared the proposed method with greyscale image based
(c)
(d)
Fig. 5. The results of fiber tracking, based on the two regions-of-interest using the diffusion tensor imaging data normalized using either DARTEL or our approach. (a) A triplanar view (viewed from top left of the brain) of the tracked fibers in the original template brain that passed through the two ROIs, viewed from above the right ear. The fibers were superimposed on the corresponding FA image. (b)The same as that shown in (b), but only the tracked fibers are displayed here, and the view was generated from top of the brain. (c)The tracked fibers in the data normalized using DARTEL. (d) The tracked fibers in the data normalized using our method.
Please cite this article as: Wen Y, et al, A highly accurate symmetric optical flow based high-dimensional nonlinear spatial normalization of brain images, Magn Reson Imaging (2015), http://dx.doi.org/10.1016/j.mri.2015.01.013
8
Y. Wen et al. / Magnetic Resonance Imaging xxx (2015) xxx–xxx
Table 2 A comparison of the statistics of the fiber tracking results based on the two seeding volumes defined on the genu and the splenium of the corpus callosum (The statistics in parentheses denotes the difference between the testing item and the original template). Object
Template
Fibers_bundle
Splenium
Genu
DARTEL Splenium
Genu
Splenium
Genu
Mean fiber length Std of length Mean fiber angle Std of angle Mean fiber count Std of count Mean volume Std of volume
49.84 / 3.09 / 693 / 14604 /
32.05 / 4.88 / 512 / 5384 /
44.69 1.76 3.06 0.025 699 18.19 14 760 262.21
31.86 1.54 4.98 0.21 348 27.65 4182 182.62
49.81 1.37 3.09 0.02 689 12.91 14 512 202.63
31.38 1.57 4.84 0.06 490 11.88 5039 152.36
methods, such as symmetric diffeomorphic image registration (SyN) [21], SPM5 DARTEL [22], as well as with Thirion’s Demons algorithm [23] used extensively in the medical imaging community. SyN algorithm was implemented with the extended version of the ITK deformable image registration frame, described in Yoo [33]. Thirion’s Demons algorithm was known to perform well in inter-subject deformable image registration and used an approximate elastic regularizer to solve an optical flow problem. Demons algorithm was freely available in the Insight ToolKit and has been optimized by the ITK community (www.itk.org). We considered the above three registration methods for three reasons: (1) we wanted to evaluate brain registration algorithms, not brain labeling algorithms or particular atlas-based approaches, (2) these methods are based on greyscale image but not feature points, and (3) their performances are good according to Klein’s report [34]. In this experiment, each of the 20 brain images was registered to the 19 others, resulting in 380 pairwise registrations. Each of the brains is represented 19 times as the registration source and 19 times as the target. We segmented the brain image into three tissues: white matter (WM), gray matter (GM), cerebrospinal fluid (CSF) using the software FSL [35–37]. Then, we adopted two measurements to estimate the effectiveness of the registration: Target overlap [38] and Dice coefficient [39] and their definitions are: T S T ρ ¼ O \ O O
=
ð15Þ
and T T S S ρ ¼ 2 O \ O = O þ O
ð16Þ
where O T and O S denote the ROI in the template and subject image, respectively. Target overlap and Dice coefficient measure varies in the range [0, 1]. Table 3 presented a registration performance comparison of four registration methods. The average target overlay and dice coefficient and their standard deviation of four methods were listed in Table 3.
Our method
The performance of Demons algorithm was weaker than that of others. It can be seen that the performances of DARTEL, SyN and our methods were very close to each other. The target overlap from our method was roughly 0.835 in case of WM, 0.722 for GM and 0.807 for CSF. The dice coefficient from our method was roughly 0.820 in case of WM, 0.683 for GM and 0.811 for CSF. These values were slightly higher than or close to the results using other algorithms. As a whole, our method was comparable to other methods used extensively in the medical imaging community.
4. Conclusion We have developed a new approach based on three dimensional symmetric optic flow models for deformable registration of brain images. The proposed algorithm takes FA maps to generate the point correspondence between the subject and template spaces. We designed it so as to integrate the DTI-specific characteristics into the registration framework. This framework may not be suitable for neuroimaging DTI studies that actually study FA variances in diseases or across populations, but will be fine for those studying other DAIs (e.g., D?, DO, MD, etc.), because this framework is supposed to already remove all possible FA and morphometric variances during the normalization procedure. However, the framework is ready to use other replacements instead of the FA images for generating the point correspondence. One alterative is the baseline data acquired along with the DWI data, which are actually heavily T2-weighted images. Our algorithm is a fully automated process of registration for brain image data and the software implements the approach presented in http://www.columbia.edu/~dx2103/brainimagescope. html. Compared to other advanced algorithms currently available such as SPM5 DARTEL [22], HAMMER [7] and Timer [5], our algorithm is relatively more convenient and user-friendly with comparably high accuracy and efficiency. We plan to further evaluate and validate the method in larger datasets, especially pediatric databases in the future.
Table 3 The registration results comparison of demons, SPM5 dartel, SYN and our method. Method
Demons
DARTEL
SyN
Our method
0.831 ± 0.038 0.705 ± 0.052 0.785 ± 0.127
0.837 ± 0.024 0.711 ± 0.067 0.788 ± 0.138
0.835 ± 0.025 0.722 ± 0.059 0.807 ± 0.135
Target overlap WM GM CSF
0752 ± 0.120 0666 ± 0.121 0.741 ± 0.166
Demons
DARTEL
SyN
Our Method
0.817 ± 0.058 0.661 ± 0.067 0.791 ± 0.106
0.821 ± 0.047 0.672 ± 0.087 0.790 ± 0.091
0.820 ± 0.049 0.683 ± 0.069 0.811 ± 0.097
Dice coefficient 0767 ± 0.164 0.620 ± 0.110 0.753 ± 0.131
Please cite this article as: Wen Y, et al, A highly accurate symmetric optical flow based high-dimensional nonlinear spatial normalization of brain images, Magn Reson Imaging (2015), http://dx.doi.org/10.1016/j.mri.2015.01.013
Y. Wen et al. / Magnetic Resonance Imaging xxx (2015) xxx–xxx
References [1] Du J, Goh A, Qiu A. Large deformation diffeomorphic metric mapping of orientation distribution functions. Int. Conf. Information Processing in Medical Imaging, 22; 2011. p. 448–62. [2] Geng X, Ross TJ, Gu H, Shin W, Zhan W, Chao YP, et al. Diffeomorphic image registration of diffusion MRI using spherical harmonics. IEEE Trans Med Imaging 2011;30(3):747–58. [3] Duarte-Carvajalino JM, Sapiro G, Harel N, Lenglet C. A Framework for Linear and Non-Linear Registration of Diffusion-Weighted MRIs Using Angular Interpolation. Front Neurosci 2013:7–41. [4] Bookstein FL. Principal warps: Thin-plate splines and the decomposition of deformations. IEEE Trans Pattern Anal Machine Intell 1989;11:567–85. [5] Rohr K. “Image registration based on thin plate splines and local estimates of anisotropic landmark localization uncertainties”. In Proc. Conf. Medical Image Computing and Computer-Assisted Intervention (MICCAI’98); 1998. p. 1174–83. [6] Vaillant M, Davatzikos C. Hierarchical matching of cortical features for deformable brain image registration. Proc. Int. Conf. Information Processing in Medical Imaging; 1999. p. 182–95. [7] Chui H, Win L, Schultz R, Duncan J, Rangarajan A. A unified feature registration method for brain mapping. Proc. Int. Conf. Information Processing in Medical Imaging; 2001. p. 300–14. [8] Yap PT, Wu GR. TIMER: Tensor Image Morphing for Elastic Registration. NeuroImage 2009;47:549–63. [9] Yap PT, Wu GR, Zhu HT, Lin W, Shen DG. F-TIMER: Fast Tensor Image Morphing for Elastic Registration. IEEE Trans Med Imaging 2010;29(5):1192–203. [10] Shen D, Davatzikos C. HAMMER: heirarchical attribute matching mechanism for elastic registration. IEEE Trans Med Imaging 2002;21(11):1421–39. [11] Greve DN, Fishl B. Accurate and robust brain image alignment using boundarybased registration. NeuroImage 2009;48:63–72. [12] Yang JZ, Shen DG, Misra C, Xu XY, Resnick S, Davatzikos C, et al. “Spatial normalization of diffusion tensor images based on anisotropic segmentation”. Medical Imaging 2008: Image Processing 2008;6914:69140 L1–10. [13] Xue Z, Hai L, Lei G, Stephen TC. A local fast marching-based diffusion tensor image registration algorithm by simultaneously considering spatial deformation and tensor orientation. NeuroImage 2010;52:119–30. [14] Memin E, Perez P. Hierarchical estimation and segmentation of dense motion fields. Int J Comput Vis 2002;46(2):129–55. [15] Black MJ, Anandan P. Robust dynamic motion estimation over time. Proc. 1991 IEEE Computer Society Conference on Computer Vision and Pattern Recognition. IEEE Computer Society Press; 1991. p. 292–302. [16] Black MJ, Anandan P. The robust estimation of multiple motions: parametric and piecewise smooth flow fields. Comput Vis Image Underst 1996;63(1): 75–104. [17] Christensen GE, Rabbitt RD, Miller MI. Deformable templates using large deformation kinematics. IEEE Trans Image Process 1996;5(10):1435–47. [18] Horn B, Schunck B. Determining optical flow. Artif Intell 1981;17:185–203.
9
[19] Lucas B, Kanade T. An iterative image registration technique with an application to stereo vision. Proc. Seventh International Joint Conference on Artificial Intelligence; 1981. p. 674–9 [Vancouver, Canada]. [20] Xu L, Jia J, Matsushita Y. Motion detail preserving optical flow estimation. IEEE Trans Pattern Anal Mach Intell 2012;34:1744–57. [21] Stoll M, Volz S, Bruhn A. Adaptive integration of feature matches into variational optical flow methods. Computer Vision–ACCV 2012. Springer; 2013. p. 1–14. [22] Avants BB, Epstein CL, Grossman M, Gee JC. Symmetric diffeomorphic image registration with cross-correlation: evaluating automated labeling of elderly and neurodegenerative brain. Med Image Anal 2008;12:26–41. [23] Ashburner J. A fast diffeomorphic image registration algorithm. NeuroImage 2007;38:95–113. [24] Thirion JP. Non-rigid matching using demons. IEEE Computer Vision and Pattern Recognition; 1996. p. 245–51. [25] Fox Charles. An introduction to the calculus of variations. Courier Dover Publications; 1987. [26] David M. Iterative methods for solving partial difference equations of elliptical type. (PhD thesis) Harvard University; 1950. [27] Hill S Francis. Computer Graphics. Prentice Hall; 2001. [28] Black MJ, Anandan P. “Robust dynamic motion estimation over time”. In Proc 1991 IEEE Computer Society Conference on Computer Vision and Patter Recognition; 1991. p. 292–302. [29] Yin W, Osher S, Goldfarb D, Darbon J. Bregman Iterative Algorithms for L1Minimization with Applications to Compressed Sensing, SIAM J. Imaging Sciences 2008;1(1):143–68. [30] Sun D, Roth S, Black M. A Quantitative Analysis of Current Practices in Optical Flow Estimation and the Principles Behind Them. Int J Comput Vis 2013:1–23. [31] Álvarez L, Castano CA, García M, Krissian K, Mazorra L, Salgado A, et al. Symmetric optical flow. Computer Aided Systems Theory–EUROCAST 2007. Springer; 2007. p. 676–83. [32] Alvarez L, Deriche R, Papadopoulo T, Sánchez J. Symmetrical dense optical flow estimation with occlusions detection. Int J Comput Vis 2007;75:371–85. [33] Yoo T. Insight Into Images: Principles and Practice for Segmentation, Registration and Image Analysis. Natick, MA: A.K. Peters Ltd; 2003. [34] Klein A, Andersson J, Ardekani BA, Ashburner J, et al. Evaluation of 14 nonlinear deformation algorithms applied to human brain MRI registration. NeuroImage 2009;46:786–802. [35] Smith SM, Jenkinson M, Woolrich MW, Beckmann CF, Behrens TEJ, JohansenBerg H, et al. Advances in functional and structural MR image analysis and implementation as FSL. NeuroImage 2004;23(S1):208–19. [36] Jenkinson M, Beckmann CF, Behrens TE, Woolrich MW, Smith SM. FSL. NeuroImage 2012;62:782–90. [37] Jenkinson M, Smith SM. Bayesian analysis of neuroimaging data in FSL. NeuroImage 2009;45:S173–86. [38] Wu G, Kim M, Wang Q, Shen D. S-HAMMER: Hierarchical attribute‐guided, symmetric diffeomorphic registration for MR brain images. Hum Brain Mapp 2014;35(3):1044–60. [39] Chang HH, Zhuang AH, Valentino DJ, Chu WC. Performance measure characterization for evaluating neuroimage segmentation algorithms. NeuroImage 2009; 47:122–35.
Please cite this article as: Wen Y, et al, A highly accurate symmetric optical flow based high-dimensional nonlinear spatial normalization of brain images, Magn Reson Imaging (2015), http://dx.doi.org/10.1016/j.mri.2015.01.013