==== Front ArXiv ArXiv arxiv ArXiv 2331-8422 Cornell University arXiv:2306.08610v7 2306.08610 7 preprint Article Tomographic Image Reconstruction Using an Advanced Score Function (ADSF) Cong Wenxiang Xia Wenjun Wang Ge Biomedical Imaging Center, Center for Computational Innovations, Center for Biotechnology and Interdisciplinary Studies, Department of Biomedical Engineering, Rensselaer Polytechnic Institute, Troy, NY 12180 29 2 2024 arXiv:2306.08610v7https://creativecommons.org/licenses/by/4.0/ This work is licensed under a Creative Commons Attribution 4.0 International License, which allows reusers to distribute, remix, adapt, and build upon the material in any medium or format, so long as attribution is given to the creator. The license allows for commercial use. nihpp-2306.08610v7.pdf Computed tomography (CT) reconstructs volumetric images using X-ray projection data acquired from multiple angles around an object. For low-dose or sparse-view CT scans, the classic image reconstruction algorithms often produce severe noise and artifacts. To address this issue, we develop a novel iterative image reconstruction method based on maximum a posteriori (MAP) estimation. In the MAP framework, the score function, i.e., the gradient of the logarithmic probability density distribution, plays a crucial role as an image prior in the iterative image reconstruction process. By leveraging the Gaussian mixture model, we derive a novel score matching formula to establish an advanced score function (ADSF) through deep learning. Integrating the new ADSF into the image reconstruction process, a new ADSF iterative reconstruction method is developed to improve image reconstruction quality. The convergence of the ADSF iterative reconstruction algorithm is proven through mathematical analysis. The performance of the ADSF reconstruction method is also evaluated on both public medical image datasets and clinical raw CT datasets. Our results show that the ADSF reconstruction method can achieve better denoising and deblurring effects than the state-of-the-art reconstruction methods, showing excellent generalizability and stability. Computed tomography (CT) radiation dose reduction image reconstruction maximum a posteriori (MAP) estimation Gaussian mixture model score function deep learning ==== Body pmcX-ray computed tomography (CT) is the most popular imaging modality used in various fields, including medical imaging, homeland security, and industrial applications [1]. Filtered backprojection (FBP) is an analytical method for tomographic image reconstruction [2], and often produces severe noise and artifacts in low-dose or sparse-view settings [3]. To address this challenge, iterative reconstruction methods were developed to incorporate the photons statistics and desirable image prior models especially sparsity [3–5]. The key issue in this iterative reconstruction method is the image prior modeling. A representative method based on compressed sensing (CS) converts images into sparse data and then performs image reconstruction through sparse regularization such as total variation (TV) minimization [4]. However, CS-based image reconstruction tends to over smoothen textures and eliminate subtle details in a reconstructed image. To overcome this shortcoming, low dimensional manifold model (LDMM) was developed for image reconstruction utilizing the low dimensional characteristics of the image patch manifold [6, 7]. However, these presumptions did not accurately and comprehensively reflect actual image structures, especially for sophisticated biomedical images. Deep learning is a powerful approach to perform various types of uncertainty estimation and data modeling [8]. Supervised learning was applied to perform image denoising and deburring by establishing a convolution neural network-based image processing model [9], and implement image reconstruction by replacing regularization terms with a trained convolutional neural network (CNN) within the framework of an unrolling iterative optimization scheme [10, 11]. The supervised learning requires a labeled dataset pairing each input sample with the corresponding target, which is usually unavailable in practice. Recently, diffusion models or score matching models have attracted a major attention for image generation [12] and image reconstruction [13, 14]. In an unsupervised fashion, the score matching methods learn a score function, i.e., the gradient of the logarithmic probability density function, from a training data set without any labels [13]. In contrast to conventional regularized reconstruction methods that aim to constrain images to become sparse, low-rank, and low-dimensional, the scoring function can generate an optimal image prior using the data-driven approach, allowing an effective representation of the image distribution. In this paper, we model the number of X-ray photon reaching each pixel on the detector as a Poisson distribution, and perform CT image reconstruction using maximum a posteriori (MAP) estimation. In the MAP framework, the score function plays a crucial role as an image prior in the iterative image reconstruction process. By leveraging Gaussian mixture model to characterize the discrepancy between a reconstructed image and the target image, we derive a novel score matching formula to generate an advanced score function (ADSF) through deep learning. By integrating a new ADSF into the image reconstruction process, we develop a new iterative reconstruction method for image reconstruction. The convergence of the iterative reconstruction algorithm is proven through mathematical analysis. The performance of proposed ADSF reconstruction method is also evaluated on both public medical image datasets and clinical raw datasets, demonstrating the superiority of our approach over state-of-the-art image reconstruction methods in terms of accuracy and generalizability. Results Based on the benchmark dataset used in the NIH-AAPM-Mayo Clinic low-dose CT Challenge, the advanced score function (ADSF) was established through deep learning according to the new score matching formula Eq. (16) in theMethods section. Using the obtained ADSF, we conducted extensive experiments on public medical image datasets and clinical CT raw datasets to evaluate the performance of the ADSF reconstruction method in comparison with the representative image reconstruction methods, including the filtered backprojection (FBP) [1], model-based iterative reconstruction (MBIR) with total variation minimization (MBIR-TV) [4, 15], and state-of-the-art score function based image reconstruction (SFrecon) [13, 14]. Phantoms CT images from the NIH-AAPM-Mayo Clinic CT Grand Challenge were selected as realistic digital phantoms to evaluate the performance of the ADSF reconstruction method. Specifically, the dataset contains 2,378 images from 10 patients. Image phantoms were randomly selected from 2,378 images. In the X-ray imaging simulation of the image phantoms, the scan trajectory radius was 595mm. The distance from the source to the detector was 1085.6mm. There were 736 detector elements with 1.2858mm pitch, equiangularly distributed along a curvilinear array. Fan-beam projections of X-ray imaging were generated using an industrial CT simulator called CatSim at 120kV and 100mAs. A total of 360 projection views were uniformly acquired over a 360 degree range. We performed fan-beam image reconstruction for 100 digital phantoms using FBP, MBIR-TV, SFrecon, and our ADSF reconstruction methods, respectively. The reconstructed images were quantitatively evaluated using peak signal-to-noise ratio (PSNR) and structural similarity index (SSIM) metrics. Table 1 presents the quantitative results of the reconstructed images with respect to PSNR and SSIM, showing that the ADSF reconstruction method achieved higher average PSNR and SSIM values than other competing reconstruction methods. Figure 1 displays the representative images and zoom-in patches for visualization. By comparing reconstructed images, it can be observed that the image reconstructed using the FBP is noisier than the image reconstructed by the other reconstruction methods in the low-dose and few-view scenarios, while the ADSF reconstruction method achieves excellent denoising and deblurring effects and well preserves structural information. The qualitative and quantitative results show that the proposed ADSF reconstruction method outperforms the other reconstruction methods. Clinical CT Raw Datasets Clinical raw data from the Siemens scanner: NIH-AAPM-Mayo Clinic Low-dose CT Grand Challenge contained 10 anonymized patient normal dose abdominal CT images. The dataset was acquired on the Siemens SOMATOM Definition Flash CT system with voltage 120 kV and 200 mAs. The scanner had a scanning radius of 595mm. The distance from the source to the detector was 1085.6mm. There were 736 detector elements with 1.2858mm pitch, equiangularly distributed along a curvilinear array. Over a 360 degree range, 768 projection views were uniformly acquired for image reconstruction. Here, the helical raw data were converted to flat detector fan-beam projections. Based on the fan-beam projection dataset, we performed 100 image reconstructions from one-third projection views using FBP, MBIR-TV, SFrecon, and ADSF reconstruction methods, respectively. Table 2 shows the quantitative results of the reconstructed images relative to the reference image reconstructed from the 2304 full-view projections using FBP. Figure 2 displays the representative images, showing better performance in terms of noise removal as well as preservation of weak edges and exhibit structural information. Clinical raw data from the GE scanner: Furthermore, a clinical raw dataset were also acquired from GE CT scanner. The scanner had a scanning radius of 541mm. The distance from source to detector was 949mm. There were 888 detector elements with 1.024mm pitch, equiangularly distributed along a curvilinear array. 984 projection views were uniformly acquired over a 360 degree range. Here, we uniformly selected 328 projection views for few-view image reconstruction. Table 3 gives a quantitative comparison of the reconstructed images relative to the reference image reconstructed from 984 full-view projections using FBP. From image reconstruction results shown in Figure 3, it is clear that FBP, MBIR-TV, and SFrecon methods still exhibit blurry boundaries and unremoved noise, while our ADSF reconstruction method is quite robust in denoising and deblurring. Discussion Based on X-ray imaging physics, the number of photons reaching individual pixels on the detector is model by a Poisson distribution, and CT image reconstruction is performed via maximum a posteriori (MAP) estimation. Under the MAP estimation framework, the CT images can be reconstructed by maximizing likelihood of the Poisson statistics and an image prior probability density. The logarithmic probability density distribution as image prior plays a crucial role in the MAP estimation. Hence, image prior modeling is an important task in iterative reconstruction methods. Generally, there is no explicit formula for an image prior model. Conventional iterative reconstruction methods usually utilize regularization techniques such as compressive sensing to represent the image prior model, which assumes that images can be converted into sparse, low-rank, or low-dimensional counterparts using transformation methods such as total variation, wavelet transform, or Fourier transform. However, these presumptions do not always accurately reflect structure and characteristics of actual images, especially for sophisticated biomedical images. Nevertheless, the emerging machine learning provides a powerful technique for the data modeling. Based on the score matching model [13, 16], the score function, i.e., the gradient of the logarithmic probability density distribution, can be established as the image prior through deep learning. The score function serves as the key element in this iterative process, contributing to image refinement in the each iteration. By leveraging the Gaussian mixture model to characterize noise distributions between a reconstructed image and the target image, we have developed a novel score matching formula to learn an advanced score function (ADSF) from a CT image dataset. Because the image reconstruction has a complex relationship between the reconstructed images and the true images, the Gaussian mixture model can describe a practical image noise distribution much better than the popular single Gaussian distribution. The new score function learns directly from the data distribution, and provides an advanced mechanism for effective and efficient extraction of data-driven prior information and a superior representation of underlying images for image reconstruction. In summary, we have proposed a novel score-matching formula to generate an advanced score function (ADSF) through deep learning from an image dataset. By incorporating the ADSF into the image reconstruction process, a new ADSF iterative reconstruction method has been developed in the MAP estimation framework to improve image reconstruction quality. The convergence of the ADSF iterative reconstruction algorithm has been proved through mathematical analysis. The performance of ADSF image reconstruction method has been evaluated on both public medical image datasets and clinical raw datasets. By comparing with the competing image reconstruction techniques such as the filtered backprojection (FBP), model-based iterative reconstruction (MBIR) with total variation minimization (MBIR-TV), and the state-of-the-art score function-based reconstruction (SFrecon) method, the quantitative results has shown that our proposed ADSF reconstruction method can achieve higher quality images in terms of PSNR and SSIM metric on diverse dataset. For Siemens and GE clinical CT raw datasets, the ADSF image reconstruction approach has achieved better denoising and deblurring effects than the competing methods, showing excellent generalizability and stability. The proposed new score matching formula can also be applied to image denoising, image deblurring, and image generation. Further algorithmic optimization and systematic evaluation are in progress to translate this new approach into clinical applications. Methods Image reconstruction approach In X-ray imaging, the number of X-ray photons recorded by a detector element is a random variable ξ, which can be modeled as the Poisson distribution [2]: (1) pξ=yi=y‾iyiyi!exp−y‾i, where y‾i is the expectation value of recorded X-ray photons along a path l from the X-ray source to the ith detector element, and obeys the Beer-Lambert law: (2) y‾i=niexp−∫l  μr→dl, where ni is the number of X-ray photons recorded by the ith detector element in the blank scan, without any object in the beam path, and μ(r) is the linear attenuation coefficient distribution within an object to be reconstructed. For the numerical implementation, Eq. (2) is discretized as, (3) y‾i=niexp−Aiμ, where μ is a vector of pixel values in the linear attenuation coefficient image, and Ai is weighting coefficients of the pixel values along the ith beam path. Assuming that measurements are independent, the likelihood function of the X-ray projection data can be obtained by (4) p(Y∣μ)=∏i=1m(y¯i)yiyi!exp(−y¯i), where Y=y1,y2,⋯,ymT is the number of photons measured by detectors, and m is the total number of X-ray detectors. Based on the Bayesian theorem: p(μ∣Y)p(Y)=p(Y∣μ)p(μ), the image reconstruction can be performed using the maximum a posteriori (MAP) estimation, which is equivalent to the following minimization problem [17]: (5) μmin=argminμ(∑i=1m[y¯i−yilog(y¯i)]−logp(μ)), where logp(μ) is the logarithmic probability density of an attenuation image, which expresses the prior knowledge about the underlying images. Combining Eqs. (3)–(5), we have (6) μmin=argminμ(∑i=1mniexp−Aiμ+yiAiμ−logpμ). Applying the second-order Taylor approximation, Eq. (6) can be simplified to the following quadratic optimization problem [4]: (7) μmin=argminμ12(Aμ−b)TD(Aμ−b)−logp(μ) where A is the m×n system matrix composed of the row vectors A1,A2,⋯,Am, and D is the diagonal matrix in the form diagy1,y2,⋯,ym. The optimization problem defined by Eq. (7) can be solved using the gradient-based method iteratively: (8) μk+1=μk−ωATDAμk−b−σ∇logpμk,k=1,2,⋯ where ∇logp(μ) is the score function, defined as the gradient of the logarithmic probability density distribution with respect to the current image, and ω and σ are parameters for trade-offs between data fidelity and image prior information. Convergence of score-function-based image reconstruction For the convergent analysis of the iteration scheme Eq. (8), we assume that A is a m×n system matrix of rank n. The probability density function p(μ) is assumed to be sufficiently smooth, and the Hessian matrix ∇2logp(μ) of the logarithmic probability density function has bounded eigenvalues, denoted by hi0 for any nonzero vector x∈Rn, the matrix ATDA is positive definite, denoting its smallest and largest eigenvalues as λmin and λmax respectively. From Eq. (11), it is easy to find that by choosing the parameters ω and σ to satisfy that 0<σ<ωrmin/C and σC/rmin<ω<(2−σC)/rmax, there exist a positive constant 0