A new stationary gridline artifact suppression method based on the 2D discrete wavelet transform Hui Tang, Dan Tong, Xu Dong Bao, and Jean-Louis Dillenseger Citation: Medical Physics 42, 1721 (2015); doi: 10.1118/1.4914861 View online: http://dx.doi.org/10.1118/1.4914861 View Table of Contents: http://scitation.aip.org/content/aapm/journal/medphys/42/4?ver=pdfcov Published by the American Association of Physicists in Medicine Articles you may be interested in Microcalcification detection based on wavelet domain hidden Markov tree model: Study for inclusion to computer aided diagnostic prompting system Med. Phys. 34, 2206 (2007); 10.1118/1.2733800 Performance analysis of a new semiorthogonal spline wavelet compression algorithm for tonal medical images Med. Phys. 27, 276 (2000); 10.1118/1.598830 A wavelet-based algorithm for detecting clustered microcalcifications in digital mammograms Med. Phys. 26, 1294 (1999); 10.1118/1.598624 An introduction to wavelet theory and application for the radiological physicist Med. Phys. 25, 1985 (1998); 10.1118/1.598387 Optimally weighted wavelet transform based on supervised training for detection of microcalcifications in digital mammograms Med. Phys. 25, 949 (1998); 10.1118/1.598273
A new stationary gridline artifact suppression method based on the 2D discrete wavelet transform Hui Tanga) Laboratory of Image Science and Technology, School of Computer Science and Engineering, Southeast University, Nanjing 210096, China; Key Laboratory of Computer Network and Information Integration (Southeast University), Ministry of Education, Nanjing 210000, China; Centre de Recherche en Information Biomédicale sino-français, Laboratoire International Associé, Inserm, Université de Rennes 1, Rennes 35000, France; and Southeast University, Nanjing 210000, China
Dan Tong and Xu Dong Bao Laboratory of Image Science and Technology, School of Computer Science and Engineering, Southeast University, Nanjing 210096, China
Jean-Louis Dillenseger INSERM, U1099, Rennes F-35000, France; Université de Rennes 1, LTSI, Rennes F-35000, France; Centre de Recherche en Information Biomédicale sino-français, Laboratoire International Associé, Inserm, Université de Rennes 1, Rennes 35000, France; and Southeast University, Nanjing 210000, China
(Received 28 August 2014; revised 12 February 2015; accepted for publication 17 February 2015; published 18 March 2015) Purpose: In digital x-ray radiography, an antiscatter grid is inserted between the patient and the image receptor to reduce scattered radiation. If the antiscatter grid is used in a stationary way, gridline artifacts will appear in the final image. In most of the gridline removal image processing methods, the useful information with spatial frequencies close to that of the gridline is usually lost or degraded. In this study, a new stationary gridline suppression method is designed to preserve more of the useful information. Methods: The method is as follows. The input image is first recursively decomposed into several smaller subimages using a multiscale 2D discrete wavelet transform. The decomposition process stops when the gridline signal is found to be greater than a threshold in one or several of these subimages using a gridline detection module. An automatic Gaussian band-stop filter is then applied to the detected subimages to remove the gridline signal. Finally, the restored image is achieved using the corresponding 2D inverse discrete wavelet transform. Results: The processed images show that the proposed method can remove the gridline signal efficiently while maintaining the image details. The spectra of a 1D Fourier transform of the processed images demonstrate that, compared with some existing gridline removal methods, the proposed method has better information preservation after the removal of the gridline artifacts. Additionally, the performance speed is relatively high. Conclusions: The experimental results demonstrate the efficiency of the proposed method. Compared with some existing gridline removal methods, the proposed method can preserve more information within an acceptable execution time. C 2015 American Association of Physicists in Medicine. [http://dx.doi.org/10.1118/1.4914861] Key words: x-ray image, stationary grid, wavelet transform, gridline suppression, Gaussian filter 1. INTRODUCTION X-ray radiography is one of the most common imaging modalities in the medical field for diagnosis. Compared with analog radiographic systems that use screen/film combinations, digital radiographic (DR) systems directly produce digital images with optimal density and brightness for interpretation or diagnosis, because the detected signal can be modified by electronic adjustment of the system gain.1 When the x-ray photons strike the object being imaged, some photons are scattered, which leads to contrast degradation and reduced image quality. For both analog and digital systems, the most commonly used technique for reducing the scattered radiation is to insert an antiscatter grid between the patient and the image receptor.2 The use of antiscatter 1721
Med. Phys. 42 (4), April 2015
grids has been proven to increase the image contrast and the signal-to-noise ratio in radiographic images.3–6 Two types of grid, stationary or moving, are used according to the situation. Compared with stationary grids, systems with moving grids require a more complex system architecture. If a grid is stationary during the exposure, the grid shadows may be visible in the image and are known as gridline artifacts. Numerous approaches have been undertaken to suppress the gridline artifacts in x-ray images.7–18 Some researchers have reported methods that avoid gridline artifacts by adjusting some of the physical modules during the x-ray projection and detection process in the imaging system.7–11 These methods are limited in some specific situations because of the induced system architecture complexity.
0094-2405/2015/42(4)/1721/9/$30.00
© 2015 Am. Assoc. Phys. Med.
1721
1722
Tang et al.: Gridline suppression
1722
Comparatively, image processing methods are more adaptive.12–22 Stationary gridline artifacts are expressed as regularly distributed periodic signals in the x-ray images. The different image processing based artifact suppression methods can be classified into three categories: 1. Space domain methods.12–15 In this class of methods, the gridline artifacts are analyzed according to the image gray level information and suppressed in the spatial domain. The disadvantage of these methods is that the gridline detection process is usually affected by the object image information and the convolution gridline suppression process will blur the object image. 2. Frequency domain filtering methods.16–20 In the frequency domain, the gridline signals only appear at certain frequencies related to the image acquisition situation. With prior knowledge, the gridline artifacts and the imaging object are relatively easy to separate in the frequency domain. But when the grid frequency is suppressed in the frequency domain, the object image information with a similar frequency will also be reduced. 3. Wavelet domain suppression methods.21,22 In this method, a recursive wavelet transform is applied to the image until subimages of the desired resolution are generated. The subimages that are expected to contain the gridline artifacts are set to zero, and finally the restored image is obtained using the inverse wavelet transform. These methods are proposed to preserve more of the object information, but the disadvantage is that these methods did not provide a stop condition for the recursive wavelet transform and the zero assignment operation will lead to a ringing effect in the final image. To efficiently remove gridline artifacts from x-ray images while preserving object information, we propose a new gridline artifact suppression method. First, the image with gridline artifacts is decomposed into several subimages using a recursive wavelet transform. The decomposition is stopped when the gridline signals are detected to be dominant in one or several of the subimages using an adaptive gridline detection module. An adaptive Gaussian band-stop filter is then applied to this or these subimages, which is or are replaced by the filtered version(s). Finally, an inverse wavelet transform is applied to the wavelet pyramid to obtain a restored image without gridline artifacts. 2. GRIDLINE ARTIFACT EXPRESSION IN AN X-RAY IMAGE In the stationary grid case, the x-rays emitted by the source pass through and are thus absorbed by the target object and the stationary antiscatter grid. Both will leave shadows on the final image. The grid’s ideal shadow is show in Fig. 1, denoted as Sgrid (x), multiplies the desired image on the final image. It can be expressed using the Fourier expansion Sgrid (x) =
d Tgrid
+
∞
(ai cos(2iπx f g )).
i=1
Medical Physics, Vol. 42, No. 4, April 2015
(1)
F. 1. Ideal gridline signal in the space domain.
According to Lin’s analysis,18 from this formula, we can get d/Tgrid, the DC term; hg 1 (x) = a1 cos(2πx f g ), the fundamental wave with the frequency f g ; and hg m (x), m = 2, ..., ∞, the harmonics with the frequency m f g for each term. Since the harmonics have a higher frequency but a lower magnitude than the fundamental wave, the effect of the harmonics can be considered small and can be ignored. For the following discussion, we will only take the fundamental wave into account. Equation (1) describes the gridline artifacts in the presampled image. However, the digital image is also affected by the sampling function. We will denote the sampling frequency as f s . According to the Nyquist sampling theorem, if f s ≥ 2 f g , the data at f g can be reconstructed without aliasing. Otherwise, aliasing occurs.23 So according to the relationship between the sampling frequency and the grid frequency, we can estimate the grid artifact frequency, denoted as f estimate, by if f s ≥ 2 f g , f estimate = f g ,
(2)
else, f estimate = f g − k f s , (3) ( ) ( ) fg fg if mod < 0.5 floor fs fs ( ) , (4) k≡ fg otherwise ceil fs where floor(x) is x rounded down to the nearest integer and ceil(x) is x rounded up to the nearest integer. 3. PROPOSED METHOD 3.A. Algorithm outline
The goal of our work is to remove grid patterns from the x-ray image while preserving the object information. Wavelet transformation is a useful tool to separate signals, so we endeavor to implement this tool to focus the image processing on only the gridline artifacts’ part of the signal instead of the whole image.
1723
Tang et al.: Gridline suppression
1723
F. 2. Gridline artifact suppression method based on the 2D recursive wavelet transform.
As shown in Fig. 2, when the wavelet transform is performed on the image, it will be decomposed into four parts (subimages): the approximate part of the image is denoted as A; the part with high frequency in the vertical direction and low frequency in the horizontal direction is denoted as H; the part with high frequency in the horizontal direction and low frequency in the vertical direction is denoted V ; the part with high frequency in both the horizontal and vertical directions is denoted as D. Assuming the gridlines are aligned horizontally, the proposed method outline has the following steps: 1. The image with gridline artifacts is first decomposed into four subimages: A, H, V , and D using a wavelet transform. According to the estimated gridline frequency in the image, A or H or both of them will continue to be decomposed. The decomposition operation on the subimages is the same as that on the original image. 2. The recursive decomposition is stopped when one or several subimages that contain mainly the gridline signals are detected using a gridline detection module. 3. An adaptive Gaussian band-stop filter is applied to the detected gridline dominant subimages, and then they are replaced by their filtered ones. 4. Finally, an inverse wavelet transform is applied to the wavelet pyramid and a restored image without gridline artifacts is obtained. 3.B. Magnification of the gridline artifacts signal 3.B.1. Recursive wavelet decomposition
The 2D discrete wavelet transform (DWT) is a useful tool for multiscale image analysis and can be used as a Medical Physics, Vol. 42, No. 4, April 2015
“mathematical signal magnifier.”24 We adapted the recursive wavelet decomposition framework to magnify the gridline artifacts signal into one of the subimages. The DWT can be realized using digital filtering and a down-sampling operation.25 An M × N image can be considered to be M column vectors of N elements or N row vectors of M elements. First, we apply a high-pass filter hψ (−n) and a low-pass filter hϕ (−n) to each row of the image and then downsample the two filtered images by removing every other column. Thus, we get two subimages, one is the low-pass filtered M × N/2 subimage and another is the high-pass filtered M × N/2 subimage. A high-pass filter hψ (−m) and a lowpass filter hϕ (−m) are then applied to the each column of the two downsampled images, respectively, and then the resulted four images are downsampled by removing every other row. Thus, the image is decomposed into four M/2 × N/2 parts (subimages): A, H, V , and D. The image can be reconstructed from these four parts using the inverse transform. Orthogonal wavelet bases are often used in the discrete wavelet transform. In our framework, we will use the Daubechies wavelets,26 which can be realized using fast wavelet transform algorithms. The Daubechies wavelets form a family of wavelet bases, any one of which can be used as the basis for a discrete wavelet transform. The wavelets differ by their level, each with a different frequency response. A wavelet of Level L is denoted by dbL. Figure 3 presents the frequency response of the low pass filter (the “scaling function”) and the high pass filter (the “wavelet function”) of the db10 wavelet. The frequency locations where the curves start to drop down in the low-pass filter and rise up to the maximum magnitude in the high-pass filter are denoted as r 1 and r 2, respectively. As mentioned previously, after one-step decomposition, we obtain four wavelet coefficients’ subimages: A, H, V , and
1724
Tang et al.: Gridline suppression
1724
F. 3. The frequency response function of the low-pass filter and the high-pass filter of discrete wavelet transform db10. The horizontal axis is the ratio between the spatial frequency and the Nyquist frequency.
D. Assuming that the gridline appears horizontally in the original image, because of this directional specificity, we can predict that the gridline component will appear in either A or H according to the gridline frequency in the image. Because the filters are smooth but not ideal (see Fig. 3), the gridline information can also appear in both A and H if its frequency is located in the interval (r 1,r 2). We can estimate the gridline frequency in the image according to Eqs. (2) or (3). During the gridline components localization process (as shown in Fig. 2), only the subimage or subimages that contain gridline signals will be placed in the further decomposition pipeline. DWT is used as a tool to magnify the gridline signal. Based on the above analysis, the recursive wavelet decomposition process can be summarized as follows: 1. On the input image I, calculate the sampling frequency f s according to the pixel resolution. Implement the discrete wavelet transform on I and obtain four subimages: A, H, V , and D; 2. Check if f s satisfies the Nyquist sampling frequency, if f s ≥ 2 f g , calculate the estimated gridline frequency f estimate in the image according to Eq. (2), otherwise calculate it according to Eq. (3); 3. If f estimate ∈ [0,r 1], set I = A and repeat step 1; else if f estimate ∈ [r 2,1], set I = H and repeat step 1; else if f estimate ∈ (r 1,r 2), set I = A and repeat step 1, then set I = H and repeat step 1. Note that in step 3, when the input image in step 1 is replaced by A or H or both, because of the down-sampling operation, the pixel resolution changes so that the sampling frequency f s varies accordingly. 3.B.2. Gridline detection module
The recursive wavelet decomposition process is performed until a subimage, which contains mainly the gridline signal, is found by a gridline detection module. The gridline artifacts are typically line-shaped shadows. We can use some statistical texture features to quantitatively decide if the gridline is the main signal in this subimage. Haralick27 summarized fourteen statistical texture features based on the image gray Medical Physics, Vol. 42, No. 4, April 2015
level co-occurrence matrix (GLCM). But these features have information redundancy. Ulaby et al.28 found that among the fourteen features, only four are uncorrelated (correlation, contrast, energy, and inverse moment). We choose correlation and contrast to describe the signal strength of the gridline artifacts. From Haralick’s definition,27 the spatial GLCM is created by calculating how often a pixel with the intensity (gray-level) value i occurs in a specific spatial relationship to a pixel with the value j. Each element (i, j) in the resultant GLCM is simply the sum of the number of times that the pixel with value i occurred in the specified spatial relationship to a pixel with value j in the input image. In order to describe the spatial relationship, a distance d and the angular relationship are needed. Assume the gridline artifacts are in parallel with axis in the image, we select two angular relationships, 0◦ and 90◦, to calculate the corresponding spatial GLCM, denoted as P(d,0◦) and P(d,90◦), respectively. We call P(d,0◦) as horizontal GLCM and P(d,90◦) as vertical GLCM. The element of the matrix is denoted as pij. Assume the intensity range of the image is [0,G − 1], from the corresponding GLCM, we can calculate the correlation which describes the pixel correlation in the certain direction of the image. It is calculated by G−1 G−1 pij (i − µi ) j − µ j , (5) Correlation = σi σ j i=0 j=0 where µi =
G−1 G−1
(pij ×i),
i=0 j=0
µj =
G−1 G−1
(pij × j),
i=0 j=0
σi =
G−1 G−1 (
) pij × (i − µi )2 ,
i=0 j=0
σj =
G−1 G−1 ( i=0 j=0
pij × j − µ j
2)
.
1725
Tang et al.: Gridline suppression
1725
Contrast reflects the degree of image clarity and the texture groove depth in the direction. It is calculated by Contrast =
G−1 G−1
pij|i − j|2.
(6)
i=0 j=0
Assuming that the gridline appears horizontally in the image, we perform gridline detection using the following steps: 1. Let d = 2, calculate the image horizontal GLCM P(d,0◦) and vertical GLCM P(d,90◦), respectively; 2. Calculate the following parameters from the two GLCMs: the correlation from P(d,0◦) according to Eq. (5), denoted as hCor; the correlation from P(d,90◦) according to Eq. (5), denoted as vCor; the contrast from P(d,90◦) according to Eq. (6), denoted as vCon. 3. Let Dcor = |hCor − vCor|, if Dcor > thr1 and vCon > thr2, we can make the decision that the main content of this image is composed of horizontal lines. The two threshold thr1 and thr2 are chosen according to the image grey value range. If the image grey value range is normalized to [0,1], we suggest thr1 = 0.2 and thr2 = 0.1775, respectively. 3.C. Gaussian band-stop filter
Assume that the grids are aligned horizontally, for a column in the subimage containing mainly the gridline signal Ic (x, y), a Gaussian band-stop filter is applied to remove the gridline signal. It is designed as follows: 2
B(u) = 1 − e−(1/2)((u−µ)/σ) , u = 1,...,M,
(7)
where µ = f u′ is the gridline’s actual frequency in this subimage. The value µ is chosen by finding the maximum in the power spectrum near the estimated gridline frequency f estimate. The standard deviation σ is computed from the interval f u′ − 0.5 ∗W, f u′ + 0.5 ∗W , in which W is twice or some other multiple of the wave peak width. After applying the Gaussian band-stop filter to all the columns of Ic (x, y), we obtain an image Ic′ (x, y) free of gridline artifacts. Ic (x, y) is then replaced by Ic′ (x, y) in the pyramid and an inverse discrete wavelet transform (IDWT) is applied to obtain a restored gridline free image as shown in Fig. 2. 3.D. Algorithm summary
Based on the above analysis and equations, the gridline artifact suppression method can be summarized as follows. As inputs, we have the x-ray image with gridline artifacts I (x, y), the number of lines per centimeter in the stationary grid f g , the pixel resolution Rs , and the gridline artifacts direction in image d. Step 1: Check the gridline artifacts direction in the image, if it is not horizontal, rotate the image until d → d ′, where d ′ denotes the horizontal direction; Medical Physics, Vol. 42, No. 4, April 2015
Step 2: Let f s = 1/Rs , calculate f estimate according to Eqs. (2) or (3); Step 3: Implement the 2D DWT on I (x, y) and obtain four components: Ai , Hi , Vi , and Di , where i stands for the scale of the decomposition; Step 4: Estimate the subimages that could contain the gridline signal according to Sec. 3.B.1 and denote it or them as G i , which could be Ai or Hi or both; Step 5: Calculate the texture features of G i according to Eqs. (5) and (6) and make the decision if G i contains mainly the gridline according to Sec. 3.B.2, if not, perform a recursive DWT decomposition, that is, let G i become the input image and go back to step 3, i = i + 1; if yes, go to step 6; Step 6: Implement a Gaussian band-stop filter on G i and obtain a gridline free subimage G i′, the detail of the filter is described in Sec. 3.C; Step 7: Replace G i with G i′ and perform a IDWT to obtain the restored gridline free image I ′ (x, y).
4. EXPERIMENTAL RESULTS AND DISCUSSION We evaluated the proposed method on several phantom images, taken with different stationary grids, provided by Siemens Shanghai Medical Equipment Ltd., Shanghai, China. The image receptor is Varian 4336w CsI (Varian medical systems, Palo Alto, CA) with a linear dose response. All the results were checked and qualified manually by an image expert from Siemens Shanghai Medical Equipment Ltd. In this section, we will provide an illustration of our method and compare it with some existing methods. Figure 4(a) shows an image of a chest phantom. We processed the For Processing image produced by the detector; however, in the illustrations presented here, we converted the images to a logarithmic dose response and applied an inverse lookup table to match the appearance of the clinical display. Figure 4(b) shows a selected region in Fig. 4(a). Figure 4(c) shows the corresponding portion in which the artifacts are eliminated. The image size was 2560 × 3072 pixels with a resolution of 0.139 mm/pixel. A 50-lines/cm grid was used while acquiring the image. Thus, f g was 5 lines/mm (50 lines/cm) and f s was 7.194 pixels/mm. Because f s < 2 f g , the sampling frequency did not satisfy the Nyquist sampling theorem. Therefore, the grid texture in the image was subject to aliasing. The estimated gridline artifacts frequency was estimated by Eq. (3): k = 1, we obtained f estimate = 2.194 lines/mm. For each row with 2560 elements, when we performed the fast Fourier transform (FFT), the length was expanded to N = 210 = 4096. So the estimated artifacts frequency was located around the position N∗ f estimate/ f s = 1249 on the frequency axis. Because the FFT spectrum is symmetrical, we only show the positive frequencies’ half axis in the following illustration. Figure 4(f) shows the 1D Fourier spectrum of the sum of several rows (No. 90-100) of the original image. We notice that
1726
Tang et al.: Gridline suppression
1726
F. 4. (a) A chest phantom. Logarithmic and inverse operation was performed on the image’s intensity to simulate the real clinical DR output; (b) The portion that is highlighted in white in (a); (c) The grid pattern was removed using the proposed method; (d) The gridline-mainly component calculated according to Sec. 3.B.2; (e) After applying a Gaussian band-stop filter on the image shown in (d); (f) The spectrum of a 1D Fourier transform of the 90-100 rows of the image shown in (a); (g) The spectrum after grid artifacts removed.
the actual gridline frequency peak is close to the estimated position so that it can be located by looking for a crest around the theoretical position. After applying a first DWT on the image, the sample frequency and the FFT length N needed to be divided by two. Figure 4(d) shows one of the subimages Medical Physics, Vol. 42, No. 4, April 2015
that contain mainly the gridline signal that was detected by the module described in Sec. 3.B.2. When we implemented the Gaussian band-stop filter expressed in Sec. 3.C on this subimage, we obtained Fig. 4(e). The filtered subimage was set into the pyramid and after applying an IDWT, we obtained
1727
Tang et al.: Gridline suppression
the restored gridline free image, as shown in Fig. 4(c). The 1D Fourier spectrum of the restored image in Fig. 4(g) helps to verify the proposed gridline suppression method. To compare our proposed method against some previously published artifact removal methods, we used a phantom containing lines and numbers. The original image was 3072 × 2560 pixels with resolution 0.139 mm/pixel. A 40-lines/cm grid was used while acquiring the image. We used nonclinical images because Lin et al.18 reported that using geometricshaped characters is the best way to demonstrate the effect of the applied methods. We compared our method to the frequency domain filtering (Gaussian band-stop filter) method proposed by Lin et al.18 and to the wavelet domain suppression 2D method proposed by Sasada et al.22 In Fig. 5(a), a zoomed part of the original image with the gridline artifacts can be seen. The result of the several gridline artifact removal methods can be seen on Fig. 5(b), our proposed method; Fig. 5(c), Lin’s Gaussian band-stop filter method; and Fig. 5(d), the Sasada’s 2D method. From visual inspection, both our proposed method and that of Lin’s obtained a gridlinefree image while maintaining the image details. However in Fig. 5(d), we clearly observe the ringing effect along the gridline direction in the resulting image obtained using Sasada’s method. The goal of the proposed method is to maximally preserve the information of the target object during the gridline removal process. From visual inspection of the image, we cannot
1727
clearly see the information preservation effect. To explore this advantage of our proposed method, the corresponding frequency spectrums of the four images in Fig. 5 are shown in Fig. 6. The black circle highlights the region close to the gridline information area. The highlighted regions are magnified simultaneously. The corresponding data loss after processing is also shown. We can see that compared with the other two methods, the data loss for the proposed method was the lowest so that it had the best information preservation after gridline artifacts were removed. The Lin’s method preserves some useful information. In Lin’s paper, the power spectrum is logarithmic so that the difference before and after gridline suppression is not so significant as we show here. The 2D method lost most of the information close to the grid artifacts frequency because of the replacement of the gridline subimages by zeros. The above methods were implemented on a PC with an i5 (2.80 GHz) CPU running the Windows 7 operating system (Microsoft, Redmond WA). To compare the execution time of the three previously mentioned methods, we processed several images of 3072 × 2560 pixels. The average execution times are given in Table I. We can see that Sasada’s 2D wavelet transform method was the fastest. The reason for this is that it assigns zero directly to the considered gridline containing subimages without any analysis. The direct Gaussian bandstop filter method was the slowest because it applied the relatively time consuming filter process on the original
F. 5. (a) A zoomed part of the original image with gridline artifacts. The gridline artifacts were removed using (b) our proposed method; (c) the Gaussian band-stop filter method proposed by Lin et al. (Ref. 18) and (d) the 2D method proposed by Sasada et al. (Ref. 22). Medical Physics, Vol. 42, No. 4, April 2015
1728
Tang et al.: Gridline suppression
1728
F. 6. The spectrum of a 1D Fourier transform of the CDRAD phantom image, the highlighted portion and the corresponding data loss after the grid artifacts were removed. (a) The power spectrum of the image with gridline artifacts; (b) The power spectrum after grid artifacts were removed by the proposed method; (c) The power spectrum after grid artifacts were removed using Lin’s method; (d) The power spectrum after grid artifacts were removed using Sasada’s method.
high resolution image. Our proposed method achieved a relatively high computation speed because the time consuming Gaussian band-stop filter was applied to a subimage with a reduced resolution compared with the original image one. Consequently, it was almost three times faster than the direct
T I. The average execution time of the three methods. Methods Gaussian band-stop filter method Proposed method 2D method
Medical Physics, Vol. 42, No. 4, April 2015
Average execution time (s) 12 4.3 2.3
Gaussian band-stop filter method but two times slower than Sasada’s zero-assign method. 5. CONCLUSIONS In this paper, a new automatic method was presented to remove the gridline artifacts while preserving the maximum useful information. It used the discrete wavelet transform to magnify the gridline signal in the decomposed subimages. Several subimages that contained mainly the gridline artifacts signal were detected by a detection module, which also determined at which level to stop the decomposition process. An automatic Gaussian band-stop filter was then applied to these detected subimages to remove the gridline artifacts
1729
Tang et al.: Gridline suppression
signals. Finally, the gridline free image was obtained using the inverse discrete wavelet transform. From the experimental results, we can see that the proposed method achieved a gridline-free image while preserving the signals with similar frequency to the gridline artifacts. It also achieved a good balance between the suppression effect and the computation time. We hope and believe that this framework can also be adapted to other applications.29,30 ACKNOWLEDGMENT This work was partly supported by Siemens Shanghai Medical Equipment Ltd. All the original data for the experiment were supplied by Siemens Shanghai Medical Equipment Ltd.
a)Author
to whom correspondence should be addressed. Electronic mail: [email protected] 1N. Tanaka, K. Naka, A. Saito, J. Morishita, F. Toyofuku, M. Ohki, and Y. Higashida, “Investigation of optimum anti-scatter grid selection for digital radiography: Physical imaging properties and detectability of low-contrast signals,” Radiol. Phys. Technol. 6, 54–60 (2013). 2G. Gennaro, L. Katz, H. Souchay, R. Klausz, C. Alberelli, and C. Maggio, “Grid removal and impact on population dose in full-field digital mammography,” Med. Phys. 34, 547–555 (2007). 3T. Lehnert, N. N. Naquib, S. Wutzler, R. W. Bauer, J. M. Kerl, T. Burkhard, B. Schulz, M. C. Larson, H. Ackermann, T. J. Vogl, and J. Balzer, “Comparative study between mobile computed radiography and mobile flat-panel radiography for bedside chest radiography: Impact of an antiscatter grid on the visibility of selected diagnostically relevant structures,” Invest. Radiol. 49, 1–6 (2014). 4W. J. Vedkamp, M. A. Thijssen, and N. Karssemeijer, “The value of scatter removal by a grid in full field digital mammography,” Med. Phys. 30, 1712–1718 (2003). 5C. C. Shaw, T. Wang, and D. Gur, “Effectiveness of antiscatter grids in digital radiography. A phantom study,” Invest. Radiol. 29, 636–642 (1994). 6K. A. Fetterly and B. A. Schueler, “Physical evaluation of prototype highperformance anti-scatter grids: Potential for improved digital radiographic image quality,” Phys. Med. Biol. 54, N37–N42 (2009). 7J. A. Sorenson, L. T. Niklason, S. C. Jacobsen, D. F. Knutti, and T. Johnson, “Tantalum air-interspace crossed grid: Design and performance characteristics,” Radiology 145, 485–492 (1982). 8J. D. Robinson, D. Ferlic, F. Kotula, L. Ferlic, R. A. Geise, and K. Amplatz, “Improved mammography with a reduced radiation dose,” Radiology 188, 868–871 (1993). 9D. M. Gauntt and G. T. Barnes, “Grid line artifact formation: A comprehensive theory,” Med. Phys. 33, 1668–1677 (2006). 10D. M. Gauntt and G. T. Barnes, “A novel technique to suppress grid line artifacts,” Med. Phys. 33, 1654–1667 (2006). 11J. W. Yoon, Y. G. Park, C. J. Park, D. I. Kim, J. H. Lee, N. K. Chung, B. Y. Choe, T. S. Suh, and H. Lee, “Reduction of a grid moiré pattern by
Medical Physics, Vol. 42, No. 4, April 2015
1729 integrating a carbon-interspaced high precision x-ray grid with a digital radiographic detector,” Med. Phys. 34, 4092–4097 (2007). 12L. L. Barski and X. Wang, “Characterization, detection, and suppression of stationary grids in digital projection radiography imagery,” Proc. SPIE 3658, 502–519 (1999). 13T. Maruyama and H. Yamamoto, “Elimination of gridlines in x-ray imaging for mammography,” Imagining Systems and Techniques, Minori, Italy, 29 April 2006 (IEEE, New York, NY, 2006), pp. 104–109. 14T. Maruyama and H. Yamamoto, “Elimination of gridlines in x-ray image,” Instrumentation and Measurement Technology Conference Proceedings (IEEE, Victoria, BC, 2008), pp. 1091–1096. 15T. Maruyama and H. Yamamoto, “Elimination of gridlines by using nonlinear filter in mammographic image,” IET Image Process. 5, 457–465 (2011). 16L. L. Barski and X. Wang, “Method for x-ray antiscatter grid detection and suppression in digital radiography,” U.S. patent No. US 6269176B1 (31 July 2001). 17I. N. Belykh and C. Cornelius, “Method for antiscatter stationary grid artifacts detection and attenuation in digital radiographic images,” U.S. patent No. US 7050618B2 (23 May 2006). 18C. Y. Lin, W. J. Lee, S. J. Chen, C. H. Tsai, J. H. Lee, C. H. Chang, and Y. T. Ching, “A study of grid artifacts formation and elimination in computed radiographic images,” J. Digital Imaging 19, 351–361 (2006). 19D. S. Kim and S. Lee, “Adaptive grid pattern artifact reduction in radiography imaging based on the significant-signal bandwidth,” IEEE International Conference on Image Processings, Brussels, Belgium, 2011, (IEEE, New York, NY, 2011), pp. 1473–1476. 20D. S. Kim, S. Lee, and J. K. Yoon, “Grid artifact reduction based on homomorphic filtering in digital radiography imaging,” Proc. SPIE 8688, 86682C (2013). 21H. Kaisei-machi, “Method and unit for suppressing a periodic pattern,” U.S. patent No. US 7336811B2 (Jan. 22, 2011). 22R. Sasada, M. Yamada, S. Hara, H. Takeo, and K. Shimura, “Stationary grid pattern removal using 2-dimensional technique for moiré-free radiographic image display,” Proc. SPIE 5029, 688 (2003). 23D. Linden, “A discussion of sampling theorem,” IEEE Trans. Pattern Anal. Mach. Intell. 47, 1219 (1959). 24S. G. Mallat, “A theory for multiresolution signal decomposition: The wavelet representation,” IEEE Trans. Pattern Anal. Mach. Intell. 11, 674–693 (1989). 25R. C. Gonzalez and R. Woods, Digital image processing, 2nd ed. (Prentice Hall, Upper Saddle River, NJ, 2002), pp. 124–322. 26I. Daubechies, “Orthonormal bases of compactly supported wavelets,” Commun. Pure Appl. Math. 41, 909–996 (1988). 27R. Haralick, “Statistical and structural approaches to texture,” Proc. IEEE 67, 786–804 (1979). 28F. T. Ulaby, F. Kouyate, B. Brisco, and T. Williams, “Textural infornation in SAR images,” IEEE Trans. Geosci. Remote Sens. GE-24, 235–245 (1986). 29Y. Chen, L. Shi, Q. Feng, J. Yang, H. Shu, L. Luo, J. Cortrieux, and W. Chen, “Artifact suppressed dictionary learning for low-dose CT image processing,” IEEE Trans. Med. Imaging 33, 2271–2292 (2014). 30Y. Chen, Y. Zhou, Y. Hu, G. Yang, L. Luo, W. Chen, and C. Toumoulin, “Thoracic low-dose CT image processing using an artifact suppressed largescale nonlocal means,” Phys. Med. Biol. 57, 2667–2688 (2012).
A new stationary gridline artifact suppression method based on the 2D discrete wavelet transform Hui Tanga) Laboratory of Image Science and Technology, School of Computer Science and Engineering, Southeast University, Nanjing 210096, China; Key Laboratory of Computer Network and Information Integration (Southeast University), Ministry of Education, Nanjing 210000, China; Centre de Recherche en Information Biomédicale sino-français, Laboratoire International Associé, Inserm, Université de Rennes 1, Rennes 35000, France; and Southeast University, Nanjing 210000, China
Dan Tong and Xu Dong Bao Laboratory of Image Science and Technology, School of Computer Science and Engineering, Southeast University, Nanjing 210096, China
Jean-Louis Dillenseger INSERM, U1099, Rennes F-35000, France; Université de Rennes 1, LTSI, Rennes F-35000, France; Centre de Recherche en Information Biomédicale sino-français, Laboratoire International Associé, Inserm, Université de Rennes 1, Rennes 35000, France; and Southeast University, Nanjing 210000, China
(Received 28 August 2014; revised 12 February 2015; accepted for publication 17 February 2015; published 18 March 2015) Purpose: In digital x-ray radiography, an antiscatter grid is inserted between the patient and the image receptor to reduce scattered radiation. If the antiscatter grid is used in a stationary way, gridline artifacts will appear in the final image. In most of the gridline removal image processing methods, the useful information with spatial frequencies close to that of the gridline is usually lost or degraded. In this study, a new stationary gridline suppression method is designed to preserve more of the useful information. Methods: The method is as follows. The input image is first recursively decomposed into several smaller subimages using a multiscale 2D discrete wavelet transform. The decomposition process stops when the gridline signal is found to be greater than a threshold in one or several of these subimages using a gridline detection module. An automatic Gaussian band-stop filter is then applied to the detected subimages to remove the gridline signal. Finally, the restored image is achieved using the corresponding 2D inverse discrete wavelet transform. Results: The processed images show that the proposed method can remove the gridline signal efficiently while maintaining the image details. The spectra of a 1D Fourier transform of the processed images demonstrate that, compared with some existing gridline removal methods, the proposed method has better information preservation after the removal of the gridline artifacts. Additionally, the performance speed is relatively high. Conclusions: The experimental results demonstrate the efficiency of the proposed method. Compared with some existing gridline removal methods, the proposed method can preserve more information within an acceptable execution time. C 2015 American Association of Physicists in Medicine. [http://dx.doi.org/10.1118/1.4914861] Key words: x-ray image, stationary grid, wavelet transform, gridline suppression, Gaussian filter 1. INTRODUCTION X-ray radiography is one of the most common imaging modalities in the medical field for diagnosis. Compared with analog radiographic systems that use screen/film combinations, digital radiographic (DR) systems directly produce digital images with optimal density and brightness for interpretation or diagnosis, because the detected signal can be modified by electronic adjustment of the system gain.1 When the x-ray photons strike the object being imaged, some photons are scattered, which leads to contrast degradation and reduced image quality. For both analog and digital systems, the most commonly used technique for reducing the scattered radiation is to insert an antiscatter grid between the patient and the image receptor.2 The use of antiscatter 1721
Med. Phys. 42 (4), April 2015
grids has been proven to increase the image contrast and the signal-to-noise ratio in radiographic images.3–6 Two types of grid, stationary or moving, are used according to the situation. Compared with stationary grids, systems with moving grids require a more complex system architecture. If a grid is stationary during the exposure, the grid shadows may be visible in the image and are known as gridline artifacts. Numerous approaches have been undertaken to suppress the gridline artifacts in x-ray images.7–18 Some researchers have reported methods that avoid gridline artifacts by adjusting some of the physical modules during the x-ray projection and detection process in the imaging system.7–11 These methods are limited in some specific situations because of the induced system architecture complexity.
0094-2405/2015/42(4)/1721/9/$30.00
© 2015 Am. Assoc. Phys. Med.
1721
1722
Tang et al.: Gridline suppression
1722
Comparatively, image processing methods are more adaptive.12–22 Stationary gridline artifacts are expressed as regularly distributed periodic signals in the x-ray images. The different image processing based artifact suppression methods can be classified into three categories: 1. Space domain methods.12–15 In this class of methods, the gridline artifacts are analyzed according to the image gray level information and suppressed in the spatial domain. The disadvantage of these methods is that the gridline detection process is usually affected by the object image information and the convolution gridline suppression process will blur the object image. 2. Frequency domain filtering methods.16–20 In the frequency domain, the gridline signals only appear at certain frequencies related to the image acquisition situation. With prior knowledge, the gridline artifacts and the imaging object are relatively easy to separate in the frequency domain. But when the grid frequency is suppressed in the frequency domain, the object image information with a similar frequency will also be reduced. 3. Wavelet domain suppression methods.21,22 In this method, a recursive wavelet transform is applied to the image until subimages of the desired resolution are generated. The subimages that are expected to contain the gridline artifacts are set to zero, and finally the restored image is obtained using the inverse wavelet transform. These methods are proposed to preserve more of the object information, but the disadvantage is that these methods did not provide a stop condition for the recursive wavelet transform and the zero assignment operation will lead to a ringing effect in the final image. To efficiently remove gridline artifacts from x-ray images while preserving object information, we propose a new gridline artifact suppression method. First, the image with gridline artifacts is decomposed into several subimages using a recursive wavelet transform. The decomposition is stopped when the gridline signals are detected to be dominant in one or several of the subimages using an adaptive gridline detection module. An adaptive Gaussian band-stop filter is then applied to this or these subimages, which is or are replaced by the filtered version(s). Finally, an inverse wavelet transform is applied to the wavelet pyramid to obtain a restored image without gridline artifacts. 2. GRIDLINE ARTIFACT EXPRESSION IN AN X-RAY IMAGE In the stationary grid case, the x-rays emitted by the source pass through and are thus absorbed by the target object and the stationary antiscatter grid. Both will leave shadows on the final image. The grid’s ideal shadow is show in Fig. 1, denoted as Sgrid (x), multiplies the desired image on the final image. It can be expressed using the Fourier expansion Sgrid (x) =
d Tgrid
+
∞
(ai cos(2iπx f g )).
i=1
Medical Physics, Vol. 42, No. 4, April 2015
(1)
F. 1. Ideal gridline signal in the space domain.
According to Lin’s analysis,18 from this formula, we can get d/Tgrid, the DC term; hg 1 (x) = a1 cos(2πx f g ), the fundamental wave with the frequency f g ; and hg m (x), m = 2, ..., ∞, the harmonics with the frequency m f g for each term. Since the harmonics have a higher frequency but a lower magnitude than the fundamental wave, the effect of the harmonics can be considered small and can be ignored. For the following discussion, we will only take the fundamental wave into account. Equation (1) describes the gridline artifacts in the presampled image. However, the digital image is also affected by the sampling function. We will denote the sampling frequency as f s . According to the Nyquist sampling theorem, if f s ≥ 2 f g , the data at f g can be reconstructed without aliasing. Otherwise, aliasing occurs.23 So according to the relationship between the sampling frequency and the grid frequency, we can estimate the grid artifact frequency, denoted as f estimate, by if f s ≥ 2 f g , f estimate = f g ,
(2)
else, f estimate = f g − k f s , (3) ( ) ( ) fg fg if mod < 0.5 floor fs fs ( ) , (4) k≡ fg otherwise ceil fs where floor(x) is x rounded down to the nearest integer and ceil(x) is x rounded up to the nearest integer. 3. PROPOSED METHOD 3.A. Algorithm outline
The goal of our work is to remove grid patterns from the x-ray image while preserving the object information. Wavelet transformation is a useful tool to separate signals, so we endeavor to implement this tool to focus the image processing on only the gridline artifacts’ part of the signal instead of the whole image.
1723
Tang et al.: Gridline suppression
1723
F. 2. Gridline artifact suppression method based on the 2D recursive wavelet transform.
As shown in Fig. 2, when the wavelet transform is performed on the image, it will be decomposed into four parts (subimages): the approximate part of the image is denoted as A; the part with high frequency in the vertical direction and low frequency in the horizontal direction is denoted as H; the part with high frequency in the horizontal direction and low frequency in the vertical direction is denoted V ; the part with high frequency in both the horizontal and vertical directions is denoted as D. Assuming the gridlines are aligned horizontally, the proposed method outline has the following steps: 1. The image with gridline artifacts is first decomposed into four subimages: A, H, V , and D using a wavelet transform. According to the estimated gridline frequency in the image, A or H or both of them will continue to be decomposed. The decomposition operation on the subimages is the same as that on the original image. 2. The recursive decomposition is stopped when one or several subimages that contain mainly the gridline signals are detected using a gridline detection module. 3. An adaptive Gaussian band-stop filter is applied to the detected gridline dominant subimages, and then they are replaced by their filtered ones. 4. Finally, an inverse wavelet transform is applied to the wavelet pyramid and a restored image without gridline artifacts is obtained. 3.B. Magnification of the gridline artifacts signal 3.B.1. Recursive wavelet decomposition
The 2D discrete wavelet transform (DWT) is a useful tool for multiscale image analysis and can be used as a Medical Physics, Vol. 42, No. 4, April 2015
“mathematical signal magnifier.”24 We adapted the recursive wavelet decomposition framework to magnify the gridline artifacts signal into one of the subimages. The DWT can be realized using digital filtering and a down-sampling operation.25 An M × N image can be considered to be M column vectors of N elements or N row vectors of M elements. First, we apply a high-pass filter hψ (−n) and a low-pass filter hϕ (−n) to each row of the image and then downsample the two filtered images by removing every other column. Thus, we get two subimages, one is the low-pass filtered M × N/2 subimage and another is the high-pass filtered M × N/2 subimage. A high-pass filter hψ (−m) and a lowpass filter hϕ (−m) are then applied to the each column of the two downsampled images, respectively, and then the resulted four images are downsampled by removing every other row. Thus, the image is decomposed into four M/2 × N/2 parts (subimages): A, H, V , and D. The image can be reconstructed from these four parts using the inverse transform. Orthogonal wavelet bases are often used in the discrete wavelet transform. In our framework, we will use the Daubechies wavelets,26 which can be realized using fast wavelet transform algorithms. The Daubechies wavelets form a family of wavelet bases, any one of which can be used as the basis for a discrete wavelet transform. The wavelets differ by their level, each with a different frequency response. A wavelet of Level L is denoted by dbL. Figure 3 presents the frequency response of the low pass filter (the “scaling function”) and the high pass filter (the “wavelet function”) of the db10 wavelet. The frequency locations where the curves start to drop down in the low-pass filter and rise up to the maximum magnitude in the high-pass filter are denoted as r 1 and r 2, respectively. As mentioned previously, after one-step decomposition, we obtain four wavelet coefficients’ subimages: A, H, V , and
1724
Tang et al.: Gridline suppression
1724
F. 3. The frequency response function of the low-pass filter and the high-pass filter of discrete wavelet transform db10. The horizontal axis is the ratio between the spatial frequency and the Nyquist frequency.
D. Assuming that the gridline appears horizontally in the original image, because of this directional specificity, we can predict that the gridline component will appear in either A or H according to the gridline frequency in the image. Because the filters are smooth but not ideal (see Fig. 3), the gridline information can also appear in both A and H if its frequency is located in the interval (r 1,r 2). We can estimate the gridline frequency in the image according to Eqs. (2) or (3). During the gridline components localization process (as shown in Fig. 2), only the subimage or subimages that contain gridline signals will be placed in the further decomposition pipeline. DWT is used as a tool to magnify the gridline signal. Based on the above analysis, the recursive wavelet decomposition process can be summarized as follows: 1. On the input image I, calculate the sampling frequency f s according to the pixel resolution. Implement the discrete wavelet transform on I and obtain four subimages: A, H, V , and D; 2. Check if f s satisfies the Nyquist sampling frequency, if f s ≥ 2 f g , calculate the estimated gridline frequency f estimate in the image according to Eq. (2), otherwise calculate it according to Eq. (3); 3. If f estimate ∈ [0,r 1], set I = A and repeat step 1; else if f estimate ∈ [r 2,1], set I = H and repeat step 1; else if f estimate ∈ (r 1,r 2), set I = A and repeat step 1, then set I = H and repeat step 1. Note that in step 3, when the input image in step 1 is replaced by A or H or both, because of the down-sampling operation, the pixel resolution changes so that the sampling frequency f s varies accordingly. 3.B.2. Gridline detection module
The recursive wavelet decomposition process is performed until a subimage, which contains mainly the gridline signal, is found by a gridline detection module. The gridline artifacts are typically line-shaped shadows. We can use some statistical texture features to quantitatively decide if the gridline is the main signal in this subimage. Haralick27 summarized fourteen statistical texture features based on the image gray Medical Physics, Vol. 42, No. 4, April 2015
level co-occurrence matrix (GLCM). But these features have information redundancy. Ulaby et al.28 found that among the fourteen features, only four are uncorrelated (correlation, contrast, energy, and inverse moment). We choose correlation and contrast to describe the signal strength of the gridline artifacts. From Haralick’s definition,27 the spatial GLCM is created by calculating how often a pixel with the intensity (gray-level) value i occurs in a specific spatial relationship to a pixel with the value j. Each element (i, j) in the resultant GLCM is simply the sum of the number of times that the pixel with value i occurred in the specified spatial relationship to a pixel with value j in the input image. In order to describe the spatial relationship, a distance d and the angular relationship are needed. Assume the gridline artifacts are in parallel with axis in the image, we select two angular relationships, 0◦ and 90◦, to calculate the corresponding spatial GLCM, denoted as P(d,0◦) and P(d,90◦), respectively. We call P(d,0◦) as horizontal GLCM and P(d,90◦) as vertical GLCM. The element of the matrix is denoted as pij. Assume the intensity range of the image is [0,G − 1], from the corresponding GLCM, we can calculate the correlation which describes the pixel correlation in the certain direction of the image. It is calculated by G−1 G−1 pij (i − µi ) j − µ j , (5) Correlation = σi σ j i=0 j=0 where µi =
G−1 G−1
(pij ×i),
i=0 j=0
µj =
G−1 G−1
(pij × j),
i=0 j=0
σi =
G−1 G−1 (
) pij × (i − µi )2 ,
i=0 j=0
σj =
G−1 G−1 ( i=0 j=0
pij × j − µ j
2)
.
1725
Tang et al.: Gridline suppression
1725
Contrast reflects the degree of image clarity and the texture groove depth in the direction. It is calculated by Contrast =
G−1 G−1
pij|i − j|2.
(6)
i=0 j=0
Assuming that the gridline appears horizontally in the image, we perform gridline detection using the following steps: 1. Let d = 2, calculate the image horizontal GLCM P(d,0◦) and vertical GLCM P(d,90◦), respectively; 2. Calculate the following parameters from the two GLCMs: the correlation from P(d,0◦) according to Eq. (5), denoted as hCor; the correlation from P(d,90◦) according to Eq. (5), denoted as vCor; the contrast from P(d,90◦) according to Eq. (6), denoted as vCon. 3. Let Dcor = |hCor − vCor|, if Dcor > thr1 and vCon > thr2, we can make the decision that the main content of this image is composed of horizontal lines. The two threshold thr1 and thr2 are chosen according to the image grey value range. If the image grey value range is normalized to [0,1], we suggest thr1 = 0.2 and thr2 = 0.1775, respectively. 3.C. Gaussian band-stop filter
Assume that the grids are aligned horizontally, for a column in the subimage containing mainly the gridline signal Ic (x, y), a Gaussian band-stop filter is applied to remove the gridline signal. It is designed as follows: 2
B(u) = 1 − e−(1/2)((u−µ)/σ) , u = 1,...,M,
(7)
where µ = f u′ is the gridline’s actual frequency in this subimage. The value µ is chosen by finding the maximum in the power spectrum near the estimated gridline frequency f estimate. The standard deviation σ is computed from the interval f u′ − 0.5 ∗W, f u′ + 0.5 ∗W , in which W is twice or some other multiple of the wave peak width. After applying the Gaussian band-stop filter to all the columns of Ic (x, y), we obtain an image Ic′ (x, y) free of gridline artifacts. Ic (x, y) is then replaced by Ic′ (x, y) in the pyramid and an inverse discrete wavelet transform (IDWT) is applied to obtain a restored gridline free image as shown in Fig. 2. 3.D. Algorithm summary
Based on the above analysis and equations, the gridline artifact suppression method can be summarized as follows. As inputs, we have the x-ray image with gridline artifacts I (x, y), the number of lines per centimeter in the stationary grid f g , the pixel resolution Rs , and the gridline artifacts direction in image d. Step 1: Check the gridline artifacts direction in the image, if it is not horizontal, rotate the image until d → d ′, where d ′ denotes the horizontal direction; Medical Physics, Vol. 42, No. 4, April 2015
Step 2: Let f s = 1/Rs , calculate f estimate according to Eqs. (2) or (3); Step 3: Implement the 2D DWT on I (x, y) and obtain four components: Ai , Hi , Vi , and Di , where i stands for the scale of the decomposition; Step 4: Estimate the subimages that could contain the gridline signal according to Sec. 3.B.1 and denote it or them as G i , which could be Ai or Hi or both; Step 5: Calculate the texture features of G i according to Eqs. (5) and (6) and make the decision if G i contains mainly the gridline according to Sec. 3.B.2, if not, perform a recursive DWT decomposition, that is, let G i become the input image and go back to step 3, i = i + 1; if yes, go to step 6; Step 6: Implement a Gaussian band-stop filter on G i and obtain a gridline free subimage G i′, the detail of the filter is described in Sec. 3.C; Step 7: Replace G i with G i′ and perform a IDWT to obtain the restored gridline free image I ′ (x, y).
4. EXPERIMENTAL RESULTS AND DISCUSSION We evaluated the proposed method on several phantom images, taken with different stationary grids, provided by Siemens Shanghai Medical Equipment Ltd., Shanghai, China. The image receptor is Varian 4336w CsI (Varian medical systems, Palo Alto, CA) with a linear dose response. All the results were checked and qualified manually by an image expert from Siemens Shanghai Medical Equipment Ltd. In this section, we will provide an illustration of our method and compare it with some existing methods. Figure 4(a) shows an image of a chest phantom. We processed the For Processing image produced by the detector; however, in the illustrations presented here, we converted the images to a logarithmic dose response and applied an inverse lookup table to match the appearance of the clinical display. Figure 4(b) shows a selected region in Fig. 4(a). Figure 4(c) shows the corresponding portion in which the artifacts are eliminated. The image size was 2560 × 3072 pixels with a resolution of 0.139 mm/pixel. A 50-lines/cm grid was used while acquiring the image. Thus, f g was 5 lines/mm (50 lines/cm) and f s was 7.194 pixels/mm. Because f s < 2 f g , the sampling frequency did not satisfy the Nyquist sampling theorem. Therefore, the grid texture in the image was subject to aliasing. The estimated gridline artifacts frequency was estimated by Eq. (3): k = 1, we obtained f estimate = 2.194 lines/mm. For each row with 2560 elements, when we performed the fast Fourier transform (FFT), the length was expanded to N = 210 = 4096. So the estimated artifacts frequency was located around the position N∗ f estimate/ f s = 1249 on the frequency axis. Because the FFT spectrum is symmetrical, we only show the positive frequencies’ half axis in the following illustration. Figure 4(f) shows the 1D Fourier spectrum of the sum of several rows (No. 90-100) of the original image. We notice that
1726
Tang et al.: Gridline suppression
1726
F. 4. (a) A chest phantom. Logarithmic and inverse operation was performed on the image’s intensity to simulate the real clinical DR output; (b) The portion that is highlighted in white in (a); (c) The grid pattern was removed using the proposed method; (d) The gridline-mainly component calculated according to Sec. 3.B.2; (e) After applying a Gaussian band-stop filter on the image shown in (d); (f) The spectrum of a 1D Fourier transform of the 90-100 rows of the image shown in (a); (g) The spectrum after grid artifacts removed.
the actual gridline frequency peak is close to the estimated position so that it can be located by looking for a crest around the theoretical position. After applying a first DWT on the image, the sample frequency and the FFT length N needed to be divided by two. Figure 4(d) shows one of the subimages Medical Physics, Vol. 42, No. 4, April 2015
that contain mainly the gridline signal that was detected by the module described in Sec. 3.B.2. When we implemented the Gaussian band-stop filter expressed in Sec. 3.C on this subimage, we obtained Fig. 4(e). The filtered subimage was set into the pyramid and after applying an IDWT, we obtained
1727
Tang et al.: Gridline suppression
the restored gridline free image, as shown in Fig. 4(c). The 1D Fourier spectrum of the restored image in Fig. 4(g) helps to verify the proposed gridline suppression method. To compare our proposed method against some previously published artifact removal methods, we used a phantom containing lines and numbers. The original image was 3072 × 2560 pixels with resolution 0.139 mm/pixel. A 40-lines/cm grid was used while acquiring the image. We used nonclinical images because Lin et al.18 reported that using geometricshaped characters is the best way to demonstrate the effect of the applied methods. We compared our method to the frequency domain filtering (Gaussian band-stop filter) method proposed by Lin et al.18 and to the wavelet domain suppression 2D method proposed by Sasada et al.22 In Fig. 5(a), a zoomed part of the original image with the gridline artifacts can be seen. The result of the several gridline artifact removal methods can be seen on Fig. 5(b), our proposed method; Fig. 5(c), Lin’s Gaussian band-stop filter method; and Fig. 5(d), the Sasada’s 2D method. From visual inspection, both our proposed method and that of Lin’s obtained a gridlinefree image while maintaining the image details. However in Fig. 5(d), we clearly observe the ringing effect along the gridline direction in the resulting image obtained using Sasada’s method. The goal of the proposed method is to maximally preserve the information of the target object during the gridline removal process. From visual inspection of the image, we cannot
1727
clearly see the information preservation effect. To explore this advantage of our proposed method, the corresponding frequency spectrums of the four images in Fig. 5 are shown in Fig. 6. The black circle highlights the region close to the gridline information area. The highlighted regions are magnified simultaneously. The corresponding data loss after processing is also shown. We can see that compared with the other two methods, the data loss for the proposed method was the lowest so that it had the best information preservation after gridline artifacts were removed. The Lin’s method preserves some useful information. In Lin’s paper, the power spectrum is logarithmic so that the difference before and after gridline suppression is not so significant as we show here. The 2D method lost most of the information close to the grid artifacts frequency because of the replacement of the gridline subimages by zeros. The above methods were implemented on a PC with an i5 (2.80 GHz) CPU running the Windows 7 operating system (Microsoft, Redmond WA). To compare the execution time of the three previously mentioned methods, we processed several images of 3072 × 2560 pixels. The average execution times are given in Table I. We can see that Sasada’s 2D wavelet transform method was the fastest. The reason for this is that it assigns zero directly to the considered gridline containing subimages without any analysis. The direct Gaussian bandstop filter method was the slowest because it applied the relatively time consuming filter process on the original
F. 5. (a) A zoomed part of the original image with gridline artifacts. The gridline artifacts were removed using (b) our proposed method; (c) the Gaussian band-stop filter method proposed by Lin et al. (Ref. 18) and (d) the 2D method proposed by Sasada et al. (Ref. 22). Medical Physics, Vol. 42, No. 4, April 2015
1728
Tang et al.: Gridline suppression
1728
F. 6. The spectrum of a 1D Fourier transform of the CDRAD phantom image, the highlighted portion and the corresponding data loss after the grid artifacts were removed. (a) The power spectrum of the image with gridline artifacts; (b) The power spectrum after grid artifacts were removed by the proposed method; (c) The power spectrum after grid artifacts were removed using Lin’s method; (d) The power spectrum after grid artifacts were removed using Sasada’s method.
high resolution image. Our proposed method achieved a relatively high computation speed because the time consuming Gaussian band-stop filter was applied to a subimage with a reduced resolution compared with the original image one. Consequently, it was almost three times faster than the direct
T I. The average execution time of the three methods. Methods Gaussian band-stop filter method Proposed method 2D method
Medical Physics, Vol. 42, No. 4, April 2015
Average execution time (s) 12 4.3 2.3
Gaussian band-stop filter method but two times slower than Sasada’s zero-assign method. 5. CONCLUSIONS In this paper, a new automatic method was presented to remove the gridline artifacts while preserving the maximum useful information. It used the discrete wavelet transform to magnify the gridline signal in the decomposed subimages. Several subimages that contained mainly the gridline artifacts signal were detected by a detection module, which also determined at which level to stop the decomposition process. An automatic Gaussian band-stop filter was then applied to these detected subimages to remove the gridline artifacts
1729
Tang et al.: Gridline suppression
signals. Finally, the gridline free image was obtained using the inverse discrete wavelet transform. From the experimental results, we can see that the proposed method achieved a gridline-free image while preserving the signals with similar frequency to the gridline artifacts. It also achieved a good balance between the suppression effect and the computation time. We hope and believe that this framework can also be adapted to other applications.29,30 ACKNOWLEDGMENT This work was partly supported by Siemens Shanghai Medical Equipment Ltd. All the original data for the experiment were supplied by Siemens Shanghai Medical Equipment Ltd.
a)Author
to whom correspondence should be addressed. Electronic mail: [email protected] 1N. Tanaka, K. Naka, A. Saito, J. Morishita, F. Toyofuku, M. Ohki, and Y. Higashida, “Investigation of optimum anti-scatter grid selection for digital radiography: Physical imaging properties and detectability of low-contrast signals,” Radiol. Phys. Technol. 6, 54–60 (2013). 2G. Gennaro, L. Katz, H. Souchay, R. Klausz, C. Alberelli, and C. Maggio, “Grid removal and impact on population dose in full-field digital mammography,” Med. Phys. 34, 547–555 (2007). 3T. Lehnert, N. N. Naquib, S. Wutzler, R. W. Bauer, J. M. Kerl, T. Burkhard, B. Schulz, M. C. Larson, H. Ackermann, T. J. Vogl, and J. Balzer, “Comparative study between mobile computed radiography and mobile flat-panel radiography for bedside chest radiography: Impact of an antiscatter grid on the visibility of selected diagnostically relevant structures,” Invest. Radiol. 49, 1–6 (2014). 4W. J. Vedkamp, M. A. Thijssen, and N. Karssemeijer, “The value of scatter removal by a grid in full field digital mammography,” Med. Phys. 30, 1712–1718 (2003). 5C. C. Shaw, T. Wang, and D. Gur, “Effectiveness of antiscatter grids in digital radiography. A phantom study,” Invest. Radiol. 29, 636–642 (1994). 6K. A. Fetterly and B. A. Schueler, “Physical evaluation of prototype highperformance anti-scatter grids: Potential for improved digital radiographic image quality,” Phys. Med. Biol. 54, N37–N42 (2009). 7J. A. Sorenson, L. T. Niklason, S. C. Jacobsen, D. F. Knutti, and T. Johnson, “Tantalum air-interspace crossed grid: Design and performance characteristics,” Radiology 145, 485–492 (1982). 8J. D. Robinson, D. Ferlic, F. Kotula, L. Ferlic, R. A. Geise, and K. Amplatz, “Improved mammography with a reduced radiation dose,” Radiology 188, 868–871 (1993). 9D. M. Gauntt and G. T. Barnes, “Grid line artifact formation: A comprehensive theory,” Med. Phys. 33, 1668–1677 (2006). 10D. M. Gauntt and G. T. Barnes, “A novel technique to suppress grid line artifacts,” Med. Phys. 33, 1654–1667 (2006). 11J. W. Yoon, Y. G. Park, C. J. Park, D. I. Kim, J. H. Lee, N. K. Chung, B. Y. Choe, T. S. Suh, and H. Lee, “Reduction of a grid moiré pattern by
Medical Physics, Vol. 42, No. 4, April 2015
1729 integrating a carbon-interspaced high precision x-ray grid with a digital radiographic detector,” Med. Phys. 34, 4092–4097 (2007). 12L. L. Barski and X. Wang, “Characterization, detection, and suppression of stationary grids in digital projection radiography imagery,” Proc. SPIE 3658, 502–519 (1999). 13T. Maruyama and H. Yamamoto, “Elimination of gridlines in x-ray imaging for mammography,” Imagining Systems and Techniques, Minori, Italy, 29 April 2006 (IEEE, New York, NY, 2006), pp. 104–109. 14T. Maruyama and H. Yamamoto, “Elimination of gridlines in x-ray image,” Instrumentation and Measurement Technology Conference Proceedings (IEEE, Victoria, BC, 2008), pp. 1091–1096. 15T. Maruyama and H. Yamamoto, “Elimination of gridlines by using nonlinear filter in mammographic image,” IET Image Process. 5, 457–465 (2011). 16L. L. Barski and X. Wang, “Method for x-ray antiscatter grid detection and suppression in digital radiography,” U.S. patent No. US 6269176B1 (31 July 2001). 17I. N. Belykh and C. Cornelius, “Method for antiscatter stationary grid artifacts detection and attenuation in digital radiographic images,” U.S. patent No. US 7050618B2 (23 May 2006). 18C. Y. Lin, W. J. Lee, S. J. Chen, C. H. Tsai, J. H. Lee, C. H. Chang, and Y. T. Ching, “A study of grid artifacts formation and elimination in computed radiographic images,” J. Digital Imaging 19, 351–361 (2006). 19D. S. Kim and S. Lee, “Adaptive grid pattern artifact reduction in radiography imaging based on the significant-signal bandwidth,” IEEE International Conference on Image Processings, Brussels, Belgium, 2011, (IEEE, New York, NY, 2011), pp. 1473–1476. 20D. S. Kim, S. Lee, and J. K. Yoon, “Grid artifact reduction based on homomorphic filtering in digital radiography imaging,” Proc. SPIE 8688, 86682C (2013). 21H. Kaisei-machi, “Method and unit for suppressing a periodic pattern,” U.S. patent No. US 7336811B2 (Jan. 22, 2011). 22R. Sasada, M. Yamada, S. Hara, H. Takeo, and K. Shimura, “Stationary grid pattern removal using 2-dimensional technique for moiré-free radiographic image display,” Proc. SPIE 5029, 688 (2003). 23D. Linden, “A discussion of sampling theorem,” IEEE Trans. Pattern Anal. Mach. Intell. 47, 1219 (1959). 24S. G. Mallat, “A theory for multiresolution signal decomposition: The wavelet representation,” IEEE Trans. Pattern Anal. Mach. Intell. 11, 674–693 (1989). 25R. C. Gonzalez and R. Woods, Digital image processing, 2nd ed. (Prentice Hall, Upper Saddle River, NJ, 2002), pp. 124–322. 26I. Daubechies, “Orthonormal bases of compactly supported wavelets,” Commun. Pure Appl. Math. 41, 909–996 (1988). 27R. Haralick, “Statistical and structural approaches to texture,” Proc. IEEE 67, 786–804 (1979). 28F. T. Ulaby, F. Kouyate, B. Brisco, and T. Williams, “Textural infornation in SAR images,” IEEE Trans. Geosci. Remote Sens. GE-24, 235–245 (1986). 29Y. Chen, L. Shi, Q. Feng, J. Yang, H. Shu, L. Luo, J. Cortrieux, and W. Chen, “Artifact suppressed dictionary learning for low-dose CT image processing,” IEEE Trans. Med. Imaging 33, 2271–2292 (2014). 30Y. Chen, Y. Zhou, Y. Hu, G. Yang, L. Luo, W. Chen, and C. Toumoulin, “Thoracic low-dose CT image processing using an artifact suppressed largescale nonlocal means,” Phys. Med. Biol. 57, 2667–2688 (2012).
