2s
Nuclear Ins@uments and Methods in Physics Research B 111 (1996) 337-340
__ __
k!Wil
B
Beam Interactions with Materials h Atoms
!!IJ
ELSEVIER
A fidelity check of a practical cone-beam tomography algorithm P. Munshi aT* , K.J. Machin b, S. Webb b a Nuclear and Mechanical Engineering, Indian Institute of Technology, Kanpur 208016, India b Joint Department of Physics, Institute of Cancer Research, The Royal Marsden NHS Trust. Downs Road, Sutton, Surrey SM2 5PT. England, UK Received
12 October
1995
Abstract Cone-beam tomography is currently being developed for commercial applications wherever two-dimensional X-ray tomography is being used. One popular method for reconstructing a three dimensional image of an object (from cone-beam data), is the Feldkamp-Davis-Kress (FDK) algorithm, and it is an extension of the well known convolution and backprojection (CBP) method used in the usual two dimensional cases. This study investigates the fidelity of approximations involved in the FDK algorithm. Test data is obtained by the cone-beam X-ray scanner developed at the Institute of Cancer Research (Surrey, England). The error formula for the CBP filters plays a key role in this work.
1. Introduction Direct three-dimensional nondestructive imaging has been suggested and attempted by several researchers in the past [l-8] and it has also been known as cone-beam tomography (CBT). Webb et al. [9], and Bateman et al. [lO,l l] demonstrated the applicability of the CBT concepts, and in the reconstruction process they used the inversion method [S], suggested by Feldkamp, Davis, and Kress (FDK algorithm). This method, though not exact, has demonstrated its simplicity and accuracy under certain constraints. The FDK algorithm is essentially a 3-d extension of the popular 2-d convolution and backprojection (CBP) method employed by all medical and industrial scanners. Consequently, the notion of CBP filters also gets carried through in the FDK implementation. These filters play a crucial role in the tomographic inversion process and the user is provided with a choice of filter functions [ 1,2]. The present work is an attempt towards understanding the approximate nature of the FDK algorithm [5] from the point-of-view of the user dependent filter functions. A quantification of the effect of filters is attempted following the study reported by Munshi et al. [ 12-141, and in this process the fidelity of the FDK algorithm is quantified with the help of the error theorem [12-141. Data has been collected by the CBT scanner developed at the Institute of Cancer Research (ICR) located in Surrey
* Corresponding author. Tel. +91 250260, e-mail
[email protected].
512 257243,
fax f91
(England). A brief discussion of this ICR-CBT scanner is given in Section 4 and a detailed account is given by Machin and Webb [15].
2. The filter functions Following Feldkamp et al. [5], the filter function, in the inversion algorithm is given by, g,,(Y)
=/“’
IRlexp(i2nRY)dR, -A,
g&Y ),
(1)
where, R is the Fourier frequency, Y depends explicitly on the point being reconstructed, and A, is the Fourier cut-off frequency, and it is given by, A, = 1/(2AY),
where AY is the data-ray spacing in detector-object Eq. (1) can be rewritten as, g,(Y)
=/“’
W(R)IRlexp(i2mRY)dR,
plane.
(2)
-A, where W(R) is a Fourier filter satisfying all the properties of the CBP filters for two-dimensional reconstruction [1,2]. A similar form has been adopted in the CBP error theorems [12-141. These filters have also been suggested by Webb [ 161 separately. The filter chosen for this study is the generalized Hamming filter (due to its popularity in medical/NDT situations), and it is given by,
512
0168-583X/96/$15.00 0 1996 Elsevier Science B.V. All rights reserved SSDI 0168-583X(95)01456-X
W(R)=B+(l-B)(cos(ITR/A,)),
(3)
P. Munshi
338
Table I Generalized
Hamming
Value of B
Code
IW”(O)l
0.540
h54 h75 h80 h91 h99
0.460 0.250 0.200 0.083 0.001
0.750 0.800 0.917 0.999
et al./Nucl.
Instr.
md
filters chosen for the study
Note: The equivalence (Refs. [ l2- 141).For,
(normalized)
is w.r.t.
Equivalent
to
cosine tiiter sine filter Ram-Lak filter W”(O)
as per the error theory
W(R)=B+(l-@(cos(~R/A,~)), we have, W”(O)
Fi*W( R) = ---I$-
or, normalized
W”(0)
for this study as A,
= (I-
E)(a/A,)’
x=0 = (1 - B). This normalization is fixed for all the data sets.
is convenient
where, 0.5 < B < 1.0. Table 1 gives the details of the B values, the “codes”, and the Fourier-space second derivative at the origin. We note that these filters have been chosen in earlier 2-d CBP studies involving the numerical [14] and the experimental [ 17-201 investigations of the error theory reported earlier [12,13]. The difference in Eqs. (1) and (2) leads essentially to the following notion of inherent error of reconstruction,
Meth.
in Phys. Rrs. B I I I (19961337-340
If V*f is zero, or near zero, other errors will be dominant. The Laplacian is zero for smooth regions of the object, while it does not exist for rough edges. For simulated objects, E, can be calculated from the original image, but for real objects the distribution is unknown, hence the error in reconstruction cannot be calculated directly. This fact motivates an indirect representation of error, and following Davis et al. [20], we choose *‘sharpness” for this purpose. For a single point in a plane, “sharpness” is defined as the value of the reconstructed grey-level of that point. For a general image, the sharpness parameter corresponds to the maximum grey-level (linear absorption coefficient) obtained in the reconstruction. It has been earlier reported [ 19,201 that, for a given data set, sharpness can be used as an indicator of the behaviour of E, arising due to the choice of the filter function. For this work, sharpness, along with Eq. (4) will be used as one of the parameters to represent the physical nature of the object cross section. A high value of sharpness corresponds to a low value of E, as defined by Eq. (4). Here we also mention that the present study involves 3-d reconstruction while the sharpness parameter used in the earlier work has been essentially a 2-d parameter. We know that the FDK algorithm is an approximate 3-d CBP, hence we extend the the notion of 2-d sharpness to the present 3-d case. Thus, in this work, “sharpness” is the maximum grey-level in the reconstructed 3-d object.
4. The ICR-CBT where, (r,4) is the point being reconstructed, fir,+) is the true value at CT,&), and, fir,&~) is the reconstructed value at (r,+).
3. Error
formula
It has been shown [ 12,131 that, for 2-d CBP, Et at a point (r,+) in the object cross section, is given by,
given
E,(r,dJ)
= K(W”(O))(
V2”f(r&).
(4)
where W”(0)
a*W( R)
= 7
3 R=O
V*f is the Laplacian of f,and K is a constant depending on the data-ray spacing. Eq. (4) is valid for objects having certain smoothness properties [ 121 provided the data is perfect [I]. The error 15, represents the point-wise theoretical error in reconstruction, and it is also obvious that the Laplacian of fir,41 has to exist for the predictions of the theorem to be valid. For points in the cross section, where V*f does not exist, the linearity between E, and W”(O) is disturbed.
scanner
In this section we summarise briefly the salient features of the CBT scanner as reported by Machin and Webb [15]. The CBT equipment consists of a microfocal X-ray tube (energy up to 30 kVp, current up to 0.2 mA, focal spot size < 5 pm), a rotating specimen platform and a high-resolution X-ray image intensifier coupled optically to a CCD video camera. Data acquisition and 3-d image reconstruction are performed by a desktop computer. This instrument can image samples with diameters of S-50 mm and the resulting tomograms have a geometric resolution of the order lo-100 pm. The signal-to-noise ration is better than 5: 1. Typical data collection time is 10.24 seconds per view implying approximate 1 hour for a full rotation in steps of 1”. Reconstruction time on the Sun-Sparc2 is about 2 hours of CPU for 128 X 128 X 128 voxel output.
5. Results and discussion Three small specimens have been selected for the present study. They are, I. a phantom of Plexiglas with two steel pins (29 x 29 x 29 mm),
P. Munshi et al. /Nucl. 2.
Instr. and Meth. in Phys. Res. B I11 (1996) 337-340
a rose-hip berry with seeds (22.5 X 22.5 X 22.5 mm),
and 3. a dry lady-bird (9.8 X 9.8 X 9.8 mm). These specimens were chosen since their radiological path-lengths are approximately the reciprocal of the attenuation coefficient for the radiation quality used. The detector plane discretization is 64 X 64 and the angular discretization is 2”. The lady-bird exhibits a relatively high absorbing nature, while the berry is comparatively less absorbing than the phantom. These specimens also have different contrast due to their internal structures (density distribution), thus making this combination quite useful for analyzing the effect of the filter functions. Table 1 lists the 5 filters which were selected for this study. Fig. 1 shows some sample images which are reconstructions obtained by the h99 filter. The cross sections shown are transaxial for specimen 1 and 2, and parallel to the detector plane (passing through the origin) for specimen 3. Fig. 2 shows the plot of “sharpness” for the three
specimens for the various filters. The x-axis denotes the normalized values of W”(O), the Fourier space second derivative
(at the origin, R = 0) of W(R).
We observe that, for all the three samples, the sharpness (i.e., l/E,) parameter shows an interestingly linear behaviour with respect to W”(0). This linearity is not surprising in view of the approximations involved in the FDK algorithm for small cone beam angles of less than 5”. It may be mentioned here, that some numerical experiments with synthetic data carried out by Buck [21] confirm this linearity for the FDK algorithm. One of the theorems reported earlier [12] states a precise order of E, in 2-d CBP implementation. This error has been shown to be proportional to the magnitude of W”(O). The approximate analysis [13] also details this result in a simple way without involving Sobolev spaces. We now see the significance of the FDK assumptions [5] in light of this result. Eq. (1) is essentially the bandlimiting filter (or the Ramachandran-Lakshminarayan filter) for two dimensional CBP [ 11. Modification of Eq. (1) by Eq. (2) is also identical to the existing 2-d CBP implementation [ 11. The error theorem [ 12,131 has given pleasing results with synthetic data [14]. The results with experimental data, in the absence of a reference image for error computations, are reasonably explained by the error theorem with the “sharpness” concept [17-201. If we follow the same strategy here we arrive at an interesting conclusion. Fig. 2 indicates that, for all the three samples, the sharpness decreases as W”(0) increases, as is the case with me 2-d CBP where the image becomes less sharp as W”(0) is increased (e.g. h99 is sharper than h54). Furthermore, this decrease in sharpness (or increase in error) is linear as is the case for 2-d CBP. This linearity is predicted by the 2-d error theorem [ 12,131, but the the results here are for the 3-d FDK algorithm which is not 3-d CBP but an approximation to it (e.g., for small cone beam angles).
Fig. 1. Reconstructed images for (a> phantom, (b) berry, and (c) lady-bird. The filter is h99. The cross section is transaxial for (a) and (b), and perpendicular to the detector plane for (c).
We thus arrive at the conclusion that the FDK assumptions for recovering a cross section are indeed within the acceptable limits to be called a good approximation to the
P. Munshi
340
et al./Nucl.
Instr. and Meth. in Phys. Res. B II1
1000
Fig. 2. Sharpness numbers for the all the reconstructions.
true 3-d CBP algorithm which has some practical difficulties in implementation of the data collection conditions.
Acknowledgements One of the authors (P.M.) thanks, the Alexander von Humboldt Foundation (Germany), Prof. Dr. W. Arnold (Fraunhofer Institute for Nondestructive Testing, Germany) and Prof. J.C. Elliot (London Hospital Medical College) for arranging the visit to England, during which detailed discussions took place leading to this work.
References Image Reconstruction from Projections: The Fundamentals of Computerized Tomography (Academic Press, New York, 1980). [2] F. Natterer, The Mathematics of Computerized Tomography (Wiley, New York, 1986). [l] G.T. Herman,
(1996)
337-340
[3] B.D. Smith, Comput. Bio. Med. 13 (1983) 81. [4] H.K. Tuy, SIAM J. Appl. Math. 43 (1983) 546. [5] L.A. Feldkamp, L.C. Davis and J.W. Kress, J. Opt. Sot. Am. A 1 (1984) 612. [6] D.V. Finch, SIAM J. Appl. Math. 45 (1985) 665. [7] P. Grangeat, Description of a 3D reconstruction algorithm for diverging X-ray beam, Conf. Proc. Radiology Society of America (1985). (81 F.C. Peyrin, IEEE Trans. Nucl. Sci. NS-32 (1985) 1512. [9] S. Webb, J. Sutcliffe, L. Burkinshaw and A Horsman, IEEE Trans. Med. Imag. MI-6 (1987) 67. [lo] J.E. Bateman, P. Rockett and S. Webb, J. Photographic Sci. 37 (1989) 92. [I 1] J.E. Bateman, G.R. Moss, P. Rockett and S. Webb, Nucl. Instr. and Meth. A 273 (1988) 767. [12] P. Munshi, R.K.S. Ratbore, K.S. Ram and M.S. Kalra, Inverse Problems 7 (1991) 399. [ 131 P. Munshi, NDT&E International 25 (1992) 191. [14] P. Munshi, R.K.S. Rathore, K.S. Ram and M.S. Kalra, NDT&E International 26 (1993) 235. [15] K. Machin and S. Webb, Phys. Med. Biol. 39 (1994) 1639. [16] S. Webb, Phys. Med. Biol. 27 (1982) 419. [17] P. Munshi, M. Maisl and H. Reiter, Proc. Spring meeting of the American Society for Nondestructive Testing, Nashville ( 1993) p. 53. [ 181 P. Munshi, M. Maisl and H. Reiter, Experimental aspects of the approximate errOr formula for tomographic reconstmcdon, technical report, Fraunhofer-Institute for Nondestructive Testing, Saarbruecken(1993). [ 191 P. Wells and P. Munshi, Nucl. Instr. and Meth. B 93 (1994) 87. [20] G.R. Davis, P. Munshi and J.C. Elliott, An investigation of hard biological tissues using the tomographic reconstruction error formula, J. X-ray Science and Technology (1994). in press. [21] J. Buck, Cone beam X-ray tomography, Fraunhofer-Institute of Nondestructive Testing, Saarbruecken, Germany, in progress.