Local Linear Constraint based Optimization Model for Dual Spectral CT

TL;DR

This paper introduces a novel optimization model for dual spectral CT that leverages local linear constraints to improve image quality and noise robustness, validated through simulations and real data.

Contribution

It proposes a new local linear constraint-based optimization model for DSCT, enhancing image quality and noise suppression over traditional methods.

Findings

Improved signal-to-noise ratio in decomposed images

Effective noise reduction demonstrated in experiments

Enhanced material distinguishability in dual spectral CT

Abstract

Dual spectral computed tomography (DSCT) can achieve energy- and material-selective images, and has a superior distinguishability of some materials than conventional single spectral computed tomography (SSCT). However, the decomposition process is illposed, which is sensitive with noise, thus the quality of decomposed images are usually degraded, especially the signal-to-noise ratio (SNR) is much lower than single spectra based directly reconstructions. In this work, we first establish a local linear relationship between dual spectra based decomposed results and single spectra based directly reconstructed images. Then, based on this constraint, we propose an optimization model for DSCT and develop a guided image filtering based iterative solution method. Both simulated and real experiments are provided to validate the effectiveness of the proposed approach.

Equations

$$P_{k,l}=-\ln\int_{E}S_{k}(E)\exp{\big{(}-\mathcal{P}_{l}(\mu(E,\mathbf{x}))\big{)}}\,\mathrm{d}E,l\in\mathcal{L},k=1,2,$$

$$\mu (E,\mathbf{x}))=\sum _{i=1}^2{\psi }_i(E)f_i(\mathbf{x}),$$

$$P_{k,l}=-\ln\sum_{j=1}^{J_{k}}S_{k,j}\exp{\big{(}-\sum_{i=1}^{2}\psi_{i,j}\mathcal{P}_{l}(\mathbf{f}_{i})\big{)}}\Delta E,$$

$$P_{k,l}\approx P_{k,l}(n)+\left(\begin{array}{c}\frac{{\Psi }_{k,l}^1(n)}{Q_{k,l}(n)},\frac{{\Psi }_{k,l}^2(n)}{Q_{k,l}(n)}\end{array}\right)\left(\begin{array}{c}\mathcal{P}({\mathbf{f}}_1-{\mathbf{f}}_1(n))\\ \mathcal{P}({\mathbf{f}}_2-{\mathbf{f}}_2(n))\end{array}\right),$$

$$P_{k,l}(n)$$

$$Q_{k,l}(n)$$

$${\Psi }_{k,l}^1(n)$$

$${\Psi }_{k,l}^2(n)$$

$$\left(\begin{array}{c}\mathcal{P}({\mathbf{f}}_1(n+1))\\ \mathcal{P}({\mathbf{f}}_2(n+1))\end{array}\right)=\left(\begin{array}{c}\mathcal{P}({\mathbf{f}}_1(n))\\ \mathcal{P}({\mathbf{f}}_2(n))\end{array}\right)+\frac{C_l(n)}{\det (M_l(n))}\left(\begin{array}{c}{P}_{1,l}-P_{1,l}(n)\\ P_{2,l}-P_{2,l}(n)\end{array}\right),$$

$$M_l(n)=\left(\begin{array}{c}\frac{{\Psi }_{1,l}^1(n)}{Q_{1,l}(n)},\frac{{\Psi }_{1,l}^2(n)}{Q_{1,l}(n)}\\ \frac{{\Psi }_{2,l}^1(n)}{Q_{2,l}(n)},\frac{{\Psi }_{2,l}^2(n)}{Q_{2,l}(n)}\end{array}\right),C_l(n)=\left(\begin{array}{c}\frac{{\Psi }_{2,l}^2(n)}{Q_{2,l}(n)},-\frac{{\Psi }_{1,l}^2(n)}{Q_{1,l}(n)}\\ -\frac{{\Psi }_{2,l}^1(n)}{Q_{2,l}(n)},\frac{{\Psi }_{1,l}^1(n)}{Q_{1,l}(n)}\end{array}\right).$$

$$y_i=a_k{I}_i+b_k,\forall i\in {\omega }_k,$$

$$y_i=x_i-t_i.$$

$$\min\sum_{i\in\omega_{k}}\big{(}(a_{k}I_{i}+b_{k}-x_{i})^{2}+\epsilon a_{k}^{2}\big{)},$$

$$a_k$$

$$b_k$$

$$y_i={\bar{a}}_i{I}_i+{\bar{b}}_i.$$

$${\mathbf{f}}^{(j)}=a^{(k)}{\mathbf{g}}^{(j)}+b^{(k)},\forall j\in {\omega }^{(k)},$$

$$\displaystyle\min_{(\mathbf{f}_{1},\mathbf{f}_{2})}\bigg{\{}\Big{\|}\begin{pmatrix}\mathcal{P}(\mathbf{f}_{1})-\mathcal{P}(\mathbf{f}_{1}(n))\\ \mathcal{P}(\mathbf{f}_{2})-\mathcal{P}(\mathbf{f}_{2}(n))\end{pmatrix}-\frac{C_{l}(n)}{\det(M_{l}(n))}\begin{pmatrix}P_{1,l}-P_{1,l}(n)\\ P_{2,l}-P_{2,l}(n)\end{pmatrix}\Big{\|}_{L_{2}}^{2}$$

$$\displaystyle+\sum_{i=1}^{2}\sum\limits_{k}\sum\limits_{j\in\omega_{i}^{(k)}}\frac{\xi_{i}}{|\omega_{i}|}\Big{(}\big{(}a_{i}^{(k)}\mathbf{g}_{\mathbf{f}_{i}}^{(j)}+b_{i}^{(k)}-\mathbf{f}_{i}^{(j)}(n)\big{)}^{2}+\epsilon_{i}(a_{i}^{(k)})^{2}\Big{)}\bigg{\}},$$

$$\displaystyle\big{(}\mathcal{P}(\mathbf{f}_{1}(*)),\mathcal{P}(\mathbf{f}_{2}(*))\big{)}=$$

$$\displaystyle\min_{\big{(}\mathcal{P}(\mathbf{f}_{1}),\mathcal{P}(\mathbf{f}_{2})\big{)}}\Big{\|}\begin{pmatrix}\mathcal{P}(\mathbf{f}_{1})-\mathcal{P}(\mathbf{f}_{1}(n))\\ \mathcal{P}(\mathbf{f}_{2})-\mathcal{P}(\mathbf{f}_{2}(n))\end{pmatrix}-\frac{C_{l}(n)}{\det(M_{l}(n))}\begin{pmatrix}P_{1,l}-P_{1,l}(n)\\ P_{2,l}-P_{2,l}(n)\end{pmatrix}\Big{\|}_{L_{2}}^{2},$$

$$\displaystyle\big{(}\mathbf{f}_{1}(*),\mathbf{f}_{2}(*)\big{)}=\mathcal{P}^{-1}\big{(}\mathcal{P}(\mathbf{f}_{1}(*)),\mathcal{P}(\mathbf{f}_{2}(*))\big{)},$$

$$\displaystyle\big{(}a_{i}(n+1)),b_{i}(n+1))\big{)}=$$

$$\displaystyle\min_{(a_{i},b_{i})}\sum\limits_{k}\sum\limits_{j\in\omega_{i}^{(k)}}\frac{\xi_{i}}{|\omega_{i}|}\Big{(}\big{(}a_{i}^{(k)}\mathbf{g}_{\mathbf{f}_{i}}^{(j)}+b_{i}^{(k)}-\mathbf{f}_{i}^{(j)}(*)\big{)}^{2}+\epsilon_{i}(a_{i}^{(k)})^{2}\Big{)},$$

$${\bar{a}}_i^{(k)}=\frac{1}{|{\omega }_i|}\sum _{j\in {\omega }_i^{(k)}}{a}_i^{(j)},$$

$${\bar{b}}_i^{(k)}=\frac{1}{|{\omega }_i|}\sum _{j\in {\omega }_i^{(k)}}{b}_i^{(j)},$$

$${\mathbf{f}}_i(n+1)={\bar{a}}_i{\mathbf{g}}_{{\mathbf{f}}_i}+{\bar{b}}_i,$$

Full text

1 Introduction

In X-ray dual spectral computed tomography (DSCT), a specimen is scanned with two different X-ray energy spectra, then the collected polychromatic projections from this procedure are utilized to perform energy- and material-selective reconstruction [1, 2, 3, 4, 5, 6]. Compared with conventional single spectral computed tomography (SSCT), DSCT has many advantages for contrast enhancement, artifact reduction, and material composition assessment. Thus, DSCT has wide potential applications in medical and industrial domains, such as bone mineral density and liver iron concentrations measurements, beam-hardening correction and contrast enhancement of soft tissue, positron emission tomography (PET) attenuation correction, rock core characterization for petrochemical industry, and so forth [7, 8, 9, 10, 11, 12, 13, 14].

Existing methods to perform the decomposition of DSCT can be classified into three groups: image based methods, projection based methods, and iterative methods. Image based methods treat the projection data sets as being independent until they are reconstructed. Then, images of each spectra are combined linearly to obtain two decomposed images [15]. Such methods fail to describe the real nonlinearity relationship between decomposed results and polychromatic projections, thus beam hardening artifacts cannot be removed successfully [16, 7]. Projection based methods correctly treat the available information by passing the projection data through a high order decomposition function, followed by image reconstruction [17, 18]. Generally, they can obtain better decomposition performances than image based ones, but the combination of polychromatic projections before reconstruction requires a satisfication of geometrical consistency. Several iterative methods are proposed based on statistical model and nonlinear optimization [19, 20, 21, 22]. By introducing prior knowledge or establishing nonlinear model, these methods improve the quality of decomposed results effectively, but they usually have a drawback of slow convergence and huge computational cost. An extended algebraic reconstruction technique (E-ART) for DSCT was proposed by Zhao et al. [23]. It describes the DSCT reconstruction as a nonlinear system problem, and then extends the classic ART method to solve the model iteratively. This method can produce high quality decompositions from polychromatic projections, but it also has the drawback of slow convergence. Hu et al. developed the E-ART method into an extended simultaneous algebraic reconstruction version, i.e. E-SART [24]. This method is based on the matrix inversion and has a high degree of parallelism, thus the convergence rate is improved dramatically. But the illposedness of the decomposition process render it noticeably sensitive with noise, thus the achieved quality is degraded with reduced signal-to-noise ratio (SNR).

Comparing decomposed results of DSCT with reconstructed images of SSCT, we find that although single spectra based reconstructions have weaker capability of destinguishing some materials, the achieved quality like SNR is dramatically higher. Moreover, there is an interesting relationship between them, i.e., decomposed results of DSCT can be considered as modifications of reconstructed images of SSCT by removing some components and adjusting gray values. Further, this structure-based feature can be described as a local linear relationship mathematically. By combing it as a constraint into an optimization model for DSCT, the reconstructed image of SSCT works as a guided image and contributes to improving the smoothness of decomposed results of DSCT effectively. Thus, the systematic noise are well surpressed and the quality of decomposed results are improved considerably.

In this work, a local linear model is first established to describe the relationship between decomposed results of DSCT and reconstructed images of SSCT. Then by employing it as a constraint, we propose an optimization model for DSCT. Based on the guided image filtering [25, 26], a correlative iterative solution method is developed, which produces decomposed images with high SNR and maintains a fast convergence meanwhile.

The remainder of this paper is organized as follows. In section 2, the mathematical model of DSCT is presented, then the E-SART method and the guided image filtering technique are reviewed briefly. In section 3, we first propose a local linear constraint based optimization model DSCT. Then we develop a guided image filtering based iterative solution method. In section 4, both simulated and real experiments are provided to verify the effectiveness of the proposed methods. Section 5 contains our conclusions.

2 Methods

2.1 Mathematical Model of DSCT

By considering the X-ray is polychromatic and assuming collected raw data are consistent geometrically, we can describe the physical process of DSCT as follows,

[TABLE]

where $\mu(E,\mathbf{x})$ is the linear attenuation coefficient of a specimen at a spatial position $\mathbf{x}$ and energy $E$; $\mathcal{P}_{l}(\cdot)$ represents the ray transform which is an integral transform along a ray path $l$; $S_{k}(E)$ is the $k$-th normalized emission spectrum; $P_{k,l}$ indicates the acquired information, frequently called projection data.

In DECT, the linear attenuation coefficient $\mu(E,\mathbf{x})$ is usually considered splittable with variable $E$ and $\mathbf{x}$, i.e.,

[TABLE]

where $\psi_{i}(E)$ is a function of energy; $f_{i}(\mathbf{x})$ is a function of spatial position. There are commonly two physical explanations of eq. (2): basis material based decomposition and effect based decomposition. For the former, $\psi_{i}(E)$ is the mass attenuation coefficient for material $i$, and $f_{i}(\mathbf{x})$ represents the correlative density distribution. For the latter, $\psi_{1}(E)=E^{-3}$ and $\psi_{2}(E)=KN(E)$ (Klein-Nishina function) corresponds to the photoelectric effect and Compton scattering respectively, and $f_{i}(\mathbf{x})$ represents the effect distribution correspondingly. The aim of DSCT is to reconstruct images of distrubution functions $f_{i}(\mathbf{x}),i=1,2$.

2.2 E-SART Method

By substituting eq. (2) into eq. (1) and discretizing the correlative result according to energy bin, we get

[TABLE]

where $J_{k}$ is the enery bin number of spectra $k$; $\Delta E$ represents the bin length; $S_{k,j}$ and $\psi_{i,j}$ are the samplings of $S_{k}(E)$ and $\psi_{i}(E)$ within bin $j$; $\mathbf{f}_{i}$ is a one demensional column vector representing the discretized distribution function. Then, the 1st order Taylor expansion of eq. (3) at point $(\mathbf{f}_{1}(n);\mathbf{f}_{2}(n))$ is

[TABLE]

where

[TABLE]

Along each ray path, two projection data are acquired based on different X-ray spectra, i.e., low-energy and high-energy. By solving the system of linear equation (4), $k=1,2$, we can get the projection of distribution fuction in a iteration form,

[TABLE]

where

[TABLE]

Then by using the conventional Simultaneous Algebraic Reconstruction Technique (SART), distribution fuction $\mathbf{f}_{1}$ and $\mathbf{f}_{2}$ are updated iteratively.

Comparing with E-ART, E-SART improves the convergence rate dramatically. However, the illposedness of the inverse problem render this matrix inversion based decomposition process sensitive with inevitable systematic noise, i.e., obtained results suffer from low SNR. Thus, some prior knowledge or constraints are needed to improve the robust against noise.

2.3 Guided Image Filtering

Guided filter is an edge-preserving filter with a great variety of applications [27, 28, 29, 30], of which the key assumption is a local linear model between the guidance image $I$ and the filtering output $y$, i.e.,

[TABLE]

where $\omega_{k}$ is a window centered at the pixel $k$ with radius $r$; $(a_{k},b_{k})$ are some linear coefficients constant in $\omega_{k}$. Modeling the output $y$ as the input $x$ removing some unwanted noise or textures $t$:

[TABLE]

Thus, by minimizing the difference between $y$ and $x$ within a window $\omega_{k}$ while maintaining the linear model (6), the correlative optimization model is established as follows,

[TABLE]

where $\epsilon$ is a regularization parameter penalizing large $a_{k}$. The solution is given by

[TABLE]

where $\nu_{k}$ and $\sigma_{k}^{2}$ are the mean and variance of $I$ in $\omega_{k}$; $|\omega|$ is the number of pixels in $\omega_{k}$; $\bar{x}_{k}$ is the mean of $x$ in $\omega_{k}$. Then, the filtering output $y$ can be computed according to eq. (6). Because a pixel $i$ is involved in all the covered windows, by averagig all the possibles output values, we get

[TABLE]

Here $\bar{a}_{i}$ and $\bar{b}_{i}$ are the average coefficients of all windows overlapping $i$. By using the guided image filtering, input $x$ is refined by the guided image $I$ according to the local linear relationship between them.

3 Local Linear Constraint based Optimization Model and Iterative Solution Method

3.1 Local Linear Constraint based Optimization Model

Although directly reconstructed images under single spectra scan have weak capability of distinguishing some materials, the quality is improved remarkbaly than dual spectra based decomposition results, especially when the noise level is relatively high. Moreover, there are some structure-based relationships between them. An intuitionistic character is that decomposed results of DSCT can be viewed as modifications of reconstructed images of SSCT by removing some components and adjusting gray values. When analyzing this feature in detail, as is illustrated in Fig. 1, we can find that a linear relationship is usually held in small patches, of which the discreted version reads,

[TABLE]

where $\mathbf{f}$ represnets a decomposed result of DSCT; $\mathbf{g}$ represents a reconstructed image of SSCT; $j$ is a pixel index; $\omega^{(k)}$ is a window centered at the pixel $k$ with radius $r$; $(a^{(k)},b^{(k)})$ are some linear coefficients constant in $\omega^{(k)}$.

Based on this constraint, we proposed an optimization model for DSCT as follows,

[TABLE]

where $\mathbf{g}_{\mathbf{f}_{i}}$ is a selected result of SSCT corresponding to $\mathbf{f}_{i}$; $\xi_{i}$ and $\epsilon_{i}$ are regularization parameters.

In model (3.1), for each searched-for decomposition result $\mathbf{f}_{i},i=1,2$, we employ a correlative local linear constraint. Thus, the smoothness knowledge of reconstructed image of SSCT is effectively introduced into the decomposition problem. By weakening the illposedness, the noise is surpressed noticeably, thus the quality of decomposed result, like SNR, is improved dramatically.

3.2 Iterative Solution Method

When solving the proposed optimization model, considering the data term is measured in projection domain and the regularity terms are measured in image domain, we split model (3.1) into two sub-optimization problems and develop an iterative scheme as follows,

[TABLE]

where Eqs. (8)-(8f) are implemented respectively for $i=1,2$. We use the E-SART method reviewed in 2.2 to solve Eqs. (8) and (8b). When solving Eqs. (8)-(8f), the guided image filtering reviewed in 2.3 is employed.

For one thing, the proposed method suppresses the noise effectively, thus the quality of decomposed results is improved dramatically. For another, it maintains the merit of high convergence rate of the E-SART meanwhile. By considering the relatively weaker nonlinearity of high-energy spectra based data than low-energy spectra based ones, in this study we just use the high-energy spectra based reconstructed image as the guided image for both decomposition results $f_{i},i=1,2$. More analyses and comparisons are still needed to choose a more satisfactory guided image, which will be fully studied in our continuous work.

4 Experiments

In order to verify the effectiveness of the proposed method, experiments with both simulated and real data are performed. In the simulated experiments, both noise-free case and noisy case are tested. These experiments are restricted to fan beam CT for simplicity. The generalization to cone beam case is straightforward. As a comparison, we implement the E-SART method as well, because it is a relatively new method with fast convergence rate. We employ all the methods for the case of basis material based decomposition, which can be extended to the case of effect based case easily. And the iteration number is fixed to 30.

4.1 Simulated Experiments

The phantom used in all the simulated experiments is the 2D FORBILD head phantom without ears shown in Fig. 1 [31]. Water and bone are chosen as the two basis materials, of which the mass attenuation coefficients are retrieved from the National Institute of Standard Technology (NIST) tables of X-ray mass attenuation coefficient [32]. The polychromatic spectra of a GE Maxiray 125 X-ray tube are simulated by using an open source X-ray spectra simulator, SpectrumGUI [33]. Two tube voltages, 80 kV and 140 kV, are chosen, where the latter is filtered with 1mm copper. The correlative spectra are shown in Fig. 2. The energy of photons emitted from the source is 8 MeV. The detector consists of 512 channels with length 0.03 cm. The source-object distance (SOD) is 100 cm and the source-detector distance (SDD) is 120 cm. With this configuration, images are reconstructed to a 512 $\times$ 512 digitization with a pixel size of 0.0249 cm.

Under this setting, we test both the noise-free case and the noisy case, seen Fig. 3 and Fig. 4 respectively. All the experimenal resutls demonstrate the conspicuous noise tolerance of the proposed method.

4.2 Real Experiments

In the real experiment, the measured specimen is a bone submerged in water. An X-ray source (YXLON Y.TU450 D09 tube) is operated at the tube voltage of 80 kV and 140 kV for low- and high-energy spectra scan respectively, and the tube current is 5 mA. The employed flat-pannel detector (YXLON Y.LDA detector) has $1920\times 1920$ channels with mesh size of 0.0127 cm. The SOD is 23.15 cm and the SDD is 69.67 cm. By using a collimator, the data from the central slice are obtained to validate the proposed method. The iteration number is 10.

The decomposed results by using the E-SART method and the proposed method are shown in Fig. 5. It is noticeable that the proposed method can surpress the noice effectively and improve the smoothness dramatically.

5 Conclusion

In this work, we first establish a local linear constraint to describe the structure relationship between dual spectra based decomposed results and single spectra based directly reconstructed images. Then the correlative optimization model and the iterative solution method are proposed respectively. By employing the guided image filtering, the smoothness knowledge of the single spectra based reconstructed image is effectively introduced into the decomposition process. This method reduces the illposedness of the DSCT, thus the noise in the decomposition process is surpressed successfully. Both simulated and real experiments demonstrate the effectiveness of the proposed method. Further studies are needed to select a satisfactory guided image.