# Tomographic Inverse Problem with Estimating Missing Projections

Authors: Masashi Kimura; Yusaku Yamaguchi; Omar M. Abou Al-Ola; Tetsuya Yoshinaga
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

Image reconstruction in computed tomography can be treated as an inverse problem, namely, obtaining pixel values of a tomographic image from measured projections. However, a seriously degraded image with artifacts is produced when a certain part of the projections is inaccurate or missing. A novel method for simultaneously obtaining a reconstructed image and an estimated projection by solving an initial-value problem of differential equations is proposed. A system of differential equations is constructed on the basis of optimizing a cost function of unknown variables for an image and a projection. Three systems described by nonlinear differential equations are constructed, and the stability of a set of equilibria corresponding to an optimized solution for each system is proved by using the Lyapunov stability theorem. To validate the theoretical result given by the proposed method, metal artifact reduction was numerically performed.


---

## Body

## 1. Introduction

In the field of computed tomography (CT) [1–4], image reconstruction is generally considered an inverse problem, namely, obtaining pixel values of an image from measured projections and a known projection operator. In X-ray CT, however, an image is seriously degraded by, for example, metal and ring artifacts [5–7] when a certain part of the projections received by the detector is inaccurate or missing. Interpolation projection is a well-known method for reducing such artifacts [5, 7]; that is, inaccurate projections are interpolated from their neighboring data and replaced by synthesized values. Another method for reducing metal artifacts completes inaccurate projections by subtracting the projections that cause the artifacts [6, 7]; in particular, the subtraction process uses synthesized projections reprojected as an image including an estimated and corrected metal portion.As for formulating the inverse problem,p∈R+I and A∈R+I×J denote, respectively, projection data and a projection operator, where R+ represents a set of non-negative real numbers and I and J respectively indicate the numbers of projections and pixels in a reconstructed image. Here, the case in which a part of projection p is given by inaccurate or missing measurements (for example, due to the existence of a metal object in the X-ray CT scan) is considered. It is assumed that, without loss of generality, the elements of p are sorted and divided into two parts as follows:(1)p=qr,where q∈R+M and r∈R+I-M correspond to the inaccurate (I>M>0) and the accurate projection data, respectively. It is also assumed that(2)A=BC,where B∈R+M×J and C∈R+(I-M)×J with rows corresponding to elements of p.To reconstruct an image with pixel values and to estimate certain parts of projections instead of inaccurate partq, unknown variables x∈R+J and y∈R+M are treated, respectively. It is assumed that pixel values correspond to, for example, in the case of X-ray CT, actual values of attenuation coefficients, including the minimum value zero corresponding to attenuation by air. It is also assumed that some elements of vector x satisfying r=Cx are undetermined because rankC<J due to missing measurements. Accordingly, the undetermined elements should be permitted to have arbitrary values. The problem considered in this paper is therefore defined as follows:Definition 1.
The tomographic inverse problem with an inaccurate projection is consistent if the set(3)E=e∈R+J:r=Ceis not empty.It follows that the problem is to obtain unknown variablesx and y by minimizing an appropriate cost function that takes zero at x=e and y=Be for some e∈E, respectively. To solve this problem, the following cost function of x∈R+J and y∈R+M with e∈E is considered:(4)Vx,y=∑j=1Jλj-1KLej,xj+α-1KLBe,y=∑j=1Jλj-1ejlog⁡ejxj+xj-ej+α-1∑i=1MBeilog⁡Beiyi+yi-Beiwhere α>0 is a weight parameter, (v)l denotes the lth element of vector v, KL indicates the generalized Kullback–Leibler divergence [8, 9], and λj is defined as(5)λj=∑i=1IAij-1 for j=1,2,…,J. Note that the function in Eq. (4) is nonnegative for arbitrary x∈R+J and y∈R+M and zero if and only if x=e and y=Be.This paper is organized as follows. In Section2, a novel method of simultaneously obtaining a reconstructed image and an estimated projection is proposed. It solves an initial-value problem of differential equations consisting of state variables including parts of projections as well as image pixel values. A system of differential equations based on an extension of dynamical methods [10–20] is constructed. It is an approach that optimizes a cost function of unknown variables consisting of an image and a projection. Three systems described by nonlinear differential equations with different vector fields among each other are presented. In Section 3, for each system, the stability of a set of equilibria corresponding to an optimized solution is proved by using the Lyapunov stability theorem [21]. In Section 4, the intention of the proposed dynamical system is described in detail, and how it works is shown by using a simple toy model. Additionally, to validate the theoretical results given by the proposed method in comparison with results obtained by linear interpolation and a reduced system, reduction of metal artifacts is numerically simulated by using a large-sized model.


## 2. Dynamical Systems

Three kinds of continuous-time dynamical systems with state variables(x⊤,y⊤)⊤∈R++K are proposed, where K≔J+M and R++ denotes the set of positive real numbers. In these systems, V(x(t),y(t)) is expected to decrease in time t along their solutions. The three dynamical systems are described by autonomous nonlinear differential equations. The first system is defined as(6)dxdt=XΛB⊤y-Bx+XΛC⊤r-Cx,dydt=αYBx-y,with initial states x(0)=x0∈R++J and y(0)=y0∈R++M, where X, Y, and Λ indicate the diagonal matrices of vectors x, y, and λ consisting of λj in Eq. (5), respectively. Note that variables x(t) and y(t) are mutually coupled in the vector field. The derivative dx/dt for the dynamics of image reconstruction in the case M=0 has the same form as the continuous-time image reconstruction system presented in Ref. [17].The second and third systems, which are respectively inspired by continuous analogs [19, 22] of the multiplicative algebraic reconstruction technique [23, 24] and the maximum-likelihood expectation-maximization [3, 25] algorithm, are given by(7)dxdt=XΛB⊤Logy-LogBx+XΛC⊤Logr-LogCx,dydt=αYLogBx-Logy,and(8)dxdt=XLogΛB⊤Exp⁡Logy-LogBx+ΛC⊤Exp⁡Logr-LogCx,dydt=αYLogBx-Logy,with the same initial states as in Eq. (6). Here, Log(β) and Exp(β) denote vector-valued functions Log(β)≔(log(β1),log(β2),…,log(βL))⊤ and Exp(β)≔(exp(β1),exp(β2),…,exp(βL))⊤ of each element in vector β=(β1,β2,…,βL)⊤, respectively.


## 3. Theoretical Analysis

Theoretical results for the common behavior of the solutions to the three dynamical systems are presented hereafter. It is first shown that any solution to each dynamical system is positive regardless of whether the inverse problem is consistent or not.Proposition 2.
If initial value(x0⊤,y0⊤)⊤∈R++K in each of the dynamical systems described by Eqs. (6), (7), and (8) is taken, solution φ(t,x0,y0) stays in R++K for all t>0.Proof.
Since the dynamics of thekth elements of w≔(x⊤,y⊤)⊤ can be written as dwk/dt=wkfk(w) with a function fk, it follows that, on the subspace restricted to wk=0, the solution satisfies dφk/dt≡0 for any k=1,2,…,K. Therefore, according to the uniqueness of solutions for the initial value problem [26], the subspace is invariant and trajectories cannot pass through every invariant subspace. This restriction leads to any solution φ(t,x0,y0) of any system with initial value (x0⊤,y0⊤)⊤∈R++K at t=0 being in R++K for all t>0.Under the assumption that the tomographic inverse problem minimizing cost functionV in Eq. (4) is consistent, equilibrium (e⊤,(Be)⊤)⊤∈R++K for each of the differential equations in Eqs. (6), (7), and (8) is asymptotically stable according to the Lyapunov theorem [21].Theorem 3.
If the setE in Eq. (3) is not empty, the solution (x⊤,y⊤)⊤ to the dynamical system described by Eq. (6) with positive initial values asymptotically converges to (e⊤,(Be)⊤)⊤ for some e∈E. The same property holds for the dynamical systems in Eqs. (7) and (8).Proof.
Consider the functionV(x,y)≥0 in Eq. (4) as a Lyapunov-candidate-function, which is well-defined because, from Proposition 2, the state variable (x(t)⊤,y(t)⊤)⊤ of the dynamical system in Eq. (6) belongs to R++M for any t by choosing positive initial values. The function can be rewritten as (9)Vx,y=∑j=1Jλj-1∫ejxju-ejudu+α-1∑i=1M∫Beiyiv-Beivdv,so the derivative along the solution to the system in Eq. (6) is given by (10)ddtVx,y6=-∑j=1Jxj-ejB⊤jBx-y-∑j=1Jxj-ejC⊤jCx-r-∑i=1Myi-Beiyi-Bxi=-Bx-Be⊤Bx-y-Cx-r⊤Cx-r-Be-y⊤Bx-y=-Cx-r22+Bx-y22≤0,with equality if and only if x=e and y=Be. Thus, the asymptotic convergence of solutions to some e in E is guaranteed by using the Lyapunov theorem.
In the same way, we see thatV is also Lyapunov functions for the systems described by Eqs. (7) and (8) according to (11)ddtVx,y7=-∑j=1Jxj-ejB⊤jLogBx-Logy-∑j=1Jxj-ejC⊤jLogCx-Logr-∑i=1Myi-Beilogyi-logBxi=-Bx-Be⊤LogBx-Logy-Cx-r⊤LogCx-Logr-Be-y⊤LogBx-Logy=-KLCx,r+KLr,Cx+KLBx,y+KLy,Bx≤0,and (12)ddtVx,y8=-∑j=1Jej-xjλj-1log⁡λjB⊤jExp⁡Log⁡y-LogBx+λjC⊤jExp⁡Logr-LogCx-∑i=1Myi-Beilogyi-logBxi≤-∑j=1JejB⊤jLogy-LogBx-∑j=1JejC⊤jLogr-LogCx+∑j=1JxjB⊤jExpLogy-LogBx-u+∑j=1JxjC⊤jExp⁡Log⁡r-Log⁡Cx-v-∑i=1Myilog⁡yi-log⁡Bxi+∑i=1MBeilogyi-logBxi=-∑j=1JejC⊤jLogr-LogCx+∑j=1JxjC⊤jExpLogr-LogCx-v-∑i=1Myilogyi-logBxi+∑i=1MBxiexplogyi-logBxi-1≤-∑i=1I-Mri1-explogCxi-logri+∑i=1I-MCxiexplogri-logCxi-1-∑i=1Myi1-explogBxi-logyi+∑i=1MBxiexplogyi-logBxi-1=0,respectively, where u=(1,1,…,1)⊤∈R++M and v=(1,1,…,1)⊤∈R++I-M. The derivatives are zero if and only if x=e and y=Be. Therefore, the solutions to the dynamical systems in Eqs. (7) and (8) also asymptotically converge to some e∈E.


## 4. Application

Numerical experiments were performed to illustrate the efficiency of the proposed method in applications for reducing metal artifacts.


### 4.1. Experimental Method

The systems of differential equations described by Eqs. (6), (7), and (8) with initial states x(0)=x0∈R++J and y(0)=y0∈R++M were used for the numerical experiments. In the experiments, sinograms consisting of projections described as(13)p∗=q∗r,which resulted in the effects of beam hardening and sever attenuation of the X-ray beam caused by the metal hardware [27–29], were synthesized. Here, q∗∈R+M was obtained by the following procedure. Let inaccurate projections be the outer neighborhood of the metal affected projections that are identified by thresholding [7, 28, 30]. The inaccurate measured projections, q, in Eq. (1) were linearly interpolated from their neighboring projections in r at each common projection angle and replaced by synthesized values, q∗. Because the main feature of the proposed method is that inaccurate projections are considered as variables of optimization, the synthesized sinogram is treated as an initial value for the proposed dynamical system and as a fixed projection value for the system to be compared. Then, initial state (x0⊤,y0⊤)⊤ of the proposed system is chosen as xj0=x∗0>0 for j=1,2,…,J, and y0=q∗. However, the system described by the following differential equations is defined as a linear-interpolated system compared to the proposed system in Eq. (6):(14)dxdt=XΛA⊤p∗-Ax=XΛB⊤q∗-Bx+XΛC⊤r-Cxwith initial state x(0)=x0. Note that both solutions to the differential equation in Eq. (14) with x(0)=x0 and to the equations in Eq (6), in which the time derivative of y is replaced with dy/dt=0, with x(0)=x0 and y(0)=q∗ are equivalent because y(t)≡q∗, for all t≥0, satisfies in regard to the latter system. An ODE solver ode113 in MATLAB (MathWorks, Natick, USA) was used for solving the initial-value problem of the differential equations in Eqs. (6) and (14).For comparison with the proposed system, another system, called a reduced system, is established. Namely, for reconstructing imagez∈R++J subject to minimizing r-Cz2, an iterative formula by a projected Landweber algorithm [31] is defined as(15)zn+1=gzn+,n=0,1,2,…with z(0)=x0, where(16)gzn≔zn+ρ-1C⊤r-Cznwith ρ denoting the largest eigenvalue of matrix C⊤C and vector-valued function (β)+ of vector β=(β1,β2,…,βL)⊤ indicating (β)+≔(max{0,β1},max{0,β2},…,max{0,βL})⊤. Relaxation parameter ρ-1 is set according to Theorem 3.1 in Ref. [32] so that the iterative sequence converges to a solution closest to e∈E.To evaluate the quality of the images quantitatively, a distance measure is defined by theL1-norm as(17)Ux=∑j=1j∉MJxj-ej∗where M⊂{1,2,…,J} denotes the set of indices corresponding to the metal position in the true image e∗∈R+J.


### 4.2. Simple Model

To explain how the proposed system works, a simple toy model is considered as follows. In the model, true imagee∗∈R+9 is defined as(18)e∗=0.9,1,0,0,0,0,0.7,0,0⊤.A graphical image of 3×3 pixels corresponding to e∗ is shown in Figure 1(a). It is assumed that a metal object in the phantom is located at the two pixels colored in white. Although the pixel values are undetermined, due to the high attenuation of metal in X-ray CT-scan imaging, for convenience, the values of e5∗ and e6∗ were both set to zero. A sinogram obtained by taking samples at every 120 degrees in the range of 360 and 7 detector bins per projection angle is shown in Figure 1(b). The white regions in the 3×7-sized sinogram indicate an inaccurate part q∈R+8. The correspondence relation between elements p and (q⊤,r⊤)⊤ in Eq. (1) is described by the following format, which corresponds to the sinogram in Figure 1(b).(19)p1p2p3p4p5p6p7p8p9p10p11p12p13p14p15p16p17p18p19p20p21=r1r2q1q2r3r4r5r6q3q4q5r7r8r9r10r11q6q7q8r12r13≕Sq1Sq2Sq3=Sq∈R+3×7.When using this notation, in the case of proposed system, the elements of measured projection r and estimated projection y(t), t≥0, instead of q, can be arranged as matrix S(y(t)) in the same format.Figure 1
(a) Phantom image and (b) sinogram.
(a) (b)Images were reconstructed by using solutions of three systems described by Eqs. (6), (14), and (15). The parameter α used for the proposed system, uniform initial pixel values x∗0, integration time τ, and iteration number ν were set to 0.1, 0.5, 105, and 105, respectively, unless otherwise specified. An example of reconstructed images is shown in Figure 2. Values of distance measures U(x(τ)) and U(z(ν)) defined in Eq. (17) with M={5,6} for reconstructed images (a), (b), and (c) in Figure 2 are 0.2887, 1.3773, and 2.4881, respectively. It can be seen that the proposed reconstruction has better image quality. Moreover, as indicted in Table 1, compared to the reduced method, the proposed method gives more robust performance in terms of the variation of initial values.Table 1
U ( x ( τ ) ) and U(z(ν)) for images reconstructed using the proposed method with α=0.01 and 0.1, linear-interpolated method, and reduced method with variation of initial pixel values.
x∗0 proposed method interpolated reduced α = 0.01 α = 0.1 method method 0.1 0.2255 0.2300 1.3773 0.7227 0.5 0.2308 0.2887 1.3773 2.4881 1 0.2331 0.3200 1.3773 5.0240Figure 2
Reconstructed images using (a) proposed, (b) linear-interpolated, and (c) reduced systems.
(a) (b) (c)For each of the proposed and reduced systems, it was numerically confirmed that the distance betweenr and Cx(τ) (or Cz(ν)) is zero; that is, convergent values x(τ) and z(ν) belong to set E in Eq. (3). An explicit description of this set is introduced as follows. Neglecting any trivial element ri of r∈R+13 and the corresponding row Ci of C∈R+13×9 such that ri=0 and Cij=0 for any j (excluding the case that ri=0, but Cij≠0 for some j) obtains(20)r2r6r7r12=0.54860.141300.2573 and(21)C∗≔C2C6C7C12=0.3429000.3429000.3429000.10660.04530.0075000000000000000.00430.0573000.1887000.29400.01310.The set defined in Eq. (3) can therefore be explicitly described as(22)E=e∗+Ns:s∈R5where N is the matrix whose column vectors are the orthonormal basis for the null space of matrix C∗, i.e., (23)N=00-0.2766-0.00110.0382000.65620.00390.074000-0.0302-0.0080-0.9918000.60680.1247-0.0849100000100000-0.3302-0.12360.046700-0.12150.98440.008200000.Approximate values of coefficients s in Eq. (22) used for the proposed and reduced systems are respectively given as(24)0.00440.1407-0.01630.0014-0.1194and0.50000.50000.07610.5737-0.5958.It is clear that the convergent value obtained by the proposed system gives smaller absolute values of s3, s4, and s5.A graph of the distance measuresU(x(τ)) and U(z(ν)) for comparing the three methods is shown in Figure 3. The distance given by the proposed method varies with different values of α. The figure shows that there exists a range of α at least in interval (0,0.1] in which the quantitative measure has a minimum value. Note that the system of Eq. (6) as α tends to zero is similar to the interpolated system of Eq. (14). A smaller value of α results in not only a larger weight of term KL(Be,y) in Eq. (4) but also a slower dynamics of state variable y(t). On the other hand, when parameter α tends to infinity, the system of Eq. (6) is equivalent to a reduced system optimizing the function ∑j=1Jλj-1KL(ej,xj), which is a part of the cost function V in Eq. (4), without information about q∗. The dynamics of solutions to the proposed system are explained as follows. In Figure 4, the dynamics expressed by the time evolution of state variables obtained by the three dynamical systems is illustrated. With regard to the proposed system, the behavior of Bx(t) is restricted by the dynamics of y(t) emanating from interpolated value q∗ and converges to coincide with y(t) after the passage of time. However, in the interpolated and reduced systems, Bx(t) attempts to reach q∗ without reaching it and is independent of q∗, respectively. The mutual dynamical coupling of Bx(t) and y(t) in the proposed system gives rise to a better convergence performance of the state x(t) in regard to approaching ideal value e∗. The dynamics of y(t) that eventually converges to a limit set in the neighborhood of q∗ in the state space is effective for estimating the inaccurate projection and produces a better quality image through Bx(t), which converges to y(t).Figure 3
Distances (atτ=105 and ν=105) for the proposed method with variation of α, interpolated method, and reduced method, which is independent of α, indicated by blue points, a green point, and a magenta line, respectively.Figure 4
Time evolution of state variable for pixel valuesx (upper panel), and projection values (lower panel) obtained by using (a) proposed, (b) linear-interpolated, and (c) reduced systems. In the upper graph, trajectories xj and zj with lines and ej∗ marked with colored dots at the right edge are indicated by the same color as the corresponding pixel in the phantom image. In the lower graph, blue and magenta lines denote y (or q∗) and Bx (or Bz), respectively.
(a) (b) (c)


### 4.3. Large-Sized Model

The proposed method was experimentally demonstrated by applying it to reducing metal artifacts by using a larger sized model. A phantom image with192×192 pixels, so J=36,864, was simulated as shown in Figure 5(a). An image was reconstructed by using 275 detector bins per projection angle and 360-degree scanning with sampling every two degrees, which corresponds to number of projections I=49,500. For simulating reduction of metal artifacts, it was assumed that three small metal objects exist in the phantom. The white regions in the gray-scaled sinogram shown in Figure 5(b) indicate inaccurate part q in the sinogram or projection p in Eq. (1).Figure 5
(a) Phantom image with three metal positions (small black regions) indicated by white arrows and (b) its sinogram.
(a) (b)The experimental results obtained by using the proposed system described by Eq. (6) with α=0.01, the linear-interpolated system described by Eq. (14), and the reduced system described by Eq. (15) are presented as follows. Reconstructed images obtained from solutions x(τ) at τ=100 for the continuous-time system and z(ν) at ν=3,500 for the discrete-time system, which are emanating from the initial value with x∗0=0.5, are shown in Figure 6. The values of τ and ν for adjusting the difference between the continuous-time and discrete-time systems describing the proposed and reduced methods, respectively, were chosen so that the values of objective functions defined by r-Cx(t)1 and r-Cz(n)1 coincide at 35τ=ν as shown in Figure 7. Under the same conditions for reconstruction, a qualitative reduction of artifacts in the reconstructed image obtained by the proposed method was compared with that by the other two methods. From Figure 8, which compares quantitative measures U(x(τ)) for the proposed method and U(z(ν)) for the reduced method in addition to the observation of the measure 830.6 for the interpolated method, it is clear that the proposed method with α∈[0.002,0.1] has the smallest measure. The advantage of the proposed method is due to the simultaneous dynamical behavior of a pair of variables, x(t) and y(t), which is based on mutual dynamical coupling described by Eq. (6).Figure 6
Reconstructed images by (a) proposed, (b) interpolated, and (c) reduced methods.
(a) (b) (c)Figure 7
Time evolutions of objective function by the proposed, interpolated, and reduced methods indicated by blue, green, and magenta lines, respectively.Figure 8
Distances (atτ=100 and ν=3,500) for the proposed method with variation of α and the reduced method being independent of α indicated by blue points and magenta lines, respectively.Discretization of the differential equations is required for low-cost computation in practical use. Iterative formulae derived from a multiplicative Euler discretization [33] of differential equations in Eqs. (7) and (8) are, respectively, given by (25)zn+1=ZnExp⁡ΛB⊤Logwn-LogBzn+ΛC⊤Logr-LogCzn,wn+1=WnExpαLogBzn-Logwn,and(26)zn+1=ZnΛB⊤Exp⁡Logwn-LogBzn+ΛC⊤Exp⁡Logr-LogCzn,wn+1=WnExp⁡αLogBzn-Logwn,for n=0,1,2,…, with initial states (z(0)⊤,w(0)⊤)⊤=(x0⊤,y0⊤)⊤, where Z and W indicate the diagonal matrices of z and w, respectively. Graphs of the distances measured by U(z) for reconstructed images obtained after 3,500 iterations of the proposed systems in Eqs. (25) and (26) and the reduced system in Eq. (15) are shown in Figure 9. The multiplicative Euler discretization approximates solutions to the differential equations of Eqs. (7) and (8) at the discrete points with a step size equal to one. It is noted that the course of the distances with variation α for the proposed method given by Eqs. (25) and (26) as well as Eq. (6) is qualitatively consistent and the values of these distances, in the range of α, are lower than that given by the reduced method. The system of Eq. (25) yields the best performance with respect to minimizing distance U(z), under the same conditions.Figure 9
Distances (atν=3,500) for proposed method given (a) by Eq. (25) and (b) Eq. (26) with variation of α along with reduced method being independent of α indicated by blue points and magenta lines, respectively.
(a) (b)


## 4.1. Experimental Method

The systems of differential equations described by Eqs. (6), (7), and (8) with initial states x(0)=x0∈R++J and y(0)=y0∈R++M were used for the numerical experiments. In the experiments, sinograms consisting of projections described as(13)p∗=q∗r,which resulted in the effects of beam hardening and sever attenuation of the X-ray beam caused by the metal hardware [27–29], were synthesized. Here, q∗∈R+M was obtained by the following procedure. Let inaccurate projections be the outer neighborhood of the metal affected projections that are identified by thresholding [7, 28, 30]. The inaccurate measured projections, q, in Eq. (1) were linearly interpolated from their neighboring projections in r at each common projection angle and replaced by synthesized values, q∗. Because the main feature of the proposed method is that inaccurate projections are considered as variables of optimization, the synthesized sinogram is treated as an initial value for the proposed dynamical system and as a fixed projection value for the system to be compared. Then, initial state (x0⊤,y0⊤)⊤ of the proposed system is chosen as xj0=x∗0>0 for j=1,2,…,J, and y0=q∗. However, the system described by the following differential equations is defined as a linear-interpolated system compared to the proposed system in Eq. (6):(14)dxdt=XΛA⊤p∗-Ax=XΛB⊤q∗-Bx+XΛC⊤r-Cxwith initial state x(0)=x0. Note that both solutions to the differential equation in Eq. (14) with x(0)=x0 and to the equations in Eq (6), in which the time derivative of y is replaced with dy/dt=0, with x(0)=x0 and y(0)=q∗ are equivalent because y(t)≡q∗, for all t≥0, satisfies in regard to the latter system. An ODE solver ode113 in MATLAB (MathWorks, Natick, USA) was used for solving the initial-value problem of the differential equations in Eqs. (6) and (14).For comparison with the proposed system, another system, called a reduced system, is established. Namely, for reconstructing imagez∈R++J subject to minimizing r-Cz2, an iterative formula by a projected Landweber algorithm [31] is defined as(15)zn+1=gzn+,n=0,1,2,…with z(0)=x0, where(16)gzn≔zn+ρ-1C⊤r-Cznwith ρ denoting the largest eigenvalue of matrix C⊤C and vector-valued function (β)+ of vector β=(β1,β2,…,βL)⊤ indicating (β)+≔(max{0,β1},max{0,β2},…,max{0,βL})⊤. Relaxation parameter ρ-1 is set according to Theorem 3.1 in Ref. [32] so that the iterative sequence converges to a solution closest to e∈E.To evaluate the quality of the images quantitatively, a distance measure is defined by theL1-norm as(17)Ux=∑j=1j∉MJxj-ej∗where M⊂{1,2,…,J} denotes the set of indices corresponding to the metal position in the true image e∗∈R+J.


## 4.2. Simple Model

To explain how the proposed system works, a simple toy model is considered as follows. In the model, true imagee∗∈R+9 is defined as(18)e∗=0.9,1,0,0,0,0,0.7,0,0⊤.A graphical image of 3×3 pixels corresponding to e∗ is shown in Figure 1(a). It is assumed that a metal object in the phantom is located at the two pixels colored in white. Although the pixel values are undetermined, due to the high attenuation of metal in X-ray CT-scan imaging, for convenience, the values of e5∗ and e6∗ were both set to zero. A sinogram obtained by taking samples at every 120 degrees in the range of 360 and 7 detector bins per projection angle is shown in Figure 1(b). The white regions in the 3×7-sized sinogram indicate an inaccurate part q∈R+8. The correspondence relation between elements p and (q⊤,r⊤)⊤ in Eq. (1) is described by the following format, which corresponds to the sinogram in Figure 1(b).(19)p1p2p3p4p5p6p7p8p9p10p11p12p13p14p15p16p17p18p19p20p21=r1r2q1q2r3r4r5r6q3q4q5r7r8r9r10r11q6q7q8r12r13≕Sq1Sq2Sq3=Sq∈R+3×7.When using this notation, in the case of proposed system, the elements of measured projection r and estimated projection y(t), t≥0, instead of q, can be arranged as matrix S(y(t)) in the same format.Figure 1
(a) Phantom image and (b) sinogram.
(a) (b)Images were reconstructed by using solutions of three systems described by Eqs. (6), (14), and (15). The parameter α used for the proposed system, uniform initial pixel values x∗0, integration time τ, and iteration number ν were set to 0.1, 0.5, 105, and 105, respectively, unless otherwise specified. An example of reconstructed images is shown in Figure 2. Values of distance measures U(x(τ)) and U(z(ν)) defined in Eq. (17) with M={5,6} for reconstructed images (a), (b), and (c) in Figure 2 are 0.2887, 1.3773, and 2.4881, respectively. It can be seen that the proposed reconstruction has better image quality. Moreover, as indicted in Table 1, compared to the reduced method, the proposed method gives more robust performance in terms of the variation of initial values.Table 1
U ( x ( τ ) ) and U(z(ν)) for images reconstructed using the proposed method with α=0.01 and 0.1, linear-interpolated method, and reduced method with variation of initial pixel values.
x∗0 proposed method interpolated reduced α = 0.01 α = 0.1 method method 0.1 0.2255 0.2300 1.3773 0.7227 0.5 0.2308 0.2887 1.3773 2.4881 1 0.2331 0.3200 1.3773 5.0240Figure 2
Reconstructed images using (a) proposed, (b) linear-interpolated, and (c) reduced systems.
(a) (b) (c)For each of the proposed and reduced systems, it was numerically confirmed that the distance betweenr and Cx(τ) (or Cz(ν)) is zero; that is, convergent values x(τ) and z(ν) belong to set E in Eq. (3). An explicit description of this set is introduced as follows. Neglecting any trivial element ri of r∈R+13 and the corresponding row Ci of C∈R+13×9 such that ri=0 and Cij=0 for any j (excluding the case that ri=0, but Cij≠0 for some j) obtains(20)r2r6r7r12=0.54860.141300.2573 and(21)C∗≔C2C6C7C12=0.3429000.3429000.3429000.10660.04530.0075000000000000000.00430.0573000.1887000.29400.01310.The set defined in Eq. (3) can therefore be explicitly described as(22)E=e∗+Ns:s∈R5where N is the matrix whose column vectors are the orthonormal basis for the null space of matrix C∗, i.e., (23)N=00-0.2766-0.00110.0382000.65620.00390.074000-0.0302-0.0080-0.9918000.60680.1247-0.0849100000100000-0.3302-0.12360.046700-0.12150.98440.008200000.Approximate values of coefficients s in Eq. (22) used for the proposed and reduced systems are respectively given as(24)0.00440.1407-0.01630.0014-0.1194and0.50000.50000.07610.5737-0.5958.It is clear that the convergent value obtained by the proposed system gives smaller absolute values of s3, s4, and s5.A graph of the distance measuresU(x(τ)) and U(z(ν)) for comparing the three methods is shown in Figure 3. The distance given by the proposed method varies with different values of α. The figure shows that there exists a range of α at least in interval (0,0.1] in which the quantitative measure has a minimum value. Note that the system of Eq. (6) as α tends to zero is similar to the interpolated system of Eq. (14). A smaller value of α results in not only a larger weight of term KL(Be,y) in Eq. (4) but also a slower dynamics of state variable y(t). On the other hand, when parameter α tends to infinity, the system of Eq. (6) is equivalent to a reduced system optimizing the function ∑j=1Jλj-1KL(ej,xj), which is a part of the cost function V in Eq. (4), without information about q∗. The dynamics of solutions to the proposed system are explained as follows. In Figure 4, the dynamics expressed by the time evolution of state variables obtained by the three dynamical systems is illustrated. With regard to the proposed system, the behavior of Bx(t) is restricted by the dynamics of y(t) emanating from interpolated value q∗ and converges to coincide with y(t) after the passage of time. However, in the interpolated and reduced systems, Bx(t) attempts to reach q∗ without reaching it and is independent of q∗, respectively. The mutual dynamical coupling of Bx(t) and y(t) in the proposed system gives rise to a better convergence performance of the state x(t) in regard to approaching ideal value e∗. The dynamics of y(t) that eventually converges to a limit set in the neighborhood of q∗ in the state space is effective for estimating the inaccurate projection and produces a better quality image through Bx(t), which converges to y(t).Figure 3
Distances (atτ=105 and ν=105) for the proposed method with variation of α, interpolated method, and reduced method, which is independent of α, indicated by blue points, a green point, and a magenta line, respectively.Figure 4
Time evolution of state variable for pixel valuesx (upper panel), and projection values (lower panel) obtained by using (a) proposed, (b) linear-interpolated, and (c) reduced systems. In the upper graph, trajectories xj and zj with lines and ej∗ marked with colored dots at the right edge are indicated by the same color as the corresponding pixel in the phantom image. In the lower graph, blue and magenta lines denote y (or q∗) and Bx (or Bz), respectively.
(a) (b) (c)


## 4.3. Large-Sized Model

The proposed method was experimentally demonstrated by applying it to reducing metal artifacts by using a larger sized model. A phantom image with192×192 pixels, so J=36,864, was simulated as shown in Figure 5(a). An image was reconstructed by using 275 detector bins per projection angle and 360-degree scanning with sampling every two degrees, which corresponds to number of projections I=49,500. For simulating reduction of metal artifacts, it was assumed that three small metal objects exist in the phantom. The white regions in the gray-scaled sinogram shown in Figure 5(b) indicate inaccurate part q in the sinogram or projection p in Eq. (1).Figure 5
(a) Phantom image with three metal positions (small black regions) indicated by white arrows and (b) its sinogram.
(a) (b)The experimental results obtained by using the proposed system described by Eq. (6) with α=0.01, the linear-interpolated system described by Eq. (14), and the reduced system described by Eq. (15) are presented as follows. Reconstructed images obtained from solutions x(τ) at τ=100 for the continuous-time system and z(ν) at ν=3,500 for the discrete-time system, which are emanating from the initial value with x∗0=0.5, are shown in Figure 6. The values of τ and ν for adjusting the difference between the continuous-time and discrete-time systems describing the proposed and reduced methods, respectively, were chosen so that the values of objective functions defined by r-Cx(t)1 and r-Cz(n)1 coincide at 35τ=ν as shown in Figure 7. Under the same conditions for reconstruction, a qualitative reduction of artifacts in the reconstructed image obtained by the proposed method was compared with that by the other two methods. From Figure 8, which compares quantitative measures U(x(τ)) for the proposed method and U(z(ν)) for the reduced method in addition to the observation of the measure 830.6 for the interpolated method, it is clear that the proposed method with α∈[0.002,0.1] has the smallest measure. The advantage of the proposed method is due to the simultaneous dynamical behavior of a pair of variables, x(t) and y(t), which is based on mutual dynamical coupling described by Eq. (6).Figure 6
Reconstructed images by (a) proposed, (b) interpolated, and (c) reduced methods.
(a) (b) (c)Figure 7
Time evolutions of objective function by the proposed, interpolated, and reduced methods indicated by blue, green, and magenta lines, respectively.Figure 8
Distances (atτ=100 and ν=3,500) for the proposed method with variation of α and the reduced method being independent of α indicated by blue points and magenta lines, respectively.Discretization of the differential equations is required for low-cost computation in practical use. Iterative formulae derived from a multiplicative Euler discretization [33] of differential equations in Eqs. (7) and (8) are, respectively, given by (25)zn+1=ZnExp⁡ΛB⊤Logwn-LogBzn+ΛC⊤Logr-LogCzn,wn+1=WnExpαLogBzn-Logwn,and(26)zn+1=ZnΛB⊤Exp⁡Logwn-LogBzn+ΛC⊤Exp⁡Logr-LogCzn,wn+1=WnExp⁡αLogBzn-Logwn,for n=0,1,2,…, with initial states (z(0)⊤,w(0)⊤)⊤=(x0⊤,y0⊤)⊤, where Z and W indicate the diagonal matrices of z and w, respectively. Graphs of the distances measured by U(z) for reconstructed images obtained after 3,500 iterations of the proposed systems in Eqs. (25) and (26) and the reduced system in Eq. (15) are shown in Figure 9. The multiplicative Euler discretization approximates solutions to the differential equations of Eqs. (7) and (8) at the discrete points with a step size equal to one. It is noted that the course of the distances with variation α for the proposed method given by Eqs. (25) and (26) as well as Eq. (6) is qualitatively consistent and the values of these distances, in the range of α, are lower than that given by the reduced method. The system of Eq. (25) yields the best performance with respect to minimizing distance U(z), under the same conditions.Figure 9
Distances (atν=3,500) for proposed method given (a) by Eq. (25) and (b) Eq. (26) with variation of α along with reduced method being independent of α indicated by blue points and magenta lines, respectively.
(a) (b)


## 5. Conclusion

An image-reconstruction method for optimizing a cost function of unknown variables corresponding to missing projections as well as image densities was proposed. Three dynamical systems described by nonlinear differential equations were constructed, and the stability of their equilibria was proved by using the Lyapunov theorem. The efficiency of the proposed method was evaluated through numerical experiments simulating the reduction of metal artifacts. Qualitative and quantitative evaluations using image quality and a distance measure indicate that the proposed method performs better than the conventional linear-interpolation and reduced methods. This paper shows a novel approach to optimizing the cost function of unknown variables consisting of both image and projection. However, not only further analysis of convergence rates and investigation on the influence of measured noise, but also experimental evaluation and detailed comparison with other existing methods [7, 29] of artifact reduction are required to extend the proposed approach for practical use.



---

*Source: 7932318-2019-03-13.xml*