1.

INTRODUCTION

Energy discriminating photon counting detector (PCD) enables estimating line-integrals of basis materials. It provides various potential applications in PCD-based computed tomography (CT) compared to the conventional energy integrated detector-based CT.¹ The maximum likelihood (ML) principle has been a gold standard for estimating basis line-integrals in PCD-CT.^2,3 It has the optimal statistical property-unbiased and achieves the Cramer-Rao lower bound (CRLB), the minimum noise for an unbiased estimator. However, the noise level in the estimated basis sinograms and the associated reconstructed basis images may not be satisfactory for clinical purposes.^4,5 When the dose level is low, or the attenuation is large, the estimates greatly suffer from noise.^4,6 One may apply low-pass filtering (LPF) in the estimated basis sinograms or statistical iterative reconstruction (SIR) by regularizing neighboring pixels in the image domain while exploiting noise statistics in the forward model.^4,5,7 Both approaches effectively reduce the noise in the basis images, but the degraded spatial resolution is a concern.

In this study, we propose a likelihood-based bilateral filter (LBF) to reduce noise while preserving spatial resolution. Bilateral filter (BF) is an effective edge-preserving method for image denoising.^8,9 In a region of interest (ROI; the central pixel is the one we want to correct), it gives more weight for a neighborhood pixel whose value is more similar to that of the centered one, in addition to a pixel distance weight (e.g., Gaussian shape). The method is effective when the central pixel in ROI is a reliable reference for filtering purposes. However, the estimated basis sinograms may be too noisy to apply in practice (we validate it in the simulation study later). Instead of the pixel value-dependent weighting, we apply the neighboring likelihood-dependent one. The proposed method first estimates basis sinograms using the ML principle. It then calculates log-likelihoods using the estimates in the neighboring pixels for given PCD measurements at the centered one. It gives more weight to the neighbors whose likelihoods are greater than others. It is possible that pixel values that are very different from the center one (the ML estimates) may produce likelihoods close to the ML, so the method reflects the statistical distance rather than the distance between noisy estimates.

We validate the proposed method for a two-material decomposition (water and bone) and compare it to conventional methods. The ML-based estimation for basis sinograms is the starting point for all the methods we consider here. We used Gaussian and the original bilateral filters in the estimated basis sinograms as the conventional methods. We also applied model-based iterative reconstruction with edge-preserved penalty in the image domain. Both noise-reduced images (basis images and CT images) and noise-FWHM (full width at half maximum) behavior demonstrate that the proposed method effectively reduces noise without sacrificing spatial resolution compared to the conventional ones.

2.

METHODS

2.1

The ML-based Estimation for Basis Sinograms in PCD-CT

The log-likelihood of unknown basis line-integrals for the j-th sinogram point (x_j = [x_1,j, x₂,j, ⋯, x_M,j]^T ∈ ℝ^M×1 [cm]) for a given PCD measurements y_j = y₁, j₁,y_2,j, ⋯, y_B,j]^T ∈ ℝ^B×1 [counts] is given by (without constant terms),

where y_b,j is the number of photons measured at the b-th energy bin of the PCD pixel, B is the total number of energy bins, and M is the number of basis materials. We assume that the counts in energy bins are independent Poisson random variables with the expected value (λ_b(x_j)) as follows:

for b = 1, 2, ⋯, B, where S_b(E) [counts/keV] is the bin-wise energy-dependent function, including incident spectrum, and detector response function² and ϕ(E) [cm^-1] is the energy-dependent basis function. The ML-based estimation can be described as follows:

for m = 1, 2, ⋯, M. It includes a box constraint to the line-integrals to restrain excessively large estimates. After the estimation for all detector channels and projection views, one can reconstruct the basis images by applying e.g., filtered backprojection (FBP).

2.2

Proposed Likelihood-based Bilateral Filtering (LBF)

The conventional bilateral filtration (BF) we applied to the basis sinograms is given by,

where N_j is the set of neighboring sinogram pixels centered by the j-th pixel and α(j) is a normalization factor. The first weighting function G(i,j,σ₁) is a spatial low-pass filter (LPF), and we used a Gaussian filter with a standard deviation σ₁. The second weighting function is given by , which has a control parameter σ₂, where || · || denotes the standard ℓ₂ norm (Euclidean distance).

We propose the likelihood-based bilateral filtration (LBF) as follows:

where the second weighting function is changed to and again β(j) is a normalization factor.

2.3

Simulation Settings and Conventional Methods

We used an incident x-ray spectrum (140 kVp) with a 5 mm aluminium pre-filtration.¹⁰ A numerical thorax phantom was used, as shown in Fig. 1 (512 × 512 with a pixel size 0.75 × 0.75mm).¹¹ We used water (1 g/cm³) and bone (0.925g/cm³) as the two basis functions for a two-material decomposition. We applied a fan-beam geometry (equidistance) and used the ASTRA toolbox for the projection and backprojection.¹² The source to center and center to detector were both 50 cm. The number of views per rotation was 1000. The number of PCD pixels was 998 with a detector pitch of 0.1 cm. We used four energy bins 20, 67, 77, 95keV with minimum CRLB when the basis line-integrals of water and bone were 30 and 6 cm, respectively.¹³ The total number of incident photons was 5 × 10⁵ for each ray. The minimum photon numbers measured at PCD with the phantom were around 350 (summed for all energy bins).

Figure 1.

Reconstructed basis and CT images using various methods. For each image, the left half is the bias calculated from 100 noise trials, and the right half is the result of a single noise trial.

As illustrated in Fig. 1 (the first row), we measured FWHM for the reconstructed water and CT images (indicated by arrows) and also standard deviation indicated by the rectangles (ROIs 1 and 2). The ROI2 was selected to examine the severe streak artifacts in the image. The standard deviation was calculated for a single noise trial for each method and the FWHM using an average image from 100 noise trials. We also showed the reconstructed two basis images and the synthesized CT image for a single noise trial and corresponding bias images. The ground truth was set as the FBP (Shepp-Logan) images following the ML estimates using a substantially high number of incident photons (10⁹).

We compared the proposed LBF to the ML estimates, the original BF, and the LPF given by,

For the ML estimates, we used the derivative-free Nelda-Mead method with the box constraint for water and bone, −2 ≤ x_water ≤ 40 cm and −5 ≤ x_bone ≤ 15 cm, respectively. We set σ₁ = 0.9 [pix] with a 3 × 3 window for BF and LBF, and we varied σ₂. For LPF, we also used a 3 × 3 window but with different values of σ₁. Once the filtration was applied to the basis sinograms, we also combined them to synthesize a CT image (at 65keV).

We also applied statistical iterative reconstruction (SIR) to the ML estimates as described in.⁷ We used edge preserved pseudo-Huber penalty¹⁴ with the separable surrogate-based monotonic algorithm.¹⁵ The detailed algorithm is given in⁵ (Appendix C). We set the maximum iterations as 2000 and terminated when the relative change in the cost function is less than 2 × 10^-4. There are two control parameters in SIR: δ for the shape of the pseudo-Huber penalty and λ for the regularization parameter (weight on the penalty term). We tried various combinations of δ and λ and selected the best one by observing the reconstructed images and the FWHM and noise measures.

3.

RESULTS

We present the reconstructed basis images and the CT image overlapped by the bias image in Fig. 1 for various methods. For basis images (water and bone), we selected the results showing similar noise levels in ROI1 for all the methods (except the ML estimates; see the arrows in Fig. 2(a)). For the CT image, we also selected the ones with similar noise levels in ROI2 (see the arrows in Fig. 2(b)). Thus, the CT image in Fig. 1 may not represent the linear combination of the water and bone images as shown in the same figure. The images from the ML estimates are noisier than other methods, but there is no bias at edges inside the phantom. The LPF is effective in reducing the streak artifacts in the basis images, but it shows noticeable bias at edges for all basis and CT images. The images using the BF also show substantial bias at the edges, and the streak artifact is larger than the LPF. However, the bias at the edges from the proposed LBF is smaller than those of the LPF and BF, and the reduced streak artifact is qualitatively similar to the LPF. It also clearly demonstrates the effectiveness of the proposed LBF over the original BF (statistical difference versus pixel value difference). The SIR shows substantial bias at the edges in the CT image compared to the basis images. Note that all the CT images have similar noise levels at ROI2, and it clearly shows the trade-off between spatial resolution and noise, compared to the results in LBF.

Figure 2.

Noise-FWHM plots calculated from the reconstructed water image (ROI1 and FWHM1) and CT image (ROI2 and FWHM2) for various methods. The arrows indicate the points we selected to show in Fig. 1.

We show the noise-FWHM plots measured as indicated in Fig. 1 for the water and CT images in Fig. 2(a) and (b), respectively. For the water image (FWHM1 and ROI1), the FWHM from the ML is around 0.8 mm, but it has the largest noise among the methods. The LPF can substantially reduce the noise by sacrificing spatial resolution. The tendency is similar for the BF, but it degrades the resolution more than the LPF for reducing the same amount of noise. Interestingly, the proposed LBF reduces noise without sacrificing the resolution (even slight improvement). However, for a large value of σ₂, we could see slight but visible artifacts between bones in the average image (from 100 noise trials). It corresponds to the last point in the curve, and we indicated it with a cross mark (dotted line). The SIR presents better property on the noise-FWHM than the LPF, but the proposed method shows a slightly smaller resolution than the SIR at some point (indicated by arrow). The overall tendency in the CT image (FWHM2 and ROI2-streak artifacts; Fig. 2(b)) is similar to the water image. However, the noise-FWHM relation of SIR is not as superior as shown in Fig. 1(a). We searched for all the possible combinations of regularization parameters used for water and bone basis images, but we might have missed the (unknown) optimal parameters.