Joint phase reconstruction and magnitude segmentation from velocity-encoded MRI data

TL;DR

This paper introduces a joint optimization method using non-convex Bregman iteration to simultaneously estimate velocity, magnitude, and segmentation in velocity-encoded MRI, enhancing accuracy over traditional sequential methods.

Contribution

The work presents a novel joint estimation algorithm for velocity, magnitude, and segmentation in velocity-encoded MRI data, improving upon classical approaches.

Findings

Joint model outperforms sequential methods in accuracy

Numerical experiments validate the approach on synthetic and real data

Method effectively estimates dynamic flow components in bubbly flow imaging

Abstract

Velocity-encoded MRI is an imaging technique used in different areas to assess flow motion. Some applications include medical imaging such as cardiovascular blood flow studies, and industrial settings in the areas of rheology, pipe flows, and reactor hydrodynamics, where the goal is to characterise dynamic components of some quantity of interest. The problem of estimating velocities from such measurements is a nonlinear dynamic inverse problem. To retrieve time-dependent velocity information, careful mathematical modelling and appropriate regularisation is required. In this work, we propose an optimisation algorithm based on non-convex Bregman iteration to jointly estimate velocity-, magnitude- and segmentation-information for the application of bubbly flow imaging. Furthermore, we demonstrate through numerical experiments on synthetic and real data that the joint model improves velocity, magnitude and segmentation over a classical sequential approach.

Tables

Table 1. Table 1: MSE for phase ( ϕ 1 subscript italic-ϕ 1 \phi_{1} and ϕ 2 subscript italic-ϕ 2 \phi_{2} ) and magnitude ( u 1 subscript 𝑢 1 u_{1} and u 2 subscript 𝑢 2 u_{2} ) images for the sequential and joint approaches. The error is significantly decreased using our proposed joint approach.

	$u_1$	$u_2$	${\varphi }_1$	${\varphi }_2$
\svhline Sequential	0.0019	0.0028	0.0032	0.0059
Joint	0.0011	0.0012	0.0018	0.0051

Equations

$$\displaystyle\frac{d}{dt}\left(\begin{array}[]{c}M_{x}(t)\\ M_{y}(t)\\ M_{z}(t)\end{array}\right)=\left(\begin{array}[]{ccc}-\frac{1}{T_{2}}&\gamma B_{z}(t)&-\gamma B_{y}(t)\\ -\gamma B_{z}(t)&-\frac{1}{T_{2}}&\gamma B_{x}(t)\\ \gamma B_{y}(t)&-\gamma B_{x}(t)&-\frac{1}{T_{1}}\end{array}\right)\left(\begin{array}[]{c}M_{x}(t)\\ M_{y}(t)\\ M_{z}(t)\end{array}\right)+\left(\begin{array}[]{c}0\\ 0\\ \frac{M_{0}}{T_{1}}\end{array}\right)\,\text{.}$$

$$\frac{d}{dt}{M}_{xy}(t)$$

$$\frac{d}{dt}{M}_z(t)$$

$$\frac{d}{dt}{M}_{xy}(t)$$

$$\frac{d}{dt}{M}_z(t)$$

$$M_{xy}(t)$$

$$M_z(t)$$

$$\frac{d}{dt}{M}_{xy}(t)$$

$$M_{xy}(t)$$

$$x(t)=\sum _{n=0}^{\infty }\frac{x^{(n)}(\Delta t)}{n!}{t}^n,$$

$${\int }_{\Delta t}^{t}x(\tau )\cdot g(\tau )d\tau =\sum _{n=0}^{\infty }\left[\frac{x^{(n)}(\Delta t)}{n!}\cdot {\int }_{\Delta t}^{t}g(\tau ){\tau }^{n}d\tau \right]\text{.}$$

$${\int }_{\Delta t}^{t}x(\tau )\cdot g(\tau )d\tau =x\cdot {\int }_{\Delta t}^{t}g(\tau )d\tau +\phi \cdot {\int }_{\Delta t}^{t}g(\tau )\tau d\tau \text{,}$$

$${\int }_{\Delta t}^{t}g(\tau )d\tau =0\text{and}{\int }_{\Delta t}^{t}g(\tau )\tau d\tau =\zeta (t),$$

$${\int }_{t_2}^{t}g(\tau )d\tau =\xi (t)\text{and}{\int }_{t_2}^{t}g(\tau )\tau d\tau =0,$$

$$f(t)={\int }_{{\mathbb{R}}^3}u(x)e^{-i\gamma \left(B_0(x)t+\phi (x)\cdot \zeta (t)\right)}{e}^{-i\gamma x\cdot \xi (t)}dx.$$

$$f_{+}(t)$$

$$f_{-}(t)$$

$$\phi (x)\cdot \zeta (t)=\frac{1}{2}\left({\phi }_{+}(x,t)-{\phi }_{-}(x,t)\right).$$

$$f_{{\text{noflow}}_{+}}(t)$$

$$f_{{\text{noflow}}_{-}}(t)$$

$$\phi (x)\cdot \zeta (t)=\frac{1}{2}\left(\left({\phi }_{+}(x,t)-{\phi }_{-}(x,t)\right)-\left({\phi }_{{\text{noflow}}_{+}}(x,t)-{\phi }_{{\text{noflow}}_{-}}(x,t)\right)\right).$$

$$f_j={\int }_0^{t_{\ast }}\Psi (t,t_j)\left[{\int }_{{\mathbb{R}}^3}u(x)e^{i\varphi (x,t)}{e}^{-i\gamma x\cdot \xi (t)}dx\right]dt,$$

$$f_j={\int }_{{\mathbb{R}}^3}u(x)e^{i\varphi (x)}{e}^{-i\gamma x\cdot {\xi }_j}dx,j\in \{1,\dots ,m\}.$$

$$f=SF\left(u{e}^{i\varphi }\right),$$

$$f=SF\left(u{e}^{i\varphi }\right)+\eta ,$$

$$r_z=A^{\ast }f$$

$$r_{j}\in\underset{r}{\operatorname{arg}\,\operatorname{min}}\;\Big{\{}\frac{1}{2}\|A_{j}r-f_{j}\|_{2}^{2}+\alpha J(r)\Big{\}},$$

$${\varphi }_j=\arg (r_j).$$

$$E(u,v,\varphi )=$$

$$\displaystyle+\delta\underbrace{\sum_{n}v_{nj}(c_{1}-u_{nj})^{2}+(1-v_{nj})(c_{2}-u_{nj})^{2}\Bigg{\}}}_{\text{segmentation}}\,.$$

$$u_l^{k+1}$$

$$p_l^{k+1}$$

$$v_l^{k+1}$$

$$q_l^{k+1}$$

Full text

1 Introduction

Magnetic resonance imaging (MRI) is an imaging technique that allows to visualise the chemical composition of patients or materials in a non-invasive fashion. Besides resolving in great detail the morphology of the object under consideration, MRI is intrinsically sensitive to motion, flow and diffusion Burt19821044 ; Axel1984 . This means that in a single experiment, MRI can produce both structural and functional information. By designing the acquisition protocol appropriately, MRI can provide flow and motion estimation. This technique is known as MR velocimetry or phase-encoded MR velocity imaging Callaghan ; doi:10.1146/annurev.fluid.31.1.95 ; Elkins2007 ; GLADDEN20132 . In this work, we will focus on the dynamic inverse problem involved in recovering velocities from this kind of data.

In many MRI applications, the goal is not only to extract the structure of the object of interest, but also to estimate some functional features. An example is flow imaging, in which the aim is to reconstruct the velocity of the fluid that is moving in some structure. In order to acquire the velocity information and assess flow motion, phase-encoded MR velocity imaging is widely used in different areas. In medical imaging, this is used for example in cardiovascular blood flow studies to assess the distribution and variation in flow in blood vessels around the heart Gatehouse2005 . Other industrial applications include the study the rheology of complex fluids Callaghan1999 , liquids and gases flowing through packed beds Sederman1998 ; Sankey2009 ; Holland2010 , granular flows Fukushima99 ; Holland2008 and multiphase turbulence Tayler2012 .

MRI scanners use strong magnetic fields and radio waves to excite subatomic particles (like protons) that subsequently emit radio frequency signals which can be measured by the radio frequency coils. Because the local magnetisation of the spins is a vector quantity, it is possible to derive both magnitude and phase images from the signal. Furthermore, for appropriately designed experiments, the velocity information can be estimated from the phase image. The problem of retrieving magnitude and phase (and therefore velocities) from such measurements is non-linear. Many standard approaches reduce this inverse problem to a complex but linear inverse problem, where magnitude and phase are estimated subsequently. With this strategy, however, it is impossible to impose regularity on the velocity information. In this work, we therefore propose a joint framework to simultaneously estimate magnitude and phase from undersampled velocity-encoded MRI. Based on Corona , we additionally introduce a third task, that is the segmentation on the magnitude, to improve the overall reconstruction quality. The main motivation is that by estimating edges simultaneously from the data, both magnitude and segmentation are reconstructed more accurately. By enhancing the magnitude reconstruction, we expect in turn to improve the corresponding phase image and therefore the final velocity estimation.

Contributions

In this work we consider the problem of estimating flow, magnitude and segmentation of regions of interest from undersampled velocity-encoded MRI data. The problem is of great interest in different areas including cardiovascular blood flow analysis in medical imaging and rheology of complex fluids in industrial applications. To this end, we propose a joint variational model for undersampled velocity-encoded MRI. The significance of our approach is that by tackling the phase and magnitude reconstruction jointly, we can exploit their strong correlation and finally impose regularity on the velocity component. This is further assisted by the introduction of a segmentation term as additional prior to enhance edges of the regions of interest. Our main contributions are

• •

A description of the forward and inverse problem of velocity-encoded MRI in the setting of bubbly flow estimation.

• •

A joint variational framework for the approximation of the non-linear inverse problem of velocity estimation. We show that by exploiting the strong correlation in the data, our joint method yields an accurate estimation of the underlying flow, alongside a magnitude reconstruction that preserves and enhances intrinsic structures and edges, due to a joint segmentation approach. Moreover, we achieve an accurate segmentation to discern between different areas of interest, e.g. fluid and air.

• •

An alternating Bregman iteration method for non-convex optimisation problems.

• •

Numerical experiments on synthetic and real data in which we demonstrate the suitability and potential of our approach and provide a comparison with sequential approach.

Organisation of the paper

This paper is organised as follows. In Section 2 we describe the derivation of the inverse problem of velocity-encoded MRI from the acquisition process to the spin proton density estimation. In Section 3 we present our joint variational model to jointly estimate phase and magnitude reconstruction and its segmentation. In Section 4 we propose an optimisation scheme to solve the non-convex and non-linear problem using Bregman iteration. To conclude, in Section 5 we demonstrate the performance of our proposed joint method in comparison with a sequential approach for synthetic and real MRI data.

2 Velocity-encoded MRI

In the following we will briefly describe the mathematics of the acquisition process involved in MRI velocimetry. Subsequently we are going to see that finding the unknown spin proton density basically leads to solving the inverse problem of the Fourier transform.

2.1 From the Bloch equations to the inverse problem

The magnetisation of a so-called spin isochromat can be described by the Bloch equations

[TABLE]

Here $M(t)=(M_{x}(t),M_{y}(t),M_{z}(t))$ is the nuclear magnetisation (of the spin isochromat), $\gamma$ is the gyromagnetic ratio, $B(t)=(B_{x}(t),B_{y}(t),B_{z}(t))$ denotes the magnetic field experienced by the nuclei, $T_{1}$ is the longitudinal and $T_{2}$ the transverse relaxation time and $M_{0}$ the magnetisation in thermal equilibrium. If we define $M_{xy}(t)=M_{x}(t)+iM_{y}(t)$ and $B_{xy}(t)=B_{x}(t)+iB_{y}(t)$, we can rewrite (13) to

[TABLE]

with $\overline{\cdot}$ denoting the complex conjugate of $\cdot$.

If we assume for instance that $B=(0,0,B_{0})$ is just a constant magnetic field in $z$-direction, (14) reduces to the decoupled equations

[TABLE]

It is easy to see that this system of equations (15) has the unique solution

[TABLE]

for $\omega_{0}:=\gamma B_{0}$ denoting the Lamor frequency, and $M_{xy}(\Delta t)$, $M_{z}(\Delta t)$ being the initial magnetisations at time $t=\Delta t$.

2.2 Signal recovery

The key idea to enable spatially resolved nuclear magnetic resonance spectrometry is to add a magnetic field $\hat{B}(t)$ to the constant magnetic field $B_{0}$ in $z$-direction that varies spatially over time. Then, (15a) changes to

[TABLE]

if we ensure $\hat{B}(\Delta t)=0$. If now $x(t)$ denotes the spatial location of a considered spin isochromat at time $t$, we can write $\hat{B}(t)$ as $\hat{B}(t)=x(t)\cdot g(t)$, with a function $g$ that describes the influence of the magnetic field gradient over time.

If a radio-frequency (RF) pulse that has been used to induce magnetisation in the $x$-$y$-plane is subsequently turned off at time $t=t_{\ast}$ and thus, $B_{x}(t)=0$ and $B_{y}(t)=0$ for $t>t_{\ast}>\Delta t$, the same coils that have been used to induce the RF pulse can be used to measure the $x$-$y$ magnetisation. Using (16a) and assuming $t_{\ast}<t\ll T_{2}$ for all $x\in\mathbb{R}^{3}$, this gives rise to the following model-equation:

[TABLE]

In the following we assume that $x(t)$ can be approximated reasonably well via its Taylor approximation around $t=\Delta t$, i.e.

[TABLE]

which yields

[TABLE]

It is well-known that appropriate application of gradients (i.e. appropriate design of $g$) enables the approximation of individual moments of (19). If we further assume that the system to be observed does only contain zero- and first-order moments, we can assume

[TABLE]

where $x$ is now short for $x({\Delta t})$ and $\varphi:=x^{\prime}({\Delta t})$ is the corresponding velocity information.

In order to turn (18) into a useful mathematical model we need to encode velocity information and remove the temporal dependency of $x$. In order to do so, the gradient coils need to be programmed to first enable the encoding of velocity information via a function $g$ that satisfies

[TABLE]

for time $\Delta t\leq t\leq t_{1}$ and a function $\zeta:\mathbb{R}_{\geq 0}\mapsto\mathbb{R}^{3}$. Subsequently, the gradients need to be programmed to encode spatial information by choosing $g$ with

[TABLE]

for $t_{2}>t_{1}$ and a function $\xi:\mathbb{R}_{\geq 0}\mapsto\mathbb{R}^{3}$. Since the RF-coils measure a volume of the whole $x$-$y$ net-magnetisation, the acquired signal then equals

[TABLE]

with $u(x)$ denoting the spin-proton density $M_{xy}({\Delta t})$ at a specific spatial coordinate $x\in\mathbb{R}^{3}$. Note that for $r(x):=u(x)\,e^{-i\gamma\left(B_{0}(x)t+\varphi(x)\cdot\zeta(t)\right)}$ we observe that $f$ is just the Fourier transform of the complex signal $r$ with magnitude $u$ and phase $-\gamma(B_{0}t+\varphi\cdot\zeta)$.

2.3 Removal of background magnetic field

Our goal is to recover the velocity information $\varphi$ from $f$. Assuming that we do not know $B_{0}$, we can alternatively conduct two experiments, where the setup is identical apart from the velocity-encoding gradients having opposite polarities, i.e. we measure

[TABLE]

Hence, if we denote $\varphi_{+}(x,t):=B_{0}(x)t+\varphi(x)\cdot\zeta(t)$ and $\varphi_{-}(x,t):=B_{0}(x)t-\varphi(x)\cdot\zeta(t)$, we immediately observe

[TABLE]

The inverse problem of (22) is to recover $u$ and $\varphi$ from $f_{+}$ and $f_{-}$.

2.4 Zero-flow experiment

A zero-flow experiment that allows for the removal of additional artefacts is also conducted. This experiment is to account for imperfections in the measurement system which cause an added signal between the positive and negative $\zeta$ experiments even in the absence of flow, and enables a correction that allows direct quantification of flow and tissue motion. We refer to this technique as flow compensation, which consists of acquiring a reference scan, with any flow switched-off, with vanishing zero and first gradient moments, before the actual velocity encoding scan with added bipolar gradients is performed. In this way, we obtain background phase images from the reference scan, and velocity sensitivity with the second flow-sensitive scan. In practice, this means that in addition to (22), the following two measurements are taken:

[TABLE]

so that the actual velocity information can be recovered via

[TABLE]

The inverse problem is to recover $u$ and $\varphi$ from (22) and (23) via (24). More details on phase-encoded MR velocity imaging can be found in Markl20061VE .

In other words, for a given direction of the velocity to be measured ($x$, $y$ or $z$), the corresponding component velocity map ($v_{x}$, $v_{y}$ or $v_{z}$) is acquired by applying repeatedly a pulse sequence with the velocity-encoding gradient in the respective direction ($x$, $y$ or $z$) and with alternating polarity between consecutive pulse sequences (from $\pm g$ to $\mp g$). The difference between the phase of the MRI image reconstructed from the acquired k-space data of consecutive pulse sequences, and the reference to a zero flow experiment, yields the component velocity map.

2.5 Sampling

The MRI signal is acquired by sampling the continuous signals of $f_{+}$, $f_{-}$, $f_{\text{noflow}_{+}}$ and $f_{\text{noflow}_{-}}$ at $m$ discrete points in time. Hence, for each phase $-\phi$ the data acquisition reads as

[TABLE]

for $j\in\{1,\ldots,m\}$ and where $\Psi$ denotes the sampling function or distribution. If we for example assume $\Psi(t,t_{j})=\delta(t-t_{j})$, where $\delta$ denotes the Dirac delta distribution, and that $\phi(x,t)$ is constant w.r.t. time, then (25) simplifies to

[TABLE]

We want to denote the sub-sampled Fourier transform with $SF$, and therefore rewrite (26) to

[TABLE]

where $f\in\mathbb{R}^{m}$ denotes the vector of k-space samples. Sampling strategies are very important to reduce the acquisition times and therefore to be able to image dynamic systems using velocity-encoded MRI through fast imaging techniques. The main idea is to exploit redundancy in some specific domain of the measured data. This approach is strongly related to the theory of compressed sensing (CS) Candes2004 ; Donoho2006 ; Lustig2007 and many image reconstruction techniques have been proposed Holland2010 ; Tayler2012 ; Jung2008 ; Paciok2011 ; Paulsen2010 ; Parasoglou2009 .

Depending on whether $\gamma\xi$ is sampled on a uniform or non-uniform grid, $SF$ can be realised via the Fast Fourier Transform (FFT) fft or via a non-uniform Fourier Transform such as NUFFT NUFFT .

2.6 Dynamic inverse problem

We want to highlight that every $u$ and $\phi$ in (27) implicitly depends on an initial time $\Delta t$, which becomes evident from the derivation in Section 2.2. Hence, if we take measurements for a sequence $\{\Delta t_{j}\}_{j=1}^{s}$ with $0=\Delta t_{1}<\Delta t_{2}<\ldots<\Delta t_{s}$, we are introducing a discrete temporal dimension to our inverse problem that potentially allows us to exploit any temporal correlation between frames $\{u_{j}\}_{j=1}^{s}$ and $\{\phi_{j}\}_{j=1}^{s}$. However, we will only consider the reconstruction of individual frames throughout this work for reasons that we are going to address later.

In the following we will refer to an individual frame of the dynamic inverse problem for velocity-encoded MRI in the discrete setting and under the presence of noise making use of the notation of the discrete Fourier transform operator.

3 Mathematical modelling

In this section we first present the velocity-encoded MRI reconstruction inverse problem in the presence of noise and discuss a sequential variational regularisation scheme to approximate the solution. Secondly, we introduce our joint reconstruction and segmentation approach in a Bregman iteration framework to jointly estimate phase, magnitude and segmentation.

3.1 Indirect phase-encoded MR velocity imaging

The velocity-encoded MRI image reconstruction problem is described as follows. Let $u,\phi\in\mathbb{R}^{n}$ be the proton density or magnitude image and correspondent phase image, respectively, in a discretised image domain $\Omega:=\{1,\dots,n_{1}\}\times\{1,\dots,n_{2}\}$, with $n=n_{1}n_{2}$. The vector $f=(f_{i})_{i=1}^{m}\in\mathbb{C}^{m}$ with $m\ll n$ are the measured Fourier coefficients obtained from (27). Based on (27) the forward model for noisy data is given by

[TABLE]

where $i^{2}=-1$ and $\eta$ is Gaussian noise with zero mean and standard deviation $\sigma$. For brevity we will follow the notation $A=SF$. As explained in the previous section, velocity information is encoded in the phase image. However, during the acquisition the phase is perturbed by an error due to field inhomogeneity and chemical shift. To account for this error, usually different measurements corresponding to different polarities of encoding flow gradients are acquired. Then the velocity (in one direction) at one particular time will be estimated as in (24), where $\zeta$ is a constant known from the acquisition setting.

Given the presence of noise and partial observation of the data due to undersampling, the problem described in (28) is ill-posed. A simple strategy to obtain an approximated solution is to replace with zero the missing Fourier coefficients and compute the so-called zero-filling solution

[TABLE]

where $r=ue^{i\phi}$. However, these reconstructed images will present aliasing artefacts because of the undersampling. A classical approach to solve this problem is to compute approximate solutions of (28) using a variational regularisation approach. We consider a Tikhonov-type regularisation approach that reads

[TABLE]

for $j\in\{1,\dots,4\}$ being the different measurements, where the first term is the data fidelity that imposes consistency between the reconstruction and the given measurements $f$, the second term is the regularisation, which incorporates some prior knowledge of the solution. The parameter $\alpha>0$ is a regularisation parameter that balances the two terms in the variational scheme. In this setting, the survey proposed in BenningGladdenHollandEtAl2014 describes different choices for the regularisation functional $J$, including wavelets and higher-order total variation (TV) schemes. Subsequently, the phases can be extracted from these complex images $r_{j}=u_{j}e^{i\phi_{j}}$ as

[TABLE]

More recently, other reconstruction approaches have been proposed to regularise the phase of the image Fessler2004 ; Zibetti2010 ; Zhao2012 ; valkonen2014primal ; Zibetti:2017 . All these methods rely on modelling separately prior knowledge on the magnitude and on phase images and differ on the optimisation schemes involved in the non-convex and non-linear problem. However, while it is possible to exploit information about the velocity from fluid mechanics, it is in general hard to assume specific knowledge on the individual phases. As explained in the previous section and described in (24), velocities are computed as phase differences of different MR measurements and therefore the regularisation needs to be imposed on the phase difference rather than individual phases. In this work, we step away from the approach of only regularising individual phases and propose instead to regularise the velocity as difference of phases. In the following we describe our choice of regularisation and algorithmic framework for velocity-encoded MRI.

3.2 Joint variational model

In many industrial applications, velocity-encoded MRI is used to estimate flow of different chemical species in different physical status, such as gas-liquid systems Gladden2017 . In this case, one aims at recovering a piecewise constant image or an image with sharp edges to facilitate further analysis such as identification of regions of interest. It was proposed in Corona to use a segmentation task as additional regularisation on the reconstruction to impose regularity in terms of sharp edges. It was shown there that this is highly beneficial for very low undersampling rates in MRI. In this work, we expand this idea to the phase-encoded MR velocity imaging data, where the idea is to jointly solve for magnitude, segmentation and phase improving performances on the three tasks.

Following the work in Corona , we are interested in the joint model to recover magnitude $u_{j}$ and velocity $\varphi$ components through the measured phases $\phi_{j}$ from undersampled MRI data $f_{i}$ and to estimate a segmentation $v_{j}$ on the magnitude images. As described in the previous section, we are dealing with four MRI measurements to obtain one component velocity image. Defining the shorthand notations $u:=\{u_{j}\}_{j=1}^{4}$, $v:=\{v_{j}\}_{j=1}^{4}$ and $\phi:=\{\phi_{j}\}_{j=1}^{4}$, this joint model reads as

[TABLE]

The first term in (32) describes the reconstruction fidelity term for the magnitudes $u$ and phases $\phi$ for the given data $f:=\{f_{j}\}_{j=1}^{4}$. The second term represents the segmentation problem to find partitions $v$ of the images $u$ in two disjoint regions that have mean intensity values close to the constants $c_{1}$ and $c_{2}$ Chant2001 ; Chan2006 . The parameter $\delta$ weighs the effect of the segmentation onto the reconstruction. The underlying idea is to exploit structure and redundancy in the data, estimating edges simultaneously from the data, ultimately improving the reconstruction. By incorporating prior knowledge of the regions of interest we impose additional regularity of the solution.

The joint cost function (32) is non-convex. While sub-problems in $u$ and $v$ (leaving the other parameters fixed) are convex, the sub-problems in $\phi$ are non-linear and non-convex. In the next section we present a unified framework based on non-convex Bregman iterations to solve the joint model.

4 Optimisation

There are many ways of minimising (32). We want to pursue a strategy that guarantees smooth velocity-components, piecewise-constant segmentations and magnitude images with sharp transitions in an inverse scale-space fashion. In order to achieve those features, we aim to approximate minimisers of (32) via an alternating Bregman proximal method or Bregman iteration of the form

[TABLE]

for $l=1,\ldots,d:=4$, $u:=(u_{l})_{l=1}^{d}$, $v:=(v_{l})_{l=1}^{d}$ and $\phi:=(\phi_{l})_{l=1}^{d}$. Here $J_{u}$, $J_{v}$ and $J_{\phi}$ are proper, lower semi-continuous and convex functions and $D_{J_{u}}^{p_{l}^{k}}(u,u_{l}^{k})$, $D_{J_{v}}^{q_{l}^{k}}(v,v_{l}^{k})$ and $D_{J_{\phi}}^{w^{k}}(\phi,\phi^{k})$ are the corresponding generalised Bregman distances bregman1967relaxation ; kiwiel1997proximal with arguments and corresponding subgradients $p_{l}^{k}$, $q_{l}^{k}$ and $w^{k}$. A generalised Bregman distance is the distance between a function $J$ evaluated at argument $u$ and its linearisation around argument $v$, i.e.

[TABLE]

for a subgradient $q\in\partial J(v)$. Note that algorithm (33) has update rules for the subgradients, as $J_{u}$, $J_{v}$ and $J_{\phi}$ are allowed to be non-smooth, which makes the selection of particular subgradients necessary.

The algorithm is a hybrid of the algorithms proposed in Benning2017 and Corona . For both algorithms global convergence results, motivated by xu2013block ; bolte2014proximal , have been established. Since we deal with imperfect data potentially corrupted by measurement noise and numerical errors, we will, however, use (33) in combination with an early-stopping criterion in order not to converge to a minimiser of (32) but to approximate the solution of (28) via iterative regularisation.

The crucial part for the application of (33) are the choices of the underlying functions $J_{u}$, $J_{v}$ and $J_{\phi}$ of the corresponding Bregman distances. We want both the magnitude images and the segmentations to maintain sharp discontinuities and therefore want to penalise their discretised, isotropic, total variation. On the other hand, we want to guarantee smooth components of our velocity field, which is why we penalise them with the two-norm of a discretised gradient. In particular, we choose

[TABLE]

to be the isotropic total variation with weights $\alpha>0$ and $\beta>0$, where $\nabla:\mathbb{R}^{n}\rightarrow\mathbb{R}^{2n}$ denotes a forward finite-difference approximation of the gradient, $|\cdot|$ the Euclidean vector norm and $\|\cdot\|_{1}$ the pixel-wise one-norm. Further, we choose $J_{\phi}$ in a way that allows to enable an $H^{1}$-norm-type smoothing on the difference of the phases, i.e.

[TABLE]

where $\eta>0$ denotes another weight. Note that all convex sub-optimisation-problems in (33) are solved numerically with a primal-dual hybrid gradient (PDHG) method zhu2008efficient ; pock2009algorithm ; esser2010general ; chambolle2011first . Once we have approximated the magnitudes, labels and phases with this iterative regularisation strategy, we can compute the velocity components via (24).

5 Numerical results

In this section we present numerical results of our method for the specific application of bubble burst hydrodynamics using MR velocimetry. The hydrodynamics of bursting bubbles is important in many different areas such as geophysical processes and bioreactor design. We refer to Reci for an overview on the field and the description of results on the first experimental measurement of the liquid velocity field map during the burst of a bubble at the liquid surface interface.

5.1 Case-study on simulated dataset

To quantitatively evaluate our method, we consider the simulated k-space data of a rising spherical bubble in an infinite fluid during Stoke’s flow regime. The simulated data consists of 32 time frames, but for the sake of compactness we will show some visual outputs for one time step $t=19$.

We assess the performance of our approach for velocity and magnitude estimation by comparing our solutions with respect to the groundtruth and using the mean squared error (MSE) defined as $\|x^{\text{groundtruth}}-x\|_{2}^{2}/n$, where $n$ is the number of pixels in the image.

We also present a comparison with a sequential approach, where the magnitude is obtained with a classic CS TV-regularised approach and the phase is subsequently estimated using the method proposed in Benning2017 and presented in Reci for the evaluation of bubbly flow estimation.

In Fig. 1 we can see the results for the sequential approach compared to the joint approach when sampling only 11% of the k-space data. Although visually there is not significant change, the MSE shows a big improvement for the joint approach. This confirms that using our joint model is relevant for the problem of velocity-encoded MRI. For the 32 frames, we report the average MSE for magnitude and phase in Table 5.1 where can see a drastic improvement compared to the sequential approach.