PaperlayerView exact revision on arXiv

Source-linked AI summary

An algorithm for constrained one-step inversion of spectral CT data

Rina Foygel Barber, Emil Y. Sidky, Taly Gilat Schmidt, Xiaochuan Pan

arXiv:1511.03384v1physics.med-ph

TL;DR

Spectral CT reconstruction traditionally estimates basis-map sinograms before image reconstruction, while one-step reconstruction with basis-map constraints remains challenging. This paper develops a constrained primal-dual one-step algorithm and demonstrates it in simulations, including TV-constrained data.

  • Problem

    Spectral CT workflows estimate basis-map sinograms before image reconstruction, leaving one-step reconstruction as a central challenge.

  • Method

    The paper develops a primal-dual framework for directly reconstructing basis maps from spectral CT data while enforcing non-smooth convex constraints.

  • Results

    Simulations demonstrate the one-step algorithm for TV-constrained spectral CT data.

  • Takeaways & Limitations

    The framework supports applying non-smooth convex constraints to basis maps during one-step spectral CT reconstruction.

Abstract

from arXiv · show

We develop a primal-dual algorithm that allows for one-step inversion of spectral CT transmission photon counts data to a basis map decomposition. The algorithm allows for image constraints to be enforced on the basis maps during the inversion. The derivation of the algorithm makes use of a local upper bounding quadratic approximation to generate descent steps for non-convex spectral CT data discrepancy terms, combined with a new convex-concave optimization algorithm. Convergence of the algorithm is demonstrated on simulated spectral CT data. Simulations with noise and anthropomorphic phantoms show examples of how to employ the constrained one-step algorithm for spectral CT data.

I. INTRODUCTION

Spectral CT enables quantitative imaging from energy-windowed photon measurements, but existing processing choices trade simplicity against model fidelity and artifact suppression. This work develops a constrained one-step reconstruction framework that directly fits the spectral model while supporting convex constraints on basis maps and composite images.

  • Photon-counting spectral CT can reconstruct quantitative attenuation and material information from energy-windowed transmission measurements, with potential improvements in beam-hardening, CNR, SNR, and quantitative imaging.With three or more energy windows, the approach can also image contrast agents with diagnostic-range K-edges.
  • Pre-reconstruction processing avoids beam-hardening by estimating consistent basis-map sinograms before reconstruction, but requires spectral calibration and registered projections.
  • Post-reconstruction processing is simpler and needs no projection registration, but negative-logarithm data can produce energy-window images susceptible to beam-hardening artifacts.
  • The one-step approach directly estimates basis maps from energy-windowed transmission data, treating the spectral model exactly while incorporating basis-map constraints, but requiring iterative optimization.
  • The proposed framework applies non-smooth convex constraints, including TV, non-negativity, upper bounds, and sum bounds, while optimizing least-squares or transmission-Poisson data discrepancies.The paper develops a convex-concave primal-dual algorithm for the non-convex discrepancy arising from the nonlinear spectral CT model and demonstrates it on simulated data.

II. ONE-STEP IMAGE RECONSTRUCTION FOR SPECTRAL CT · A. Spectral CT data model · B. Constrained optimization for one-step basis decomposition

The paper models spectral CT transmission counts with a low-dimensional basis-material representation of energy-dependent attenuation, then formulates one-step decomposition as constrained optimization. The framework accommodates non-convex data discrepancies and convex physical or image constraints on basis and monochromatic maps.

  • A. Spectral CT data model: Energy-dependent attenuation is represented with a low-dimensional basis-material expansion, with fractional density maps as the reconstruction targets.The expansion uses material mass attenuation coefficients and spatially varying fractional density maps.
  • A. Spectral CT data model: The discretized model uses pixel or voxel material maps and a line-intersection matrix to represent ray integration.Matrix entries measure the intersection length between each transmission ray and pixel or voxel.
  • A. Spectral CT data model: Photon-counting measurements are modeled through incident and transmitted spectral fluence, with transmitted fluence interpreted as an expected value.This interpretation resolves the non-integer value produced by the forward model despite photon counts being integer-valued observations.
  • B. Constrained optimization for one-step basis decomposition: One-step basis decomposition minimizes a non-convex data-discrepancy objective subject to convex constraints encoded by indicator functions.The constraints can represent requirements such as nonnegativity of material maps.
  • B. Constrained optimization for one-step basis decomposition: The framework considers transmission Poisson likelihood and least-squares discrepancies, which weight measurements differently and are non-convex in the basis maps.Least squares treats measurements equally, whereas transmission Poisson likelihood gives greater weight to higher-count measurements.
  • B. Constrained optimization for one-step basis decomposition: Convex constraints can enforce material-map bounds, sum-of-fractions restrictions, monochromatic-image nonnegativity, and total-variation limits.Total-variation constraints promote gradient-magnitude sparsity in basis maps or monochromatic images, while other convex constraints can also be incorporated.
  • B. Constrained optimization for one-step basis decomposition: Bounds and sum constraints may need loosening when computational basis materials differ from those present in the scanned object, reducing their effectiveness.The mismatch arises because the assumed basis may not represent all tissue or density combinations in the subject.

III. A FIRST-ORDER ALGORITHM FOR SPECTRAL CT CONSTRAINED OPTIMIZATION

The section generalizes Chambolle–Pock’s primal-dual algorithm to constrained spectral CT optimization with non-convex data-discrepancy terms. It uses local convex quadratic upper bounds and a mirrored convex-concave formulation to generate descent steps for the reconstruction algorithm.

  • Algorithm development: The main contribution generalizes CP’s primal-dual algorithm to minimize the spectral CT data-discrepancy term while handling convex image constraints.The data-fidelity derivation focuses specifically on DTPL and DLSQ.
  • Algorithm development: The algorithm uses nested outer and inner loops: each outer iteration builds a local convex quadratic upper bound, while the inner loop takes descent steps for that bound.This construction targets non-convex data-discrepancy terms.
  • Implementation: In practice, one inner-loop iteration is used, so each iteration re-expands the data-discrepancy term and performs a single effective update loop.The formal nested-loop construction is therefore implemented as a single loop.
  • Convex-concave optimization: Because the local quadratic upper bound does not fit CP’s generic form, the method uses a convex-concave generalization called the mirrored convex-concave algorithm, or MOCCA.The section applies MOCCA to one-step spectral CT reconstruction along with preconditioning and convergence checks.

A. A local convex quadratic upper bound to the spectral CT data discrepancy terms · 1) Quadratic expansion:

The method forms local quadratic approximations to spectral CT data discrepancies using their gradients and Hessians, then enforces convexity by retaining positive-semidefinite curvature components. The resulting convex function Q is combined with convex constraints to generate descent steps for the non-convex optimization problem, although it bounds only the quadratic expansion rather than necessarily the discrepancy itself.

  • 1) Quadratic expansion:: The algorithm derives local quadratic approximations of the spectral CT data discrepancy terms around the material maps using explicit gradients and Hessians.The derivations are carried out for both DLSQ and DTPL, with gradients expressed through count or log-count residuals.
  • 1) Quadratic expansion:: The resulting Taylor quadratic is generally non-convex because the Hessian expressions can contain negative curvature.This motivates modifying the Hessian before using the quadratic in optimization.
  • 1) Quadratic expansion:: The method splits each Hessian into positive and negative components and sets the negative components to zero, producing non-negative curvature.The split is formulated so the resulting Hessian components are positive semidefinite.
  • 1) Quadratic expansion:: A tighter positive-semidefinite Hessian could retain only non-negative eigenvalue directions, but eigenvalue decomposition is impractical for realistic tomography configurations.The implemented algebraic split avoids this computational burden.
  • 1) Quadratic expansion:: The convex local quadratic Q uses the original discrepancy value and gradient at expansion center f0, with the positive Hessian component in its quadratic term.Q depends on the counts data c through L and on the expansion center f0.
  • 1) Quadratic expansion:: The constructed Q is not necessarily an upper bound for the data discrepancy, even locally, because only the quadratic expansion is bounded.The method uses Q as a convex surrogate rather than asserting a bound on the original discrepancy functions.
  • 1) Quadratic expansion:: Combining Q with convex constraints generates descent steps for the generic non-convex optimization problem.This supplies the optimization step used after constructing the convex local approximation.

B. The motivation and application of MOCCA · 1) Summary of the Chambolle-Pock (CP) primal-dual framework: · 2) The need to generalize the CP primal-dual framework:

The paper presents the CP primal-dual framework for minimizing convex objectives with potentially nonsmooth terms, then explains why the quadratic surrogate Q requires a generalized framework. The key obstacle is that Q is convex in f but not in z, preventing the standard convex-conjugate CP formulation.

  • 1) Summary of the Chambolle-Pock (CP) primal-dual framework:: The CP framework minimizes F(Kx) + G(x), where F and G are convex, possibly nonsmooth functions and K is a matrix acting on x.Its ability to handle nonsmooth convex functions is important for enforcing the paper’s convex image constraints.
  • 1) Summary of the Chambolle-Pock (CP) primal-dual framework:: The CP algorithm solves a saddle-point problem while simultaneously obtaining the solution to the corresponding primal minimization problem.The formulation uses convex conjugation to embed the primal problem in a larger primal-dual saddle-point problem.
  • 1) Summary of the Chambolle-Pock (CP) primal-dual framework:: The primal-dual gap tends to zero as iterations increase, providing a convergence check because the primal and dual objectives are solved simultaneously.The CP step sizes must satisfy a norm-based condition, and diagonal step-size matrices can replace scalar step sizes when appropriate.
  • 2) The need to generalize the CP primal-dual framework:: Applying CP to Q for fixed f0 requires rewriting Q in the form F(Kf) + G(f), with D and E depending on c and f0, b depending on c and f0, and K depending only on f0.Both terms of Q depend on Kf, so Q is identified with F in the CP objective.
  • 2) The need to generalize the CP primal-dual framework:: Although Q is convex in f, its associated function FQ is not convex in z and instead is a difference of convex functions.The decomposition uses convex functions FQ+(z) and FQ−(z), with nonnegative matrices D and E enabling this difference-of-convex representation.
  • 2) The need to generalize the CP primal-dual framework:: Because FQ is not convex in z, it cannot be written as the convex conjugate of F*, so the standard CP update no longer yields the original objective.This failure motivates generalizing the CP primal-dual framework for the quadratic surrogate Q.

3) Heuristic derivation of MOCCA: · 4) Application of MOCCA to optimization of the spectral CT data fidelity:

MOCCA heuristically extends convex-concave primal-dual optimization by locally approximating the convex component around a dual-mirrored primal point. Applied to spectral CT fidelity, repeated local re-expansion enables a constrained one-step algorithm without nested loops, while introducing a convergence–stability control parameter.

  • 3) Heuristic derivation of MOCCA:: MOCCA applies a convex approximation around a dual-mirrored primal expansion point, then iterates a convex primal-dual update while relocating that point.The expansion point is selected to mirror the current dual variables through the Legendre-transform maximizer.
  • 3) Heuristic derivation of MOCCA:: The resulting MOCCA iteration combines the mirrored expansion-point update and convex-conjugate dual step with the remaining primal-dual steps inherited from CP.The simulations initialize all variables to zero, while convergence analysis is referenced in an accompanying paper.
  • 4) Application of MOCCA to optimization of the spectral CT data fidelity:: For spectral CT fidelity, the method repeatedly re-expands a local quadratic data-discrepancy approximation at the current material maps and applies MOCCA.Core MOCCA operations are called inner iterations, whereas re-expansion and MOCCA application form outer iterations; convex image constraints can be incorporated.
  • 4) Application of MOCCA to optimization of the spectral CT data fidelity:: A single inner iteration is sufficient, eliminating nested iteration loops and yielding an algorithm that re-expands the spectral CT fidelity at every iteration.The resulting one-step procedure initializes the material-map, expansion-point, and dual variables to zero vectors.
  • 4) Application of MOCCA to optimization of the spectral CT data fidelity:: Re-expansion reduces the system size because only the first matrix block contributes to the updates, while the second block becomes trivial.Vector step sizes are computed using matrix-vector products rather than a power method, making repeated re-expansion practical.
  • 4) Application of MOCCA to optimization of the spectral CT data fidelity:: Because re-expansion at every step is not guaranteed to converge, λ balances convergence rate against potentially unstable iterations.λ is introduced while maintaining the relevant preconditioned step-size product, and its effect on convergence is demonstrated in Section IV.
  • 4) Application of MOCCA to optimization of the spectral CT data fidelity:: The presented update equations analytically evaluate the dual and primal minimizations for spectral CT fidelity, with primal and prediction steps matching the corresponding CP steps.This formulation addresses data-fidelity optimization alone; the full version with TV constraints is provided separately in Appendix C.

C. One-step algorithm µ-preconditioning

Spectral CT attenuation curves can be nearly linearly dependent, making material separation sensitive to inconsistent data and degrading reconstruction conditioning. The proposed µ-preconditioning orthogonalizes these curves before applying the constrained one-step algorithm.

  • Motivation: Near-linear dependence of the attenuation matrix makes material-map separation sensitive to inconsistent spectral CT data and worsens reconstruction conditioning.The poor conditioning also negatively impacts algorithm efficiency.
  • Preconditioning: µ-preconditioning orthogonalizes the attenuation curves to address the conditioning-related efficiency problem.The step is distinguished from preconditioning used in the CP algorithm.
  • Preconditioning: The method constructs a preconditioning matrix from an eigenvalue decomposition of a matrix formed from the attenuation matrix, with ordered eigenvalues defining its singular values.The supplied passage identifies the singular values through the ordered eigenvalues and their ratio.
  • Implementation: Implementation transforms the material maps and attenuation matrix, substitutes them into the one-step algorithm, and reformulates constraints in terms of transformed maps.The vector stepsizes require care because µ-preconditioning can introduce negative entries into the transformed linear operator.

D. Convergence checks · IV. RESULTS · A. Head phantom studies with material map TV-constraints

The studies use primal-dual gap and TV metrics to assess convergence, showing that constrained one-step spectral CT recovers head-phantom material maps under ideal and Poisson-noisy data. Performance depends on λ and TV constraints, with TPL-TV and LSQ-TV exhibiting different RMSE behavior.

  • D. Convergence checks: Convergence is checked through the primal-dual gap, defined from the convex quadratic approximation and its Fenchel dual objective.When additional constraints are included, the primal-dual gap is modified and material-map estimates should be checked against those constraints.
  • IV. RESULTS: The experiments use simulated ideal photon-counting transmission data and include a two-material FORBILD head phantom plus a more realistic anthropomorphic chest-phantom study.The results section defines separate studies for convergence/material-map recovery and anthropomorphic application.
  • A. Head phantom studies with material map TV-constraints: For ideal noiseless data, the cPD gap and data discrepancy tend to zero while material-map TVs converge to their designed values.These convergence metrics demonstrate convergence for the particular problem studied, but the data discrepancy is not generally a convergence check with inconsistent data.
  • A. Head phantom studies with material map TV-constraints: Progress toward convergence depends on λ, so the algorithm requires a search over λ, and no theoretical convergence guarantee follows from the re-expansion step.The noisy simulations nevertheless show convergence for both TPL-TV and LSQ-TV, with convergence rates similar to the noiseless results.
  • A. Head phantom studies with material map TV-constraints: The material-map estimates converge to the phantom maps despite projection-view undersampling, indicating that TV constraints combat undersampling artifacts.Difference images also show that convergence can be spatially non-uniform.
  • A. Head phantom studies with material map TV-constraints: In noisy simulations, TPL-TV yields lower image RMSE than LSQ-TV, with a larger gap under looser TV constraints.The data are generated using a transmission Poisson model, although RMSE may not fully represent task-based image quality.

B. Chest phantom studies with a mono-energetic image TV constraint

Chest phantom simulations apply the one-step TPL-monoTV algorithm to reconstruct water, bone, and Gadolinium basis maps from noisy spectral CT data. The resulting 70 keV monoenergetic images show high-contrast structures with few artifacts across TV-constraint settings.

  • Chest phantom studies with a mono-energetic image TV constraint: The study used an anthropomorphic phantom with 24 material/density combinations and reconstructed water, bone, and Gadolinium basis maps.The phantom was created by segmenting an actual CT chest image, and the transmitted counts followed a Poisson model with 4×106 photons per detector pixel.
  • Chest phantom studies with a mono-energetic image TV constraint: TPL-monoTV was selected because it requires one monoenergetic-image TV constraint parameter rather than three parameters for TPL-TV.The selected optimization problem was intended to demonstrate a simpler constrained one-step reconstruction on the realistic phantom.
  • Chest phantom studies with a mono-energetic image TV constraint: Unregularized TPL reconstructions demonstrated the impact of simulated noise on the 70 keV monoenergetic image.The phantom and unregularized reconstruction were displayed at 70 keV for comparison.
  • Chest phantom studies with a mono-energetic image TV constraint: TPL-monoTV produced 70 keV monoenergetic chest images showing high-contrast structures with few artifacts under different TV constraints.The images include global views and a right-lung ROI highlighting soft-tissue detail.

V. CONCLUSION

The paper develops a constrained one-step algorithm that directly inverts spectral CT transmission data into basis material maps using local quadratic majorization and convex-concave optimization. Simulations demonstrate TV-constrained inversion, while future work targets robustness to data inconsistencies and spectral-response auto-calibration.

  • V. CONCLUSION: The constrained one-step algorithm directly inverts spectral CT transmission data into basis material maps while addressing non-convex discrepancy terms with a local convex quadratic upper bound.The method is derived for TPL and LSQ data discrepancies, with the strategy applicable to other data fidelities.
  • V. CONCLUSION: Simulations demonstrate TV-constrained data-discrepancy minimization with constraints applied either to individual basis maps or to an estimated monochromatic X-ray attenuation map.The algorithm derives from MOCCA to address an intermediate problem arising from the local convex quadratic approximation.
  • V. CONCLUSION: Basis-map constraints should help control artifacts caused by data inconsistencies from spectral miscalibration, X-ray scatter, and physical photon-counting detection processes.Robustness to these inconsistencies is identified as future work.
  • V. CONCLUSION: Future work also pursues generalization of the algorithm to enable auto-calibration of the CT system’s spectral response.

APPENDIX A GRADIENT OF LTPL · APPENDIX B

Appendix A derives the gradient of LTPL and motivates the linear transform A. Appendix B simplifies A’s indexing and proves inequalities establishing positive-semidefinite Hessians through vector, Cauchy–Schwarz, and nonnegative-coefficient arguments.

  • APPENDIX A GRADIENT OF LTPL: Appendix A derives the LTPL gradient in Eq. (12), motivating the definition of the linear transform A.The derivation recalls Eqs. (3) and (16) and explicitly states the LTPL gradient.
  • APPENDIX A GRADIENT OF LTPL: The LTPL gradient derivation continues by inserting I_ℓℓ′ into the algebraic manipulation.The remaining gradient and Hessian calculations are described as following by similar manipulations.
  • APPENDIX B: Appendix B collapses A’s double indices into single indices, with s = w · N_ℓ + ℓ and t = i · N_ℓ + ℓ′.This reindexing gives A_s,t(f) = A_wℓ,ℓ′i(f).
  • APPENDIX B: An inequality is used to prove that the Hessians in Eqs. (20)–(23) are positive semidefinite by setting b to the relevant residual vectors.The specified choices are r−(f), r+(f), and r(log).
  • APPENDIX B: The proof of Eq. (61) expands b in unit vectors and reduces the argument to showing the required inequality for any vector u.The derivation fixes s and defines v using components A(f)_s,t.
  • APPENDIX B: Because v_t ≥ 0, the left-hand side of Eq. (62) can be bounded using the Cauchy–Schwarz inequality.The appendix derives the left-hand side from v’s definition and then applies Cauchy–Schwarz.
  • APPENDIX B: The appendix proves Eq. (62), then uses it to establish Eq. (61).The latter step relies on b_s ≥ 0, making the expression a linear combination of positive definite matrices with non-negative coefficients.

APPENDIX C

Appendix C derives the constrained one-step TPL-TV and LSQ-TV algorithms by expressing the TV-constrained problem in primal-dual form and deriving its dual and proximity updates. It concludes with pseudocode using μ-preconditioning and inverse transformation of the final material maps.

  • Writing constrained TV optimization in the general form F(Kx)+G(x): The derivation formulates TV-constrained optimization as F(Kx)+G(x), enabling construction of the primal-dual algorithm through convex conjugates and dual functions.The data discrepancy is augmented with convex constraints on material-map TV values and written using indicator functions.
  • Convex conjugate of F2: The TV-gradient conjugate reduces to a magnitude maximization, with the maximizing dual direction aligned to the largest component of the input gradient vector.The largest component is identified with the infinity norm, yielding the corresponding maximization expression.
  • The pseudocode for TPL-TV and LSQ-TV: The complete TPL-TV and LSQ-TV pseudocode uses μ-preconditioning, and the final material maps are recovered by applying the inverse preconditioner after N iterations.For the reported results, all variables are initialized to zero.
Loading 1511.03384v1…