Dawson, Christopher A

Hsieh, Jiang

Johnson, Roger H.

Molthen, Robert C.

English

epublications@Marquette

Marquette Universitye-Publications@MarquetteBiomedical Engineering Faculty Research andPublications Biomedical Engineering, Department of1-1-1999An Iterative Approach to the Beam HardeningCorrection in Cone Beam CT (Proceedings)Jiang HsiehGE Medical SystemsRobert C. MolthenMarquette University, robert.molthen@marquette.eduChristopher A. DawsonMedical College of WisconsinRoger H. JohnsonPublished version. Published as part of the proceedings of the conference, 1999 International Meetingon Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine, 1999: 85-95.Publisher Link. © 1999 Fully 3-D. Used with permission.88 1999 International Meeting on Fully Three-Dimensional Image Reconstruction In Radiology and Nuclear Medicine AN ITERATIVE APPROACH TO THE BEAM HARDENING CORRECTION IN CONE BEAM CT Jiang Hsieh l, Robert C. Molthcnl.3, Christopher A. Dawson3, Roger H. Johnsonv I Applied Science Laboratory, OE Medical Systems, Milwaukee, WI 53201 ZDepmtmcnt of Diomedical Engineering, Marquette University, Milwaukee, WI 53201 'Zablocki VA Medical Center, Milwaukee, WI 53201 Introduction It is well known that most materials preferentially absorb low energy x-ray photons much more than the high energy photons. For most of the commercially available CT scrumers, the output of the x-ray tube exhibits a very broad energy spectrum. As a result, the measured projections. deviate from the true line integral of the object, which often leads to artifacts in the reconstructed images [1,2]. In the past, many studies have been performed to understand the cause of these artifacts; and vru.-ious correction schemes have been proposed to combat the artifacts. For example, when an object is made of a single material, a polynomial correction can be applied to the detected signal to arrive at a modifIed representation based on the incident x-ray beam spectrum and the attenuation characteristics of the material [2]. Although this method performs well in most cases, it fails to remove shading artifacts caused by multiple high density materials inside the scanned field of view. An alternative approach is to pre-filter the spectrum of the x-ray beam so that the low energy x-ray photons are removed from the output of the x-ray tube [3]. This filtering approach can only achieve a compromise between the low contrast detectability and the reduction of image artifacts. To overcome these difficulties, an iterative method was proposed to estimate and correct for the projection error in the reconstructed images [4]. This method has been shown to be effective in dealing with image artifact in a fan beam CT system. The trade-off, however, is that two fan beam reconstruction processes need to be perfonned to produce a single image. For a cone beam CT geometry, the computational complexity becomes even more significant. In this paper, we present an iterative beam hardening correction method for cone beam CT. We make use of a tilted parallel beam geometry for an approximation of the beam hardening effects in cone beam geometry. \Ve show that the amount of error introduced by this approximation is negligible. Following a similar derivation of the Feldkamp reconstruction formula, we present a closed form solution for reconstruction in the tilted parallel beam sampling geometry. In the last section, we present a set of experiments to demonstrate the effectiveness of our approach. Error due to beam hardening To fully understand the correction algorithm, a brief discussion of beam hardening artifact is in order. Let us denote by 10 the incident x-ray beam intensity, h the transmitted beam intensity, and PL the line integral of attenuation coefficient along path L. For a homogeneous object, the following relationship exist: N (J j 10 PL:=:ILailn -j::ol I L (1) In this equation, N represents the order of polynomial to model beam hardening effect. (Xi depends not only on the nature of the object material, but also on the incident x-ray photon energy spectrum. Under monoenergetic x-ray beam conditions, a j == 0, Vi ¢ 1. When two different materials are present (as is often the case for typical CT scans), the above relationship is no longer valid. If we denote the transmitted beam intensity after the fIrst (low density) material by Is, the true line integral estimation can be represented by: p L ~ t. a, In '( ~: J + t. Po In' U J (2) Similar to Uj, ~i depends on the second material and the x-ray beam spectrum. If we treat the object as a single material in the beam hardening correction, the error in the projection, M L, will be: N , ML == I (ai - f3z r' == (al -.Bt t+ (a2 - P2 ~2+.. (3) i=l where 5' == In(J s / J L ). The logarithm terms in this equation are closely related to the line integral of the second (high density) material. In fact, to a first order approximation, it is proportional to second material's line integral. The linear error term in this equation causes a low frequency CT number shift (close to a DC error in the reconstructed images). The nonlinear terms in this equation are responsible for the production of the streaking and shading artifacts. When ignoring all the higher order tenus, we could conclude that the error introduced in the projection is roughly proportional to the line integral of the high density object along the x-ray beam path. The nonlinear beam hardening artifact can be illustrated by the following experiment. We scanned an excised rat lung with a micro-focal cone beam CT scanner. The set of projections were processed and reconstructed using the 1999 International Meeting on Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine 89 cone beam reconstruction algorithm as described in reference [5]. For the data acquisition, there are 480 detectors in a row and 480 rows in each projection. For the cone beam sampling, 360 projections were acquired. The data were reconstructed on a 457x457x457 matrix. To separate the effect of the potential artifacts introduced by cone beam reconstruction, we selected a reconstructed slice that is close to the center cone plane (at this location, the projection data set is essentially in a fan beam geometry). The resulting image is shown in Fig. 1 (a). Note that dark streaks are clearly visible in the regions that connect several denser objects. (a) (b) Fig. 1 Reconstructed rat lung image at slice location 199 (display window width=1500, window level=100) (a) original image (b) corrected image Iterative correction algorithm Based on Eq. (3), the logarithm term in the projection error can be estimated by forward projecting the reconstructed images containing only the high density objects along the same path as the data collection. This would require a cone beam forward projection operation followed by a cone beam backprojection reconstruction process to produce the corrected images. This process is quite computational intensive. To overcome this difficulty, we propose the projection error estimation based on a tilted parallel beam geometry. Instead of modifying the original projections for the beam hardening error, the error will be removed in image space by reconstructing the projection error derived from the tilted parallel beam geometry. To analyze the error caused by approximating the cone beam by the tilted parallel beam geometry in the beam hardening error estimation, let us define the coordinates for both systems. Let us consider a point, P, inside the scan field in which the beam hardening error needs to be estimated. There is a unique ray for each cone beam projection that intersects P. We will use two angles to define such a ray. The first angle, ~, is the angle in the x-y plane formed with the y-axis by the plane passing through P, containing the source, and parallel to the z-axis, as shown in Fig.2(a). The second angle, S, is the angle formed with the x-y plane by the ray passing through P. Note that in this notation, ~ equals the projection angle only at the iso-ray (ray intersecting z-axis). Therefore, for the same projection angle a., the ray that interests P will be significantly different in the tilted parallel beam (Fig. 2(b)) than in cone beam geometry. As a result, the estimated beam path through the dense object will be in error. However, if we compare the parallel projection whose projection angle is identical to ~ (since many views are generated, we can always find a view that eithe~ satisfies or nearly satisfies this condition), the difference in the parallel ray and the cone beam ray will be the tilt angle, S. Since the geometry is rotational symmetric with respect to the z-axis, we could simplify our analysis by selecting P to be a point in the x-z plane. Let us denote by Scone and Stilt as the tilt angle for cone beam and the parallel beam, respectively. It can be shown that the difference tilt angle, .6.S, for a point, P(x, z), can be expressed by the following equation: .6.0 = Ocone - Btilt = tan - tan I -1 ( Z J -1 ( z 'I ~ x2 + R 2 _ 2xR cos(j3')- x2 sin 2 (.s') ~ x2 + R2 - 2xR cos(.s') ) (4) 90 1999 International Meeting on Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine where B '=B-rc/2 and R is the x-ray source to iso-center distance. In general, this error increases with an increase in x and z. To estimate the impact on a cone beam system, we calculate the average error (over all ~), .6.8.v ' for R=200mm, and let x and z vary from Orum to 30l11m (which represents a cone angle of 16.7°). The maximum average en'or for this case is less than 0.019°, which clearly justifies the approximation. y y (a) (b) . Fig. 2 different sampling geometry (a): cone beam (b) tilted parallel beam Once the projection error is generated by a set of tilted parallel beam projections from only the high density images, the error in the reconstructed image can be estimated by reconstructing the error projections. The derivation for the tilted parallel beam reconstruction fonnula can be carried out in a similar fashion as the derivation used in the Feldkamp cone beam reconstruction algorithm. The resulting equation for the estimated error, e(x, y, z), is: e(x, y, z) = r .. ..j d [f" S p({J),Z)e i21{{JX uxl{J)]df3 (5) o d 2 +Z2 0 where Sp(fJ),Z) = f: Pp (t,Z)e- i2troX dt In this equation, P~(t, Z) is the projection intersecting the point (x, y, z). Unlike the parallel beam reconstruction formula, the integration limits are from 0 to 2n. This results primarily from the fact that opposing projections (separated by 180°) no longer contain complementary rays for the tilted parallel case. The fmal reconstructed image is then simply the summation of the fIrst pass reconstruction, q(x, y, z), with the scaled error image: f(x,y,z) =q(x,y,z)+A·e(x,y,z) (6) Experimental results and discussion To investigate the effectiveness of our beam hardening correction approach, we applied the correction to the cone beam data set shown previously. Figure l(b) shows the same scan reconstructed with the iterative beam hardening correction algorithm. It is evident that majority of the shading artifacts have been removed or signifIcantly suppressed. . In the derivation of our correction algorithm, we have utilized a tilted parallel beam approximation for the cone beam geometry. It can be shown that the error, .6.9, increases with an increase in z. To demonstrate the accuracy of our approximation, we performed the same correction on a slice that is near the top of the field of view in z (slice number 90). The original cone beam reconstruction and the beam hardening corrected images are shown in Figs. 3(a) and (b), respectively. Again, the effectiveness of the correction algorithm is clearly demonstrated. When comparing the correction performed on the central slice (Fig. 1) with that of a peripheral slice (Fig. 3), no clear degradation can be observed due to the modified tilted parallel beam geometry. The same conclusion can be drawn from the sagittal images generated from the reconstructed volume, as shown in Fig. 4. We would like to point out some considerations in the selection of several parameters. The first is the number of views used for the projection error estimation. Note that Eq. (4) is based on the assumption that the projection angle, ~, can be matched or made nearly identical to the angle between any cone beam ray and a titled parallel projection. This condition can only be met when the number of projections generated for the error estimation is large. If the condition can not be met, additional error for P has to be added to the error estimation presented in Eq. (4). For the implementation of our correction algorithm, we have selected the number of views to be 500. Another consideration for the correction is the selection of the reconstruction kernel in Eq. (5). If the cutoff frequency is selected too high; the error estimation will be quite noisy and its addition to the original image will 1999 International Meeting on Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine 91 degrade the signal to noise ratio of the original image. On the other hand, if the cutoff frequency is too low, many beam hardening features will be suppressed during the reconstruction and significant residual error will remain. Close examination of the beam hardening artifacts have indicated that the majority of the streaking and shading artifacts are low frequency in nature. Conclusion In this paper, we have presented an iterative con'ection algorithm for beam hardening artifact in cone beam CT. Base on the error analysis, we use a tilted parallel beam geometry for the estimation of the cone beam projection errOr. Analysis has indicated that this approximation is adequate for this application. The advantage of this approximation is a significant reduction ill the computational complexity of error estimation and reconstruction. A reconstruction formula for the tilted parallel beam geometry is also presented in the paper. Experiments with x-ray cone beam data have confirmed the robustness and effectiveness of our approach. References [1] R. A. Brooks and G. DiChiro, "Beam hardening in reconstructive tomography," Physics in Medicine and Biology 21(3), pp. 390-398,1976. [2] J. Hsieh, "Image aJ.1ifacts, causes, and correction," in kledical CT and ultrasound: current technology and applications, Advanced Medical Publishing, Madison, 1995. [3] R. J. Jennings, "A method for comparing beam-hardening filter materials for diagnostic radiology," Med. Phys. 15(4), pp. 588-599, 1988. [4] P. M. Joseph and R. D. Spital, "A method for correcting bone induced artifacts in computed tomography scanners," J. Compo Assis. Tomgr., vol. 2, pp. 100-108, 1978. [5] R. H. Johnson, H. Hu, S. T. Haworth, P. S. Cho, C. A. Dawson, and J. H. Linehan, "Feldkamp and circ1e-and-line cone-beam reconstruction for 3D micro-CT of vascular networks," Phys. Med. BioI. 43, pp. 929-940, 1998. [6] 1. A. Feldkamp, L. C. Davis, and J. W. Kress, "Practical cone-beam algorithm," J. Opt. Soc. Am. A., vol. 1, no. 6, pp. 612-619,1984. (a) (b) Fig.3 Reconstructed rat lung image at slice location 90 (display window width=1500, window level=100) (a) original image (b) corrected image (a) (b) Fig. 4 Sagittal image of the lung scan (slice 80-slice 199) (display window width=1500, window level=lOO) (a) original image (b) corrected image 

An Iterative Approach to the Beam Hardening Correction in Cone Beam CT (Proceedings)

https://epublications.marquette.edu/cgi/viewcontent.cgi?article=1124&amp;context=bioengin_fac

An Iterative Approach to the Beam Hardening Correction in Cone Beam CT (Proceedings)

Abstract

Similar works

Full text

Available Versions

epublications@Marquette