Limited data X-ray tomography using accretive operators

We prove that a nonlinear evolution equation which gives a novel approach to the X-ray tomography problem (see Kolehmainen et al., SIAM J. Sci. Comput. 30(3):1413–1429, 2008) has a solution. To this end, we list some of our results on theory of accretive operators and then we apply them to this concrete context.


Introduction
In medical X-ray tomography, the inner structure of a patient is reconstructed from a collection of projection images. The widely used computerized tomography (CT) imaging uses an extensive set of projections acquired from all around the body. This type of reconstruction is well understood, the most popular method being filtered back-projection (FBP).
In mammography and intraoral dental imaging, the X-ray detector is in a fixed position behind the tissue, and the X-ray source moves with respect to the detector. In these cases the projections can be taken from a view angle significantly less than 180 • , leading to a limited angle tomography problem. In some applications, such as the radiation dose to the patient is minimized by keeping the number of projection small. In addition, the projections are typically truncated to detector size, yielding a local tomography problem. We refer to the above types of incomplete data as sparse projection data.
Sparse projection data do not contain enough information to describe completely the tissue, and therefore successful reconstruction requires some form of regularization or priori information. It is well known that methods such as (FBP) are not suited for sparse projection data (see [25,26]). More promising approaches include, among others, total variation methods (see [7,8]), variational methods (see [20]), and deformable models (see [27]). Kolehmainen et al. in [19] study a variant of the level set method, where the X-ray attenuation coefficient is modeled as the function max{ (x), 0}, with a smooth function. Thus they make use of the natural "a priori" information that the X-ray attenuation coefficient is always nonnegative (the intensity of X-ray does not increase inside the tissue).
They assume that the attenuation coefficient v ∈ L 2 ( ) for a bounded subset ⊆ R 2 and use the following model for the direct problem: where A is a linear operator on L 2 ( ) with appropriate target space and ε is a measurement of the noise. To reconstruct v approximately from m, Kolehmainen et al. in [19] solve numerically the evolution equation with a suitable initial condition φ(x, 0) = φ 0 and r ≥ 0, β > 0. Here ν is the interior normal vector of the boundary. The cutoff function f : R → R is given by and consider the function as the reconstructed attenuation coefficient. The theory of accretive operators will allow us to show that the evolution equation (1) has a solution when the measurement equation (2) comes from one of the two most popular models for X-ray tomography: the pencil beam model or the Radon Transform. Moreover, we prove that the limit (4) exists when β is large enough and r > 0.

X-ray measurement models
In medical X-ray imaging, an X-ray source is placed on one side of the target tissue. Radiation passes through the tissue, and the attenuated signal is detected on the other side.
A two-dimensional model slice through the target tissue given by a rectangle ⊆ R 2 and a nonnegative attenuation coefficient v : The tissue is contained in a subset 1 ⊂ , and v(x) = 0 for x ∈ \ 1 . This yields the linear model where L is the line of the X-ray, I 0 is the initial intensity of the X-ray beam when entering , and I 1 is the attenuated intensity at the detector. Below, we present two usual ways to organize and interpret collections of measured line integrals (6) in the following form (1): • The Radon transform, • Pencil beam model.

Pencil beam model
Suppose we have N 1 projection images with a digital detector consisting of N 2 pixels. Then our data consist of integral of v over N = N 1 N 2 different lines L in (6). Accordingly, the operator in (1) is defined as the measurement is a vector m ∈ R N , and noise is modeled by a Gaussian zero-centered random vector ε taking values in R N .

Classical level set method for inverse obstacle problems
Consider a physical parameter of the form σ = σ 0 + cχ 1 , where σ 0 (x) is known background, c is a constant, and the characteristic function χ 1 (x) is discontinuous at the boundary of 1 . In inverse obstacle problems one aims to recover the set 1 from indirect measurement of σ. For example, the parameter σ may be sound speed or electrical impedance, and one may measure scattered waves or voltage-to-current boundary maps, respectively.
In the classical level set approach the obstacle is represented as H ( ), where H is the Heaviside function and is smooth. The boundary ∂ 1 of the obstacle is given by the zero level set of . The measurement is written in the form m = A(H ( )) =: Q( ).
In the classical level set method the function is found as the limit where φ is the solution of the evolution equation Here θ is a nonnegative function, (D Q| φ )ρ is the Gateaux derivative of Q at the point φ in direction ρ ∈ C ∞ 0 ( ), and (D Q| φ ) * is the adjoint operator. The intuition behind this approach is the following: define a cost functional where ., . is the inner product of L 2 (D). If we compute, Now, since (D Q) * is the adjoint operator, we obtain that Consequently, if φ is a solution of Eq. (7), we derive formally The above inequality means that the function F 0 (φ(x, t)) is decreasing with respect to t and therefore there exists F 0 ( ) = lim t→∞ F 0 (φ(x, t)) and moreover F 0 ( ) ≤ F 0 (φ). [19] Consider the measurement given by Eq. (1) in the case that the X-ray attenuation coefficient v is smooth and nonnull only inside a subset 1 ⊂ . Now, the operator A * A arising from Radon transform, or pencil beam model, is nonlocal, and mathematical justification of the classical level set approach described in such section does not seem easy.

Motivation for the method developed in
Therefore, we design an algorithm that (i) constructs an approximation 2 for the subset 1 , and (ii) with given approximation 2 produces a reconstruction w that solves the Tikhonov regularization problem where β > 0 is a parameter and the minimum is taken over all u satisfying

Formulation of this method
We approximate the X-ray attenuation coefficient v by w = f ( ), where f is given by (3) and is smooth. Notice that (a) is achieved naturally with ∂ 2 given by the zero level of . The measurement of X-ray projection images is now modeled by m = A( f ( )).
In this method the function is found as the limit ( with ν the interior normal of ∂ , β > 0 a regularization parameter, and r ≥ 0. Compare Eq. (9) to (7) with the choice θ ≡ 1.
The function w in requirements (a) and (b) satisfies The solution of Eq. (9) converges to the solution of the above equation and simultaneously produces a useful approximation 2 for 1 .
How did they come up with such a formulation? Tikhonov regularization yields the cost functional Computing the derivative ∂ t F(u) in a similar way to the process for F 0 , we obtain By Green's formula and using the definition of φ t , we have Therefore, which suggests the evolution equation However, the last equation is numerically unstable. Outside the level set 2 (t) := {x : φ(x, t) = 0} the evolution is driven by the term −β φ alone, pushing towards constant value zero in \ 2 .
Thus we remove the Heaviside function in the last equation and arrive at Eq. (9). Numerical tests show that such an equation is stable and gives much better reconstructions than the last one.

Preliminaries
A mapping A : X → 2 X will be called an operator on X. The domain of A is denoted by D(A) and its range by R( A). Sometimes, we will identify an operator by its graph and we will write (x, y) ∈ A instead of x ∈ D(A) and y ∈ Ax. An operator A on X is said to be accretive if the inequality If, in addition, R(I + λA) is for one, hence for all, λ > 0, precisely X , then A is called m-accretive. Accretive operators were introduced by Browder [6] and Kato [17] independently.
Those accretive operators which are m-accretive play an important role in the study of nonlinear partial differential equations.
Consider the Cauchy problem where A is m-accretive on X and f ∈ L 1 (0, T, X ). It is well known that (13) has a unique integral solution in the sense of Bénilan [4], i.e., there exists a unique continuous function u : [0, T ] → D(A) such that u(0) = x 0 , and moreover, for each (x, y) ∈ A and 0 ≤ s ≤ t ≤ T, we have Here the function ·, · + : We now recall some important facts regarding accretive operators which will be used in our paper (see for example [10]). Proposition 3.1 Let A : D(A) → 2 X be an operator on X. The following conditions are equivalent: A strong solution of Problem (13) is a function u ∈ W 1,∞ (0, T ; X ), i.e., u is locally absolutely continuous and almost differentiable everywhere, u ∈ L ∞ (0, T ; X ), and u (t) Concerning the existence of strong solutions, the following theorem is known (see page 133 of [5]): We refer to [9], where the reader will find a deep study on the Radon-Nikodym property and in particular the reader may pay attention to two very important classes of Banach spaces that enjoy this property: reflexive Banach spaces and Banach spaces whose duals are separable.
On the other hand, we say that u ∈ C(0, T ; X ) is a weak solution of Problem (13)  With respect to the existence of weak solutions, the following result, which is an easy consequence of Theorem 3.2 (see page 134 of [5] for the case of reflexive Banach spaces), is important: Theorem 3.3 Let X be a Banach space with the Radon-Nikodym property. Then Problem (13) admits a unique weak solution which is the unique integral solution of this problem.

General theory
has a unique integral solution.
Proof Let T be a positve number. Consider the subset K := {u ∈ C(0, T ; X ) : u(0) = x 0 }. Given v ∈ K, the problem has a unique integral solution, namely S(v) ∈ K. This fact allows us to introduce a mapping S : K → K defined as follows: given v ∈ K, S(v) is the unique integral solution of the aforementioned problem.

From the definition of an integral solution, it follows that
Using an inductive process, we deduce that, for each n ∈ N, This means that there exists n 0 ∈ N such that S n 0 is a strict contraction on K. Since K is closed, S has a unique fixed point in K by Banach's fixed point theorem. This fixed point is the unique integral solution of Eq. (16).
Finally, we can define, given t > 0, u(t) := u T (t), where u T is the unique integral solution of Problem (16) with T > t. It is clear that u is the unique integral solution of Problem (15).

Remark 3.5
This type of results has been studied for instance in [14].

Theorem 3.6 Let E be a Banach space with Radon-Nikodym property (RN for short). Under the assumptions of Theorem 3.4, if we define the Cauchy problem
then it has a unique weak solution.
Proof By Theorem 3.4, we know that Problem (17) has a unique integral solution, say w. This means that w is the unique solution of the problem Since B(w(.)) ∈ L 1 (0, T ; E) for every T > 0, by Theorem 3.3, w has to be a weak solution of the Problem.

Remark 3.9
The main result of [11] and remark 3.8 of [16] establish that if X is a Banach space and A : D(A) → 2 X is an m-accretive and φ-expansive operator, then A is surjective.

Definition 3.10
Let E be a Banach space, let φ : E → [0, ∞) be a continuous function such that φ(0) = 0, φ(x) > 0, for x = 0, and which satisfies the following condition: For every sequence (x n ) in E such that ( x n ) is decreasing and φ(x n ) → 0 as n → ∞, then x n → 0. An accretive operator A : D(A) → 2 E , with 0 ∈ Az is said to be φ-accretive at zero whenever the inequality holds.

Remark 3.11
The uniqueness of a zero for an operator either φ-expansive or φ-accretive at zero is an immediate consequence of (19) or (20), respectively. On the other hand, Proposition 3.4 and Remark 4.5 of [13] prove that every m-ψ-strongly accretive operator is both ψ-expansive and φ-accretive at zero with φ = ψ • · . Finally, it is proved in [13] that there is not any relationship between being φ-expansive and being φ-accretive at zero.

Proposition 3.12 If P is m-ψ-strongly accretive and h ∈ E, then the operator
Proof Since P is ψ-strongly accretive, then by the above remark we know that P is φ-expansive. Since P is also m-accretive then, by Remark 3.9, P is surjective. Then, there exists z ∈ D(P) such that 0 = H(z).
For each x ∈ D(P) let u = H(x) we have the following: Consequently, Finally, if we define φ : E → [0, ∞) by φ(u) := ψ( u ) u , we conclude that H is an operator φ-accretive at zero.

Theorem 3.13 Let E be a Banach space with the RN property. Consider P : D(P) ⊆ E → E an m-ψstrongly accretive operator on E. Assume that u 0 is an element of D(P), and h ∈ E. If u : [0, ∞) → D(P) is the unique weak solution of the Cauchy problem
where H = P − h, then lim t→+∞ u(t) = z, z being the unique element in D(P) such that h = P(z).
Proof We know that H is φ-accretive at zero and 0 = H(z). First, consider that the initial data u 0 ∈ D(P). In this case by Theorem 3.
where lim h→0 ξ(t, h) = 0. Since j (t) = u(t) − z , elementary calculus yields On the other hand, since the mapping t → u(t) − z is lipschitzian, it is also differentiable almost everywhere. Consequently, Moreover, since t → u(t) − z is decreasing, the function t → 1 2 d dt u(t) − z 2 is Lebesgue integrable on [0, ∞). Hence by (22) we know that the function t → −u (t), j (t) is also Lebesgue integrable on [0, ∞). Then lim inf t→∞ −u (t), j (t) = 0, which means that there exists a sequence (t n ) with t n → ∞ such that Since H is φ-accretive at zero, we know that and, since the sequence ( u(t n ) − z ) is decreasing, by (23) we derive lim n→∞ u(t n ) − z = 0.
Finally, since the function t → u(t) − z is decreasing, we conclude that Finally, if we suppose that u 0 ∈ D(P), then there exists a sequence (x n ) ⊆ D(P) such that x n → u 0 . If we call u n (t) = S(t)(x n ), by the above argument we have that lim t→∞ u n (t) = z. Now, let us see that lim t→∞ u(t) − z = 0. Indeed, since S(t) is a nonexpansive mapping we have Remark 3.14 The above result is a particular case of Corollary 9 of [12].

Existence of solution for the method
In this section we will see, using the results of the above section, that Eq. (9) has a strong solution in L 2 ( ).
We remark that similar analysis fails for the classical level set approach because using H instead of f leads to a heat equation with very singular source terms.
Consider the initial boundary value Problem (9). If we denote g = A * (m), such Problem is ⎧ ⎨ It is well known, see for instance [3], that if β > 0, and we consider the function j : R → R given by j (s) = r 2 s 2 and define the function φ : L 2 ( ) →] − ∞, +∞] by then ϕ is a proper lower semi-continuous convex function in L 2 ( ) such that D(ϕ) = W 1,2 ( ) and moreover, its subdifferential is given by where D(∂ϕ) = {u ∈ W 2,2 ( ) : ∂ ∂ν u − ru = 0, a.e. on ∂ }. In this case, it is also well known that D(∂ϕ) ⊆ D(ϕ) and D(∂ϕ) We interpret and rewrite Problem (24) as follows: Next we will study the existence of a strong solution of Problem (27) when φ 0 ∈ L 2 ( ).  ( f (u))) + g, then Problem (27) has a unique strong solution whenever φ 0 ∈ L 2 ( ). Proof It is clear that A * A : L 2 ( ) → L 2 ( ) is a continuous linear operator; hence we choose the value of μ such that Since for all s, r ∈ R we have that | f (s) − f (r )| ≤ |s − r |; then it is clear that B : On the other hand, by the above comments, the operator −β is m-accretive on L 2 ( ) when its domain is given by {u ∈ W 2,2 ( ) : ∂ ν u = ru a. e. on ∂ }.
Under these conditions Theorem 3.6 allows us to conclude that Problem (27) has a unique weak solution.
Let us see that such weak solution is in fact a strong solution. Indeed, let w be the solution of the problem, in this case we can consider the function B(w(·)), since B is k-lipschitzian and B(0) = 0, it is clear that B(w(·)) ∈ L 2 (0, T ; L 2 ( )) for all T > 0.
On the other hand, by the above argument, we have that A = ∂ϕ, and thus Theorem 3.6 of [5] yields that w is a strong solution in [0, T ] for all T > 0.

Remark 3.16
The assumptions of Theorem 3.15 are satisfied when A is the Radon transform (including limited angle and local tomography cases). Moreover, R N can be identified with the space L 2 (D) when D = {1, 2, . . . , N } equipped with the counting measure. Thus Theorem 3.15 also covers the pencil beam model. Moreover, this theorem proves the existence of a strong solution for every initial data in L 2 ( ) which allows us to improve Theorem 4.1 of [19], since in this theorem we can only assume the existence of a strong solution when the initial data belongs to W 1,2 ( ).
Suppose that φ(x, t) is the unique strong solution of Problem (9) with φ 0 ∈ L 2 ( ); next we are going to study under what conditions the limit lim t→∞ φ(x, t) exists for r > 0 and β > 0.
The following result tells us that operator −β + A * A f is m-accretive on L 2 ( ) for big enough β :

Theorem 3.17
Let an open bounded subset of R 2 with its boundary ∂ smooth. The operator B := −β + A * A f is m-accretive whenever β is large enough.
Proof First let us see that B is an accretive operator. Indeed, it is well known that under this boundary condition satisfies where λ 0 > 0 is the smallest eigenvalue of − in D(A). Thus, since we can assume that A * A : L 2 ( ) → L 2 ( ) is a continuous linear operator, we choose the value of β such that If we define B(u) = A * (A( f (u))), clearly B is a k-lipschitzian mapping with k = A * A < βλ 0 . Let u, v be two elements of D(A); then we have This proves that B is φ-strongly accretive with φ(t) = (βλ 0 − k)t. In order to show that B is m-accretive, now we are going to prove that −β : D(A) → L 2 ( ) is bijective. It is well known that −β : D(A) → L 2 ( ) is m-accretive; moreover since such operator is linear, by Inequality (29), we derive that this operator is also ψ-expansive with ψ(t) = βλ 0 t. Thus by Theorem 8 of [11] we obtain that the operator −β : D(A) → L 2 ( ) is bijective.
Let Q : L 2 ( ) → D(A) be the inverse operator of −β : D(A) → L 2 ( ); from inequality (29) it is clear that Q is continuous.
To prove that B is m-accretive, we have to see that given h ∈ L 2 ( ) there exists u ∈ D(A) such that u + B(u) = h. This means that we have to solve the equation In order to find a solution of Problem (31) it will be enough to show that the operator K : L 2 ( ) → L 2 ( ) defined by K (u) = −Q(u) + B(Q(u)) − h has a fixed point. Since, in this case, if u is a fixed point of K , then v = Q(u) is a solution of Problem 31.
In order to obtain the result it will be enough to see that there exists R > 0 such that K (S R ) ⊆ B R . It is not difficult to see that If we assume that u = R, since taking β such that 1+ A * A βλ 0 < 1 we achieve the result.

Remark 3.18
In [16] one may find results on the existence of fixed points for the sum of two nonlinear operators which are used to study similar problems to (31). Proof Under the hypothesis of this theorem, Theorem 3.17 proves that the operator B := −β + A * A f is m-accretive. Further, if u, v ∈ D(B) since − (u), u ≥ λ 0 u 2 2 , we have which means that B is φ-expansive and m-accretive, and therefore it is surjective. Then there exists ∈ D(B) such that On the other hand, if βλ 0 > A * A + 1 it is clear that B is ψ-strongly accretive (see the proof of Theorem 3.17) and thus by Proposition 3.12 the operator H = B − g is m-φ-accretive at zero. Moreover, is the unique element of D(B) such that H( ) = 0.
By Theorem 3.15 we know that Problem (9) admits a strong solution whenever the initial data belongs to L 2 ( ). Nevertheless, when βλ 0 > A * A + 1, it is not difficult to see that such strong solution is given by Finally, if we apply Theorem 3.13 we derive that u(t) → as t → ∞.