Next Article in Journal
A Method for Accessing the Non-Slip Function of Socks in an Acute Maneuver
Previous Article in Journal
A Wi-Fi-Based Passive Indoor Positioning System via Entropy-Enhanced Deployment of Wi-Fi Sniffers
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Joint-Parameter Estimation and Bayesian Reconstruction Approach to Low-Dose CT †

1
Department of Radiology, Stony Brook University, Stony Brook, NY 11794, USA
2
Department of Medical Physics, Memorial Sloan Kettering Cancer Center, New York, NY 10065, USA
3
Department of Preventive Medicine, Stony Brook University, Stony Brook, NY 11794, USA
4
Department of Engineering Science and Physics, CUNY/CSI, Staten Island, NY 10314, USA
*
Author to whom correspondence should be addressed.
This paper is an extended version of our conference paper published in Gao, Y.; Liang, Z.; Lu, S.; Shi, Y.; Chang, S.; Hou, W. Hyperparameter Selection for Bayesian Image Reconstruction by Mimicking Physical Crystallization. In Proceedings of the IEEE Symposium on Nuclear Science (NSS/MIC), Boston, MA, USA, 31 October 2020–7 November 2020. (https://ieeexplore.ieee.org/document/9508084, accessed on 28 December 2022).
Sensors 2023, 23(3), 1374; https://doi.org/10.3390/s23031374
Submission received: 29 December 2022 / Revised: 18 January 2023 / Accepted: 23 January 2023 / Published: 26 January 2023
(This article belongs to the Section Biosensors)

Abstract

:
Most penalized maximum likelihood methods for tomographic image reconstruction based on Bayes’ law include a freely adjustable hyperparameter to balance the data fidelity term and the prior/penalty term for a specific noise–resolution tradeoff. The hyperparameter is determined empirically via a trial-and-error fashion in many applications, which then selects the optimal result from multiple iterative reconstructions. These penalized methods are not only time-consuming by their iterative nature, but also require manual adjustment. This study aims to investigate a theory-based strategy for Bayesian image reconstruction without a freely adjustable hyperparameter, to substantially save time and computational resources. The Bayesian image reconstruction problem is formulated by two probability density functions (PDFs), one for the data fidelity term and the other for the prior term. When formulating these PDFs, we introduce two parameters. While these two parameters ensure the PDFs completely describe the data and prior terms, they cannot be determined by the acquired data; thus, they are called complete but unobservable parameters. Estimating these two parameters becomes possible under the conditional expectation and maximization for the image reconstruction, given the acquired data and the PDFs. This leads to an iterative algorithm, which jointly estimates the two parameters and computes the to-be reconstructed image by maximizing a posteriori probability, denoted as joint-parameter-Bayes. In addition to the theoretical formulation, comprehensive simulation experiments are performed to analyze the stopping criterion of the iterative joint-parameter-Bayes method. Finally, given the data, an optimal reconstruction is obtained without any freely adjustable hyperparameter by satisfying the PDF condition for both the data likelihood and the prior probability, and by satisfying the stopping criterion. Moreover, the stability of joint-parameter-Bayes is investigated through factors such as initialization, the PDF specification, and renormalization in an iterative manner. Both phantom simulation and clinical patient data results show that joint-parameter-Bayes can provide comparable reconstructed image quality compared to the conventional methods, but with much less reconstruction time. To see the response of the algorithm to different types of noise, three common noise models are introduced to the simulation data, including white Gaussian noise to post-log sinogram data, Poisson-like signal-dependent noise to post-log sinogram data and Poisson noise to the pre-log transmission data. The experimental outcomes of the white Gaussian noise reveal that the two parameters estimated by the joint-parameter-Bayes method agree well with simulations. It is observed that the parameter introduced to satisfy the prior’s PDF is more sensitive to stopping the iteration process for all three noise models. A stability investigation showed that the initial image by filtered back projection is very robust. Clinical patient data demonstrated the effectiveness of the proposed joint-parameter-Bayes and stopping criterion.

1. Introduction

Tomography, as one non-invasive imaging technique, has a wide application in many scientific fields, e.g., physics, chemistry, astronomy, medicine, etc. [1]. There are a variety of tomographic modalities, such as computed tomography (CT), nuclear medicine, ultrasound diffraction tomography and so on. Under certain conditions, the measurements in each modality can be converted into samples of the Radon transform distribution from which we wish to reconstruct the corresponding image [2]. The inversion of Radon transform is well-established by deterministic methods: for example, the filtered back-projection (FBP) [3]. However, there are many scenarios in which the measurements are incomplete data, which will significantly degrade the performance of the deterministic method. For instance, noise artifacts will appear in low-dose CT imaging due to the incomplete data. In such a scenario, the Bayesian theorem-based penalized maximum likelihood methods [4], which incorporate the data fidelity and prior knowledge into one unified objective function, provide an alternative way to deal with the imperfect data acquisitions.
Most penalized maximum likelihood reconstruction methods include one freely adjustable hyperparameter to balance the strength between data fidelity and prior terms. Specifically, the hyperparameter controls the noise–resolution tradeoff in the reconstructed image. When the hyperparameter is too small, the resulting reconstruction will be noisy; when the hyperparameter is too large, the resulting reconstruction will be oversmoothed. In practice, this hyperparameter is determined empirically in a trial-and-error fashion for an optimal result, which is very time-consuming. In the past decades, several methods have been proposed to determine the optimal hyperparameter, such as the discrepancy principle [5], generalized cross validation [6], the L-curve method [7], the balance principle [8] and the Z-index method [9]. These methods directly or indirectly quantify the relationship between the image quality and the hyperparameter value. For example, the idea of the discrepancy principle is to select the hyperparameter so that the fidelity residual norm is close to a prior upper bound. Therefore, it requires the estimation of the noise level. The generalized cross-validation method is based on the hypothesis that a good hyperparameter should predict missing data values. It constructs an objective function of a hyperparameter and minimizes it. With the objective function, the L-curve plots the curve of fidelity residual norm vs. prior norm to determine the hyperparameter. Then, the optimal hyperparameter is at the corner of the “L” shape curve. The Z-index method plots the sparsity prior vs. the hyperparameter. The optimal parameter is selected at the second corner of the “Z” shape curve. Some iterative algorithms are also implemented to find the optimal hyperparameter [10,11]. However, all the aforementioned methods are conducted trial-and-error style and very time-consuming, even though the trials can sometimes be performed automatically by programs. To avoid this difficulty, a theory-based strategy is desired to determine the free hyperparameter, which can help to reduce the computational load.
Therefore, this study aims to investigate a theory-based strategy for parameter estimation in Bayesian image reconstruction without a freely adjustable hyperparameter. The main purpose of this study is the introduction of parameters to ensure that the condition of probability density function (PDF) in Bayes’ law is satisfied. Specifically, the Bayesian image reconstruction problem is formulated by two PDFs, one for the data fidelity term and the other for the prior term. Each PDF contains a parameter. While these two parameters ensure the PDFs completely describe the data and prior terms, they cannot be determined by the acquired data; thus, they are called complete but unobservable parameters. To differentiate from the conventional hyperparameter, we use the introduced parameters, which makes the PDF completely specified, and these are referred to as unobservable parameters hereafter. Estimating these two unobservable parameters becomes possible under the conditional expectation and maximization for the image reconstruction, given the acquired data and the PDFs. This leads to an iterative algorithm, which jointly estimates the unobservable parameters and computes the to-be reconstructed image by maximizing a posteriori probability (MAP), denoted as joint-parameter-Bayes. Although the joint Bayesian approach has been used in some studies [12,13,14], comprehensive studies are lacking that address the related issues, such as explanation, algorithm stability and so on. For instance, Hsiao et al. proposed one deterministic annealing method for transmission tomography based on Bayesian image reconstruction, where joint-MAP is used to estimate the parameters of mixed Gamma distribution in its prior model [12]. Cai et al. used joint-MAP to estimate the decomposition fractions and observation variance in Bayesian based dual-energy CT [13]. Zhang et al. applied this joint-MAP into cone beam X-ray luminescence computed tomography [14]. However, there are no extensive theory-based studies and discussions for parameter estimation in the mentioned studies. In this work, we investigate the stability of the proposed approach through factors including the initialization of the iterative process, the PDF specification, and renormalization in an iterative manner.
Additionally, by monitoring the unobservable parameters, another contribution of this work is to provide a sensitive metric to stop the iteration. Currently, most reconstruction terminates when the convergence criterion is met. It has already been reported that a better image can be obtained by stopping an iterative maximum likelihood algorithm early, well before convergence [15,16,17]. Veklerov and Llacer proposed a criterion based on a probabilistic interpretation [15,16]. Hebert reported that one explanation of the benefit of early termination is that the image degradation is one’s perception of the convergence process from the smooth initial image to the non-smooth maximum likelihood solution [17]. All these works inspired us to investigate the stopping criterion for this joint-parameter-Bayes method. The unobservable parameter of prior is found to have a turning point coincident with the optimal image quality in their curves along with iteration times. It is also found that the stopping criterion should be part of the joint-parameter-Bayes method. Ignoring the stopping criterion, however, may result in unsatisfying results.
Finally, given the acquired data, an optimized reconstructed image is obtained without any freely adjustable hyperparameters by satisfying the PDF condition for both data likelihood and a priori probability of Bayes’ law, and by satisfying the stopping criterion. Finally, although limited to the CT application in this paper, the proposed Bayesian image reconstruction can be generalized to other modalities and applications. This work is an extension of our previous conference paper of the 2020 IEEE Symposium on Nuclear Science [18].

2. Materials and Methods

2.1. Bayesian Image Reconstruction for CT

For an ill-posed problem such as low-dose CT (LdCT) image reconstruction, Bayesian reconstruction provides an alternative method for an improved solution compared to the deterministic method. According to Bayes’ law, we can have:
p ( μ | y ) = p ( y | μ ) p ( μ ) p ( y ) p ( y | μ ) p ( μ )
where the vector y = { y i | i = 1 , , I } I × 1 represents the measured LdCT post-log sinogram data (with I data elements) and vector μ = { μ j | j = 1 , , J } J × 1 represents the linear attenuation coefficient distribution (with J image elements). The acquired data can be written as:
y = y ¯ + ζ = A μ + ζ
where y ¯ I × 1 is the expectation of y , A = { a i j | i = 1 , , I ; j = 1 , , J } I × J is called the system matrix, and ζ I × 1 indicates the noise vector corresponding to the measured sinogram data.

2.1.1. Multivariate Gaussian Distribution Data Model

The post-log sinogram data of the i t h X-ray, y i , is the logarithm of transmission data N i and the incident photon flux N 0 , which is calculated as:
y i = l o g   { N i N 0 } .
In the transmission domain, the noise is usually modeled as the sum of Poisson and Gaussian distribution [19] as:
N i   ~   P o i s s o n   { N ¯ i } + G a u s s i a n   { m e ,   σ e 2 }
where N ¯ i is the mean transmission measurement of i t h ray, and m e and σ e 2 are the mean and variance value of the Gaussian type electronic noise. After logarithm, the non-linearity will be introduced, which can be approximated by Taylor expansion [20]. When only taking the leading order terms, the variance of post-log sinogram measurement is inversely proportional to the mean transmission data. Thus, the variance of i t h sinogram measurement can be expressed as:
V a r   ( y i )     1 / N ¯ i ,
where is the proportional operator. Under the condition that all measurements { y i } are independent to each other, the covariance matrix is diagonal. If we denote the measured transmission counts as ρ 2 = { ρ i 2 | i = 1 , , I } I × 1 , we have:
Σ y = d i a g { 1 / N ¯ i } = d i a g { 1 / ρ i 2 }
Given this knowledge, the PDF of the acquired sinogram can be expressed as:
p ( y | μ , s ) = i = 1 I 1 2 π s / ρ i 2 e x p { ( y i y ¯ i ) 2 2 s / ρ i 2 } .
where s is a parameter to completely define the variance.

2.1.2. Multivariate Gaussian MRF Prior Model

While the linear attenuation coefficients are not random variables, their distribution patterns across the image can be assumed as random variables and the patterns are frequently modeled as the Markov random field (MRF). In this paper, the well-established Gaussian MRF penalty is used [21]. We assume that the variance of each pixel in the image field is σ 2 , which is independent between two pixels. We can obtain the variance of the pair-wise pixel difference following the analysis for a pair-wise neighbor system as below:
V a r ( μ j μ k ) = V a r ( μ j ) + V a r ( μ k ) 2 C o v ( μ j , μ k ) = σ 2 + σ 2 0 = 2 σ 2
where V a r (.) represents the variance operator and C o v (.) denotes the covariance operator. If we denote t = σ 2 , the Gaussian MRF a priori PDF can be expressed as:
p ( μ | t ) = j = 1 J 1 2 π t   e x p { k Ω j ω j k ( μ j μ k ) 2 2 t } ,
where ω j k represents the normalized pair-wise MRF coefficients among the neighbors Ω j around the center j, k Ω j ω j k = 1 , and can be estimated based on some established knowledge or principles [22,23,24]. The scale factor is introduced as an unobservable parameter to ensure that the a priori PDF is completely specified for the pattern distribution to estimate the solution μ . One established MRF weight is used. In two-dimensional presentation, the eight image pixels weights in a 3 × 3 window can be given as:
[ 0.104 0.146 0.104 0.146 0 0.146 0.104 0.146 0.104 ]

2.2. Unobservable Parameter Estimation by Joint-Parameter-Bayes

Introducing the PDF of the two terms of Equations (7) and (9) into Equation (1), we can obtain the solution by MAP estimation, expressed as:
μ ^ M A P = argmax | μ   { log p ( μ | y ) }                                                                     = a r g   m a x | μ   { l o g   p ( ( y | μ ) ) + l o g   p ( μ ) }                                                                                         = a r g   m a x | μ     { l o g   p ( ( y | μ , s   ) ) + l o g   p ( μ | t ) }  
However, the s , t are two unobservable parameters to completely obtain the PDF of data and penalty distribution. Therefore, for a given image μ , we use the joint-parameter-Bayes estimation of the two unobservable parameters. We can estimate s , t by:
s ^ = a r g   m a x | s   { l o g   p ( ( y | μ , s ) )   } ,
t ^ = a r g   m a x | t     { l o g   p ( μ | t )   } .
According to the definition, for the nth iteration, with the current image estimate μ ( n ) , the data variance parameter s ^ ( n ) can be estimated by:
s ^ ( n ) = 1 I i = 1 I 1 N i ( y i [ A μ ( n ) ] i ) 2   .
The prior variance parameter t ^ ( n ) can be obtained by:
t ^ ( n ) = 1 J j = 1 J k Ω j ω j k ( μ j ( n ) μ k ( n ) ) 2 .
To solve the problem of Equation (11), we use the Newton method to iteratively find the solution. The pseudo code is summarized in Algorithm 1. For initialization, we first have the initial FBP image μ ^ from the measured post-log sinogram y . Then, a forward projection is applied on μ ^ with q = A   μ ^ . The residual r ^ = y     q is obtained to describe the difference between estimation and measurement. ρ 2 describes the covariance factor obtained from measurement. The two unobservable parameters are initialized by Equations (14) and (15). For each iteration, we first update the attenuation map given the value of s ^ and t ^ . Then, the residual r and the covariance factor ρ 2 are obtained based on the new attenuation map. Then, the s ^ and t ^ are updated orderly given the known attenuation. A j denotes the j t h column of the system matrix A .
Algorithm 1. Joint-MAP-Bayes
Initialization:
                    μ ^ = F B P { y } ;   q = A   μ ^   ;   r ^ = y     q
                    ρ 2 = d i a g { N 0 e q i }  
                    λ j = A j T ρ 2 A j ,   j
          Initialize s ^ and t ^ by Equations (14) and (15) with μ ^
For each iteration:
While (Stopping criterion is not met)
          For each voxel j:
                    μ ^ j o l d = μ ^ j
                    μ ^ j n e w = A j T ρ 2 r ^ + λ j μ ^ j o l d + s ^ t ^ k Ω j ω j k μ k λ j + s ^ t ^
                    μ ^ j = m a x { 0 ,   μ ^ j n e w }
                    r ^ = r ^ + A j ( μ ^ j o l d μ ^ j )
          end
                    ρ 2 = d i a g { N 0 e q i }
                    λ j = A j T ρ 2 A j
          Update s ^ and t ^ by Equations (14) and (15)

2.3. Stopping Criterion Investigation

According to the pseudo-code above, the stopping criterion is the only remaining part in our proposed framework. Conventionally, the mean squared difference between two successive reconstructed images is often employed as an indicator for the reconstruction convergency. It is also popular to assign a threshold of iteration times, e.g., 100 times based on experience. However, these methods always focus on the convergency from one indicator and pay little attention to image quality. As mentioned in the introduction, early stopping may provide better reconstruction results in the maximum likelihood algorithm. All these works inspired us to investigate the stopping criterion for this joint-parameter-Bayes method.
Based on the analysis from Section 2.2, we can see that there are two unobservable parameters in addition to the unknown attenuation map. Therefore, we can investigate the tendency of these two parameters along with their iterations. At the same time, we can quantitatively record the image quality, e.g., the mean squared error (MSE) compared to the ground truth. By doing so, we can investigate whether there is a correlation between the two parameters and image quality in the simulation, which always contains the ground truth. Based on this correlation, the real sinogram data experiment can also be evaluated.

2.4. Stability Investigation

Figure 1 compares the difference between the conventional MAP (left) and joint-parameter-Bayes (right). The orange region is one starting point. The star indicates the position of the optimal solution. For the conventional MAP method, we manually set one smoothness strength (hyperparameter) each time. After several trials, a set of solutions can be obtained. Then, the optimal solution is determined by comparing the final image quality. The proposed joint-parameter-Bayes method is more like that shown on the right hand. At every stage, it gives an estimation of the solution distribution and gradually approaches the optimal solution. Comparing the two methods, the conventional MAP is seen to be very time-consuming. The joint-parameter-Bayes method is much more convenient. However, the joint-parameter-Bayes method might give the “too far away” estimation and may lead to unsatisfying results. Therefore, the stability of the proposed joint-parameter-Bayes method should be explored from different factors, such as initials, variance and so on.

2.5. Phantom Simulation and Patient Data Acquisition

One advantage of using numerical phantom is that the simulation is fully controlled and the ground truth is available. Thus, we used the classical Shepp–Logan phantom to explore the effectiveness and stability of the proposed joint-parameter-Bayes method. The Shepp–Logan phantom includes 512 × 512 pixels with a resolution of 0.74 mm × 0.74 mm. This simulation setup is the same as that for the clinical patient data used in this work. Specifically, the distance from the source to the system center was 570 mm, and the distance from the source to the detector was 1040 mm. The sinogram data has 672 detector elements with a width of 1.4 mm per element, and 1160 projection views over a 360° range by forward projection. A clinical patient was also recruited to this study under informed consent after approval by the Institutional Review Board (IRB). The patient was scanned with normal dose settings with a tube voltage of 120 kVp and a tube current of 100 mAs. The geometric settings are consistent with the simulation study. Since the sinogram is from a normal-dose scan, we just used the reconstruction by FBP as reference for comparison.

3. Results

3.1. Results of Numerical Simulation Data

3.1.1. Reconstruction Comparison

In Section 2.1, the multivariate Gaussian distribution data model is described in detail. Note that this variance model analysis is based on the realistic CT signal generation model, i.e., the compound of the Gaussian and Poisson model. In a simplified simulation case, the variance can be degenerated to be a scaler when directly applying white Gaussian noise to the post-log sinogram. With simulations, we can manipulate the noise type to perform comprehensive evaluations.
The simple noise model is used to add the white Gaussian noise directly to the post-log sinogram data. This provides a prefect evaluation model, since we know exactly the value of s . If we add noise ζ ~ G a u s s i a n ( 0 ,   σ G 2 ) , s will be of the value σ G 2 by assigning ρ i 2 = 1 . We can evaluate the estimation of s ^ by the proposed method. In the experiment, five noise levels are simulated with σ G 2 = { 0.1 2 , 0.3 2 , 0.5 2 , 0.7 2 , 0.9 2 , 1.1 2 } . For a comparison of the reconstructed images, the conventional MAP method is regarded as the baseline. For the baseline, we swept its hyperparameter within a wide range by sampling three points at each order of magnitude. For each sampled hyperparameter, the reconstruction is stopped after 1000 iterations. If we denote the reconstruction for one sampled hyperparameter as one trial, the computing time for the conventional method is T i m e M A P = t i m e / i t e r a t i o n   * # i t e r a t i o n s *   t r i a l   n u m b e r . It is ensured that the reconstructed images are from too noisy to oversmoothed. Visual judgement and quantitative measures are used to determine the optimal images as the results of the baseline.
Figure 2 shows the reconstructed images both by the conventional MAP (baseline) and joint-parameter-Bayes with different noise levels. It can be observed that the results are comparable. The quantitative comparisons are summarized in Table 1 which are also consistent with our visual inspection. However, their computation time is quite different based on the iterations they use in Table 1. The proposed joint-parameter-Bayes saves approximately 66~88% of computation time in comparison with the baseline method. The number of iterations for joint-parameter-Bayes is further discussed in the stopping criterion subsection. To test the robustness, we performed the experiments on two other noise types: (1) Poisson-like signal-dependent noise to post-log sinogram data; (2) Poisson noise to the pre-log transmission data. We repeated the above evaluations and yielded very similar results, which can be found in the Appendix A. Some details of the two noise models are also be discussed in the following subsections.
To evaluate the estimation of joint-parameter-Bayes for the unobservable parameters, we also recorded the value of s with the iteration times for a different simulation. The results are shown in Figure 3. The red line is the simulated value or true value of s . The blue curve is the estimated value by the joint-parameter-Bayes method. It is found that the estimated s agrees well with the true value, which demonstrates the effectiveness of the estimation method. It is also noted there is slightly higher disagreement when the iteration passes a certain number, especially in the high-noise scenarios, e.g., σ G 2 = 1.1 2 . This observation also indicates that it might be better to terminate the iteration early.

3.1.2. Stopping Criterion Investigation

As described in Section 2.3, the stopping position (i.e., different number of iterations) may result in different image quality. In this section, we investigate the stopping criterion for the proposed joint-parameter-Bayes by exploring the correlation between the two unobservable parameters and image quality. To see the response of the stopping criterion to different noises, three noise types mentioned above are added to the simulation data, including white Gaussian noise to post-log sinogram data, Poisson-like signal-dependent noise to post-log sinogram data, and Poisson noise to the pre-log transmission data. As introduced in Section 3.1.1, we can manipulate the variance σ G 2 to mimic different noise levels for white Gaussian noise. For Poisson noise to sinogram data, we scale the values of Shepp–Logan phantom by different factors and then add Poisson noise to its sinogram data from forward projection. By doing so, we can mimic different noise levels by varying the factor value. Three scaling factors are used, which are 1, 5 and 10. For the Poisson noise to transmission data, it is more like the realistic noise model introduced in Section 2.1.1. By varying the incident flux N 0 , we simulated different noise levels. In this work, we chose five incident flux, N 0 = { 1 × 10 4 , 5 × 10 4 , 1 × 10 5 , 5 × 10 5 , 1 × 10 6 } .
Figure 4 shows the change of two unobservable parameters along with their iteration. The image quality measure RMSE is also plotted in the same figure to find some correlations. It can be found the unobservable parameter t has a coincident turning point with the optimal image quality. For example, in the case of σ G 2 = 0.7 2 , the image quality (red line) has a valley around 250 iterations. The s curve has a very fast drop at several iterations and then becomes flat. The t curve has a relatively slow drop and then becomes flat. If we call the largest slope rate change point the turning point, the turning point of image quality and t are almost identical. Therefore, we can monitor the transformation of t along with iteration and terminate the reconstruction at its turning point, since the image quality measure RMSE cannot be obtained in real cases. For the other two noise types, the above observations are found in the case of scale factor = 1. The change of s remains the same, but the t curve decreases first and then increases. Interestingly, the image quality also deteriorates first and then improves. After the changing point, both image quality and t tend to converge. One possible reason for the initial un-monotonic tendency could be the initial condition. This effect can be removed, however, if we start to monitor the value of t several iterations (e.g., 10 iterations) late. After the removal of the initial effect, the regulations can be summarized as if t decreases monotonically, the reconstruction can stop at the turning point, and if t increases, the reconstruction can stop when it is converged. It is interesting to find that the prior variance parameter t has such an explicit correlation with the image quality. Data variance parameter s has no such clear tendency. One possibility could be that the prior variance t is directly related to the image domain and data variance parameter s is directly related to the measurement. Therefore, the relationship between image quality and t is not explicit. It is also noted that the t curves have different patterns under different noise conditions, regardless of the initial effect. This is expected because the data-independent, data-dependent and logarithm non-linearity type noises are introduced on purpose. The three noise types cover the most popular models used in the image reconstruction. Based on the observation, we find some semi-experimental regulation of the early stop criteria for optimal image reconstruction. To fully understand the curve patterns under different noise conditions is one of our future research interests.
In the experiments, it is noticed that the stopping criterion should be part of the joint-parameter-Bayes method. Ignoring the stopping criterion may result in unsatisfying results. We presented the reconstructed images of early stop criterion and hard threshold in Figure 5. We can see that the early stop can produce better reconstructed images.

3.1.3. Stability Investigation

(a)
Effect of Initialization
We further investigated the effect of initialization on the proposed joint-parameter-Bayes method. Three different initials were tested: uniform initial by averaging the FBP results, FBP with ramp filter, and FBP with hann50 filter. We investigated the effect on the Poisson-like signal-dependent noise to post-log sinogram data. This simulation provides the opportunity to deal with the signal-dependent noise data, where the signal variance is known to equal its signal mean. This means that the true value of ρ 2 is y . First, we used the true variance y for ρ 2 , instead of the estimated y ^ and just varied the initials. The reconstructed images and their initials are presented in Figure 6. The estimated parameters along the iterations are presented in Figure 7. According to Figure 7, the pattern of s does not change too much for different initials, but the pattern of t varies. Their ratio s / t follows a similar pattern with a different convergence rate. According to Figure 6, the reconstructed results are stable with all three different initials. This implies that joint-parameter-Bayes can effectively estimate the parameter values starting with different initials and provide stable solutions anyway when knowing the true variance.
We then tested the initials by only knowing that the variance distribution is Poisson. Specifically, ρ 2 is estimated by the A μ ^ and updates every iteration. Stable reconstructions are still observed with FBP initials. However, something interesting is observed for the uniform initials. Figure 8a presents the reconstructed images at different iterations with the same uniform initial used in Figure 6. Comparing Figure 6 and Figure 8a, the only difference is the variance value ρ 2 . Therefore, we recorded the change of ρ 2 estimated by A μ ^ at different iterations, shown as Figure 8c. In the simulation, we know that the ground truth of variance should have the pattern of the first image in Figure 8d. However, the estimated variance in Figure 8c cannot produce the right estimation. To further investigate the effect, we directly use the ground truth variance as the initial variance and still use A μ ^ to estimate the variance for the following iterations. The tracked reconstructions and variance estimation are shown Figure 8b,d. According to Figure 8b, joint-parameter-Bayes can somehow give us a reasonable reconstruction. From Figure 8d, although the variance is not estimated very accurately during the iteration, finally, it can return by giving a good initial. Therefore, this joint-parameter-Bayes depends on the initial values of the parameters. In the experiments, it is found that the initials using FBP are quite robust, and they are recommended for joint-parameter-Bayes.
(b)
Effects of Variance Normalization
As we know, ρ 2 is the variance factor from the measurement. The element value range might vary from different modalities. For example, in the simulation, ρ i 2 = y i for the Poisson-like signal-dependent noise to post-log sinogram data, and ρ i 2 = N i for the Poisson noise to the pre-log transmission data. Therefore, it is useful to know whether the normalization of matrix ρ 2 will affect the performance of the joint-parameter-Bayes method. Figure 9 presents the reconstructed images for ρ 2 without (upper) and with normalization (lower) by s u m ( ρ 2 ) = 1 on the Poisson-like signal-dependent noise to post-log sinogram data. The scale factor is 1. Comparable reconstructed images are obtained. The estimated parameters follow the same pattern with linear scaling strength. This observation agrees with our expectation that even though the magnitude of s   a n d   t are scaled due to normalization, the fidelity strength and prior strength will also be scaled by the same degree. Therefore, the results would not be affected.

3.2. Results of Clinical Patient Data

Figure 10 shows the results of the reconstructed images using different algorithms. The pattern of t along with iteration is also shown in the top right panel. The early stopping criterion used is the turning point of t, as described in Section 3.1.2. From the t curve, we can see that it has a rapid decrease and then becomes flat. According to the experimental study in Section 3.1.2, for this type of t curve pattern, the turning point of t—which is the largest slope rate change point—will give the best image quality. This reconstruction should stop early at this iteration number. The reconstruction with early stopping criteria is shown in the bottom left panel. For comparison, the reconstruction with 1000 iterations is shown in the bottom right panel. According to Figure 10, joint-parameter-Bayes can give a comparable image with the reference using visual judgement. With the convergence stopping criterion, the reconstructed image of joint-parameter-Bayes is blurred. Early stopping is necessary to obtain satisfying results. The stopping criterion obtained from the simulation model has also shown its effectiveness in clinical patient data.

4. Discussion

We performed extensive studies to investigate the joint-parameter-Bayes method without a freely adjustable hyperparameter. It is found that an early stopping criterion is necessary for optimal image quality. The stopping criterion constructed from one of the unobservable parameters has shown its effectiveness and robustness through three different noise models and evaluation on the clinical patient data. However, this criterion observed is based on the initial image being FBP-ramp. More evaluations are needed to investigate the stopping criterion under different conditions. Currently, we construct the criterion based on our observations in the experiments. It is still unclear why the two unobservable parameters have the patterns they do and why the parameter of the prior is more sensitive.
In the conventional MAP method, the hyperparameter is used to balance the data fidelity term and the prior/penalty term, which is a constant through the reconstruction process. Comparing the updating equation of the attenuation map between joint-parameter-Bayes (see the pseudo code) and the conventional MAP method [23], we can see that the hyperparameter of the conventional method is proportional to the ratio of unobservable parameters s and t, where there might be some scale factor considering the normalization effect. As we know, s is introduced to completely specify the data fidelity PDF, and t is introduced to completely specify the prior term PDF. s is related to the data fidelity variance and t is related to the prior variance. This might be one interpretation of the physical meaning of the hyperparameter, which is related to the ratio of two unobservable parameters reflecting the data fidelity variance and the prior variance. This implies that, given the count-limited imaging, such as LdCT, we desire high confidence or low variance on the prior penalty. In other words, we desire the medical truth to derive prior and train machine learning. Additionally, in the reconstruction of joint-parameter-Bayes, the s/t is no longer a constant, which will be self-adjusted by the current reconstruction estimation. In the future, we will compare the proposed method with the deep learning-based approaches of estimating the optimal hyperparameter for the MAP methods.
There is another remaining issue to be addressed. This work used the normal-dose clinical patient data to evaluate the performance of the joint-parameter-Bayes algorithm. In the future, evaluations on the LdCT and other modalities are needed to validate the model.

5. Conclusions

A theory-based joint-parameter-Bayes method is presented in this work to avoid the freely adjustable hyperparameter for penalized maximum likelihood tomographic image reconstruction. The strategy is based on the PDF requirements of Bayes’ law. For both fidelity and prior term in Bayesian reconstruction, their multivariate Gaussian PDF distributions are completely described by two unobservable parameters, which are jointly estimated by maximizing a posteriori probability for the proposed Bayesian reconstruction. The optimal image can be reconstructed by joint-parameter-Bayes without any freely adjustable hyperparameter for the given data, by satisfying the PDF condition for both the data likelihood and the prior probability of Bayes’ law, and by satisfying the stopping criterion. Compared to the conventional MAP method, the proposed scheme can save at least 66% computation time. It has great potential for application to joint-parameter-Bayes to improve the reconstructed image quality with incomplete data, such as LdCT and limited angle CT. It can also be generalized to other tomographic modalities and applications, such as positron emission tomography and ultrasound imaging reconstruction.

Author Contributions

Conceptualization, Z.L.; Methodology, Y.G., Y.S., S.C. and Z.L.; Validation, Y.G., W.H. and L.L.; Formal analysis, W.H. and L.L.; Investigation, Y.G., S.L., Y.S., S.C., H.Z., L.L. and Z.L.; Resources, Z.L.; Data curation, H.Z. and Z.L.; Writing—original draft, Y.G. and Z.L.; Writing—review & editing, S.L., Y.S., S.C., H.Z., W.H., L.L. and Z.L. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the NIH/NCI grant #CA206171.

Institutional Review Board Statement

The study was conducted in accordance with the Declaration of Helsinki, and approved by the Institutional Review Board (or Ethics Committee) of Stony Brook University.

Informed Consent Statement

Informed consent was obtgained from all subjects involved in the study.

Data Availability Statement

Patient sinogram data is unavailable due to privacy or ethical restrictions.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Appendix A

Figure A1 and Table A1 present the comparison results between the baseline and the proposed joint-parameter-Bayes on Poisson noise in post-log sinogram data. Figure A2 and Table A2 present the comparison results between the baseline and the proposed joint-parameter-Bayes on Poisson noise to pre-log transmission data.
Figure A1. Comparison of reconstructed images between the conventional MAP (baseline) and joint-parameter-Bayes with Poisson noise to post-log sinogram data.
Figure A1. Comparison of reconstructed images between the conventional MAP (baseline) and joint-parameter-Bayes with Poisson noise to post-log sinogram data.
Sensors 23 01374 g0a1
Table A1. Quantitative comparison of joint-parameter-Bayes with conventional MAP (baseline) with Poisson noise to post-log sinogram data.
Table A1. Quantitative comparison of joint-parameter-Bayes with conventional MAP (baseline) with Poisson noise to post-log sinogram data.
Conventional MAP (Baseline)Joint-Parameter-Bayes
Scale FactorRMSESSIMPSNRIterationsRMSESSIMPSNRIterations
1 0.06430.910923.841000 * trial number0.05470.824925.2481
5 0.18440.644514.681000 * trial number0.18480.635814.671000
10 0.29940.609110.481000 * trial number0.30730.590410.251000
Figure A2. Comparison of reconstructed images between the conventional MAP (baseline) and joint-parameter-Bayes with Poisson noise to pre-log transmission data.
Figure A2. Comparison of reconstructed images between the conventional MAP (baseline) and joint-parameter-Bayes with Poisson noise to pre-log transmission data.
Sensors 23 01374 g0a2
Table A2. Quantitative comparison of joint-parameter-Bayes with conventional MAP (baseline) with Poisson noise to pre-log transmission data.
Table A2. Quantitative comparison of joint-parameter-Bayes with conventional MAP (baseline) with Poisson noise to pre-log transmission data.
Conventional MAP (Baseline)Joint-Parameter-Bayes
Incident FluxRMSESSIMPSNRIterationsRMSESSIMPSNRIterations
1 × 10 4 0.00650.974443.691000 * trial number0.00670.973143.531000
5 × 10 4 0.00480.986446.361000 * trial number0.00510.985345.831000
1 × 10 5 0.00420.989847.461000 * trial number0.00460.988546.731000
5 × 10 5 0.00330.994149.681000 * trial number0.00380.993148.251000
1 × 10 6 0.00310.994650.041000 * trial number0.00370.993748.531000

References

  1. Schwenzer, N.F.; Springer, F.; Schraml, C.; Stefan, N.; Machann, J.; Schick, F. Non-invasive assessment and quantification of liver steatosis by ultrasound, computed tomography and magnetic resonance. J. Hepatol. 2009, 51, 433–445. [Google Scholar] [CrossRef] [PubMed]
  2. Deans, S.R. The Radon Transform and Some of Its Applications; Courier Corporation: Chelmsford, MA, USA, 2007. [Google Scholar]
  3. Shepp, L.A.; Logan, B.F. The Fourier reconstruction of a head section. IEEE Trans. Nucl. Sci. 1974, 21, 21–43. [Google Scholar] [CrossRef]
  4. Liang, Z.; La Riviere, P.; El Fakhri, G.; Glick, S.; Siewerdsen, J. Low-Dose CT: What Has Been Done, and What Challenges Remain. IEEE Trans. Med. Imaging 2017, 36, 2409–2416. [Google Scholar] [CrossRef]
  5. Jin, B.; Zhao, Y.; Zou, J. Iterative parameter choice by discrepancy principle. IMA J. Numer. Anal. 2012, 32, 1714–1732. [Google Scholar] [CrossRef]
  6. Golub, G.H.; Wahba, G. Generalized cross-validation as a method for choosing a good ridge parameter. Technometrics 1979, 21, 215–223. [Google Scholar] [CrossRef]
  7. Hansen, P.C. Analysis of discrete ill-posed problems by means of the L-curve. SIAM Rev. 1992, 34, 561–580. [Google Scholar] [CrossRef]
  8. Kazufumi, I.; Jin, B.; Takeuchi, T. A regularization parameter for nonsmooth Tikhonov regularization. SIAM J. Sci. Comput. 2011, 33, 1415–1438. [Google Scholar]
  9. Bai, T.; Yan, H.; Jia, X.; Jiang, S.; Wang, G.; Mou, X. Z-index parameterization for volumetric CT image reconstruction via 3-D dictionary learning. IEEE Trans. Med. Imaging 2017, 36, 2466–2478. [Google Scholar] [CrossRef]
  10. Mou, X.; Bai, T.; Chen, X.; Yu, H.; Yang, Q.; Wang, G. Optimal Selection for Regularization Parameter in Iterative CT Reconstruction Based on the Property of Natural Image Statistics. Fully 3D 2015, 13, 383–386. [Google Scholar]
  11. Johnson, V.E.; Wong, W.H.; Hu, X.; Chen, C.T. Image restoration using Gibbs priors: Boundary modeling, treatment of blurring, and selection of hyperparameter. IEEE Trans. Pattern Anal. Mach. Intell. 1991, 5, 413–425. [Google Scholar] [CrossRef]
  12. Hsiao, T.; Rangarajan, A.; Gindi, G.R. Bayesian image reconstruction for transmission tomography using deterministic annealing. J. Electron. Imaging 2003, 12, 7–17. [Google Scholar]
  13. Cai, C.; Rodet, T.; Legoupil, S.; Mohammad-Djafari, A. A full-spectral Bayesian reconstruction approach based on the material decomposition model applied in dual-energy computed tomography. Med. Phys. 2013, 40, 111916. [Google Scholar] [CrossRef]
  14. Zhang, G.; Liu, F.; Liu, J.; Luo, J.; Xie, Y.; Bai, J.; Xing, L. Cone beam X-ray luminescence computed tomography based on Bayesian method. IEEE Trans. Med. Imaging 2016, 36, 225–235. [Google Scholar] [CrossRef]
  15. Veklerov, E.; Llacer, J. Stopping rule for the MLE algorithm based on statistical hypothesis testing. IEEE Trans. Med. Imaging 1987, 6, 313–319. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  16. Herbert, T.J. Statistical stopping criteria for iterative maximum likelihood reconstruction of emission images. Phys. Med. Biol. 1990, 35, 1221. [Google Scholar] [CrossRef]
  17. Gaitanis, A.; Kontaxakis, G.; Spyrou, G.; Panayiotakis, G.; Tzanakos, G. PET image reconstruction: A stopping rule for the MLEM algorithm based on properties of the updating coefficients. Comput. Med. Imaging Graph. 2010, 34, 131–141. [Google Scholar] [CrossRef] [PubMed]
  18. Gao, Y.; Liang, Z.; Lu, S.; Shi, Y.; Chang, S.; Hou, W. Hyperparameter Selection for Bayesian Image Reconstruction by Mimicking Physical Crystallization. In Proceedings of the 2020 IEEE Symposium on Nuclear Science (NSS/MIC), Boston, MA, USA, 31 October–7 November 2020. [Google Scholar]
  19. Gao, Y.; Liang, Z.; Xing, Y.; Zhang, H.; Pomeroy, M.; Lu, S.; Ma, J.; Lu, H.; Moore, W. Characterization of tissue-specific pre-log Bayesian CT reconstruction by texture–dose relationship. Med. Phys. 2020, 47, 5032–5047. [Google Scholar] [CrossRef]
  20. Ma, J.; Liang, Z.; Fan, Y.; Liu, Y.; Huang, J.; Chen, W.; Lu, H. Variance Analysis of X-ray CT Sinograms in the Presence of Electronic Noise Background. Med. Phys. 2012, 39, 4051–4065. [Google Scholar] [CrossRef] [Green Version]
  21. Wang, J.; Li, T.; Lu, H.; Liang, Z. Penalized Weighted Least-squares Approach to Sinogram Noise Reduction and Image Reconstruction for Low-dose X-ray CT. IEEE Trans. Med. Imaging 2006, 25, 1272–1283. [Google Scholar] [CrossRef] [PubMed]
  22. Zhang, H.; Han, H.; Wang, J.; Ma, J.; Liu, Y.; Moore, W.; Liang, Z. Deriving Adaptive MRF Coefficients from Previous Normal-dose CT Scan for Low-dose Image Reconstruction via Penalized Weighted Least-squares Minimization. Med. Phys. 2014, 41, 041916. [Google Scholar] [CrossRef] [Green Version]
  23. Zhang, H.; Han, H.; Liang, Z.; Hu, Y.; Moore, W.; Ma, J.; Lu, H. Extracting Information from Previous Full-dose CT Scan for Knowledge-based Bayesian Reconstruction of Current Low-dose CT Images. IEEE Trans. Med. Imaging 2016, 35, 860–870. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  24. Gao, Y.; Liang, Z.; Moore, W.; Zhang, H.; Pomeroy, M.J.; Ferretti, J.A.; Bilfinger, T.V.; Ma, J.; Lu, H. A Feasibility Study of Extracting Tissue Textures from a Previous Full-Dose CT Database as Prior Knowledge for Bayesian Reconstruction of Current Low-dose CT Images. IEEE Trans. Med. Imaging 2019, 38, 1981–1992. [Google Scholar] [CrossRef] [PubMed]
Figure 1. Comparison between the conventional MAP (left) and the proposed joint-parameter-Bayes (right).
Figure 1. Comparison between the conventional MAP (left) and the proposed joint-parameter-Bayes (right).
Sensors 23 01374 g001
Figure 2. Comparison of reconstructed images between the conventional MAP (top row) and joint-parameter-Bayes (bottom row). The display window is [0, 0.035] mm−1.
Figure 2. Comparison of reconstructed images between the conventional MAP (top row) and joint-parameter-Bayes (bottom row). The display window is [0, 0.035] mm−1.
Sensors 23 01374 g002
Figure 3. Comparison of estimated s by the joint-parameter-Bayes method (blue curve) and the true value of s (red curve). They agree well with each other, demonstrating the effectiveness of our proposed estimation method.
Figure 3. Comparison of estimated s by the joint-parameter-Bayes method (blue curve) and the true value of s (red curve). They agree well with each other, demonstrating the effectiveness of our proposed estimation method.
Sensors 23 01374 g003
Figure 4. Stopping criterion under three different noise types. The change of two unobservable parameters ( s ,   t ) and change of image quality (RMSE) along with iteration are plotted.
Figure 4. Stopping criterion under three different noise types. The change of two unobservable parameters ( s ,   t ) and change of image quality (RMSE) along with iteration are plotted.
Sensors 23 01374 g004
Figure 5. Reconstructed images by FBP (left column) and the joint-parameter-Bayes method with early stop (middle column) and hard threshold stop (right column). The display window is [0, 0.035] mm−1.
Figure 5. Reconstructed images by FBP (left column) and the joint-parameter-Bayes method with early stop (middle column) and hard threshold stop (right column). The display window is [0, 0.035] mm−1.
Sensors 23 01374 g005
Figure 6. Performance of joint-parameter-Bayes with different initials. From top to bottom, the initials (left column) are uniform averaged from FBP, FBP with ramp filter, and FBP with hann50 filter, while the corresponding joint-parameter-Bayes reconstructions are shown on the right (right column).
Figure 6. Performance of joint-parameter-Bayes with different initials. From top to bottom, the initials (left column) are uniform averaged from FBP, FBP with ramp filter, and FBP with hann50 filter, while the corresponding joint-parameter-Bayes reconstructions are shown on the right (right column).
Sensors 23 01374 g006
Figure 7. Estimated free parameters by joint-parameter-Bayes along with iterations with three different initials in accordance with Figure 6.
Figure 7. Estimated free parameters by joint-parameter-Bayes along with iterations with three different initials in accordance with Figure 6.
Sensors 23 01374 g007
Figure 8. (a,b): Reconstructed images at different iterations with uniform initials but different initial variance. (c,d): The change of variance ρ 2 at different iterations according to (a,b).
Figure 8. (a,b): Reconstructed images at different iterations with uniform initials but different initial variance. (c,d): The change of variance ρ 2 at different iterations according to (a,b).
Sensors 23 01374 g008
Figure 9. Performance of joint-parameter-Bayes images for ρ 2 without (upper) and with normalization (lower) by s u m ( ρ 2 ) = 1 . The CT image display window is [0, 0.035] mm−1.
Figure 9. Performance of joint-parameter-Bayes images for ρ 2 without (upper) and with normalization (lower) by s u m ( ρ 2 ) = 1 . The CT image display window is [0, 0.035] mm−1.
Sensors 23 01374 g009
Figure 10. Comparison of reconstructed images between the conventional MAP and joint-parameter-Bayes. The CT image display window is [0, 0.035] mm−1.
Figure 10. Comparison of reconstructed images between the conventional MAP and joint-parameter-Bayes. The CT image display window is [0, 0.035] mm−1.
Sensors 23 01374 g010
Table 1. Quantitative comparison of the conventional MAP (baseline) and joint-parameter-Bayes.
Table 1. Quantitative comparison of the conventional MAP (baseline) and joint-parameter-Bayes.
Conventional MAP (Baseline)Joint-Parameter-Bayes
σ G 2 RMSESSIMPSNRIterationsRMSESSIMPSNRIterations
0.1 2 0.0023970.997052.411000 * trial number0.0022640.997252.901000
0.3 2 0.0035810.991148.921000 * trial number0.0041250.991847.69339
0.5 2 0.0044400.989447.051000 * trial number0.0047780.988146.42200
0.7 2 0.0049950.986746.031000 * trial number0.005180.985345.71157
0.9 2 0.0051900.984445.671000 * trial number0.0055010.982945.19132
1.1 2 0.0054260.981545.311000 * trial number0.0057860.980544.75115
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Gao, Y.; Lu, S.; Shi, Y.; Chang, S.; Zhang, H.; Hou, W.; Li, L.; Liang, Z. A Joint-Parameter Estimation and Bayesian Reconstruction Approach to Low-Dose CT. Sensors 2023, 23, 1374. https://doi.org/10.3390/s23031374

AMA Style

Gao Y, Lu S, Shi Y, Chang S, Zhang H, Hou W, Li L, Liang Z. A Joint-Parameter Estimation and Bayesian Reconstruction Approach to Low-Dose CT. Sensors. 2023; 23(3):1374. https://doi.org/10.3390/s23031374

Chicago/Turabian Style

Gao, Yongfeng, Siming Lu, Yongyi Shi, Shaojie Chang, Hao Zhang, Wei Hou, Lihong Li, and Zhengrong Liang. 2023. "A Joint-Parameter Estimation and Bayesian Reconstruction Approach to Low-Dose CT" Sensors 23, no. 3: 1374. https://doi.org/10.3390/s23031374

APA Style

Gao, Y., Lu, S., Shi, Y., Chang, S., Zhang, H., Hou, W., Li, L., & Liang, Z. (2023). A Joint-Parameter Estimation and Bayesian Reconstruction Approach to Low-Dose CT. Sensors, 23(3), 1374. https://doi.org/10.3390/s23031374

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop