Diffusion tensor imaging with deterministic error bounds

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 nonlinear 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.


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 [41]. Partial order in Banach lattices has therefore recently been investigated in [30,32,31] 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 u R(u) subject to u 0, g l j A j u g u j , L n -a.e. on Ω, (j = 1, . . . , N ).
That is, we want to find a field of symmetric 2-tensors u : Ω → Sym 2 (R 3 ) on the domain Ω ⊂ R 3 , minimising the value of the regulariser R on the feasible set. The tensor field u is our unknown image. It is subject 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 g l j := log(ŝ l j /ŝ u 0 ) and g u j := log(ŝ u j /ŝ l 0 ). These model, in terms of error intervals after logarithmic transformation, the Stejskal-Tanner equation s j (x) = s 0 (x) exp(− b j , u(x)b j ), (j = 1, . . . , N ), (1.1) central to the diffusion tensor imaging process.
To shed more light on u and the equation (1.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 b i 's in (1.1). Eventually, multiple DWI images {s j } are related through the Stejskal-Tanner equation (1.1) to the symmetric positivedefinite diffusion-tensor field u : Ω → Sym 2 (R 3 ) [4,27]. 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 b i '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 [23]. However, as the histogram of an in vivo measurement in Figure 1 illustrates, this may not be the case for practical data sets from black-box devices. Moreover, the DWI process is prone to eddy-current distortions [50], and due to the slowness of it, it is very sensitive to patient motion [26,1]. 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 [54]. Special Perona-Malik type constructions on Riemannian manifolds can also be used to This divides the pixels into bins of 50 different noise levels. However, we only find approximately 10 noise levels to have non-zero pixel count. As the Rician distribution is continuous, we see that the noise cannot be Rician, some bins of the 50-bin histogram being empty. The measurement setup of the data used here is described in Section 5.3. maintain the structure of the tensor field [14,51]. Such anisotropic diffusion is however severely ill-posed [?]. Recently manifold-valued discretedomain 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 [45]. Namely, we follow up on the work in [54,56,55,52] 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 [42] 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 [23]. This is not modelled by the L 2 fidelity.
Numerical implementation of Rician noise modelling has been studied in [37,21]. 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 Section 3. We study the practical performance in Section 5 using the numerical method presented in Section 4. First we however start with the general error modelling theory in Section 2. Readers who are not familiar with notation for Banach lattices or symmetric tensor fields are advised to start with the Appendix, where we introduce our mathematical notation and techniques.

Mathematical basis
We now briefly outline the theoretical framework [30] that is the basis for our approach. Consider a linear operator equation where U and F are Banach lattices, A : U → F is a regular injective operator. The inaccuracies in the right-hand side f and the operator A are represented as bounds by means of appropriate partial orders, i.e.
where the symbol F stands for the partial order in F and L ∼ (U,F ) for the partial order for regular operators induced by partial orders in U and F . Further, we will drop the subscripts at inequality signs where it will not cause confusion. The exact right-hand side f and operator A are not available. Given the approximate data (f l , f u , A l , A u ), we need to find an approximate solution u that converges to the exact solutionū as the inaccuracies in the data diminish. This statement needs to be formalised. We consider monotone convergent sequences of lower and upper bounds We are looking for an approximate solution u n such that u n −ū U → 0 as n → ∞. Let us ask the following question. What are the elements u ∈ U that could have produced data within the tolerances (2.3)? Obviously, the exact solution is one of such elements. Let us call the set containing all such elements the feasible set U n ⊂ U .
Suppose that we know a priori that the exact solution is positive (by means of the appropriate partial order in U ). Then it is easy to verify that the following inequalities hold for all n ∈ N This observation motivates our choice of the feasible set: It is clear that the exact solutionū belongs to the sets U n for all n ∈ N. Our goal is to define a rule that will choose for any n an element u n of the set U n such that the sequence u n ∈ U n will strongly converge to the exact solutionū. We do so by minimising an appropriate regularisation functional R(u) on U n : u n = arg min u∈Un R(u). (2.4) This method, in fact, is a lattice analogue of the well-known residual method [22,49]. The convergence result is as follows [30].
C} (C = const) are sequentially compact in U (in the strong topology induced by the norm).
Then the sequence defined in (2.4) strongly converges to the exact solutionū and R(u n ) → R(ū).
Examples of regularisation functionals that satisfy the conditions of Theorem 2.1 are as follows. Total Variation in L 1 (Ω), where Ω is a subset of R n , assures strong convergence in L 1 , given that the L 1 -norm of the solution is bounded. The Sobolev norm u W 1,q (Ω) yields strong convergence in the spaces L p (Ω), where p 1, q > np p+n . The latter fact follows from the compact embedding of the corresponding Sobolev W 1,q (Ω) space into L p (Ω) [17].
However, the assumption that the sets {u : R(u) C} are strong compacts in U is quite strong. It can be replaced by the assumption of weak compactness, provided that the regularisation functional possesses the socalled Radon-Riesz property. Definition 2.1. A functional F : U → R has the Radon-Riesz property (sometimes referred to as the H-property), if for any sequence u n ∈ U weak convergence u n u 0 and simultaneous convergence of the values F (u n ) → F (u 0 ) imply strong convergence u n → u 0 . Then the sequence defined in (2.4) strongly converges to the exact solutionū and R(u n ) → R(ū).
It is easy to verify that the norm in any Hilbert space possesses the Radon-Riesz property. Moreover, this holds for the norm in any reflexive Banach space [17].
As we explain in the Appendix, the spaces L p (Ω; Sym 2 (R m )) are not Banach lattices, therefore, Theorems 2.1 and 2.2 cannot be applied directly. Further theoretical work will be undertaken to extend the framework to the non-lattice case. For the moment, however, we will prove that if there are no errors in the operator A in (2.1), the requirement that the solution space U is a lattice can be dropped. Proof. Let us define an approximate right-hand side and its approximation error in the following way One can easily verify, that the inequality f − f δn δ n holds. Indeed, we have The first two inequalities are consequences of the conditions (2.3), the third one holds by the definition of supremum and the equality |f | = f ∨ (−f ) that holds for all f ∈ F , and the last inequality is due to the monotonicity of the norm in a Banach lattice. Similarly, one can show that for any u ∈ U n , we have Therefore, the inclusion U n ⊂ {u : Au − f δn δ n } holds. Now we will proceed with the proof of convergence u n −ū → 0. Will prove it for the case when the regulariser R(u) satisfies conditions of Theorem 2.1. Suppose that the sequence u n does not converge to the exact solutionū. Then it contains a subsequence u n k such that u n k −ū ε for any k ∈ N and some fixed ε > 0.
Since the inclusionū ∈ U n holds for all n ∈ N, we have R(u n ) R(ū) for all n ∈ N. Since the level set {u : R(u) R(ū)} is a compact set by assumptions of the theorem, the sequence u n k contains a strongly convergent subsequence. With no loss of generality, let us assume that u n k → u 0 . We will now show that u 0 =ū. Indeed, we have On the other hand, we have Au n k − Aū → Au 0 − Aū due to continuity of A and · . Therefore, Au 0 = Aū and u 0 =ū, since A is an injective operator. By contradiction, we get u n −ū → 0. Finally, since the regulariser R(u) is lower semi-continuous, we get that lim inf R(u n ) = R(ū). However, for any n we have R(u n ) R(ū), therefore, we get the convergence R(u n ) → R(ū) as n → ∞.

Philosophical discussion and statistical interpretation
In practice, our data is discrete. So let us momentarily switch to measurementsf = (f 1 , . . . ,f n ) ∈ R n of a true data f ∈ R n . If we actually knew the pointwise noise model of the data, then one way to obtain potentially useful upper and lower bounds for the deterministic model is by means of statistical interval estimates: confidence intervals. Roughly, the idea is to find for each true signal f individual random upper and lower boundsf u andf l such that Iff u andf l are computed based on multiple experiments (i.e., multiple noisy samplesf 1 , . . . ,f m , of the true data f ), the interval [f u,i ,f l,i ] will converge in probability to the true dataf i , as the number of experiments m increases. Thus we obtain a probabilistic version of the convergences in (2.3). Let us try to see, how such intervals might work in practice. For the purpose of the present discussion, assume that the noise is additive and normal-distributed with variance σ and zero mean-an assumption that does not hold in practice, as we have already seen in Figure 1, but will suffice for the next thought experiments. That is,f j = f + ν j for the noise ν j . Let the sample mean of {f j } m j=1 bef m = (f 1 m , . . .f n m ), and pointwise sample variance σ = (σ 1 , . . . , σ n ). The product of the pointwise confidence intervals I i with confidence 1 − θ is [15,46] for Φ the cumulative normal distribution function. For θ = 0.05, i.e., the 95% confidence interval, Φ −1 (0.05/2) = 1.96. Now, the probability that I i covers the true f i is 1 − θ, e.g. 95%. If we have only a single sample m = 1, the intervals stay large, but the joint probability, (1 − θ) n goes to zero as n increases. As an example, for a rather typical single 128 × 128 slice of a DTI measurement, the probability that exactly φ = 5% (to the closest discrete value possible) of the 1 − θ = 95% confidence intervals do not cover the true parameter would be about 1.4%, or The probability of at least φ = 5% of the pointwise 95% confidence intervals not covering the true parameter is in this setting approximately 49%. This can be verified by summing the above estimates over m = φn , . . . , n.
In summary, unless θ simultaneously goes to 1, the product intervals are very unlikely to cover the true parameter. Based on a single experiment, the deterministic approach as interpreted statistically through confidence intervals, is therefore very likely to fail to discover the true solution as the data size n increases unless the pointwise confidence is very low. But, if we let the pointwise confidences be arbitrarily high, such that the intervals are very large, the discovered solution in our applications of interest would be just a constant! Asymptotically, the situation is more encouraging. Indeed, if we could perform more experiments to compute the confidence intervals, then for any fixed n and θ, it is easy to see that the solution of the "deterministic" error model is an asymptotically consistent and hence asymptotically unbiased estimator of the true f . That is, the estimates converge in probability to f as the experiment count m increases. Indeed, the error-bounds based estimatorf m , based on m experiments, by definition satisfiesf m ∈ n i=1 I i . Therefore, we have Thus f P →f in probability. Since by the law of large numbers alsof m P → f , this proves the claim, and to some extent justifies our approach from the statistical viewpoint.
It should be noted that this is roughly the most that has previously been known of the maximum a posteriori estimate (MAP), corresponding to the Tikhonov models In particular, the MAP is not the Bayes estimator for the typical squared cost functional. This means that it does not minimisef The minimiser in this case is the conditional mean (CM) estimate, which is why it has been preferred by Bayesian statisticians despite its increased computational cost. The MAP estimate is merely an asymptotic Bayes estimator for the uniform cost function. In a very recent work [12], it has however been proved that the MAP estimate is the Bayes estimator for certain Bregman distances. One possible critique of the result is that these distances are not universal and do depend on the regulariser R, unlike the squared distance for CM. The CM estimate however has other problems in the setting of total variation and its discretisation [34,33].

Application to diffusion tensor imaging
We now build our model for applying the deterministic error modelling theory to diffusion tensor imaging. We start by building our forward model based on the Stejskal-Tanner equation, and then briefly introduce the regularisers we use.

The forward model
For u : Ω → Sym 2 (R 3 ), Ω ⊂ R 3 , a mapping from Ω to symmetric second order tensors, let us introduce non-linear operators T j , defined by Their role is to model the so-called Stejskal-Tanner equation [4] s Each tensor u(x) models the covariance of a Gaussian probability distribution at x for the diffusion of water molecules. The data s j ∈ L 2 (Ω), (j = 1, . . . , N ), are the diffusion-weighted MRI images. Each of them is obtained by performing the MRI scan with a different non-zero diffusion sensitising gradient b j , while s 0 is obtained with a zero gradient. After correcting the original k-space data for coil sensitivities, each s j is assumed real. As a consequence, any measurementŝ j of s j has-in theory-Rician noise distribution [23].
Our goal is to reconstruct u with simultaneous denoising. Following [52,29], we consider using a suitable regulariser R the Tikhonov model The constraint u 0 is to be understood in the sense that u(x) is positive semidefinite for L n -a.e. x ∈ Ω (see Appendix for more details). Due to the Rician noise ofŝ j , the Gaussian noise model implied by the L 2 -norm in (3.2) is not entirely correct. However, in some cases the L 2 model may be accurate enough, as for suitable parameters the Rician distribution is not too far from a Gaussian distribution. If one were to model the problem correctly, one should either modify the fidelity term to model Rician noise, or include the (unit magnitude complex number) coil sensitivities in the model. The Rician noise model is highly nonlinear due to the Bessel functional logarithms involved. Its approximations have been studied in [5,21,37] for single MR images and DTI. Coil sensitivities could be included either by knowing them in advance, or by simultaneous estimation as in [28]. Either way, significant complexity is introduced into the model, and for the present work, we are content with the simple L 2 model.
We may also consider, as is often the case, and as was done with TGV in [54], the linearised model where, for each x ∈ Ω, f (x) is solved by regression for u(x) from the system of equations (3.1) with s j (x) =ŝ j (x). Further, as in [55], we may also consider ). In both of these linearised models, the assumption of Gaussian noise is in principle even more remote from the truth than in the nonlinear model (3.2). We will employ (3.3) and (3.2) as benchmark models. We want to further simplify the model, and forgo with accurate noise modelling. After all, we often do not know the real noise model for the data available in practice. It can be corrupted by process artefacts from blackbox algorithms in the MRI devices. This problem of black box devices has been discussed extensively in [41], in the context of Computed Tomography. Moreover, as we have discussed above, even without such artefacts, the correct model may be difficult to realise numerically. So we might be best off choosing the least assuming model of all -that of error bounds. This is what we propose in the reconstruction model (3.5) Here g l j := log(ŝ l j /ŝ u 0 ) and g u j := log(ŝ u j /ŝ l 0 ), g l j , g u j ∈ L 2 (Ω), are our upper and lower bounds on g j that we derive from the data.

Choice of the regulariser R
A prototypical regulariser in image processing is the total variation, first studied in this context in [42]. It can be defined for a symmetric tensor field u ∈ L 1 (Ω; Sym k (R m )) as Observe that for scalar or vector fields, i.e., the cases k = 0, 1, we have Therefore, for scalars in particular, this gives the usual isotropic total variation Total generalised variation was introduced in [9] as a higher-order extension of TV. Following [54], the second-order variant may be defined using the differentiation cascade formulation for symmetric tensor fields u ∈ L 1 (Ω; Sym k (R m )) as the marginal It turns out that the standard BV norm and the "BGV norm" [9] u : are topologically equivalent norms [10,11] on BV(Ω; Sym k (R m )), yielding the same convergence results for TGV regularisation as for TV regularisation. The geometrical regularisation behaviour is however different, and TGV tends to avoid the staircasing observed in TV regularisation.
Regarding topologies, we say that a sequence {u i } in BV(Ω; Sym k (R m )) converges weakly* to u, if u i → u strongly in L 1 , and Eu i * Eu weakly* as Radon measures [2,48,54]. The latter is characterised as
holds, we may then bound By the weak* lower semicontinuity of the BV-norm, and the weak* compactness of the unit ball in BV(Ω; Sym k (R m ))-we refer to [2] for these and other basic properties of BV-spaces-we may thus find a representativeũ in the BV 0,R (Ω; Sym k (R m )) equivalence class of u, satisfying Again using the weak* compactness of the unit ball in BV(Ω; Sym k (R m )), and the weak* lower semicontinuity of R, it follows that the sets lev a R := {u ∈ BV 0,R (Ω; Sym k (R m )) | R(u) ≤ a}, (a > 0), are weak* compact in BV 0,R (Ω; Sym k (R m )), in the topology inherited form BV(Ω; Sym k (R m )). Consequently, they are strongly compact subsets of L 1 (Ω; Sym k (R m )). This feature is crucial for the application of the regularisation theory in Banach lattices above. On a connected domain Ω, in particular That is, the space consists of zero-mean functions. Then u → Du M(Ω;R m ) is a norm on BV 0,TV (Ω) [39], and this space is weak* compact. In particular, the sets lev a TV are compact in L 1 (Ω).
More generally, we know from [8] that on a connected domain Ω, ker TV consists of Sym k (R m )-valued polynomials of maximal degree k. By extension, ker TGV 2 consists of Sym k (R m )-valued polynomials of maximal degree k + 1. In both cases, (3.7), weak* lower semicontinuity of R, and the equivalence of · to · BV(Ω;Sym k (R m )) hold by the results in [8,11,48]. Therefore, we have proved the following.

Lemma 3.1.
Let Ω ⊂ R m and k ≥ 0. Then the sets lev a TV and lev a TGV 2 are weak* compact in BV(Ω; Sym k (R m )) and strongly compact in L 1 (Ω; Sym k (R m )). Now, in the above cases, ker R is finite-dimensional, and we may write BV(Ω; Sym k (R m )) BV 0,R (Ω; Sym k (R m )) ⊕ ker R.

Denoting by
the closed ball of radius r in a normed space X, we obtain by the finitedimensionality of ker R the following result.
for some constant a > 0, then for R = TV and R = TGV 2 , the sequence u n = arg min u∈Un R(u) converges strongly in L 1 (Ω; Sym k (R m )) to the exact solutionū and R(u n ) → R(ū).
Proof. With the decomposition u n = u 0,n + u ⊥ n , where u ⊥ n ∈ ker R, we have u 0,n ∈ lev a R for suitably large a > 0 through The assumption (3.8) bounds u ⊥ n ≤ a. Thus u n ∈ V for V as in Proposition 3.1. The proposition thus implies the necessary compactness in U = L 1 (Ω; Sym k (R m )) for the application of Theorem 2.3.
Remark 3.1. The condition (3.8) simply says for R = TV that the data has to bound the solution in mean. This is very reasonable to expect for practical data; anything else would be very non-degenerate. For R = TGV 2 we also need that the data bounds the entire affine part of the solution. Again, this is very likely for real data. Indeed, in DTI practice, with at least 6 independent diffusion sensiting gradients, A is an invertible or even overdetermined linear operator. In that typical case, the bounds f l n and f u n will be translated into U n being a bounded set.

Solving the optimisation problem 4.1 The Chambolle-Pock method
The Chambolle-Pock algorithm is an inertial primal-dual backward-backward splitting method, classified in [19] as the modified primal-dual hybrid gradient method (PDHGM). It solves the minimax problem where G : X → R and F * : Y → R are convex, proper, lower semicontinuous functionals on (finite-dimensional) Hilbert spaces X and Y . The operator K : X → Y is linear, although an extension of the method to nonlinear K has recently been derived [52]. The PDHGM can also be seen as a preconditioned ADMM (alternating directions method of multipliers); we refer to [19,44,53] for reviews of optimisation methods popular in image processing. For step sizes τ, σ > 0, and an over-relaxation parameter ω > 0, each iteration of the algorithm consists of the updates We should remark that the order of the primal (u) and dual (y) updates here is reversed from the original presentation in [13]. The reason is that when reordered, the updates can, as discovered in [25], be easily written in a proximal point form.
The first and last update are the backward (proximal) steps for the primal (x) and dual (y) variables, respectively, keeping the other fixed. However, the dual step includes some "inertia" or over-relaxation, as specified by the parameter ω. Usually ω = 1, which is required for convergence proofs of the method. If G or F * is uniformly convex, by smartly choosing for each iteration the step length parameters τ, σ, and the inertia ω, the method can be shown to have convergence rate O(1/N 2 ). This is similar to Nesterov's optimal gradient method [40]. In the general case the rate is O(1/N ). In practice the method produces visually pleasing solutions in rather few iterations, when applied to image processing problems.
In implementation of the method, it is crucial that the resolvents (I + τ ∂G) −1 and (I + σ∂F * ) −1 can be computed quickly. We recall that they may be written as Usually in applications, these computations turn out to be simple projections or linear operations -or the soft-thresholding operation for the L 1norm.
As a further implementation note, since the algorithm (4.2) is formulated in Hilbert spaces (see however [?]), and we work in the Banach space BV(Ω; Sym 2 (R 3 )), we have to discretise our problems before application of the algorithm. We do this by simple forward-differences discretisation of the operator E with cell width h = 1 on a regular rectangular grid corresponding to the image voxels.

Implementation of deterministic constraints
We now reformulate the problem (3.5) of DTI imaging with deterministic error bounds in the form (4.1). Suppose we have some upper and lower bounds s l j s j s u j on the DWI signals s j , (j = 0, . . . , N ). Then the bounds for g j = log(s j /s 0 ) are because g j is monotone in regards to s j . We are thus trying to solve For the ease of notation, we write g = g 1 , . . . , g N , and so that the Stejskal-Tanner equation is satisfied with The problem (4.4) may be rewritten as for F 0 (y) = δ(g l y g r ) with δ(g l y g u ) denoting the indicator function of the convex set {y : g l y g r }. Solving the conjugate F * 0 (y) = g l , y , y < 0 g u , y , y ≥ 0. , and also writing R(u) = R 0 (K 0 u) for some R 0 and K 0 , the problem can further be written in the saddle point form (4.1) with G(u) = 0, F * (y, ψ) = F * 0 (y) + R * 0 (ψ). To apply algorithm (4.2), we need to compute the resolvents of G * 0 and R * 0 . For details regarding R * 0 for R = TGV 2 (β,α) and R = α TV in the discretised setting, we refer to [54,55]; here it suffices to note that for R = α TV, we have K 0 = αE and R * 0 (φ) is the indicator function of the dual ball {φ | sup x φ(x) F ≤ 1}. Thus the resolvent (I + τ ∂R * 0 ) −1 becomes a projection to the dual ball. The case of R = TGV 2 (β,α) is similar. For F * 0 we have which resolves pointwise at each ξ ∈ Ω into the expression Finally, we note that the saddle point system (4.1) has to have a solution for the Chambolle-Pock algorithm to converge. In our setting, in particular, we need to find error bounds g l and g u , for which there exists a solution u to g l Au g u . If one uses at most six independent diffusion directions (N = 6), as we will, then, for any g, there in fact exists a solution to g = Au. The condition (4.6) becomes g l g u , immediately guaranteed through the monotonicity of (4.3), and the trivial conditions s l j s u j . We are thus ready to apply the algorithm (4.2) to diffusion tensor imaging with deterministic error bounds. For the realisation of (4.2) for models (3.3), (3.4), and (3.2), we refer to [54,55,56,52].

Experimental results
We now study the efficiency of the proposed reconstruction model in comparison to the different L 2 -squared models, i.e., ones with Gaussian error assumption. This is based on a synthetic data, for which a ground-truth is available, as well as a real in vivo DTI data set. First we, however, have to describe in detail the procedure for obtaining upper and lower bounds for real data, when we do not know the true noise model, and are unable to perform multiple experiments as required by the theory of Section 2.2.

Estimating lower and upper bounds from real data
As we have already discussed, in practice the noise in the measurement signalsŝ j is not Gaussian or Rician; in fact we do not know the true noise distribution and other corruptions. Therefore, we have to estimate the noise distribution from the image background. To do this, we require a known correspondence between the measurement, the noise, and the true value. As we have no better assumptions available, the standard one that we use is that of additive noise. Continuing in the statistical setting of Section 2.2, we now describe the procedure, working on discrete images expressed as vectorsf = s j ∈ R n for some fixed j ∈ {0, 1, . . . , N }. We use superscripts to denote the voxel indices, that isf = (f 1 , . . . , f n ).
In the i-th voxel, the measured valuef i is the sum of the true value f i and additive noise ν i : All ν i are assumed independent and identically distributed (i.i.d.), but their distribution is unknown. If we did know the true underlying cumulative distribution function F of the noise, we could choose a confidence parameter θ ∈ (0, 1) and use the cumulative distribution function to calculate ν θ/2 , ν 1−θ/2 such that 1 Let us instead proceed non-parametrically, and divide the whole image into two groups of voxels -the ones belonging to the background region and the rest. For simplicity, let the indices i = 1, . . . , k, (k < n), stand for the background voxels. In this region, we have f i = 0 andf i = ν i . Therefore, the background region provides us with a number of independent samples from the unknown distribution of ν. The Dvoretzky-Kiefer-Wolfowitz inequality [18,38,35] states that for the empirical cumulative distribution function Therefore, computing ν θ/2 and ν 1−θ/2 such that we also have We may therefore use the valueŝ as lower and upper bounds for the true values f i outside the background region. The Dvoretzky-Kiefer-Wolfowitz inequality implies that the interval estimates converge to the true intervals, determined by (5.1), as the number of background pixels k increases with the image size n. This procedure, with large k, will therefore provide an estimate of a single-experiment (m = 1) confidence interval for f i . We note that this procedure will, however, not yield the convergence of the interval estimate [f l,i ,f u,i ] to the true data; for that we would need multiple experiments, i.e., multiple sample images (m > 1), not just agglomeration of the background voxels into a single noise distribution estimate. In practice, however, we can only afford a single experiment (m = 1), and cannot go to the limit.

Verification of the approach with synthetic data
To verify the effectiveness of the considered approach and to compare it to the other models, we use synthetic data. For the ground-truth tensor field u g.t. we take a helix region in a 3D box 100 × 100 × 30, and choose the tensor in each point inside the helix in such a way that the principal eigenvector coincides with the helix direction ( Figure 2). The helix region is described by the following equations: x = (R + r cos(θ)) cos(φ), y = (R + r cos(θ)) sin(φ), The vector direction in every point coincides with helix direction: We take the parameters R = 0.3, φ max = 4π, r max = 0.07 in this numerical experiment.
We apply several models for solving the inverse problem of reconstructing u: the linear and non-linear L 2 approaches (3.3) and (3.2), and the constrained problem (3.5). As the regulariser we use R = TGV 2 (0.9α,α) , where the choice β = 0.9α was made somewhat arbitrarily, however yielding good results for all the models. This is slightly lower than the range [1, 1.5]α discovered in comprehensive experiments for other imaging modalities [7,16].
For the linear and non-linear L 2 models (3.3) and (3.2), respectively, the regularisation parameter α is chosen either by a version of the discrepancy principle for inconsistent problems [49] or optimally with regard to the · F,2 distance between the solution and the ground-truth. In case of the discrepancy principle, such an α was chosen that the following equality holds: We find α by solving this equation numerically using bisection method. We start by finding such α 1 , α 2 that ∆ρ(α 1 ) > 0 and ∆ρ(α 2 ) < 0. We calculate ∆ρ(α 3 ), α 3 = α 1 +α 2 2 and depending on its sign replace either α 1 or α 2 with α 3 . We repeat this procedure until the stopping criteria is reached.  Figure 2: Helix vector field for the principal eigenvector of the ground-truth tensor field As stopping criteria we use |f (α)| < . We use τ = 1.05, = 0.01 for linear and τ = 1.2, = 0.0001 for non-linear L 2 solution. A value of τ yielding a reasonable degree of smoothness has been chosen by trial and error, and is different for the non-linear model, reflecting a different non-linear objective in the discrepancy principle. For the constrained problem we calculate θ = 90%, 95%, and 99% confidence intervals to generate the upper and lower bounds. We however digress a little bit from the approach of Section 2.2. Minding that we do not know the true underlying distribution, which fails to be Rician as illustrated in Figure 1, we do not use it to calculate the confidence intervals, but use the estimation procedure described in Section 5.1. We stress that we only have a single sample of each signal s j , so are unable to verify any asymptotic estimation properties.
The numerical results are in Table 1 and Figures 3-5, with the first of the figures showing the colour-coded principal eigenvector of the reconstruction, the second showing the fractional anisotropy and principal eigenvectors, and the last one the errors in the latter two, in a colour-coded manner. All plots are masked to represent only the non-zero region. The field of Table 1: Numerical results for the synthetic data. For the linear and nonlinear L 2 the free parameter chosen by the parameter choice criterion is the regularisation parameter α, and for the constrained problem it is the confidence interval.

Method
Parameter fractional anisotropy is defined for a field u of 2-tensors on Ω ⊂ R m as with λ 1 , . . . , λ m denoting the eigenvalues of u(x). It measures how far the ellipsoid prescribed by the eigenvalues and eigenvectors is from a sphere, with FA u (x) = 1 corresponding a full sphere, and FA u (x) = 0 corresponding to a degenerate object not having full dimension. As we can see, the non-linear approach (3.2) performs overall the best by a wide margin, in terms of the pointwise Frobenius error, i.e., error in · F,2 . This is expressed as a PSNR in Table 1. What is, however, interesting, is that the constraint-based approach (3.5) has a much better reconstruction of the principal eigenvector angle, and a comparable reconstruction of its magnitude. Indeed, the 95% confidence interval in Figure 3(g) and Figure  4(g) suggests a nearly perfect reconstruction in terms of smoothness. But, the Frobenius PSNR in Table 1 for this approach is worse than the simple unregularised inversion by regression. The problem is revealed by Figure  5(f): the large white cloudy areas indicate huge fractional anisotropy errors, while at the same time, the principal eigenvector angle errors expressed in colour are much lower than for other approaches. Good reconstruction of the principal eigenvector is important for the process of tractography, i.e., the reconstruction of neural pathways in a brain. One explanation for our good results is that the regulariser completely governs the solution in areas where the error bounds are inactive due to generally low errors. This results in very smooth reconstructions, which is in the present case desirable as our synthetic tensor field is also smooth within the helix.

Results with in vivo brain imaging data
We now wish to study the proposed regularisation model on a real invivo diffusion tensor image. Our data is that of a human brain, with the measurements of a volunteer performed on a clinical 3T system (Siemens Magnetom TIM Trio, Erlangen, Germany), with a 32 channel head coil. A 2D diffusion weighted single shot EPI sequence with diffusion sensitising gradients applied in 12 independent directions (b = 1000s/mm 2 ). An additional reference scan without diffusion was used with the parameters: TR = 7900ms, TE = 94ms, flip angle 90 • . Each slice of the 3D data set has plane resolution 1.95mm × 1.95mm, with a total of 128 × 128 pixels. The total number of slices is 60 with a slice thickness of 2mm. The data set consists of 4 repeated measurements. The GRAPPA acceleration factor is 2.
Prior to the reconstruction of the diffusion tensor, eddy current correction was performed with FSL [47]. Written informed consent was obtained from the volunteer before the examination. For error bounds calculation according to the procedure of Section 5.1, to avoid systematic bias near the brain, we only use about 0.6% of the total volume near the borders, or roughly k ≈ 6000 voxels.
To estimate errors for the all the considered reconstruction models, for each gradient direction b i we use only one out of the four duplicate measurements. We then calculate the errors using a somewhat less than ideal pseudo-ground-truth, which is the linear regression reconstruction from all the available measurements.
The results are in Table 2 and Figures 6-8, again with the first of the figures showing the colour-coded principal eigenvector of the reconstruc- Table 2: Numerical results for the in-vivo brain data data. For the L 2 and non-linear L 2 reconstruction models the free parameter chosen by the parameter choice criterion is the regularisation parameter α, and for the constrained problem it is the confidence interval.  Table 2. This time, the linear L 2 approach (3.3) has best overall reconstruction (Frobenius PSNR), while the nonlinear L 2 approach (3.2) has clearly the best principal eigenvector angle reconstruction besides the regression, which does not seem entirely reliable regarding our regression-based pseudo-groundtruth. The constraints based approach (3.5), with 95% confidence intervals is, however, not far behind in terms of numbers. More detailed study of the corpus callosum in Figure 8 (small picture in picture) and Figure 7 however indicates a better reconstruction of this important region by the nonlinear approach. The constrained approach has some very short vectors there in the white region. Naturally, however, these results on the in vivo data should be taken with a grain of salt, as we have only a somewhat unreliable pseudo-ground-truth available for comparison purposes.

Conclusions from the numerical experiments
Our conclusion is that the error bounds based approach is a feasible alternative to standard modelling with incorrect Gaussian assumptions. It can produce good reconstructions, although the non-linear L 2 approach of [52] is possibly slightly more reliable. The latter does, however, in principle depend on a good initialisation of the optimisation method, unlike the convex bounds based approach.
Further theoretical work will be undertaken to extend the partial-orderbased approach to modelling errors in linear operators to the non-lattice case of the semidefinite partial order for symmetric matrices, which will allow us to consider problems of diffusion MRI with errors in the forward operator.
It shall also be investigated whether the error bounds approach needs to be combined with an alternative, novel, regulariser that would ameliorate the fractional anisotropy errors that the approach exhibits. It is important to note, however, that from the practical point of view, of using the reconstruction tensor field for basic tractography methods based solely on principal eigenvectors, these are not that critical. As pointed out by one of the reviewers, the situation could differ with more recent geodesic tractography methods [?, ?, ?] employing the full tensor. We provide basic principal eigenvector tractography results for reference in Figure 9, without attempting to extensively interpret the results. It suffices to say that the results look comparable. With this in sight, the error bounds approach produces a very good reconstruction of the direction of the principal eigenvectors, although we saw some problems with the magnitude within the corpus callosum. ages, and Florian Knoll for many inspiring discussions.

A data statement for the EPSRC
The MRI scans used for this publication have been used by the courtesy of the producer. Data that is similar for all intents and purposes, is widespread, and can be easily produced by making a measurement of a human subject with an MRI scanner. Our source codes are archived at https://www. repository.cam.ac.uk/handle/1810/253422.

Appendix: Notation and techniques
We recall some basic, not completely standard, mathematical notation and concepts in this appendix. We begin with partially ordered vector spaces, following the book [43]. Then we proceed to tensor calculus and tensor fields of bounded variation and of bounded deformation. These are also covered in more detail for the diffusion tensor imaging application in [54].

Banach lattices
A linear space X, endowed with a partial order relation is called an ordered vector space if the partial order agrees with linear operations in the following way: x y =⇒ x + z y + z ∀x, y, z ∈ X, x y =⇒ λx λy ∀x, y ∈ X and λ ∈ R + .
An ordered vector space is called a vector lattice if each pair of elements x, y ∈ X have a supremum x ∨ y ∈ X and infimum x ∧ y ∈ X. Supremum of two elements x, y of a Banach lattice X is the element z = x ∨ y with the following properties: z x, z y and ∀z ∈ X such thatz x andz y we havez z.
For any x ∈ X, the element x + = x ∨ 0 is called its positive part, the element x − = (−x)∨0 = (−x) + is called its negative part, the element |x| = x + + x − is its absolute value. The equalities x = x + − x − and |x| = x ∨ (−x) hold for any x ∈ X.
It is obvious that suprema and infima exist for any finite number of elements of a vector lattice. A vector lattice X is said to be order complete if any bounded from above set in X has a supremum.
Let X and Y be ordered vector spaces. A linear operator U : X → Y is called positive, if x X 0 implies U x Y 0. An operator U is called regular, if it can be written as U = U 1 − U 2 , where U 1 and U 2 are positive operators.
Denote the linear space of all regular operators X → Y by L ∼ (X, Y ). A partial order can be introduced in L ∼ (X, Y ) in the following way: U 1 U 2 , if U 1 −U 2 is a positive operator. If X and Y are vector lattices and Y is order complete, then L ∼ (X, Y ) is also an order complete vector lattice.
A norm · defined in a vector lattice X is called monotone if |x| |y| implies x y . A vector lattice endowed with a monotone norm is called a Banach lattice if it is norm complete. If X and Y are Banach lattices, then all operators in L ∼ (X, Y ) are continuous.
Let us list some examples of Banach lattices. The space of continuous functions C(Ω), where Ω ⊂ R n , is a Banach lattice under the natural pointwise ordering: f C g if and only if f (x) g(x) for all x ∈ Ω. The spaces L p (Ω), 1 p ∞, are also Banach lattices under the following partial ordering: f L p g if and only if f (x) g(x) almost everywhere in Ω. With this partial order, L p (Ω), 1 p ∞, are order complete Banach lattices. The Banach lattice of continuous functions C(Ω) is not order complete.

Tensors in the Euclidean setting
We call a k-linear mapping A : R m × · · · × R m → R a k-tensor, denoted A ∈ T k (R m ). This is a simplification from the full differential-geometric definition, sufficient for our finite-dimensional setting. We say that A is symmetric, denoted A ∈ Sym k (R m ), if it satisfies for any permutation π of {1, . . . , k} that A(c π1 , . . . , c πk ) = A(c 1 , . . . , c k ).
With e 1 , . . . , e m the standard basis of R m , we define on T k (R m ) the inner product A, B := p∈{1,...,m} k A(e p 1 , . . . , e p k )B(e p 1 , . . . , e p k ), and the Frobenius norm The Frobenius norm is rotationally invariant in a sense crucial for DTI. We refer to [54] for a detailed discussion of this, as well of alternative rotationally invariant norms.

Example .1 (Vectors). Vectors
A ∈ R m are of course symmetric 1-tensors, The inner product is the usual inner product in R m , and the Frobenius norm is the two-norm, A F = A 2 .
Example .2 (Matrices). Matrices are 2-tensors: A(x, y) = Ax, y , while symmetric matrices A = A T are symmetric 2-tensors. The inner product is A, B = i,j A ij B ij and A F is the matrix Frobenius norm. We use the notation A ≥ 0 for positive-semidefinite matrices A. One can verify that this relation indeed defines a partial order in the space of symmetric matrices: With this partial order, the space of all symmetric matrices becomes an ordered vector space, but not a vector lattice. However, it enjoys some properties similar to those of vector lattices: for example, any directed upwards subset 2 has a supremum [36,Ch.8].

Symmetric tensor fields of bounded deformation
Let u : Ω → Sym k (R m ) for a domain Ω ⊂ R m . We set and u F,∞ := ess sup x∈Ω u(x) F , The spaces L p (Ω; Sym k (R m )) are defined in the natural way using these norms, and clearly satisfy all the usual properties of L p spaces.
In the particular case of matrices (k = 2), partial order can be introduced in the space L p (Ω; Sym 2 (R m )) in the following way: u v iff u(x) ≥ v(x) almost everywhere in Ω. (.5) Recall that the symbol ≥ stands for the positive semidefinite order (.4) in the space of symmetric matrices. Since the positive semidefinite order is not a lattice, neither is the order (.5).
Example .3. The space BD(Ω; Sym 0 (R m )) agrees with the space BV(Ω) of functions of bounded variation. The space BD(Ω; Sym 1 (R m )) = BD(Ω) is the space of functions of bounded deformation of [48].