## 2 Materials and Methods  
### CT Image Blind Restoration and Reconstruction  
A standard formulation of the CT image blind restoration and reconstruction is given in the following vector-matrix form:  
\[\mathbf{P}\ =\ \mathbf{A}\mathbf{H}\mathbf{\psi}\boldsymbol{\alpha}+\mathbf{n} \tag{1}\]  
where \(\mathbf{P}\in\mathbb{R}^{m\times 1}\) is the known noisy projection data; \(\mathbf{A}\in\mathbb{R}^{m\times m}\) is the known measurement matrix of the X-ray CT system; \(\mathbf{H}\in\mathbb{R}^{m\times m}\) is the unknown blurring matrix constructed by the point spread function (\(\mathbf{h}\in\mathbb{R}^{m\times 1}\) denotes the point spread function (PSF) in vector form, \(\mathbf{H}\) performs convolution with \(\mathbf{h}\)); \(\mathbf{\psi}\in\mathbb{R}^{m\times k}\) (\(k<<m\)) is a sparse representation transform matrix and is known. \(\mathbf{f}\) denotes the CT image. \(\mathbf{f}\in\mathbb{R}^{m\times 1}\) can be changed into another form \(\boldsymbol{\alpha}\in\mathbb{R}^{k\times 1}\) by the sparse representation transform matrix, \(\mathbf{f}\ =\ \mathbf{\psi}\boldsymbol{\alpha}\). \(\boldsymbol{\alpha}\) denotes the sparse representation vector because it can describe the image and include fewer data. \(\mathbf{n}\in\mathbb{R}^{m\times 1}\) is the additive noise vector, \(\mathbf{n}\sim N(0,\ \sum_{n})\), \(\sum_{n}\in\mathbb{R}^{m\times m}\) is the covariance matrix. For simplicity, we assume that \(\sum_{n}\ =\ \sigma_{n}^{2}\mathbf{I}\). The problem in (1) is designed to estimate \(\boldsymbol{\alpha}\) and \(\mathbf{H}\) when the projection data \(\mathbf{P}\) is as few as possible. Here, \(\boldsymbol{\alpha}\) is for the estimation of CT image. It doesn't matter what is \(\mathbf{H}\) itself. The purpose of estimating \(\mathbf{H}\) is to eliminate the degradation caused by PSF. As to estimate \(\boldsymbol{\alpha}\in\mathbb{R}^{k\times 1}\) with the assumption, the problem is to solve the following problem  
\[\hat{\boldsymbol{\alpha}}\ =\ \operatorname*{argmax}_{\boldsymbol{\alpha}} \log p(\boldsymbol{\alpha}|\mathbf{P}) \tag{2}\]  
By applying the Bayesian formula, \(p(x|y)\ =\ \frac{p(y|x)p(x)}{p(y)}\), (2) can be written as  
\[\hat{\boldsymbol{\alpha}}\ =\ \operatorname*{argmax}_{\boldsymbol{\alpha}} \{\log p(\mathbf{P}|\boldsymbol{\alpha})+\log p(\boldsymbol{\alpha})\} \tag{3}\]  
where \(p(\boldsymbol{\alpha})\) is the prior probability distribution of \(\boldsymbol{\alpha}\). \(p(\boldsymbol{\alpha})\) provides the prior statistical information of the sparse representation. The closer the estimated prior distribution is to the actual distribution, the more accurate the estimated value of \(\boldsymbol{\alpha}\) will be. To this end, the statistical prior of \(\boldsymbol{\alpha}\) was described by a Gaussian compound distribution as did in [22]. In this case, the probability density is given by  
\[p(\boldsymbol{\alpha})\ =\ \int\limits_{0}^{+\infty}\frac{1}{(2\pi)^{\frac{1}{2}} \big{\lvert}2^{2}\sum_{v}\big{\rvert}^{\frac{1}{2}}}\exp\left(-\frac{\boldsymbol {\alpha}^{T}\sum_{v}^{-1}\boldsymbol{\alpha}}{2z^{2}}\right)\phi_{z}(z)dz \tag{4}\]  
where the mixing variable \(\mathbf{z}\geq 0\) is a random vector, \(\mathbf{z}\in\mathbb{R}^{k\times 1}\). \(\phi_{z}(z)\) defines the probability density of the scale variable and results in the different distribution in different applications. The Gaussian compound distribution \(\boldsymbol{\alpha}\) is not easy to solve directly, and the hierarchical Bayesian strategy is applied to convert \(\boldsymbol{\alpha}\) into the synthesis of two variables with the following form. \[\boldsymbol{\alpha}\ =\ \mathbf{z}\odot\mathbf{V} \tag{5}\]  
where \(\odot\) defines the Hadamard product. \(\mathbf{V}\sim N(0,\ \sum_{v})\) denotes a Gaussian random vector, \(\mathbf{V}\in\mathbb{R}^{k\times 1}\). \(\sum_{v}\in\mathbb{R}^{k\times k}\) is the covariance matrix of \(\mathbf{V}\). For the sake of simplicity, we assume that \(\mathbf{z}\) and \(\mathbf{V}\) are independent.