Abstract
In whole-core CT imaging, scanned data corresponding to the central portion of a cylindrical core often suffer from photon starvation, because increasing photon flux will cause overflow on some detector units under the restriction of detector dynamic range. Either photon starvation or data overflow will lead to increased noise or severe artifacts in the reconstructed CT image. In addition, cupping shaped beam hardening artifacts also appear in the whole-core CT image. In this paper, we present a method to design an attenuator for cone beam whole-core CT, which not only reduces the dynamic range requirement for high SNR data scanning, but also corrects beam hardening artifacts. Both simulation and real data are employed to verify our design method.
© 2016 Optical Society of America
1. Introduction
Digital core analysis has greatly advanced in recent years, showing high potential to explore tight (or shale) oil and gas. Computed Tomography (CT) is one of the powerful techniques for building a digital core model, since it can reconstruct internal 3D fine structure of a core sample [1–4]. Image reconstruction and image segmentation are the two necessary processes of this technique. The noise and the artifacts in the reconstructed CT images may lead an error segmentation, which will further affect the accuracy of the digital core model.
A typical standard core is a 100 mm diameter cylindrical sample. Generally, the X-ray photons penetrating the central portion of the sample are far less than those penetrating the peripheral portion of the sample because of its high density. Therefore, the scanned data corresponding to the central portion of the sample often suffer from the photon starvation [5,6], which may degrade the SNR of the CT image. On the other hand, an increased photon flux may cause detector overflow, due to the restricted dynamic range of the detector units, which will also conduct artifacts in the reconstructed CT image. In order to scan the data with adequate SNR under the restriction that all detector units do not overflow, an attenuator is usually placed between the X-ray source and the detector. An attenuator is also called a bowtie filter, as it is usually shaped like a bow tie [7–10]. Inappropriately designed attenuator will distort the CT values of the reconstructed CT image [11]. Moreover, the X-ray photons employed in whole-core CT imaging are polychromatic, while the common reconstruction algorithms, such as FBP (Filtered Back Projection) and ART (Algebraic Reconstruction Technique), are based on monochromatic X-ray assumption. An image directly reconstructed by the mentioned algorithms from polychromatic projections will have beam hardening artifacts. The aim of this study is to design an attenuator which is conducive to increasing the SNR in CT image and alleviating beam hardening artifacts simultaneously.
Most studies on attenuator design are involved in medical CT [12–21], which focus on the X-ray flux adjustment for decreasing the requirement of the detector dynamic range, reducing the absorbed dose, alleviating the scattered artifacts, or correcting beam hardening artifacts. In 2004, Boone et al [15] designed an attenuator to deliver ideal beam flattening for 50 keV photons in 50% fibroglandular breast tissue. Another simulation study by Luck et al [16] optimized the attenuators for breast CT to create either uniform dose distribution in the breast or uniform noise distribution in the reconstructed image. In 2015, Kontson et al [17,18] designed three different bowtie filters for a cylindrical 14-cm diameter phantom with a uniform composition of 40/60 breast tissue for different purpose: to transmit the same spectral shape as breast tissue, to produce both constant X-ray intensity and the same spectral shape transmitted through the breast, and to generate a constant effective attenuation coefficient across the field-of-view of the breast. The attenuators above are stationary during data scanning. There are a few studies on dynamic attenuators, which change the dynamic range as scanning angle [19–21]. Photon counting detectors (PCDs) with energy discrimination capabilities are of many advantages for compositional analysis, contrast enhancement, and dose reduction [22]. However, when such detectors are operated under a high X-ray flux, coincident pulses distort the recorded energy spectrum. These distortions are called pulse pileup effects. To avoid pulse pileup effects, attenuators are also used to reduce the X-ray flux of photons reaching the detector [23].
In this paper, we propose an attenuator design method for whole-core CT. The main idea of our method is that the thickness of the attenuator is defined appropriately along each X-ray path so that the polychromatic projections can be equivalent to certain monochromatic ones. Therefore, the image directly reconstructed by common algorithms from the projections modulated by the attenuator is close to a monochromatic image, but other existing methods on attenuator design mainly focus on the modulation of the distribution of X-ray flux, and have not considered the polychromaticity of the X-ray beam. Besides, the attenuator effectively reduces dynamic range requirement on the detector so that the SNR of the scanning data are improved significantly. Both simulation and real data are employed to verify our design method. The preliminary results are desirable.
2. Method
In this section, we first introduce the main idea of the attenuator design and its mathematical model, and then analyze its feasibility and provide the numerical solving algorithm. The function representation of the attenuator is given at last.
2.1 Design Idea and Mathematical Model
A dedicated whole-core CT system is composed of X-ray source, detector, mechanical system, attenuator, controller, and computer system. Figure 1 shows a schematic diagram of a whole-core CT scanner. After penetrating both the core sample and the attenuator, the detector counts along the X-ray path can be expressed as
where denotes the linear attenuation coefficient of the attenuator at energy and the linear attenuation coefficient of the core sample at energy and position ; denotes the spectrum of the X-ray beam; and are the minimum and maximum photon energy of the X-ray beam, respectively; denotes the detector response, the thickness of the attenuator, and the scattered photons detected by the detector unit along the X-ray path .In medical CT, water phantom is usually used to pre-correct CT data so that the measured data are converted into line integrals through the patient. In this study, we try to directly convert the photon flux passing through the core sample into the approximation of its line integrals by compensating attenuation of a well-designed attenuator. Let be the mean value of the linear attenuation coefficient of a scanned core. Imagine a cylindrical reference sample with attenuation coefficient . If omitting scattered photons, the polychromatic projection of the reference sample can be expressed as
where and are the thicknesses of the rock core and the attenuator along the X-ray path, respectively; is the effective X-ray spectrum. It is noted that the polychromatic projection in Eq. (2) can be regarded as the negative logarithm of the weighted average of , where the weight isThe idea of the attenuator design is to determine the thickness of the attenuator for a given , so that the obtained polychromatic projection equals to certain monochromatic projection at some given energy, i.e.,
where the left-hand side of Eq. (3) is the polychromatic projection, and the right-hand side is the monochromatic projection at energy . Figure 2 shows a geometrical illustration of the Eq. (3). In Fig. 2, the dashed black and blue curve denote the polychromatic projections and , satisfying and respectively, and the red line denotes the monochromatic projection of the whole-core at energy . The Eq. (3) means that we should find for each core thickness so that the polychromatic projection curve intersects with the monochromatic projection line. In the following, we will prove the feasibility of finding such under some conditions.2.2 Feasibility and algorithm
To demonstrate the feasibility of our idea, we just need to prove the existence of the solution of Eq. (3), which is equivalent to prove that the following implicit equation has a solution .
Let be the diameter of a cylindrical core sample, be the thickness of the attenuator corresponding to . Setting , then we have
Suppose that and are monotone decreasing (note that the effect of K edge of most common materials can be ignored within the energy range of the detector). As , , are nonnegative functions, it can be proven (see the Appendix) that
where , or.So, according to the condition (5), (6), and the implicit function theorem, the solution of Eq. (4) exists, in other words, there exist a function such that for all .
According to the derivative of implicit function, we have
where So, to solve Eq. (3) is equivalent to solve the following ordinary differential equation (ODE):The ODE can be solved by Runge-Kutta method numerically. It can also be proved that is a singular point of the equation.It should be noted that, there is a certain degree of freedom on the selection of , which will affect the shape of the attenuator and the ratio of low energy photons in the X-ray beam. In practice, is also restricted by the manufacture process of the attenuator.
2.3 The function representation of the attenuator
In this subsection, we give the function representation of the attenuator in the attenuator coordinate system which is used to manufacture the attenuator. In a dedicated whole-core CT system, the cone beam formed by the X-ray focal spot and the flat panel detector rotates relatively around the core sample. The attenuator can be placed between the X-ray source and the sample, or between the sample and the detector. In this study, the attenuator is placed between the sample and the detector. For simplicity, the attenuator is manufactured with only one curved surface and the inconsistency of the detector is ignored. We establish a right handed coordinate system , where the turn-table is placed at the origin , the X-ray source focal spot is at , the center line of the cylindrical sample coincides with the axis, the left side of the attenuator is in the plane , its center is at , and the right side is a curved surface. The attenuator coordinate system is the translation of the coordinate system with relation , as shown in Fig. 1.
Let the intersection point of the X-ray path with the left plane of the attenuator be . For a given cylindrical sample, the interception length of the X-ray path and the sample can be calculated as
where,, and are the horizontal and vertical cone angles, respectively.Now we represent the curved surface of the attenuator with parameter and in the and coordinate system, respectively. Substitute Eq. (11) into the function , then we have , which denotes the thickness of the attenuator corresponding to the X-ray determined by coordinate . Thereby, in coordinate system, the intersection point of the X-ray path with the curved surface of the attenuator can be expressed as
After coordinate transform, the curved surface can be expressed in the coordinate system as follows
where ,.3. Experiment results
In this section, we use both simulated data and real data to verify whether an attenuator designed by our method can reduce the dynamic range requirement for high SNR data scanning and correct beam hardening artifacts simultaneously.
According to Eq. (2) and Eq. (3), we need to know the information of the CT system and the core sample, including X-ray effective spectrum, scanning geometrical parameters, the attenuation coefficients of attenuator material, as well as the attenuation coefficients of the reference sample.
For simplicity, all the parameters and setups we use for simulation are the same as our whole-core CT, as shown in Fig. 3. The X-ray source is YXLON Y.TU-225-D04, the flat panel detector is a Varian PS2520V, which is composed of detector units, and the pitch of each unit is 0.127 mm. The distance from the X-ray focal spot to the rotation center is 357 mm, and the distance from the X-ray source to the panel detector is 500 mm. For CT imaging, 1080 projections are scanned per rotation.
Through a series of experiments, we found that 140 kV is an optimal voltage for CT imaging on 100 mm diameter cylindrical core samples. So, we simulated the X-ray tube with tube voltage 140 kV. Figure 4 shows the spectrum which is obtained with open source software SpectrumGUI. While for real CT system, the effective spectrum is estimated by the method proposed in [24] from the projection data of some standard phantoms scanned at the 140 kV.
Generally, the rock core is composed of several components. It is known from [25] that the major components are (about 61.5%) and (about 15.1%). Their mass attenuation coefficients are similar to that of , as shown in Fig. 5. Therefore, we use as a reference sample. We also choose as the attenuator material, since it easy to be manufactured and relative cheap. The attenuation coefficients mentioned above is obtained from NIST [26].
3.1 Numerical simulation
In this subsection, we simulate an attenuator based on a cylindrical reference sample. From this simulation, we mainly evaluate: (1) the dependency on initial thickness ; (2) the reduction of the detector dynamic range; (3) correction of beam hardening.
Setting the diameter of cylindrical reference sample mm, we simulate with initial conditions mm and mm, respectively, as shown in Fig. 6. Based on the two curves, we design two attenuators. Figures 7 and 8 show the orthogonal sections of the attenuators corresponding to the two initial values. Clearly, the initial condition dominates the shape of the attenuator, as well as mean energy of the X-ray beam. Besides, choosing should also consider the manufacture technique of attenuator.
Next we analyze the effect of the attenuator on the dynamic range of the detector count. For comparison, the detection counts before and after penetrating the attenuator are calculated, and are expressed as a function of the thickness of the attenuator. Let the number of incident photon along each X-ray path be . Figure 9 shows the number of detected photons with the two different conditions. The quantitative results are presented in Table 1, where the percentage is defined by dividing the interval length of detected photon counts by , i.e., the number of incident photon along each X-ray path. These results indicate that there is significant potential of using this attenuator to reduce the dynamic range. Moreover, we can adjust the initial condition of algorithm to design the attenuator for meeting the practical dynamic range requirement.
In order to verify the effect of the attenuator on beam hardening correction, we simulate polychromatic projections with and without attenuator, and reconstruct images from the simulated projections with the FBP algorithm. Figure 10 shows the reconstructed images, where the images from left to right are reconstructed from data simulated without attenuator, with mm attenuator, and mm attenuator. Figure 11 shows the profiles of lines in Fig. 10. We can see from the images and the profiles that, beam hardening artifacts are effectively corrected. In addition, the reconstructed CT values are quite similar to their theoretical ones.
3.2 Real data experiments
The digital model and the object of the manufactured attenuator are shown in Fig. 12, which is calculated with the initial value mm.
To validate the designed attenuator, a series of cores are scanned and reconstructed, among which two cores are illustrates here. One is an artificial core made of cement and cobblestones, and the other is a sedimentary rock core. Please note that cement and cobblestone have different attenuation coefficients, so this sample can be regarded as a heterogeneous object; while the other sample is relatively homogenous. Each sample is scanned twice, with and without attenuator, respectively.
Since the samples are relatively heavy and are of high absorption, a slow scanning speed 1.667 degree/second is chosen. Without attenuator, detector units corresponding to the peripheral portions of the samples are about to overflow when the tube current reaches to 0.9 mA, while with the attenuator, all measured data of detector units are just within half of the detector’s dynamic range, even the tube current reaches to 3.5 mA, furthermore the measurement corresponding to the central potations of the samples are increased about 4 times in this situation. If the power of X-ray source available, then one can increase the tube current and decrease the integral time for each angle, to improve the scanning efficiency.
All images are reconstructed by FBP algorithm, as shown in Figs. 13-16, where subfigures (a) and (c) of Figs. 13 and 15 are two orthogonal sections reconstructed from the data scanned without attenuator, and subfigures (b) and (d) of Figs. 13 and 15 are two orthogonal sections reconstructed from the data scanned with the attenuator. Figures 14 and 16 show the profiles of the corresponding lines in Figs. 13 and 15. The SNR [27] of the region-of-interests (ROI) marked with squares in the cross section images are also calculated, which are listed in Table 2.
We can see from the images and the profiles that the images reconstructed from the data set acquired without attenuator are contaminated with severe beam hardening artifacts; on the contrary, no obvious beam hardening artifacts are shown in the CT image reconstructed from the data set scanned with attenuator. We can also see from Table 2 that the SNRs of the images reconstructed from the data scanned with attenuator are significantly increased.
4. Discussion
A rock core is heterogeneous in microscopic scale. But it is generally assumed to be homogeneous, because beam hardening artifact is an impact in macroscopic scale. In the following, we test two extreme cases to illustrate that the homogeneous assumption is reasonable. In one case, the sample contains a large hollow; in the other case, the sample contains a large copper block. The reconstruction results including the profiles are shown in Figs. 17-20. From these results, we can see that in both case, the parts similar to reference material are calibrated well, and the reconstructed CT value of the hollow part is bigger than its ideal value, while the reconstructed CT value of the copper part is smaller than its ideal value. We take the first case as an example to explain such phenomenon. Actually, for a pathpenetrating both reference material (with length) and air (with length), the measured polychromatic projection is
while the ideal polychromatic projection isAs, we have, which results in that the reconstructed CT value of the hollow part is bigger than its ideal value. It should be noted that these are two extreme cases, most rock cores can be regarded as homogenous objects.5. Conclusion
In this paper, we have proposed an attenuator design method for the cone beam dedicated whole-core CT. Compared with existing software artifact correction methods, the proposed method not only significantly reduces dynamic range requirement on the detector so that the SNR of the scanning data are improved, but also effectively corrects the beam hardening artifacts in the images reconstructed by the common algorithms. Moreover, the proposed method needs no extra computation cost. Consequently, it is quite fit for fast CT imaging on the whole-core samples.
Appendix
In this appendix, we prove that under some condition, for all.
Lemma: Let , and be non-negative integral functions on the interval , and for all , then
Proof: Let Then we have thatBy exchange the order of integral in the second term of the right side, we get the following formulaThereforeSince is non-negative on and. Then, i.e., the lemma holds.Theorem: if and are non-negative integral function and decreasing. is non-negative integral, then for all .
Proof: By the definition of , we have
Let , . Then, the derivative of is as follows:Since and are monotonically decreasing. We have for all . According to the above lemma, we haveThen . In addition, by formula (20), and hence for all . Thus for all .Funding
Financial Support of Beijing Education Committee; Project of National Natural Science Foundation of China (NSFC)(61401289).
References and links
1. R. A. Ketcham and W. D. Carlson, “Acquisition, optimization and interpretation of X-ray computed tomographic imagery: applications to the geosciences,” Comput. Geosci. 27(4), 381–400 (2001). [CrossRef]
2. O. P. Wennberg, L. Rennan, and R. Basquet, “Computed tomography scan imaging of natural open fractures in a porous rock; geometry and fluid flow,” Geophys. Prospect. 57(2), 239–249 (2009). [CrossRef]
3. V. Cnudde and M. N. Boone, “High-resolution X-ray computed tomography in geosciences: A review of the current technology and applications,” Earth Sci. Rev. 123, 1–17 (2013). [CrossRef]
4. Y. Liu, A. M. Kiss, D. H. Larsson, F. Yang, and P. Pianetta, “To get the most out of high resolution X-ray tomography: A review of the post-reconstruction analysis,” Spectrochim. Acta B At. Spectrosc. 117, 29–41 (2016). [CrossRef]
5. J. F. Barrett and N. Keat, “Artifacts in CT: Recognition and Avoidance,” Radiographics 24(6), 1679–1691 (2004). [CrossRef] [PubMed]
6. I. Mori, Y. Machida, M. Osanai, and K. Iinuma, “Photon starvation artifacts of X-ray CT: their true cause and a solution,” Radiological Phys. Technol. 6(1), 130–141 (2013). [CrossRef] [PubMed]
7. J. M. Boone, A. L. C. Kwan, K. Yang, G. W. Burkett, K. K. Lindfors, and T. R. Nelson, “Computed tomography for imaging the breast,” J. Mammary Gland Biol. Neoplasia 11(2), 103–111 (2006). [CrossRef] [PubMed]
8. N. Mail, D. J. Moseley, J. H. Siewerdsen, and D. A. Jaffray, “The influence of bowtie filtration on cone-beam CT image quality,” Med. Phys. 36(1), 22–32 (2009). [CrossRef] [PubMed]
9. G. J. Bootsma, F. Verhaegen, and D. A. Jaffray, “The effects of compensator and imaging geometry on the distribution of x-ray scatter in CBCT,” Med. Phys. 38(2), 897–914 (2011). [CrossRef] [PubMed]
10. W. A. Kalender, M. Beister, J. M. Boone, D. Kolditz, S. V. Vollmar, and M. C. C. Weigel, “High-resolution spiral CT of the breast at very low dose: concept and feasibility considerations,” Eur. Radiol. 22(1), 1–8 (2012). [CrossRef] [PubMed]
11. R. B. Benítez and R. Ning, “Development of a beam hardening filter and characterization of the spatial resolution for a cone beam CT imaging system,” Proc. SPIE 91, 391–396 (2011).
12. M. Blessing, M. S. Bhagwat, Y. Lyatskaya, J. R. Bellon, J. Hesser, and P. Zygmanski, “Kilovoltage beam model for flat panel imaging system with bow-tie filter for scatter prediction and correction,” Phys. Med. 28(2), 134–143 (2012). [CrossRef] [PubMed]
13. S. Bartolac, S. Graham, J. Siewerdsen, and D. Jaffray, “Fluence field optimization for noise and dose objectives in CT,” Med. Phys. 38(S1), S2–S17 (2011). [CrossRef] [PubMed]
14. S. E. McKenney, A. Nosratieh, D. Gelskey, K. Yang, S. Y. Huang, L. Chen, and J. M. Boone, “Experimental validation of a method characterizing bow tie filters in CT scanners using a real-time dose probe,” Med. Phys. 38(3), 1406–1415 (2011). [CrossRef] [PubMed]
15. J. M. Boone, N. Shah, and T. R. Nelson, “A comprehensive analysis of DgN(CT) coefficients for pendant-geometry cone-beam breast computed tomography,” Med. Phys. 31(2), 226–235 (2004). [CrossRef] [PubMed]
16. F. Lück, D. Kolditz, M. Hupfer, and W. A. Kalender, “Effect of shaped filter design on dose and image quality in breast CT,” Phys. Med. Biol. 58(12), 4205–4223 (2013). [CrossRef] [PubMed]
17. K. Kontson and R. J. Jennings, “Bowtie filters for dedicated breast CT: theory and computational implementation,” Med. Phys. 42(3), 1453–1462 (2015). [CrossRef] [PubMed]
18. K. Kontson and R. J. Jennings, “Bowtie filters for dedicated breast CT: Analysis of bowtie filter material selection,” Med. Phys. 42(9), 5270–5277 (2015). [CrossRef] [PubMed]
19. S. S. Hsieh and N. J. Pelc, “The feasibility of a piecewise-linear dynamic bowtie filter,” Med. Phys. 40(3), 031910 (2013). [CrossRef] [PubMed]
20. F. Liu, G. Wang, W. Cong, S. S. Hsieh, and N. J. Pelc, “Dynamic bowtie for fan-beam CT,” J. XRay Sci. Technol. 21(4), 579–590 (2013). [PubMed]
21. F. Liu, Q. Yang, W. Cong, and G. Wang, “Dynamic bowtie filter for cone-beam/multi-slice CT,” PLoS One 9(7), e103054 (2014). [CrossRef] [PubMed]
22. K. Taguchi and J. S. Iwanczyk, “Vision 20/20: Single photon counting X-ray detectors in medical imaging,” Med. Phys. 40(10), 100901 (2013). [CrossRef] [PubMed]
23. K. Taguchi, E. C. Frey, X. Wang, J. S. Iwanczyk, and W. C. Barber, “An analytical model of the effects of pulse pileup on the energy spectrum recorded by energy resolved photon counting x-ray detectors,” Med. Phys. 37(8), 3957–3969 (2010). [CrossRef] [PubMed]
24. E. Y. Sidky, L. Yu, X. Pan, Y. Zou, and M. Vannier, “A robust method of X-ray source spectrum estimation from transmission measurements: Demonstrated on computer simulated, scatter-free transmission data,” J. Appl. Phys. 97(12), 124701 (2005). [CrossRef]
25. K. H. Wedepohl, “The composition of the continental crust,” Geochim. Cosmochim. Acta 59(7), 1217–1232 (1995). [CrossRef]
26. J. H. Hubbell and S. M. Seltzer, Radiation and biomolecular physics division, PML NIST [Online]. Available: http://physics.nist.gov/xaamdi.
27. B. M. Gramer, D. Muenzel, V. Leber, A. K. von Thaden, H. Feussner, A. Schneider, M. Vembar, N. Soni, E. J. Rummeny, and A. M. Huber, “Impact of iterative reconstruction on CNR and SNR in dynamic myocardial perfusion imaging in an animal model,” Eur. Radiol. 22(12), 2654–2661 (2012). [CrossRef] [PubMed]