Diffusion Tensor Imaging with Deterministic Error Bounds

J Math Imaging Vis (2016) 56:137–157 DOI 10.1007/s10851-016-0639-7 Diffusion Tensor Imaging with Deterministic Error Bounds Artur Gorokh1 · Yury Korolev2 · Tuomo Valkonen3 Received: 7 September 2015 / Accepted: 8 February 2016 / Published online: 29 February 2016 © The Author(s) 2016. This article is published with open access at Springerlink.com Abstract Errors in the data and the forward operator of an inverse problem can be handily modelled using partial order in Banach lattices. We present some existing results of the theory of regularisation in this novel framework, where errors are represented as bounds by means of the appropriate partial order. We apply the theory to diffusion tensor imaging, where correct noise modelling is challenging: it involves the Rician distribution and the non-linear Stejskal–Tanner equation. Linearisation of the latter in the statistical framework would complicate the noise model even further. We avoid this using the error bounds approach, which preserves simple error structure under monotone transformations. Keywords Diffusion tensor imaging · Noise modelling · Total generalised variation · Error bounds · Deterministic Mathematics Subject Classification 92C55 · 94A08 1 Introduction Often in inverse problems, we have only very rough knowledge of noise models, or the exact model is too difficult to realise in a numerical reconstruction method. The data may also contain process artefacts from black-box devices [44]. Partial order in Banach lattices has therefore recently been investigated in [32–34] as a less-assuming error modelling approach for inverse problems. This framework allows the representation of errors in the data as well as in the forward operator of an inverse problem by means of order intervals (i.e. lower and upper bounds by means of appropriate partial orders). An important advantage of this approach vs. statistical noise modelling is that deterministic error bounds preserve their simple structure under monotone transformations. We apply partial order in Banach lattices to diffusion tensor imaging (DTI). We will in due course explain the diffusion tensor imaging progress, as well as the theory of inverse problems in Banach lattices, but start by introducing our model min R(u) s.t. u  0, u B Tuomo Valkonen Artur Gorokh Yury Korolev 1 Faculty of Physics, Lomonosov Moscow State University, Moscow, Russia 2 School of Engineering and Materials Science, Queen Mary University of London, London, UK 3 Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Cambridge, UK glj  A j u  g uj , ( j = 1, . . . , N ). That is, we want to find a field of symmetric 2-tensors u : Ω → Sym2 (R3 ) on the domain Ω ⊂ R3 , minimising the value of the regulariser R on the feasible set. The tensor field u is our unknown image. It is subject to a positivity constraint, as well as partial order constraints imposed through the operators [A j u](x) := −b j , u(x)b j , and the upper and lower bounds glj := log(ŝ lj /ŝ0u ) and g uj := log(ŝ uj /ŝ0l ). These arise from the linearisation (via the monotone logarithmic transformation) of the Stejskal-Tanner equation s j (x) = s0 (x) exp(−b j , u(x)b j ), ( j = 1, . . . , N ), (1) central to the diffusion tensor imaging process. 123 138 J Math Imaging Vis (2016) 56:137–157 To shed more light on u and the equation (1), let us briefly outline the diffusion tensor imaging process. As a first step towards DTI, diffusion-weighted magnetic resonance imaging (DWI) is performed. This process measures the anisotropic diffusion of water molecules. To capture the diffusion information, the magnetic resonance images have to be measured with diffusion-sensitising gradients in multiple directions. These are the different bi ’s in (1). Eventually, multiple DWI images {s j } are related through the Stejskal–Tanner equation (1) to the symmetric positive-definite diffusion tensor field u : Ω → Sym2 (R3 ) [4,29]. At each point x ∈ Ω, the tensor u(x) is the covariance matrix of a normal distribution for the probability of water diffusing in different spatial directions. The fact that multiple bi ’s are needed to recover u leads to very long acquisition times, even with ultra fast sequences like echo planar imaging (EPI). Therefore, DTI is inherently a low-resolution and low-SNR method. In theory, the amplitude DWI images exhibit Rician noise [24]. However, as the histogram of an in vivo measurement in Fig. 1 illustrates, this may not be the case for practical datasets from black-box devices. Moreover, the DWI process is prone to eddy-current distortions [53], and due to the slowness of it, it is very sensitive to patient motion [1,27]. We therefore have to use techniques that remove these artefacts in solving for u(x). We also need to ensure the positivity u, as non-positive-definite diffusion tensor are non-physical. One proposed approach for the satisfaction of this constraint is that of log-Euclidean metrics [3]. This approach has several theoretically desirable aspects, but some practical shortcomings [57]. Special Perona–Malik-type constructions on Riemannian manifolds can also be used to maintain the structure of the tensor field [14,54]. Such anisotropic diffusion is, however, severely illposed [60]. Recently manifold-valued discrete-domain total variation models have also been applied to diffusion tensor imaging [6]. Our approach is also in the total variation family, first considered for diffusion tensor imaging in [48]. Namely, we follow up on the work in [55,57–59] on the application of total generalised variation regularisation [9] to DTI. We should note that in all of these works, the fidelity function was the ROF-type [45] L 2 fidelity. This would only be correct, according to the assumption that noise of MRI measurements is Gaussian, if we had access to the original complex k-space MRI data. The noise of the inverse Fourier-transformed magnitude data s j , that we have in practice access to, is however Rician under the Gaussian assumption on the original complex data [24]. This is not modelled by the L 2 fidelity. Numerical implementation of Rician noise modelling has been studied in [22,40]. As already discussed, in this work, we take the other direction. Instead of modelling the errors in a statistically accurate fashion, not assuming to know an exact noise model, we represent them by means of pointwise bounds. The details of the model are presented in Sect. 3. We study the practical performance in Sect. 5 using the numerical method presented in Sect. 4. First we, however, start with the general error modelling theory in Sect. 2. Readers who are not familiar with notation for Banach lattices or symmetric tensor fields are advised to start with the Appendix 2 where we introduce our mathematical notation and techniques. Fig. 1 The noise in the absolute values of complex MRI data should be Rician. Here we have taken a 50-bin histogram of the noise in (...truncated)