Derivative-Free Iterative One-Step Reconstruction for Multispectral CT

Abstract

Image reconstruction in multispectral computed tomography (MSCT) requires solving a challenging nonlinear inverse problem, commonly tackled via iterative optimization algorithms. Existing methods necessitate computing the derivative of the forward map and potentially its regularized inverse. In this work, we present a simple yet highly effective algorithm for MSCT image reconstruction, utilizing iterative update mechanisms that leverage the full forward model in the forward step and a derivative-free adjoint problem. Our approach demonstrates both fast convergence and superior performance compared to existing algorithms, making it an interesting candidate for future work. We also discuss further generalizations of our method and its combination with additional regularization and other data discrepancy terms.

1. Introduction

Classical computed tomography (CT) is based on the inversion of the linear Radon transform, where a scalar-valued attenuation map $\mu :\mathcal{X}\to \mathbb{R}$ of the patient is recovered from observation of its Radon transform $\mathcal{R}\mu :\mathcal{L}\to \mathbb{R}$ derived from projection data. Here and below $\mathcal{X}\subseteq {\mathbb{R}}^d$ is the image domain in $d=2,3$ dimensions, and $\mathcal{L}$ is the set of integration lines. While sufficient in many applications, the linear problem ignores the polychromatic nature of the X-rays and the energy-dependent absorption characteristics of real-world objects. The sample is more accurately represented by a family of attenuation maps $\mu \left(e\right):\mathcal{X}\to \mathbb{R}$ dependent on the photon energy $e\in (0,\infty )$. Recovering a single $\mu$ from projection data using a single energy bin results in a mixture of density maps from different energies resulting in severe nonuniqueness. Additionally, the nonlinearity results in beam hardening artifacts that may be partially accounted for by iterative algorithms or analytic modeling [1,2,3,4,5,6]. In order to overcome such weaknesses, the idea of multispectral CT (MSCT) is to measure projection data for different energy bands, which are then used to reconstruct multiple attenuation maps. The reconstruction problem, however, becomes nonlinear and much more challenging than pure Radon inversion [7,8,9,10,11,12,13]. In this work, we develop a simple and efficient strategy for tackling the nonlinearity.

1.1. Multispectral CT

Specifically, in this work, we use the material decomposition paradigm in MSCT. In the material decomposition approach, it is assumed that the energy-dependent attenuation maps $\mu \left(e\right)$ can be written as $\mu (e,x)={\sum }_{m=1}^M{\mu }_m\left(e\right)f_m\left(x\right)$, where $f_1,f_2,\dots ,f_M:\mathcal{X}\to \mathbb{R}$ are the densities of M separate materials (with $f_m\left(x\right)\in [0,1]$) to be recovered and ${\mu }_m\left(e\right)$ are known and tabled absorption characteristic of the m-th material. The aim is to recover the material densities from data collected for several energy bands. This does not only allow to improve image quality, but also offers a broad range of applications, as it reconstructs multiple images encoding different characteristics of specific regions, which enables a deeper understanding of the object under examination. Recent significant advancements in the manufacturing of energy-sensitive sensors [14,15] has considerably increased the interest in and applicability of MSCT.

Assuming B spectral measurements $Y_1,\dots ,Y_B$, the material decomposition problem in MSCT can be written as the problem of recovering $f_1,\dots ,f_M$ based on the modeling equation ${\left(Y_b\right)}_{b=1}^B=\Phi \left({\left(\mathcal{R}{f}_m\right)}_{m=1}^M\right)$. Here, $\mathcal{R}{f}_m$ is the Radon transform applied to the m-th material density map $f_m$, and $\Phi$ is a nonlinear map. Classical CT would correspond to the unrealistic case where the application of the pointwise logarithm in the modeling equation results in a linear inverse problem. In MSCT, one accounts for the full nonlinear problem for the unknown density maps $f_m$.

1.2. Two-Step and One-Step Algorithms

Various algorithms have been developed for solving the reconstruction problem in MSCT. They can be broadly classified into two categories: two-step methods and one-step algorithms. The idea of earlier two-step methods is to perform Radon inversion and material decomposition in two separate steps. Material decomposition can be performed either in the projection domain $\mathcal{L}$ (before Radon inversion) or in the image domain $\mathcal{X}$ (after Radon inversion). Both methods have their specific advantages and disadvantages. The image-domain decomposition approach allows incorporating prior information about the objects that is naturally contained in the image domain $\mathcal{X}$. However, the nonlinear nature of the problem leads to approximate linear models that introduce severe reconstruction artifacts. The projection-domain decomposition approach, on the other hand, allows working with the correct nonlinear model. However, the prior structure in the Radon domain is not directly available. See the works [12,13,16,17,18] and references therein for proposed two-step approaches.

One-step methods reconstruct the material densities $f_1,\dots ,f_M$ through iterative minimization techniques and thus overcome the drawbacks of both two-step methods. For some one-step algorithms, we refer to [7,8,11,19,20,21,22,23,24]. Despite their superior performance, such one-step iterative algorithms are computationally expensive. Existing methods require many iterative steps due to poor conditioning of the problem with computationally expensive iterative steps. The algorithms proposed in this paper are specific one-step algorithms that address these two drawbacks of existing one-step methods. The structure of our proposed algorithm is shown in Figure 1.

1.3. Our Contributions

The following list summarizes the main contributions of this paper:

• We present a novel derivative-free algorithm designed to combine the advantages of one-step and two-step approaches. To achieve this, we introduce a simple and computationally efficient iterative update that incorporates appropriate preconditioning.
• Image reconstruction is performed in the image domain, which naturally allows for the inclusion of an image smoothness prior. Our method can be combined with additional regularization. However, in order to show the method in its pure form, we will not include such a modification.
• Our methods integrate benefits of two-step approaches by separating iterative updates into two parts. Moreover, the main ingredient that makes the algorithm efficient is the use of the full nonlinear forward model but linearization around zero for the adjoint problem. In addition to avoiding computation and evaluation of the derivative of the forward map, this also allows for including simple but efficient channel preconditioning.

2. Mathematical Modelling of MSCT

We assume that the object to be imaged lies in some domain $\mathcal{X}\subseteq {\mathbb{R}}^2$ and consists of a combination of M different materials with densities $f_m:\mathcal{X}\to \mathbb{R}$ with $m=1,\dots ,M$. Each material has a separate mass attenuation coefficient ${\mu }_m:[0,\infty )\to [0,\infty )$, which is a known function of the X-ray energy e. The total energy dependent (linear) X-ray attenuation coefficient is then given by:

$$\mu (x,e)=\sum _{m=1}^M{\mu }_m\left(e\right)f_m\left(x\right)\mathrm{for}(x,e)\in \mathcal{X}\times \mathbb{R}.$$

(1)

Assuming that the material specific attenuation functions ${\mu }_m\left(e\right)$ are known, the goal is to recover densities $f_m$ from indirect x-ray measurements using different energy bins.

2.1. Continuous Model

We start with continuous modeling, where the quantities involved are functions on continuous domains that will be discretized later. Suppose that X-ray energy with a known incident spectral density $I_0\left(e\right)$ is sent along a line $\ell \in \mathcal{L}$ from the source position to the detector position. While propagating along ℓ, the X-rays are attenuated according to $\mu (x,e)$ defined in (1). This results in an outgoing spectral density $I_1\left(e\right)=I_0\left(e\right)exp(-{\int }_{\ell }\mu (x,e)d\ell \left(x\right))$ at the detector. The energy-sensitive detector with the spectral profile $\widehat{s}\left(e\right)$ records the integral ${\int }_0^{\infty }\widehat{s}\left(e\right)I_1\left(e\right)de$. Denoting the product of the incident spectral density of the source and the spectral sensitivity of the detector by $s\left(e\right)\triangleq I_0\left(e\right)\widehat{s}\left(e\right)$, referred to as the effective spectrum, the recorded data are given by:

$$\begin{array}{cc}\hfill Y(s,\ell )& ={\int }_0^{\infty }s\left(e\right)exp\left(-{\int }_{\ell }\mu (x,e)d\ell \left(x\right))\right)de\hfill \\ \hfill & ={\int }_0^{\infty }s\left(e\right)exp\left(-\sum _{m=1}^M{\mu }_m\left(e\right){\int }_{\ell }{f}_{m}d\ell \left(x\right)\right)de.\hfill \end{array}$$

(2)

The data in Equation (2) represent a single measurement in MSCT. The goal of material decomposition in MSCT is to determine the density distributions $f_m$ from multispectral X-ray measurements obtained by varying the line ℓ and the effective spectra s.

For simplicity of presentation, we consider parallel beam mode, where any line $\ell =\ell (\theta ,r)$ is parametrized by its normal vector $\theta$ and its distance r from the origin. In this case, ${\int }_{\ell }{f}_{m}d\ell \left(x\right)=\mathcal{R}{f}_m(\theta ,x)$ is given by the classical Radon transform of $f_m$. Assuming further a total number of B different effective spectra and writing $f=(f_1,\dots ,f_M)$, we obtain the continuous MSCT forward model:

$$Y={\left[{\int }_0^{\infty }{s}_b\left(e\right)exp\left(-\sum _{m=1}^M{\mu }_m\left(e\right)\mathcal{R}\left(f_m\right)\right)de\right]}_{b=1,\dots B}.$$

(3)

Equation (3) gives the complete continuous forward model in material decomposition in multispectral CT. The unknown f consists of M functions $f_1,\dots ,f_M$ defined on the image domain $\mathcal{X}$ and the data of B functions $Y_1,\dots ,Y_B$ defined on the projection domain $\mathcal{L}={\mathbb{S}}^1\times \mathbb{R}$. The methods that we describe, however, would also work with a three-dimensional image domain and a general projection domain $\mathcal{L}$ of lines in ${\mathbb{R}}^3$.

2.2. Discretization

In order to avoid technical details and to concentrate on the main ideas, we derive the algorithm for the discrete forward model throughout this paper. For that purpose, we represent the material densities via a discrete column vector $X_1,\dots ,X_M\in {\mathbb{R}}^{N_x}$ and the Radon transform via a matrix $\mathbf{A}\in {\mathbb{R}}^{N_y\times N_x}$, where $N_x$ is the number of discretization points in the image domain and $N_y$ is the number of lines used in the projection domain. Furthermore, we discretize the effective energy spectra by vectors ${\mathbf{S}}_1,\dots ,{\mathbf{S}}_B\in {\mathbb{R}}^E$ and the known material attenuations by vectors ${\mu }_1,\dots ,{\mu }_M\in {\mathbb{R}}^E$, where E is number of discretization points of the energy variable. The discretization of (3) yields the following discrete image reconstruction problem.

Problem 1 (Discrete MSCT image reconstruction problem).

Recover the unknown $X\in {\mathbb{R}}^{N_x\times M}$ from data $Y=\mathbf{F}\left(X\right)+\delta \in {\mathbb{R}}^{N_y\times B}$ where:

$$\mathbf{F}\left(X\right)\triangleq {(\mathbf{S}·{exp}_{E\times N_y}(-\mathbf{M}·{\left(\mathbf{A}X\right)}^{⊺}))}^{⊺}$$

(4)

Here and below, we use the convention that the boldface notation ${exp}_{E\times N_y}$ indicates that the scalar function exp is applied pointwise to a vector in ${\mathbb{R}}^{E\times N_y}$. Furthermore, in (4):

• The columns of $X=[X_1,\dots ,X_M]$ are the discrete material images;
• $\mathbf{A}\in {\mathbb{R}}^{N_y\times N_x}$ is the discretized Radon transform;
• The columns of $\mathbf{M}=[{\mathbf{M}}_1,\dots ,{\mathbf{M}}_M]\in {\mathbb{R}}^{E\times M}$ are the discretized material attenuations;
• The columns of ${\mathbf{S}}^{⊺}=[{\mathbf{S}}_1^{⊺},\dots ,{\mathbf{S}}_B^{⊺}]\in {\mathbb{R}}^{E\times B}$ the discretized effective spectra;
• The columns of $Y=[Y_1,\dots ,Y_B]$ are the observed spectral data.

Note that we included the transpose operations in (4) such that all involved linear operations can be written as matrix multiplications from left. Alternatively, we can also write $\mathbf{F}\left(X\right)={exp}_{N_y\times E}(-\mathbf{A}X{\mathbf{M}}^{⊺}){\mathbf{S}}^{⊺}$, reflecting the fact that the discrete Radon transform $\mathbf{A}$ operates on the columns and the matrices $\mathbf{M}$ and $\mathbf{S}$ on the rows of $X\in {\mathbb{R}}^{N_y\times M}$. At some places, we will denote operation of $\mathbf{M}$ on X from the right by $\mathbf{M}\diamond X\triangleq X·{\mathbf{M}}^{⊺}$ such that we have $\mathbf{F}\left(X\right)=\mathbf{S}\diamond {exp}_{N_y\times E}(-\mathbf{M}\diamond \mathbf{A}X)$. The structure of the discrete forward model is illustrated in Figure 2.

Remark 1 (Recalibration).

With $\langle {\mathbf{S}}_b,\mathbf{1}\rangle ={\sum }_{e=1}^E{\mathbf{S}}_{b,e}$ we obtain:

$$\mathbf{F}{\left(0\right)}^{⊺}=\mathbf{S}·{\mathbf{exp}}_{E\times N_y}\left(0\right)=\mathbf{S}·{\mathbf{1}}_{E\times N_y}=\left(\begin{array}{ccc}\langle \mathbf{1},{\mathbf{S}}_1\rangle & \dots & \langle \mathbf{1},{\mathbf{S}}_1\rangle \\ ⋮& & ⋮\\ \langle \mathbf{1},{\mathbf{S}}_B\rangle & \dots & \langle \mathbf{1},{\mathbf{S}}_B\rangle \end{array}\right).$$

(5)

Thus, with ⊙ denoting pointwise multiplication (also known as Hadamard product) and ⊘ the pointwise division, we get:

$$\begin{array}{c}\mathbf{F}{\left(X\right)}^{⊺}\oslash \mathbf{F}{\left(0\right)}^{⊺}=\left(\begin{array}{ccc}1/\langle \mathbf{1},{\mathbf{S}}_1\rangle & \dots & 1/\langle \mathbf{1},{\mathbf{S}}_1\rangle \\ ⋮& & ⋮\\ 1/\langle {\mathbf{S}}_b,\mathbf{1}\rangle & \dots & 1/\langle {\mathbf{S}}_b,\mathbf{1}\rangle \end{array}\right)\odot (\mathbf{S}·{\mathbf{exp}}_{E\times N_y}(-\mathbf{M}·{\left(\mathbf{A}X\right)}^{⊺}))\hfill \\ \hfill =\left(\begin{array}{c}{\mathbf{S}}_1/\langle \mathbf{1},{\mathbf{S}}_1\rangle \\ ⋮\\ {\mathbf{S}}_B/\langle \mathbf{1},{\mathbf{S}}_B\rangle \end{array}\right)·{\mathbf{exp}}_{E\times N_y}(-\mathbf{M}·{\left(\mathbf{A}X\right)}^{⊺}).\end{array}$$

(6)

This means that the recalibrated forward model $\mathbf{F}{\left(X\right)}^{⊺}\oslash \mathbf{F}{\left(0\right)}^{⊺}$ is the same as the original forward model with normalized effective spectra ${\mathbf{S}}_b/\langle {\mathbf{S}}_b,\mathbf{1}\rangle$. The matrix with normalized spectra can be written as $\mathbf{S}\oslash (\mathbf{S}·{\mathbf{1}}_{E\times E})$.

In this work we, use rescaled data, to which we apply the pointwise logarithm ${log}_{N_y\times B}$ and the corresponding least squares (LSQ) functional.

Definition 1 (Forward model and LSQ functional).

We define the logarithmic forward model and the LSQ functional by:

$$\begin{array}{cc}\hfill & \mathbf{H}:{\mathbb{R}}^{N_x\times M}\to {\mathbb{R}}^{N_y\times B}:Y\mapsto {\mathbf{log}}_{N_y\times B}(\mathbf{F}\left(X\right)\oslash \mathbf{F}\left(0\right))\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill & \mathcal{D}:{\mathbb{R}}^{N_x\times M}\to \mathbb{R}:X\mapsto {∥\mathbf{H}\left(X\right)-Y^{\mathbf{H}}∥}_2^2/2.\hfill \end{array}$$

(8)

where $\mathbf{F}\left(X\right)={\mathbf{exp}}_{N_y\times E}(-\mathbf{A}X{\mathbf{M}}^{⊺}){\mathbf{S}}^{⊺}$ is defined in (4), $Y\in {\mathbb{R}}^{N_y\times B}$ are the given data, $Y_{\mathbf{H}}={\mathbf{log}}_{N_y\times B}(Y\oslash \mathbf{F}\left(0\right))$ the modified data, and ${\mathbf{log}}_{N_y\times B}(·)$ the pointwise logarithm.

Using the notations of Definition 1, material decomposition in MSCT amounts to the near solution of $\mathbf{H}\left(X\right)=Y^{\mathbf{H}}$ or the near minimization of $\mathcal{D}$.

Remark 2 (Noise modelling).

In the statistical context, LSQ minimization derives from maximum likelihood estimation using a Gaussian noise model on ${\mathbf{Y}}_{\mathbf{H}}$. From a statistical perspective, maximum likelihood estimation for Poisson noise on Y might be more reasonable, resulting in $\mathcal{L}\left(X\right)=\langle \mathbf{F}\left(X\right),\mathbf{1}\rangle -\langle log\left(\mathbf{F}\left(X\right)\right),Y\rangle \to {min}_X$. As the focus in this paper is the derivation of an efficient reconstruction algorithm rather than statistical optimality, we work with $\mathcal{D}$. However, we expect that our strategy can also be applied to $\mathcal{L}$ instead of $\mathcal{D}$.

Remark 3 (Regularization).

Due to the ill-conditioning of $\mathbf{H}\left(X\right)=Y^{\mathbf{H}}$, the reconstruction problem has to be regularized [25,26]. In the context of LSQ minimization, a natural approach is variational regularization, where one considers ${∥\mathbf{H}\left(X\right)-Y^{\mathbf{H}}∥}_2^2/2+\alpha \mathcal{R}\left(X\right)$ instead of $\mathcal{D}$ with a regularization functional $\mathcal{R}\left(X\right)$. Recently, in [27], the plug-and-play method has been identified as a regularization technique where the regularization is incorporated by a denoiser. Another class is given by iterative regularization [28], where regularization comes from early stopping. All these regularization methods require the gradient of $\mathcal{D}$, which we compute below.

3. Algorithm Development

In this section, we derive the proposed algorithms for MSCT based on channel preconditioning (CP). The first one (CP-full) integrates channel preconditioning into a gradient scheme. In the second algorithm (CP-fast), the derivative of the forward map is replaced by the derivative at zero. Both methods greatly reduce the number of iterations compared to standard gradient methods and the numerical cost per iteration compared to Newton type methods. Before presenting our algorithms, we start by computing derivatives and gradients and recall existing gradient and Newton type methods.

3.1. Derivatives Computation

Standard algorithms for minimizing (8) require the derivative of the forward map and the gradient of the LSQ functional that we compute next. Recall the original and logarithmic MSCT forward operator $\mathbf{F},\mathbf{H}:{\mathbb{R}}^{N_x\times M}\to {\mathbb{R}}^{N_y\times B}$ and the LSQ functional $\mathcal{D}$ defined by (4), (7), and (8). We equip ${\mathbb{R}}^{N_x\times M}$ and ${\mathbb{R}}^{N_y\times B}$ with the Hilbert space structure induced by the standard inner product $\langle {\xi }_1,{\xi }_2\rangle ={\sum }_{i,j}{\xi }_1[i,j]{\xi }_2[i,j]$. Furthermore, we use ${\mathbf{H}}^{\prime }\left[X\right]:{\mathbb{R}}^{N_x\times M}\to {\mathbb{R}}^{N_y\times B}$ to denote the derivative of $\mathbf{H}$ at location $X\in {\mathbb{R}}^{N_x\times M}$ and ${\mathcal{D}}^{\prime }\left[X\right]:{\mathbb{R}}^{N_x\times M}\to \mathbb{R}$ and $\nabla \mathcal{D}\left[X\right]\in {\mathbb{R}}^{N_x\times M}$ to denote the derivative and the gradient of $\mathcal{D}$ at X, respectively.

Remark 4 (Gradients, inner products, and preconditioning).

By the definition of gradients, we have $\langle \nabla \mathcal{D}\left[X\right],\xi \rangle ={\left({\mathcal{D}}^{\prime }\left[X\right]\right)}^{*}\left(\xi \right)$, where ${(·)}^{*}$ denotes the adjoint of a linear operator. Furthermore, by the chain rule, $\nabla \mathcal{D}\left[X\right]={\left({\mathbf{H}}^{\prime }\left[X\right]\right)}^{·}(\mathbf{H}\left(X\right)-Y^{\mathbf{H}})$. Gradients and adjoints depend on the chosen inner product. For example, the inner product $\langle {\xi }_1,{\mathbf{U}}^{-1}{\xi }_2\rangle$ on the image space with a positive-definite matrix $\mathbf{U}$ yields the modified gradient $\mathbf{U}·\nabla \mathcal{D}\left[X\right]$. Choosing $\mathbf{U}$ such that $\mathbf{U}·\nabla \mathcal{D}\left[X\right]$ has improved condition significantly improves gradient based methods for minimizing $\mathcal{D}$.

Remark 5 (Some calculus rules).

For the following computation we use some elementary calculus rules listed next. Let ${\mathbf{G}}_1,{\mathbf{G}}_2:{\mathbb{R}}^{N_x}\to {\mathbb{R}}^{N_y}$ be vector valued functions and $\varphi :\mathbb{R}\to \mathbb{R}$ a scalar function with derivative $\psi :\mathbb{R}\to \mathbb{R}$. Then, for $X,\xi \in {\mathbb{R}}^{N_x}$, we have:

$$\begin{array}{cc}\hfill {({\mathbf{G}}_1\odot {\mathbf{G}}_2)}^{\prime }\left[X\right]\left(\xi \right)& =\left({\mathbf{G}}_1^{\prime }\left[X\right]\left(\xi \right)\right)\odot \left({\mathbf{G}}_2\left(X\right)\right)+\left({\mathbf{G}}_1\left(X\right)\right)\odot \left({\mathbf{G}}^{\prime }\left[X\right]\left(\xi \right)\right)\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill {\mathit{\varphi }}_N^{\prime }\left[X\right]\left(\xi \right)& ={\mathit{\psi }}_N\left(X\right)\odot \xi .\hfill \end{array}$$

(10)

As usual, we define the vector value functions ${\mathit{\varphi }}_N,{\mathit{\psi }}_N:{\mathbb{R}}^{N_x}\to {\mathbb{R}}^{N_x}$ by pointwise application ${\mathit{\varphi }}_N\left(X\right)={\left(\varphi \left(X_i\right)\right)}_i$ and ${\mathit{\psi }}_N\left(X\right)={\left(\psi \left(X_i\right)\right)}_i$.

We have the following explicit expressions for derivatives, adjoints, and gradients.

Theorem 1 (Derivatives computation).

Let $\mathbf{F},\mathbf{H}:{\mathbb{R}}^{N_x\times M}\to {\mathbb{R}}^{N_y\times B}$ and $\mathcal{D}$ be defined by (4), (7) and (8). The derivative of $\mathbf{F}$ and $\mathbf{H}$, as well as the adjoint and the gradient of $\mathcal{D}$ are given by:

$$\begin{array}{cc}\hfill {\mathbf{F}}^{\prime }\left[X\right]{\left(\xi \right)}^{⊺}& =-\mathbf{S}·(Q_X\odot \left(\mathbf{M}{\left(\mathbf{A}\xi \right)}^{⊺}\right)\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill {\mathbf{H}}^{\prime }\left[X\right]{\left(\xi \right)}^{⊺}& =-\mathbf{F}{\left(X\right)}^{-T}\odot (\mathbf{S}·(Q_X\odot \left(\mathbf{M}{\left(\mathbf{A}\xi \right)}^{⊺}\right)))\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill {\mathbf{F}}^{\prime }{\left[X\right]}^{*}\left(\eta \right)& =-{\mathbf{A}}^{⊺}·{({\mathbf{M}}^{⊺}·(Q_X\odot \left({\mathbf{S}}^{⊺}{\eta }^{⊺}\right)))}^{⊺}\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill {\mathbf{H}}^{\prime }{\left[X\right]}^{*}\left(\eta \right)& =-{\mathbf{A}}^{⊺}·{({\mathbf{M}}^{⊺}·(Q_X\odot \left({\mathbf{S}}^{⊺}(\mathbf{F}{\left(X\right)}^{-T}\odot {\eta }^{⊺})\right)))}^{⊺}\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill \nabla \mathcal{D}\left[X\right]& =-{\mathbf{A}}^{⊺}·{({\mathbf{M}}^{⊺}·(Q_X\odot \left({\mathbf{S}}^{⊺}(\mathbf{F}{\left(X\right)}^{-T}\odot {(\mathbf{H}\left(X\right)-Y^{\mathbf{H}})}^{⊺})\right)))}^{⊺},\hfill \end{array}$$

(15)

where $Q_X\triangleq {\mathbf{exp}}_{E\times N_y}(-\mathbf{M}{\left(\mathbf{A}X\right)}^{⊺})$ denotes virtual spectrally resolved data.

Proof. 

The proof is given in the Appendix A. □

For CP-fast, the derivative of $\mathbf{H}$ at zero plays a central role.

Remark 6 (Derivative at zero).

Let us consider the derivative at the zero image $X=0$. In this case, we have $Q_0:={\mathbf{exp}}_{E\times N_y}(-\mathbf{M}{\left(\mathbf{A}0\right)}^{⊺})={\mathbf{1}}_{E\times N_y}$, and therefore, ${\mathbf{H}}^{\prime }\left[0\right]{\left(\xi \right)}^{⊺}=-\left(\mathbf{S}\mathbf{M}{\left(\mathbf{A}\xi \right)}^{⊺}\right)\oslash \mathbf{F}{\left(0\right)}^{⊺}$ and ${\mathbf{H}}^{\prime }{\left[0\right]}^{*}\left(\eta \right)=-{\mathbf{A}}^{⊺}(\eta \oslash \mathbf{F}\left(0\right))·\mathbf{S}·\mathbf{M}$. Using that $\mathbf{F}{\left(0\right)}^{⊺}=\mathbf{S}·{\mathbf{1}}_{E\times N_y}$, we get:

$$\begin{array}{cc}\hfill {\mathbf{H}}^{\prime }\left[0\right]{\left(\xi \right)}^{⊺}& =-\mathbf{U}{\left(\mathbf{A}\xi \right)}^{⊺}\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill {\mathbf{H}}^{\prime }{\left[0\right]}^{*}\left(\eta \right)& =-{\mathbf{A}}^{⊺}{\eta }^{⊺}\mathbf{U}\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill \mathbf{U}& \triangleq \left(\mathbf{S}\mathbf{M}\right)\oslash (\mathbf{S}·{\mathbf{1}}_{E\times M}).\hfill \end{array}$$

(18)

The derivative ${\mathbf{H}}^{\prime }\left[0\right]\left(\xi \right)=-\mathbf{A}\xi {\mathbf{U}}^{⊺}$ may also be seen as the linearization of $\mathbf{H}$ around zero. It has been used previously in MSCT and can be simply derived by first order Taylor series approximation [12]. In fact, with $\mathbf{F}{\left(0\right)}_{\ell ,b}=\langle {\mathbf{S}}_b,{\mathbf{1}}_{E\times 1}\rangle$, we find:

$$\begin{array}{cc}\hfill log\left(\frac{\langle {\mathbf{S}}_b,exp(-\mathbf{M}{\left(\mathbf{A}\xi \right)}_{\ell }^{⊺})\rangle }{\langle {\mathbf{S}}_b,{\mathbf{1}}_{E\times 1}\rangle }\right)& \simeq log\left(\frac{\langle {\mathbf{S}}_b,(1-\mathbf{M}{\left(\mathbf{A}\xi \right)}_{\ell }^{⊺})\rangle }{\langle {\mathbf{S}}_b,{\mathbf{1}}_{E\times 1}\rangle }\right)\hfill \\ & =log\left(1-\frac{\langle {\mathbf{S}}_b,\mathbf{M}{\left(\mathbf{A}\xi \right)}_{\ell }^{⊺}\rangle }{\langle {\mathbf{S}}_b,{\mathbf{1}}_{E\times 1}\rangle }\right)\simeq -〈\frac{{\mathbf{S}}_b}{\langle {\mathbf{S}}_b,{\mathbf{1}}_{E\times 1}\rangle },\mathbf{M}{\left(\mathbf{A}\xi \right)}_{\ell }^{⊺}〉.\hfill \end{array}$$

(19)

The final expression in (19) in matrix notation is (16) and (18). For dual energy CT (the case where $M=B=2$), the use of the inverse of $\mathbf{U}$ has been proposed in [29]. We emphasize that while we utilize the linearization ${\mathbf{H}}^{\prime }\left[0\right]\left(\xi \right)=-\mathbf{A}\xi {\mathbf{U}}^{⊺}$ as an auxiliary tool, we actually solve the full nonlinear problem. However, the linearized LSQ problem ${∥\mathbf{A}\xi {\mathbf{U}}^{⊺}-Y_{\mathbf{H}}∥}_2^2/2\to {min}_{\xi }$ is also of interest in its own. Theoretically proven convergent algorithms for such problems can be found in [30,31].

The derivatives in Theorem 1 have a clear composite structure that we exploit in our algorithms. This is discussed next.

Remark 7 (Composite structure of derivatives).

Consider $X\in {\left({\mathbb{R}}^M\right)}^{N_x}$ as a signal of size $N_x$ with M channels (each channel is a material) and data $Y_{\mathbf{H}}\in {\left({\mathbb{R}}^B\right)}^{N_y}$ of size $N_y$ with B channels (each channel represents an energy bin). Then, we can write $\mathbf{H}\left(X\right)=\Phi \left(\mathbf{A}X\right)$, where the following nonlinear function:

$$\Phi \left(Z\right)={\mathbf{log}}_{N_y\times B}({\mathbf{exp}}_{N_y\times E}(-Z{\mathbf{M}}^{⊺})·{\underline{\mathbf{S}}}^{⊺})$$

(20)

operates on the multichannel sinogram $Z=\mathbf{A}X$ along the (horizontal) channel dimension and $\underline{\mathbf{S}}=\mathbf{S}\oslash (\mathbf{S}·{\mathbf{1}}_{E\times N_y})$ are the normalized effective spectra. Application of the chain rule and some computations results in:

$$\begin{array}{cc}\hfill {\mathbf{H}}^{\prime }\left[X\right]\left(\xi \right)& ={\mathbf{\Phi }}^{\prime }\left[\mathbf{A}X\right]\left(\mathbf{A}\xi \right)\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill {\mathbf{\Phi }}^{\prime }\left[Z\right]\left(\zeta \right)& =-({\mathbf{exp}}_{N_y\times E}(-Z{\mathbf{M}}^{⊺})\odot (\zeta ·{\mathbf{M}}^{⊺}))·{\mathbf{S}}^{⊺})\oslash \left({\mathbf{exp}}_{N_y\times E}(-Z{\mathbf{M}}^{⊺}){\mathbf{S}}^{⊺}\right)\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill {\mathbf{\Phi }}^{\prime }{\left[Z\right]}^{*}\left(\eta \right)& =-({\mathbf{exp}}_{N_y\times E}(-Z{\mathbf{M}}^{⊺})\odot ((\eta \oslash \left({\mathbf{exp}}_{N_y\times E}(-Z{\mathbf{M}}^{⊺}){\mathbf{S}}^{⊺}\right))·\mathbf{S}))·\mathbf{M}.\hfill \end{array}$$

(23)

from which we recover (12) and (14). Furthermore, for the zero material sinogram $Z=0$, we get ${\mathbf{\Phi }}^{\prime }\left[0\right]\left(\zeta \right)=-\zeta {\mathbf{U}}^{⊺}$ and ${\mathbf{\Phi }}^{\prime }{\left[0\right]}^{*}\left(\eta \right)=-\eta \mathbf{U}$ with $\mathbf{U}=(\mathbf{S}·\mathbf{M})\oslash (\mathbf{S}·{\mathbf{1}}_{N_y\times M})$ as in Remark 6. Equations (21)–(23) factorize the derivative ${\mathbf{\Phi }}^{\prime }\left[Z\right]$ and its adjoint into two separate parts: a high-dimensional ill-posed but linear part $\mathbf{A}$ operating in the pixel dimension and small size well-posed but nonlinear part operating in the channel dimension. Our algorithms will target this structure for fast and effective iterative updates.

3.2. Gradient and Newton—One-Step Algorithms

It is helpful to start with gradient based method for minimizing the LSQ functional (8). Our first method, CP-full, can be seen as a modified version of a hybrid between the standard gradient iteration (or Landweber’s method) and the Gauss–Newton method. Our second method, CP-fast, involves a simplification based on linearization around zero.

Gradient-based one-step algorithms for solving the MSCT problem $\mathbf{H}X=Y_{\mathbf{H}}$ using the LSQ functional $\mathcal{D}\left(X\right)$ start with the optimality condition $\nabla \mathcal{D}\left[X\right]={\mathbf{H}}^{\prime }{\left[X\right]}^{*}(\mathbf{H}\left(X\right)-Y^{\mathbf{H}})=0$ and fixed point equations derived from it. Applying a nonstationary positive-definite preconditioner ${\mathbf{Q}}_k$ and a step size ${\omega }_k>0$ results in:

$$X_{k+1}=X_k-{\omega }_k·{\mathbf{Q}}_k{\mathbf{H}}^{\prime }{\left[X_k\right]}^{*}\left(\mathbf{H}\left(X_k\right)-Y^{\mathbf{H}}\right).$$

(24)

Explicit expressions for $\mathbf{H}$ and ${\mathbf{H}}^{\prime }{\left[X_k\right]}^{*}$ are given by (7) and (14). Particular choices for the preconditioner and the step size yield various iterative solution methods. including Landweber’s iteration, the steepest descent method, Gauss–Newton iteration, Newton-CG iterations, or Quasi-Newton methods. To motivate our algorithm, it is most educational to discuss the Landweber and the Gauss–Newton iteration.

Landweber method:

In the context of inverse problems, the standard gradient method with a constant step size $\omega$ is known as the (nonlinear) Landweber iteration $X_{k+1}=X_k-\omega ·{\mathbf{H}}^{\prime }{\left[X_k\right]}^{*}(\mathbf{H}\left(X_k\right)-Y^{\mathbf{H}})$, which is (24) for the case where ${\mathbf{Q}}_k$ is the identity and ${\omega }_k=\omega$. Landweber’s iteration is stable, robust, and easy to implement. It is even applicable in ill-posed cases, where, with an appropriate stopping criterion, it serves as a regularization method [32]. On the other hand, it is also known to be slow in the sense that many iterative steps are required. In our case, this is due to the ill-conditioning of the forward operator.

Gauss–Newton method:

Several potential accelerations of Landweber’s method exist, and preconditioning seems one of the most natural ones. In the context of nonlinear least squares, the Gauss–Newton method and its variants are well-established and effective. In this case, one chooses the preconditioner ${\mathbf{Q}}_k={\left({\mathbf{H}}^{\prime }{\left[X_k\right]}^{*}{\mathbf{H}}^{\prime }\left[X_k\right]\right)}^{-1}$ in (24), which results in:

$$X_{k+1}=X_k-{\omega }_k·{\left({\mathbf{H}}^{\prime }{\left[X_k\right]}^{*}{\mathbf{H}}^{\prime }\left[X_k\right]\right)}^{-1}{\mathbf{H}}^{\prime }{\left[X_k\right]}^{*}\left(\mathbf{H}\left(X_k\right)-Y^{\mathbf{H}}\right).$$

(25)

The Gauss–Newton method (25) has the potential to significantly reduce the required number of iterations. On the other hand, each one of these iterations is numerically costly, as it requires inversion of the nonstationary normal operator ${\mathbf{H}}^{\prime }{\left[X_k\right]}^{*}{\mathbf{H}}^{\prime }\left[X_k\right]$. Moreover, due to ill conditioning, the inversion needs to be regularized [28,33,34]. The algorithms proposed in this paper use simplifications that do not need to be regularized and avoid the costly inversion of the normal operator.

3.3. Proposed Algorithms

Now we move on to the proposed iterative algorithms for MSCT. We start with CP-full, which is a gradient-based algorithm with channel preconditioning. We then derive CP-fast, which is a derivative-free iterative algorithm using a stationary adjoint problem.

CP-full: 

The first proposed algorithm is an instance of (24). Instead of no preconditioning, as in Landweber’s method, or the costly preconditioning ${\mathbf{Q}}_k={\left({\mathbf{H}}^{\prime }{\left[X_k\right]}^{*}{\mathbf{H}}^{\prime }\left[X_k\right]\right)}^{-1}$, as in the Gauss–Newton method, we propose preconditioning with the channel mixing term $\mathbf{\Phi }$ only. That is, we exploit the factorization $\mathbf{H}\left(X\right)=\mathbf{\Phi }\left(\mathbf{A}X\right)$ and propose the choice ${\mathbf{Q}}_k={\left({\mathbf{\Phi }}^{\prime }{\left[\mathbf{A}{X}_k\right]}^{*}{\mathbf{\Phi }}^{\prime }\left[\mathbf{A}{X}_k\right]\right)}^{-1}$ for the preconditioner. This results in the following CP-full iteration:

$$X_{k+1}=X_k-{\omega }_k·{\mathbf{A}}^{⊺}·({\mathbf{\Phi }}^{\prime }{\left[\mathbf{A}{X}_k\right]}^{*}{\left({\mathbf{\Phi }}^{\prime }\left[\mathbf{A}{X}_k\right]\right)}^{-1}{\mathbf{\Phi }}^{\prime }{\left[\mathbf{A}{X}_k\right]}^{*}(\mathbf{H}\left(X_k\right)-Y^{\mathbf{H}}).$$

(26)

While efficiently addressing the nonlinearity via a Gauss–Newton-type preconditioner in the channel dimension, it is computationally much less costly than the full Gauss–Newton update. Instead of inverting ${\mathbf{H}}^{\prime }{\left[X_k\right]}^{*}{\mathbf{H}}^{\prime }\left[X_k\right]$, which in matrix form has size $\left(N_{x}M\right)\times \left(N_{x}M\right)$ in the Gauss–Newton method, it requires inversion of the smaller $M\times M$ matrices ${\mathbf{\Phi }}^{\prime }{\left[\mathbf{A}{X}_k\right]}^{*}{\mathbf{\Phi }}^{\prime }\left[\mathbf{A}{X}_k\right]$ only, which can be done separately for each pixel in the projection domain. Assuming $M,B=\mathcal{O}\left(1\right)$ and $N_y=\mathcal{O}\left(N_x\right)$, this dramatically reduces the cost of preconditioning from $\mathcal{O}\left(N_x^3\right)$ to $\mathcal{O}\left(N_x\right)$ per iterative update.

CP-fast: 

In the derivative-free version, we go one step further and completely avoid the derivative ${\mathbf{H}}^{\prime }\left[X_k\right]$. For that purpose, we replace the derivative ${\mathbf{\Phi }}^{\prime }\left[\mathbf{A}{X}_k\right]$ in (26) by the derivative at zero. According to Remark 6, we have ${\mathbf{\Phi }}^{\prime }\left[0\right]\left(\zeta \right)=-\zeta {\mathbf{U}}^{⊺}$ with $\mathbf{U}=(\mathbf{S}·\mathbf{M})\oslash (\mathbf{S}·{\mathbf{1}}_{E\times M})$. Now, with ${\mathbf{U}}^{‡}={\left({\mathbf{U}}^{⊺}\mathbf{U}\right)}^{-1}{\mathbf{U}}^{⊺}$ denoting the pseudoinverse of $\mathbf{U}$, we arrive at the iterative update:

$$X_{k+1}=X_k-{\omega }_k·{\mathbf{A}}^{⊺}·(\mathbf{H}\left(X_k\right)-Y^{\mathbf{H}})·{\left({\mathbf{U}}^{‡}\right)}^{⊺}.$$

(27)

We refer to (27) as the derivative-free fast channel-preconditioned (cp-fast) iteration. It only involves the derivative at zero, which can be computed once before the actual iteration. In this sense, it is actually derivative-free and fast. It can be interpreted as using the full nonlinear model for the forward problem, the linearization at zero for the adjoint problem, and including channel preconditioning.

Both iterations (26) and (27) are of fixed-point type and we, therefore, expect convergence for sufficiently small step sizes. Theoretically, proving convergence seems possible, but this is beyond the scope of this paper. As (26) is of gradient type, it seems easier to derive convergence for CP-full, while for the derivative free version CP-fast, such a proof seems challenging. Note further that for the results presented below, we integrated a positivity constraint by alternating iterative updates with the orthogonal projection onto the cone of nonnegative images.

4. Numerical Simulations

We compared our algorithms CP-full and CP-fast to existing iterative one-step algorithms in MSCT. Our evaluation builds on [11], which compares five such algorithms and provides open source code (https://github.com/SimonRit/OneStepSpectralCT, accessed on 13 July 2022.) that is used for our results. We compare CP-full and CP-fast with the best performing one of [11], and further with a two-step method.

4.1. Comparison Methods

The work [11] compares the following iterative one-step algorithms for MSCT in terms of memory usage and convergence speed to reach a fixed image quality threshold:

• Ref. [19] derives a nonlinear CG method for a weighted LSQ term;
• Refs. [20,22] derive surrogate approaches for Poisson maximum likelihood;
• Ref. [21] extends [22] by including Nesterov’s momentum acceleration;
• Ref. [23] generalizes the Chambolle–Pock algorithm [35] to nonconvex functionals.

Specifically, ref. [11] found the algorithm of [21] (referred to as Mechlem2018) to be significantly faster than the other four methods, and thus, we use it for comparison.

In addition, we compare with the algorithm [17] (referred to as Niu2014) as a prime example of an image domain two-step method. They use a penalized weighted least squares estimation technique applied to an empirical linear model. Note that more recently, data-driven methods based on neural networks and deep learning have also been proposed. Such methods are beyond the scope of this manuscript and we refer the interested reader to the review articles [10,36].

4.2. Numerical Implementation

For the presented results, we build on the Matlab code of [11], which we extend with our algorithms. In particular, we work with $M=3$ base materials (water, iodine, and gadolinium) and $B=5$ energy bins. The energy variable is discretized using $E=150$ uniform nodes between 0 and 150 keV. The attenuation functions and energy spectra used are shown in Figure 3. We use $N_x=256\times 256$ image pixels and $N_y=262450$ line integrals for the Radon transform. In particular the code https://github.com/SimonRit/OneStepSpectralCT (accessed on 13 July 2022) creates matrices:

• $\mathbf{M}\in {\mathbb{R}}^{150\times 3}$ for the base materials;
• $\mathbf{S}\in {\mathbb{R}}^{5\times 150}$ for the effective energy spectra;
• $\mathbf{A}\in {\mathbb{R}}^{N_y\times N_x}$ for the Radon transform.

After row normalizing $\underline{\mathbf{S}}=\mathbf{S}\oslash (\mathbf{S}·{\mathbf{1}}_{E\times N_y})$, we have $\mathbf{H}\left(X\right)=log(exp(-\mathbf{A}·X·{\mathbf{M}}^{⊺})·{\underline{\mathbf{S}}}^{⊺})$ for the MSCT forward model. Furthermore, noisy data Y are created with a different realistic forward model and Poisson noise added.

Besides $\mathbf{M}$, $\mathbf{S}$, $\mathbf{A}$, $\mathbf{H}$, and Y, we require implementations of ${\mathbf{U}}^{‡}={\left({\mathbf{U}}^{*}\mathbf{U}\right)}^{-1}{\mathbf{U}}^{*}$ for CP-fast and ${\mathbf{U}}_k^{‡}={\left({\mathbf{\Phi }}^{\prime }{\left[\mathbf{A}{X}_k\right]}^{*}{\mathbf{\Phi }}^{\prime }\left[\mathbf{A}{X}_k\right]\right)}^{-1}{\mathbf{\Phi }}^{\prime }{\left[\mathbf{A}{X}_k\right]}^{*}$ for CP-full. Computing ${\mathbf{U}}^{‡}$ is trivial and can be done in advance; ${\mathbf{U}}_k^{‡}$ are computed in each step of CP-fast using (22) and (23).

4.3. Results

Reconstruction results using the proposed algorithms CP-fast (top row) and CP-full (second row) and the two comparison methods Mechlem2018 [21] (row three) and Niu2014 [17] (bottom row) can be seen in Figure 4. The phantom shown in the top row is made out of iodine (left), gadolinium (middle), and water (right). In all cases, we use noisy data and plot the iteration with minimal ${\ell }^2$-reconstruction error ${∥X_k-X_{★}∥}_2/{∥X_{★}∥}_2$, where $X★$ is the ground truth. All calculations are performed on the same standard laptop, where one iteration of CP-full takes around six seconds, one iteration of CP-fast takes around one second, and one iteration of Mechlem2018 [21] about four seconds. These times are comparable to Niu2014 [17], which takes around 5 to 8 s when using 500 and 1000 CG iterations. Figure 5 shows the evolution of the relative ${\ell }^2$-reconstruction error for various one-step methods. Note that for CP-fast, the minimum error in Iodine and Gadolinium is reached at approximately the same number of iterations, which shows efficient preconditioning and is important in application. Furthermore, we do not enforce that the sum over the three density images is one, and in the example, it indeed does not hold. The proposed algorithms turned out to be more stable than Mechlem2018 [21] and produce better results. In particular, CP-full gives the best results, while CP-fast is fastest.

Additional results for a different phantom are shown in Figure 6 and Figure 7. The concentration of iodine, gadolinium, and water has been computed as in [11].

5. Conclusions and Outlook

Image reconstruction in MSCT requires the solution of a nonlinear ill-posed problem. Iterative one-step methods are known to be accurate for this purpose. In this work, we propose two generic algorithms named CP-full (channel-preconditioned full gradient iteration) and CP-fast (channel-preconditioned fast iteration). Both algorithms use preconditioning in the channel dimen•sion only, which considerably accelerates the updates compared to Newton-type methods that require solving numerically costly linear inverse problems at each iteration. CP-fast replaces the derivative in the channel nonlinearity with linearization at zero, making it even more efficient. Both algorithms turn out to be fast and robust.

There are several future directions emerging from our work. First, proving the convergence of the two algorithms and demonstrating their regularization properties is important. Second, we will combine them with more realistic noise priors, such as Poisson noise, resulting in the maximum likelihood estimation (MLE) functional. Additionally, we will integrate explicit image priors, use plug-and-play strategies, and incorporate learned components.