{'fraction_non_alphanumeric': 0.030072463768115943, 'fraction_numerical': 0.0029790660225442835, 'mean_word_length': 4.951365596550072, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Preconditioned Algorithm for Model-Based Iterative CT Reconstruction and Material Decomposition from Spectral CT Data

Model-based material decomposition is a statisticaliterative reconstruction framework where basis material densityimages are estimated directly from spectral CT data. This methoduses a physical model for polyenergetic x-ray transmission andattenuation and therefore it does not typically suffer frombeam-hardening artifacts. However, this estimation is a poorly-conditioned inverse problem due to the strong anticorrelationbetween basis materials. In this work we propose an precondi-tioned optimization algorithm for a nonlinear penalized weightedleast-squares objective function.


I. INTRODUCTION
Model-based iterative reconstruction (MBIR) was proposed for x-ray computed tomography (CT) over 20 years ago [4] [3], yet they are only recently starting to be implemented on commercial CT systems due, in large part, to the computational cost [7]. These methods are generally based on maximum a posteriori (MAP) estimation or an approximation thereof (e.g. penalized likelihood). Through a combination of advanced physical models, statistical weighting of the data, and the incorporation of prior knowledge via regularization terms, MBIR has consistently demonstrated the potential for improved image quality for the same radiation dose, or the same image quality for a reduced radiation dose. However, to achieve these performance improvements, the algorithm generally involves running many iterations of an image update routine to numerically optimize the MAP objective.
One such optimization algorithm is the separable parabolic surrogates (SPS) algorithm, originally proposed by [4]. It was recently extended an advanced physical and statistical model for the realistic non-idealities of data acquisition with modern multi-slice CT or cone-beam CT systems [9]. By incorporating effects such as system blur (e.g. focal spot, scintillating detectors), noise correlations, and non-linear polyenergetic effects (e.g. beam hardening) into the reconstruction model, MBIR can result in an inversion of some of these effects (e.g. focal spot deconvolution) to reconstruct high-fidelity images.
This model can also be extended to material density quantization. If the polyenergetic attenuation from the object is parameterized by a finite number of basis material density images, and the projection measurements have varied sensitivity spectra, as in spectral CT systems, those density images can be estimated in a one-step model-based material decomposition (MBMD) algorithm [10] [11].
Due to the large size of the data from multi-slice or conebeam CT scanners, the various complex geometric and physical models used by MBIR, and the large number of iterations necessary to optimize the MAP objective function, the computational cost of these algorithms is high. For this reason, it is an ongoing area of interest to investigate methods to accelerate the optimization routine. Gradient-based iterative optimization algorithms typically converge slowly for large-scale or poorlyconditioned inverse problems such as the MAP estimator used in MBIR. For that reason, preconditioning techniques can have a dramatically positive impact on convergence.
The simplest preconditioner is a diagonal scaling matrix. Several algorithms including Expectation Maximization (EM) [6], Scaled Gradient Descent (SGD), [4], and SPS itself can be viewed as diagonal preconditioners applied to the standard gradient descent algorithm. While these diagonal preconditioners are effective and reliable, they are considered to be relatively conservative approximations of the ideal preconditioner which would be the inverse of the hessian of the objective function.
Non-diagonal preconditioners based on a Fourier basis have have also been applied to a conjugate-gradient-based reconstruction algorithm assuming a linear model. These include Fourier preconditioners which are closely connected to the ramp filter [1] and an extension to shift-variant systems using channelized frequency bands [5]. For MBMD, the conditioning is even worse than conventional CT image reconstruction with MBIR due to the strong correlations between different estimated basis materials. A block-diagonal preconditioner for the primal-dual optimization algorithm was proposed by [8]. This strategy dramatically improves the conditioning of the material decomposition portion of the problem.
In this work we propose a general preconditioned version of the SPS algorithm. We define the conditions for the optimal preconditioner and describe a process for finding approximations thereof. We also apply this general preconditioned SPS algorithm to three specific cases: conventional image reconstruction from single-energy CT data, projection-domain decomposition from spectral CT data, and one-step basis material density estimation from spectral CT data in a full preconditioned MBMD algorithm. For the single-energy CT case, we derive a shift-invariant preconditioner based on a fourier basis approximation of the hessian, as well as a shiftvariant version based on a haar-wavelet basis approximation. For the projection-domain decomposition case, we describe how to same model can be used to estimate material line integrals (without spatial reconstruction) and we propose a cross-material preconditioner. Finally, for the one-step MBMD case, we propose a block-diagonal cross-material preconditioner, and a cross-material cross-voxel preconditioner based on the fourier basis for each material, as well as a shift-variant version using a haar-wavelet basis.

A. Generalized Models for Data Acquisition and Estimation
For a generalized physical model for x-ray transmission data acquisition, we assume that the projection data, y, are normally distributed with known covariance, Σ y , and mean, y(x) given by (1).ȳ where, x, is an unknown parameter vector, and, A A A and B B B are placeholders for linear operators inside and outside the exponential operator, respectively. They are used to define the expected non-linear relationship between x and the projection measurements, y. This general formulation can be used to model physical phenomena such as system blur or beam hardening for single-energy CT systems. We will show that, based on the definition of x, this model can also be used for projection-domain decomposition or direct estimation of material density distributions via model-based material decomposition from spectral CT data.
A maximum a posteriori estimator for this statistical model can be approximated in terms of the negative log-likelihood of the data L(y|x) and an approximation of the negative logprior probability of x, given by the penalty R(x) as shown in the penalized likelihood objective function shown in (2).
For the purposes of this work we will assume the penalty function takes the quadratic form, R(x) = 1 2 x T Rx. The objective function, Φ(x), is globally convex, but it has no closed-form analytical solution for the global minimum. Instead, estimatex with an iterative optimization algorithm.

B. The Separable Parabolic Surrogates Algorithm
The concept of operation for the SPS algorithm involves defining a surrogate objective function,Φ(x; x (n) ) which is separable (i.e. diagonal Hessian) with respect to the elements of x, and majorizes the original objective, Φ(x). Therefore, we can iteratively update x (n+1) as the minimum of the surrogate function for which a closed-form solution exists. The full algorithm is given by the following set of formulas.
whereL(x; x (n) ) andR(x) are separable surrogates for the data likelihood and regularization terms of the objective function, respectively. The quantities η, γ, ρ (n) , and c (n) are defined below.
where f c (x (n) ) is the maximum curvature calculation described in [9]. While this algorithm guarantees that the cost will monotonically decrease with each iteration, the separable approximation of the objective function will be very poor if the estimates, x, are strongly inter-correlated, leading to a small step size and slow convergence. Therefore, we propose to derive a preconditioned version of the algorithm to find a new basis for estimation for which separability is a more effective approximation.

C. Preconditioned Optimization
The Hessian of the objective function evaluated at the solution,x, is given by For the penalized likelihood objective in (2), this hessian is Note that for many cases, an effective approximation for D is available without dependence onx. In many cases, there are strong cross-estimate correlations that make this a poorly-conditioned inverse problem to estimatex given y. We propose to find a linear transformation to a new basis, x ′ , which improves the conditioning of the estimation problem. This transformation is defined by the relationship, x = M x ′ . Therefore, the new hessian, H ′ , associated with the preconditioned estimation ofx ′ is Therefore, effective preconditioner will satisfy, M T HM ≈ I, or at least the condition number of M T HM will be much less than that of H.
Substituting Mx ′ into (1) gives Based on this form we propose to use a preconditioned version of the SPS algorithm. This involves replacing the linear operator, A A A, with the preconditioned version, A A AM, estimating x ′ via the SPS algorithm, and finally applying x = M x ′ to reproject the estimates to the original basis. Since H is symmetric by construction and assumed to be full-rank, we will aim to find a precondtitioner which approximates the unique symmetric matrix square root of the inverse Hessian (i.e. M = H − 1 2 ). The eigen decomposition of H is given by where U H is an orthogonal matrix containing the eigenvectors of the hessian, H, and Λ H is a diagonal matrix containing the eigenvalues. Therefore, the ideal preconditioner is given by In the following sections we will propose preconditioner designs for a few specific applications by finding an orthogonal basisŨ H which is anapproximation of the eigenbasis of the hessian. Then we will find the corresponding approximate eigenvalues,Λ H according to the following formula.
which will result in an approximation of the hessian given byH =Ũ T HΛ HŨH which matches the response of H for some test function, δx test . Note that δx test should be designed appropriately to excite all of the eigenvectors such that D{Ũ H ∆x test }] −1 is defined. Finally, we can substitute the approximationsŨ H andΛ H into (17) to establish a establish a preconditioned optimization algorithm.

D. Standard CT Reconstruction
For standard CT reconstruction, we aim to estimate an image composed of attenuation coefficients, µ, given standard single-energy CT data, y. Therefore, we have the following definitions.
where A is the forward projector, capturing the system's spatial sampling geometry, B models the projection-domain blur (e.g. focal spot, scintilating detectors), and G is the gain. Note that this particular model does not include a polyenergetic absorption model and would therefore be subject to beamhardening artifacts if applied to measurements acquired with a polyenergetic source. This leads to the following definition of the Hessian.
An effective approximation for D is available directly from the gain-corrected measured data. The vector µ is sized N j × 1 and therefore the Fisher information matrix F and the quadratic regularization matrix R are sized N j × N j . The following is a quadratic smoothness regularizer, for example.
where N j is a neighborhood around voxel j (not including j) and |N j | is the cardinality of that neighborhood.

1) Fourier Preconditioner:
To begin the derivation of a cross-voxel preconditioner, consider the case where A represents the discrete radon transform, corresponding to a parallel beam imaging geometry, and D, B, G, and Σ y are all diagonal. For that special case, F −1 = (A T A) −1 is the ramp filter used in filtered backprojection, so F takes the form shown below.
where U DFT is the unitary discrete fourier transform, which is the eigenbasis of F, and Λ F is diagonal and contains the eigenvalues, or the frequency coefficients for the inverse ramp filter. This form implies that F is a shift-invariant operator.
For the example in (26), the matrix operator R is also shiftinvariant and can therefore also be diagonalized by the discrete fourier transform. Therefore, the eigenbasis of H is also the a fourier basis or U H = U DFT , and we can solve for Λ H exactly using the formula by applying (18) as shown below.
where δ j is an impulse with a small value only at voxel j.
Although we understand that for the general case where A may represent a non-parallel geometry or for non-diagonal matrices D, B, G, and Σ y , the response of F will be shiftinvariant, we can still make the shift-invariant approximation U H ≈ U DFT and apply (28) to establish an approximate set of eigenvalues (or filter coefficients) for a shift-invariant preconditioner. The result will be an approximate model for the hessian that perfectly matches the frequency response of H at position j, but it may be a poor approximation for positions j ′ which are far from j. This Fourier preconditioner is therefore constructed as follows.
This is a desirable preconditioner because it is essentially a linear filter which can be applied using fast fourier transform algorithm, a diagonal scaling matrix, and the inverse fourier transform. This is very similar to the preconditioner proposed in [1]. The pitfall of this approach is that for cases where the operator H is strongly shift-invariant, this approximation of the hessian may be very poor.
2) Modified Fourier Preconditioner: The Fourier preconditioner described above assumes the same frequency response at every position in the image domain. One straight-forward modification to this algorithm is to first normalize the zero-frequency response before assuming shift invariance. This modification extends the model to an operator with relative frequency weights which are shift invariant with a shift-variant overall scale.
The zero-frequency response normalization is given by the following diagonal preconditioner.
Then, we assume the normalized hessian is a Fourier system.
and solve for the eigenvalues accordingly.
The modified Fourier preconditioner is therefore given by the following formula.
This combination of a diagonal preconditioner and a Fourier preconditioner is capable of approximating the shift-variant scale but still assumes a shift-invariant relative frequency response.

3) Wavelet Precondtioner:
For a shift-variant preconditioner, we propose to approximate the eigenbasis of H with a unitary discrete wavelet transform, that is, U H ≈ U DWT . Examples of unitary wavelet transforms include the Haar transform and other unitary transforms based on Daubechies wavelets [2].
We have chosen this particular basis because wavelets are capable of encoding information about both frequency response and shift-variant behaviors.
To find the eigenvalues, we use a test function δx test = U DWT T 1 which is excites all eigenvectors of the wavelet basis . This is analogous to the impulse function, δ 1 = U DFT T 1, which excites all eigenvectors of the Fourier basis.
Following a similar procedure as described in the previous section, the zero-frequency response is normalized via the diagonal preconditioner shown in (30). Then, assuming a wavelet eigenbasis, the approximate eigenvalues of the normalized hessian, M 0 TΛ H M 0 , can be established using (18) and the preconditioner can be formulated as In general, the pitfalls of the wavelet transform are that there are relatively few number of eigenvectors used to describe the shift-variance of the low-frequency response, and there is also an upper limit on frequencies that are captured.
This combination of a diagonal preconditioner and a wavelet preconditioner is capable of modeling a shift-variant local frequency response.

E. Projection-Domain Decomposition
For model-based projection-domain material decomposition, we aim to estimate material line integrals, ℓ, given spectral CT data, y which is composed of projections with varied spectral sensitivity arranged into channels. Typically there are at least as many channels as materials. Therefore, we have the following definitions.
where Q maps material basis-material line integrals to projection-domain attenuation spectra using a weighted sum of the attenuation spectra for each basis material, S is the projection-dependent system spectral sensitivity for each channel of the spectral CT system, B models the channeldependent projection-domain blur, and G is the gain. This leads to the following definition of the Hessian.
An effective approximation for D can be obtained by assuming all of the attenuation was due to water and conducting an inexpensive single-material line integral estimation. The vector ℓ is sized N i N k × 1, so the hessian is sized Here i is used to index projections and k is used to index materials.

1) Cross-Material Preconditioner:
Since the estimation of ℓ occurs in the projection-domain, different positions can be effectively approximated as separable. That is, the elements of the cross-projection secondderivatives in the hessian can be approximated as zero, or H ik,i ′ k ≈ 0 for i = i ′ . Therefore, an approximate preconditioner can be constructed from N k ×N k blocks, denoted asH i for each projection. There are N i N k values which compose column k of those blocks, denoted as the N i N k × 1 vector H k , which can be computed as where the vector 1 k has a flat response of ones for material k and zeros for other materials. After N k applications of (42), the block-diagonal hessian can be reorganized into the blocks,H i , and the preconditioner for can be constructed as M i =H i . This leads to a blockdiagonal preconditioner M which assumes separability across projections, but applies a non-diagonal preconditioning accross materials.

F. One-Step Reconstruction and Material Decomposition
In direct model-based material decomposition, we aim to estimate material density maps, ρ, given spectral CT data, y which is composed of projections with varied spectral sensitivity. Therefore, we have the following definitions.
where A is the forward projector, capturing the system's spatial sampling geometry, Q contains the mass attenuation spectra for each basis material, S is the projection-dependent system spectral sensitivity, B models the projection-domain blur (e.g. focal spot, scintilating detectors), and G is the gain. This leads to the following definition of the Hessian.
We assume that R is the hessian of a cross-material quadratic smoothness penalty, R(ρ), which we define as (49) where ρ jk is the denisty of basis material k at voxel j, the set N j is the set of voxels neighboring j (not including j) which has cardinality |N j |, and the parameters, β k,k ′ > 0 are the cross-material regularization weights. We can write R in terms of the two symmetric matrices β, modeling the regularization strength, and Γ modeling the kernel for the spatial smoothness penalty as follows.

1) Cross-Material Precondtioner:
The correlations between estimates for different materials severely impacts the the conditioning of the estimation problem for model-based material decomposition. Therefore we will defineH as shown below in order to define a shift-variant cross-material preconditioner.
Therefore, H can be constructed from N j blocks denoted byH j which are sized N k × N k . Column j of all of those blocks can be organized into a N k N j × 1 vector, denoted bỹ H k , which can be computed as where 1 k is an image of ones for material k and zero for the other materials. Then, the hessian can be reorganized into N j blocks,H j . Then, a block-diagonal preconditioner, M 0 , which has a block for each voxel equal toH − 1 2 j . This block-diagonal preconditioner matches the zerofrequency response of the hessian including the cross-material zero-frequency response. However, this format does not handle the cross-voxel correlations in the estimates or the more abstract cross-voxel-cross-material correlations.
2) Cross-Material Fourier Precondtioner: In this section we seek to derive a cross-material Fourier preconditioner to handle both cross-material and cross-voxel correlations. We begin with the approximation that after normalizing the zero-frequency response, the modified hessian operator is approximately shift-invariant. Note that this shiftinvariance describes both the in-basis response as well as the cross-basis response of the normalized hessian. We express this in the formula below.
U DFT represents the unitary discrete Fourier transform applied individually and identically to all materials. The matrix L contains weights applied to spatial the N j spatial frequencies.
If all spatial frequencies follow the same cross-material relationship, then the spatial correlations and material correlations would be entirely separable and L would be a diagonal matrix. In general, this may not be the case so we assume L is block diagonal contain N j blocks, one fore each spatial frequency, denoted as L j , each sized N k × N k . We can compute the column k for all N j blocks organized into the N j N k × 1 vector, L k , as follows.
The formula above describes the Fourier transform of the zerofrequency-normalized hessian (M 0 T HM 0 ) applied to δ jk , an impulse in material k at position j, divided by the frequency response of δ jk to correct for the phase. Therefore, the vector L k represents the local response of all materials and spatial frequencies to an impulse in material k at position j.
After rearranging the columns, L k , into the N j blocks, L j , the cross-material Fourier preconditioner can be constructed by taking the symmetric matrix negative square root of each of the blocks as written below.
This is the cross-material analog of the modified Fourier preconditioner defined in a previous section for standard CT reconstructions. It makes the assumption that after normalization of the zero-frequency response (including cross-material), the modified hessian is approximately shift-invariant. This form will perfectly match the local impulse response (including cross-material) for at least one position.

3) Cross-Material Wavelet Precondtioner:
Instead of assuming a Fourier basis for each material as in the previous section, we can assume a wavelet basis of the zero-frequency-normalized hessian as shown below.
we can solve for the block-diagonal matrix K sequentially for each material as formulated below.
where U DWT represents a unitary discrete wavelet transform (such as the Haar transform) applied individually and identically for each material, and the test function ∆x is designed to excite all eigenvectors of the wavelet basis just as an impulse function was used to excite all eigenvectors of a Fourier basis. After rearranging the block-diagonal matrix L, the crossmaterial wavelet preconditioner can be written as This form is capable of modeling complicated crossmaterial and cross-voxel correlations with relatively few total paramters. The cross-material preconditioner M 0 inverts the zero-frequency response of the hessian accross all positions, and the wavelet-domain operator paramterized by the block diagonal matrix K can be described as a model for the residual shift-variant cross-material spatial frequency response.

III. CONCLUSION
The slow computational speed of MBIR is a major factor preventing widespread implementation of clinical systems. The SPS algorithm is able to monotonically decrease the objective function in an iterative optimization algorithm. However, if the estimates are highly correlated, the separable surrogate objective function will be a poor match to the truth, leading to very small step sizes.
Effective preconditioners are an efficient way to improve the conditioning of the inverse problem. In this work, we have proposed a preconditioned version of the SPS algorithm. We have also derived specific preconditioners for applications including standard CT reconstruction, projection-domain material decomposition, and one-step direct model-based material decomposition and image reconstruction. Some of the nondiagonal preconditioners presented in this paper (e.g. Fourier preconditioner) have been previously proposed in a similar form. We have also presented novel preconditioning strategies such as the wavelet preconditioner and cross-material preconditioners.
With effective preconditioning, there is a potential to the computation time necessary for MBIR. These preconditioned optimization algorithms could accelerate reconstructions of standard CT and spectral CT data in a way that makes MBIR more practical for clinical implementation on commercial systems.