The Conjugate Gradient Regularization Method in ... - Semantic Scholar
The Conjugate Gradient Regularization Method in Computed Tomography Problems ELENA LOLI PICCOLOMINI, and FABIANA ZAMA
Department of Mathematics, University of Bologna, Italy e{mail: [email protected], [email protected] ABSTRACT In this work we solve inverse problems coming from the area of Computed Tomography by means of regularization methods based on conjugate gradient iterations. We develop a stopping criterion which is ecient for the computation of a regularized solution for the least{squares normal equations. The stopping rule can be suitably applied also to the Tikhonov regularization method. We report computational experiences based on di erent physical models and with di erent degrees of noise. We compare the results obtained with those computed by other currently used methods such as Algebraic reconstruction Techinques (ART) and Backprojection.
1. Introduction The problems arising in the Computed Tomography area are well known for their high dimensions and ill{posedness. There are di erent kinds of Computed Tomography problems (TAC, SPECT, PET), which are all modeled by a rst kind Fredholm integral equation: Z1 Z1 g(; t) = K(x; y; ; t)f(x; y)dxdy (1) ?1 ?1
In practical problems the unknown function f(x; y) represents di erent physical phenomena (attenuation in TAC, emission in SPECT and PET) and the kernel K(x; y; ; t) may have di erent expressions. In the easiest Department of Mathematics, Piazza di Porta San Donato,5 40127 Bologna (ITALY)
APPLIED MATHEMATICS AND COMPUTATION xx:1{xxx (1994) 1
c Elsevier Science Inc., 1994 655 Avenue of the Americas, New York, NY 10010 0096-3003/94/$6.00
2 models it only accounts for the geometrical characteristics of the problem, i.e.,[1]: K(x; y; ; t) = (t ? xcos() ? ysin()) (2) and the value of g(; t) represents the integral of f(x; y) along the line r of orientation + =2 at displacement t from the center of the spatial coordinate system. In more realistic models, in which other physical phenomena are considered, it represents a gaussian response to an impulse. For example, in the case of SPECT, the following kernel has been used[13]: z(x; y; ; t) 2 p K(x; y; ; t) = 1=( 2(x; y; ; t)) exp(?0:5 (x; y; ; t) ) (3) where is the standard deviation of the Gaussian system response and z is a function of the distance of the point (x; y) from the line r. It is well known that (1) is an ill{posed problem and the solution of the integral equation in the discrete domain leads to the following linear system: Kf = g (4) where the matrix K is usually neither square nor full rank. For this reason and for the presence of noise in g, the following least{squares problem is solved: min k Kf ? g k2 (5) The ill{conditioning of the matrix K makes the least{squares solution useless and requires the computation of a regularized solution. Both direct and iterative methods can be utilized [8], but iterative methods based on the conjugate gradient iterations are more ecient when applied to large size problems [14]. In this work we compute least{squares regularized solutions by means of regularization methods based on conjugate gradient iterations. In particular, in section 2 we show the properties of the Conjugate Gradient method and derive a suitable stopping criterion to compute eciently a regularized solution. In section 3 the same criterion is extended to the Tikhonov regularization method. Computational experiments are reported in section 4. 2. The Conjugate Gradient Method The matrix K of (5) is usually not full rank, severely ill{conditioned and has a few very small singular values . Although the Conjugate Gradient method applied to the least{squares normal equations (CGNE) does not necessarily converge, it has been observed that the approximations seem
3 to converge in the rst iterations, before they start oscillating and nally diverge. This phenomenon is known as semiconvergence and the best approximation obtained is called a regularized solution of the problem ([8], [9]). When the number of iterations exceeds a threshold value, the solution diverges from the regularized one and the method must be stopped. Hence the number of iterations plays the role of regularization parameter for the method. The merits of the CGNE as a regularization method depend on the determination of the optimal stopping iteration index k , that is generally related to the conditioning of K and to the speci c problem. The ill{conditioning of K can be analysed in terms of the Singular Values Decomposition (SVD) of the matrix ([12], [11]). If q singular values are comparatively large with respect to the other singular values, i.e. the problem is extremely ill{posed (its singular values decay very fast as ri with 0 < r < 1 and i > 1), then after about q iterations the solution will be in uenced by noise and will diverge. Alternatively if the singular values are almost uniformly distributed over the interval the regularized solution will be reached in a much large number of iterations. In our problems, only the last few singular values decay very fast ( gure 4) and hence the iterations must be stopped for k < q. We have determined a stopping criterion on the basis of the behaviour of the following successions computed during the CGNE iterations: k = k xk ? xk?1 k2 k = k rk k2
(6)
where xk is the approximate solution and rk is the residual at the k?th step. We compare the behaviour of these successions with that of the euclidean norm of the absolute error: E2k =k xk ? f k2. In the case of a positive de nite well conditioned matrix both the successions (6) decrease monotonically, but in the case of an ill{posed problem the successions k and k have a typical oscillating behaviour while E2k reaches a minimum and then start to diverge. We report, as an example, the plots of the singular values of K representing 16 20 projections of a 16 16 image f ( gure 4). In gures 1, 2 and 3 the successions k , k and E2k are reported for a noisy g when a kernel of the kind (2) is used, while the same plots for a kernel (3) are not reported because they are very similar. We used these observations to derive two almost equivalent stopping rules: if k > k?1 (1 < 10) then stop. if k > k?1 (1 < 4) then stop. The values of have been sperimentally determined for these kind of problems; they essentially depend on the matrix K and on the data g.
4 || r ||
3
10
2
10
1
10
0
10
−1
10
−2
10
−3
10
−4
10
0
100
200
300 iteration
400
500
600
500
600
. Succession .
Fig. 1
k
|| x_k − x_(k−1) ||
1
10
0
10
−1
10
−2
10
−3
10
−4
10
−5
10
0
100
Fig. 2
200
300 iteration
. Succession
400
k
3. Tikhonov Method Another highly regarded method for solving (5) is the Tikhonov regularization method, which corresponds to the solution of: min k K f ? g k2 + k Qf k2 where Q is a suitable regularization matrix and is the regularization parameter. The regularized solution f is then obtained by the normal equations: (K t K + Qt Q)f = K t g (7) In practical applications the regularization matrix Q is the discretization of some di erential operator which expresses the smoothness of the solution.
5 E2
1
10
0
10
0
100
200
Fig. 3
300 iteration
400
500
600
200
250
300
. E2
singular values
2
10
0
10
−2
10
−4
10
−6
10
−8
10
−10
10
−12
10
−14
10
−16
10
0
50
Fig. 4
100
150
. Singular values.
One choice commonly used in tomographic reconstruction [3] is to de ne the regularization matrix by means of the discrepancy between the value of each image point and its nearest neighbours: (8) (f ) = 21 f T(I + C tC)f where y = C f has the following expression: yj = fj ? 18 (fj ?n?1 + fj ?n + fj ?n+1 + fj ?1 + fj +1 + fj +n?1 + fj +n + fj +n+1) and is a positive weight factor. Since the matrix (I + C tC) results positive de nite we can determine Q such that: Qt Q = (I + C tC). The system (7) can be solved with direct or iterative methods. If the matrix is sparse, as is the case of kernel (2), the conjugate gradient method is more ecient, since it accounts for the sparsity of the matrix and has a
6 lower computational complexity if compared with a matrix factorization. Moreover, even in this case a suitable stopping criteria can be adopted in order to obtain a regularized solution before the convergence of the method. The behaviour of the iterations depend strongly on the regularization parameter and on the function used. When the regularization parameter is small the successions k and k behave like in the CGNE previously described, while for larger values of the regularization parameter they decrease smoothly ( gures 5, 6). We stopped the iterations of the conjugate gradient in Tikhonov method using the same CGNE criterion for small values of , while the the usual convergence criterion for conjugate gradients is used when is large and the problem conditioning is far better. For the computation of the optimal regularization parameter there are several proposals in literature, such as the L{curve method, the Generalized Cross Validation method or the discrepancy principle. They usually requires the solutions f for di erent values of and a criterion to select the best solution obtained, hence they are computationally very expensive and some counterexamples show that the methods are not always precises. 4. Numerical Results In this section we give some numerical results using simulations of brain sections (phantoms). We used the following test images: P1 95 95 samples at uniform distance 1:6 obtained by using snark93 test problem b.7 ( gure 7). P2 32 32 simulated brain section (Ho man phantom) ( gure 10). The projection data are obtained by applying (4) where K has the values given in (2) for problem P1 and (3) for problem P2. When an additive Gaussian white noise is introduced in the projection data it is measured by means of the noise level parameter: NL = k gk ?g kg~ k2 2 where g~ is the noisy signal. The pictures ~f reconstructed with di erent methods are compared with the original images f by means of the following error parameters [1]: v u P P jf ? ~f j P P (f ? ~f )2 u u =1 =1 =1 =1 NRMSE = u u t P P (f ? f )2 NMAE = P P jf j u;v
n
v
u;v
u;v
v
u
v
n
n
n
n
n
=1 u=1
u;v
u
n
u;v
v
n
=1 u=1
u;v
7
1
10
0
10
−1
10
−2
10
0
5
10
15
20
25
30
35
40
45
Fig. 5. Euclidean norm of the residuals ( = 0:06). 3
10
2
10
1
10
0
10
−1
10
−2
10
−3
10
−4
10
−5
10
0
5
10
15
20
25
30
Fig. 6. Euclidean norm of the residuals ( = 256).
35
40
8 where f is the average of the values in the original image f. In order to evaluate the reconstruction quality we report in table 1 the results obtained for P1 with four methods: ART, Backprojection, CGNE and Tikhonov methods. We observe that the errors increase when the number of projections decrease, and in particular ART and Backprojection are more sensitive to the number of projections. Noise also increases the errors, expecially for ART and Backprojection, while CGNE and Tikhonov (TIKH) methods give better results in more unfavourable situations (high noise, few projections). In table 2 we report the results obtained with CGNE algorithm for di erent values of the parameter . By increasing the value of the number of iterations increases while the errors reach a minimum value and then start to increase (as shown in section 2). The results obtained with Tikhonov method are reported in table 3 with di erent values of the regularization parameter . We observe that the number of iterations is greater with respect to CGNE and in general each iteration is computationally more expensive. In the case of test P1 the errors are a little bit better with respect to CGNE but the number of iterations increases signi cantly. Some reconstructed images are shown in gures 8, 9, 11 and 12. 5. Conclusions We can conclude that CGNE is the preferable regularization method for the problems in the area of Computed Tomography not only for the quality of the results obtained but also for its eciency (number of iterations) and the easy parameter setting. In fact, while the setting of the regularization Tikhonov parameter is dicult and costly (section 3), only a few values of
are sucient to obtain the best results with CGNE. REFERENCES 1 G.T Herman, Image Reconstruction from Projections, Academic Press, 1980. 2 F. Natterer, The Mathematics of Computerized Tomography, John Wiley & Sons 1989. 3 R.L Kayshap and C. Mittal, IEEE Trans. on Computers 24, 915 (1975) 4 P.C.Hansen, Numerical Algorithms 6, 1 (1994) 5 J.M.Varah, SIAM J. Numer. Anal. 10, 257 (1973)
9
7. 95 95 Phantom
Fig.
Fig. 8. = 50 Backprojection. k
9. 50proj: CGNE 25 it.
Fig.
10
10. 32 32 Ho man Phantom
Fig.
Fig. 11. 50proj: CGNE 29 it.
Fig. 12. 50proj: Tikhonov
11 6 7 8 9 10 11 12 13 14
P.C.Hansen, SIAM J. Sci. Stat. Comp. 11, 503 (1990) P.C.Hansen, BIT 27, 543 (1987) M. Hanke, P.C.Hansen, Surveys on Mathematics for Industry 3, 543 (1993) M. Hanke Pitman Reserch Notes in Mathematics , 327 (1995) C. C. Paige, M. A. Saunders,Anal Comp Math 8, 43 (1982) R. Plato, G. Vainikko, Numer. Funct. Anal. Optim. 11, 111 (1990) A. van der Sluis, H. A. Van der Vorst, Linear Algebra Appl. 130, 257 (1990) Formiconi et alt. Phys. Med. Biol. 37, (1992) E. Loli Piccolomini, F. Zama, Regularization Algorithms for Image Reconstruction from Projections, Proc. of the Conf. "Advanced Mathematical Tools in Metrology III", Berlin, 1996.
12
. Test Image 95 95.
Table 1
Proj. NL 100 0.
50
25
2. 2. 2. 2. 0. 0. 0. 0. 2. 2. 2. 2. 0. 0. 0. 0. 2. 2. 2. 2.
Method NRMSE NMAE Iter ART 0.61 0.42 11 BPRJ 0.14 9.3(-2) CGNE 0.97(-1) 0.52(-1) 16 TIKH(1.e-3) 0.60(-1) 0.1 18 ART 0.88 0.59 96 BPRJ 0.63 0.42 CGNE 0.29 0.19 9 TIKH (1.e+2) 0.22 0.14 17 ART 0.44 0.29 32 BPRJ 0.24 0.16 CGNE 0.15 9.1(-2) 9 TIKH(1.e-3) 0.15 9.1(-2) 9 ART 0.84 0.55 28 BPRJ 0.74 0.49 CGNE 0.33 0.33 11 TIKH(0.25e+2) 0.22 0.13 11 ART 0.46 0.31 25 BPRJ 0.42 0.27 CGNE 0.21 0.13 9 TIKH(1.e-3) 0.21 0.13 11 ART 0.85 0.56 18 BPRJ 0.83 0.55 CGNE 0.34 0.21 9 TIKH(0.25e+2) 0.26 0.18 11
13
. Conjugate Gradient Iterations.
Table 2
Test proj P1 100
NL 0
P1
100 5.9e-2
P1
100
0.59
P2
50
0
P2
50
5.0e-3
P2
50
5.0e-2 5.0e-2
P2
50
0.5
0.10E+01 0.22E+01 0.72E+01 0.10E+02 0.10E+01 0.32E+01 0.52E+01 0.10E+02 0.10E+01 0.12E+01 0.32E+01 0.12E+01 0.22E+01 0.32E+01 0.12E+01 0.26E+01 0.32E+01 0.12E+01 0.22E+01 0.32E+01 0.12E+01 0.22E+01 0.32E+01
iter 9 16 84 100 9 29 81 100 9 16 100 5 29 100 5 73 99 5 29 95 5 33 127
Emax 0.733E+02 0.650E+02 0.434E+02 0.398E+02 0.728E+02 0.581E+02 0.444E+02 0.389E+02 0.769E+02 0.748E+02 0.183E+03 0.120E+03 0.112E+03 0.109E+03 0.120E+03 0.109E+03 0.109E+03 0.120E+03 0.111E+03 0.108E+03 0.120E+03 0.107E+03 0.180E+03
NRMSE 0.117E+00 0.970E-01 0.660E-01 0.637E-01 0.117E+00 0.844E-01 0.731E-01 0.738E-01 0.124E+00 0.115E+00 0.466E+00 0.681E+00 0.393E+00 0.219E+00 0.681E+00 0.256E+00 0.220E+00 0.681E+00 0.393E+00 0.228E+00 0.681E+00 0.404E+00 0.555E+00
NMAE 0.615E-01 0.519E-01 0.362E-01 0.349E-01 0.615E-01 0.454E-01 0.405E-01 0.412E-01 0.665E-01 0.638E-01 0.262E+00 0.453E+00 0.203E+00 0.884E-01 0.453E+00 0.112E+00 0.891E-01 0.453E+00 0.203E+00 0.966E-01 0.453E+00 0.227E+00 0.278E+00
14
. Tikhonov Method
Table 3
Test proj. NL P1 100 5.9e-2 0.10E-02 0.10E-01 0.10E+00 P1 100 0.59 0.10E-02 0.10E-01 0.10E+00 0.10E-02 0.10E+00 0.10E+01 0.10E+02 0.15E+02 0.50E+02 P2 50 0.5 0.10E-01 0.10E+00 0.10E+01 0.10E+02
0.32E+01 0.32E+01 0.32E+01 0.10E+01 0.10E+01 0.10E+01 0.32E+01 0.32E+01 0.32E+01 0.32E+01 0.32E+01 0.32E+01 0.22E+01 0.22E+01 0.32E+01 0.32E+01
iter 68 68 17 9 9 9 16 29 30 69 54 29 33 33 56 23
Emax 0.492E+02 0.491E+02 0.628E+02 0.769E+02 0.769E+02 0.769E+02 0.748E+02 0.735E+02 0.717E+02 0.676E+02 0.681E+02 0.715E+02 0.107E+03 0.103E+03 0.106E+03 0.117E+03
NRMSE 0.739E-01 0.737E-01 0.943E-01 0.124E+00 0.124E+00 0.124E+00 0.115E+00 0.155E+00 0.143E+00 0.110E+00 0.111E+00 0.137E+00 0.405E+00 0.424E+00 0.546E+00 0.723E+00
NMAE 0.407E-01 0.406E-01 0.505E-01 0.665E-01 0.665E-01 0.666E-01 0.638E-01 0.880E-01 0.815E-01 0.614E-01 0.613E-01 0.748E-01 0.226E+00 0.238E+00 0.339E+00 0.493E+00
Department of Mathematics, University of Bologna, Italy e{mail: [email protected], [email protected] ABSTRACT In this work we solve inverse problems coming from the area of Computed Tomography by means of regularization methods based on conjugate gradient iterations. We develop a stopping criterion which is ecient for the computation of a regularized solution for the least{squares normal equations. The stopping rule can be suitably applied also to the Tikhonov regularization method. We report computational experiences based on di erent physical models and with di erent degrees of noise. We compare the results obtained with those computed by other currently used methods such as Algebraic reconstruction Techinques (ART) and Backprojection.
1. Introduction The problems arising in the Computed Tomography area are well known for their high dimensions and ill{posedness. There are di erent kinds of Computed Tomography problems (TAC, SPECT, PET), which are all modeled by a rst kind Fredholm integral equation: Z1 Z1 g(; t) = K(x; y; ; t)f(x; y)dxdy (1) ?1 ?1
In practical problems the unknown function f(x; y) represents di erent physical phenomena (attenuation in TAC, emission in SPECT and PET) and the kernel K(x; y; ; t) may have di erent expressions. In the easiest Department of Mathematics, Piazza di Porta San Donato,5 40127 Bologna (ITALY)
APPLIED MATHEMATICS AND COMPUTATION xx:1{xxx (1994) 1
c Elsevier Science Inc., 1994 655 Avenue of the Americas, New York, NY 10010 0096-3003/94/$6.00
2 models it only accounts for the geometrical characteristics of the problem, i.e.,[1]: K(x; y; ; t) = (t ? xcos() ? ysin()) (2) and the value of g(; t) represents the integral of f(x; y) along the line r of orientation + =2 at displacement t from the center of the spatial coordinate system. In more realistic models, in which other physical phenomena are considered, it represents a gaussian response to an impulse. For example, in the case of SPECT, the following kernel has been used[13]: z(x; y; ; t) 2 p K(x; y; ; t) = 1=( 2(x; y; ; t)) exp(?0:5 (x; y; ; t) ) (3) where is the standard deviation of the Gaussian system response and z is a function of the distance of the point (x; y) from the line r. It is well known that (1) is an ill{posed problem and the solution of the integral equation in the discrete domain leads to the following linear system: Kf = g (4) where the matrix K is usually neither square nor full rank. For this reason and for the presence of noise in g, the following least{squares problem is solved: min k Kf ? g k2 (5) The ill{conditioning of the matrix K makes the least{squares solution useless and requires the computation of a regularized solution. Both direct and iterative methods can be utilized [8], but iterative methods based on the conjugate gradient iterations are more ecient when applied to large size problems [14]. In this work we compute least{squares regularized solutions by means of regularization methods based on conjugate gradient iterations. In particular, in section 2 we show the properties of the Conjugate Gradient method and derive a suitable stopping criterion to compute eciently a regularized solution. In section 3 the same criterion is extended to the Tikhonov regularization method. Computational experiments are reported in section 4. 2. The Conjugate Gradient Method The matrix K of (5) is usually not full rank, severely ill{conditioned and has a few very small singular values . Although the Conjugate Gradient method applied to the least{squares normal equations (CGNE) does not necessarily converge, it has been observed that the approximations seem
3 to converge in the rst iterations, before they start oscillating and nally diverge. This phenomenon is known as semiconvergence and the best approximation obtained is called a regularized solution of the problem ([8], [9]). When the number of iterations exceeds a threshold value, the solution diverges from the regularized one and the method must be stopped. Hence the number of iterations plays the role of regularization parameter for the method. The merits of the CGNE as a regularization method depend on the determination of the optimal stopping iteration index k , that is generally related to the conditioning of K and to the speci c problem. The ill{conditioning of K can be analysed in terms of the Singular Values Decomposition (SVD) of the matrix ([12], [11]). If q singular values are comparatively large with respect to the other singular values, i.e. the problem is extremely ill{posed (its singular values decay very fast as ri with 0 < r < 1 and i > 1), then after about q iterations the solution will be in uenced by noise and will diverge. Alternatively if the singular values are almost uniformly distributed over the interval the regularized solution will be reached in a much large number of iterations. In our problems, only the last few singular values decay very fast ( gure 4) and hence the iterations must be stopped for k < q. We have determined a stopping criterion on the basis of the behaviour of the following successions computed during the CGNE iterations: k = k xk ? xk?1 k2 k = k rk k2
(6)
where xk is the approximate solution and rk is the residual at the k?th step. We compare the behaviour of these successions with that of the euclidean norm of the absolute error: E2k =k xk ? f k2. In the case of a positive de nite well conditioned matrix both the successions (6) decrease monotonically, but in the case of an ill{posed problem the successions k and k have a typical oscillating behaviour while E2k reaches a minimum and then start to diverge. We report, as an example, the plots of the singular values of K representing 16 20 projections of a 16 16 image f ( gure 4). In gures 1, 2 and 3 the successions k , k and E2k are reported for a noisy g when a kernel of the kind (2) is used, while the same plots for a kernel (3) are not reported because they are very similar. We used these observations to derive two almost equivalent stopping rules: if k > k?1 (1 < 10) then stop. if k > k?1 (1 < 4) then stop. The values of have been sperimentally determined for these kind of problems; they essentially depend on the matrix K and on the data g.
4 || r ||
3
10
2
10
1
10
0
10
−1
10
−2
10
−3
10
−4
10
0
100
200
300 iteration
400
500
600
500
600
. Succession .
Fig. 1
k
|| x_k − x_(k−1) ||
1
10
0
10
−1
10
−2
10
−3
10
−4
10
−5
10
0
100
Fig. 2
200
300 iteration
. Succession
400
k
3. Tikhonov Method Another highly regarded method for solving (5) is the Tikhonov regularization method, which corresponds to the solution of: min k K f ? g k2 + k Qf k2 where Q is a suitable regularization matrix and is the regularization parameter. The regularized solution f is then obtained by the normal equations: (K t K + Qt Q)f = K t g (7) In practical applications the regularization matrix Q is the discretization of some di erential operator which expresses the smoothness of the solution.
5 E2
1
10
0
10
0
100
200
Fig. 3
300 iteration
400
500
600
200
250
300
. E2
singular values
2
10
0
10
−2
10
−4
10
−6
10
−8
10
−10
10
−12
10
−14
10
−16
10
0
50
Fig. 4
100
150
. Singular values.
One choice commonly used in tomographic reconstruction [3] is to de ne the regularization matrix by means of the discrepancy between the value of each image point and its nearest neighbours: (8) (f ) = 21 f T(I + C tC)f where y = C f has the following expression: yj = fj ? 18 (fj ?n?1 + fj ?n + fj ?n+1 + fj ?1 + fj +1 + fj +n?1 + fj +n + fj +n+1) and is a positive weight factor. Since the matrix (I + C tC) results positive de nite we can determine Q such that: Qt Q = (I + C tC). The system (7) can be solved with direct or iterative methods. If the matrix is sparse, as is the case of kernel (2), the conjugate gradient method is more ecient, since it accounts for the sparsity of the matrix and has a
6 lower computational complexity if compared with a matrix factorization. Moreover, even in this case a suitable stopping criteria can be adopted in order to obtain a regularized solution before the convergence of the method. The behaviour of the iterations depend strongly on the regularization parameter and on the function used. When the regularization parameter is small the successions k and k behave like in the CGNE previously described, while for larger values of the regularization parameter they decrease smoothly ( gures 5, 6). We stopped the iterations of the conjugate gradient in Tikhonov method using the same CGNE criterion for small values of , while the the usual convergence criterion for conjugate gradients is used when is large and the problem conditioning is far better. For the computation of the optimal regularization parameter there are several proposals in literature, such as the L{curve method, the Generalized Cross Validation method or the discrepancy principle. They usually requires the solutions f for di erent values of and a criterion to select the best solution obtained, hence they are computationally very expensive and some counterexamples show that the methods are not always precises. 4. Numerical Results In this section we give some numerical results using simulations of brain sections (phantoms). We used the following test images: P1 95 95 samples at uniform distance 1:6 obtained by using snark93 test problem b.7 ( gure 7). P2 32 32 simulated brain section (Ho man phantom) ( gure 10). The projection data are obtained by applying (4) where K has the values given in (2) for problem P1 and (3) for problem P2. When an additive Gaussian white noise is introduced in the projection data it is measured by means of the noise level parameter: NL = k gk ?g kg~ k2 2 where g~ is the noisy signal. The pictures ~f reconstructed with di erent methods are compared with the original images f by means of the following error parameters [1]: v u P P jf ? ~f j P P (f ? ~f )2 u u =1 =1 =1 =1 NRMSE = u u t P P (f ? f )2 NMAE = P P jf j u;v
n
v
u;v
u;v
v
u
v
n
n
n
n
n
=1 u=1
u;v
u
n
u;v
v
n
=1 u=1
u;v
7
1
10
0
10
−1
10
−2
10
0
5
10
15
20
25
30
35
40
45
Fig. 5. Euclidean norm of the residuals ( = 0:06). 3
10
2
10
1
10
0
10
−1
10
−2
10
−3
10
−4
10
−5
10
0
5
10
15
20
25
30
Fig. 6. Euclidean norm of the residuals ( = 256).
35
40
8 where f is the average of the values in the original image f. In order to evaluate the reconstruction quality we report in table 1 the results obtained for P1 with four methods: ART, Backprojection, CGNE and Tikhonov methods. We observe that the errors increase when the number of projections decrease, and in particular ART and Backprojection are more sensitive to the number of projections. Noise also increases the errors, expecially for ART and Backprojection, while CGNE and Tikhonov (TIKH) methods give better results in more unfavourable situations (high noise, few projections). In table 2 we report the results obtained with CGNE algorithm for di erent values of the parameter . By increasing the value of the number of iterations increases while the errors reach a minimum value and then start to increase (as shown in section 2). The results obtained with Tikhonov method are reported in table 3 with di erent values of the regularization parameter . We observe that the number of iterations is greater with respect to CGNE and in general each iteration is computationally more expensive. In the case of test P1 the errors are a little bit better with respect to CGNE but the number of iterations increases signi cantly. Some reconstructed images are shown in gures 8, 9, 11 and 12. 5. Conclusions We can conclude that CGNE is the preferable regularization method for the problems in the area of Computed Tomography not only for the quality of the results obtained but also for its eciency (number of iterations) and the easy parameter setting. In fact, while the setting of the regularization Tikhonov parameter is dicult and costly (section 3), only a few values of
are sucient to obtain the best results with CGNE. REFERENCES 1 G.T Herman, Image Reconstruction from Projections, Academic Press, 1980. 2 F. Natterer, The Mathematics of Computerized Tomography, John Wiley & Sons 1989. 3 R.L Kayshap and C. Mittal, IEEE Trans. on Computers 24, 915 (1975) 4 P.C.Hansen, Numerical Algorithms 6, 1 (1994) 5 J.M.Varah, SIAM J. Numer. Anal. 10, 257 (1973)
9
7. 95 95 Phantom
Fig.
Fig. 8. = 50 Backprojection. k
9. 50proj: CGNE 25 it.
Fig.
10
10. 32 32 Ho man Phantom
Fig.
Fig. 11. 50proj: CGNE 29 it.
Fig. 12. 50proj: Tikhonov
11 6 7 8 9 10 11 12 13 14
P.C.Hansen, SIAM J. Sci. Stat. Comp. 11, 503 (1990) P.C.Hansen, BIT 27, 543 (1987) M. Hanke, P.C.Hansen, Surveys on Mathematics for Industry 3, 543 (1993) M. Hanke Pitman Reserch Notes in Mathematics , 327 (1995) C. C. Paige, M. A. Saunders,Anal Comp Math 8, 43 (1982) R. Plato, G. Vainikko, Numer. Funct. Anal. Optim. 11, 111 (1990) A. van der Sluis, H. A. Van der Vorst, Linear Algebra Appl. 130, 257 (1990) Formiconi et alt. Phys. Med. Biol. 37, (1992) E. Loli Piccolomini, F. Zama, Regularization Algorithms for Image Reconstruction from Projections, Proc. of the Conf. "Advanced Mathematical Tools in Metrology III", Berlin, 1996.
12
. Test Image 95 95.
Table 1
Proj. NL 100 0.
50
25
2. 2. 2. 2. 0. 0. 0. 0. 2. 2. 2. 2. 0. 0. 0. 0. 2. 2. 2. 2.
Method NRMSE NMAE Iter ART 0.61 0.42 11 BPRJ 0.14 9.3(-2) CGNE 0.97(-1) 0.52(-1) 16 TIKH(1.e-3) 0.60(-1) 0.1 18 ART 0.88 0.59 96 BPRJ 0.63 0.42 CGNE 0.29 0.19 9 TIKH (1.e+2) 0.22 0.14 17 ART 0.44 0.29 32 BPRJ 0.24 0.16 CGNE 0.15 9.1(-2) 9 TIKH(1.e-3) 0.15 9.1(-2) 9 ART 0.84 0.55 28 BPRJ 0.74 0.49 CGNE 0.33 0.33 11 TIKH(0.25e+2) 0.22 0.13 11 ART 0.46 0.31 25 BPRJ 0.42 0.27 CGNE 0.21 0.13 9 TIKH(1.e-3) 0.21 0.13 11 ART 0.85 0.56 18 BPRJ 0.83 0.55 CGNE 0.34 0.21 9 TIKH(0.25e+2) 0.26 0.18 11
13
. Conjugate Gradient Iterations.
Table 2
Test proj P1 100
NL 0
P1
100 5.9e-2
P1
100
0.59
P2
50
0
P2
50
5.0e-3
P2
50
5.0e-2 5.0e-2
P2
50
0.5
0.10E+01 0.22E+01 0.72E+01 0.10E+02 0.10E+01 0.32E+01 0.52E+01 0.10E+02 0.10E+01 0.12E+01 0.32E+01 0.12E+01 0.22E+01 0.32E+01 0.12E+01 0.26E+01 0.32E+01 0.12E+01 0.22E+01 0.32E+01 0.12E+01 0.22E+01 0.32E+01
iter 9 16 84 100 9 29 81 100 9 16 100 5 29 100 5 73 99 5 29 95 5 33 127
Emax 0.733E+02 0.650E+02 0.434E+02 0.398E+02 0.728E+02 0.581E+02 0.444E+02 0.389E+02 0.769E+02 0.748E+02 0.183E+03 0.120E+03 0.112E+03 0.109E+03 0.120E+03 0.109E+03 0.109E+03 0.120E+03 0.111E+03 0.108E+03 0.120E+03 0.107E+03 0.180E+03
NRMSE 0.117E+00 0.970E-01 0.660E-01 0.637E-01 0.117E+00 0.844E-01 0.731E-01 0.738E-01 0.124E+00 0.115E+00 0.466E+00 0.681E+00 0.393E+00 0.219E+00 0.681E+00 0.256E+00 0.220E+00 0.681E+00 0.393E+00 0.228E+00 0.681E+00 0.404E+00 0.555E+00
NMAE 0.615E-01 0.519E-01 0.362E-01 0.349E-01 0.615E-01 0.454E-01 0.405E-01 0.412E-01 0.665E-01 0.638E-01 0.262E+00 0.453E+00 0.203E+00 0.884E-01 0.453E+00 0.112E+00 0.891E-01 0.453E+00 0.203E+00 0.966E-01 0.453E+00 0.227E+00 0.278E+00
14
. Tikhonov Method
Table 3
Test proj. NL P1 100 5.9e-2 0.10E-02 0.10E-01 0.10E+00 P1 100 0.59 0.10E-02 0.10E-01 0.10E+00 0.10E-02 0.10E+00 0.10E+01 0.10E+02 0.15E+02 0.50E+02 P2 50 0.5 0.10E-01 0.10E+00 0.10E+01 0.10E+02
0.32E+01 0.32E+01 0.32E+01 0.10E+01 0.10E+01 0.10E+01 0.32E+01 0.32E+01 0.32E+01 0.32E+01 0.32E+01 0.32E+01 0.22E+01 0.22E+01 0.32E+01 0.32E+01
iter 68 68 17 9 9 9 16 29 30 69 54 29 33 33 56 23
Emax 0.492E+02 0.491E+02 0.628E+02 0.769E+02 0.769E+02 0.769E+02 0.748E+02 0.735E+02 0.717E+02 0.676E+02 0.681E+02 0.715E+02 0.107E+03 0.103E+03 0.106E+03 0.117E+03
NRMSE 0.739E-01 0.737E-01 0.943E-01 0.124E+00 0.124E+00 0.124E+00 0.115E+00 0.155E+00 0.143E+00 0.110E+00 0.111E+00 0.137E+00 0.405E+00 0.424E+00 0.546E+00 0.723E+00
NMAE 0.407E-01 0.406E-01 0.505E-01 0.665E-01 0.665E-01 0.666E-01 0.638E-01 0.880E-01 0.815E-01 0.614E-01 0.613E-01 0.748E-01 0.226E+00 0.238E+00 0.339E+00 0.493E+00
