Joint Reconstruction in Low Dose Multi-Energy CT

TL;DR

This paper introduces a joint regularization approach combined with a low-dose, non-overlapping projection protocol to improve multi-energy CT reconstruction, reducing radiation dose while maintaining image quality.

Contribution

It proposes a novel joint regularization model based on the structural similarity index and demonstrates its effectiveness with a new low-dose data acquisition protocol.

Findings

Joint regularization outperforms individual reconstructions.

The combination reduces radiation dose significantly.

The new model improves multi-energy CT image quality.

Abstract

Multi-energy CT takes advantage of the non-linearly varying attenuation properties of elemental media with respect to energy, enabling more precise material identification than single-energy CT. The increased precision comes with the cost of a higher radiation dose. A straightforward way to lower the dose is to reduce the number of projections per energy, but this makes tomographic reconstruction more ill-posed. In this paper, we propose how this problem can be overcome with a combination of a regularization method that promotes structural similarity between images at different energies and a suitably selected low-dose data acquisition protocol using non-overlapping projections. The performance of various joint regularization models is assessed with both simulated and experimental data, using the novel low-dose data acquisition protocol. Three of the models are well-established, namely the joint total variation, the linear parallel level sets and the spectral smoothness promoting regularization models. Furthermore, one new joint regularization model is introduced for multi-energy CT: a regularization based on the structure function from the structural similarity index. The findings show that joint regularization outperforms individual channel-by-channel reconstruction. Furthermore, the proposed combination of joint reconstruction and non-overlapping projection geometry enables significant reduction of radiation dose.

Tables

Table 1. Table 1: Imaging geometry used for collecting the experimental data.

Parameter	Value
Focus-center distance	252 mm
Focus-detector distance	420 mm
Magnification	5/2
Detector pixel size	0.200 mm
Effective pixel size	0.120 mm
Projection size	552 $\times$ 576 pixels
Angular range	$360°$
#projections	720

Table 2. Table 2: Energy-specific settings used for collecting the experimental data.

Energy	$U$ (kV)	Filtration	$I$ ($\mu A$)	Exposure time (ms)	Frame averaging
$E_1$	50	None	300	125	4
$E_2$	80	1 mm Al	180	125	4
$E_3$	120	0.5 mm Cu	120	250	4

Table 3. Table 3: Selected α 𝛼 \alpha and γ k subscript 𝛾 𝑘 \gamma_{k} parameters for the different regularization terms for the simulated data test case ( α SIM , γ k , SIM subscript 𝛼 SIM subscript 𝛾 𝑘 SIM \alpha_{\textrm{SIM}},\gamma_{k,\textrm{SIM}} ) and for the experimental data test case ( α EXP , γ k , EXP subscript 𝛼 EXP subscript 𝛾 𝑘 EXP \alpha_{\textrm{EXP}},\gamma_{k,\textrm{EXP}} ). The three values of γ k subscript 𝛾 𝑘 \gamma_{k} are for the three different tube energies with indices k = 1 , 2 , 3 𝑘 1 2 3 k=1,2,3 .

Prior	Acronym	${\alpha }_{\text{SIM}}$	${\gamma }_{k,\text{SIM}}$	${\alpha }_{\text{EXP}}$	${\gamma }_{k,\text{EXP}}$
No prior	No prior	-	-	-	-
Total variation	TV	-	2.25, 0.8, 0.7	-	1.0, 0.5, 0.25
Joint total variation	JTV	2.25	0	1	0
Linear parallel level sets	LPLS	4000	0	6000	0
First difference	D1	400	0	100	0
Structural	S	4e6	0	1e6	0
First difference + TV	D1+TV	50	1.125, 0.4, 0.35	70	0.5, 0.25, 0.125
Structural + TV	S+TV	0.75e5	1.125, 0.4, 0.35	0.5e6	0.5, 0.25, 0.125

Table 4. Table 4: Measurement angles used for energies E 1 subscript 𝐸 1 E_{1} , E 2 subscript 𝐸 2 E_{2} and E 3 subscript 𝐸 3 E_{3} in the simulated data case (left) and experimental data case (right).

$E_1$	$E_2$	$E_3$		$E_1$	$E_2$	$E_3$
0	2	4		0	4	8
6	8	10		12	16	20
12	14	16		24	28	32
$⋮$	$⋮$	$⋮$		$⋮$	$⋮$	$⋮$
162	164	166		324	328	332
168	170	172		336	340	344
174	176	178		348	352	356

Equations

$$I=I_0(E_k)e^{-{\int }_{\text{ray}}\mu (r,E_k)ds},$$

$${\int }_{\text{ray}}\mu (r,E_k)ds=-\ln \frac{I}{I_0(E_k)}.$$

$$I=\underset{0}{\overset{E_{\text{max}}}{\int }}{I}_0(E)e^{-{\int }_{\text{ray}}\mu (r,E)ds}dE.$$

$$y_k=A_k{x}_k+e_k$$

$$x_k=\arg \underset{x_k\ge 0}{\min }\left\{||y_k-A_k{x}_k|{|}^2\right\}\text{(No prior)}$$

$$x_k=\arg \underset{x_k\ge 0}{\min }\left\{||y_k-A_k{x}_k|{|}^2+{\gamma }_{k}TV(x_k)\right\}\text{(TV)}$$

$$TV(x_k)={\int }_{\Omega }{\left(||\nabla x_k|{|}^2+{\beta }^2\right)}^{1/2}\mathrm{d}r$$

$$(x_1,\dots ,x_n)=\arg \underset{x_1\ge 0,\dots ,x_n\ge 0}{\min }\left\{\sum _{k=1}^n||y_k-A_k{x}_k|{|}^2+R(x_1,x_2,\dots ,x_n)\right\}$$

$$R(x_1,x_2,\dots ,x_n)=\sum _{k=1}^n{\gamma }_{k}U(x_k)+\alpha W(x_1,x_2,\dots ,x_n)$$

$$JTV(x_1,x_2,\dots ,x_n)={\int }_{\Omega }{\left(||\nabla x_1|{|}^2+||\nabla x_2|{|}^2+\dots +||\nabla x_n|{|}^2+{\beta }^2\right)}^{1/2}\mathrm{d}r$$

$$R(x_1,x_2,\dots ,x_n)=\alpha JTV(x_1,x_2,\dots ,x_n)\text{(JTV)}$$

$$LPLS(x_k,x_j)={\int }_{\Omega }({\left(||\nabla x_k|{|}^2+{\beta }^2\right)}^{1/2}{\left(||\nabla x_j|{|}^2+{\beta }^2\right)}^{1/2}$$

$$-{\left(|⟨\nabla x_k,\nabla x_j⟩{|}^2+{\beta }^4\right)}^{1/2})\mathrm{d}r$$

$$LPLS(x_1,x_2,\dots ,x_n)=LPLS(x_1,x_2)+LPLS(x_2,x_3)+$$

$$\dots +LPLS(x_n,x_1).$$

$$R(x_1,x_2,\dots ,x_n)=\alpha LPLS(x_1,x_2,\dots ,x_n)\text{(LPLS)}$$

$$D1(x_k,x_{k+1})=w_k{\int }_{\Omega }||x_{k+1}-x_k|{|}^2\mathrm{d}r$$

$$D1(x_1,x_2,\dots ,x_n)=D1(x_1,x_2)+D1(x_2,x_3)+\dots +D1(x_{n-1},x_n).$$

$$R(x_1,x_2,\dots ,x_n)=\alpha D1(x_1,x_2,\dots ,x_n)\text{(D1)}$$

$$R(x_1,x_2,\dots ,x_n)=\sum _{k=1}^n{\gamma }_{k}TV(x_k)+\alpha D1(x_1,x_2,\dots ,x_n)\text{(D1+TV)}$$

$$\hat{S}(x_1,x_2)=\frac{{\sigma }_{x_1{x}_2}+C}{{\sigma }_{x_1}{\sigma }_{x_2}+C}$$

$$({\sigma }_{x_k}{)}_i=$$

$$({\mu }_{x_k}{)}_i=$$

$$({\sigma }_{x_k{x}_p}{)}_i=$$

$$\bar{S}(x_k,x_p)=\frac{1}{N}\sum _{i=1}^N(\hat{S}(x_k,x_p){)}_i$$

$$S(x_1,x_2,\dots ,x_n)=\frac{1}{\bar{S}(x_1,x_2,\dots ,x_n)}$$

$$\bar{S}(x_1,x_2,\dots ,x_n)=\bar{S}(x_1,x_2)+\bar{S}(x_2,x_3)+\dots +\bar{S}(x_n,x_1)$$

$$R(x_1,x_2,\dots ,x_n)=\alpha S(x_1,x_2,\dots ,x_n)\text{(S)}$$

$$R(x_1,x_2,\dots ,x_n)=\sum _{k=1}^n{\gamma }_{k}TV(x_k)+\alpha S(x_1,x_2,\dots ,x_n)\text{(S+TV)}$$

$$\text{RMSE}(x_k)=\sqrt{\frac{1}{N}\sum _{i=1}^N((x_k{)}_i-(x_{\text{REF,k}}{)}_i{)}^2}$$

$$=l(x_1,x_2)c(x_1,x_2)s(x_1,x_2)$$

$$=\left(\frac{2{\mu }_{x_1}{\mu }_{x_2}+C_1}{{\mu }_{x_1}^2+{\mu }_{x_2}^2+C_1}\right)\left(\frac{2{\sigma }_{x_1}{\sigma }_{x_2}+C_2}{{\sigma }_{x_1}^2+{\sigma }_{x_2}^2+C_2}\right)\left(\frac{{\sigma }_{x_1{x}_2}+C_3}{{\sigma }_{x_1}{\sigma }_{x_2}+C_3}\right)$$

$$\text{MSSIM}(x_k,x_{\text{REF},k})=\frac{1}{N}\sum _{i=1}^N(\text{SSIM}(x_k,x_{\text{REF},k}){)}_i$$

Full text

Joint Reconstruction in Low Dose

Multi-Energy CT

Jussi Toivanen

Department of Applied Physics, University of Eastern Finland, POB 1627, FI-70211 Kuopio, Finland

Alexander Meaney

Department of Mathematics and Statistics, University of Helsinki, Finland

Samuli Siltanen

Department of Mathematics and Statistics, University of Helsinki, Finland

Ville Kolehmainen

Department of Applied Physics, University of Eastern Finland, POB 1627, FI-70211 Kuopio, Finland

Abstract

Multi-energy CT takes advantage of the non-linearly varying attenuation properties of elemental media with respect to energy, enabling more precise material identification than single-energy CT. The increased precision comes with the cost of a higher radiation dose. A straightforward way to lower the dose is to reduce the number of projections per energy, but this makes tomographic reconstruction more ill-posed. In this paper, we propose how this problem can be overcome with a combination of a regularization method that promotes structural similarity between images at different energies and a suitably selected low-dose data acquisition protocol using non-overlapping projections. The performance of various joint regularization models is assessed with both simulated and experimental data, using the novel low-dose data acquisition protocol. Three of the models are well-established, namely the joint total variation, the linear parallel level sets and the spectral smoothness promoting regularization models. Furthermore, one new joint regularization model is introduced for multi-energy CT: a regularization based on the structure function from the structural similarity index. The findings show that joint regularization outperforms individual channel-by-channel reconstruction. Furthermore, the proposed combination of joint reconstruction and non-overlapping projection geometry enables significant reduction of radiation dose.

1 Introduction

Computed tomography (CT) is an imaging modality in which the interior structure of an object is reconstructed from X-ray transmission images recorded along different directions. The reconstructed object is typically presented as cross-sectional images or 3D volumetric data. CT imaging plays an ever-diversifying role in medical imaging and industrial applications, with continuing methodological developments being motivated by the desire for increased image quality and information, and for reductions in radiation dose [1, 2]. Advances in CT imaging have been driven both by hardware improvements, such as better detectors and X-ray tubes, and by algorithmic developments in reconstruction methods [3, 4].

A CT image is a greyscale image where the value of a pixel or voxel in a given location is proportional to the X-ray attenuation coefficient of the medium in the corresponding location. A significant limitation of conventional CT imaging is that no bijective relationship exists between the reconstructed greyscale value and the material properties, i.e., elemental composition and mass density. Therefore, two tissues of entirely differing properties may be reconstructed as the same value, a typical example being bone and iodine contrast agent [5]. More information on the material composition can be obtained by using dual- or multi-energy CT imaging, where the object is imaged using more than one X-ray energy, and the non-linearly varying attenuation properties of different elemental media are exploited for material identification [1, 5]. Applications of multi-energy CT, also known as spectral CT, include differentiation and quantification of materials, tissue characterization, virtual monoenergetic imaging, automated bone removal in CT angiography, cardiovascular imaging, multiple contrast agent imaging, and mapping of effective atomic number [5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17].

Multi-energy CT imaging was first proposed in the 1970s [18], and today there exist many data acquisition schemes for spectral X-ray imaging: sequential scanning, rapid tube voltage switching, dual source scanning, multilayer detectors, and photon-counting detectors [5]. Compared to traditional radiography, CT is a high-dose modality [19]. Although the risks of radiation doses in the range of a few millisieverts remains a controversial topic [20, 21], an underlying philosophy in CT development is to reduce the dose as much as possible, in accordance with the ALARA (As Low As Reasonably Achievable) principle [22, 23]. There is a concerted effort to push the radiation dose from routine CT examinations to less than 1 millisievert, which is well below the annual dose due to background radiation [3]. Both image quality and radiation dose depend fundamentally on the number of photons used in imaging, and an acceptable balance needs to be found between the two factors [24]. Radiation dose can be reduced by reducing the number of photons used, and one efficient way of achieving this is to reduce the number of X-ray projections [3, 25, 26, 27]. However, image reconstruction from sparse projections requires algorithms more advanced than the standard approach where attenuation at each measured energy is reconstructed independently using classical techniques, such as the filtered backprojection (FBP) algorithm.

A highly feasible approach for low dose multi-energy CT is to utilize a multi-channel joint reconstruction approach where all the unknown images are reconstructed simultaneously by solving one combined inverse problem. The central idea is to combine all the projection data into a single image reconstruction problem and to utilize regularization models that promote some prior information on the unknown images within and across the energies. In multi-energy CT, a reasonable prior assumption for the attenuation images at different energies is that they can be expected to be structurally similar in the sense that an edge (e.g. an organ boundary) that is present at one energy, is likely to be at same location and alignment with the other energies as well, even though the contrast between materials will be different at each energy. There are many models in the literature designed for promoting structural similarity of images. These models can work either in a one-sided way, where the model is used to promote structural similarity of the unknown image with a fixed reference image (see e.g. [28, 29]), or two-sided way where the model promotes structural similarity between two (or more) images in a joint reconstruction formulation. These include vectorial total variation [30, 31], spectral patch-based penalty for the maximum likelihood method [32], tensor-based dictionary learning [33, 34], parallel level sets [35], the prior rank, intensity and sparsity model (PRISM) [36, 37], tensor-based nuclear norm regularization [38], nonlocal low-rank and sparse matrix decomposition [39], and total variation regularization using non-convex optimization [40].

In this paper, we study the performance of various structural similarity promoting joint reconstruction models using a novel low dose data acquisition protocol. We consider three established joint reconstruction models: joint total variation which was originally developed for denoising of multichannel images [41], linear parallel level sets [42] and a smoothness regularization similar to the models that have been employed to dynamic MRI in [43] and fluoresence x-ray tomography in [44]. Furthermore, we introduce a new joint regularization model based on the structural similarity index [45].

We combine the joint reconstruction approach with a low-dose data acquisition protocol. We collect data using a combination of a sparse angular sampling and a method of using a different X-ray energy in consecutive steps similar to [46]. See Figure 1 for an illustration. For example with three different energies this protocol drops the radiation dose down to one third compared to the option of measuring all the energies at all angles.

The attenuation images at different energies are structurally similar. Therefore, the combination of non-overlapping sampling and joint regularization model should not compromise image quality. In this paper, the proposed combination of joint reconstruction and non-overlapping sampling is evaluated with simulated multi-energy CT data and experimental data from a biological specimen.

The rest of this paper is organized as follows. The forward modeling of multi-energy CT is explained in Section 2 and the reconstruction methods used in this paper in Section 3. The generation of simulated data and measuring of the experimental data are explained in Section 4. The results with simulated and experimental measurement data are shown and discussed in Section 5. Section 6 gives the conclusions.

2 Forward Problem for Multi-Energy CT

In multi-energy CT, X-ray images of an object are acquired at two or more X-ray energies to obtain information on the energy-dependent attenuation properties of the object for differentiation and identification of materials. Multiple acquisition schemes exist for obtaining spectral CT data, detailed in, e.g., [5]. Regardless of how the multi-energy data is acquired, the forward model for multi-energy CT is typically based on the linear projection model of conventional CT.

To define our forward model, we consider a multi-energy CT experiment with $n$ X-ray energies $\{E_{k},\ k=1,\ldots,n\}$. Let $\Omega\subset\mathbb{R}^{m}$, $m\in\{2,3\}$ denote the image domain (2D area or 3D volume) and $\mu(r,E_{k}):\Omega\mapsto\mathbb{R}_{+}$ denote the attenuation function for energy $E_{k}$ with X-ray source intensity $I_{0}(E_{k})$, where $r\in\Omega$ are the spatial coordinates. With these notations, the X-ray intensity at a detector element placed behind the object is modelled by the Beer-Lambert law:

[TABLE]

where the integration path is along the straight ray from the X-ray source to the detector element [24]. By normalizing with the source intensity and taking the logarithm, we obtain the conventional linear attenuation model

[TABLE]

We remark that the linear attenuation model (2) is based on several physical approximations, namely neglecting any scattering effects and the energy dependencies of the attenuation coefficients across the X-ray spectrum for polychromatic sources. While approaches such as [47, 48, 49] have been proposed for the polychromatic model, these require solving computationally intensive nonlinear optimization problems. Figure 2 shows on the left examples of the energy-dependencies of the X-ray coefficients of three materials of medical relevance, and on the right a realistic example of an energy spectrum of an X-ray source. In principle, when only the scattering effects are neglected, an energy dependent ray propagation based attenuation model would be

[TABLE]

where $E_{\text{max}}$ is the maximum energy in the X-ray source spectrum, leading to a non-linear problem for the X-ray tomography. While the non-linear model could in principle be used, it would require accurate knowledge of the input spectrum $I_{0}(E)$ and either a large number of different measured energies (input spectra) or a spectrally resolved measurement in an attempt to recover $\mu(r,E)$ with a high energy resolution. In practice, however, the projection measurements are typically taken only with a few energies (X-ray spectra) and the measurement for each energy is approximated by the linear model (2) where $E_{k}$ is thought of as an effective energy for the respective X-ray source spectrum. The effective energy of a polychromatic source is defined as the energy of a monochromatic beam that has the same half-value layer (HVL) as the polychromatic x-ray beam in a given material, usually aluminum or copper [50]. The use of the linear approximation, however, often results in problems in image reconstruction, such as beam-hardening and metal artifacts [24].

For a computational implementation, the forward model (2) that we use must be appropriately discretized. By concatenating all the X-ray projection data that is measured for the (effective) energy $E_{k}$ into a vector $y_{k}\in\mathbb{R}^{M_{k}}$, where $M_{k}$ is the number of data points for energy $k$, the forward model for the projection data becomes

[TABLE]

where $x_{k}\in\mathbb{R}^{N}$ is the vectorized representation of the (spatial) discretization of the unknown attenuation image $\mu(r,E_{k})$ with energy $E_{k}$ and $N$ is the number of pixels/voxels, $e_{k}$ models the measurement noise, and $A_{k}$ is the matrix which implements the discretization of the projection model (2) for the particular projection geometry and projection angles used with energy $E_{k}$. The whole multi-energy CT experiment consist of a set $\{y_{k},k=1,\ldots,n\}$ data vectors of the form (4).

We remark that the number $M_{k}$ of data points for different energies in the joint reconstruction can be freely chosen and could be different for each measured energy (e.g. using fewer measurements for the higher energies).

3 Inverse Problem in Multi-Energy CT

In this section, we discuss the joint reconstruction approach and the joint regularization models that are employed in this study. However, first we will briefly review unregularized reconstruction and total variation (TV) regularized reconstruction for (single energy) CT since they will be employed as references for the joint reconstructions in the numerical examples.

3.1 Least Squares (LS) Reconstruction

When densely collected projection data for each energy is available, the image reconstructions are typically carried out by using classical methods such as the filtered backprojection (FBP) [53], Kaczmarz iterations such as the algebraic reconstruction technique (ART) [53], or iterative techniques which minimize the least squares (LS) data misfit

[TABLE]

where $y_{k}$ is the projection data (using single energy $k$), $A_{k}$ is the system matrix (or projection matrix) and $x_{k}$ is the unknown (vectorized) image.

In this paper, the LS estimate is used as one reference estimate for the joint reconstruction, giving a reference of an estimate which is obtained from the given data $y_{k}$ without prior information.

3.2 Sparsity Promoting Reconstruction in CT

The ALARA principle for X-ray dose in CT imaging has led to significant efforts in development of reconstruction methods that can provide high quality reconstructions from sparsely collected projection data. In the spirit of the theory of compressed sensing [54], the methods are usually based on some type of regularization which promotes sparsity of the unknown image by an $l_{1}$ norm regularization functional on the image on some domain. The $l_{1}$ regularization can be based, for example, on wavelets [26], shearlets [55], or the total variation (TV) which promotes piecewise regular solutions by posing the $l_{1}$ regularization for the gradient of the unknown image [56, 57, 58, 59].

Using the framework of regularized least squares minimization, the TV regularized reconstruction can be obtained as

[TABLE]

where

[TABLE]

is the (smooth) total variation functional [60], $\beta>0$ is a smoothing parameter, $\gamma_{k}>0$ is a regularization parameter, and $r$ is the spatial coordinate vector.

In multi-energy CT, one can reconstruct the images at different energies using the TV regularization separately for each energy. In this paper, this approach is used as the other reference for the joint reconstructions, giving a reference of sparsity promoting reconstructions which are computed independently for each energy without any prior model for the correlation of the images between the energies.

3.3 Joint Reconstruction in Multi-Energy CT

In joint reconstruction, the central idea is to combine all the data sets $\{y_{1},\ldots,y_{n}\}$, where $n$ is the number of energies, together with a joint regularization functional into a single problem of reconstructing all the unknown images $\{x_{1},\ x_{2},\ \ldots,\ x_{n}\}$ at once. In the regularized least squares framework, the reconstructions can be obtained as

[TABLE]

where $y_{k}$ is the measurement vector and $A_{k}$ is the system matrix corresponding to energy $k$, and

[TABLE]

is the regularization functional where $U(x_{k})$ is a (spatial) regularization functional acting separately on each of the images $x_{k}$, $W(x_{1},x_{2},\ldots,x_{n})$ is a joint regularization functional which is used to incorporate prior information across the different energies, and $\gamma_{k}$ and $\alpha$ are regularization parameters.

The regularization functional $W(x_{1},x_{2},\ldots,x_{n})$ should be formulated to promote feasible a priori information between the unknown images at different energies. In multi-energy CT, it is reasonable to assume that the attenuation images at different energies are structurally similar in the sense that if an edge is present with one energy, it is likely to be at the same location and alignment with other energies as well.

The following subsections introduce the joint regularization functionals used in this paper.

3.3.1 Joint Total Variation

The joint total variation (JTV) functional was originally developed for denoising and deblurring of the red, green, and blue channels in RGB images [41]. It can be defined as

[TABLE]

and it promotes sparsity of the joint gradient, leading to edge preserving regularization within each image and favoring the locations of edges to be the same in the images $x_{k}$.

In the JTV regularization, in (9) we set $U(x_{i})=0$ and

$W(x_{1},x_{2},\ldots,x_{n})=JTV(x_{1},x_{2},\ldots,x_{n})$ leading to

[TABLE]

in (8).

3.3.2 Linear Parallel Level Sets

The linear parallel level sets (LPLS) prior was originally used for denoising and demosaicking RGB images [42]. In tomographic context, it has been used in joint reconstruction of PET and MRI images [61], reconstruction of PET and EIT images with MRI images as side information [62, 63] and in multi-energy computed tomography [35].

The LPLS prior is based on the idea that level sets can be used to identify the structure in images and that the gradients of an image are always perpendicular to the level sets. Therefore, structural similarity can be measured by measuring the parallelism of the gradients’ orientations in two images. In the context of a minimization problem, this can be formulated as [62]

[TABLE]

and in this paper, the computation is done in pairs, so that

[TABLE]

In the LPLS regularization, we set $U(x_{i})=0$ and

$W(x_{1},x_{2},\ldots,x_{n})=LPLS(x_{1},x_{2},\ldots,x_{n})$ in (9), leading to the regularization functional

[TABLE]

3.3.3 Spectral Smoothness Regularization

A straightforward way to promote structural similarity is to use smoothness regularization in the spectral (energy) dimension. Spectral smoothness regularization can be obtained simply by penalizing differences between the images $\{x_{k}\}$. In this paper, we consider a spectral smoothness regularization term based on the first finite difference but also higher order finite differences could be used. Somewhat similar smoothness penalties have been used for example in dynamic magnetic resonance imaging [43] and X-ray fluorescence tomography [44].

The first finite difference can be used for regularization with

[TABLE]

and it favors images that are close to being identical. The parameter $w_{k}$ in (15) basically corresponds to the reciprocal of the (spectral) distance (energy difference) between the images $x_{k+1}$ and $x_{k}$. However, in cases where one might have prior information about the spatial variations in the expected magnitudes of the changes between the images $x_{k+1}$ and $x_{k}$, the parameter $w_{k}$ could be made spatially distributed for construction of a spatially weighted regularization model for the spectral smoothness. For more than two images $x_{i}$ we can write

[TABLE]

In this paper, we set the weights $w_{k}\equiv 1$, and use the regularization functional

[TABLE]

for D1 regularization.

Notice that the D1 regularizer imposes a penalty only for the differences of the pixel values in the energy direction. When using highly sparse projection data, this lack of spatial constraint can result in spatially noisy images. In such cases, an additional spatial regularization term, such as total variation (7), can be used to promote regularity within each of the images. This can be done by choosing $U(x_{i})=TV(x_{i})$ in (9), leading to

[TABLE]

for the D1+TV regularization.

3.3.4 Structural Similarity Index Based Regularization

A new joint regularization model which favors structural similarity can be formulated by considering the structural similarity (SSIM) index [45, 64] which is an image metric that assesses the similarity of two images by using three characteristics of the input images: luminance, contrast, and structure. Here we consider only the structure part which associates the unit vectors $(x_{1}-\mu_{x_{1}})/\sigma_{x_{1}}$ and $(x_{2}-\mu_{x_{2}})/\sigma_{x_{2}}$ with the structure of the two images and uses their inner product to measure the structural similarity, resulting in

[TABLE]

where $\sigma_{x_{1}x_{2}}$ is the cross correlation, $\sigma_{x_{1}}$ and $\sigma_{x_{2}}$ are the standard deviations, and $\mu_{x_{1}}$ and $\mu_{x_{2}}$ are the expected values of $x_{1}$ and $x_{2}$, and $C$ is a small constant introduced to avoid division by zero. Geometrically, the cross correlation $\sigma_{x_{1}x_{2}}$ corresponds to the cosine of the angle between the vectors $x_{1}-\mu_{x_{1}}$ and $x_{2}-\mu_{x_{2}}$, and thus (19) is maximized when these vectors point in the same direction and minimized when they point in opposite directions.

As with the SSIM index, it is possible to compute an image showing the similarity in structure of two images. Following the procedure for computing a SSIM map [45], but only for the structure part (19), the cross correlations and standard deviations, as well as $\hat{S}$ itself, are computed locally using a 11x11 window that moves pixel by pixel over the entire image. If we let $\mathcal{N}_{i}$ denote the set of indices for the 11x11 window for pixel $i$ in the image lattice, and let $w=\{w_{k}|k\in\mathcal{N}_{i}\}$ denote the local weights obtained from a circular symmetric Gaussian weighting function with standard deviation of 1.5, local weighted statistics for window $i$ can be computed as

[TABLE]

where $(x_{k})_{j}$ is the $j$th element of $x_{k}$. The obtained pixel-by-pixel local $\hat{S}$ form the structural similarity map with larger values indicating similarity and smaller values indicating structural differences. Now we can use the mean of the local $\hat{S}$ values

[TABLE]

where $N$ is the number of pixels in the images and $(\hat{S}(x_{k,i},x_{p,i}))_{i}$ is the value of $\hat{S}$ in pixel $i$, as a scalar measure of the structural similarity which is maximized with the value 1 when the images $x_{k}$ and $x_{p}$ are identical. Now the regularization term in (8) can be chosen to be some decreasing function of all the $\bar{S}(x_{k},x_{k+1})$ which has minimum when the image pairs are identical. In this paper, we use a simple choice

[TABLE]

where

[TABLE]

computes the structure parts in pairs.

When employing the structural similarity regularization alone, we set $U(x_{i})=0$ and use

[TABLE]

as the regularization functional in (8) for S regularization.

Similarly as with the spectral smoothness regularization, the lack of penalty for spatial regularity within each channel in the structural regularization term (24) can result in spatially noisy reconstructions when using sparse projection data. In such cases, an additional spatial regularization term, such as total variation (7), can be used to promote spatial regularity of the reconstructions. In this paper, we use $U(x_{i})=TV(x_{i})$ in (9), leading to

[TABLE]

for the S+TV regularization.

4 Materials and Methods

In this section, we present the simulated and measured X-ray projection data that were used to test the different regularization models utilized in multi-energy CT reconstruction. The datasets were generated using dense angle sampling, and could then later be appropriately subsampled to simulate sparse-angle acquisition schemes.

4.1 Simulated Data

Simulated X-ray projection data was generated using idealized, monochromatic sources in a parallel beam geometry. For this, a computational anthropomorphic 2D phantom was created using the XCAT software [65, 66]. The phantom consisted of a single cross-sectional slice of an adult human male chest. A calcified coronary artery was placed in the heart of the phantom. We simulated three sets of noise-free projections using 40 keV, 80 keV, and 120 keV X-rays. Each set consisted of 360 projections taken evenly across [math] at [math] intervals. The X-ray detector was a line detector consisting of 729 elements, with a pixel size of 0.875 mm. To avoid committing inverse crime, the projections were generated using XCAT’s own CT Projector tool.

The noisy realizations of the simulated measurements for each energy were obtained by adding Gaussian zero mean random noise with standard deviation equal to one percent of the maximum amplitude of the measurements with that energy to the simulated noise-free measurements.

The reconstruction grid size was $512\times 512$ pixels. The forward operators for the model $Ax_{k}=y_{k}$ were created using the ASTRA Toolbox [67, 68, 69].

4.2 Experimental Data

Experimental data was measured at the Department of Physics, University of Helsinki, using a cone-beam micro-CT scanner with an end-window tube and a tungsten target (GE Phoenix nanotom 180 NF). The chest of a small bird (the common quail, Coturnix coturnix) was used as a test phantom, as it contains multiple different tissue types as well as fine details arising from the bone structure. The bird phantom was imaged using three different X-ray tube settings in the same geometry. Multiple frames were averaged for each projection in order to increase signal-to-noise ratio. The imaging geometry is specified in Table 1 and the settings for each effective energy are specified in Table 2. 2D sinograms were created using the central plane of the cone-beam projections, in which the geometry reduces to a fan-beam geometry.

The reconstruction grid size was $512\times 512$ pixels with a pixel size of 120 $\mathrm{\SIUnitSymbolMicro}\mathrm{m}$, and the forward operators for the imaging model were created using the ASTRA Toolbox.

4.3 Image Fidelity Measures

To quantify the quality of reconstructions in the simulation test case, we use the root mean square error (RMSE) and the mean structural similarity (MSSIM) index [45].

The root mean square error can be written as

[TABLE]

where $(x_{k})_{i}$ denotes the $i$th pixel of the solution $x$ corresponding to energy $k$ and $x_{\textrm{REF},k}$ is the reference image corresponding to energy $k$.

The structural similarity (SSIM) index [45] aims to measure the perceived similarity in two images, $x_{1}$ and $x_{2}$. It assesses the luminance $l(x_{1},x_{2})$, contrast $c(x_{1},x_{2})$ and structure $s(x_{1},x_{2})$ and can be written as

[TABLE]

where $\mu$ is the mean, $\sigma$ is the standard deviation, $\sigma_{x_{1}x_{2}}$ is the cross-correlation and $C_{1}$, $C_{2}$ and $C_{3}$ are small constants used to avoid division by zero. The mean SSIM index is obtained by first computing the SSIM index locally for each pixel using a 11x11 pixel window whilst also using a Gaussian weighting function for the statistics, as in (20 – 22), and then the final index is obtained as the mean of these local SSIM indices. For a reconstructed image $x_{k}$ and a reference image $x_{\textrm{REF},k}$, this can be written as

[TABLE]

where $(\textrm{SSIM}(x_{k},x_{\textrm{REF},k}))_{i}$ is the local value of SSIM for pixel $i$ in the image lattice.

As the anthropomorphic XCAT target phantom that was used as the basis for the projection data does not translate well into attenuation images that could be used as ground truth reference images, the ground truth reference images $x_{\textrm{REF},k}$ for the simulation test case were obtained by computing full noiseless data (360 angles for each energy) filtered backprojection (FPB) reconstructions using the ASTRA Toolbox [67, 68, 69] FBP algorithm.

For the test case with experimental data, the reference images $x_{\textrm{REF},k}$ were obtained by the ’No Prior’ method in (5) applied to 720 projections for each measured energy.

4.4 Selection of the Regularization Parameters

For the simulation test case, the regularization parameters $\gamma_{k}$ and $\alpha$ in (6) and (8) and $\beta$ for each regularization term were selected by computing a series of reconstructions with varying $\gamma_{k}$, $\alpha$ and $\beta$. Then the mean structural similarity index (30) was computed for each reconstructed image using the full data FPB reconstructions as the reference images $x_{\textrm{REF},i}$, and the $\gamma_{k}$, $\alpha$ and $\beta$ corresponding to the highest structural similarity were selected for the test cases. A similar procedure was used in the experimental data test case but the identification of best reconstructions was done with visual comparison. The selected $\gamma_{k}$ and $\alpha$ for each regularization term are given in Table 3. For $\beta$, the value of 1e-6 was used for the simulated data test case and 1e-9 for the experimental data test case. These values are small compared to the magnitudes of the gradients in Equations (7), (10) and (3.3.2), resulting to an acceptable level of smoothing in the reconstructions.

4.5 Solution of the image reconstruction problems

The minimization problems (5), (6) and (8) were solved using a Polak-Ribière nonlinear conjugate gradient method [70] where the search direction was reset to the deepest descent direction whenever the Polak-Ribière $\beta^{PR}$ parameter went negative. The step lengths were chosen using a golden section line search and the non-negativity constraints were implemented using projection. The iteration was terminated when the change in the minimized function or the norm of the change in the parameters was less than 1e-9, or after 512 iterations. The maximum number of iteration limit was reached only when computing the no prior solution (5). Using a smaller termination criterion or a larger maximum iteration limit did not lead to noticeable changes in the reconstructions. We remark that the nonlinear Polak-Ribière conjugate gradient algorithm does not possess global convergence results under the present setup where we minimize different non-linear functionals under non-negativity constraints for the solution. However, the algorithm implementation for the different regularization models was tested with simulated data with various targets, with and without simulated measurement errors, and with full and sparse measurement data. The minimization algorithm converged to a (possibly local) minimum and produced feasible reconstructions in all the test cases.

5 Results

In this section, reconstructions computed using the joint reconstruction approaches in Section 3.3 are presented and compared to the full data reference images and to the single-energy reference methods (5) (no prior) and (6) (TV) where attenuation at each energy is reconstructed separately. The first subsection introduces the multi-energy projection geometry, the second subsection shows results computed with simulated data and the third subsection shows results computed with experimental measurement data.

5.1 Multi-Energy Projection Geometry

The principle of the low dose multi-energy projection geometry is illustrated in Figure 1 for an experiment with three different energies, denoted by the X-ray tube voltages $\mathrm{kV}_{1}$, $\mathrm{kV}_{2}$ and $\mathrm{kV}_{3}$ in the picture. Instead of sampling all the projection directions with all the energies, the directions are divided to non-overlapping sets of interleaved projection directions making the projection geometry different for each of the used source energies. This type of imaging setup could be realistically implemented in a clinical CT scanner using a rapid-kVp switching scheme.

In this study, the non-overlapping low dose data sets were obtained by retrospective subsampling of the full angle data sets with three different source energies. In both the simulation and experimental data cases, a subset of 90 equally spaced projection angles were chosen from the whole angular span of the projection geometry. This subset was then further subsampled to three subsets of 30 angles, leading to 30 equally spaced projection directions for each energy, see Figure 1 and Table 4. For each energy, the projection data corresponding to the set of 30 angles was used in the computations.

5.2 Simulated Data

The reconstructed images for the simulated data test case are shown in Figure 3. The first three columns show the reconstructions for energies 40, 80 and 120 keV, respectively, and the next three columns show a region of interest from the corresponding reconstructions. Figure 4 shows the RMSE and MSSIM fidelity metrics for the reconstructions.

All the joint reconstructions (JTV, LPLS, D1, S, D1+TV, S+TV) in Figure 3 are visually better than the independent single energy reconstructions (No prior, TV) in the sense that more of the details that are present in the full data reference images can be detected in the reconstructions. This could be expected since the sub-problems with different energies and projection directions are coupled to each other in the joint reconstruction approaches via the joint regularization functional. The (S+TV) reconstructions have the best fidelity in the RMSE measure and the (LPLS) reconstructions in the MSSIM measure, see Fig 4. One interesting feature in the reconstructions is the difference in the amount of details in the reconstructions with the (JTV, LPLS) models compared to the reconstructions with the (D1, S) models. The JTV and LPLS models promote spatially piecewise regular solutions by posing different penalties for the sparsity of the image gradients. Consequently, the reconstructions lack some of the fine details but on the other hand have a very low level of spatial noise and very good RMSE and MSSIM metrics. In contrast, the reconstructions with the (D1, S) models show more details but are somewhat noisy, which is also reflected in the fact that the (D1, S) models have worse RMSE and MSSIM metrics than the (TV) model. This behavior is due to the fact that the (D1, S) models are based on promoting structural similarity in the energy direction without any penalty on spatial regularity. A good combination between the spatial details and noise level is reached when the (D1, S) models are used in combination with the TV regularization applied to each of the images for promoting spatial regularity (D1+TV, S+TV).

The results in Figures 3 and 4 reveal that the joint reconstruction models outperform the independent reconstruction models (No prior, TV) when using the same projection data of a total of 90 projection images. Figure 5 shows a comparison of TV and S+TV reconstructions using the same measurement data of 90 projection images with the 30+30+30 direction non-overlapping geometry that was used in Fig. 3 and also reconstructions using measurement data of 270 projection images using the same 90 directions for each energy. Error metrics for the reconstructions are shown in Figure 6. By comparison of the reconstructions S+TV(90) and TV(90), the S+TV outperforms the TV reconstructions also in this case when using classical overlapping projection geometry with the same 90 directions measured at each energy. More importantly, the comparison of TV(90) to S+TV(30W) reveals that they both perform approximately equally well even though S+TV(30W) uses only 90 projection images instead of the 270 used by TV(90). This finding suggests that the proposed combination of joint reconstruction and non-overlapping projection direction geometry can facilitate significant reduction of the x-ray dose.

5.3 Experimental Data

The reconstructed images for the experimental data test case are shown in Figure 7. The first three columns show the reconstructions for $E_{1}$, $E_{2}$, and $E_{3}$, respectively, and the next three columns show a region of interest from the corresponding reconstructions. For reference, the first row (Ref) shows ground truth reconstructions computed with the least squares (No prior) model using the full 720 projection angle data for each energy.

The reconstructions of the experimental data test case support the findings of the simulation test case. The joint reconstructions (JTV, LPLS, D1, S, D1+TV, S+TV) are visually better than the reference single energy reconstructions (No prior, TV) and there is a difference in the amount of details between (TV, JTV, LPLS) and (D1, S), visible for example from the outer edge of the target in the zoomed in images. The combined regularization approaches (D1+TV, S+TV) provide again a good combination of the spatial details and noise level, noticeable for example in the visibility of the two small low value areas on the right hand side of the zoomed in images.

6 Discussion and Conclusions

In this paper, we studied the performance of various structural similarity promoting joint reconstruction models using a novel low dose data acquisition protocol. The models were tested using both simulated and experimental measurement data and in both cases the linear parallel level sets (LPLS), joint total variation (JTV) and the combined structural similarity index based and total variation regularization (S+TV) performed the best. The regularization models that utilized only spectral regularization (S, D1) performed worse with respect to the RMSE and MSSIM fidelity metrics because they resulted in noisier reconstructions but nevertheless showed more details than the non-joint reconstruction approaches (No prior, TV). The combination of the joint reconstruction and the novel low dose projection protocol with non-overlapping projection directions was also compared to TV reconstructions using classical projection geometry, and was found to tolerate a clear reduction in the number of projection images, suggesting that the proposed approach can facilitate an effective way to reduce x-ray dose in multi-energy CT.

An implicit assumption in the proposed approach is that the target remains stationary during the imaging with different energies. However, with a data acquisition protocol that uses mutually non-overlapping projection angles and mid-scan energy switching, these motion issues should not, at least in theory, be any more problematic than with regular CT scanning. In addition, since the proposed combination of joint reconstruction and non-overlapping projection directions can operate on a smaller number of projection images than the combination of conventional projection geometries and reconstruction methods, it could potentially also facilitate a shorter overall scanning time to mitigate the risks of motion artifacts.

In this study, we employed smoothed versions of the total variation and linear parallel level set regularization functionals and computed the optimization solutions using a Polak-Ribière non-linear conjugate gradient method. However, using the non-smooth versions of the regularization functionals can lead to sharper edges in the reconstructions and the optimization could be handled using well-known algorithms such as the Chambolle-Pock algorithm [71, 72] or the primal dual hybrid gradient algorithm [73]. The approach used here, on the other hand, has the benefit that the optimization using the smoothed functionals requires significantly fewer iterations (and thus less forward and back projection operations) than the minimization of the non-smooth functionals, implying significantly faster computation times, especially in large-scale reconstruction problems with three spatial dimensions or a temporal dimension, or both.

The joint reconstruction approach used in this study can be applied directly to multi-energy CT data obtained with energy discriminating detectors. In such a case, the target would be imaged for each angle using a single exposure and the detector would distribute the measured photons to different energy bins, implying that each energy channel is measured at all projection angles. However, as indicated by comparison of S+TV(90) and TV(90) in Figures 5 and 6, the proposed joint reconstruction approach would be beneficial also when working with data with overlapping projection angles.

In this study, we employed sparse data which contained 30 projections at each energy, leading to an overall of 90 projection directions. In practical applications, the critical number of angles needed could be sought for by a similar simulated phantom based study that is proposed for limited angle tomography of the breast in [27].

Acknowledgments

This work was supported by the Academy of Finland (Project 312343, Finnish Centre of Excellence in Inverse Modelling and Imaging), the Jane and Aatos Erkko Foundation, and Business Finland project 6614/31/2016.