Journal of X-Ray Science and Technology 22 (2014) 653–671 DOI 10.3233/XST-140452 IOS Press
653
A scatter correction method for dual-energy digital mammography: Monte Carlo simulation Kai Aia , Yanhua Gaob and Gang Yua,∗ a School
of Geosciences and Info-Physics, Central South University, Changsha, Hunan, China Imaging Center, Shaanxi Provincial People’s Hospital, Xi’an, Shaanxi, China
b Medical
Received 21 March 2014 Revised 6 June 2014 Accepted 7 July 2014 Abstract. PURPOSE: To develop a novel scatter correction method without additional patient dose for dual-energy digital mammography (DEDM) to reduce scatter’s impacts and enhance microcalcification detectability in dual-energy X-ray subtraction image. METHODS: Combining scatter radiation is lower spatial frequency component and calcifications are sparsely distributed in digital mammogram, we develop a new scatter correction strategy. First, an adaptive sampling scheme is presented to find possible noncalcification (zero calcification) pixels. Then the maximum likelihood expectation maximization (MLEM) algorithm is applied to evaluate initial scatter surface. The accurate scatter radiation of sampling pixels is obtained by solving dual-energy computational formula with zero calcification constraint and scatter surface constraint. RESULTS: After scatter correction, the scatter-to-primary ratio (SPR) of wedge phantom is reduced from ∼36.0% to ∼3.1% for low-energy (LE) image and ∼29.6% to ∼0.6% for high-energy (HE) image. For step phantom, the SPR is reduced from ∼42.1% and ∼30.3% to ∼3.9% and ∼0.9% for LE and HE image, respectively. The calcification contrast-to-noise ratio is improved by two orders of magnitudes in calcification images. CONCLUSIONS: The proposed method shows an excellent performance on scatter reduction and calcification detection. Compared with hardware based scatter correction strategy, our method need no extra exposure and is easy to implementation. Keywords: Scatter correction, dual-energy mammography, calcification, monte carlo
1. Introduction Breast cancer is a major carcinoma that affects women’s health. Intramammary microcalcifications have significant association with early breast cancer [1,2]. Therefore, the efficient detection and visualization of these calcifications in mammograms can provide a good basis for clinical diagnosis. However, conventional screen-film and single-energy digital mammography can only generate the overlapped calcification image, which could be obscured by overlying the fibroglandular structures. ∗ Corresponding author: Gang Yu, School of Geosciences and Info-Physics, Central South University, Changsha 41003, Hunan, China. Tel.: +86 13077327285; E-mail: [email protected].
c 2014 – IOS Press and the authors. All rights reserved 0895-3996/14/$27.50
654
K. Ai et al. / A scatter correction method for dual-energy digital mammography
Dual-energy digital mammography (DEDM) can overcome the inherent drawbacks of conventional single-energy mammography and is viewed as a promising method to enhance the calcification detectability. The DEDM makes use of two different X-ray spectrum to expose the same breast and generate a low-energy (LE) mammographic image and a high-energy (HE) mammographic image. The subtraction mammogram (calcification image especially) is then obtained by calculating dual-energy (DE) computational formula [3–5]. Theoretically, the DE formula can cancel unconcerned tissue, improve image contrast and increase the visualization of very subtle calcifications. The above theory is based on ideal conditions. That is to say, DE subtraction calculation needs accurate data without deteriorating effects. Unfortunately, there are many deleterious factors in practice such as noise interference, anode heel effect and scatter impact [6], which will affect the computational precision. Among them, scatter is one of the most serious effects. According to Cooper’s study, the measured scatter-to-primary ratio (SPR) of 4 cm breast phantom at 28 kVp were ∼39%, ∼35% and ∼42% for 0%, 43% and 100% glandular ratios respectively [7]. And it has little dependence on X-ray spectra range from 22 to 40 kVp [8]. Considering the size of calcification, these quantities of scatter radiation can lead to an enormous calculation error which may eliminate the calcification information in DE subtraction image and destroy the advantages of DEDM [9,10]. So we mainly take scatter correction into account in this paper. Many previous literatures have investigated the performance of different scatter correction methods in DEDM. These methodologies have good performances on scatter suppression and contrast enhancement, but they all have their own limitations. The anti-scatter grid can effectively remove most scatter radiation, but it also wipes off some valuable primary radiation, and this will lead to the reduction of signalto-noise ratio [11]. Pinhole-array interpolation method is used to fit scatter surface and then subtract estimated surface from original image [12]. However, this strategy would increase patient dose by taking an extra exposure. Bornefalk et al. utilize scanned multislit system [13] to eliminate deteriorating effects caused by scatter and other factors. This method is inconvenient to operate, because it needs to precisely align detectors and collimator slits. Saito applies air-gap [14] to achieve the same goal. Nevertheless, air-gap often requires an additional radiation dose, and leads to a reduction of the effective field of view [15]. Chen et al. present an algorithm without extra exposure, but this method easily falls into local optimal solutions and it’s a numerical instability method because of empirical choice of initial scatter values. Besides, the uniform sampling pixels are maybe calcification pixels lead to constraint errors of computational formulae [16]. The main work of this paper is to present a novel DEDM algorithmic scatter correction method. This software based technique exploits adaptive sampling scheme to establish DE equations and then the Bayesian image estimation (BIE) and maximum likelihood expectation maximization (MLEM) algorithm [17–22] is applied to estimate the initial scatter radiation. Finally, the accurate scatter value of sampling pixels is calculated based on noncalcification (zero calcification) constraint and scatter surface (scatter surface is a slowly varying surface) constraint. The whole scatter surface can be obtained by surface fitting. The following sections are organized as below. In Section 2, we introduce the basic DE theory that used in our work. The main idea of our scatter correction method and corresponding dataset are described in detail in Section 3. Section 4 shows experiments and results. At last, we present our discussions and draw a conclusion.
K. Ai et al. / A scatter correction method for dual-energy digital mammography
655
2. Theory According to the principles of X-ray imaging and Lemacks et al.’s work [23], human breast mainly have three different materials: adipose tissue (thickness ta ), glandular tissue (thickness tg ) and calcifications (thickness tc ). The total thickness of breast can be written as T = ta + tg + tc . Due to the negligible thickness (size) of calcifications, the three above variables can be rewritten as two unknowns: glandular ratio (gr = tg /T ≈ tg /(ta + tg )) and tc , where T is a known quantity because the compression device is used in mammography. The incident photon fluence I0 (E) (per unit energy, E ) transmits the compressed breast and then reach detector. In clinical diagnosis, the X-ray spectrum contains different photon energy, so the transmitted signal is obtained by a mathematical integration. With all the energy levels in the high and low X-ray spectrum, the response signals can be described as Elmax Rl = Elmin I0l (E)exp[−μa (E)T − gr (μg (E) − μa (E))T − μc (E)tc ]Q(E)dE , (1) Ehmax Rh = Ehmin I0h (E)exp[−μa (E)T − gr (μg (E) − μa (E))T − μc (E)tc ]Q(E)dE where Rl and Rh are the detector signals produced by the photon spectrum of low and high energy. μa (E), μg (E) and μc (E) are linear attenuation coefficients of adipose tissue, glandular tissue and calcifications at photon energy E , respectively. Elmax and Eimin (i = l, h) are energy range corresponding to LE or HE X-ray spectrum. Q(E) is the detector response. In order to increase the dynamic ranges, two reference signals Rlref and Rhref are applied. Then logarithmic signal can be written as ⎛E ⎞ E lmax lmax Di = ln ⎝ Iref (E)Q(E)dE I0 (E) exp[−μa (E)T−gr (μg (E)−μa (E))T−μc (E)tc ]Q(E)dE ⎠ Eimin
Eimin
= ln(Riref /Ri ) , i = l, h .
(2)
In DEDM, the calcification thickness tc is generated by calculating a cubic inverse-mapping function pixel-by-pixel [24] tc = c0 + c1 Dl + c2 Dh + c3 Dl2 + c4 Dh2 + c5 Dl Dh + c6 Dl3 + c7 Dh3 ,
(3)
where c0 to c7 are constants determined by calibration data. If a certain pixel is zero calcification pixel, tc equals to zero and this pixel is regarded as noncalcification. On the contrary, it is a calcification pixel. Equations (2) and (3) only consider the primary radiation, R. The actual signal received by detector of jth pixel also contains scatter radiation (if other nuisance factors are suppressed), Sj , i.e. Rij = Rij + Sij , i = l, h .
(4)
This may lead to an incorrect thickness of calcifications during DE subtraction calculation. Therefore, to generated accurate tc at jth pixel, we have to estimate and eliminate the contribution of Sij from . Then the estimated D is written as received signal Rij ij
ˆ ij = ln(Riref j (Rij − Sij )) , i = l, h . (5) D For jth pixel, substituting Eq. (5) in Eq. (3) yields an equation with three unknowns: tcj , Slj and Shj . ˆ lj + c2 D ˆ hj + c3 D ˆ 2 + c4 D ˆ 2 + c5 D ˆ lj D ˆ hj + c6 D ˆ 3 + c7 D ˆ3 . tcj = c0 + c1 D lj hj lj hj
(6)
656
K. Ai et al. / A scatter correction method for dual-energy digital mammography
3. Materials and methods To calculate tc , Sl and Sh should be firstly estimated. However, it’s hard to know the accurate scatter surface pixel-to-pixel through a common method. Fortunately, scatter radiation is lower spatial frequency component in X-ray images [25]. Wagner et al. use a polynomial with lower-order (order 2) coefficients to represent it [26]. Kappadath and Shaw assume that scatter can be describe as a quadratic surface function [12] Sl (x, y) = (a0 + a1 x + a2 x2 )(a3 + a4 y + a5 y 2 ) , (7) Sh (x, y) = (b0 + b1 x + b2 x2 )(b3 + b4 y + b5 y 2 ) where (x, y) is spatial location, a0 to a5 and b0 to b5 are constants which can be determined by a leastsquares fit of some known scatter values. In breast cancer, calcifications are sparsely distributed [16]. That is to say, almost all the pixels are zero calcification pixels in LE and HE image. Then, at zero calcification pixels, Eq. (6) can be rewritten as ˆ lj + c2 D ˆ hj + c3 D ˆ 2 + c4 D ˆ 2 + c5 D ˆ lj D ˆ hj + c6 D ˆ 3 + c7 D ˆ 3 = 0. c0 + c1 D lj hj lj hj
(8)
The proposed method utilizes the Eqs (7) and (8) to estimate LE and HE scatter values. It includes three parts. The first is to sample zero calcification pixels. The second is to estimate initial scatter surfaces of LE and HE images. And the final is to calculate scatter values at sampling pixels and then generate the whole scatter surface using quadratic B-spline function fitting method. 3.1. Sampling scheme To calculate scatter value from Eq. (8), the uniform sampling is applied in LE and HE image firstly. The adaptive sampling scheme is then performed to analyze whether the sampling pixels are calcification pixels or not. 3.1.1. Uniform sampling scheme The Eq. (8) of all sampling pixels (circles in Fig. 1(a)), can build a system of nonlinear equations. Because there are two unknowns for Eq. (8) of jth sampling pixel, the scatter values, Slj and Shj , can’t be obtained directly. However, as the scatter radiation is lower spatial frequency component, it can be visualized as a slowly changing variable within a small spatial area. That is to say, one pixel’s scatter value could be expressed as the scatter of its neighboring pixels through Eq. (7) (Fig. 1(b)). Let A1 to A4 are four uniform sampling pixels, and B1 to B4 are their intermediate pixels, i.e. additional sampling pixels. Due to scatter surface constraint, i.e. Eq. (7), the scatter value of B1 (horizontal intermediate point between A1 and A2) can be described as SB1 = (SA1 + SA2 )/2 + Δ,
(9)
where Δ = −a2 ((xA1 − xA2 )/2)2 (a3 + a4 y + a5 y 2 ). It is the error term of scatter value of additional pixel. xA1 , xA2 are horizontal ordinates of SA1 and SA2 , y is vertical ordinates of SA1 , SA2 and SB1 . Others (SB2 , SB3 and SB4 ) are described similarly. The additional sampling pixel Bi is then introduced into the system of nonlinear equations, which can be written as D × C = [ 0 0 · · · 0 · · · 0 0 ]T ,
(10)
K. Ai et al. / A scatter correction method for dual-energy digital mammography
657
(a)
(b)
(c)
Fig. 1. (a) Uniform sampling scheme for LE or HE image (circles). (b) Additional sampling pixels (black squares) in dashed box from Fig. 1(a). A1 to A4 are uniform sampling pixels, and B1 to B4 are additional selected pixels of which the scatter values can be represented by A1 to A4. (c) Adaptive sampling on a n-by-n mask. There are two cases in our adaptive scheme. For a certain sampling pixel j, if the location of pixel j is different from k, then replaces Slj and Shj with Slk and Shk (left), else remain unchanged (right).
where D = [ DA1 DA2 · · · DAn DB1 DB2 · · · DBn ]T , and
T DAi = 1 Dl (SAi ) Dh (SAi ) Dl2 (SAi ) Dh2 (SAi ) Dl Dh (SAi ) Dl3 (SAi ) Dh3 (SAi ) , i = 1, 2, · · · , n ,
T DBi = 1 Dl (SBi ) Dh (SBi ) Dl2 (SBi ) Dh2 (SBi ) Dl Dh (SBi ) Dl3 (SBi ) Dh3 (SBi ) , i = 1, 2, · · · , n , C = [ c0 c1 c2 c3 c4 c5 c6 c7 ]T , n is the number of sampling pixels.
The Eq. (10) is the system of nonlinear equations based on zero calcification constraint and scatter surface constraint. 3.1.2. Adaptive sampling scheme Obviously, the unknowns are equal to the number of equations, and the Slj and Shj can be obtained by solving Eq. (10). However, the above uniform sampling scheme has an inherent drawback. Although the probability of the appearance of calcifications is rather small in early breast cancer mammogram, we can’t make sure that all sampling pixels’ tc in nonlinear equations are equal to zeros (zero calcifications). This deficiency may lead to calculation error in DE subtraction image. Therefore, for each pixel in A sets and B sets, it should be further checked whether it is zero calcification or not. For a certain sampling pixel j , we search for its neighboring pixels with n-by-n mask. Because the calcification pixels are sparsely distributed and exist in isolation, most pixels are zero calcifications under the mask. One pixel k that has most neighbors with similar signal intensity should be considered as zero
658
K. Ai et al. / A scatter correction method for dual-energy digital mammography
calcification. As scatter surface varies slowly, Slj ≈ Slk and Shj ≈ Shk . Therefore, Eq. (10) is also workable when pixel k replaces pixel j . Obviously, if the location of pixel k doesn’t equal to that of pixel j , the pixel k should be used to generate one equation in Eq. (10) instead of pixel j , else pixel j is used. This method is so called adaptive sampling scheme proposed by us (Fig. 1(c)). 3.2. Initial scatter surface estimation As known to all, the convergence of a nonlinear equation highly depends on the choice of initial value [27] and an inappropriate starting point may get an error result (locally optimal solution) or leads to no convergence. The initial value of a0 to a5 and b0 to b5 in Eq. (7) are also needed to be known to fix Δ in Eq. (9). BIE and MLEM algorithm are utilized to estimate initial scatter. The general BIE used in our image estimation can be described as ln f (Ri |Ri ) = ln f (Ri |Ri ) + ln f (Ri ) − ln f (Ri ) , i = l, h ,
(11)
where f (Ri |Ri ) represents a posterior probability of the unknown primary radiation Ri , f (Ri |Ri ) represents the probability related to likelihood function, f (Ri ) is a priori probability of the unknown Ri and f (Ri ) is a constant. According to MLEM algorithm, Eq. (11) can be rewritten as ˆ i = arg max ln f (Ri |R ) = arg max(ln f (R |Ri ) + ln f (Ri )) , i = l, h , R i i Ri
Ri
ˆ i is the estimation of Ri . To calculate the R ˆ i , we use iterative formula presented by [17] where R ˆ (n) · R (T ∗ ∗R ˆ (n) ) , i = l, h , ˆ (n+1) = R R i i i i
(12)
(13)
where, T represents the scatter kernel K (a 2-D exponential function with two parameters: full-width at half-maximum and magnitude) plus a Dirac delta functional, i.e. T = δ + K . The denominator ˆ (0) can be easily set to R . Finally, the estimated scatter represents a convolution computation. And R i i ˆ surface is derived by subtract Ri from Ri . 3.3. Scatter values calculation and surface fitting The accurate scatter values at adaptive sampling points will be calculated by solving system of nonlinear Eq. (10). Finally, the whole accuracy scatter surface is produced by fitting method. In this paper, a quadratic B-spline function [28,29] is used. After subtract the accuracy scatter surfaces from original images, we can obtain a precisely DE calcification image. 3.4. Algorithm process The presented software based scatter correction method for DEDM consists of six successive steps. The brief overview of these steps list as follows: Step 1. Initial scatter surface estimation: The BIE and MLEM algorithm is used to estimate initial Sl and Sh . And then by using least-squares approximation, a0 to a5 and b0 to b5 in Eq. (7) can be obtained through Sl and Sh . Step 2. Adaptive sampling scheme: The uniform sampling scheme is used to sample some pixels. The n-by-n mask is then applied for adaptive sampling to avoid sampling calcification pixels.
K. Ai et al. / A scatter correction method for dual-energy digital mammography
659
Fig. 2. The flow chart of the scatter correction method.
Step 3. Establishment of system of nonlinear equations: We have made sure that the Eq. (8) is workable for each sampling pixel and additional pixel. Therefore, the system of nonlinear equations is built with the unknowns equal to the number of equations just like Eq. (10). Step 4. Estimation for initial solutions: Due to the BIE and MLEM algorithm, the initial scatter values of sampled pixels, Slj and Shj , are gained from Sl and Sh . And the initial Δ value for each additional pixel can also be fixed. old (for the first iteration, S old = Step 5. Equations solving: By using the known values Sljold and Shj lj old = S ), a new group of S new and S new can be obtained by solving Eq. (10) Slj and Shj hj lj hj iteratively.
660
K. Ai et al. / A scatter correction method for dual-energy digital mammography Table 1 Data of adipose and glandular for human breast-tissue-equivalent materials [30] Materials CIRS adipose CIRS glandular
(a)
Density (g/cm3 ) 0.924 1.04
H 11.8 10.9
Elemental compositions (%) C N O Cl Ca 76.0 1.2 9.8 1.2 − 70.2 1.2 12.5 1.1 0.6
Al − 3.5
(b)
Fig. 3. Schematic drawing of designed phantoms. (a) The lateral view of wedge phantom. (b) The lateral view of step phantom. The positive direction of y-axis pointed into paper and the x-axis represented chest wall side. new old new old Step 6. Convergence determination: If the difference between Slj and Slj , Shj and Shj don’t new old ξ , a to meet a certain criterion ξ (a constant), e.g. Sljnew − Sljold ξ and Shj − Shj 0 new and S new , as well as the Δ value. Then a5 and b0 to b5 in Eq. (7) are reinitialized by Slj hj turn back to Step 5. new old ξ or the iteration time is enough, the algorithm is Step 7. If Sljnew − Sljold ξ and Shj − Shj convergent. Then quadratic B-spline function is utilized to generate accurate scatter surface Sl new . And then S and S is subtracted from R and R , we can obtain and Sh by Sljnew and Shj l h l h a precisely primary radiation Rl and Rh . The flow chart of the algorithm is shown in Fig. 2.
3.5. Phantom design A cuboid object with breast-tissue-equivalent materials was used to simulate breast phantom. The length, width and height of the cuboid were 9 cm, 4 cm and 5 cm, respectively. The elemental compositions data (Table 1) of the adipose and glandular tissue-equivalent materials were provided by CIRS (Computerized Imaging Reference Systems, Inc., Norfolk, VA). In our experiment, two phantoms with different inner structures were designed. The first wedge phantom was used to simulate the slowly varying adipose-glandular ratio in breast. The second step phantom was used to simulate the rapid varying adipose-glandular ratio. Usually, 5 cm was the common thickness of compressed breast and the microcalcifications in breast were mainly composed of CaCO3 and reached to several hundred microns. So we set the thickness of the phantoms to 5 cm and used CaCO3 cubes to simulate calcifications. The wedge phantom was combined by two reversed wedges (Fig. 3(a)) with glandular ratio ranging from 0% to 100%, and the step phantom was made up of nine steps (Fig. 3(b)). The length of each
K. Ai et al. / A scatter correction method for dual-energy digital mammography
661
(a)
(b) Fig. 4. (a) The HE step phantom primary image. The calcifications size varied along the vertical direction and stayed fixed along the horizontal direction. And glandular ratio ranged from 10% to 90% from left to right. (b) The HE step phantom scatter image.
step was 1 cm and the height was 0.5 cm, i.e. the glandular ratio ranged from 10% to 90% with each step interval is 10%. And then three different sizes of CaCO3 cubes (200 μm, 300 μm, 400 μm and the mass density is 2.93g/cm3) were placed on the middle top of two phantoms to simulate calcifications. Calcifications were arranged to a 3×3 pattern with the size varying along the y-axis and stayed fixed along the x-axis. All the geometries above were generated by egspp [31] (a C++ class library provided by EGSnrc package [32,33]). 3.6. Image acquisition Monte Carlo (MC) method has been widely used for simulating x-ray scatter radiation and X-ray imaging [34–36]. And the accuracy of many MC simulation packages (e.g. EGSnrc, GEANT and MCNP) has been proved. In our experiment, we developed a new EGSnrc C++ user code to generate LE and HE primary and scatter data based on Lippuner et al.’s work [37]. Besides, we rewrote the phantom design code to generate wedge phantom and step phantom. During image acquisition, incident X-ray photons irradiated the designed phantom from top, and a virtual detector with 400 × 900 pixels (i.e. 100 × 100 μm2 for each pixel) was placed right after the bottom of phantom. For easy implementation, we made use of the monochromatic x-ray energies [38, 39] with 18 keV for LE imaging and 36 keV for HE imaging. To save simulating time, we applied Richardson-Lucy fitting algorithm [40] to process the raw scatter data. Figure 4(a) and (b) shown the HE step phantom primary and scatter image that acquired from our EGSnrc based C++ program (LE step phantom, LE and HE wedge phantoms were just similar to HE step phantom). The nine black dots were calcifications with different sizes.
662
K. Ai et al. / A scatter correction method for dual-energy digital mammography
Fig. 5. Schematic drawing of calibration phantom (top view). The dash lines were used to distinguish different glandular ratios, and the dotted squares represented 0 µm calcifications. The positive direction of z-axis pointed out of paper and the x-axis represented chest wall side.
3.7. Calibration procedure The aim of calibration procedure was to confirm coefficients in Eq. (3). To this end, we designed a calibration phantom of 5 cm thickness with five steps (the glandular ratios were 0%, 30%, 50%, 80% and 100%), and on each step we placed six different thickness of CaCO3 cubes (0, 100, 350, 500, 750, 1200 μm) orthogonal to x-axis (Fig. 5). The value of c0 to c7 were determined by least squares fit of thirty acquired data. And the reference signals, Rlref and Rhref were measured through calibration phantom with 100% adipose ratio. According to Kappadath and Shaw’s experiment [23], the mean fitting error of Eq. (6) was about 50 μm and the max fitting error was 150 μm for microcalcification thickness. As our experiment is a simulation program, the mean and max errors were suffered a further reduction. This meant the Eq. (6) is valid for modeling calcification thickness in this paper. 3.8. Evaluation metrics The presented method was used to cleanup scatter radiation. Scatter-to-primary ratio (SPR) was used to estimate the number of scatter components in original LE and HE image. In contrast, residual scatterto-primary ratio (rSPR) represented how much scatter components remained after scatter correction. It could be written as Sij SPR = , i = l, h , (14) Rij Sij − Sij rSPR = , i = l, h , (15) Rij represents the estimated scatter value, j is the pixel index. where Sij To prove our method could enhance calcifications contrast effectively, we also employed calcification contrast-to-noise ratio (cCNR) which was described as |Vc − Vb | cCNR = , (16) σb where Vb and σb are average signal intensity and standard deviation of background. Vc is average signal intensity of calcification.
K. Ai et al. / A scatter correction method for dual-energy digital mammography
(a)
(b)
(c)
(d)
663
Fig. 6. (a) to (b) The SPR image of LE and HE wedge phantom before scatter correction. (c) to (d) The rSPR image of LE and HE wedge phantom after scatter correction.
3.9. Scatter correction settings In this paper, we sampled 256 pixels, including 128 original sampling pixels and 128 additional sampling pixels. A 5-by-5 mask for adaptive sampling was proved to have a good performance in the following experiment. In Eq. (13), that two parameters of scatter kernel K were fixed to 50 pixels for full-width at halfmaximum and 0.029, 0.032 for magnitude of LE and HE scatter estimation respectively can work well. After 10 times iterations, Eq. (10) became convergent. 4. Results The LE and HE scatter surfaces are estimated by the presented method and subtracted from original LE and HE images, Rl and Rh . Then Eq. (6) is used to generate calcification image based on the corrected primary radiation. The SPR, rSPR and cCNR are calculated to evaluate the performance of our method.
664
K. Ai et al. / A scatter correction method for dual-energy digital mammography
(a)
(b)
(c)
(d)
Fig. 7. (a) to (b) The SPR image of LE and HE step phantom before scatter correction. (c) to (d) The rSPR image of LE and HE step phantom after scatter correction.
4.1. SPR and rSPR of two phantoms Figure 6(a) and (b) show the SPR of LE and HE wedge phantom of sampling pixels. Figure 6(c) and (d) show the rSPR of sampling pixels after scatter correction. The negative value of rSPR represents our method has overestimated scatter radiation. In Fig. 6 (a) and (b), the mean SPR of LE and HE wedge phantom of sampling pixels are ∼36.0% and ∼29.6% before scatter correction, respectively. By contrast, the mean absolute value of rSPR is suppressed to ∼3.1% and 0.6% after apply the proposed method. Our correction strategy has a very good performance on scatter reduction of wedge phantom. Figure 7 shows the corresponding result of step phantom. As can be seen from Fig. 7(a), the SPR of LE step phantom varies quickly due to the shape of phantom. The mean SPR of LE and HE step phantom of sampling pixels are ∼42.1% and ∼30.3% before scatter correction, respectively. And after scatter correction, the mean absolute value of rSPR becomes to ∼3.9% and ∼0.9%. That is to say, our method is also suitable for steep tissue in breast.
K. Ai et al. / A scatter correction method for dual-energy digital mammography
665
Table 2 The statistical analysis of SPR and absolute value of rSPR of sampling pixels Phantom shape Wedge (LE) Wedge (HE) Step (LE) Step (HE)
Mean 35.96 29.60 42.09 30.33
SPR (%) Max Min 40.00 27.25 33.26 21.95 53.75 26.63 34.93 21.05
Std. 2.92 2.63 6.24 3.21
Mean 3.09 0.61 3.86 0.85
|rSPR| (%) Max Min 8.21 0.02 1.63 0.01 11.05 0.15 3.02 0.00
Std. 2.23 0.44 3.04 0.83
Table 3 The statistical analysis of the absolute value of rSPR of sampling pixels using Chen et al.’s method Phantom shape Wedge (LE) Wedge (HE) Step (LE) Step (HE)
Mean 7.31 1.32 11.95 2.11
|rSPR |(%) Max Min 11.63 0.33 2.07 0.06 22.13 0.10 3.80 0.03
Std. 2.91 0.52 5.55 0.95
Mean 12.60 2.50 4.80 0.90
|rSPR| (%) Max Min 22.25 7.33 4.59 1.41 13.51 0.04 2.69 0.00
Std. 3.55 0.75 3.20 0.60
Table 4 The statistical analysis of absolute value of rSPR of fitting scatter surface Phantom shape Wedge (LE) Wedge (HE) Step (LE) Step (HE)
Mean 2.51 0.49 3.38 0.72
Fitting data |rSPR| (%) Max Min 7.91 0.00 1.58 0.00 11.14 0.00 3.17 0.00
Std. 1.63 0.32 2.96 0.78
Table 2 shows the maximum, minimum, mean and standard deviation of SPR and absolute value of rSPR of sampling pixels. In contrast, Chen et al. just use empirical method to select initial scatter estimation. It is a numerical instability method. According to the SPR data of Table 2, we empirically choose two different sets of initial scatter values (35% and 45% of the total received signal, Ri ) to test Chen et al.’s method and some statistical values are list below. The left and right half of Table 3 shows the statistical values with initial scatter estimation equals to 0.35×Ri and 0.45×Ri . Obviously, with initial scatter estimation changes, the absolute value of statistical values of rSPR varies seriously. This indicates that Chen et al.’s scatter correction method completely depends on the choice of initial scatter estimation. In practice, we have little prior information about SPR. So it is a hard work to set suitable initial scatter surfaces by empirical method. 4.2. Surface fitting To generate DE calcification image, we have to cleanup the whole scatter surface. Here the quadratic B-spline function and least-squares approximation are utilized to fit the smooth scatter surface through estimated scatter values at sampling pixels. Figure 8 shows the central row profiles of real scatter surface and fitting surface (normalized). The mean fitting relative error of LE and HE wedge phantom of central rows are ∼4.2% and ∼0.9%, respectively. For LE and HE step phantom, the values are ∼9.2% and ∼2.6%, respectively. The simulation experiment has demonstrated that the chose of function and order are entirely satisfied our expecta-
666
K. Ai et al. / A scatter correction method for dual-energy digital mammography
(a)
(b)
(c)
(d)
Fig. 8. The central row profiles of real scatter data (solid line), fitting data (dash-dot line) and error data (dotted line). (a) to (b) are wedge phantom of LE and HE image, respectively. (c) to (d) are step phantom of LE and HE image, respectively.
tion. Table 4 shows the maximum, minimum, mean and standard deviation of absolute value of rSPR of the whole fitting scatter surface. 4.3. Calcification image Figure 9 shows local calcification image of wedge and step phantom (2 × magnification for clear visual inspection) which generated by calculating Eq. (6). In Fig. 9, we can clearly see calcifications with the size of 300 μm and 400 μm. Due to the small size of 200 μm calcifications, the visibility is weaker than 300 μm and 400 μm. In these two local images, the glandular tissue is removed by DE algorithm and we can faintly see some residual signal (‘white’ background), but this interference signal hardly affects the visibility of calcifications.
K. Ai et al. / A scatter correction method for dual-energy digital mammography
(a)
667
(b)
Fig. 9. (a) The DE calcification image of wedge phantom. (b) The DE calcification image of step phantom.
(a)
(b)
Fig. 10. The cCNR of DE calcification images. (a) The cCNR of nine calcification pixels of wedge phantom. (b) The cCNR of nine calcification pixels of step phantom.
Figure 10 shows the cCNR of nine calcification pixels. We choose nine 20 × 20 region of interests (ROIs) as background to calculate the cCNR of the nine pixels with 200 μm, 300 μm, 400 μm calcifications on different glandular ratios. The ROIs are just nearby the corresponding calcifications. For wedge phantom, the mean cCNR of 200 μm calcifications raises from 6.45 to 124.31. For 300 μm and 400 μm calcifications, the cCNR raises from 17.97, 2.40 to 188.06 and 296.00, respectively. For step phantom, the mean cCNR for 200 μm, 300 μm, 400 μm calcifications raise from 5.43, 5.28, 4.85 to 131.44, 158.21 and 282.95, respectively. Obviously, the cCNR of DE calcification image is significantly improved after scatter correction. In Table 5, the thickness of nine calcification pixels after scatter correction is shown outside the parentheses and the uncorrected values are shown in parentheses (use Eq. (3) to obtain them). The calculated thickness after scatter correction is close to the real thickness. And the max errors of 200 μm, 300 μm and 400 μm calcification are only 28 μm, 65 μm and 74 μm for wedge phantom and 35 μm, 68 μm and 106 μm for step phantom.
668
K. Ai et al. / A scatter correction method for dual-energy digital mammography Table 5 The calcification thickness (µm) in DE image
Real thickness ( µm) 200 300 400
Glandular ratio of wedge phantom increases from left to right 215 (3347) 191 (3420) 172 (3450) 313 (3394) 276 (3478) 235 (3504) 474 (3202) 427 (3263) 432 (3302)
Glandular ratio of step phantom 40% 50% 60% 235 (3276) 221 (3337) 191 (3359) 294 (3315) 265 (3347) 232 (3388) 397 (3108) 350 (3150) 294 (3155)
5. Discussions As compared with single-energy mammography, DEDM can remove unwanted tissue to improve the detectibility of microcalcifications. Therefore, it has a good prospect for breast imaging. Due to the deteriorating effects of scatter radiation, however, the thickness of calcifications may become incorrect. Based on the two characteristics of DEDM (i.e. X-ray scatter is a lower spatial frequency component and calcifications are sparsely distributed), this paper present a novel software based scatter correction method. From a theoretical perspective, the presented method includes some key points. Because the noncalcification tissue plays a dominant role in breast, so the thickness of calcifications can be regarded as zero in most mammogram regions. That is to say, most pixels are satisfied with Eq. (8). By uniformly sampling some pixels in LE and HE image and substituting these pixels in Eq. (8), a system of nonlinear equations can be built. Nevertheless, this system of nonlinear equations can’t be solved, because the number of the equations is smaller than the unknowns. Fortunately, Eq. (7) is used to represent the low-frequency variation of the scatter surface and the scatter value of intermediate pixel between two sampling pixels can be expressed by Eq. (9). The introduction of these additional pixels makes the number of the equations equivalent to the unknowns and then a solvable Eq. (10) is obtained. Meanwhile, the error term of scatter value of additional pixel is also included in Eq. (9) to increase the computation accuracy. It’s worth noting that Eq. (10) based on uniform sampling scheme has an inherent shortcoming where we can’t ensure Eq. (8) is workable for all the sampling pixels and additional pixels. Thus the adaptive sampling scheme is necessarily utilized. This innovative point can ensure that all the selected pixels are noncalcification pixels. The initial solution estimation is also important for Eq. (10). If the inappropriate initial scatter values are used such as reference [16], Eq. (10) may be fall into local optimal solution or directly not convergence. BIE and MLEM algorithm is employed to generate a good initial scatter surface estimation and fix initial values of a0 to a5 and b0 to b5 at the same time. It plays an important role in the presented method. The results indicated that our method has a good performance on scatter correction. It can improve the visibility of microcalcifications with the size of ∼200 μm which is unachievable in reference [12] and other hardware based strategies. Compared with other algorithmic correction, our method is adaptive process and need no empirical set of initial scatter values. Despite all this, there are still some useful information can be get from Figs 6, 7 and Table 2. Most important among them is that the correction effectiveness for different shape of phantoms is not the same. The absolute value of rSPR of LE wedge phantom image is 3.09%, 0.76% less than step phantom. And for HE wedge phantom image, the absolute value of rSPR is 0.61%, 0.24% less than step phantom. This phenomenon may be caused by two reasons. Firstly, compared with the slowly varying glandular ratio in wedge phantom, the step phantom contains rapid change of glandular ratio. The boundary effect can constitute a handicap for MLEM algorithm to estimate initial scatter surface, which will further affect the solution of Eq. (10). In addition, as can be
K. Ai et al. / A scatter correction method for dual-energy digital mammography
669
seen from Table 2, the mean SPR of step LE and HE phantom image is greater than the corresponding wedge phantom image. Under the same simulation conditions and system parameters, the original scatter radiation of step phantom image may be more than that of wedge phantom, which reduces the performance of scatter correction. For the same phantom, the scatter correction effectiveness is also different between LE and HE image. The performance of scatter correction in HE image is much better than LE image no matter which phantom is used. For step phantom, scatter fraction reduces from 35.96% to 3.09% in LE image and 29.60% to 0.61% in HE image. For wedge phantom, the fraction reduces from 42.09% to 3.86% in LE image and 30.33% to 0.85% in HE image. The observed difference may be caused by the amount of scatter radiation in LE image. Compared with the scatter fraction in HE image, the larger scatter fraction in LE image leads to the reduction of primary fraction. This will affect complex calculation in Eq. (10) and lead to the error. As can be seen from Table 4, the error between the fitting scatter surface and real scatter surface suffers from a further reduction in both LE and HE image after a quadratic B-spline function is applied. And still this observation is in agreement with the findings in which the correction performance is better in wedge phantom image than step phantom image, and HE image than LE image. However, the calcification images shown in Fig. 10 and Table 5 demonstrate that the above difference has not serious impact on the visibility of calcifications. Besides the scatter radiation, X-ray radiation contains many other interference factors. Noise is another important effect among them. For simplicity, the presented method assumes that the noise in the received radiation has been removed or depressed by many known noise reduction methods [41]. Future work is underway to integrate both scatter and noise in received radiation and a new methodology which can fully correct the two effects will be carried out. 6. Conclusions In this paper, we present a novel software based scatter correction used for DEDM. This method can reduce scatter fraction within 4% and 1% for LE and HE image, and enhance the microcalcifications detectability at the same time. The cCNR is improved by two orders of magnitudes. Compared with hardware based scatter correction, the proposed method is easy to implementation and need no extra exposure. Our next work will take other deleterious effects into account (e.g. noise) and make a further investigation in calculation error to develop a complete, accurate and clinical DEDM processing application. Acknowledgments The authors would like to thank the following funding organizations for support: Hunan Provincial Natural Science Foundation of China (No. 11JJ6507) and China Postdoctoral Science Foundation (No. 20090461303), and the Fundamental Research Funds for the Central Universities of Central South University (No. 2013zzts251). References [1]
R.W. Powell, M.B. McSweeney and C.E. Wilson, X-ray calcifications as the only basis for breast biopsy, Ann Surg 197 (1983), 555–559.
670 [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30]
K. Ai et al. / A scatter correction method for dual-energy digital mammography A. Fandos-Morera, M. Prats-Esteve, J.M. Tura-Soteras and A. Traveria-Cros, Breast tumors: composition of microcalcifications, Radiology 169 (1988), 325–327. L.A. Lehmann, R.E. Alvarez, A. Macovski and W.R. Brody, Generalized image combination in dual KVP digital radiography, Med Phys 8 (1981), 659–667. A. Taibi, S. Fabbri, P. Baldelli, C. di Maggio, G. Gennaro, M. Marziani, A. Tuffanelli and M. Gambaccini, Dual-energy imaging in full-field digital mammography: A phantom study, Phys Med Biol 48 (2003), 1945–1956. S.C. Kappadath and C.C. Shaw, Quantitative evaluation of dual-energy digital mammography for calcification imaging, Phys Med Biol 49 (2004), 2563–2576. J.S. Maltz, B. Gangadharan, S. Bose, D.H. Hristov, B. A. Faddegon, A. Paodi and A.R. Bani-Hashemi, Algorithm for X-ray Scatter, Beam-Hardening, and Beam Profile Correction in Diagnostic (Kilovoltage) and Treatment (Megavoltage) Cone Beam CT, IEEE Trans. on Med Imaging 27 (2008), 1791–1810. V.N. Cooper III, J.M. Boone, J.A. Seibert and C.J. Pellot-Barakat, An edge spread technique for measurement of the scatter-to-primary ratio in mammography, Med Phys 27 (2000), 845–853. J.M. Boone, K.K. Lindfors, V.N. Cooper III and J.A. Seibert, Scatter/primary in mammography Comprehensive results, Med Phys 27 (2000), 2408–2416. T.B. Hamed and R.D. Speller, The effect of scattered radiation in dual-energy analysis, Phys Med Biol 40 (1995), 1619– 1632. F.C. Wagner, A. Macovski and D.G. Nishimura, Effects of Scatter in Dual-Energy Imaging: An Alternative Analysis, IEEE Trans. on Med Imaging 8 (1989), 236–244. S.Z. Shen, A.K. Bloomquist, G.E. Mawdsley and M.J. Yaffe, Effect of scatter and an antiscatter grid on the performance of a slot-scanning digital mammography system, Med Phys 33 (2006), 1108–1115. S.C. Kappadath and C.C. Shaw, Dual-energy digital mammography for calcification imaging: Scatter and nonuniformity corrections, Med Phys 32 (2005), 3395–3408. H. Bornefalk, M. Hemmendorff and T. Hjärn, Contrast-enhanced dual-energy mammography using a scanned multislit system: evaluation of a differential beam filtering technique, Journal of Electronic Imaging 16 (2007), 0230061– 0230067. M. Saito, Dual-energy approach to contrast-enhanced mammography using the balanced filter method: Spectral optimization and preliminary phantom measurement, Med Phys 34 (2007), 4236–4246. J.M. Boone, J.A. Seibert, C.M. Tang and S.M. Lane, Grid and Slot Scan Scatter Reduction in Mammography Comparison by Using Monte Carlo Techniques, Radiology 222 (2002), 519–527. X. Chen, R.M. Nishikawa, S.T. Chan, B.A. Lau, L. Zhang and X.Q. Mou, Algorithmic scatter correction in dual-energy digital mammography, Med Phys 40 (2013), 111919. A.H. Baydush and C.E. Floyd, Jr., Improved image quality in digital mammography with image processing, Med Phys 27 (2000), 1503–1508. H. Wieczorek, The image quality of FBP and MLEM reconstruction, Phys Med Biol 55 (2010), 3161–3176. J.Q. Xia, G.D. Tourassi, J.Y. Lo and C.E. Floyd, Jr., On the Development of a Gaussian Noise Model for Scatter, Proc. of SPIE 6510 (2007), 2M1–10. B. Guérin and G.E. Fakhri, Novel Scatter Compensation of List-Mode PET Data Using Spatial and Energy Dependent Corrections, IEEE Trans. on Med Imaging 30 (2011), 759–773. S.J. Feng and I. Sechopoulos, A software-based x-ray scatter correction method for breast tomosynthesis, Med Phys 38 (2011), 6643–6652. A.H. Baydush and C.E. Floyd Jr., Decreasing patient dose in digital mammography with Bayesian image estimation, J X-ray Sci Technol 10 (2001), 87–93. M.R. Lemacks, S.C. Kappadath, C.C. Shaw and X.M. Liu, A dual-energy subtraction technique for microcalcification imaging in digital mammography – A signal-to-noise analysis, Med Phys 29 (2002), 1739–1751. S.C. Kappadath and C.C. Shaw, Dual-energy digital mammography: Calibration and inverse-mapping techniques to estimate calcification thickness and glandular-tissue ratio, Med Phys 30 (2003), 1110–1117. M. Petersilka, K. Stierstorfer and H. Bruder, Strategies for scatter correction in dual source CT, Med Phys 37 (2010), 5971–5992. F.C. Wagner, A. Macovski and D.G. Nishimura, Dual-energy x-ray projection imaging: Two sampling schemes for the correction of scattered radiation, Med Phys 15 (1988), 732–748. R.L. Burden and J.D. Faires, [NUMERICAL ANALYSIS, Ninth Edition], Boston, MA: Brooks Cole (2010), 630–636. F.H. Cheng, X.F. Wang and B.A. Barsky, Quadratic B-Spline Curve Interpolation, Computer Math Appl 41 (2001), 39–50. A. Gálvez, A. Iglesias and J. Puig-Pey, Iterative two-step genetic-algorithm-based method for efficient polynomial Bspline surface reconstruction, Information Sciences 182 (2012), 56–76. X.Q. Mou, X. Chen, L.J. Sun, H.Y. Yu, Z. Ji and L. Zhang, The impact of calibration phantom errors on dual-energy digital mammography, Phys Med Biol 53 (2008), 6321–6336.
K. Ai et al. / A scatter correction method for dual-energy digital mammography [31] [32]
671
I. Kawrakow, egspp: the EGSnrc C++ class library, NRCC Report PIRS-899, Ottawa, Canada, 2005. I. Kawrakow, E. Mainegra-Hing, F. Tessier and B.R.B. Walters, The EGSnrc Code System: Monte Carlo Simulation of Electron and Photon Transport, NRCC Report PIRS-701, Ottawa, Canada, 2011. [33] I. Kawrakow, E. Mainegra-Hing, F. Tessier and B.R.B. Walters, The EGSnrc C++ class library, NRC Report PIRS-898 (rev A), Ottawa, Canada, 2009. [34] G. Poludniowski, P.M. Evans, V. N. Hansen and S. Webb, An efficient Monte Carlo-based algorithm for scatter correction in keV cone-beam CT, Phys Med Biol 54 (2009), 3847–3864. [35] Y. Chen, B. Liu, J.M. O’Connor, C.S. Didier and S.J. Glick, Characterization of scatter in cone-beam CT breast imaging: Comparison of experimental measurements and Monte Carlo simulation, Med Phys 36 (2009), 857–869. [36] S. Kim, H.J. Song, B. Movsas and I.J. Chetty, Characteristics of x-ray beams in two commercial multidetector computed tomography simulators: Monte Carlo simulations, Med Phys 39 (2012), 320–329. [37] J. Lippuner, I.A. Elbakri, C.W. Cui and H.R. Ingleby, Epp: A C++ EGSnrc user code for x-ray imaging, Med Phys 38 (2011), 1705–1708. [38] S. Puong, R. Iordache, X. Bouchevreau and S. Muller, Dual-Energy Contrast Enhanced Digital Mammography: theoretical and experimental study of optimal monoenergetic beam parameters using synchrotron radiation, Proc. of SPIE 7258 (2009), 3U1–10. [39] M. Marziani, A. Taibi, A. Tuffanelli and M. Gambaccini, Dual-energy tissue cancellation in mammography with quasimonochromatic x-rays, Phys Med Biol 47 (2002), 305–313. [40] A.P. Colijn and F.J. Beekman, Accelerated Simulation of Cone Beam X-Ray Scatter Projections, IEEE Trans on Med. Imaging 23 (2004), 584–590. [41] S.C. Kappadath and C.C. Shaw, Dual-energy digital mammography for calcification imaging: noise reduction techniques, Phys Med Biol 53 (2008), 5421–5443.
Copyright of Journal of X-Ray Science & Technology is the property of IOS Press and its content may not be copied or emailed to multiple sites or posted to a listserv without the copyright holder's express written permission. However, users may print, download, or email articles for individual use.
653
A scatter correction method for dual-energy digital mammography: Monte Carlo simulation Kai Aia , Yanhua Gaob and Gang Yua,∗ a School
of Geosciences and Info-Physics, Central South University, Changsha, Hunan, China Imaging Center, Shaanxi Provincial People’s Hospital, Xi’an, Shaanxi, China
b Medical
Received 21 March 2014 Revised 6 June 2014 Accepted 7 July 2014 Abstract. PURPOSE: To develop a novel scatter correction method without additional patient dose for dual-energy digital mammography (DEDM) to reduce scatter’s impacts and enhance microcalcification detectability in dual-energy X-ray subtraction image. METHODS: Combining scatter radiation is lower spatial frequency component and calcifications are sparsely distributed in digital mammogram, we develop a new scatter correction strategy. First, an adaptive sampling scheme is presented to find possible noncalcification (zero calcification) pixels. Then the maximum likelihood expectation maximization (MLEM) algorithm is applied to evaluate initial scatter surface. The accurate scatter radiation of sampling pixels is obtained by solving dual-energy computational formula with zero calcification constraint and scatter surface constraint. RESULTS: After scatter correction, the scatter-to-primary ratio (SPR) of wedge phantom is reduced from ∼36.0% to ∼3.1% for low-energy (LE) image and ∼29.6% to ∼0.6% for high-energy (HE) image. For step phantom, the SPR is reduced from ∼42.1% and ∼30.3% to ∼3.9% and ∼0.9% for LE and HE image, respectively. The calcification contrast-to-noise ratio is improved by two orders of magnitudes in calcification images. CONCLUSIONS: The proposed method shows an excellent performance on scatter reduction and calcification detection. Compared with hardware based scatter correction strategy, our method need no extra exposure and is easy to implementation. Keywords: Scatter correction, dual-energy mammography, calcification, monte carlo
1. Introduction Breast cancer is a major carcinoma that affects women’s health. Intramammary microcalcifications have significant association with early breast cancer [1,2]. Therefore, the efficient detection and visualization of these calcifications in mammograms can provide a good basis for clinical diagnosis. However, conventional screen-film and single-energy digital mammography can only generate the overlapped calcification image, which could be obscured by overlying the fibroglandular structures. ∗ Corresponding author: Gang Yu, School of Geosciences and Info-Physics, Central South University, Changsha 41003, Hunan, China. Tel.: +86 13077327285; E-mail: [email protected].
c 2014 – IOS Press and the authors. All rights reserved 0895-3996/14/$27.50
654
K. Ai et al. / A scatter correction method for dual-energy digital mammography
Dual-energy digital mammography (DEDM) can overcome the inherent drawbacks of conventional single-energy mammography and is viewed as a promising method to enhance the calcification detectability. The DEDM makes use of two different X-ray spectrum to expose the same breast and generate a low-energy (LE) mammographic image and a high-energy (HE) mammographic image. The subtraction mammogram (calcification image especially) is then obtained by calculating dual-energy (DE) computational formula [3–5]. Theoretically, the DE formula can cancel unconcerned tissue, improve image contrast and increase the visualization of very subtle calcifications. The above theory is based on ideal conditions. That is to say, DE subtraction calculation needs accurate data without deteriorating effects. Unfortunately, there are many deleterious factors in practice such as noise interference, anode heel effect and scatter impact [6], which will affect the computational precision. Among them, scatter is one of the most serious effects. According to Cooper’s study, the measured scatter-to-primary ratio (SPR) of 4 cm breast phantom at 28 kVp were ∼39%, ∼35% and ∼42% for 0%, 43% and 100% glandular ratios respectively [7]. And it has little dependence on X-ray spectra range from 22 to 40 kVp [8]. Considering the size of calcification, these quantities of scatter radiation can lead to an enormous calculation error which may eliminate the calcification information in DE subtraction image and destroy the advantages of DEDM [9,10]. So we mainly take scatter correction into account in this paper. Many previous literatures have investigated the performance of different scatter correction methods in DEDM. These methodologies have good performances on scatter suppression and contrast enhancement, but they all have their own limitations. The anti-scatter grid can effectively remove most scatter radiation, but it also wipes off some valuable primary radiation, and this will lead to the reduction of signalto-noise ratio [11]. Pinhole-array interpolation method is used to fit scatter surface and then subtract estimated surface from original image [12]. However, this strategy would increase patient dose by taking an extra exposure. Bornefalk et al. utilize scanned multislit system [13] to eliminate deteriorating effects caused by scatter and other factors. This method is inconvenient to operate, because it needs to precisely align detectors and collimator slits. Saito applies air-gap [14] to achieve the same goal. Nevertheless, air-gap often requires an additional radiation dose, and leads to a reduction of the effective field of view [15]. Chen et al. present an algorithm without extra exposure, but this method easily falls into local optimal solutions and it’s a numerical instability method because of empirical choice of initial scatter values. Besides, the uniform sampling pixels are maybe calcification pixels lead to constraint errors of computational formulae [16]. The main work of this paper is to present a novel DEDM algorithmic scatter correction method. This software based technique exploits adaptive sampling scheme to establish DE equations and then the Bayesian image estimation (BIE) and maximum likelihood expectation maximization (MLEM) algorithm [17–22] is applied to estimate the initial scatter radiation. Finally, the accurate scatter value of sampling pixels is calculated based on noncalcification (zero calcification) constraint and scatter surface (scatter surface is a slowly varying surface) constraint. The whole scatter surface can be obtained by surface fitting. The following sections are organized as below. In Section 2, we introduce the basic DE theory that used in our work. The main idea of our scatter correction method and corresponding dataset are described in detail in Section 3. Section 4 shows experiments and results. At last, we present our discussions and draw a conclusion.
K. Ai et al. / A scatter correction method for dual-energy digital mammography
655
2. Theory According to the principles of X-ray imaging and Lemacks et al.’s work [23], human breast mainly have three different materials: adipose tissue (thickness ta ), glandular tissue (thickness tg ) and calcifications (thickness tc ). The total thickness of breast can be written as T = ta + tg + tc . Due to the negligible thickness (size) of calcifications, the three above variables can be rewritten as two unknowns: glandular ratio (gr = tg /T ≈ tg /(ta + tg )) and tc , where T is a known quantity because the compression device is used in mammography. The incident photon fluence I0 (E) (per unit energy, E ) transmits the compressed breast and then reach detector. In clinical diagnosis, the X-ray spectrum contains different photon energy, so the transmitted signal is obtained by a mathematical integration. With all the energy levels in the high and low X-ray spectrum, the response signals can be described as Elmax Rl = Elmin I0l (E)exp[−μa (E)T − gr (μg (E) − μa (E))T − μc (E)tc ]Q(E)dE , (1) Ehmax Rh = Ehmin I0h (E)exp[−μa (E)T − gr (μg (E) − μa (E))T − μc (E)tc ]Q(E)dE where Rl and Rh are the detector signals produced by the photon spectrum of low and high energy. μa (E), μg (E) and μc (E) are linear attenuation coefficients of adipose tissue, glandular tissue and calcifications at photon energy E , respectively. Elmax and Eimin (i = l, h) are energy range corresponding to LE or HE X-ray spectrum. Q(E) is the detector response. In order to increase the dynamic ranges, two reference signals Rlref and Rhref are applied. Then logarithmic signal can be written as ⎛E ⎞ E lmax lmax Di = ln ⎝ Iref (E)Q(E)dE I0 (E) exp[−μa (E)T−gr (μg (E)−μa (E))T−μc (E)tc ]Q(E)dE ⎠ Eimin
Eimin
= ln(Riref /Ri ) , i = l, h .
(2)
In DEDM, the calcification thickness tc is generated by calculating a cubic inverse-mapping function pixel-by-pixel [24] tc = c0 + c1 Dl + c2 Dh + c3 Dl2 + c4 Dh2 + c5 Dl Dh + c6 Dl3 + c7 Dh3 ,
(3)
where c0 to c7 are constants determined by calibration data. If a certain pixel is zero calcification pixel, tc equals to zero and this pixel is regarded as noncalcification. On the contrary, it is a calcification pixel. Equations (2) and (3) only consider the primary radiation, R. The actual signal received by detector of jth pixel also contains scatter radiation (if other nuisance factors are suppressed), Sj , i.e. Rij = Rij + Sij , i = l, h .
(4)
This may lead to an incorrect thickness of calcifications during DE subtraction calculation. Therefore, to generated accurate tc at jth pixel, we have to estimate and eliminate the contribution of Sij from . Then the estimated D is written as received signal Rij ij
ˆ ij = ln(Riref j (Rij − Sij )) , i = l, h . (5) D For jth pixel, substituting Eq. (5) in Eq. (3) yields an equation with three unknowns: tcj , Slj and Shj . ˆ lj + c2 D ˆ hj + c3 D ˆ 2 + c4 D ˆ 2 + c5 D ˆ lj D ˆ hj + c6 D ˆ 3 + c7 D ˆ3 . tcj = c0 + c1 D lj hj lj hj
(6)
656
K. Ai et al. / A scatter correction method for dual-energy digital mammography
3. Materials and methods To calculate tc , Sl and Sh should be firstly estimated. However, it’s hard to know the accurate scatter surface pixel-to-pixel through a common method. Fortunately, scatter radiation is lower spatial frequency component in X-ray images [25]. Wagner et al. use a polynomial with lower-order (order 2) coefficients to represent it [26]. Kappadath and Shaw assume that scatter can be describe as a quadratic surface function [12] Sl (x, y) = (a0 + a1 x + a2 x2 )(a3 + a4 y + a5 y 2 ) , (7) Sh (x, y) = (b0 + b1 x + b2 x2 )(b3 + b4 y + b5 y 2 ) where (x, y) is spatial location, a0 to a5 and b0 to b5 are constants which can be determined by a leastsquares fit of some known scatter values. In breast cancer, calcifications are sparsely distributed [16]. That is to say, almost all the pixels are zero calcification pixels in LE and HE image. Then, at zero calcification pixels, Eq. (6) can be rewritten as ˆ lj + c2 D ˆ hj + c3 D ˆ 2 + c4 D ˆ 2 + c5 D ˆ lj D ˆ hj + c6 D ˆ 3 + c7 D ˆ 3 = 0. c0 + c1 D lj hj lj hj
(8)
The proposed method utilizes the Eqs (7) and (8) to estimate LE and HE scatter values. It includes three parts. The first is to sample zero calcification pixels. The second is to estimate initial scatter surfaces of LE and HE images. And the final is to calculate scatter values at sampling pixels and then generate the whole scatter surface using quadratic B-spline function fitting method. 3.1. Sampling scheme To calculate scatter value from Eq. (8), the uniform sampling is applied in LE and HE image firstly. The adaptive sampling scheme is then performed to analyze whether the sampling pixels are calcification pixels or not. 3.1.1. Uniform sampling scheme The Eq. (8) of all sampling pixels (circles in Fig. 1(a)), can build a system of nonlinear equations. Because there are two unknowns for Eq. (8) of jth sampling pixel, the scatter values, Slj and Shj , can’t be obtained directly. However, as the scatter radiation is lower spatial frequency component, it can be visualized as a slowly changing variable within a small spatial area. That is to say, one pixel’s scatter value could be expressed as the scatter of its neighboring pixels through Eq. (7) (Fig. 1(b)). Let A1 to A4 are four uniform sampling pixels, and B1 to B4 are their intermediate pixels, i.e. additional sampling pixels. Due to scatter surface constraint, i.e. Eq. (7), the scatter value of B1 (horizontal intermediate point between A1 and A2) can be described as SB1 = (SA1 + SA2 )/2 + Δ,
(9)
where Δ = −a2 ((xA1 − xA2 )/2)2 (a3 + a4 y + a5 y 2 ). It is the error term of scatter value of additional pixel. xA1 , xA2 are horizontal ordinates of SA1 and SA2 , y is vertical ordinates of SA1 , SA2 and SB1 . Others (SB2 , SB3 and SB4 ) are described similarly. The additional sampling pixel Bi is then introduced into the system of nonlinear equations, which can be written as D × C = [ 0 0 · · · 0 · · · 0 0 ]T ,
(10)
K. Ai et al. / A scatter correction method for dual-energy digital mammography
657
(a)
(b)
(c)
Fig. 1. (a) Uniform sampling scheme for LE or HE image (circles). (b) Additional sampling pixels (black squares) in dashed box from Fig. 1(a). A1 to A4 are uniform sampling pixels, and B1 to B4 are additional selected pixels of which the scatter values can be represented by A1 to A4. (c) Adaptive sampling on a n-by-n mask. There are two cases in our adaptive scheme. For a certain sampling pixel j, if the location of pixel j is different from k, then replaces Slj and Shj with Slk and Shk (left), else remain unchanged (right).
where D = [ DA1 DA2 · · · DAn DB1 DB2 · · · DBn ]T , and
T DAi = 1 Dl (SAi ) Dh (SAi ) Dl2 (SAi ) Dh2 (SAi ) Dl Dh (SAi ) Dl3 (SAi ) Dh3 (SAi ) , i = 1, 2, · · · , n ,
T DBi = 1 Dl (SBi ) Dh (SBi ) Dl2 (SBi ) Dh2 (SBi ) Dl Dh (SBi ) Dl3 (SBi ) Dh3 (SBi ) , i = 1, 2, · · · , n , C = [ c0 c1 c2 c3 c4 c5 c6 c7 ]T , n is the number of sampling pixels.
The Eq. (10) is the system of nonlinear equations based on zero calcification constraint and scatter surface constraint. 3.1.2. Adaptive sampling scheme Obviously, the unknowns are equal to the number of equations, and the Slj and Shj can be obtained by solving Eq. (10). However, the above uniform sampling scheme has an inherent drawback. Although the probability of the appearance of calcifications is rather small in early breast cancer mammogram, we can’t make sure that all sampling pixels’ tc in nonlinear equations are equal to zeros (zero calcifications). This deficiency may lead to calculation error in DE subtraction image. Therefore, for each pixel in A sets and B sets, it should be further checked whether it is zero calcification or not. For a certain sampling pixel j , we search for its neighboring pixels with n-by-n mask. Because the calcification pixels are sparsely distributed and exist in isolation, most pixels are zero calcifications under the mask. One pixel k that has most neighbors with similar signal intensity should be considered as zero
658
K. Ai et al. / A scatter correction method for dual-energy digital mammography
calcification. As scatter surface varies slowly, Slj ≈ Slk and Shj ≈ Shk . Therefore, Eq. (10) is also workable when pixel k replaces pixel j . Obviously, if the location of pixel k doesn’t equal to that of pixel j , the pixel k should be used to generate one equation in Eq. (10) instead of pixel j , else pixel j is used. This method is so called adaptive sampling scheme proposed by us (Fig. 1(c)). 3.2. Initial scatter surface estimation As known to all, the convergence of a nonlinear equation highly depends on the choice of initial value [27] and an inappropriate starting point may get an error result (locally optimal solution) or leads to no convergence. The initial value of a0 to a5 and b0 to b5 in Eq. (7) are also needed to be known to fix Δ in Eq. (9). BIE and MLEM algorithm are utilized to estimate initial scatter. The general BIE used in our image estimation can be described as ln f (Ri |Ri ) = ln f (Ri |Ri ) + ln f (Ri ) − ln f (Ri ) , i = l, h ,
(11)
where f (Ri |Ri ) represents a posterior probability of the unknown primary radiation Ri , f (Ri |Ri ) represents the probability related to likelihood function, f (Ri ) is a priori probability of the unknown Ri and f (Ri ) is a constant. According to MLEM algorithm, Eq. (11) can be rewritten as ˆ i = arg max ln f (Ri |R ) = arg max(ln f (R |Ri ) + ln f (Ri )) , i = l, h , R i i Ri
Ri
ˆ i is the estimation of Ri . To calculate the R ˆ i , we use iterative formula presented by [17] where R ˆ (n) · R (T ∗ ∗R ˆ (n) ) , i = l, h , ˆ (n+1) = R R i i i i
(12)
(13)
where, T represents the scatter kernel K (a 2-D exponential function with two parameters: full-width at half-maximum and magnitude) plus a Dirac delta functional, i.e. T = δ + K . The denominator ˆ (0) can be easily set to R . Finally, the estimated scatter represents a convolution computation. And R i i ˆ surface is derived by subtract Ri from Ri . 3.3. Scatter values calculation and surface fitting The accurate scatter values at adaptive sampling points will be calculated by solving system of nonlinear Eq. (10). Finally, the whole accuracy scatter surface is produced by fitting method. In this paper, a quadratic B-spline function [28,29] is used. After subtract the accuracy scatter surfaces from original images, we can obtain a precisely DE calcification image. 3.4. Algorithm process The presented software based scatter correction method for DEDM consists of six successive steps. The brief overview of these steps list as follows: Step 1. Initial scatter surface estimation: The BIE and MLEM algorithm is used to estimate initial Sl and Sh . And then by using least-squares approximation, a0 to a5 and b0 to b5 in Eq. (7) can be obtained through Sl and Sh . Step 2. Adaptive sampling scheme: The uniform sampling scheme is used to sample some pixels. The n-by-n mask is then applied for adaptive sampling to avoid sampling calcification pixels.
K. Ai et al. / A scatter correction method for dual-energy digital mammography
659
Fig. 2. The flow chart of the scatter correction method.
Step 3. Establishment of system of nonlinear equations: We have made sure that the Eq. (8) is workable for each sampling pixel and additional pixel. Therefore, the system of nonlinear equations is built with the unknowns equal to the number of equations just like Eq. (10). Step 4. Estimation for initial solutions: Due to the BIE and MLEM algorithm, the initial scatter values of sampled pixels, Slj and Shj , are gained from Sl and Sh . And the initial Δ value for each additional pixel can also be fixed. old (for the first iteration, S old = Step 5. Equations solving: By using the known values Sljold and Shj lj old = S ), a new group of S new and S new can be obtained by solving Eq. (10) Slj and Shj hj lj hj iteratively.
660
K. Ai et al. / A scatter correction method for dual-energy digital mammography Table 1 Data of adipose and glandular for human breast-tissue-equivalent materials [30] Materials CIRS adipose CIRS glandular
(a)
Density (g/cm3 ) 0.924 1.04
H 11.8 10.9
Elemental compositions (%) C N O Cl Ca 76.0 1.2 9.8 1.2 − 70.2 1.2 12.5 1.1 0.6
Al − 3.5
(b)
Fig. 3. Schematic drawing of designed phantoms. (a) The lateral view of wedge phantom. (b) The lateral view of step phantom. The positive direction of y-axis pointed into paper and the x-axis represented chest wall side. new old new old Step 6. Convergence determination: If the difference between Slj and Slj , Shj and Shj don’t new old ξ , a to meet a certain criterion ξ (a constant), e.g. Sljnew − Sljold ξ and Shj − Shj 0 new and S new , as well as the Δ value. Then a5 and b0 to b5 in Eq. (7) are reinitialized by Slj hj turn back to Step 5. new old ξ or the iteration time is enough, the algorithm is Step 7. If Sljnew − Sljold ξ and Shj − Shj convergent. Then quadratic B-spline function is utilized to generate accurate scatter surface Sl new . And then S and S is subtracted from R and R , we can obtain and Sh by Sljnew and Shj l h l h a precisely primary radiation Rl and Rh . The flow chart of the algorithm is shown in Fig. 2.
3.5. Phantom design A cuboid object with breast-tissue-equivalent materials was used to simulate breast phantom. The length, width and height of the cuboid were 9 cm, 4 cm and 5 cm, respectively. The elemental compositions data (Table 1) of the adipose and glandular tissue-equivalent materials were provided by CIRS (Computerized Imaging Reference Systems, Inc., Norfolk, VA). In our experiment, two phantoms with different inner structures were designed. The first wedge phantom was used to simulate the slowly varying adipose-glandular ratio in breast. The second step phantom was used to simulate the rapid varying adipose-glandular ratio. Usually, 5 cm was the common thickness of compressed breast and the microcalcifications in breast were mainly composed of CaCO3 and reached to several hundred microns. So we set the thickness of the phantoms to 5 cm and used CaCO3 cubes to simulate calcifications. The wedge phantom was combined by two reversed wedges (Fig. 3(a)) with glandular ratio ranging from 0% to 100%, and the step phantom was made up of nine steps (Fig. 3(b)). The length of each
K. Ai et al. / A scatter correction method for dual-energy digital mammography
661
(a)
(b) Fig. 4. (a) The HE step phantom primary image. The calcifications size varied along the vertical direction and stayed fixed along the horizontal direction. And glandular ratio ranged from 10% to 90% from left to right. (b) The HE step phantom scatter image.
step was 1 cm and the height was 0.5 cm, i.e. the glandular ratio ranged from 10% to 90% with each step interval is 10%. And then three different sizes of CaCO3 cubes (200 μm, 300 μm, 400 μm and the mass density is 2.93g/cm3) were placed on the middle top of two phantoms to simulate calcifications. Calcifications were arranged to a 3×3 pattern with the size varying along the y-axis and stayed fixed along the x-axis. All the geometries above were generated by egspp [31] (a C++ class library provided by EGSnrc package [32,33]). 3.6. Image acquisition Monte Carlo (MC) method has been widely used for simulating x-ray scatter radiation and X-ray imaging [34–36]. And the accuracy of many MC simulation packages (e.g. EGSnrc, GEANT and MCNP) has been proved. In our experiment, we developed a new EGSnrc C++ user code to generate LE and HE primary and scatter data based on Lippuner et al.’s work [37]. Besides, we rewrote the phantom design code to generate wedge phantom and step phantom. During image acquisition, incident X-ray photons irradiated the designed phantom from top, and a virtual detector with 400 × 900 pixels (i.e. 100 × 100 μm2 for each pixel) was placed right after the bottom of phantom. For easy implementation, we made use of the monochromatic x-ray energies [38, 39] with 18 keV for LE imaging and 36 keV for HE imaging. To save simulating time, we applied Richardson-Lucy fitting algorithm [40] to process the raw scatter data. Figure 4(a) and (b) shown the HE step phantom primary and scatter image that acquired from our EGSnrc based C++ program (LE step phantom, LE and HE wedge phantoms were just similar to HE step phantom). The nine black dots were calcifications with different sizes.
662
K. Ai et al. / A scatter correction method for dual-energy digital mammography
Fig. 5. Schematic drawing of calibration phantom (top view). The dash lines were used to distinguish different glandular ratios, and the dotted squares represented 0 µm calcifications. The positive direction of z-axis pointed out of paper and the x-axis represented chest wall side.
3.7. Calibration procedure The aim of calibration procedure was to confirm coefficients in Eq. (3). To this end, we designed a calibration phantom of 5 cm thickness with five steps (the glandular ratios were 0%, 30%, 50%, 80% and 100%), and on each step we placed six different thickness of CaCO3 cubes (0, 100, 350, 500, 750, 1200 μm) orthogonal to x-axis (Fig. 5). The value of c0 to c7 were determined by least squares fit of thirty acquired data. And the reference signals, Rlref and Rhref were measured through calibration phantom with 100% adipose ratio. According to Kappadath and Shaw’s experiment [23], the mean fitting error of Eq. (6) was about 50 μm and the max fitting error was 150 μm for microcalcification thickness. As our experiment is a simulation program, the mean and max errors were suffered a further reduction. This meant the Eq. (6) is valid for modeling calcification thickness in this paper. 3.8. Evaluation metrics The presented method was used to cleanup scatter radiation. Scatter-to-primary ratio (SPR) was used to estimate the number of scatter components in original LE and HE image. In contrast, residual scatterto-primary ratio (rSPR) represented how much scatter components remained after scatter correction. It could be written as Sij SPR = , i = l, h , (14) Rij Sij − Sij rSPR = , i = l, h , (15) Rij represents the estimated scatter value, j is the pixel index. where Sij To prove our method could enhance calcifications contrast effectively, we also employed calcification contrast-to-noise ratio (cCNR) which was described as |Vc − Vb | cCNR = , (16) σb where Vb and σb are average signal intensity and standard deviation of background. Vc is average signal intensity of calcification.
K. Ai et al. / A scatter correction method for dual-energy digital mammography
(a)
(b)
(c)
(d)
663
Fig. 6. (a) to (b) The SPR image of LE and HE wedge phantom before scatter correction. (c) to (d) The rSPR image of LE and HE wedge phantom after scatter correction.
3.9. Scatter correction settings In this paper, we sampled 256 pixels, including 128 original sampling pixels and 128 additional sampling pixels. A 5-by-5 mask for adaptive sampling was proved to have a good performance in the following experiment. In Eq. (13), that two parameters of scatter kernel K were fixed to 50 pixels for full-width at halfmaximum and 0.029, 0.032 for magnitude of LE and HE scatter estimation respectively can work well. After 10 times iterations, Eq. (10) became convergent. 4. Results The LE and HE scatter surfaces are estimated by the presented method and subtracted from original LE and HE images, Rl and Rh . Then Eq. (6) is used to generate calcification image based on the corrected primary radiation. The SPR, rSPR and cCNR are calculated to evaluate the performance of our method.
664
K. Ai et al. / A scatter correction method for dual-energy digital mammography
(a)
(b)
(c)
(d)
Fig. 7. (a) to (b) The SPR image of LE and HE step phantom before scatter correction. (c) to (d) The rSPR image of LE and HE step phantom after scatter correction.
4.1. SPR and rSPR of two phantoms Figure 6(a) and (b) show the SPR of LE and HE wedge phantom of sampling pixels. Figure 6(c) and (d) show the rSPR of sampling pixels after scatter correction. The negative value of rSPR represents our method has overestimated scatter radiation. In Fig. 6 (a) and (b), the mean SPR of LE and HE wedge phantom of sampling pixels are ∼36.0% and ∼29.6% before scatter correction, respectively. By contrast, the mean absolute value of rSPR is suppressed to ∼3.1% and 0.6% after apply the proposed method. Our correction strategy has a very good performance on scatter reduction of wedge phantom. Figure 7 shows the corresponding result of step phantom. As can be seen from Fig. 7(a), the SPR of LE step phantom varies quickly due to the shape of phantom. The mean SPR of LE and HE step phantom of sampling pixels are ∼42.1% and ∼30.3% before scatter correction, respectively. And after scatter correction, the mean absolute value of rSPR becomes to ∼3.9% and ∼0.9%. That is to say, our method is also suitable for steep tissue in breast.
K. Ai et al. / A scatter correction method for dual-energy digital mammography
665
Table 2 The statistical analysis of SPR and absolute value of rSPR of sampling pixels Phantom shape Wedge (LE) Wedge (HE) Step (LE) Step (HE)
Mean 35.96 29.60 42.09 30.33
SPR (%) Max Min 40.00 27.25 33.26 21.95 53.75 26.63 34.93 21.05
Std. 2.92 2.63 6.24 3.21
Mean 3.09 0.61 3.86 0.85
|rSPR| (%) Max Min 8.21 0.02 1.63 0.01 11.05 0.15 3.02 0.00
Std. 2.23 0.44 3.04 0.83
Table 3 The statistical analysis of the absolute value of rSPR of sampling pixels using Chen et al.’s method Phantom shape Wedge (LE) Wedge (HE) Step (LE) Step (HE)
Mean 7.31 1.32 11.95 2.11
|rSPR |(%) Max Min 11.63 0.33 2.07 0.06 22.13 0.10 3.80 0.03
Std. 2.91 0.52 5.55 0.95
Mean 12.60 2.50 4.80 0.90
|rSPR| (%) Max Min 22.25 7.33 4.59 1.41 13.51 0.04 2.69 0.00
Std. 3.55 0.75 3.20 0.60
Table 4 The statistical analysis of absolute value of rSPR of fitting scatter surface Phantom shape Wedge (LE) Wedge (HE) Step (LE) Step (HE)
Mean 2.51 0.49 3.38 0.72
Fitting data |rSPR| (%) Max Min 7.91 0.00 1.58 0.00 11.14 0.00 3.17 0.00
Std. 1.63 0.32 2.96 0.78
Table 2 shows the maximum, minimum, mean and standard deviation of SPR and absolute value of rSPR of sampling pixels. In contrast, Chen et al. just use empirical method to select initial scatter estimation. It is a numerical instability method. According to the SPR data of Table 2, we empirically choose two different sets of initial scatter values (35% and 45% of the total received signal, Ri ) to test Chen et al.’s method and some statistical values are list below. The left and right half of Table 3 shows the statistical values with initial scatter estimation equals to 0.35×Ri and 0.45×Ri . Obviously, with initial scatter estimation changes, the absolute value of statistical values of rSPR varies seriously. This indicates that Chen et al.’s scatter correction method completely depends on the choice of initial scatter estimation. In practice, we have little prior information about SPR. So it is a hard work to set suitable initial scatter surfaces by empirical method. 4.2. Surface fitting To generate DE calcification image, we have to cleanup the whole scatter surface. Here the quadratic B-spline function and least-squares approximation are utilized to fit the smooth scatter surface through estimated scatter values at sampling pixels. Figure 8 shows the central row profiles of real scatter surface and fitting surface (normalized). The mean fitting relative error of LE and HE wedge phantom of central rows are ∼4.2% and ∼0.9%, respectively. For LE and HE step phantom, the values are ∼9.2% and ∼2.6%, respectively. The simulation experiment has demonstrated that the chose of function and order are entirely satisfied our expecta-
666
K. Ai et al. / A scatter correction method for dual-energy digital mammography
(a)
(b)
(c)
(d)
Fig. 8. The central row profiles of real scatter data (solid line), fitting data (dash-dot line) and error data (dotted line). (a) to (b) are wedge phantom of LE and HE image, respectively. (c) to (d) are step phantom of LE and HE image, respectively.
tion. Table 4 shows the maximum, minimum, mean and standard deviation of absolute value of rSPR of the whole fitting scatter surface. 4.3. Calcification image Figure 9 shows local calcification image of wedge and step phantom (2 × magnification for clear visual inspection) which generated by calculating Eq. (6). In Fig. 9, we can clearly see calcifications with the size of 300 μm and 400 μm. Due to the small size of 200 μm calcifications, the visibility is weaker than 300 μm and 400 μm. In these two local images, the glandular tissue is removed by DE algorithm and we can faintly see some residual signal (‘white’ background), but this interference signal hardly affects the visibility of calcifications.
K. Ai et al. / A scatter correction method for dual-energy digital mammography
(a)
667
(b)
Fig. 9. (a) The DE calcification image of wedge phantom. (b) The DE calcification image of step phantom.
(a)
(b)
Fig. 10. The cCNR of DE calcification images. (a) The cCNR of nine calcification pixels of wedge phantom. (b) The cCNR of nine calcification pixels of step phantom.
Figure 10 shows the cCNR of nine calcification pixels. We choose nine 20 × 20 region of interests (ROIs) as background to calculate the cCNR of the nine pixels with 200 μm, 300 μm, 400 μm calcifications on different glandular ratios. The ROIs are just nearby the corresponding calcifications. For wedge phantom, the mean cCNR of 200 μm calcifications raises from 6.45 to 124.31. For 300 μm and 400 μm calcifications, the cCNR raises from 17.97, 2.40 to 188.06 and 296.00, respectively. For step phantom, the mean cCNR for 200 μm, 300 μm, 400 μm calcifications raise from 5.43, 5.28, 4.85 to 131.44, 158.21 and 282.95, respectively. Obviously, the cCNR of DE calcification image is significantly improved after scatter correction. In Table 5, the thickness of nine calcification pixels after scatter correction is shown outside the parentheses and the uncorrected values are shown in parentheses (use Eq. (3) to obtain them). The calculated thickness after scatter correction is close to the real thickness. And the max errors of 200 μm, 300 μm and 400 μm calcification are only 28 μm, 65 μm and 74 μm for wedge phantom and 35 μm, 68 μm and 106 μm for step phantom.
668
K. Ai et al. / A scatter correction method for dual-energy digital mammography Table 5 The calcification thickness (µm) in DE image
Real thickness ( µm) 200 300 400
Glandular ratio of wedge phantom increases from left to right 215 (3347) 191 (3420) 172 (3450) 313 (3394) 276 (3478) 235 (3504) 474 (3202) 427 (3263) 432 (3302)
Glandular ratio of step phantom 40% 50% 60% 235 (3276) 221 (3337) 191 (3359) 294 (3315) 265 (3347) 232 (3388) 397 (3108) 350 (3150) 294 (3155)
5. Discussions As compared with single-energy mammography, DEDM can remove unwanted tissue to improve the detectibility of microcalcifications. Therefore, it has a good prospect for breast imaging. Due to the deteriorating effects of scatter radiation, however, the thickness of calcifications may become incorrect. Based on the two characteristics of DEDM (i.e. X-ray scatter is a lower spatial frequency component and calcifications are sparsely distributed), this paper present a novel software based scatter correction method. From a theoretical perspective, the presented method includes some key points. Because the noncalcification tissue plays a dominant role in breast, so the thickness of calcifications can be regarded as zero in most mammogram regions. That is to say, most pixels are satisfied with Eq. (8). By uniformly sampling some pixels in LE and HE image and substituting these pixels in Eq. (8), a system of nonlinear equations can be built. Nevertheless, this system of nonlinear equations can’t be solved, because the number of the equations is smaller than the unknowns. Fortunately, Eq. (7) is used to represent the low-frequency variation of the scatter surface and the scatter value of intermediate pixel between two sampling pixels can be expressed by Eq. (9). The introduction of these additional pixels makes the number of the equations equivalent to the unknowns and then a solvable Eq. (10) is obtained. Meanwhile, the error term of scatter value of additional pixel is also included in Eq. (9) to increase the computation accuracy. It’s worth noting that Eq. (10) based on uniform sampling scheme has an inherent shortcoming where we can’t ensure Eq. (8) is workable for all the sampling pixels and additional pixels. Thus the adaptive sampling scheme is necessarily utilized. This innovative point can ensure that all the selected pixels are noncalcification pixels. The initial solution estimation is also important for Eq. (10). If the inappropriate initial scatter values are used such as reference [16], Eq. (10) may be fall into local optimal solution or directly not convergence. BIE and MLEM algorithm is employed to generate a good initial scatter surface estimation and fix initial values of a0 to a5 and b0 to b5 at the same time. It plays an important role in the presented method. The results indicated that our method has a good performance on scatter correction. It can improve the visibility of microcalcifications with the size of ∼200 μm which is unachievable in reference [12] and other hardware based strategies. Compared with other algorithmic correction, our method is adaptive process and need no empirical set of initial scatter values. Despite all this, there are still some useful information can be get from Figs 6, 7 and Table 2. Most important among them is that the correction effectiveness for different shape of phantoms is not the same. The absolute value of rSPR of LE wedge phantom image is 3.09%, 0.76% less than step phantom. And for HE wedge phantom image, the absolute value of rSPR is 0.61%, 0.24% less than step phantom. This phenomenon may be caused by two reasons. Firstly, compared with the slowly varying glandular ratio in wedge phantom, the step phantom contains rapid change of glandular ratio. The boundary effect can constitute a handicap for MLEM algorithm to estimate initial scatter surface, which will further affect the solution of Eq. (10). In addition, as can be
K. Ai et al. / A scatter correction method for dual-energy digital mammography
669
seen from Table 2, the mean SPR of step LE and HE phantom image is greater than the corresponding wedge phantom image. Under the same simulation conditions and system parameters, the original scatter radiation of step phantom image may be more than that of wedge phantom, which reduces the performance of scatter correction. For the same phantom, the scatter correction effectiveness is also different between LE and HE image. The performance of scatter correction in HE image is much better than LE image no matter which phantom is used. For step phantom, scatter fraction reduces from 35.96% to 3.09% in LE image and 29.60% to 0.61% in HE image. For wedge phantom, the fraction reduces from 42.09% to 3.86% in LE image and 30.33% to 0.85% in HE image. The observed difference may be caused by the amount of scatter radiation in LE image. Compared with the scatter fraction in HE image, the larger scatter fraction in LE image leads to the reduction of primary fraction. This will affect complex calculation in Eq. (10) and lead to the error. As can be seen from Table 4, the error between the fitting scatter surface and real scatter surface suffers from a further reduction in both LE and HE image after a quadratic B-spline function is applied. And still this observation is in agreement with the findings in which the correction performance is better in wedge phantom image than step phantom image, and HE image than LE image. However, the calcification images shown in Fig. 10 and Table 5 demonstrate that the above difference has not serious impact on the visibility of calcifications. Besides the scatter radiation, X-ray radiation contains many other interference factors. Noise is another important effect among them. For simplicity, the presented method assumes that the noise in the received radiation has been removed or depressed by many known noise reduction methods [41]. Future work is underway to integrate both scatter and noise in received radiation and a new methodology which can fully correct the two effects will be carried out. 6. Conclusions In this paper, we present a novel software based scatter correction used for DEDM. This method can reduce scatter fraction within 4% and 1% for LE and HE image, and enhance the microcalcifications detectability at the same time. The cCNR is improved by two orders of magnitudes. Compared with hardware based scatter correction, the proposed method is easy to implementation and need no extra exposure. Our next work will take other deleterious effects into account (e.g. noise) and make a further investigation in calculation error to develop a complete, accurate and clinical DEDM processing application. Acknowledgments The authors would like to thank the following funding organizations for support: Hunan Provincial Natural Science Foundation of China (No. 11JJ6507) and China Postdoctoral Science Foundation (No. 20090461303), and the Fundamental Research Funds for the Central Universities of Central South University (No. 2013zzts251). References [1]
R.W. Powell, M.B. McSweeney and C.E. Wilson, X-ray calcifications as the only basis for breast biopsy, Ann Surg 197 (1983), 555–559.
670 [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30]
K. Ai et al. / A scatter correction method for dual-energy digital mammography A. Fandos-Morera, M. Prats-Esteve, J.M. Tura-Soteras and A. Traveria-Cros, Breast tumors: composition of microcalcifications, Radiology 169 (1988), 325–327. L.A. Lehmann, R.E. Alvarez, A. Macovski and W.R. Brody, Generalized image combination in dual KVP digital radiography, Med Phys 8 (1981), 659–667. A. Taibi, S. Fabbri, P. Baldelli, C. di Maggio, G. Gennaro, M. Marziani, A. Tuffanelli and M. Gambaccini, Dual-energy imaging in full-field digital mammography: A phantom study, Phys Med Biol 48 (2003), 1945–1956. S.C. Kappadath and C.C. Shaw, Quantitative evaluation of dual-energy digital mammography for calcification imaging, Phys Med Biol 49 (2004), 2563–2576. J.S. Maltz, B. Gangadharan, S. Bose, D.H. Hristov, B. A. Faddegon, A. Paodi and A.R. Bani-Hashemi, Algorithm for X-ray Scatter, Beam-Hardening, and Beam Profile Correction in Diagnostic (Kilovoltage) and Treatment (Megavoltage) Cone Beam CT, IEEE Trans. on Med Imaging 27 (2008), 1791–1810. V.N. Cooper III, J.M. Boone, J.A. Seibert and C.J. Pellot-Barakat, An edge spread technique for measurement of the scatter-to-primary ratio in mammography, Med Phys 27 (2000), 845–853. J.M. Boone, K.K. Lindfors, V.N. Cooper III and J.A. Seibert, Scatter/primary in mammography Comprehensive results, Med Phys 27 (2000), 2408–2416. T.B. Hamed and R.D. Speller, The effect of scattered radiation in dual-energy analysis, Phys Med Biol 40 (1995), 1619– 1632. F.C. Wagner, A. Macovski and D.G. Nishimura, Effects of Scatter in Dual-Energy Imaging: An Alternative Analysis, IEEE Trans. on Med Imaging 8 (1989), 236–244. S.Z. Shen, A.K. Bloomquist, G.E. Mawdsley and M.J. Yaffe, Effect of scatter and an antiscatter grid on the performance of a slot-scanning digital mammography system, Med Phys 33 (2006), 1108–1115. S.C. Kappadath and C.C. Shaw, Dual-energy digital mammography for calcification imaging: Scatter and nonuniformity corrections, Med Phys 32 (2005), 3395–3408. H. Bornefalk, M. Hemmendorff and T. Hjärn, Contrast-enhanced dual-energy mammography using a scanned multislit system: evaluation of a differential beam filtering technique, Journal of Electronic Imaging 16 (2007), 0230061– 0230067. M. Saito, Dual-energy approach to contrast-enhanced mammography using the balanced filter method: Spectral optimization and preliminary phantom measurement, Med Phys 34 (2007), 4236–4246. J.M. Boone, J.A. Seibert, C.M. Tang and S.M. Lane, Grid and Slot Scan Scatter Reduction in Mammography Comparison by Using Monte Carlo Techniques, Radiology 222 (2002), 519–527. X. Chen, R.M. Nishikawa, S.T. Chan, B.A. Lau, L. Zhang and X.Q. Mou, Algorithmic scatter correction in dual-energy digital mammography, Med Phys 40 (2013), 111919. A.H. Baydush and C.E. Floyd, Jr., Improved image quality in digital mammography with image processing, Med Phys 27 (2000), 1503–1508. H. Wieczorek, The image quality of FBP and MLEM reconstruction, Phys Med Biol 55 (2010), 3161–3176. J.Q. Xia, G.D. Tourassi, J.Y. Lo and C.E. Floyd, Jr., On the Development of a Gaussian Noise Model for Scatter, Proc. of SPIE 6510 (2007), 2M1–10. B. Guérin and G.E. Fakhri, Novel Scatter Compensation of List-Mode PET Data Using Spatial and Energy Dependent Corrections, IEEE Trans. on Med Imaging 30 (2011), 759–773. S.J. Feng and I. Sechopoulos, A software-based x-ray scatter correction method for breast tomosynthesis, Med Phys 38 (2011), 6643–6652. A.H. Baydush and C.E. Floyd Jr., Decreasing patient dose in digital mammography with Bayesian image estimation, J X-ray Sci Technol 10 (2001), 87–93. M.R. Lemacks, S.C. Kappadath, C.C. Shaw and X.M. Liu, A dual-energy subtraction technique for microcalcification imaging in digital mammography – A signal-to-noise analysis, Med Phys 29 (2002), 1739–1751. S.C. Kappadath and C.C. Shaw, Dual-energy digital mammography: Calibration and inverse-mapping techniques to estimate calcification thickness and glandular-tissue ratio, Med Phys 30 (2003), 1110–1117. M. Petersilka, K. Stierstorfer and H. Bruder, Strategies for scatter correction in dual source CT, Med Phys 37 (2010), 5971–5992. F.C. Wagner, A. Macovski and D.G. Nishimura, Dual-energy x-ray projection imaging: Two sampling schemes for the correction of scattered radiation, Med Phys 15 (1988), 732–748. R.L. Burden and J.D. Faires, [NUMERICAL ANALYSIS, Ninth Edition], Boston, MA: Brooks Cole (2010), 630–636. F.H. Cheng, X.F. Wang and B.A. Barsky, Quadratic B-Spline Curve Interpolation, Computer Math Appl 41 (2001), 39–50. A. Gálvez, A. Iglesias and J. Puig-Pey, Iterative two-step genetic-algorithm-based method for efficient polynomial Bspline surface reconstruction, Information Sciences 182 (2012), 56–76. X.Q. Mou, X. Chen, L.J. Sun, H.Y. Yu, Z. Ji and L. Zhang, The impact of calibration phantom errors on dual-energy digital mammography, Phys Med Biol 53 (2008), 6321–6336.
K. Ai et al. / A scatter correction method for dual-energy digital mammography [31] [32]
671
I. Kawrakow, egspp: the EGSnrc C++ class library, NRCC Report PIRS-899, Ottawa, Canada, 2005. I. Kawrakow, E. Mainegra-Hing, F. Tessier and B.R.B. Walters, The EGSnrc Code System: Monte Carlo Simulation of Electron and Photon Transport, NRCC Report PIRS-701, Ottawa, Canada, 2011. [33] I. Kawrakow, E. Mainegra-Hing, F. Tessier and B.R.B. Walters, The EGSnrc C++ class library, NRC Report PIRS-898 (rev A), Ottawa, Canada, 2009. [34] G. Poludniowski, P.M. Evans, V. N. Hansen and S. Webb, An efficient Monte Carlo-based algorithm for scatter correction in keV cone-beam CT, Phys Med Biol 54 (2009), 3847–3864. [35] Y. Chen, B. Liu, J.M. O’Connor, C.S. Didier and S.J. Glick, Characterization of scatter in cone-beam CT breast imaging: Comparison of experimental measurements and Monte Carlo simulation, Med Phys 36 (2009), 857–869. [36] S. Kim, H.J. Song, B. Movsas and I.J. Chetty, Characteristics of x-ray beams in two commercial multidetector computed tomography simulators: Monte Carlo simulations, Med Phys 39 (2012), 320–329. [37] J. Lippuner, I.A. Elbakri, C.W. Cui and H.R. Ingleby, Epp: A C++ EGSnrc user code for x-ray imaging, Med Phys 38 (2011), 1705–1708. [38] S. Puong, R. Iordache, X. Bouchevreau and S. Muller, Dual-Energy Contrast Enhanced Digital Mammography: theoretical and experimental study of optimal monoenergetic beam parameters using synchrotron radiation, Proc. of SPIE 7258 (2009), 3U1–10. [39] M. Marziani, A. Taibi, A. Tuffanelli and M. Gambaccini, Dual-energy tissue cancellation in mammography with quasimonochromatic x-rays, Phys Med Biol 47 (2002), 305–313. [40] A.P. Colijn and F.J. Beekman, Accelerated Simulation of Cone Beam X-Ray Scatter Projections, IEEE Trans on Med. Imaging 23 (2004), 584–590. [41] S.C. Kappadath and C.C. Shaw, Dual-energy digital mammography for calcification imaging: noise reduction techniques, Phys Med Biol 53 (2008), 5421–5443.
Copyright of Journal of X-Ray Science & Technology is the property of IOS Press and its content may not be copied or emailed to multiple sites or posted to a listserv without the copyright holder's express written permission. However, users may print, download, or email articles for individual use.
