Tomographic Inverse Problem with Estimating Missing Projections

Kimura, Masashi; Yamaguchi, Yusaku; Abou Al-Ola, Omar M.; Yoshinaga, Tetsuya

doi:https://doi.org/10.1155/2019/7932318

Mathematical Problems in Engineering

On this page

Abstract Introduction Conclusion Data Availability Conflicts of Interest Acknowledgments References Copyright Related Articles

Research Article | Open Access

Volume 2019 | Article ID 7932318 | https://doi.org/10.1155/2019/7932318

Tomographic Inverse Problem with Estimating Missing Projections

Masashi Kimura,¹Yusaku Yamaguchi,²Omar M. Abou Al-Ola,³and Tetsuya Yoshinaga⁴

Academic Editor: Higinio Ramos

Received17 Dec 2018

Accepted10 Feb 2019

Published13 Mar 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.

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, and denote, respectively, projection data and a projection operator, where represents a set of non-negative real numbers and and respectively indicate the numbers of projections and pixels in a reconstructed image. Here, the case in which a part of projection 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 are sorted and divided into two parts as follows:where and correspond to the inaccurate () and the accurate projection data, respectively. It is also assumed thatwhere and with rows corresponding to elements of .

To reconstruct an image with pixel values and to estimate certain parts of projections instead of inaccurate part , unknown variables and 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 satisfying are undetermined because 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 setis not empty.

It follows that the problem is to obtain unknown variables and by minimizing an appropriate cost function that takes zero at and for some , respectively. To solve this problem, the following cost function of and with is considered:where is a weight parameter, denotes the th element of vector , indicates the generalized Kullback–Leibler divergence [8, 9], and is defined as for . Note that the function in Eq. (4) is nonnegative for arbitrary and and zero if and only if and .

This paper is organized as follows. In Section 2, 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 are proposed, where and denotes the set of positive real numbers. In these systems, is expected to decrease in time along their solutions. The three dynamical systems are described by autonomous nonlinear differential equations. The first system is defined aswith initial states and , where , , and indicate the diagonal matrices of vectors , , and consisting of in Eq. (5), respectively. Note that variables and are mutually coupled in the vector field. The derivative for the dynamics of image reconstruction in the case 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 byandwith the same initial states as in Eq. (6). Here, and denote vector-valued functions and of each element in vector , 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 in each of the dynamical systems described by Eqs. (6), (7), and (8) is taken, solution stays in for all .

Proof. Since the dynamics of the th elements of can be written as with a function , it follows that, on the subspace restricted to , the solution satisfies for any . 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 of any system with initial value at being in for all .

Under the assumption that the tomographic inverse problem minimizing cost function in Eq. (4) is consistent, equilibrium 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 set in Eq. (3) is not empty, the solution to the dynamical system described by Eq. (6) with positive initial values asymptotically converges to for some . The same property holds for the dynamical systems in Eqs. (7) and (8).

Proof. Consider the function in Eq. (4) as a Lyapunov-candidate-function, which is well-defined because, from Proposition 2, the state variable of the dynamical system in Eq. (6) belongs to for any by choosing positive initial values. The function can be rewritten as so the derivative along the solution to the system in Eq. (6) is given by with equality if and only if and . Thus, the asymptotic convergence of solutions to some in is guaranteed by using the Lyapunov theorem.
In the same way, we see that is also Lyapunov functions for the systems described by Eqs. (7) and (8) according to and respectively, where and . The derivatives are zero if and only if and . Therefore, the solutions to the dynamical systems in Eqs. (7) and (8) also asymptotically converge to some .

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 and were used for the numerical experiments. In the experiments, sinograms consisting of projections described aswhich 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, 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, , in Eq. (1) were linearly interpolated from their neighboring projections in at each common projection angle and replaced by synthesized values, . 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 of the proposed system is chosen as for , and . However, the system described by the following differential equations is defined as a linear-interpolated system compared to the proposed system in Eq. (6):with initial state . Note that both solutions to the differential equation in Eq. (14) with and to the equations in Eq (6), in which the time derivative of is replaced with , with and are equivalent because , for all , 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 image subject to minimizing , an iterative formula by a projected Landweber algorithm [31] is defined aswith , wherewith denoting the largest eigenvalue of matrix and vector-valued function of vector indicating . Relaxation parameter is set according to Theorem 3.1 in Ref. [32] so that the iterative sequence converges to a solution closest to .

To evaluate the quality of the images quantitatively, a distance measure is defined by the -norm aswhere denotes the set of indices corresponding to the metal position in the true image .

4.2. Simple Model

To explain how the proposed system works, a simple toy model is considered as follows. In the model, true image is defined asA graphical image of pixels corresponding to 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 and 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 -sized sinogram indicate an inaccurate part . The correspondence relation between elements and in Eq. (1) is described by the following format, which corresponds to the sinogram in Figure 1(b).When using this notation, in the case of proposed system, the elements of measured projection and estimated projection , , instead of , can be arranged as matrix in the same format.

(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 , integration time , and iteration number were set to , , , and , respectively, unless otherwise specified. An example of reconstructed images is shown in Figure 2. Values of distance measures and defined in Eq. (17) with 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.

(a)

(b)

(c)

For each of the proposed and reduced systems, it was numerically confirmed that the distance between and (or ) is zero; that is, convergent values and belong to set in Eq. (3). An explicit description of this set is introduced as follows. Neglecting any trivial element of and the corresponding row of such that and for any (excluding the case that , but for some ) obtains andThe set defined in Eq. (3) can therefore be explicitly described aswhere is the matrix whose column vectors are the orthonormal basis for the null space of matrix , i.e., Approximate values of coefficients in Eq. (22) used for the proposed and reduced systems are respectively given asIt is clear that the convergent value obtained by the proposed system gives smaller absolute values of , , and .

A graph of the distance measures and 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 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 in Eq. (4) but also a slower dynamics of state variable . On the other hand, when parameter tends to infinity, the system of Eq. (6) is equivalent to a reduced system optimizing the function , which is a part of the cost function in Eq. (4), without information about . 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 is restricted by the dynamics of emanating from interpolated value and converges to coincide with after the passage of time. However, in the interpolated and reduced systems, attempts to reach without reaching it and is independent of , respectively. The mutual dynamical coupling of and in the proposed system gives rise to a better convergence performance of the state in regard to approaching ideal value . The dynamics of that eventually converges to a limit set in the neighborhood of in the state space is effective for estimating the inaccurate projection and produces a better quality image through , which converges to .

(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 with pixels, so , 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 . 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 in the sinogram or projection in Eq. (1).

(a)

(b)

The experimental results obtained by using the proposed system described by Eq. (6) with , 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 at for the continuous-time system and at for the discrete-time system, which are emanating from the initial value with , 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 and coincide at 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 for the proposed method and for the reduced method in addition to the observation of the measure for the interpolated method, it is clear that the proposed method with has the smallest measure. The advantage of the proposed method is due to the simultaneous dynamical behavior of a pair of variables, and , which is based on mutual dynamical coupling described by Eq. (6).

(a)

(b)

(c)

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 andfor , with initial states , where and indicate the diagonal matrices of and , respectively. Graphs of the distances measured by for reconstructed images obtained after 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 , under the same conditions.

(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.

Data Availability

All data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

Acknowledgments

This research was partially supported by JSPS KAKENHI Grant Number 18K04169.

References

A. C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging, IEEE Service Center, Piscataway, NJ, USA, 1988.
View at: MathSciNet
H. Stark, Image Recorvery Theorem And Application, Academic Press, Orlando, Fla, USA, 1987.
View at: MathSciNet
L. A. Shepp and Y. Vardi, “Maximum likelihood reconstruction for emission tomography,” IEEE Transactions on Medical Imaging, vol. 1, no. 2, pp. 113–122, 1982.
View at: Publisher Site | Google Scholar
C. L. Byrne, “Iterative image reconstruction algorithms based on cross-entropy minimization,” IEEE Transactions on Image Processing, vol. 2, no. 1, pp. 96–103, 1993.
View at: Publisher Site | Google Scholar
W. A. Kalender, R. Hebel, and J. Ebersberger, “Reduction of CT artifacts caused by metallic implants,” Radiology, vol. 164, no. 2, Article ID 3602406, pp. 576-577, 1987.
View at: Publisher Site | Google Scholar
C. R. Crawford, J. G. Colsher, N. J. Pelc, and A. H. R. Lonn, “High speed reprojection and its applications,” Proceedings of SPIE - The International Society for Optical Engineering, vol. 914, pp. 311–318, 1988.
View at: Google Scholar
L. Gjesteby, B. De Man, Y. Jin et al., “Metal artifact reduction in CT: where are we after four decades?” IEEE Access, vol. 4, pp. 5826–5849, 2016.
View at: Publisher Site | Google Scholar
S. Kullback and R. A. Leibler, “On information and sufficiency,” Annals of Mathematical Statistics, vol. 22, pp. 79–86, 1951.
View at: Publisher Site | Google Scholar | MathSciNet
I. Csiszár, “Why least squares and maximum entropy? an axiomatic approach to inference for linear inverse problems,” The Annals of Statistics, vol. 19, no. 4, pp. 2032–2066, 1991.
View at: Publisher Site | Google Scholar | MathSciNet
M. K. Gavurin, “Nonlinear functional equations and continuous analogues of iteration methods,” Izvestiya Vysshikh Uchebnykh Zavedenii. Matematika, vol. 5, no. 6, pp. 18–31, 1958.
View at: Google Scholar | MathSciNet
J. Schropp, “Using dynamical systems methods to solve minimization problems,” Applied Numerical Mathematics, vol. 18, no. 1, pp. 321–335, 1995.
View at: Publisher Site | Google Scholar | MathSciNet
R. G. Airapetyan, A. G. Ramm, and A. B. Smirnova, “Continuous analog of the Gauss-Newton method,” Mathematical Models and Methods in Applied Sciences, vol. 9, no. 3, pp. 463–474, 1999.
View at: Publisher Site | Google Scholar | MathSciNet
R. Airapetyan, “Continuous Newton method and its modification,” Applicable Analysis: An International Journal, vol. 73, no. 3-4, pp. 463–484, 1999.
View at: Publisher Site | Google Scholar | MathSciNet
R. G. Airapetyan and A. G. Ramm, “Dynamical systems and discrete methods for solving nonlinear ill-posed problems,” in Applied Mathematics Reviews, G. A. Anastassiou, Ed., vol. 1, pp. 491–536, World Scientific Publishing Company, Singapore, 2000.
View at: Publisher Site | Google Scholar | MathSciNet
A. G. Ramm, “Dynamical systems method for solving operator equations,” Communications in Nonlinear Science and Numerical Simulation, vol. 9, no. 4, pp. 383–402, 2004.
View at: Publisher Site | Google Scholar | MathSciNet
K. Fujimoto, O. Abou Al-Ola, and T. Yoshinaga, “Continuous-time image reconstruction using differential equations for computed tomography,” Communications in Nonlinear Science and Numerical Simulation, vol. 15, no. 6, pp. 1648–1654, 2010.
View at: Publisher Site | Google Scholar
O. M. Abou Al-Ola, K. Fujimoto, and T. Yoshinaga, “Common Lyapunov function based on Kullback–Leibler divergence for a switched nonlinear system,” Mathematical Problems in Engineering, vol. 2011, Article ID 723509, 12 pages, 2011.
View at: Publisher Site | Google Scholar | MathSciNet
Y. Yamaguchi, K. Fujimoto, O. M. Abou Al-Ola, and T. Yoshinaga, “Continuous-time image reconstruction for binary tomography,” Communications in Nonlinear Science and Numerical Simulation, vol. 18, no. 8, pp. 2081–2087, 2013.
View at: Publisher Site | Google Scholar | MathSciNet
K. Tateishi, Y. Yamaguchi, O. M. Abou Al-Ola, and T. Yoshinaga, “Continuous analog of accelerated OS-EM algorithm for computed tomography,” Mathematical Problems in Engineering, vol. 2017, Article ID 1564123, 8 pages, 2017.
View at: Publisher Site | Google Scholar | MathSciNet
R. Kasai, Y. Yamaguchi, T. Kojima, and T. Yoshinaga, “Tomographic image reconstruction based on minimization of symmetrized Kullback-Leibler divergence,” Mathematical Problems in Engineering, vol. 2018, Article ID 8973131, 9 pages, 2018.
View at: Google Scholar
A. M. Lyapunov, Stability of Motion, Academic Press, New York, NY, USA, 1966.
View at: MathSciNet
K. Tateishi, Y. Yamaguchi, O. M. Abou Al-Ola, T. Kojima, and T. Yoshinaga, “Continuous analog of multiplicative algebraic reconstruction technique for computed tomography,” in Proceedings of the SPIE, Medical Imaging 2016, March 2016.
View at: Google Scholar
R. Gordon, R. Bender, and G. T. Herman, “Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and X-ray photography,” Journal of Theoretical Biology, vol. 29, no. 3, pp. 471–481, 1970.
View at: Publisher Site | Google Scholar
C. Badea and R. Gordon, “Experiments with the nonlinear and chaotic behaviour of the multiplicative algebraic reconstruction technique (MART) algorithm for computed tomography,” Physics in Medicine and Biology, vol. 49, no. 8, pp. 1455–1474, 2004.
View at: Publisher Site | Google Scholar
S. Vishampayan, J. Stamos, R. Mayans, K. Kenneth, N. Clinthorne, and W. Leslie Rogers, “Maximum likelihood image reconstruction for SPECT,” Journal of Nuclear Medicine (Abs. Book), vol. 26, no. 5, p. 20, 1985.
View at: Google Scholar
L. S. Pontryagin, Ordinary Differential Equations, Addison-Wesley, Pergamon, Turkey, 1962.
B. De Man, J. Nuyts, P. Dupont, G. Marchai, and P. Suetens, “Reduction of metal streak artifacts in x-ray computed tomography using a transmission maximum a posteriori algorithm,” IEEE Transactions on Nuclear Science, vol. 47, no. 3, pp. 977–981, 2000.
View at: Publisher Site | Google Scholar
S. Karimi, P. Cosman, C. Wald, and H. Martz, “Segmentation of artifacts and anatomy in CT metal artifact reduction,” Medical Physics, vol. 39, no. 10, pp. 5857–5868, 2012.
View at: Publisher Site | Google Scholar
J. Y. Huang, J. R. Kerns, J. L. Nute et al., “An evaluation of three commercially available metal artifact reduction methods for CT imaging,” Physics in Medicine and Biology, vol. 60, no. 3, pp. 1047–1067, 2015.
View at: Publisher Site | Google Scholar
H. S. Park, J. K. Choi, K.-R. Park et al., “Metal artifact reduction in CT by identifying missing data hidden in metals,” Journal of X-Ray Science and Technology, vol. 21, no. 3, pp. 357–372, 2013.
View at: Google Scholar
C. L. Byrne, “Iterative Algorithms in Inverse Problems,” Lecture note, 2006.
View at: Google Scholar
C. Byrne, “Block-iterative algorithms,” International Transactions in Operational Research, vol. 16, no. 4, pp. 427–463, 2009.
View at: Publisher Site | Google Scholar | MathSciNet
D. Aniszewska, “Multiplicative Runge-Kutta methods,” Nonlinear Dynamics, vol. 50, no. 1-2, pp. 265–272, 2007.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2019 Masashi Kimura et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies