Effect of Filtered Back-Projection Filters to Low-Contrast Object Imaging in Ultra-High-Resolution (UHR) Cone-Beam Computed Tomography (CBCT)

In this study, the effect of filter schemes on several low-contrast materials was compared using standard and ultra-high-resolution (UHR) cone-beam computed tomography (CBCT) imaging. The performance of the UHR-CBCT was quantified by measuring the modulation transfer function (MTF) and the noise power spectrum (NPS). The MTF was measured at the radial location around the cylindrical phantom, whereas the NPS was measured in the eight different homogeneous regions of interest. Six different filter schemes were designed and implemented in the CT sinogram from each imaging configuration. The experimental results indicated that the filter with smaller smoothing window preserved the MTF up to the highest spatial frequency, but larger NPS. In addition, the UHR imaging protocol provided 1.77 times better spatial resolution than the standard acquisition by comparing the specific spatial frequency (f50) under the same conditions. The f50s with the flat-top window in UHR mode was 1.86, 0.94, 2.52, 2.05, and 1.86 lp/mm for Polyethylene (Material 1, M1), Polystyrene (M2), Nylon (M3), Acrylic (M4), and Polycarbonate (M5), respectively. The smoothing window in the UHR protocol showed a clearer performance in the MTF according to the low-contrast objects, showing agreement with the relative contrast of materials in order of M3, M4, M1, M5, and M2. In conclusion, although the UHR-CBCT showed the disadvantages of acquisition time and radiation dose, it could provide greater spatial resolution with smaller noise property compared to standard imaging; moreover, the optimal window function should be considered in advance for the best UHR performance.


Ultra-High-Resolution Cone-Beam Computed Tomography System and Imaging Configurations
A photograph and a specification of a prototype CBCT system are provided in Figure 1 and Table 1. Our system was mounted with an amorphous silicon (aSi)-based thin-film transistor (TFT) array FPD (PaxScan 4030CB, Varian Imaging Products, Palo Alto, CA) and was operated in full and binning acquisition modes. As shown in Table 2, the imaging configuration was categorized into four Sensors 2020, 20, 6416 3 of 14 subsections according to the two acquisition resolution setups and two reconstructed image resolutions. Each configuration was named depending on the row and column number of the matrix.
Sensors 2020, 20, x 3 of 14 mm slice thicknesses were chosen based on previous studies [2][3][4]. The readout time of FPD with a 2 × 2 binning mode acquisition was four times faster than that of a full mode acquisition; therefore, a lower total acquisition time and lower radiation exposure were achievable owing to the higher framerate in the binning mode. All FBP reconstruction algorithms were self-programmed and coded in C++ with the CUDA toolkit version 10.0 using a single GPU card (GTX Titan-Xp, NVIDIA Co., Ltd., Santa Clara, CA, USA).

Figure 1.
Photograph of the prototype cone-beam computed tomography (CBCT) system capable of both standard and ultra-high resolution (UHR) acquisition.   The center of rotation of the system was registered using the calibration phantom while rotating a full 360 • with a 1 • angle step for projection view image acquisition of 361 images. The 0.25 and Sensors 2020, 20, 6416 4 of 14 0.5 mm slice thicknesses were chosen based on previous studies [2][3][4]. The readout time of FPD with a 2 × 2 binning mode acquisition was four times faster than that of a full mode acquisition; therefore, a lower total acquisition time and lower radiation exposure were achievable owing to the higher framerate in the binning mode. All FBP reconstruction algorithms were self-programmed and coded in C++ with the CUDA toolkit version 10.0 using a single GPU card (GTX Titan-Xp, NVIDIA Co., Ltd., Santa Clara, CA, USA).

CT Performance Phantom
We used the CIRS Model 610 American Association of Physicists in Medicine (AAPM) CT performance phantom to measure spatial resolution and noise property. The CT number linearity insert (Part No. 610-02), which includes five cylinders with different densities, was a targeted imaging object for resolution measurement. The detailed specifications of the inserted cylinders are given in Table 3. Each cylinder has the same size and shape and has a low contrast against the background material, thus presenting a small absolute signal difference between the two materials. The larger the material index, the smaller the difference between the background and target material densities. Note that a small absolute difference between the densities of two materials does not always guarantee a small image contrast because the CT numbers are represented by the linear attenuation coefficients which are dependent on both X-ray energy and density. The insert (Part No. 610-01-05) is comprised of a uniform material with an aluminum pin at the center, and is a good candidate for measuring noise power. We assumed that the noise behaviors were the same for all materials because quantum and electronic noise, which are both stochastic events, are dominant over the entire area.

Ramp Filter Design in Spatial Domain and Six Different Window Functions
Linear filtering can be categorized into two methods: applying the convolution kernel in the spatial domain and linear multiplication of a transfer function in the Fourier domain. A band-limited ramp filter constructed in the Fourier domain is defined as follows: where ω is the discretized spatial frequency by considering the Nyquist frequency. However, the ramp filter in Equation (1) has a zero at ω = 0 lp/mm such that the signals at the DC offset (zero frequency Sensors 2020, 20, 6416 5 of 14 component) after linear multiplication go to zero. The Fourier transform of the ramp convolution kernel constructed in the spatial domain can be defined as follows [17]: if n = 0 if n is an even number if n is an odd number The ramp filter in Equation (2) does not include a zero, as shown in the comparison of the two shapes in Figure 2. Filtering with non-zero conditions avoids the zero signals that might have occurred if the filters were used with zero conditions. The ramp filter in Equation (2) does not include a zero, as shown in the comparison of the two shapes in Figure 2. Filtering with non-zero conditions avoids the zero signals that might have occurred if the filters were used with zero conditions. Many window functions have been introduced depending on the strength of noise suppression at different cutoff frequencies for each purpose [18]. However, the reduction of the critical signal is inevitable during noise suppression; therefore, the optimal window function is often heuristically chosen after multiple reconstruction trials. Six different smoothing windows were implemented herein in the ramp filter. Each window function was followed by the equation summarized in Table 4, where a term L in (b), (c), (d), and (f) indicates the length of the window.  Many window functions have been introduced depending on the strength of noise suppression at different cutoff frequencies for each purpose [18]. However, the reduction of the critical signal is inevitable during noise suppression; therefore, the optimal window function is often heuristically chosen after multiple reconstruction trials. Six different smoothing windows were implemented herein in the ramp filter. Each window function was followed by the equation summarized in Table 4, where a term L in (b), (c), (d), and (f) indicates the length of the window.

Modulation Transfer Function (MTF)
Spatial resolution for each imaging configuration and each filter scheme was evaluated by the MTF measurement of the cylindrical materials as conducted by Richard et al. [19]. After subtracting the two-dimensional planar fit from the original region of interest (ROI) of each targeted cylinder, the radial pixel values around the edge of the circular shape were rearranged to yield a one-dimensional edge spread function (ESF). When converting the image grid from a Cartesian to polar map, the center of each disk was measured on a binary image through a gray-level threshold. The ESF, which is equivalent to the radial profile of the circle, was resampled with one-tenth of the reconstructed pixel size to reduce the non-uniformly distributed pixel noise [20]. The final ESF was derived by averaging the ESFs measured from consecutive axial slices. The MTF was the Fourier amplitude of the derivative of the ensemble-averaged ESF. In addition, the high-frequency noise of the ESF derivative was relieved through a Hanning window having the same length as the ESF size. The overall process of radial MTF measurement is depicted in Figure 3. Table 4. Description of each window that was implemented with the ramp filter.

Modulation Transfer Function (MTF)
Spatial resolution for each imaging configuration and each filter scheme was evaluated by the MTF measurement of the cylindrical materials as conducted by Richard et al. [19]. After subtracting the two-dimensional planar fit from the original region of interest (ROI) of each targeted cylinder, the radial pixel values around the edge of the circular shape were rearranged to yield a onedimensional edge spread function (ESF). When converting the image grid from a Cartesian to polar map, the center of each disk was measured on a binary image through a gray-level threshold. The ESF, which is equivalent to the radial profile of the circle, was resampled with one-tenth of the reconstructed pixel size to reduce the non-uniformly distributed pixel noise [20]. The final ESF was derived by averaging the ESFs measured from consecutive axial slices. The MTF was the Fourier amplitude of the derivative of the ensemble-averaged ESF. In addition, the high-frequency noise of the ESF derivative was relieved through a Hanning window having the same length as the ESF size. The overall process of radial MTF measurement is depicted in Figure 3.

Normalized Noise Power Spectrum
The normalized noise power spectrum (NNPS) was measured to quantify the noise level in the homogeneous volume of interest (VOI) of the poly methyl methacrylate (PMMA) background. The three-dimensional (3D) NPS was measured as described in Figure 4. The eight different VOIs without interference of any structure with the size of 150 × 150 × 45 (300 × 300 × 90 for high-resolution reconstruction) were selected for measuring the 3D NPS. Each sub-volume overlapped with others to evaluate the radially and symmetrically distributed noise property (location independent noise pattern) [21].
Each mean subtracted sub-volume patch was Fourier transformed, absolute squared, and ensemble averaged to yield the power spectrum as follows [22]: where , , and are spatial frequencies (mm −1 ), , , and are pixel sizes (mm), , , and are the numbers of voxels in the sub-volume patch, ℱ[•] is the fast Fourier transform operator, and ( , , ) and ̅ indicate each voxel value and the mean intensity of the sub-volume patch, respectively. The 1D NNPS can be derived by radially averaging the 3D NPS [23].

Normalized Noise Power Spectrum
The normalized noise power spectrum (NNPS) was measured to quantify the noise level in the homogeneous volume of interest (VOI) of the poly methyl methacrylate (PMMA) background. The three-dimensional (3D) NPS was measured as described in Figure 4. The eight different VOIs without interference of any structure with the size of 150 × 150 × 45 (300 × 300 × 90 for high-resolution reconstruction) were selected for measuring the 3D NPS. Each sub-volume overlapped with others to evaluate the radially and symmetrically distributed noise property (location independent noise pattern) [21].
Each mean subtracted sub-volume patch was Fourier transformed, absolute squared, and ensemble averaged to yield the power spectrum as follows [22]: where f x , f y , and f z are spatial frequencies (mm −1 ), d x , d y , and d z are pixel sizes (mm), N x , N y , and N z are the numbers of voxels in the sub-volume patch, F [·] is the fast Fourier transform operator, and S(i, j, k) and S indicate each voxel value and the mean intensity of the sub-volume patch, respectively. The 1D NNPS can be derived by radially averaging the 3D NPS [23].

Filter Shape
Six different Fourier transformed and band-limited filters designed in the spatial domain with regard to the frequency response are depicted in Figure 5. Because the Fourier transformed sinograms were forced to be band limited with a band width of 0.5, the signals outside of the band frequency range went to zero, as shown in Figure 5. Similarly, each window function was also band limited and multiplied by the band-limited ramp filter. The magnitude of the filter at high frequencies was rejected when going from scheme (a) to (f) in Table 4, which is generally interpreted as noise suppression. Unlike other filters, some of the value of the flat-top window are negative.  Figure 6 shows the reconstructed images with configuration (1, 1) using the Hanning window and its cropped ROI images around the centers of five different materials. The relative contrast between each material and background with standard deviation error are plotted in Figure 6f. All five materials showed a low contrast, showing a small relative contrast below 0.15 (maximum contrast is 1). As mentioned above, the higher density material does not always represent the higher

Filter Shape
Six different Fourier transformed and band-limited filters designed in the spatial domain with regard to the frequency response are depicted in Figure 5. Because the Fourier transformed sinograms were forced to be band limited with a band width of 0.5, the signals outside of the band frequency range went to zero, as shown in Figure 5. Similarly, each window function was also band limited and multiplied by the band-limited ramp filter. The magnitude of the filter at high frequencies was rejected when going from scheme (a) to (f) in Table 4, which is generally interpreted as noise suppression. Unlike other filters, some of the value of the flat-top window are negative.

Filter Shape
Six different Fourier transformed and band-limited filters designed in the spatial domain with regard to the frequency response are depicted in Figure 5. Because the Fourier transformed sinograms were forced to be band limited with a band width of 0.5, the signals outside of the band frequency range went to zero, as shown in Figure 5. Similarly, each window function was also band limited and multiplied by the band-limited ramp filter. The magnitude of the filter at high frequencies was rejected when going from scheme (a) to (f) in Table 4, which is generally interpreted as noise suppression. Unlike other filters, some of the value of the flat-top window are negative.  Figure 6 shows the reconstructed images with configuration (1, 1) using the Hanning window and its cropped ROI images around the centers of five different materials. The relative contrast between each material and background with standard deviation error are plotted in Figure 6f. All five materials showed a low contrast, showing a small relative contrast below 0.15 (maximum contrast is 1). As mentioned above, the higher density material does not always represent the higher  Figure 6 shows the reconstructed images with configuration (1, 1) using the Hanning window and its cropped ROI images around the centers of five different materials. The relative contrast between each material and background with standard deviation error are plotted in Figure 6f. All five materials showed a low contrast, showing a small relative contrast below 0.15 (maximum contrast is 1).

Reconstructed Images with Different Filters and Configurations
Sensors 2020, 20, 6416 8 of 14 As mentioned above, the higher density material does not always represent the higher CT number when we measure the contrast between each material and the background (PMMA). M2 and M5 showed the lowest contrast among the five materials.
Sensors 2020, 20, x 8 of 14 CT number when we measure the contrast between each material and the background (PMMA). M2 and M5 showed the lowest contrast among the five materials. The reconstructed cropped images of M1 with different configurations using the Hanning window are shown in Figure 7a. Figure 7b shows the radial profiles of each image grid in Figure 7a. The images reconstructed using standard detector resolution (configuration (1, 1) and (1, 2)) showed an unstable fluctuation in their radial profiles at the initial radial location. On the contrary, the images of configurations (2, 1) and (2, 2) showed relatively flat signals. To understand the effect of filter schemes on image quality, the reconstructed images of M1 with different filters using configurations (1, 1) and (2, 2) are shown in Figure 8a. The radial profiles in Figure 8b correspond to the bottom row images in Figure 8a) (configuration (2, 2)). The fluctuations of the radial profiles are gradually smoothed with an increase in the index number of filter schemes, demonstrating that the high-frequency noise was rejected by using the smoothing windows. The more oscillations in the signal, the coarser the MTF curve, as shown in Figure 9a. The reconstructed cropped images of M1 with different configurations using the Hanning window are shown in Figure 7a. Figure 7b shows the radial profiles of each image grid in Figure 7a. The images reconstructed using standard detector resolution (configuration (1, 1) and (1, 2)) showed an unstable fluctuation in their radial profiles at the initial radial location. On the contrary, the images of configurations (2, 1) and (2,2) showed relatively flat signals.  The reconstructed cropped images of M1 with different configurations using the Hanning window are shown in Figure 7a. Figure 7b shows the radial profiles of each image grid in Figure 7a. The images reconstructed using standard detector resolution (configuration (1, 1) and (1, 2)) showed an unstable fluctuation in their radial profiles at the initial radial location. On the contrary, the images of configurations (2, 1) and (2,2) showed relatively flat signals. To understand the effect of filter schemes on image quality, the reconstructed images of M1 with different filters using configurations (1, 1) and (2,2) are shown in Figure 8a. The radial profiles in Figure 8b correspond to the bottom row images in Figure 8a) (configuration (2, 2)). The fluctuations of the radial profiles are gradually smoothed with an increase in the index number of filter schemes, demonstrating that the high-frequency noise was rejected by using the smoothing windows. The more oscillations in the signal, the coarser the MTF curve, as shown in Figure 9a. To understand the effect of filter schemes on image quality, the reconstructed images of M1 with different filters using configurations (1, 1) and (2,2) are shown in Figure 8a. The radial profiles in Figure 8b correspond to the bottom row images in Figure 8a (configuration (2, 2)). The fluctuations of the radial profiles are gradually smoothed with an increase in the index number of filter schemes, demonstrating that the high-frequency noise was rejected by using the smoothing windows. The more oscillations in the signal, the coarser the MTF curve, as shown in Figure 9a.

Modulation Transfer Function
Six different MTFs for each filter scheme measured in the reconstructed images of M1 are shown in Figure 9. The higher the resolution of the reconstructed images, the better the MTF is preserved up to the Sensors 2020, 20, 6416 9 of 14 high frequencies. In Figure 9f, f 50 , which indicates the specific spatial frequency when the MTF is dropped to 0.5, was 1.39, 1.40, 2.52, and 2.57 lp/mm for the configurations (1, 1), (1,2), (2,1), and (2, 2), respectively. The effect of detector resolution on the reconstruction image resolution was minor when we compared the curves between configurations (1, 1) and (1, 2) (or (2, 1) and (2,2)).     The MTF curves measured in the reconstructed images of each material using configurations (1, 1) and (2,2) are shown in Figures 10 and 11. As shown in Figure 11f, the f 50 s were 0.94, 1.86, 2.05, and 2.52 and 1.86 lp/mm from M1 to M5, respectively, which demonstrates that MTFs were preserved up to high frequencies of the order of M3, M4, M1, M5, and M2; that is, in the order of the relative contrast in Figure 6f. In contrast, the imaging configuration (1, 1) not only did not follow the order of contrast, but also presented different orders of f 50 s for the different filter schemes.

Modulation Transfer Function
Six different MTFs for each filter scheme measured in the reconstructed images of M1 are shown in Figure 9. The higher the resolution of the reconstructed images, the better the MTF is preserved up to the high frequencies. In Figure 9f, f50, which indicates the specific spatial frequency when the MTF is dropped to 0.5, was 1.39, 1.40, 2.52, and 2.57 lp/mm for the configurations (1, 1), (1, 2), (2, 1), and (2, 2), respectively. The effect of detector resolution on the reconstruction image resolution was minor when we compared the curves between configurations (1, 1) and (1, 2) (or (2,1) and (2,2)).
The MTF curves measured in the reconstructed images of each material using configurations (1,1) and (2,2) are shown in Figures 10 and 11. As shown in Figure 11f, the f50s were 0.94, 1.86, 2.05, and 2.52 and 1.86 lp/mm from M1 to M5, respectively, which demonstrates that MTFs were preserved up to high frequencies of the order of M3, M4, M1, M5, and M2; that is, in the order of the relative contrast in Figure 6f. In contrast, the imaging configuration (1, 1) not only did not follow the order of contrast, but also presented different orders of f50s for the different filter schemes.   Figure 12 shows the radially averaged 1D NPS for each configuration with different filter  Figure 12 shows the radially averaged 1D NPS for each configuration with different filter schemes. The standard reconstructed image resolution (configuration (1, 1) and (2, 1)) gave higher noise properties compared to the high-resolution images (configuration (1,2) and (2,2)). We also observed that the peak of the 1D NNPSs from the higher detector resolution was at larger spatial frequencies, which demonstrates that the noise was distributed up to a higher frequency when the smaller pixels were used in the detector. The NPSs decreased as the intensity of high-frequency smoothing increased.  Figure 12 shows the radially averaged 1D NPS for each configuration with different filter schemes. The standard reconstructed image resolution (configuration (1, 1) and (2, 1)) gave higher noise properties compared to the high-resolution images (configuration (1,2) and (2,2)). We also observed that the peak of the 1D NNPSs from the higher detector resolution was at larger spatial frequencies, which demonstrates that the noise was distributed up to a higher frequency when the smaller pixels were used in the detector. The NPSs decreased as the intensity of high-frequency smoothing increased.

Discussion
We herein designed band-limited filters for all schemes. These can effectively retrieve the sampled projections because the projections are discretized into each detector pixel so that it is band limited in the Fourier domain [24]. As a result, band-limited filters lead to the removal of unnecessary noise signals at high frequencies.
There is no universal filter in CT imaging; therefore, the user should select an optimal smoothing window to observe the detailed internal structure with a purpose. Selecting an optimal window function is often based on experience rather than theory because we do not have a high level of knowledge about whether the imaging object is lying under a low-, mid-, or high-frequency range [25]. Thus, comparing the initial imaging performance of different filters and choosing the best solution for one's purpose is a good approach [25]. The most important factor when selecting the filter scheme is the manner in which the filter removes as many of the unnecessary components as possible in the frequency domain. In this experimental study, the signals near the edge of each material that we aimed to observe mostly lie in the low-frequency range, and show severe MTF distortion in the images applied with a high-pass filter, such as Butterworth A in Figure 9a. In contrast, the results in Figure 11f indicate that the flat-top window preserved the MTF up to a high frequency without an aliasing among the six filter schemes in our experiment. This is because the reconstructed images applied with the flat-top window not only resulted in uniform pixel values but also showed small oscillations (less noise) in both the target and background, as shown in the radial profiles in Figure 8b.
The flat-top window is used for cases in which a frequency component is required to be measured with great accuracy, e.g., a fixed-sine source [26]. Measuring the MTFs in the frequency domain could be interpreted as a discrimination of the signals spreading near the circular edge region. If a much larger signal difference exists between the target and the background, such as the tungsten edge, filter selection would not have been significant. However, we measured the MTFs for materials having no significant signal difference against the background material (low-contrast imaging); therefore, the amplitude accuracy was a key factor because the principal components in the Fourier domain were largely positioned in the low-frequency area [27].
The MTFs were preserved well at higher frequencies from the images reconstructed with a higher resolution. We observed that there was an MTF preservation loss up to 1.77 times by comparing the f 50 between configurations (1, 1) and (2,2) in Figures 10f and 11f when using the same target material and detector resolution. Therefore, using a UHR imaging protocol rather than a standard imaging configuration is recommended to understand the fine sharpness of low-contrast material if the detector is available to be operated at a higher resolution. However, the high-level smoothing window is not recommended for standard resolution imaging configuration, as shown by the disagreement in the order of relative contrast in Figure 10. As shown in Figure 10, the flat-top window provided little difference in f 50 s for different materials even though there was a clear discrimination in UHR imaging protocol. This was because the flat-top window overly smoothed the low-contrast object in the standard imaging, whereas the smoothing was still effective in UHR mode.
The trend of 1D NNPS in the configuration (2, 1) showed that the noise was distributed over all of the spatial frequencies. This demonstrates the back-projection from the high-resolution to small-image array would largely reduce the quantum noise and result in uniformly distributed noise.
The main drawback of this study is that all materials used to measure the MTFs had low contrast against the background PMMA intensity. This limits the study of higher-object-contrast materials such as bone and contrast-enhanced imaging. Our future study will be directed toward the effect of various filter setups on higher-object-contrast materials.

Conclusions
In summary, we observed the effect of filter schemes on several low-contrast materials using standard and UHR imaging protocols. Although UHR image acquisition requires a higher acquisition time and greater radiation exposure, we obtained spatial resolution up to 1.77 times higher than that of standard acquisition. In addition, the performance of UHR was affected by the FBP filter schemes, showing different f 50 values and different noise patterns for different filters. Therefore, one should consider the optimal window function that can provide the best performance when observing the fine structure of the imaging object before UHR acquisition while comparing both the MTF and NPS.