|Home | About | Journals | Submit | Contact Us | Français|
In digital X-ray radiography, an anti-scatter grid is inserted between the patient and the image receptor to reduce scattered radiation. If the anti-scatter 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.
The method is as follows. The input image is first recursively decomposed into several smaller sub-images using a multi-scale 2D discrete wavelet transform (DWT). The decomposition process stops when the gridline signal is found to be greater than a threshold in one or several of these sub-images using a gridline detection module. An automatic Gaussian band-stop filter is then applied to the detected sub-images to remove the gridline signal. Finally, the restored image is achieved using the corresponding 2D inverse discrete wavelet transform (IDWT).
The processed images show that the proposed method can remove the gridline signal efficiently while maintaining the image details. The spectra of a 1-dimensional 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.
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.
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 gain1.
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 anti-scatter grid between the patient and the image receptor2. The use of anti-scatter grids has been proven to increase the image contrast and the signal-to-noise ratio in radiographic images3–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 images7–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 system7–11. These methods are limited in some specific situations because of the induced system architecture complexity.
Comparatively, image processing methods are more adaptive12–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:
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 sub-images using a recursive wavelet transform. The decomposition is stopped when the gridline signals are detected to be dominant in one or several of the sub-images using an adaptive gridline detection module. An adaptive Gaussian band-stop filter is then applied to this or these sub-images, 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.
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 anti-scatter 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:
According to Lin’s analysis18, from this formula, we can get: d/Tgrid, the DC term; hg1(x) = a1 cos(2πxfg), the fundamental wave with the frequency fg; and hgm(x),m = 2, …, ∞, the harmonics with the frequency mfg 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.
Eq. (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 fs. According to the Nyquist sampling theorem, if fs ≥ 2fg, the data at fg can be reconstructed without aliasing. Otherwise, aliasing occurs23. So according to the relationship between the sampling frequency and the grid frequency, we can estimate the grid artifact frequency, denoted as festimate, by :
If fs ≥ 2fg,
where floor(x) is x rounded down to the nearest integer, and ceil(x) is x rounded up to the nearest integer.
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.
As shown in Fig. 2, when the wavelet transform is performed on the image, it will be decomposed into four parts (sub-images): 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:
The two-dimensional discrete wavelet transform (DWT) is a useful tool for multi-scale image analysis and can be used as a “mathematical signal magnifier”24. We adapted the recursive wavelet decomposition framework to magnify the gridline artifacts signal into one of the sub-images. The DWT can be realized using digital filtering and a down-sampling operation25. 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 sub-images, one is the low-pass filtered M × N/2 sub-image and another is the high-pass filtered M × N/2 sub-image. A high-pass filter hψ(−m) and A low-pass filter h(−m) is 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 (sub-images): 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 wavelets26, 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 r1 and r2 respectively.
As mentioned previously, after one-step decomposition, we obtain four wavelet coefficients’ sub-images: A, H, V and 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 (r1, r2). We can estimate the gridline frequency in the image according to Eq. (2) or (3). During the gridline components localization process (as shown in Fig. 2), only the sub-image or sub-images 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:
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 fs varies accordingly.
The recursive wavelet decomposition process is performed until a sub-image, 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 sub-image. Haralick27 summarized fourteen statistical texture features based on the image gray 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 definition27, 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:
Contrast reflects the degree of image clarity and the texture groove depth in the direction. It is calculated by:
Assuming that the gridline appears horizontally in the image, we perform gridline detection using the following steps:
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.
Assume that the grids are aligned horizontally, for a column in the sub-image 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:
where , is the gridline’s actual frequency in this sub-image. The value μ is chosen by finding the maximum in the power spectrum near the estimated gridline frequency festimate. The standard deviation σ is computed from the interval [ ], 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 free of gridline artifacts. Ic(x,y) is then replaced by in the pyramid and an inverse discrete wavelet transform (IDWT) is applied to obtain a restored gridline free image as shown in Fig. 2.
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 fg, the pixel resolution Rs and the gridline artifacts direction in image d.
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.
Fig. 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. Fig. 4(b) shows a selected region in Fig. 4(a). Fig. 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, fg was 5 lines/mm (50 lines/cm), and fs was 7.194 pixels/mm. Because fs < 2fg, 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): let k1= 1, we obtained festimate = 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 * festimate/fs =1249 on the frequency axis. Because the FFT spectrum is symmetrical, we only show the positive frequencies’ half axis in the following illustration. Fig. 4(f) shows the 1-dimensional Fourier spectrum of the sum of several rows (No. 90–100) of the original image. We notice that 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. Fig. 4(d) shows one of the sub-images that contain mainly the gridline signal that was detected by the module described in Section 3.2.2. When we implemented the Gaussian band-stop filter expressed in Section 3.3 on this sub-image, we obtained Fig. 4(e). The filtered sub-image was set into the pyramid and after applying an IDWT, we obtained the restored gridline free image, as shown in Fig. 4(c). The 1-dimensional 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 non-clinical images because Lin et al.18 reported that using geometric-shaped 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 gridline-free 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 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 sub-images 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 1. 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 sub-images without any analysis. The direct Gaussian band-stop filter method was the slowest because it applied the relatively time consuming filter process on the original high resolution image. Our proposed method achieved a relatively high computation speed because the time consuming Gaussian band-stop filter was applied to a sub-image with a reduced resolution compared with the original image one. Consequently it was almost three times faster than the direct Gaussian band-stop filter method but two times slower than Sasada’s zero-assign method.
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 sub-images. Several sub-images 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 sub-images to remove the gridline artifacts 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
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.