MUMFORD–SHAH-TV FUNCTIONAL WITH APPLICATION IN X-RAY INTERIOR TOMOGRAPHY

,


1.
Introduction. The primary motivation of this research came from our attempt to further explore the interior problem of X-ray CT. This problem, which arises from technical limitations of the scanning apparatus and requirement of X-ray dose reduction, amounts to reconstructing the object image in a region of interest (ROI) from truncated projection data (interior data) pertaining to the ROI. Conventional wisdom asserts that the interior problem is severely ill-posed since it is not uniquely solvable [18,28]. Yet microlocal analysis indicates that the singularities of the object image in the ROI are visible in a stable way [32] from interior data by singular value decompositions [25,26] or Lambda tomography [15]. With recent advancement of CT theory, one can pursue the accurate reconstruction of the ROI with interior data by incorporating some priori knowledge of the object image. Specifically, if the object image in a subregion of the ROI is known a priori, then the ROI can be exactly reconstructed via analytic continuation method [24,41]. However, this condition is just satisfied on few occasions so as to narrow down the application of this method. A universal priori knowledge is the sparsity of image representation, which leads to the total variation (TV) minimization method for the X-ray interior problem. The ROI can be reconstructed exactly by TV minimization method under the assumption that the object image is piecewise constant in ROI [42,43]. As an extension to TV minimization method, high order total variation (HOT) minimization method was proposed to reconstruct a piecewise polynomial ROI image [40].
Since TV was introduced for image restoration [35], TV minimization method has gradually become a common tool for imaging and image processing problems in the compressed sensing scheme. For piecewise smooth image, the total variation includes two parts: edge and smooth terms, where |u + − u − | denotes the image jump across the discontinuity set S u of u . TV regularization method can be represented as the following unconstrained optimization problem: (2) min u αTV(u) + Au − g 2 .
Various algorithms, such as gradient projection method [35], second-order method [38] and split Bregman iteration [17], were proposed to resolve this problem. TV regularization is well known for recovering sharp edges of object image. But the edge set is not a direct result of this method. Edge information is always helpful in practical applications. The goal of this investigation is to present an image reconstruction scheme to separate the edge from the smooth part of the image while taking advantage of TV regularization, and then apply it to X-ray interior problem. As a successful approach to deal with the image edge, Mumford-Shah (MS) functional was originally proposed for image denoising and segmentation [27], where A is the identity operator, SBV(Ω) is the space of all special functions of bounded variation [2]. Since for any u ∈ SBV(Ω), S u is L N -negligible, the integral domain of ∇u is written as Ω in MS(u). MS functional aims to approximate the object image by a piecewise smooth function. With A being a forward operator to the measurement data g, MS functional and its variants have been applied to some other image processing and imaging problems, such as image inpainting [14], image deblurring [4], electric impedance tomography [34], X-ray tomography [33], electron tomography [22] and SPECT [23]. The regularizing properties of MS functional have been well studied in previous works [13,20]. Since MS functional in its original setting is quite sophisticated, some approximate formulations, such as level set [11], graph cut [8] and convex relaxation [31], were employed in various problems for the purpose of numerical realization. A more popular approach is to approximate the MS functional by Blake-Zisserman energy [7,10] or Ambrosio-Tortorelli elliptic functional [3] in the sense of Γ-convergence [9]. For the object image of piecewise constant, we endeavor to recover the image and its edge simultaneously. Incorporating MS functional with TV, we propose a new functional named as Mumford-Shah-TV (MSTV), where SBV c a,b (Ω) is a subset of SBV(Ω). Comparing with the results about the MS functional in [20], we introduce an extra constraint Ω |∇u| 2 dx ≤ c in SBV c a,b (Ω) to guarantee the existence of solutions to MSTV functional, where c > 0 is a constant. We systematically study the existence and stability of the minimizers of MSTV functional. Accordingly, we present an Ambrosio-Tortorelli type approximation for MSTV functional in the sense of Γ-convergence. Some related work were reported to directly apply L 1 regularization to the Blake-Zisserman [16,19] and Ambrosio-Tortorelli [36] type approximations of MS functional. In a different way, we apply L 1 regularization in MS functional directly. MSTV regularization is prone to reconstruct a piecewise constant image and preserve the sharp edges. But the reconstructed image is not restricted to piecewise constant as the levelset method for MS regularization. We apply the MSTV regularization method to the X-ray interior problem and develop an algorithm based on split Bregman and OS-SART iterations.
The remaining part of this paper is organized as follows. In §2, we define the MSTV functional and study the existence and stability of the minimizers of MSTV functional. We present an Ambrosio-Tortorelli type approximation for it in the sense of Γ-convergence and apply it to the X-ray interior problem. In §3, physical experiments are used to validate the MSTV regularization method and the algorithm we proposed. Finally, we discussed some related issues and conclude in §4.

2.
Mumford-Shah-TV functional and its application in X-ray interior problem.
2.1. Mumford-Shah functional. We will always assume that Ω is a bounded open subset of R N . A function u ∈ L 1 (Ω) is bounded variation if the distributional derivative Du is representable by a finite Radon measure in Ω [2]. The space of functions of bounded variation is denoted by BV(Ω). We say a function u ∈ L 1 (Ω) is approximately continuous at x ∈ Ω if there existsũ(x) ∈ R such that where B ρ (x) = {y ∈ Ω : |y −x| < ρ}, andũ is called the approximate limit of u. The set of points without an approximate limit is L N -negligible and is denoted by S u . If u ∈ BV(Ω), for H N −1 -a.e. x ∈ S u , there exist u + (x), u − (x) ∈ R and ν ∈ S N −1 such that ∇u is called the approximate gradient of u. For any u ∈ BV(Ω), the distributional derivative Du can be decomposed into three parts Here, ∇uL N and (u + − u − )ν u H N −1 S u are the Lebesgue integrable and jump parts respectively, and D c u is the so-called Cantor part. A function u ∈ BV(Ω) is called a special function of bounded variation if D c u = 0. The space of all special functions of bounded variation is denoted by SBV(Ω) and has closure and compactness theorem described as follows.
Theorem 2.1 (Closure and compactness theorem in SBV(Ω), [2]). Let p > 1 and {u k } k∈N be a sequence in SBV(Ω) satisfying that Then there exist a subsequence, still denoted by {u k } k∈N , and some u ∈ SBV(Ω) such that Let Θ be an open subset of R M , A : L 2 (Ω) → L 2 (Θ) a continuous operator, g ∈ L ∞ (Θ) and The MS functional for imaging and image processing problems is defined as [20] where α and β are two positive parameters, · denotes the L 2 -norm throughout this paper. One can obtain a piecewise smooth approximation of the object image by minimizing MS(u) in SBV(Ω) ∩ X b a (Ω). MS regularization is algorithmically challenging because of the different nature of u and S u . A popular approach is to approximate MS(u) by a sequence of simple functionals in the sense of Γ-convergence.
Definition 2.2 (Γ-convergence, [9]). Let X be a metric space, we say that F ε : X → [−∞, +∞] Γ−converges to F : X → [−∞, +∞] in X as ε → 0, if for any sequence {ε k } k∈N converging to 0 and u ∈ X, the following two conditions are satisfied, 1. for every sequence {u k } k∈N ⊂ X converging to u, An approximation by Ambrosio-Tortorelli elliptic functional is given as follows: otherwise. and otherwise. Then is usually minimized by an alternating scheme: It is worth noting that there is a diffusion term Ω (v k |∇u|) 2 dx in (18). The edge may be blurred after some iterations in the numerical realization [21].

2.2.
Mumford-Shah-TV functional. TV regularization is widely used in imaging and image processing problems particularly for the object image of piecewise constant. Our goal is to recover a piecewise constant image while taking into account the edge information as well. Incorporating MS functional with TV, we define a new functional named as Mumford-Shah-TV, where c is a positive constant and Without loss of generality, we will always assume that a, b and c are three fixed positive constants. The condition of the closure and compactness theorem 2.1 should be satisfied by at least one minimizing sequence of MSTV(u). Otherwise, the minimizer of MSTV(u) may not exist. Hence, we introduce an extra constraint Ω |∇u| 2 dx ≤ c to the image u, and this constraint can always be satisfied in practical situations. Then, we obtain the existence of one minimizer of MSTV(u) in SBV c a,b (Ω). Theorem 2.4. There exists at least one minimizer of MSTV(u) in SBV c a,b (Ω).
Proof. Let {u k } k∈N be a minimizing sequence of MSTV(u) in SBV c a,b (Ω), namely, (Ω), by Theorem 2.1, there exist a subsequence, still denoted by {u k } k∈N , and some u * ∈ SBV(Ω) such that Since Ω is bounded and The weak convergence ∇u k ∇u * in [L 2 (Ω)] N also implies weak convergence in [L 1 (Ω)] N . By the semicontinuity of weak convergence, we obtain that It follows that u * ∈ SBV c a,b (Ω) and (29) MSTV Thus, u * is a minimizer of MSTV(u) in SBV c a,b (Ω).
In order to study the stability of the minimizers of MSTV(u), we introduce the weak convergence in SBV(Ω). Definition 2.5 (Weak convergence in SBV(Ω), [12]). We say that a sequence {u k } k∈N ⊂ SBV(Ω) weakly converges to u in SBV(Ω) if u ∈ SBV(Ω) and Let D 1 and D 2 be two subsets of Ω, we denote D 1⊂ D 2 if H N −1 (D 1 \D 2 ) = 0 and denote D 1= D 2 if D 1⊂ D 2 and D 2⊂ D 1 . Then, we introduce a notion of convergence of sets based on the weak convergence in SBV(Ω). Definition 2.6 (σ-convergence, [12]). Let {E k } k∈N be a sequence satisfying that E k ⊂ Ω for every k ∈ N and sup k∈N H N −1 (E k ) < +∞. We say that {E k } k∈N σ-converges to E ⊆ Ω in Ω if the following conditions are satisfied, 1. if {u i } i∈N weakly converges to u in SBV(Ω) and S ui⊂ E ki for some subsequence {E ki } i∈N , then S u⊂ E; 2. there exists a sequence {u k } k∈N ⊂ SBV(Ω) weakly converging to some u in SBV(Ω) such that S u= E and S u k⊂ E k for every k ∈ N.
By the same argument as used in [20], we obtain the stability of the minimizers of MSTV(u) based on the weak convergence in SBV(Ω) and σ-convergence of sets.
Theorem 2.7. Let {g k } k∈N ⊂ L ∞ (Θ) be a sequence converging to g in L 2 (Θ), and u k a minimizer of Then there exist a subsequence, still denoted by {u k } k∈N , and some u * ∈ SBV c a,b (Ω) such that {u k } k∈N weakly converges to u * in SBV(Ω), {S u k } k∈N σ-converges to S u * , and u * minimizes MSTV(u; g) in SBV c a,b (Ω). Proof. The proof is analogous to the proof of Lemma 8 of [20] and is omitted here.
An approximation formulation of MSTV(u) is necessary for the minimization of MSTV functional to address the problem due to the edge term H N −1 (S u ).
Proof. Let {η k } k∈N be a positive sequence converging to 0 and u k a solution of the following problem (34) min By the regularizing properties of MS functional [20], there exists at least one solution for problem (34) and H N −1 (Ω ∩ (S u k \ S u k )) = 0. Using the theory of constrained problems [5], we have u k → u in L 2 (Ω). For every k ∈ N, By Theorem 2.1, we can extract a subsequence, still denoted by {u k } k∈N , such that ∇u k ∇u weakly in [L 2 (Ω)] N and Together with (35), we deduce that If there exists a subsequence, still denoted by {u k } k∈N , such that Ω |∇u k | 2 dx ≤ c, then by (39), (41) and the continuity of A, we conclude that u k ∈ SBV c a,b (Ω) and MSTV(u k ) → MSTV(u) as k → ∞. Otherwise, we may assume that Ω |∇u k | 2 dx > c for every k ∈ N. By (38), we have Ω |∇u k | 2 dx → c. Letũ k = ζ k u k + (1 − ζ k )a with (42) ζ , and ζ k → 1. It is easy to check thatũ k ∈ SBV c a,b (Ω) for every k ∈ N. Then, we have By (43) and (44), we haveũ k → u in L 2 (Ω) and Ω |∇ũ k | dx → Ω |∇u| dx. It follows that MSTV(ũ k ) → MSTV(u) as k → ∞ and the proof is completed.
Then, we present an Ambrosio-Tortorelli-type approximation for the MSTV functional in the sense of Γ-convergence.  and Then, G ε Γ-converges to G in [L 2 (Ω)] 2 as ε → 0.
Proof. For any ε k → 0, we may assume that sup k∈N G ε k (u k , v k ) < +∞. Since Moreover, let ξ ∈ S N −1 and D be an open subset of Ω, for L N −1 -a.e. z ∈ ξ ⊥ and any δ > 0, there exists a closed set K δ z,ξ ⊂ R such that L 1 (K δ z,ξ ) < δ anḋ where u z,ξ k (t) = u k (z + tξ),u z,ξ k denotes the approximate gradient of u z,ξ k , and D z,ξ = {t ∈ R : z + tξ ∈ D} with D z,ξ = ∅. By (50) and (51), we have for L N −1a.e. z ∈ ξ ⊥ , both {(v z,ξ k ) 2u z,ξ k } k∈N and {v z,ξ ku z,ξ k } k∈N weakly converge tou z,ξ in L 2 (D z,ξ \ K δ z,ξ ), which implies that Let δ → 0, we obtain that Using Fatou's lemma and (54), we have Using the same argument to (55), we obtain that By the continuity of A, we have Au k − g 2 → Au − g 2 . Then, we conclude that u ∈ SBV c a,b (Ω) and obtain the lim inf-inequality. We now prove the lim sup-inequality and assume that G(u, v) < +∞. We first consider the case that H N −1 (Ω ∩ (S u \ S u )) = 0. Let d(x, S u ) = inf{|x − y| : y ∈ S u } and define the following two functions, otherwise. (59) It is easy to verify that which follows that (u k , v k ) ∈ H c a,b (Ω). By the same argument as used in the proof of Proposition 5.2 of [3], we have that Using again the continuity of A, we have Au k −g 2 → Au−g 2 . By (61) and (63), we obtain the lim sup-inequality. Finally, we complete the proof by applying Lemma 2.8 and standard diagonal argument for the case that H N −1 (Ω ∩ (S u ∩ S u )) > 0.

2.3.
Application in X-ray interior problem. We now apply MSTV functional to the X-ray interior problem in R 2 . The object image u 0 is assumed to be compactly supported and the ROI is a disk Ω r = {x ∈ R 2 : |x| < r}. The X-ray interior problem is to reconstruct the ROI from the truncated projection data (64) Ru 0 (θ, s) = θ·x=s u 0 (x) dx, θ ∈ S 1 , |s| < r.
By Theorem 2.9, we reconstruct the ROI by solving the following problem: Here, v is the edge image of u and indicates the edge set S u . In discrete settings, u, v and g are usually represented as vectors of two indexes: where ∆ and ∆ d denote the space sampling interval of the image and the detector respectively, and ∆ a denotes the angular sampling interval of the scanning orbit. Correspondingly, the discrete projection operator is a matrix A = {a (k,l),(i,j) } (N d Na)×(NxNy) . The discrete representations of X b a (Ω) and H c a,b (Ω) are given byX where ∇u = (δ x u, δ y u) is the first order finite difference with Then, the discrete representation of (66) can be written as Combining with the algorithm we proposed to solve the L 1 -regularization problem in CT imaging in [45], we resolve (73) by the following scheme: where 0 < η < 1, and S(u k ) denotes the result of one OS-SART iteration [39] to solve Au = g. P K denotes the projection operator onto a convex set K, and (78)V k+1 = v ∈ R NxNy : v∇u k+1 2 ≤ c .
We solve problems (74) and (76) using split Bregman iteration and conjugate gradient (CG) method respectively. The projection ontoX b a can be performed by simple point-wise truncation: The projection ontoV k+1 can be represented as the following problem, By the method of Lagrange duality [6], (80) is equivalent to the dual problem, is the Lagrange dual function of (80). Problem (81) can be solve using Newton method [29].
Together with the result in Theorem 2.3, and using the same scheme as (74)-(77), MS functional can be minimized as follows, where (83) where u 0 is the object image, µ u and σ u are the mean value and standard deviation of image u within the ROI, respectively, and E rec and E SSIM are used to evaluate the difference and similarity between u and u 0 .
In practical applications, parameters a and b can usually be estimated using prior knowledge, and an appropriate choice of α and β is made in an empirical way. However, an accurate estimate of parameter c is usually difficult and a big value is preferred. Particularly, c is set to be +∞ in our experiments. The regularization parameter settings are listed in Table.3.
Reconstruction results by OS-SART iteration, TV, MS and MSTV regularization methods are shown in Fig. 4. The reconstructed images of TV and MSTV regularization are piecewise constant approximately and preserve the sharp edges of the object image. The dot matrix within the phantom implies that TV and MSTV regularization methods have the best performance in spatial resolution for object image of piecewise constant. The edge of the object image within the ROI was captured exactly by both MS and MSTV regularization method. Since the edge is smoothed while solving the diffusion problem (83), the edge captured by MS regularization is broadened. The curves of E rec and E SSIM are plotted in Fig.  5, which show that the image reconstructed by MSTV regularziation is the most accurate one in these two criterions.  Physical experiments are also conducted to test our method. A frozen chicken is scanned using cone beam geometry, and the middle slice is reconstructed using fan beam projection data. 720 views are measured in angle range [0, 2π), and the other geometry parameters of scanning are the same as the numerical simulation. The object image can be modeled as a piecewise constant function approximately, and are scanned for twice in 130 kVp, 3.5 A and 130 kVp, 0.5 A, respectively.  Table 1. Parameter settings of numerical and physical experiments.
The reconstructed images are represented by a 512 × 512 matrix with width of 140 mm. The ROI is a disk of radius 48.45 mm. The reconstructed images using non-truncated projection data and 10 OS-SART iterations are shown in Fig. 3.
While reconstructing the ROI using data of normal dose, we further reduce the scale of data to 180 × 254, namely, 180 views and 254 detector units are used to reconstruct the ROI. It makes the reconstruction more challenging. Reconstruction results of these four regularization methods are shown in Fig. 4. The edge image of MSTV regularization contains more details of the object image compared with that of MS regularization in which some edge with small jump size is missed. The reconstructed images of low-dose data are shown in Fig. 5. Compared with the result of normal dose, the reconstruction quality of low-dose data decreases significantly due to strong noise. It is well known that both TV and MS regularization have property of noise reduction. Since MSTV functional shares many common properties with them, it is understandable that MSTV regularization possesses the property of noise reduction.
4. Discussion and conclusion. MSTV regularization can also be incorporated with some other methods by applying different data term, such as Poisson-ROF denoising model for CT imaging [46] and L 1 -norm data term [37]. Blake-Zisserman type approximation and more regularization properties of MSTV functional will also be considered in our future work.
In this paper, we defined the MSTV functional and apply it to the interior problem of X-ray CT. We study the existence and stability of minimizers to the MSTV functional. We present an Ambrosio-Tortorelli type approximation for the MSTV  functional in the sense of Γ-convergence for the purpose of numerical realization. An algorithm based on split Bregman iteration and OS-SART is developed for the MSTV regularization method. Numerical and physical experiments demonstrate that both the image and its edge can be reconstructed with high-quality.