Denoising of diffusion magnetic resonance images using a modified and differentiable Monge–Kantorovich distance for measure-valued functions

Abstract

We report on the implementation of a novel total-variation denoising method for diffusion spectrum images (DSI). Our method works on the raw signal obtained from dMRI. From the Stejskal-Tanner equation [6], the signals S(x, s_k), 1 ≤ k ≤ K, at a given voxel location x may be considered as samplings of a measure supported on the unit sphere ${\mathbb{S}}^2\in {\mathbb{R}}^3$ at locations $s_k=({\theta }_k,{\varphi }_k)\in {\mathbb{S}}^2$ which quantify the ease/difficulty of diffusion in these directions. We consider the entire signal S as a measure-valued function in a complete metric space which employs the Monge–Kantorovich (MK) metric. A total variation (TV) for measures and measure-valued functions is also defined. A major advance in this paper is the use of a modification of the standard MK distance which permits rapid computation in higher dimensions. An added bonus is that this modified metric is differentiable. The resulting TV-based denoising problem is a convex optimization problem.

Introduction

In [11], a novel framework for the representation of images obtained from diffusion magnetic resonance imaging (dMRI) [8] and their denoising using total variation was presented. In this framework, the (digital) images obtained from dMRI are represented by measure-valued functions (MVFs): If we let X ⊂ Rⁿ denote the base space, i.e., the support of an image (e.g., a human brain), then the dMR image is a mapping $\mu :X\to \mathcal{M}\left(Y\right),$ where $Y\subset {\mathbb{R}}^m,$ $m=1,2$ or 3, denotes an appropriate range space and $\mathcal{M}\left(Y\right)$ the set of Borel probability measures on Y. For example, if $Y={\mathbb{S}}^2=\{u\in {\mathbb{R}}^3{|\parallel u\parallel }_2=1\},$ the unit sphere in ${\mathbb{R}}^3,$ then the measure $\mu \left(x\right)\in \mathcal{M}\left({\mathbb{S}}^2\right)$ can characterize the relative ease or difficulty of diffusion of water molecules in the voxel located at x ∈ X in the direction $s=\Omega =(\theta ,\varphi )\in {\mathbb{S}}^2$.

Let us contrast this framework with the conventional way in which images from dMRI, in particular, high angular diffusion resolution imaging (HARDI), are represented, namely, as vector-valued functions, i.e., $u\left(x\right)=(u_1\left(x\right),u_2\left(x\right),\cdots u_K\left(x\right)),$ composed of signals, u_k(x), which are obtained by aligning the linear diffusion gradient vector with gradient $g\in {\mathbb{R}}^3$ in K different directions, $s_k=({\theta }_k,{\varphi }_k)\in S^2,$ 1 ≤ k ≤ K. Most image processing algorithms simply work with such vector-valued functions in a standard manner, treating them as “cubes”. What we are proposing is to work with dMRI signals in a manner that attempts to transcend such approaches in which the net signal u is simply regarded as a “stack” of component signals u_k.

In an even earlier paper [13], we examined the use of function-valued mappings (FVMs) for HARDI. In this case, a HARDI signal is represented by the function-valued mapping $u:X\to L^2\left({\mathbb{S}}^2\right)$. In all of these works – including this paper – our philosophy has been, as described in our more recent paper on FVMs and their use in hyperspectral imaging [16], to investigate if “vector-valued images be better understood, and perhaps better algorithms be developed, if their range were a space of continuously defined functions instead of ${\mathbb{R}}^N$.” Indeed, with continued improvements in technology, the number of components K in dMRI, as well as in hyperspectral imaging for remote sensing [1], is steadily increasing, essentially approaching the “continuum limit” of FVM and MVF.

In [11], our total variation-based denoising problem was formulated as follows: Given a noisy (measure-valued) image $\tilde{\mu },$ find a solution to the following optimization problem,

$$\underset{\nu \in \mathcal{F}\left(X\right)}{\min }{d}_{\mathcal{F}\left(X\right)}(\tilde{\mu },\nu )+{\parallel \nu \parallel }_{TV}.$$

Here, $\mathcal{F}\left(X\right)$ denotes an appropriate space of measure-valued functions suppported on X, with metric $d_{\mathcal{F}}\left(X\right),$ and ‖ν‖_TV denotes the total variation of the measure-valued function $\nu \in \mathcal{F}\left(X\right)$. Both of these will be reviewed in the next section.

As discussed in [11], the metric $d_{\mathcal{F}\left(X\right)}$ requires an appropriate metric between measures in the space $\mathcal{M}\left(Y\right)$. A natural choice is the so-called Monge-Kantorovich (MK) metric for measures, also to be reviewed below. Unfortunately, the computation of the traditional MK distance between (discrete) probability measures is computationally inexpensive only in one dimension, i.e., measures on $\mathbb{R}$. The determination of efficient algorithms to compute MK distances in ${\mathbb{R}}^2$ or ${\mathbb{S}}^2$ and ${\mathbb{R}}^3$ is still an open problem. In [11], we showed how to solve the optimization problem of denoising using total variation minimization for the rather artificial – but still mathematically interesting – one-dimensional case, exploiting the fact that the MK distance between two one-dimensional measures is the L¹ distance between their respective cumulative distribution functions. For the much more difficult higher-dimensional case, as we wrote in [11], “an alternative is to replace the Monge-Kantorovich metric on measures with another metric.” This is indeed one of the major accomplishments of this paper, representing a major advance over [11]: We propose a modification of the MK metric which is computationally inexpensive not only in one dimension but also in higher-dimensional spaces. An added bonus is that this modified metric is differentiable. The resulting total variation denoising scheme for measure-valued functions becomes a convex optimization problem which can be applied to real data, as we show below.

In closing this section, we review a few basic ideas from dMRI which provide the basis of our measure-valued approach. The reader may recall that in standard magnetic resonance imaging (MRI), a linear three-dimensional magnetic field is used so that the MRI signal is a Fourier transform. Inversion of this signal produces an “image” of the object being examined. In dMRI, a linear magnetic field is also used to define the direction of diffusion which is being examined. By means of an appropriate pulsing of a radio-frequency field in a spin-echo measurement, the relative ease/difficulty of the diffusion of “active” molecules, principally water, is measured. The strength of the signal obtained at a voxel located at x ∈ X in the direction $s\in {\mathbb{S}}^2$ is given by the following form of the so-called Stejskal-Tanner equation [8],

$$S(x,s)=S_0\left(x\right)\exp (-b{\widehat{g}}_s^{T}D\left(x\right){\widehat{g}}_s),$$

where

• S₀(x) is the strength of the signal obtained in the absence of diffusion,

• b > 0 is the so-called b-value which is determined by the physics behind the measurement, i.e., the gyromagnetic ratio, γ, along with the pulsing time T₁,

• ${\widehat{g}}_s$ is the unit vector with direction $s\in {\mathbb{S}}^2,$ the direction of the gradient of the linear diffusion gradient field,

• D(x) is the 3 × 3 diffusion tensor at x. If the medium is isotropic in a neighbourhood of x, then $D=dI,$ where d is a constant (the diffusion coefficient) and I is the identity matrix.

In HARDI imaging, the signal attenuation can also be modelled as follows [9],

$$S(x,s)=S_0\left(x\right)\exp ((-bd(x,s\left)\right),$$

where d(x, s), the spherical Apparent Diffusion Coefficient(sADC), characterizes the rate of diffusion in the direction $s\in {\mathbb{S}}^2$. The sADC d(x, s) may be viewed as a function defined on the unit sphere ${\mathbb{S}}^2,$ which provides a local map of the rates of diffusion of water molecules from location x in all (or at least experimentally realizable) directions. A related quantity is the diffusion orientation probability distribution(dODF) O(x, s), which is defined to be the probability that a diffusing water molecule at x moves in direction s [8]. (As is well known – and we simply mention it here to keep our discussion brief – the method of tractography [2], [12] uses this information to determine the anatomical structure of white matter in the brain by means of connectomes, visual representations of neural fibre connections. This information can then be used for diagnostic purposes [4].)

In this study, we choose to work with the raw experimental data S(x, s_k), 1 ≤ k ≤ K, which is contaminated – usually quite highly – with instrument noise. Our goal is to denoise this data. From Eqs. (1.2) and (1.3), we see that diffusion attenuates the (theoretically noiseless) signal S. For a given x ∈ X, a direction for which d(x, s) (or O(x, s)) is large (small), i.e., a large diffusion rate, corresponds to a small (large) value of the signal strength S(x, s). This implies that the signal strength S(x, s) is (exponentially) inversely related to the diffusion rate at x in the direction s – essentially a measure of the resistance to diffusion. Even though, for each x ∈ X, the signal S(x, s), considered as a function of $s\in {\mathbb{S}}^2$ does not define a measure over ${\mathbb{S}}^2$ (in contrast to the ODF defined earlier), we propose to treat is as a probability measure (after suitable normalization). The motivation is that measure-based distance functions, such as the Monge-Kantorovich metric, should work better than simple metrics for functions such as L^p to characterize differences in information – in this case anisotropic diffusion – that is contained in the dMRI data. Indeed, there has been interest in the use of metrics derived from information theory (e.g., Fisher-Rao metric, von Mises-Fisher distribution) to measure the overlap/difference between ODFs – see the discussion and references contained in [5], [9] – but, to the best of our knowledge, no use of such metrics in any denoising algorithm has yet been performed, apart from [11] and this paper.