- Email: [email protected]
A modified conjugate gradient method based on the Tikhonov system for computerized tomography (CT)
ISA Transactions 50 (2011) 256–261
Contents lists available at ScienceDirect
ISA Transactions journal homepage: www.elsevier.com/locate/isatrans
A modified conjugate gradient method based on the Tikhonov system for computerized tomography (CT) Qi Wang ∗ , Huaxiang Wang School of Electrical Engineering & Automation, Tianjin University, Tianjin 300072, China
article
info
Article history: Received 20 April 2010 Received in revised form 8 November 2010 Accepted 8 November 2010 Available online 3 December 2010 Keywords: CT image reconstruction Multi-phase flow Conjugate gradient Tikhonov method
abstract During the past few decades, computerized tomography (CT) was widely used for non-destructive testing (NDT) and non-destructive examination (NDE) in the industrial area because of its characteristics of noninvasiveness and visibility. Recently, CT technology has been applied to multi-phase flow measurement. Using the principle of radiation attenuation measurements along different directions through the investigated object with a special reconstruction algorithm, cross-sectional information of the scanned object can be worked out. It is a typical inverse problem and has always been a challenge for its nonlinearity and ill-conditions. The Tikhonov regulation method is widely used for similar ill-posed problems. However, the conventional Tikhonov method does not provide reconstructions with qualities good enough, the relative errors between the reconstructed images and the real distribution should be further reduced. In this paper, a modified conjugate gradient (CG) method is applied to a Tikhonov system (MCGT method) for reconstructing CT images. The computational load is dominated by the number of independent measurements m, and a preconditioner is imported to lower the condition number of the Tikhonov system. Both simulation and experiment results indicate that the proposed method can reduce the computational time and improve the quality of image reconstruction. © 2010 ISA. Published by Elsevier Ltd. All rights reserved.
1. Introduction Computerized tomography (CT) imaging is used to reconstruct an image of the interior of an object after collection of projection data from different directions at its exterior. Because of the characteristics of non-invasiveness and visibility, CT is used in the industrial area for non-destructive testing (NDT) inspections, measurements and metrology, etc. [1]. Compared with other tomography techniques, e.g. electrical tomography, optical tomography and ultrasonic tomography, CT provides higher spatial resolution up to 1% of a column diameter and even more details at the scale of 200–400 µm. In recent years, some attempts have been made to adapt such techniques to the needs of multiphase flow measurement [2–4]. In CT imaging, rays emitted from the radiation sources are able to pass through an object without bending, although with their radiation intensity attenuated. The degree of attenuation is dependent on the material composition and the size of the object [5]. Using the principle of radiation attenuation, measurements along different directions through the investigated object can reconstruct internal cross-sectional information of the scanned object.
∗
Corresponding author. Tel.: +86 022 2740 5724. E-mail addresses: [email protected] (Q. Wang), [email protected] (H. Wang).
CT problems can be modeled by a first-kind Fredholm integral Equation [6], g (θ , t ) =
+∞
∫
−∞
+∞
∫
K (x, y, θ , t )f (x, y)dxdy
(1)
−∞
where g (θ , t ) represents the integral of Beer-Lambert expression f (x, y) along the line r of orientation (θ + π /2) at displacement t from the center of the spatial coordinate system, and K is a kernel function. For simplicity here, K only accounts for the geometrical characteristics of the problem, i.e. K (x, y, θ , t ) = δ(t − x cos(θ ) − y sin(θ ))
(2)
where δ(x) is the Dirac delta function which ‘‘picks’’ out rectilinear paths that are parallel to the set of projections. It is well known that the solution of Eq. (1) in the discrete domain leads to the following linear system, which can be expressed as m×n
p = Cx + e,
i.e.,
pi =
−
ci,j xj + ei ,
i = 1, 2, . . . , m.
(3)
j=1
The vector p ∈ R m×1 contains the measured or observed projection data, while the vector x ∈ R n×1 represents the density distribution of the object, which is going to be resolved. e ∈ R m×1 is the random noise vector. The matrix C ∈ R m×n models the projection process from pixel space to projection space. The entry of C , i.e. ci,j is replaced by 1 or 0, depending upon whether the
0019-0578/$ – see front matter © 2010 ISA. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.isatra.2010.11.001
Q. Wang, H. Wang / ISA Transactions 50 (2011) 256–261
257
Fig. 1. A discrete model of the image reconstruction for CT.
jth image cell is within the ith ray. They depend only on the geometry of the problem and therefore need to be calculated only once. In this paper, as only 80 projection data can be got for reconstructing an image with the resolution of 64 × 64, i.e., m ≪ n, the reconstructed problem is an ill-posed problem, in which a solution either does not exist or is not unique or unstable. The model of image reconstruction for CT is illustrated in Fig. 1. As the matrix C is usually neither square nor full rank, the following least-squares problem is solved min ‖Cx − p‖2
(4)
where ‖µ‖2 = ( i µ2i )1/2 denotes the l2 norm of µ. Both direct methods, e.g. the Linear Back Projection (LBP) method, the Tikhonov method etc. [7,8] and iterative methods, e.g. the Algebraic Reconstruction Techniques (ART) method, the Landweber method etc. [9,10] can be utilized to solve Eq. (4), but iterative methods based on the CG iterations are more efficient when applied to large size problems [11]. However, the illconditioning of the matrix C makes the least-squares solution sensitive to noise and requires the computation of a regularized solution. In this paper, a modified CG method is applied to a Tikhonov system for flow regime identification of CT system. The reconstructed results show much improvement in the accuracy and real time performance.
∑
2. Method 2.1. Tikhonov–Arsenin method
−1 L1 =
1
..
.
..
(n−1)×n , ∈R
. −1
−1 L2 =
2
..
.
1 (6)
−1 .. .
..
−1
2
.
(n−2)×n . ∈R −1
If L is the identity matrix, i.e., the zero order derivative operator, the Tikhonov problem is said to be in standard form. The regularized solution x is then obtained by the normal equations, x = (C T C + λI )−1 C T p.
(7)
2.2. The CG method The CG method of Hestenes [12] only works with a linear system of equations with a symmetric positive definite (SPD) coefficient matrix. To apply the CG method to linearized reconstruction, Eq. (5) is transformed into C T Cx = C T p.
(8) T
As it appears from Eq. (8) the coefficient matrix is C C ∈ R n×n , where the pixel number n dominates the computational complexity of the system. Since n is always much larger than m, a modified version of the CG algorithm can be employed by introducing the dummy variable y ∈ R m×1 , the computationally expensive System (8) is transformed as,
The standard Tikhonov regularization is an effective method to deal with ill-posed problems [8], the essence of the method is to transform the solving of Eq. (3) into an optimization problem
Cx = p CC y = p
(10)
min J (x) = ‖Cx − p‖2 + λ‖Lx‖2
x = C y.
(11)
(5)
where λ is a positive number, which is called the regularization parameter, L is the matrix corresponding to a certain linear operator, which always represents the first or second derivative operator, as shown in Eq. (6)
T
T
(9)
In this case, a modified conjugate gradient (MCG) method is used in system (9) to solve y and the associated value of x can be calculated by Eq. (11) [13]. In contrast to Eq. (8), the computational load for solving Eq. (10) is dominated by the coefficient matrix
258
Q. Wang, H. Wang / ISA Transactions 50 (2011) 256–261
CC T ∈ R m×m or the number of independent measurements m, which is usually much smaller than the number of pixels n in the reconstructed image. The conventional CG method has a low convergence rate due to the ill-posed nature of matrix C . Although the MCG method mentioned above can reduce the computation time per step, it leads to large iteration number and semi-convergence in early steps. Therefore, it is necessary to seek a method which has a high convergence rate or computation efficiency, and also strong stability.
As a result, the quality of reconstructions, which is evaluated by the relative error defined in Eq. (17) of Section 3, can be improved by
2.3. The MCG method for the Tikhonov system
3.1. The CT system for multiphase flow
The Tikhonov method in Eq. (5) can be solved with direct or iterative methods. If the matrix is sparse, as is the case of the kernel of Eq. (2), a more efficient method is to apply the CG method to a Tikhonov system, which is called the MCGT method in this paper,
The system consists of a sensor head equipped with five 500 mCi 241 Am gamma-ray sources at a principal energy of 59.5 keV, as shown in Fig. 2. The five radiation sources correspond to five detector modules, each consisting of 16 CdZnTe detectors. The gamma-ray sensor head therefore comprises a total of 80 gammaray detector elements placed around an 160 mm Perspex test pipe diameter. For both simulation and experiment, oil–gas two phase flows were investigated in this paper.
D−1 (CC T + λLL T )y = p
(12)
where D = diag(CC + λLL ) is ‘‘a precondition’’ to reduce the condition number of the Tikhonov system [14], thus the illposedness of the reconstruction problem is weakened. The MCGT method has two advantages. First, given that the condition number of the system in Eq. (12) is much smaller than the system in Eq. (10), the MCGT method has a higher convergence rate than the MCG method. Second, prior information can be explicitly incorporated into the inverse computation while the CG iterations possess some intrinsic regularizing properties. T
T
2.4. Regularization parameter and regularization matrix The choice of the regularization parameter value is very important because an inappropriate parameter value can lead to a useless result. There are several methods available for choosing the regularization parameter value λ, such as the L-curve method and the generalized cross-validation (GCV) method [15], but they usually acquire a considerable amount of computation. In our application, we find from tests that with a properly chosen λ of 0.006, our algorithm can provide satisfactory reconstruction images from noisy data. The regularization matrix L, which expresses the smoothness of the solution, contributes to the quality of the reconstructed image. For simplicity, the identity matrix is adopted. 2.5. Threshold operator
xˆ i =
pi pmax − pmin xi xmax − xmin
.
(13)
We define
m ∑ pˆ i i =1 p¯ = m
n ∑
α=
j =1
(14)
xˆ j
.
n Then the threshold is xˆ t = α(1 − p¯ ).
xj < xˆ t
xj ≥ 1.
xˆ t ≤ xj < 1
xj
1
(15)
(16)
3. Results and discussion
3.2. Simulation results Simulation was implemented using Geant4 [16], which is an open source simulation software toolkit based on a Monte Carlo method. Given the geometric structure of the CT system and physical process of the radiation attenuation, the information of particle distribution and sensor response (projections) can be obtained from the detectors. The phantoms of the tested flow regime and reconstructed images are shown in Fig. 3. Fig. 3 shows that the MCGT method is superior to the other two algorithms. The iteration number is 60 for both the CG and the MCGT methods. The Tikhonov method is solved based on Eq. (7). For a simple flow regime, e.g. stratified flow and core flow, all of the three methods have similar reconstructed results. For some special or complex regime, e.g. multiple bubbles or annular flow in a tube, the MCGT method can produce images with the fewest artifacts and best image quality. The relative error and correlation coefficient between the test phantoms and reconstructions described in Eqs. (17) and (18) are used to compare reconstructions based on the Tikhonov, the CG and the MCGT algorithms. error = ‖x − xˆ ‖2 /‖x‖2
To enable the best possible reconstruction with the least relative errors, the reconstruction data are further processed by a threshold operator, which is a data processing method applied to the solution to improve the resolution of the reconstructed images. We force the range of the measured data and the solution between 0 and 1, then the measured projection vector p and the discrete CT image vector x should be normalized as,
pˆ i =
xˆ j =
0
N ∑
r =
(ˆxi − x¯ˆ )(xi − x¯ )
i =1 N ∑
(17)
(18)
N ∑ (ˆxi − x¯ˆ )2 (xi − x¯ )2
i=1
i=1 n×1
where x ∈ R , the elements of x represent the real pixel values of reconstructed images and xi is the ith component of x. xˆ is the reconstructed result and xˆ i is the ith component of xˆ . x¯ and x¯ˆ are mean values of x and xˆ based on spatial average, respectively. The minimum relative errors and maximum correlation coefficients are demonstrated in Table 1. The projection data are sampled 20 times to obtain better averaged results. As can be seen the reconstructed images of the MCGT method have a smaller relative error and larger correlation coefficient (which are shown in bold) than the Tikhonov and the CG methods. To compare the performance of the CG method and the MCGT method more quantitatively, the relative errors of reconstructed images according to each iteration step, using Eq. (15), are calculated and shown in Fig. 4, indicating that the MCGT method converges faster than the CG method. In practice, since the true images are unknown beforehand, the cut off criterion cannot be determined. A method is to stop the iteration when the residual
Q. Wang, H. Wang / ISA Transactions 50 (2011) 256–261
(a) Cross sectional view of the CT system.
259
(b) Photograph of the CT system. Fig. 2. Structure of CT system.
Fig. 3. Phantoms and reconstructed images for simulation results. Table 1 Minimum relative errors and maximum correlation coefficients related to the reconstructions in Fig. 3. Phantoms Phantom (I) Phantom (II) Phantom (III) Phantom (IV)
Tikhonov
CG
MCGT
Error = 0.4296 r = 0.5155 Error = 0.3854 r = 0.5227 Error = 0.4723 r = 0.4523 Error = 0.4452 r = 0.4833
Error = 0.2683 r = 0.6819 Error = 0.2492 r = 0.6875 Error = 0.2426 r = 0.72033 Error = 0.3098 r = 0.6213
Error = 0.0872 r = 0.8292 Error = 0.0923 r = 0.8052 Error = 0.1123 r = 0.8011 Error = 0.0824 r = 0.8322
at the kth iteration ek = x − C xk < ε, where ε is the predefined tolerance corresponding to the measurement noise level. However, it is often time consuming. Under the same experimental condition, the variation trends of relative error for reconstructing different images based on the same algorithm are similar. Thus, some specific flow regimes can be tested for choosing an iteration number with a stopping criterion. This can be done before the real experiments. From Fig. 4, we can see that, the relative error plots
Table 2 Computational time using the CG method and the MCGT method (Iteration 60). Time (s)
CG
MCGT
Phantom (I) Phantom (II) Phantom (III) Phantom (IV)
14.7134 12.5560 14.3112 15.6250
1.8281 1.7240 1.7367 1.7262
for the MCGT method decrease slowly and gradually become flat after 60 iteration steps, so the iteration number is chosen as 60 for both simulation and experiments analysis. In the previous analysis, only the iteration numbers are taken into account for evaluating the convergence rate. The convergence rate can also be measured by computational time. The time spent on iteration steps of 60 is shown in Table 2. All the computations were carried out on an Intel Pentium D PC, 2.8 GHz CPU with 1G r RAM running MATLAB⃝ version 7.0.4 on a Microsoft Windows r ⃝ XP platform. Since the computational complexity of the MCGT is much less than the CG method, the MCGT method can reduce the computation time significantly.
260
Q. Wang, H. Wang / ISA Transactions 50 (2011) 256–261
Phantom (I)
Phantom (II)
Phantom (III)
Phantom (IV) Fig. 4. Relative errors according to each iteration step.
4. Experimental results In practical measurements, the projection data would be different from the numerical ones mentioned above due to the existence of the systematic and random errors together. Following the above validation in numerical simulations, this part is to testify the MCGT method using the measured projection data in real experiments. With the CT system in Fig. 2, experiments are designed for different flow regimes such as stratified flow and core flow. Limited by the hardware in our measurement system, only stratified flow and core flow regime static experiments were carried out. The regularization matrix is the identity matrix, and the regularization parameter λ is 0.006 for both the Tikhonov and the MCGT methods. The images were reconstructed with 60 iterations for both the CG and the MCGT methods. All the true reconstructed images are shown in Fig. 5, which verifies the effectiveness of the proposed MCGT problem.
standard Tikhonov method respectively, but also reduces the computational complexity. For a better understanding of the dynamics in multiphase flows, there is obviously a current need in improving the response time of the measurement while keeping a high resolution for instantaneous observation. In dynamic systems, the reconstruction may be distorted and lags behind the actual situation, especially for the highly dynamic and deformable gas–liquid interface. We will focus on the MCGT reconstruction for dynamic flow regime in future work. Acknowledgements The authors gratefully acknowledge the National Natural Science Foundation of China (60820106002, 60532020 and 50937005) for supporting this research. The authors would also like to thank the anonymous reviewers for their careful reading and valuable comments to this paper.
5. Conclusions Appendix In this paper, a modified CG method is applied to the Tikhonov system (MCGT method) for reconstructing oil–gas two phase flow regimes using the CT technique. Both simulation and experimental results show that the MCGT method not only provides reconstructions better than the CG method and the
The CG method. 1: S = C T C Form symmetric positive definite (SPD) 2: z = S T p
Q. Wang, H. Wang / ISA Transactions 50 (2011) 256–261
261
Fig. 5. Phantoms and reconstructed images for experiment results.
3: x(0) = S −1 z Set the initial search direction and residual 4: d(0) = r(0) = Sz − Sz (0) 5: ρ(0) = r(T0) r(0) 6: i ← 0 While ‖r(i) ‖2 > ε do 7: β(i) = ρ(i) /d(Ti) Sd (i) Update the step length 8: x(i+1) = x(i) + β(i) d(i) Update the solution 9: r(i+1) = r(i) − β(i) Sd (i) Update the residual 10: ρ(i+1) = r(Ti+1) r(i+1) 11: γ(i+1) = ρ(i+1) /ρ(i) 12: d(i+1) = r(i+1) + γ(i+1) di Update the direction 13: i ← i + 1 End while The MCGT method 1: D = diag(C T C + λL T L ) 2: A = D−1 (CC T + LL T ) Precondition 3: y(0) = A−1 p Set the initial search direction and residual 4: d(0) = r(0) = p − Ay (0) 5: ρ(0) = r(T0) r(0) 6: i ← 0 While ‖r(i) ‖2 > ε do 7: β(i) = ρ(i) /d(Ti) Ad (i) Update the step length 8: y(i+1) = y(i) + β(i) d(i) Update the solution 9: r(i+1) = r(i) − β(i) Ad (i) Update the residual 10: ρ(i+1) = r(Ti+1) r(i+1) 11: γ(i+1) = ρ(i+1) /ρ(i) 12: d(i+1) = r(i+1) + γ(i+1) di Update the direction 13: i ← i + 1 End while 14: x = C T y Solve the associated value x The Tikhonov method x = (C T C + λI )−1 C T p.
References [1] Rapaport MS, Gayer A, Iszak E, Goresnic C, Baran A, Polak E. A dual-mode industrial CT. Nuclear Instruments and Methods in Physics Research Section A 1995;352:652–8. [2] Hu B, Stevart C, Hale CP, Lawrence CJ, Hall ARW, Zwiens H, Hewitt GF. Development of an X -ray computed tomography (CT) system with sparse sources: application to three-phase pipe flow visualization. Experiments Fluids 2005;39:667–78. [3] Wu CC, Cheng Y, Ding Y, Wei F, Jin Y. A novel X -ray computed tomography method for fast measurement of multiphase flow. Chemical Engineering Science 2007;62:4325–35. [4] Johansen GA, Froystein T, Hjertaker BT, Olsen O. A dual sensor flow imaging tomographic system. Measurement Science and Technology 2008;7: 297–307. [5] Rapaport MS, Gayer A, Iszak E, Goresnic C, Baranand A, Polak E. A dual-mode industrial CT. Nuclear Instruments and Methods in Physics Research Section A 1995;352:652–8. [6] Hsieh J. Computed Tomography Principle, Design, Artifacts and Recent Advances, SPIE, 2003. [7] Herman GT. Image Reconstruction from Projections: The Fundamentals of Computerized Tomography. San Francisco: Academic Press; 1980. [8] Tikhonov AN, Arsenin VY. Solutions of Ill-Posed Problems. Washington DC: Winston & Sons; 1977. [9] Kak AC, Slaney M. Principles of Computerized Topographic Imaging. New York: IEEE press; 1988. [10] Landweber L. An iteration formula for Fredholm integral Equations of the first kind. American Journal of Mathematics 1951;73:615–24. [11] Piccolomini EL, Zama F. Regularization algorithms for image reconstruction from projections, in: Proceedings of the Conference on Advanced Mathematical Tools in Metrology III, Berlin, 1996. [12] Hestenes MR. Multiplier and gradient methods. Journal of Optimization Theory and Applications 1969;4(5):303–20. [13] Vogel CR, Hanke M. Two-level preconditioners for regularized inverse problems I: theory. Numerical Mathematics 1999;83:385–402. [14] Peng LH, Lu G, Sun N. A hybrid image reconstruction algorithm for electrical capacitance tomography. AIP Conference Proceedings 2007;914(1): 807–811. [15] Vauhkonen M. Electrical impedance tomography and prior information, Ph.D. Thesies, Department of Physics, University of Kuopio, Finland, 1997. [16] Sylvia G, Gerhard L, Vladimir M, Herwig GP, Werner R. GEANT4 transport calculations for neutrons and photons below 15Mev. IEEE Transactions on Nuclear Science 2009;56: 2392–2396.
Contents lists available at ScienceDirect
ISA Transactions journal homepage: www.elsevier.com/locate/isatrans
A modified conjugate gradient method based on the Tikhonov system for computerized tomography (CT) Qi Wang ∗ , Huaxiang Wang School of Electrical Engineering & Automation, Tianjin University, Tianjin 300072, China
article
info
Article history: Received 20 April 2010 Received in revised form 8 November 2010 Accepted 8 November 2010 Available online 3 December 2010 Keywords: CT image reconstruction Multi-phase flow Conjugate gradient Tikhonov method
abstract During the past few decades, computerized tomography (CT) was widely used for non-destructive testing (NDT) and non-destructive examination (NDE) in the industrial area because of its characteristics of noninvasiveness and visibility. Recently, CT technology has been applied to multi-phase flow measurement. Using the principle of radiation attenuation measurements along different directions through the investigated object with a special reconstruction algorithm, cross-sectional information of the scanned object can be worked out. It is a typical inverse problem and has always been a challenge for its nonlinearity and ill-conditions. The Tikhonov regulation method is widely used for similar ill-posed problems. However, the conventional Tikhonov method does not provide reconstructions with qualities good enough, the relative errors between the reconstructed images and the real distribution should be further reduced. In this paper, a modified conjugate gradient (CG) method is applied to a Tikhonov system (MCGT method) for reconstructing CT images. The computational load is dominated by the number of independent measurements m, and a preconditioner is imported to lower the condition number of the Tikhonov system. Both simulation and experiment results indicate that the proposed method can reduce the computational time and improve the quality of image reconstruction. © 2010 ISA. Published by Elsevier Ltd. All rights reserved.
1. Introduction Computerized tomography (CT) imaging is used to reconstruct an image of the interior of an object after collection of projection data from different directions at its exterior. Because of the characteristics of non-invasiveness and visibility, CT is used in the industrial area for non-destructive testing (NDT) inspections, measurements and metrology, etc. [1]. Compared with other tomography techniques, e.g. electrical tomography, optical tomography and ultrasonic tomography, CT provides higher spatial resolution up to 1% of a column diameter and even more details at the scale of 200–400 µm. In recent years, some attempts have been made to adapt such techniques to the needs of multiphase flow measurement [2–4]. In CT imaging, rays emitted from the radiation sources are able to pass through an object without bending, although with their radiation intensity attenuated. The degree of attenuation is dependent on the material composition and the size of the object [5]. Using the principle of radiation attenuation, measurements along different directions through the investigated object can reconstruct internal cross-sectional information of the scanned object.
∗
Corresponding author. Tel.: +86 022 2740 5724. E-mail addresses: [email protected] (Q. Wang), [email protected] (H. Wang).
CT problems can be modeled by a first-kind Fredholm integral Equation [6], g (θ , t ) =
+∞
∫
−∞
+∞
∫
K (x, y, θ , t )f (x, y)dxdy
(1)
−∞
where g (θ , t ) represents the integral of Beer-Lambert expression f (x, y) along the line r of orientation (θ + π /2) at displacement t from the center of the spatial coordinate system, and K is a kernel function. For simplicity here, K only accounts for the geometrical characteristics of the problem, i.e. K (x, y, θ , t ) = δ(t − x cos(θ ) − y sin(θ ))
(2)
where δ(x) is the Dirac delta function which ‘‘picks’’ out rectilinear paths that are parallel to the set of projections. It is well known that the solution of Eq. (1) in the discrete domain leads to the following linear system, which can be expressed as m×n
p = Cx + e,
i.e.,
pi =
−
ci,j xj + ei ,
i = 1, 2, . . . , m.
(3)
j=1
The vector p ∈ R m×1 contains the measured or observed projection data, while the vector x ∈ R n×1 represents the density distribution of the object, which is going to be resolved. e ∈ R m×1 is the random noise vector. The matrix C ∈ R m×n models the projection process from pixel space to projection space. The entry of C , i.e. ci,j is replaced by 1 or 0, depending upon whether the
0019-0578/$ – see front matter © 2010 ISA. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.isatra.2010.11.001
Q. Wang, H. Wang / ISA Transactions 50 (2011) 256–261
257
Fig. 1. A discrete model of the image reconstruction for CT.
jth image cell is within the ith ray. They depend only on the geometry of the problem and therefore need to be calculated only once. In this paper, as only 80 projection data can be got for reconstructing an image with the resolution of 64 × 64, i.e., m ≪ n, the reconstructed problem is an ill-posed problem, in which a solution either does not exist or is not unique or unstable. The model of image reconstruction for CT is illustrated in Fig. 1. As the matrix C is usually neither square nor full rank, the following least-squares problem is solved min ‖Cx − p‖2
(4)
where ‖µ‖2 = ( i µ2i )1/2 denotes the l2 norm of µ. Both direct methods, e.g. the Linear Back Projection (LBP) method, the Tikhonov method etc. [7,8] and iterative methods, e.g. the Algebraic Reconstruction Techniques (ART) method, the Landweber method etc. [9,10] can be utilized to solve Eq. (4), but iterative methods based on the CG iterations are more efficient when applied to large size problems [11]. However, the illconditioning of the matrix C makes the least-squares solution sensitive to noise and requires the computation of a regularized solution. In this paper, a modified CG method is applied to a Tikhonov system for flow regime identification of CT system. The reconstructed results show much improvement in the accuracy and real time performance.
∑
2. Method 2.1. Tikhonov–Arsenin method
−1 L1 =
1
..
.
..
(n−1)×n , ∈R
. −1
−1 L2 =
2
..
.
1 (6)
−1 .. .
..
−1
2
.
(n−2)×n . ∈R −1
If L is the identity matrix, i.e., the zero order derivative operator, the Tikhonov problem is said to be in standard form. The regularized solution x is then obtained by the normal equations, x = (C T C + λI )−1 C T p.
(7)
2.2. The CG method The CG method of Hestenes [12] only works with a linear system of equations with a symmetric positive definite (SPD) coefficient matrix. To apply the CG method to linearized reconstruction, Eq. (5) is transformed into C T Cx = C T p.
(8) T
As it appears from Eq. (8) the coefficient matrix is C C ∈ R n×n , where the pixel number n dominates the computational complexity of the system. Since n is always much larger than m, a modified version of the CG algorithm can be employed by introducing the dummy variable y ∈ R m×1 , the computationally expensive System (8) is transformed as,
The standard Tikhonov regularization is an effective method to deal with ill-posed problems [8], the essence of the method is to transform the solving of Eq. (3) into an optimization problem
Cx = p CC y = p
(10)
min J (x) = ‖Cx − p‖2 + λ‖Lx‖2
x = C y.
(11)
(5)
where λ is a positive number, which is called the regularization parameter, L is the matrix corresponding to a certain linear operator, which always represents the first or second derivative operator, as shown in Eq. (6)
T
T
(9)
In this case, a modified conjugate gradient (MCG) method is used in system (9) to solve y and the associated value of x can be calculated by Eq. (11) [13]. In contrast to Eq. (8), the computational load for solving Eq. (10) is dominated by the coefficient matrix
258
Q. Wang, H. Wang / ISA Transactions 50 (2011) 256–261
CC T ∈ R m×m or the number of independent measurements m, which is usually much smaller than the number of pixels n in the reconstructed image. The conventional CG method has a low convergence rate due to the ill-posed nature of matrix C . Although the MCG method mentioned above can reduce the computation time per step, it leads to large iteration number and semi-convergence in early steps. Therefore, it is necessary to seek a method which has a high convergence rate or computation efficiency, and also strong stability.
As a result, the quality of reconstructions, which is evaluated by the relative error defined in Eq. (17) of Section 3, can be improved by
2.3. The MCG method for the Tikhonov system
3.1. The CT system for multiphase flow
The Tikhonov method in Eq. (5) can be solved with direct or iterative methods. If the matrix is sparse, as is the case of the kernel of Eq. (2), a more efficient method is to apply the CG method to a Tikhonov system, which is called the MCGT method in this paper,
The system consists of a sensor head equipped with five 500 mCi 241 Am gamma-ray sources at a principal energy of 59.5 keV, as shown in Fig. 2. The five radiation sources correspond to five detector modules, each consisting of 16 CdZnTe detectors. The gamma-ray sensor head therefore comprises a total of 80 gammaray detector elements placed around an 160 mm Perspex test pipe diameter. For both simulation and experiment, oil–gas two phase flows were investigated in this paper.
D−1 (CC T + λLL T )y = p
(12)
where D = diag(CC + λLL ) is ‘‘a precondition’’ to reduce the condition number of the Tikhonov system [14], thus the illposedness of the reconstruction problem is weakened. The MCGT method has two advantages. First, given that the condition number of the system in Eq. (12) is much smaller than the system in Eq. (10), the MCGT method has a higher convergence rate than the MCG method. Second, prior information can be explicitly incorporated into the inverse computation while the CG iterations possess some intrinsic regularizing properties. T
T
2.4. Regularization parameter and regularization matrix The choice of the regularization parameter value is very important because an inappropriate parameter value can lead to a useless result. There are several methods available for choosing the regularization parameter value λ, such as the L-curve method and the generalized cross-validation (GCV) method [15], but they usually acquire a considerable amount of computation. In our application, we find from tests that with a properly chosen λ of 0.006, our algorithm can provide satisfactory reconstruction images from noisy data. The regularization matrix L, which expresses the smoothness of the solution, contributes to the quality of the reconstructed image. For simplicity, the identity matrix is adopted. 2.5. Threshold operator
xˆ i =
pi pmax − pmin xi xmax − xmin
.
(13)
We define
m ∑ pˆ i i =1 p¯ = m
n ∑
α=
j =1
(14)
xˆ j
.
n Then the threshold is xˆ t = α(1 − p¯ ).
xj < xˆ t
xj ≥ 1.
xˆ t ≤ xj < 1
xj
1
(15)
(16)
3. Results and discussion
3.2. Simulation results Simulation was implemented using Geant4 [16], which is an open source simulation software toolkit based on a Monte Carlo method. Given the geometric structure of the CT system and physical process of the radiation attenuation, the information of particle distribution and sensor response (projections) can be obtained from the detectors. The phantoms of the tested flow regime and reconstructed images are shown in Fig. 3. Fig. 3 shows that the MCGT method is superior to the other two algorithms. The iteration number is 60 for both the CG and the MCGT methods. The Tikhonov method is solved based on Eq. (7). For a simple flow regime, e.g. stratified flow and core flow, all of the three methods have similar reconstructed results. For some special or complex regime, e.g. multiple bubbles or annular flow in a tube, the MCGT method can produce images with the fewest artifacts and best image quality. The relative error and correlation coefficient between the test phantoms and reconstructions described in Eqs. (17) and (18) are used to compare reconstructions based on the Tikhonov, the CG and the MCGT algorithms. error = ‖x − xˆ ‖2 /‖x‖2
To enable the best possible reconstruction with the least relative errors, the reconstruction data are further processed by a threshold operator, which is a data processing method applied to the solution to improve the resolution of the reconstructed images. We force the range of the measured data and the solution between 0 and 1, then the measured projection vector p and the discrete CT image vector x should be normalized as,
pˆ i =
xˆ j =
0
N ∑
r =
(ˆxi − x¯ˆ )(xi − x¯ )
i =1 N ∑
(17)
(18)
N ∑ (ˆxi − x¯ˆ )2 (xi − x¯ )2
i=1
i=1 n×1
where x ∈ R , the elements of x represent the real pixel values of reconstructed images and xi is the ith component of x. xˆ is the reconstructed result and xˆ i is the ith component of xˆ . x¯ and x¯ˆ are mean values of x and xˆ based on spatial average, respectively. The minimum relative errors and maximum correlation coefficients are demonstrated in Table 1. The projection data are sampled 20 times to obtain better averaged results. As can be seen the reconstructed images of the MCGT method have a smaller relative error and larger correlation coefficient (which are shown in bold) than the Tikhonov and the CG methods. To compare the performance of the CG method and the MCGT method more quantitatively, the relative errors of reconstructed images according to each iteration step, using Eq. (15), are calculated and shown in Fig. 4, indicating that the MCGT method converges faster than the CG method. In practice, since the true images are unknown beforehand, the cut off criterion cannot be determined. A method is to stop the iteration when the residual
Q. Wang, H. Wang / ISA Transactions 50 (2011) 256–261
(a) Cross sectional view of the CT system.
259
(b) Photograph of the CT system. Fig. 2. Structure of CT system.
Fig. 3. Phantoms and reconstructed images for simulation results. Table 1 Minimum relative errors and maximum correlation coefficients related to the reconstructions in Fig. 3. Phantoms Phantom (I) Phantom (II) Phantom (III) Phantom (IV)
Tikhonov
CG
MCGT
Error = 0.4296 r = 0.5155 Error = 0.3854 r = 0.5227 Error = 0.4723 r = 0.4523 Error = 0.4452 r = 0.4833
Error = 0.2683 r = 0.6819 Error = 0.2492 r = 0.6875 Error = 0.2426 r = 0.72033 Error = 0.3098 r = 0.6213
Error = 0.0872 r = 0.8292 Error = 0.0923 r = 0.8052 Error = 0.1123 r = 0.8011 Error = 0.0824 r = 0.8322
at the kth iteration ek = x − C xk < ε, where ε is the predefined tolerance corresponding to the measurement noise level. However, it is often time consuming. Under the same experimental condition, the variation trends of relative error for reconstructing different images based on the same algorithm are similar. Thus, some specific flow regimes can be tested for choosing an iteration number with a stopping criterion. This can be done before the real experiments. From Fig. 4, we can see that, the relative error plots
Table 2 Computational time using the CG method and the MCGT method (Iteration 60). Time (s)
CG
MCGT
Phantom (I) Phantom (II) Phantom (III) Phantom (IV)
14.7134 12.5560 14.3112 15.6250
1.8281 1.7240 1.7367 1.7262
for the MCGT method decrease slowly and gradually become flat after 60 iteration steps, so the iteration number is chosen as 60 for both simulation and experiments analysis. In the previous analysis, only the iteration numbers are taken into account for evaluating the convergence rate. The convergence rate can also be measured by computational time. The time spent on iteration steps of 60 is shown in Table 2. All the computations were carried out on an Intel Pentium D PC, 2.8 GHz CPU with 1G r RAM running MATLAB⃝ version 7.0.4 on a Microsoft Windows r ⃝ XP platform. Since the computational complexity of the MCGT is much less than the CG method, the MCGT method can reduce the computation time significantly.
260
Q. Wang, H. Wang / ISA Transactions 50 (2011) 256–261
Phantom (I)
Phantom (II)
Phantom (III)
Phantom (IV) Fig. 4. Relative errors according to each iteration step.
4. Experimental results In practical measurements, the projection data would be different from the numerical ones mentioned above due to the existence of the systematic and random errors together. Following the above validation in numerical simulations, this part is to testify the MCGT method using the measured projection data in real experiments. With the CT system in Fig. 2, experiments are designed for different flow regimes such as stratified flow and core flow. Limited by the hardware in our measurement system, only stratified flow and core flow regime static experiments were carried out. The regularization matrix is the identity matrix, and the regularization parameter λ is 0.006 for both the Tikhonov and the MCGT methods. The images were reconstructed with 60 iterations for both the CG and the MCGT methods. All the true reconstructed images are shown in Fig. 5, which verifies the effectiveness of the proposed MCGT problem.
standard Tikhonov method respectively, but also reduces the computational complexity. For a better understanding of the dynamics in multiphase flows, there is obviously a current need in improving the response time of the measurement while keeping a high resolution for instantaneous observation. In dynamic systems, the reconstruction may be distorted and lags behind the actual situation, especially for the highly dynamic and deformable gas–liquid interface. We will focus on the MCGT reconstruction for dynamic flow regime in future work. Acknowledgements The authors gratefully acknowledge the National Natural Science Foundation of China (60820106002, 60532020 and 50937005) for supporting this research. The authors would also like to thank the anonymous reviewers for their careful reading and valuable comments to this paper.
5. Conclusions Appendix In this paper, a modified CG method is applied to the Tikhonov system (MCGT method) for reconstructing oil–gas two phase flow regimes using the CT technique. Both simulation and experimental results show that the MCGT method not only provides reconstructions better than the CG method and the
The CG method. 1: S = C T C Form symmetric positive definite (SPD) 2: z = S T p
Q. Wang, H. Wang / ISA Transactions 50 (2011) 256–261
261
Fig. 5. Phantoms and reconstructed images for experiment results.
3: x(0) = S −1 z Set the initial search direction and residual 4: d(0) = r(0) = Sz − Sz (0) 5: ρ(0) = r(T0) r(0) 6: i ← 0 While ‖r(i) ‖2 > ε do 7: β(i) = ρ(i) /d(Ti) Sd (i) Update the step length 8: x(i+1) = x(i) + β(i) d(i) Update the solution 9: r(i+1) = r(i) − β(i) Sd (i) Update the residual 10: ρ(i+1) = r(Ti+1) r(i+1) 11: γ(i+1) = ρ(i+1) /ρ(i) 12: d(i+1) = r(i+1) + γ(i+1) di Update the direction 13: i ← i + 1 End while The MCGT method 1: D = diag(C T C + λL T L ) 2: A = D−1 (CC T + LL T ) Precondition 3: y(0) = A−1 p Set the initial search direction and residual 4: d(0) = r(0) = p − Ay (0) 5: ρ(0) = r(T0) r(0) 6: i ← 0 While ‖r(i) ‖2 > ε do 7: β(i) = ρ(i) /d(Ti) Ad (i) Update the step length 8: y(i+1) = y(i) + β(i) d(i) Update the solution 9: r(i+1) = r(i) − β(i) Ad (i) Update the residual 10: ρ(i+1) = r(Ti+1) r(i+1) 11: γ(i+1) = ρ(i+1) /ρ(i) 12: d(i+1) = r(i+1) + γ(i+1) di Update the direction 13: i ← i + 1 End while 14: x = C T y Solve the associated value x The Tikhonov method x = (C T C + λI )−1 C T p.
References [1] Rapaport MS, Gayer A, Iszak E, Goresnic C, Baran A, Polak E. A dual-mode industrial CT. Nuclear Instruments and Methods in Physics Research Section A 1995;352:652–8. [2] Hu B, Stevart C, Hale CP, Lawrence CJ, Hall ARW, Zwiens H, Hewitt GF. Development of an X -ray computed tomography (CT) system with sparse sources: application to three-phase pipe flow visualization. Experiments Fluids 2005;39:667–78. [3] Wu CC, Cheng Y, Ding Y, Wei F, Jin Y. A novel X -ray computed tomography method for fast measurement of multiphase flow. Chemical Engineering Science 2007;62:4325–35. [4] Johansen GA, Froystein T, Hjertaker BT, Olsen O. A dual sensor flow imaging tomographic system. Measurement Science and Technology 2008;7: 297–307. [5] Rapaport MS, Gayer A, Iszak E, Goresnic C, Baranand A, Polak E. A dual-mode industrial CT. Nuclear Instruments and Methods in Physics Research Section A 1995;352:652–8. [6] Hsieh J. Computed Tomography Principle, Design, Artifacts and Recent Advances, SPIE, 2003. [7] Herman GT. Image Reconstruction from Projections: The Fundamentals of Computerized Tomography. San Francisco: Academic Press; 1980. [8] Tikhonov AN, Arsenin VY. Solutions of Ill-Posed Problems. Washington DC: Winston & Sons; 1977. [9] Kak AC, Slaney M. Principles of Computerized Topographic Imaging. New York: IEEE press; 1988. [10] Landweber L. An iteration formula for Fredholm integral Equations of the first kind. American Journal of Mathematics 1951;73:615–24. [11] Piccolomini EL, Zama F. Regularization algorithms for image reconstruction from projections, in: Proceedings of the Conference on Advanced Mathematical Tools in Metrology III, Berlin, 1996. [12] Hestenes MR. Multiplier and gradient methods. Journal of Optimization Theory and Applications 1969;4(5):303–20. [13] Vogel CR, Hanke M. Two-level preconditioners for regularized inverse problems I: theory. Numerical Mathematics 1999;83:385–402. [14] Peng LH, Lu G, Sun N. A hybrid image reconstruction algorithm for electrical capacitance tomography. AIP Conference Proceedings 2007;914(1): 807–811. [15] Vauhkonen M. Electrical impedance tomography and prior information, Ph.D. Thesies, Department of Physics, University of Kuopio, Finland, 1997. [16] Sylvia G, Gerhard L, Vladimir M, Herwig GP, Werner R. GEANT4 transport calculations for neutrons and photons below 15Mev. IEEE Transactions on Nuclear Science 2009;56: 2392–2396.