In this paper, we propose a variational regularization method for denoising and inpainting of diffusion tensor magnetic resonance images. We consider these images as manifold-valued Sobolev functions, i.e. in an infinite dimensional setting, which are defined appropriately. The regularization functionals are defined as double integrals, which are equivalent to Sobolev semi-norms in the Euclidean setting. We extend the analysis of  concerning stability and convergence of the variational regularization methods by a uniqueness result, apply them to diffusion tensor processing, and validate our model in numerical examples with synthetic and real data.
In this paper, we investigate denoising and inpainting of diffusion tensor (magnetic resonance) images (DTMRI) with a derivative-free, non-local variational regularization technique proposed, implemented and analyzed first in .
The proposed regularization functionals generalize equivalent definitions of the Sobolev semi-norms, which have been derived in the fundamental work of  and follow-up work [16, 40]. These papers provide a derivative-free representation of Sobolev semi-norms for intensity and vector-valued functions. The beauty of this representation is that it allows for a straightforward definition of Sobolev energies of manifold-valued data (see ), without requiring difficult differential geometric concepts (as for instance in [21, 15]).
Diffusion tensor images are considered to be re-presentable as functions from an image domain , with , respectively, into the manifold of symmetric, positive definite matrices in , denoted by 𝐾 in the following – for DTMRI images . Therefore, they are ideal objects to check the efficiency of the proposed regularization techniques. A measured diffusion tensor is often very noisy, and post-processing steps for noise removal are important. Even more, due to the noise, it is possible that a measured tensor has negative eigenvalues, which is not physical, and thus often the whole tensor at this point is omitted, leading to incomplete data. Then the missing information has to be inpainted before visualization.
Variational regularization of vector and matrix manifold-valued functions has been considered before, for instance in [49, 44, 46, 5, 8, 30, 50, 36] and [9, 12, 37, 34]. Non-local regularization formulations were studied for example in [10, 28] for filtering tensor-valued functions, and also [22, 29] for filtering intensity images. An overview of diffusion and regularization techniques for vector- and matrix-valued data is given in .
Variational methods for denoising and inpainting attempt to find a good compromise between matching some given noisy, tensor-valued data and prior information on the desired solution , also called noise-free or ideal solution. The choice of prior knowledge on is
that it is an element of the set , with a metric 𝑑 on 𝐾, the set of positive definite matrices, which is a subset of the fractional Sobolev space , with and ,
is relatively small. The function 𝜌 is a non-negative and radially symmetric mollifier with an on-off indicator denoting whether the mollifier is used or not. Note that, in case that is the Euclidean metric and if we choose in addition , becomes the fractional Sobolev semi-norm.
where the parameter determines the preference of staying close to the given function in and a small energy . One should not confuse the energy term with double integral representations approximating semi-norms on manifolds (see for instance [24, 27, 19]).
The indicator function of ,
used in (1.2) allows us to consider the two tasks of denoising ( ) and inpainting ( ) within one analysis.
The paper is organized as follows. In Section 2, we constitute our notation and setting used to analyze variational methods for DTMRI data processing. We review regularization results from  in Section 3. In Section 4, we verify that these results from Section 3 are applicable in the context of diffusion tensor imaging, meaning that we show that the functional defined in (1.2) attains a minimizer and fulfills a stability as well as a convergence result. Furthermore, we extend the analysis and give a uniqueness result using differential geometric properties of symmetric, positive definite matrices, where it is of particular importance that these matrices endowed with the log-Euclidean metric form a flat Hadamard manifold. In Section 5, we give more details on the numerical minimization of the regularization functional, and discuss different variants. In the last section, Section 6, we show numerical results for denoising and inpainting problems of synthetic and real DTMRI data.
2 Notation and setting
In the beginning, we summarize basic notation and assumptions used throughout the paper. In the theoretical part, we work with general dimensions , while we consider the particular case , , that is two-dimensional slices of a three-dimensional DTMRI image, in the numerical examples in Section 6.
We assume the following.
is a nonempty, bounded and connected open set with Lipschitz boundary, and is measurable.
, and .
is a nonempty and closed subset of .
denotes the Euclidean distance induced by the (Frobenius norm) on .
denotes an arbitrary metric on 𝐾 which is equivalent to .
Moreover, we need the definition of a mollifier which appears in the regularizer of the functional in (1.2).
Definition 2.2 (Mollifier)
We call a mollifier if
𝜌 is a non-negative, radially symmetric function,
there exists some and such that .
The last condition holds for instance if 𝜌 is radially decreasing satisfying .
2.1 Diffusion tensors
It is commonly assumed that the recorded diffusion tensor images are functions with values which are symmetric, positive definite matrices. Hence we make the assumption that
where is the set of symmetric, positive definite, real matrices defined below in (2.1). When working with data from MRI measurements, .
In the following definition, we summarize sets of matrices and associated norms on the sets.
The vector space of symmetric matrices is
Additionally, we define the set of symmetric, positive definite matrices
The set of symmetric, positive definite matrices with bounded spectrum is
where denotes the spectrum of a given matrix. For diffusion tensors, the spectrum is real.
The set of symmetric, positive definite matrices with bounded logarithm is
From now on and whenever possible, we omit the space dimension and write and instead of and .
2.2 Fractional Sobolev spaces
Moreover, we need the definition of fractional Sobolev spaces and associated subsets.
Definition 2.4 (Sobolev spaces of fractional order)
Let Assumption 2.1 hold.
We denote by the Lebesgue space of matrix-valued functions.
The Sobolev space consists of all weakly differentiable functions in for which
where is the Jacobian of 𝑤 and is the Sobolev semi-norm.
The fractional Sobolev space of order 𝑠 is defined (cf. ) as the set
equipped with the norm
where is the semi-norm on , defined by
We define the fractional Sobolev set of order 𝑠 with data in 𝐾 as(2.5)
The Lebesgue set with data in 𝐾 is defined as
Note that and are sets and not linear spaces because summation of elements in 𝐾 is typically not closed in 𝐾.
3 Metric double integral regularization on closed subsets
We start this section by stating conditions under which the regularization functional in (1.2) attains a minimizer and fulfills a stability as well as a convergence result. Therefore, we recall results established in . There the authors define a regularization functional inspired by the work of Bourgain, Brézis and Mironescu [11, 40, 16]. The analysis in turn is based on . We apply these results to diffusion tensor image denoising and inpainting in the next section.
We start by stating general conditions on the exact data , the noisy data and the functional , defined in (1.2).
for every and , the level sets
are weakly sequentially pre-compact in .
There exists such that is nonempty.
If Assumption 2.1 is fulfilled and in particular when performing image denoising ( ) or inpainting ( ) of functions with values in 𝐾, then the functional (1.2) with as in (1.1) defined on satisfies Assumption 3.1 (cf. ).
Existence: For every and , the functional attains a minimizer in .
Stability: Let be fixed, , and let be a sequence in such that
Then every sequence satisfying
has a converging subsequence with respect to the weak topology of . The limit of every such converging subsequence is a minimizer of . Moreover, converges to .
Convergence: Let be a function satisfying and for . Let be a sequence of positive real numbers converging to 0. Moreover, let be a sequence in with , and set . Then every sequence defined as
has a weakly converging subsequence as (with respect to the topology of ). In addition, . Moreover, it follows that even weakly (with respect to the topology of ) and .
In the theorem above, stability (with respect to the -norm) ensures that the minimizers of depend continuously on the given data . We emphasize that, in an Euclidean setting (that is on , for and , one could expect convergence in even stronger norms. However, here, we have to make sure that the traces into are well-defined in appropriate Sobolev spaces, which requires additional compactness assumptions, or in other words, stronger regularization.
In the next section, we apply Theorem 3.3 to diffusion tensor images, i.e. when choosing 𝐾 as a closed subset of the symmetric, positive definite matrices.
4 Diffusion tensor regularization
The goal of this section is to define appropriate fractional order Sobolev sets as defined in (2.5) of functions which can represent diffusion tensor images. To this end, we use the set of symmetric, positive definite matrices with bounded logarithm (defined in (2.3))
Below, we show that Theorem 3.3 applies to the regularization functional in (1.2) with the particular choice . In addition to what follows from the general theory from  in a straightforward manner, we present a uniqueness result for the minimizer of the regularization functional.
We begin by defining needed concepts from matrix calculus. When working with symmetric, positive definite matrices, many of the operations below reduce to their scalar counterpart applied to the eigenvalues.
4.1 Matrix calculus
We start this section by repeating basic definitions known from matrix calculus (see for instance ). Especially the matrix logarithm is needed to define the log-Euclidean metric on the symmetric, positive definite matrices.
Lemma 4.1 (Matrix properties)
The following statements hold.
Eigendecomposition: Let with eigenvalues . Then we can write
where is the orthonormal matrix whose 𝑖-th column consists of the 𝑖-th normalized eigenvector of 𝐴. Hence we have that , where denotes the identity matrix in . Further, Λ is the diagonal matrix whose diagonal entries are the corresponding eigenvalues, .
If are both unitary, then so are , , and .
Let with corresponding eigendecompositions and , where unitary and , diagonal.
Exponential map: The exponential map is defined as
It holds that
where denotes the (scalar) exponential function, and is a diffeomorphism [3, Theorem 2.8].
Logarithm: If , then A is the matrix logarithm of 𝐵. It is defined as
where is the (scalar) natural logarithm, i.e. .
When restricting to symmetric, positive definite matrices, is a diffeomorphism [3, Theorem 2.8].
The previous definition, Definition 4.2, shows that the exponential and logarithm of a symmetric (positive definite) matrix can be computed easily due to the eigendecomposition (see Remark 4.1) by calculating the scalar exponential map and logarithm of the eigenvalues.
Remark 4.3 (Matrix logarithm)
For a general matrix in , the matrix logarithm is not unique. Matrices with positive eigenvalues have a unique, symmetric logarithm, called the principal logarithm .
The next lemma states properties of the Frobenius norm (recall (2.4)).
Lemma 4.4 (Properties of Frobenius norm)
We have the following properties.
Let be symmetric and skew-symmetric, respectively, i.e. , . Then(4.2)
The Frobenius norm is unitary invariant, i.e.(4.3)
for and unitary.
If with (positive) eigenvalues , then(4.4)
The last lemma of this subsection deals with the set , the set of symmetric, positive definite matrices with bounded spectrum in the interval , defined in (2.2). We need this result later in the numerical implementation for defining a suitable projection.
Given an arbitrary matrix , there always exists a unique matrix which is closest in the Frobenius norm, i.e.
The minimizing matrix 𝑀 can be computed explicitly as stated in the following lemma. The proof is done in a way similar to [25, Theorem 2.1] and is included here for completeness.
Let . Define and as the symmetric and skew-symmetric parts of 𝐴, respectively. Let be the eigenvalues of 𝐵 which can be decomposed into , where 𝑍 is a unitary matrix, i.e. , and . Then the unique minimizer of
where is defined in (2.2), is with , where
By definition of 𝐵 and 𝐶, we can write and thus
The matrix 𝐵 is symmetric, and thus we can write , where is a unitary matrix whose columns are the eigenvectors of 𝐵 and is a diagonal matrix whose entries are the eigenvalues of 𝐵, i.e. . Let be similar to 𝑋 so that . Then we obtain, by using (4.3),
Thus the lower bound is uniquely attained for with
After giving the needed definitions from matrix calculus, the goal of this subsection is now to apply Theorem 3.3 to the regularization functional defined in (1.2) with the set defined in (4.1) and associated log-Euclidean metric defined below in (4.6). Therefore, we need to prove the equivalence of the log-Euclidean and Euclidean metric to guarantee in particular that Assumption 2.1 (v) is fulfilled. Then Assumption 3.1 holds true as stated in Remark 3.2, and therefore Theorem 3.3 is applicable.
We start by defining and reviewing some properties of the log-Euclidean metric.
Definition 4.6 (Log-Euclidean metric)
Let . The log-Euclidean metric is defined as
The log-Euclidean metric satisfies the metric axioms on .
This follows directly because is a norm and restricted to is a diffeomorphism. ∎
The reasons for choosing this measure of distance is stated in the following remark.
The log-Euclidean metric arises when considering not just as convex cone in the vector space of matrices but as a Riemannian manifold. Thus it can be endowed with a Riemannian metric defined by an inner product on the tangent space, see for example [18, 38, 3]. Two widely used geodesic distances are the affine-invariant metric
and the log-Euclidean metric as stated above. These measures of dissimilarity are more adequate in DTMRI as pointed out in  because zero or negative eigenvalues induce an infinite distance.
The affine-invariant distance measure is computationally much more demanding, which is a major drawback. This is not the case for the log-Euclidean distance, which leads to Euclidean distance computations in the matrix logarithmic domain.
The following statement can be found in [33, Sections 2.4 & 2.5].
Let denote the log-Euclidean metric (as defined in (4.6)), the affine-invariant metric (as defined in (4.7)) and the standard Euclidean distance. Then as well as form a complete metric space. This is not the case for .
Different metrics induce different properties on the corresponding regularizer. We compare and additionally a Sobolev semi-norm regularizer in the following.
Remark 4.10 (Invariances)
For the log-Euclidean metric as defined in (4.6), the following holds true.
Scale invariance: Let and , and denote by the identity matrix in . Then
Invariance under inversion: Let . Because , we directly get that
Unitary invariance: Let and 𝑈 unitary. Because of the unitary invariance of the Frobenius norm, then
Although we only work with fractional derivatives of order , we consider for comparison purposes the regularization functional (see also (6.1))
None of the invariances above, i.e. scale invariance, invariance under inversion and unitary invariance, is valid for Θ.
for some constant matrix , i.e. it is translation and reflection invariant. This, in turn, does not hold (or is not even well-defined) for our regularizer with the log-Euclidean metric but as well when considering the standard Euclidean distance, i.e. it does hold for . A comparison is shown in Table 1.
In order to show that Theorem 3.3 is applicable for defined in (1.2) with and associated log-Euclidean metric defined in (4.6), we have to show that Assumption 3.1 and therefore Assumption 2.1, in particular the equivalence stated in (v), is valid. In order to prove that, we need the following corollary.
Let (defined in (2.3)) with eigenvalues . Then, for each ,
i.e. (for the definition of the latter set, see (2.2)).
If , it holds that . Using (4.4), this is equivalent to , so the claim follows. ∎
Note that the reverse embedding in the previous lemma does not hold true. If such that, for each eigenvalue , we have that , then .
Now we can prove that the Euclidean and the log-Euclidean metric are equivalent on . In particular, we show that is bi-Lipschitz on and calculate the constant explicitly. Without explicit computation, this would follow from the fact that is a diffeomorphism on symmetric, positive definite matrices and that is a compact subset.
Let be defined in (2.3). Then
where are unitary matrices and are diagonal matrices whose entries are the corresponding positive eigenvalues of 𝐴 and of 𝐵, respectively. By Corollary 4.11, it holds that for all . We consider two cases.
Case 1. We assume that all eigenvalues of 𝐴 and 𝐵 are equal, i.e. they have the same one-dimensional spectrum , meaning that . This in turn gives that
Case 2. We now assume that there exist at least two different eigenvalues , of 𝐴 and 𝐵.
We show the lower inequality in (4.8). The upper inequality can be done analogously.
where . Using the definition of the Frobenius norm in (2.4), we obtain further that
Indices for which do not contribute to the sum in (4.10) (and do not change the following calculation), so we define as the set of such indices for which we have .
From the mean value theorem, it follows that, for every , there exists some
Further, we can write
which finishes the proof. ∎
So far we showed that the functional as defined in (1.2) over using the log-Euclidean metric as in (4.6) attains a minimizer. In this subsection, we prove that the minimum is unique. To this end, we consider the symmetric, positive definite matrices from a differential geometric point of view.
The space , where denotes the log-Euclidean metric as defined in (4.6), is a complete, connected Riemannian manifold with zero sectional curvature.
In other words, is a flat Hadamard manifold and therefore in particular a Hadamard space. The last property guarantees that the metric 𝑑 is geodesically convex [48, Corollary 2.5], i.e. let be two geodesics; then
The following result states that connecting geodesics between two points in stay in this set.
as a connecting geodesic between , and and
as the evaluation of the geodesic between and at time 𝑡 for . Then .
We split the proof into two parts. First we show that maps indeed into . Afterwards, we prove that it actually lies in .
which is equivalent to
because , i.e. . This shows that maps into .
Next we need to prove that actually , i.e. that
We obtain by Jensen’s inequality that
Using (4.8), it follows that
The last expression is finite because of the assumption that . ∎
Now we can state the uniqueness result.
Existence of a minimizer is guaranteed by Theorem 4.13. Now, let us assume that there exist two minimizers of the functional .
Because and are both minimizers, we have that
In particular, for , we obtain by the above equality and by (4.15) that
which is a contradiction to the minimizing property of and for . Hence and must be equal forcing equality in (4.15) and thus giving that the minimum is unique. ∎
Existence and uniqueness in the case
If , then existence and uniqueness of the minimizer of the functional even holds on the larger set rather than on , where is associated with the log-Euclidean distance as defined in (4.6). Existence in Theorem 4.13 and uniqueness in Theorem 4.16 (with ) are based on the theory provided in  (see Theorem 3.3), where it is a necessary assumption that the set 𝐾 is closed which is not the case for the set .
Nevertheless, it is possible to get existence and uniqueness on this set because, for every minimizing sequence , , we automatically get that so that it takes values on the closed subset . Then existence of a unique minimizer on follows by the proofs already given, see [14, Theorem 3.6] and Theorem 4.16.
We now sketch the proof of the assertion. Throughout this sketch, 𝐶 denotes a finite generic constant which, however, can be different from line to line.
Sketch of assertion
Denote by the log-Euclidean metric (as defined in (4.6)). Let us take a minimizing sequence , , of so that we can assume that for all , .
Computing the log-Euclidean metric leads to evaluations of the Euclidean metric in the matrix logarithmic domain, cf. (4.6), meaning that
Because of this and the fact that , we get that
Because of [32, Lemma 2.7], we can thus bound the -norm of ,
yielding in particular that .
By the definition of in (2.3), we thus obtain that for all . Hence every minimizing sequence , , of is automatically a minimizing sequence in .
In this section, we go into more detail on the minimization of the regularization functional defined in (1.2) with the log-Euclidean metric as defined in (4.6) (see ) over the set , the set of symmetric, positive definite matrices with bounded logarithm, as defined in (4.1) for denoising and inpainting of DTMRI images.
To optimize , we use a projected gradient descent algorithm. The implementation is done in Matlab. The gradient step is performed by using Matlab’s built-in function fminunc, where the gradient is approximated with a finite difference scheme (central differences in the interior and one-sided differences at the boundary). Therefore, after each step, we project the data which are elements of the larger space back onto by applying the following projection. We remark that, by projection, we here refer to an idempotent mapping.
Definition 5.1 (Projection operators)
We define the following projections.
Projection of onto : Let be a symmetric matrix with eigendecomposition with . Then the projection of 𝑀 onto the set is given by(5.1)
Projection of onto : Let with eigenvalues and eigendecomposition . Define
as the squared Frobenius norm of . Then the pr