Enhancing joint reconstruction and segmentation with non-convex Bregman iteration

All imaging modalities such as computed tomography (CT), emission tomography and magnetic resonance imaging (MRI) require a reconstruction approach to produce an image. A common image processing task for applications that utilise those modalities is image segmentation, typically performed posterior to the reconstruction. We explore a new approach that combines reconstruction and segmentation in a unified framework. We derive a variational model that consists of a total variation regularised reconstruction from undersampled measurements and a Chan-Vese based segmentation. We extend the variational regularisation scheme to a Bregman iteration framework to improve the reconstruction and therefore the segmentation. We develop a novel alternating minimisation scheme that solves the non-convex optimisation problem with provable convergence guarantees. Our results for synthetic and real data show that both reconstruction and segmentation are improved compared to the classical sequential approach.


Introduction
Image reconstruction plays a central role in many imaging modalities for medical and non-medical applications. The majority of imaging techniques deal with incomplete data and noise, making the inverse problem of reconstruction severely ill-posed. Based on compressed sensing (CS) it is possible to tackle this problem by exploiting prior knowledge of the signal [1][2][3]. Nevertheless, reconstructions from very noisy and undersampled data will present some errors that will be propagated into further analysis, e.g. image segmentation. Segmentation is an image processing task used to partition the image into meaningful regions. Its goal is to identify objects of interest, based on contours or similarities in the interior. Typically segmentation is performed after reconstruction, hence its result strongly depends on the quality of the reconstruction. Recently the idea of combining reconstruction and segmentation has become more popular. The main motivation is to avoid error propagations that occur in the sequential approach by estimating edges simultaneously from the data, ultimately improving the reconstruction. In this paper, we propose a new model for joint reconstruction and segmentation from undersampled MRI data. The underlying idea is to incorporate prior knowledge about the objects that we want to segment in the reconstruction step, thus introducing additional regularity in our solution. In this unified framework, we expect that the segmentation will also benefit from sharper reconstructions. We demonstrate that our joint approach improves the reconstruction quality and yields better segmentations compared to sequential approaches. In Figure 1, we consider a brain phantom from which we simulated the undersampled k -space data and added Gaussian noise. . Figure 1b and 1e present reconstructions and segmentations obtained with the sequential approaches, while Figure 1c and 1f show the results for our joint approach. The reconstruction using our method shows clearly more details and it is able to detect finer structures that are not recovered with the classical separate approach. As a consequence, the joint segmentation is also improved. In the following section we present the mathematical models that we used in our comparison. We investigated the performance of our model for two different applications: bubbly flow and cancer imaging. We show that both reconstruction and segmentation benefit from this method, compared to the traditional sequential approaches, suggesting that error propagation is reduced.
Our contribution. In our proposed joint method, we obtain an image reconstruction that preserves its intrinsic structures and edges, possibly enhancing them, thanks to the joint segmentation, and simultaneously we achieve an accurate segmentation. We consider the edge-preserving total variation regularisation for both the reconstruction and segmentation term using Bregman distances. In this unified Bregman iteration framework, we have the advantage of improving the reconstruction by reducing the contrast bias in the TV formulation, which leads to more accurate segmentation. In addition, the segmentation constitutes another prior for the reconstruction by enhancing edges of the regions of interest. Furthermore, we propose a non-convex alternating direction algorithm in a Bregman iteration scheme for which we prove global convergence.
The paper is organised as follows. In Section 2 we describe the problems of MRI reconstruction and region-based segmentation. We then introduce our joint reconstruction and segmentation approach in a Bregman iteration framework. This section also contains a detailed comparison of other joint models in the literature. In Section 3 we study the non-convex optimisation problem and present the convergence analysis for  this class of problems. Finally in Section 4 we present numerical results for MRI data for different applications. Here we investigate the robustness of our model by testing the undersampling rate up to its limit and by considering different noise levels.

MRI reconstruction and segmentation
In the following section we introduce the mathematical tools to perform image reconstruction and image segmentation. In this work, we focus on the specific MRI application; however, our proposed joint method can be applied to other imaging problems in which the measured data is connected to the image via a linear and bounded forward operator, cf. Subsection 2.1. Finally we present our model that combines the two tasks of reconstruction and segmentation in a unified framework.

Reconstruction
In image reconstruction problems, we have the general setting where A : X → Y is a bounded and linear operator mapping between two vector-spaces.
The measured data f ∈ Y is usually corrupted by some noise η and often only observed partially. In this formulation we are interested in recovering the image u given the data f .
In this work, we focus on the application of MRI and we refer to the measurements f as the k -space data. In standard MRI acquisitions, the Fourier coefficients are collected in the k -space by radio-frequency (RF) coils. Because the k -space data is acquired sequentially, the scanning time is constrained by physical limitations of the imaging system. One of the most common ways to perform fast imaging consists of undersampling the k -space; this, however, only yields satisfactory results if the dimension of the parameter space can implicitly be reduced, for example by exploiting sparsity in certain domains. In the reconstruction, this assumption is incorporated in the regularisation term. Let Ω := {1, . . . , n 1 } × {1, . . . , n 2 } with n 1 , n 2 ∈ N be a discrete image domain. Let f = (f i ) m i=1 ∈ C m with m n = n 1 n 2 be our given undersampled k -space data, where f i ∈ C are the measured Fourier coefficients that fulfil the relationship (1) with A = SF . The operator A is now composed by S : C n → C m , which is a sampling operator that selects m measurements from the F u data according to the locations provided by a binary sampling matrix (see e.g. Figure 1d), where F is the discrete Fourier transform. In MRI, the noise η is drawn from a complex-valued Gaussian distribution with zero mean and standard deviation σ [4]. In problem (1) for MRI, the aim is to recover the image u ∈ C n from the data. However, in this work we follow the standard assumption that in many applications we have negligible phase, i.e. we are working with real valued, non-negative images. Therefore, we are only interested in u ∈ R n ; hence we consider the MRI forward operator as A : R n → C m and its adjoint A * : C m → R n as modelled in [5]. Problem (1) is ill-posed due to noise and incomplete measurements. The easiest approach to approximate (1) is to compute the solution, for which the missing entries are replaced with zero However, images reconstructed with this approach will suffer from aliasing artifacts because undersampling the k -space violates the Nyquist-Shannon sampling theorem. Therefore, we consider a mathematical model that incorporates prior knowledge by using a variational regularisation approach. A popular model is to find an approximate solution for u as a minimiser of the Tikhonov-type regularisation approach where the first term is the data fidelity that forces the reconstruction to be close to the measurements and the second term is the regularisation, which imposes some regularity on the solution. The parameter α > 0 is a regularisation parameter that balances the two terms in the variational scheme. In this setting, different regularisation functionals J can be chosen (see [6] for a survey of variational regularisation approaches). Although problems of the form (2) are very effective, they also lead to a systematic loss of contrast [7][8][9]. This is typically observed for common choices of the regulariser J, i.e., convex functional. To overcome this problem, [10] proposed an iterative regularisation method based on the generalised Bregman distance [11,12]. The Bregman distance with respect to J is defined as with p k ∈ ∂J(u k ), where ∂J(u k ) is called sub-differential and it is a generalisation of the classical differential for convex functions. We replace problem (2) with a sequence of minimisation problems The update on the subgradient can be conveniently computed by the optimality condition of (4) In this work, we will focus on one particular choice for J, namely the total variation. The total variation (TV) regularisation is a well-known edge-preserving approach, first introduced by Rudin, Osher and Fatemi in [13] for image denoising. The TV regularisation, i.e., the 1-norm penalty on a discrete finite difference approximation of the two-dimensional gradient ∇ : for the isotropic case. We then consider the Bregman iteration scheme in (4) for J(u) = TV(u). This approach is usually carried on by initialising the regularisation parameter α with a large value, producing overregularised initial solutions. At every step k, finer details are added. A suitable criterion to stop iterations (4) and (5) (see [6]), is the Morozov's discrepancy principle [14]. The discrepancy principle suggests to choose the smallest k ∈ N such that u k+1 satisfies where m is the number of samples and σ is the standard deviation of the noise in the data. Note that using Bregman iterations, the contrast is improved and in some cases even recovered exactly, compared to the variational regularisation model. In addition, it makes the regularisation parameter choice less challenging. Note that for different choices of J in (2), e.g., the Mumford-Shah/Potts model [15][16][17][18][19], we do not have loss of contrast, but we deal with a non-convex NP hard problem, algorithmically more challenging.

Segmentation
Image segmentation refers to the process of automatically dividing the image into meaningful regions. Mathematically, one is interested in finding a partition of the image domain Ω subject to ∪ l i=1 Ω i = Ω and ∩ l i=1 Ω i = ∅. One way to do this is to use region-based segmentation models, which identify regions based on similarities of their pixels. The segmentation model we are considering was originally proposed by Chan and Vese in [20] and it is a particular case of the piecewise-constant Mumford-Shah model [15]. Given an image function u : Ω → R, the goal is to divide the image domain Ω in two separated regions Ω 1 and Ω 2 = Ω \ Ω 1 by minimising the following energy function where C is the desired contour separating Ω 1 and Ω 2 , and the constants c 1 and c 2 represents the average intensity value of u inside C and outside C, respectively. The parameter β penalises the length of the contour C, controlling the scale of the objects in the segmentation. From this formulation we can make two observations: first, the regions Ω 1 and Ω \ Ω 1 can be represented by the characteristic function second, the perimeter of the contour identified by the the characteristic function corresponds to its total variation, as shown by the Coarea formula [21]. This leads to the new formulation .
Even assuming fixed constants c 1 , c 2 the problem is non-convex due to the binary constraint. In [22] the authors proposed to relax the constraint, allowing v(x) to assume values in the interval [0, 1]. They showed that for fixed constants c 1 , c 2 , global minimisers can be obtained by minimising the following energy followed by thresholding, setting Σ = {x : v(x) ≥ µ} for a.e. µ ∈ [0, 1]. As the problem is convex but not strictly convex, the global minimiser may not be unique. In practice we obtain solutions which are almost binary, hence the choice of µ is not crucial. Setting the energy (8) can be written in a more general form as 1] .
In this paper, we are interested in the extension of the two-class problem to the multiclass formulation [23]. Following the simplex-constrained vector function representation for multiple regions and its convex relaxation proposed in [24], we obtain as a special case a convex relaxation of the Chan-Vese model for arbitrary number of regions, which reads where C : to lie in the standard probability simplex. As in the binary case, the constants c i describe the average intensity value inside region i. In this case we consider the vector-valued formulation of TV TV(v) = Ω ∇v 1 2 + · · · + ∇v 2 dx.

Joint reconstruction and segmentation
MRI reconstructions from highly undersampled data are subject to errors, even when prior knowledge about the underlying object is incorporated in the mathematical model. It is often required to find a trade-off between filtering out the noise and retrieving the intrinsic structures while preserving the intensity configuration and small details. As a consequence, segmentations in the presence of artifacts are likely to fail. In this paper, we propose to solve the two image processing tasks of reconstruction and segmentation in a unified framework. The underlying idea is to inform the reconstruction with prior knowledge of the regions of interest, and simultaneously update this belief according to the actual measurements. Mathematically, given the undersampled and noisy k -space data f , we want to recover the image u : Ω → R and compute its segmentation v in disjoint regions, by solving the following problem where ı C (v) is the characteristic function over C : . . , n}}, and α, β, δ > 0 are some regularisation parameters. However, instead of solving (10), we consider the iterative regularisation procedure using Bregman distances. The main motivation is to exploit the contrast enhancement aspect for the reconstruction thanks to the Bregman iterative scheme. By improving the reconstruction, the segmentation is in turn refined. Therefore, we replace (10) with the following sequence of minimisation problems for k = 0, 1, 2, . . .
Note that (11) solves a problem different from (10). Assuming that a minimiser exists, the model (11) converges to a minimiser of as we will show in Subsection 3.1. In case of noisy data f this is not desirable, so that we combine the iteration with a stopping criterion in order to form a regularisation method. This model combines the reconstruction approach described in (4) and the discretised multi-class segmentation in (9) with a variation in the regularisation term, which is now embedded in the Bregman iteration scheme. In [25] the authors used Bregman distances for the Chan-Vese formulation (8), combined with spectral analysis, to produce multiscale segmentations.
As described in the previous subsection, the parameters α and β describe the scale of the details in u and the scale of the segmented objects in v. By integrating the two regularisations into the same Bregman iteration framework, we obtain that these scales are now determined by the iteration k + 1. At the first Bregman iteration k = 0, when α is very large, we obtain an over-smoothed u 1 , and the value of β is not very important. Intuitively, u 1 is almost piecewise constant with small total variation and a broad range of values of β may lead to very similar segmentations v 1 . However, at every iteration k + 1, finer scales are added to the solution with the update p k+1 . Accordingly, with the update q k+1 , which is independent of v k+1 , the segmentation keeps up with the scale in the reconstructed image u k+1 .
The novelty of this approach is also represented by the role of the parameter δ > 0. This parameter weighs the effect of the segmentation in the reconstruction, imposing regularity in u in terms of sharp edges in the regions of interest. In Section 4 we show how different ranges of δ affects the reconstruction (see Figure 12). Intuitively, large values of δ force the solution u to be close to the piecewise constant solution described by the constants c i . This is beneficial in applications where MRI is a means to extract shapes and sizes of underlying objects, (e.g. bubbly flow in Subsection 4.1). On the other hand, with very small δ, the segmentation has little impact and the solutions for u are close to the ones obtained by solving the individual problem (4). Instead, intermediate values of δ impose sharper boundaries in the reconstruction while preserving the texture.
Obviously, we need to stop the iteration before the residual brings back noise from the data f . As we cannot use Morozov discrepancy principle in this case (due to the fact that Au k − f 2 will rather increase due to the effect of the coupling term controlled by the parameter δ), we stop when two consecutive iterates in v are smaller than a certain tolerance, v k+1 − v k < tol, following the observation that the rate at which u k+1 changes close to the optimal solution is low, in contrary to more abrupt changes at the beginning of the Bregman iteration and later on when it starts to add noise.
Clearly, problem (11) is non-convex in the joint argument (u, v) due to the coupling term. However, it is convex in each individual variable. We propose to solve the joint problem by iteratively alternating the minimisation with respect to u and to v (see Section 3 for numerical optimisation and convergence analysis).

Comparison to other joint reconstruction and segmentation approaches
In this section we will provide an overview of some existing simultaneous reconstruction and segmentation (SRS) approaches with respect to different imaging applications.
CT/SPECT. Ramlau and Ring [26] first proposed a simultaneous reconstruction and segmentation model for CT, that was later extended to SPECT in [27] and to limited data tomography [28]. In these work, the authors aim to simultaneously reconstruct and segment the data acquired from SPECT and CT. CT measures the mass density distribution µ, that represents the attenuation of x-rays through the material; SPECT measures the activity distribution f as the concentration of the radio tracer injected in the material. Given the two measurements z δ and y δ , from CT and SPECT, they consider the following energy functional They propose a joint model based on a Mumford-Shah like functional, in which the reconstructions of µ and f and the given data are embedded in the data term in a least squares sense. The operators A and R are the attenuated Radon transform (SPECT operator) and the Radon transform (CT operator), respectively. The penalty term is considered to be a multiple of the lengths of the contours of µ, Γ µ and the contours of f , Γ f . These boundaries are modelled using level set functions. In these segmented partitions of the domain, µ and f are assumed to be piecewise constant. The optimisation problem is then solved alternatively with respect to the functional variables f and µ with fixed geometric variables Γ µ and Γ f and the other way around.
In [29] the simultaneous reconstruction and segmentation is applied to dynamic SPECT imaging, which solves a variational framework consisting of a Kullback-Leibler (KL) data fidelity and different regulariser terms to enforce sharp edges and sparsity for the segmentation and smoothness for the reconstruction. The cost function is Given the data g, they want to retreive the concentration curves c k (t) in time for K disjoint regions and their indication functions u k (x) in space. The optimisation is carried out alternating the minimisation over u having c fixed and then over c having u fixed.
In [30] they propose a variational approach for reconstruction and segmentation of CT images, with limited field of view and occluded geometry. The cost function s.t. a box constraint on the image values x and the simplex constraint on the labelling function v. The operator A is the undersampled Radon transform modelling the occluded geometry and y is the given data. The second term is the edge-preserving regularisation term for u, the third term is the segmentation term which aims at finding regions in u that are close to the value c k in region k. The operator D is the finite difference approximation of the gradient. The non-convex problem is solved by alternating minimisation between updates of u, v, c.
PET and Transmission Tomography. In [31], the authors propose a maximum likelihood reconstruction and doubly stochastic segmentation for emission and transmission tomography. In their model they use a Hidden Markov Measure Field Model (HMMFM) to estimate the different classes of objects from the given data r. They want to maximise the following cost function E(u, p, θ) = log P (r|u) + log P (u|p, θ) + log P (p).
The first term is the data likelihood which will be modelled differently for emission and transmission tomography. The second term is the conditional probability or class fitting term, for which they use HMMFM. The third term is the regularisation on the HMMFM. The optimisation is carried out in three steps, where first they solve for u (image update) fixing p, θ, then for p, holding u, θ (measure field update) and finally for θ (parameter update) having u, p fixed.
A variant of this method has been presented in [32], in which they incorporate prior information about the segmentation classes through a HMMFM. Here, the reconstruction is the minimisation over a constrained Bayesian formulation that involves a data fidelity term as a classical least squares fitting term, a class fitting term as a Gaussian mixture for each pixel given K classes and dependent of the class probabilities defined by the HMMFM, and a regulariser also dependent of the class probabilities. The model to minimise is The operator A will be modelled as the Radon transform in case of CT and b represents the measured data; N is the number of pixel in the image; λ noise and λ class are the regularisation parameters; µ k , σ k are the class parameters. The cost function is non convex and they solve the problem in an alternating scheme where they either update the pixel values or the class probabilities for each pixel.
Storath and others [33] model the joint reconstruction and segmentation using the Potts Model with application to PET imaging and CT. They consider the variational formulation of the Potts model for the reconstruction. Since the solution is piecewise constant, this directly induces a partition of the image domain, thus a segmentation.
Given the data f and an operator A (e.g. Radon transform), the energy functional is in the following form where the first term is the jump penalty enforcing piecewise constant solutions and the second term is the data fidelity. As the Potts model is NP hard, they propose a discretisation scheme that allows to split the Potts problem into subproblems that can be solved efficiently and exactly.

MRI.
In [34], the authors proposed a joint model with application to MRI. Their reconstruction-segmentation model consists of a fitting term and a patch-based dictionary to sparsely represent the image, and a term that models the segmentation as a mixture of Gaussian distributions with mean, standard deviation and mixture weights µ, σ, π. Their model is where A is the undersampled Fourier transform, y is the given data, R n is a patch extraction operator, λ is a weighting parameter, T is the sparsity threshold, and γ n is the sparse representation of patch R n u organised as column n of the matrix Γ. The problem is highly non-convex and it is solved iteratively using conjugate gradient on u, orthogonal matching pursuit on Γ and Expectation-Maximisation algorithm on (µ, σ, π).
Summary. Recently, the idea to solve the problems of reconstruction and segmentation simultaneously has become more popular. The majority of these joint methods have been proposed for CT, SPECT and PET data. Mainly they differ in the way they encode prior information in terms of regularisers and how they link the reconstruction and segmentation in the coupling term. Some imposes smoothness in the reconstruction [29], others sparsity in the gradient [26,30,33], other consider a patch-dictionary sparsifying approach [34]. In [33] they do not explicitly obtain a segmentation, but they force the reconstruction to be piecewise constant. Depending on the application, the coupling term is the data fitting term itself (e.g. SPECT), or the segmentation term. In [31,32,34] the authors model the segmentation as a mixture of Gaussian distribution, while [30] has a a region-based segmentation approach similar to what we propose. However, [30] penalises the squared 2-norm of segmentation, imposing spatial smoothness. In our proposed joint approach, we perform reconstruction and segmentation in a unified Bregman iteration scheme, exploiting the advantage of improving the reconstruction, which results in a more accurate segmentation. Furthermore, the segmentation constitutes another prior imposing regularity in the reconstruction in terms of sharp edges in the regions of interest. We propose a novel numerical optimisation problem in a non-convex Bregman iteration framework for which we present a rigorous convergence result in the following section.

Optimisation
The cost function (11) is non-convex in the joint argument (u, v), but it is convex in each individual variable. To solve this problem we derive a splitting approach where we solve the two minimisation problems in an alternating fashion with respect to u and v. We present the general algorithm and its convergence analysis in the next subsection. First, we describe the solution of each subproblem.
We now solve a variant of the primaldual method [35] as suggested in [38,39]. They consider the general problem including pointwise linear terms of the form min x∈C max y∈B Kx, y + g, x − h, y where C ⊆ X, B ⊆ Y are closed, convex sets. Setting K = ∇ and h = 0, θ = 1 and step sizes σ = τ = 0.99/ ∇ , the updates are At the end, we set v k+1 = v n+1 and obtain the update q k+1 as (11d).

Convergence analysis
The proposed joint approach (11) is an optimisation problem of the form In this work we consider a finite dimensional setting and we refer to the next section for the required definitons. To prove global convergence of (12), we consider functions that satisfy the Kurdika-Lojasiewicz property, defined below, and we make the following assumptions. holds.
• If F satisfies the KL property at each point of dom(∂F ), F is called a KL function.

Lemma 1
The j=1 v ij (c j − u i ) 2 in our joint problem (11) satisfies the KL property over R n × R m . Proof. It has been proved in [40] that real-analytic functions satisfy the KL property. The function E(u, v) is polynomial and therefore it is a real-analytic function.
Likewise for any fixed u, the function v → E(u, v) is C 1 L 2 (u) . We want to study the convergence properties of the alternating scheme for initial values (u 0 , v 0 ), p 0 ∈ ∂J 1 (u 0 ) and q 0 ∈ ∂J 2 (v 0 ). We want to show that the whole sequence generated by (13) converges to a critical point of E.
In order for the updates (13a) and (13c) to exist, we want J to be of the form J = R + εG (e.g. R = ∇u 1 and G = u 2 2 , see [41]) where R and G fulfil the following assumptions. In practice, we verify that G does not significantly change the reconstruction and segmentation performance for the examples we consider in the next section, for sufficiently small parameter (e.g. ε = 10 −3 ). Therefore, in our model (11) and in the numerical results we omit it.

Algorithm 1 Alternating splitting method with Bregman iterations for two blocks.
Initialization: Assumption 2 (i) The functions G 1 : R n → R and G 2 : R m → R are strongly convex with constants γ 1 and γ 2 , respectively. They have Lipschitz continuous gradient ∇G 1 and ∇G 2 with Lipschitz constant δ 1 and δ 2 , respectively.
(ii) The functions R 1 : R n → R and R 2 : R m → R are proper, l.s.c. and convex. For Theorem 1 (Global convergence). Suppose E is a KL function for any z k = (u k , v k ) ∈ R n ×R m and r k = (p k , q k ) with p k ∈ ∂R 1 (u k ), q k ∈ ∂R 2 (v k ). Assume Assumptions 1 and 2 hold. Let {z k } k∈N and {r k } k∈N be sequences generated by (14), which are assumed to be bounded. Then (ii) The sequence {z k } k∈N converges to a critical pointz of E.

Proof of Theorem 1
In the following we are going to show global convergence of this algorithm. The first step in our convergence analysis is to show a sufficient decrease property of a surrogate of the energy function (12) and a subgradient bound of the norm of the iterates gap. We first recall the following definitions.

Definition 2 (Convex Conjugate).
Let G be a proper, l.s.c. and convex function. Then its convex conjugate G * : R n → R ∪ {∞} is defined as Lemma 2 Let G be a proper, l.s.c. and convex function and G * its convex conjugate. Then for all arguments u ∈ R n with corresponding subgradients p ∈ ∂G(u) we know • p ∈ ∂G(u) is equivalent to u ∈ ∂G * (p).

From Lemma 2 we can rewrite the Bregman distance in (3) as follows
where we can see that now it does not depend on u k anymore, but it can be defined as a function of u and p k only, D J (u, p k ).

Definition 3 (Strong convexity).
Let G be a proper, l.s.c. and convex function. Then G is said to be γ-strongly convex if there exists a constant γ such that holds true for all u, v ∈ dom(G) and q ∈ ∂G(v).

Definition 4 (Symmetric Bregman distance).
Let G be a proper, l.s.c. and convex function. Then the symmetric generalised Bregman distance D symm for u, v ∈ dom(G) with p ∈ ∂G(u) and q ∈ ∂G(v). We also observe that in case G is γ-strongly convex we have Definition 5 (Lipschitz continuity). A function G : R n → R is (globally) Lipschitzcontinuous if there exists a constant L > 0 such that Before we show global convergence, we first define the surrogate functions.
Definition 6 (Surrogate objective). Let E, R i , G i , i ∈ {1, 2} satisfy Assumption 1 and Assumption 2, respectively. For any (u k , v k ) ∈ R n × R m and subgradients p k ∈ ∂R 1 (u k ) and q k ∈ ∂R 2 (v k ), we define the following surrogate objectives F , F 1 and F 2 For convenience we will use the following notations The surrogate function F will then read F (z k+1 , r k ) = F (u k+1 , v k+1 , p k , q k ).
We can now show the sufficient decrease property of (17) for subsequent iterates.
Lemma 3 (Sufficient decrease property). The iterates generated by (14) satisfy the descent estimate In addition we observe Proof. From (12) we consider the following step for Computing the optimality condition we obtain Taking the dual product with u k+1 − u k yields Using the convexity estimate Adding α 1 D p k−1 R 1 (u k , u k−1 ) to both sides, using the strong convexity of G 1 and the surrogate function notation, we get Using the trivial estimate for the Bregman distances, we get the decrease property Similarly for v, we obtain Summing up these estimates, we verify the sufficient decrease property (20), with positive ρ 2 = max{ε 1 γ 1 , ε 2 γ 2 }. We also observe Summing over k = 0, . . . , In order to show that the sequences generated by (14) approach the set of critical point we first estimate a bound for the subgradients of the surrogate functions and verify some properties of the limit point set. We first write the subdifferential of the surrogate function as with p k ∈ ∂R 1 (u k ) and q k ∈ ∂R 2 (v k ) being equivalent to u k ∈ ∂R * 1 (p k ) and v k ∈ ∂R * 2 (q k ), respectively.
Following [41,42], we verify some properties of the limit point set. Let {z k } k∈N and {r k } k∈N be sequences generated by (14). The set of limit points is defined as ω(z 0 , r 0 ) := (z,r) ∈ R n × R n : ∃ an increasing sequence of integers {k j } j∈N such that lim j→∞ z k j =z and lim j→∞ r k j =r .
As in [41,Definition 5.4,Proposition 5.5], we are going to assume that R i , i = 1, 2 has locally bounded subgradients.
Lemma 5 Suppose Assumptions 1 and 2 hold. Let {z k } k∈N be a sequence generated by (14) which is assumed to be bounded. Let (z,r) ∈ ω(z 0 , r 0 ). Then the following assertion holds lim Proof. Since (z,r) is a limit point of {(z k , r k )} k∈N , {(z k , r k )} k∈N , there exist subsequences {z k j } j∈N and {r k j } j∈N such that lim j→∞ z k j =z and lim j→∞ r k j =r, respectively. We immediately obtain due to the continuity of E and lim j→∞ D p k j −1 R 1 (u k j , u k j −1 ) = 0 and lim j→∞ D q k j −1 R 2 (v k j , v k j −1 ) = 0. From the sufficient decrease property we conclude (23). Lemma 6 (Properties of limit point set). The limit point set w(z 0 ) is a non empty, compact and connected set, the objective function E is constant on w(z 0 ) and we have lim k→∞ dist(z k , w(z 0 )) = 0. Proof. This follows steps as in [42,Lemma 5].
To finally prove global convergence of (14), we will use the following Kurdyka-Lojasiewicz property defined and the result from [42]. Before recalling the definition, we introduce the notion of distance between any subset S ⊂ R d and any point x ∈ R d defined as where · denotes the Euclidean norm.

Lemma 7 (Uniformised KL property).
Let Ω be a compact set and let E : R n ×R m → R ∪ {∞} be a proper and l.s.c. function. Assume that E is constant on Ω and satisfy the KL property at each point in Ω. Then there exists ε > 0, η > 0 and ϕ ∈ C 1 ((0, η)) that satisfies the same conditions as in Definition KL, such that for allū ∈ Ω and all u in condition KL is satisfied. Proof. Follows from [42].
With these results we can now show global convergence of (14). Proof of Theorem 1.
By the boundedness assumption on {(z k , r k )} k∈N , there exist converging subsequences {z k j } j∈N and {r k j } j∈N such that lim j→∞ z k j =z and lim j→∞ r k j =r, respectively. We know from Lemma 5 that (23) is satisfied.
(i) KL property holds for E and therefore for E k and we write and from the concavity of ϕ we know that Thus, we obtain From (20) with Lemma 3 and using the abbreviation Multiplying by z k − z k−1 and using Young's inequality (2 Summing up from k = 1, . . . , N we get In addition we observe that the finite length property implies that the sequence {z k } k∈N is a Cauchy sequence and hence is a convergent sequence. For each z r and z s with s > r > l we have (ii) The proof follows in a similar fashion as in [41,Lemma 5.9] Remark 2 (Extension to d blocks). The analysis described above holds for the general setting of d blocks The update for each of the d blocks then reads

Numerical results
In this section we present numerical results for our joint reconstruction and segmentation model described in (11). We demonstrate its advantages and limitations, as well as a discussion on the parameter choice. In the first part, we focus on bubbly flow segmentation for simulated data. In the second part, we show results for real data acquired at the Cancer Research UK, Cambridge Institute, for tumour segmentation.
Quality measure. To assess the performance of the reconstruction we will compare our solutions u with respect to the groundtruth u gt . As quality measure we use the relative reconstuction error (RRE) and the peak signal to noise ratio (PSNR) defined as • PSNR(u, u gt ) = 10 log 10 max(u) u gt −u 2 /N For the segmentation quality, we will use the relative segmentation error (RSE) to compare our segmentations v with respect to the true segmentations v gt where N is the number of pixels in the image, δ is the Kronecker delta function that will count the number of mis-classified pixels.
Before we present our two applications, we show a more detailed result of the phantom brain in Figure 1. In this example, we show the TV reconstruction 2b, where the parameter α has been optimised with respect to PSNR and its sequential segmentation 2f with optimal β with respect to RSE. In 2c and 2g we present Bregman reconstruction and sequential segmentation where the Bregman iteration has been stopped according to the discrepancy principle Equation 7 and β has been optimised with respect to RSE. These parameter choices for the sequential approaches will be used  Figure 2: We consider 15% of the simulated k -space for the brain phantom, where Gaussian noise (σ = 0.25) was added. We compare results for the total variation reconstruction and total variation based Bregman iterative reconstruction and their segmentation in a sequential approach with our joint model. We show that both reconstruction and segmentation are improved.
in the whole paper.
In this first result, we clearly see that the joint approach performs much better compared to the separate steps in Figures 2b, 2f and 2c, 2g. Both reconstruction and segmentation are improved and more details are recovered. We refer to Appendix A for more simulated examples.

Bubbly flow
The first application considered is the characterisation of bubbly flows using MRI. Bubbly flows are two-phase flow systems of liquid and gas trapped in bubbles, which are common in industrial applications such as bioreactors [43] and hydrocarbon processing units [44]. MRI has been successfully used to characterise the bubble size distribution [45,46] and the liquid velocity field of bubbly flows [47,48]; these properties govern the heat and mass transfer between the bubbles and the liquid which ultimately determine the efficiency of these industrial systems. However, when studying fast flowing systems, the acquisition time for fully sample k -space is too long to resolve the temporal changes; the most common method of breaking the temporal resolution barrier is through under-sampling. It is therefore critical to develop reconstruction techniques for highly under-sampled k -space data for the accurate reconstruction of the MRI images which would be subsequently used in calculating the bubble size distribution or in  studying the hydrodynamics of the system. We apply our joint reconstruction and segmentation approach to simulated bubbly flow imaging. In Figure 3 we present some results for synthetic data, where Figure 3a represents the groundtruth magnitude image, from which we simulate its k -space following the forward model described in (1). From the full k -space we collect 8% of the samples using the sampling matrix in Figure 3e and we corrupt the data with Gaussian noise of standard deviation σ = 0.35. In Figure 3b and 3f we show the results for the total variation regularised reconstruction and its segmentation performed sequentially. In the same sequential way, we show the results for the Bregman iterative regularization in Figure 3c and 3g. In the last column in Figure 3d and 3h, we finally show the results for our joint approach. Although the TV and the Bregman approaches are already quite good, we can see that both RRE and PSNR are improved using our model in the reconstruction and the segmentation. Smaller details, such as the top right bubble contour, are better detected when solving the joint problem. As the goal of the bubbly flow application is to detect bubble size distribution, this improvement is really advantageous.
We tested the robustness of our approach by corrupting the data with different signal to noise ratio (SNR) and by considering different amount of sampling. In this information, we show in Figure 6 how the PSNR, RRE and RSE are affected, for the joint approach (blue lines) and for the separate approaches, TV (red dotted lines) and Bregman TV (green dotted lines). As expected, with the SNR increasing the error decreases. We can see that the joint approach performs better than the sequential approach for any SNR. The improvement is even more significant for very noisy data. As in practice we often observe high levels of noise, the joint approach is able to takle this problem better than the traditional sequential approaches. It is also interesting to investigate how the joint approach performs with very low undersampling rates. In Figure 3e we show joint reconstructions (top row) and corresponding segmentations (bottom row) for decreasing sampling rates. We can see that up to 5% results are still very good. Using 3 and 2% of the samples the results are less clean but it is possible to identify the main structures. In contrast, 1% sampling is not enough to retrieve a good image reconstruction and consequently its segmentation. In Figure 7, we plot PSNR , RRE and RSE for different sampling rates. The blue lines represent the error for our joint approach, while the red and green dotted lines are for the sequential TV and sequential Bregman TV approaches. We can see that up to 5% sampling the error measures do not change significantly. However, for lower rates, the improvement is more significant. This is highly beneficial for the bubbly flow application as increasing the temporal resolution is really important to keep track of the gas flowing in the pipe.  From left to right, we show the PSNR, RRE and RSE, respectively, for different levels of noise in the measurements. The blue lines represent the error for our joint approach, while the red and green dotted lines are for the sequential TV and sequential Bregman TV approaches. For each SNR, the joint model performs better than the separate methods. This improvement is even more significant for noisier data, which is highly advantageous as in practice we often observe lower SNR.

Cancer imaging
In this subsection, we illustrate the performance of the joint model for real cancer data. At the Cancer Research UK, Cambridge Institute, researchers acquire every day a huge amount of MRI scans to assess tumour progression and response to therapy [49]. For this reason, it is very convenient to have fast sampling through compressed sensing, and automatic segmentation methods. Furthermore, reconstructions with enhanced edges are advantageous to facilitate clinical analysis.
Here we show our results for MRI data of a rat bearing a glioblastoma. The MR image represents the rat head where the brain is the gray area in the top half of the image. Inside this gray region, a tumour is clearly visible appearing as a brighter area. For this The blue lines represent the error for our joint approach, while the red and green dotted lines are for the sequential TV and sequential Bregman TV approaches. The joint appraoch performs better than the sequential cases. The gain is not very significant for higher sampling rates, but it becomes more important for lower rates, starting from 3%. We select 15% of the samples using a spiral mask. The image show a rat brain bearing a tumour (brighter region). The zero-filled reconstruction 8a and the TV regularised reconstruction 8b are shown together with their sequential segmentation 8e and 8f respectively. In the last column 8d and 8h we show the results for our model. The parameter α for the TV reconstruction and for the joint reconstruction has been chosen such that it achieves visually optimal in the sense that it resolve all the details (e.g. the darker line cutting the tumour transversally).
experiment, we acquired the full k -space and present the zero-filling reconstruction in Figure 8a and the sequential segmentation in Figure 8e. As discussed already in the previous section, the zero-filled reconstruction presents noise and artefact which may complicate the segmentation. We want to show that the compressed sensing approach and in particular the joint model can improve this reconstruction. Given the full k -  space, we select 15% of the samples using a spiral mask. In Figure 8b, 8f and Figure 8c, 8g we show the results for the sequential approaches. In Figure 8d and 8h we show the joint reconstruction and the joint segmentation obtained for the same data. The regularised approaches perform better that the zero-filled reconstruction, producing less noisy results. However, our joint model is able to produce a cleaner reconstruction where the edges that defines the tumour and the brain are very well detected. In Figure 9, we show a zoomed section where it is easy to assess that the joint model tackle the noise and detect the region of interest. We can see that we are able to improve the reconstruction and automatically identify the tumour in the brain. The degree of enhancement of the edges in the reconstruction is controllable by the parameter δ in the model (11). In the next subsection we present a discussion on how to tune this parameter.

Parameter choice rule
In the model proposed in (11), the parameters that we need to choose are α, β and δ.
In this section we discuss a rule to choose them depending on the desired results. Some examples will clarify these empirical choices.
• α balances the total variation regularization term in the reconstruction for the magnitude. The higher the α, the more piecewise constant the reconstruction will be. See Figure 10.
• β defines the scale of the objects that will be detected in the segmentation. Smaller values of β will allow for smaller objects. See Figure 11.
• δ is the parameter linking the reconstruction and the segmentation. To better illustrate its role, let us consider a zero-filling like reconstruction and segmentation: (a) α = 0.001 (b) α = 0.01 (c) α = 0.1 Figure 10: The parameter α balances the data fidelity term and the total variation regularisation for the reconstruction. Smaller values of α produce a reconstruction closer to the data fitting term, hence less smooth as in 10a. As α increases in 10b the solution gets smoother and less noisy. Finally for large values it tends to become more piecewise constant as in 10c. where . This problem is solving the zero-filled reconstruction and segmentation jointly. For δ = 0, the reconstruction is the zerofilling solution. In Figure 12 we can see the impact of the segmentation term on the reconstruction for increasing values of δ. We can see that for very smallδ the result is close to the zero-filling solution. For δ = 1 the noise from the model is present as expected but in addition the boundaries are enhanced. For large δ the boundaries are still very pronounced and the noise is also amplified.

Comparison with another joint approach
We present a comparison of our joint model with another non-convex method, namely the Potts model approach by [33], described in Subsection 2.4. The major advantage of the joint reconstruction and segmentation using the Potts model is that it does not require to select explicitely the number of regions to segment, although this depends on Figure 12: We show the reconstructions obtained solving (26) for different values of δ. For δ = 0 we get the zero-filling solution. For small δ we expect the solution to be similar to the zero-filling reconstruction. For δ = 1 we see the effect of the joint term on the reconstruction. The solution presents the same noise artefacts but having in addition very sharp boundaries. Finally, for very large δ we still have enhanced boundaries but we also amplify the noise.
the choice of the regularisation parameter. However, by definition, it only produces a piecewise constant image, therefore a segmentation, and not a reconstruction. This is useful in some applications where one is only interested in the segmentation. In contrast, our model produces both reconstruction and segmentation. In Figure 13, we show the results for some examples. Note that because the results of the Potts model are in the range of the groundtruth image, while our segmentation are in label space, we can not directly use the RSE as before, or common metrics that compare actual intensities such as PSNR and structure similarity index measure (SSIM), for comparison. For example, for some tissue in class 1, to label it class 2 is as wrong as to label it class 3. However in this case, the SSIM and PSNR will favour the label class 2.
We therefore focus on a visual assessment and show the results of the Potts model for two different choices of the regularisation parameter γ. We recall that the proposed model requires to determine the number of classes in advance, while the model for comparison estimates the number of regions but this depends on the choice of the regularisation parameter. In the top row, we can see that the Potts model, although it retrieves the shape of the main structures for the brain phantom example, it overestimates the number of classes. By increasing the parameter γ, this issue is not resolved as it assigns different intensities to objects of the same category. In contrast, our approach is able to identify the desired classes as in the groundtruth. For the bubble case (middle row), we can see that our method works better and our segmentation is more accurate, while the Potts model fails to capture shape details (e.g., outer circle is distorted) and again overestimates the number of regions. We can also see that, when slightly decreasing γ, the Potts model is very sensitive to artefacts. For the real MR data (bottom row), we see that both methods identify the tumour quite well. Because we were only interested in identifying three classes as tumour, brain and background, we do not segment the outer region (rat's head), captured insted by the Potts model. However, the Potts model only produces the segmentation, while our method, as shown in Figure 8, also produces an enhanced reconstruction with sharp edges.

Conclusion
In this paper, we have investigated a novel mathametical approach to perform simultaneously reconstruction and segmentation from undersampled MRI data. Our motivation was to include in the reconstruction prior knowledge of the objects we are interested in. By interconnecting the reconstruction and the segmentation terms, we can achieve sharper reconstructions and more accurate segmentations. We derived a variational model based on Bregman iteration and we have verified its convergence properties. With our approach we show that by solving the more complicated joint model, we are able to improve both reconstruction and segmentation compared to the traditional sequential approach. This suggests that with the joint model it is possible to reduce error propagations that occur in sequential analysis, when the segmentation is separate and posterior to the reconstruction. We have tested our method for two different application, which are bubbly flow and cancer imaging. In both cases, the reconstructions are sharper and finer structures are detected. Additionally, the segmentations also benefit from the improvement in the reconstructions. We have found that the joint model outperforms the sequential approach by exploiting prior information on the objects that we want to segment. In addition, we also show that our method performs better than the well-known Potts model. We also presented a discussion on the parameter choice rule that offer some insight on how to tune the parameters according to the desired result. It is interesting to notice that, with our model, we are able to control the segmentation effect on the reconstruction. Furthermore, when the final analysis of the MR image is indeed the segmentation, it is possible to bias the reconstruction towards the piecewise constant solution, yet preserving finer details in the structure. In our set-up, we have specified the intensity constants characteristic of the region of interests, which were known a priori for our applications. However, it is possible to also include the optimisation with respect to c j in our joint model, where the same convergence guarantees hold (see Remark 2). Nevertheless, one limitation of the model is the need to specify the number of regions to be segmented.
In our future research, we would like to study the extension of this model for the bubbly flow to the reconstruction of the magnitude image as well as the phase image. The goal is not only to extract the structure of the bubble, but also to estimate velocity information, which is encoded in the phase image. As the problem is non-convex in the joint argument but also non-convex with respect to the phase, we need to derive a different convergence analysis.