PaperlayerView exact revision on arXiv

Source-linked AI summary

Unified Projection-Domain Material Decomposition for Energy-Integrating and Photon-Counting CT

Dale Black, Surya Naidu, Hamidreza Khodajou-Chokami, Sabee Molloi

arXiv:2610.30440v1eess.IVphysics.med-ph

TL;DR

Projection-domain material decomposition must invert a nonlinear polychromatic model, while reducing photon-counting measurements to two can discard spectral information. This paper uses nested one-dimensional searches and a Poisson likelihood to decompose two materials from arbitrary numbers of measurements without rebinning. In simulation, retaining all four photon-counting measurements reduced noise at every VMI energy and improved quantitative accuracy at 70 keV and above compared with two merged measurements.

  • Problem

    Projection-domain material decomposition requires inversion of a nonlinear polychromatic model, and reducing photon-counting measurements to two discards spectral information and makes results depend on how windows are merged.

    Source: §1 · pp. 2–3

  • Method

    A nested optimization method estimates two-material line integrals from arbitrary numbers of spectral measurements using Poisson-likelihood profile minimization and two one-dimensional searches without rebinning.

    Source: §1 · p. 3

  • Results

    Using all four photon-counting measurements reduced noise at every VMI energy and improved quantitative accuracy at 70 keV and above compared with two merged measurements.

    Source: §5 · p. 24

  • Takeaways & Limitations

    One projection-domain method was evaluated across rapid kVp-switching, dual-source, and photon-counting acquisitions, providing a basis for decomposition across systems with more energy-resolved measurements.

    Source: §5 · p. 24

    The estimator was evaluated in simulation rather than on measured projections, and the simulated systems only approximate physical systems because vendor spectra, calibrations, and processing were unavailable.

    Source: §4 · p. 24

Abstract

from arXiv · show

Objective: To develop a unified and efficient method for projection-domain material decomposition in spectral computed tomography (CT) that incorporates an arbitrary number of spectral measurements, independent of the acquisition mechanism, without reducing the measurements to two effective energies. Approach: A nested optimization method was developed for two-material basis decomposition from K>=2 spectral measurements. It extends an existing two-measurement univariate method by replacing its polynomial approximation with profile minimization of the Poisson negative log-likelihood, reducing the inversion to two one-dimensional searches that use all measurements without spectral rebinning. The method was evaluated in simulation for energy-integrating (GE Revolution Apex Elite, Siemens SOMATOM Flash) and photon-counting (Siemens NAEOTOM Alpha) CT. Main results: The method recovered line integrals with noiseless iodine errors below 3e-6 g/cm^2 and noise within 3% of the Cramer-Rao lower bound (CRLB). Virtual monoenergetic images (VMIs) agreed, at 70 keV and above, with reference CT numbers and, on all three systems, with vendor VMIs of physical scans (mean absolute error 6-22 HU). For photon counting, four measurements instead of two reduced the iodine CRLB by 19% and image noise by 8-16% at every VMI energy, and lowered the regression error at 70 keV and above. Iteration and tolerance settings changed CT numbers by less than 0.04 HU, and decomposition took 0.34-0.67 s per million rays on an A100 GPU. Significance: The method provides fast, accurate decomposition for any spectral CT acquisition mechanism and number of measurements, so one method can serve every scanner in a virtual imaging trial.

1 Introduction

Spectral CT estimates material-specific attenuation from measurements at multiple energies, supporting applications such as material-specific and virtual monoenergetic imaging. The paper extends projection-domain two-material decomposition to arbitrary numbers of spectral measurements, avoiding information loss from reducing photon-counting measurements to two.

1 Introduction

  • Two-material basis decomposition estimates line integrals for two materials, such as water and iodine, which can be combined into material maps or virtual monoenergetic images.Additional basis functions may be needed when K-edge contrast materials are present.Source: §1 · p. 2
  • Projection-domain decomposition estimates basis line integrals before reconstruction, reducing beam-hardening artifacts while avoiding the spatial regularization and parameter tuning of one-step methods.It still requires accurate spectral and detector models or calibration measurements.Source: §1 · p. 2
  • Reducing photon-counting energy windows to two effective measurements discards spectral information and makes decomposition dependent on how windows are merged.Although narrower windows have lower individual signal-to-noise ratios, all measurements contribute jointly to decomposition.Source: §1 · p. 3
  • The method extends univariate two-material decomposition to any number of spectral measurements through two coordinated one-dimensional searches, using a Poisson-likelihood profile minimization without rebinning.It replaces the earlier polynomial relation and uses each measurement through its own effective spectral response.Source: §1 · p. 3
  • The method was evaluated on rapid kVp-switching, dual-source, and four-window photon-counting acquisitions using noiseless and noisy simulations compared with physical scans.Evaluation covered line integrals, synthesized VMIs, accuracy, noise, convergence, and computational efficiency against Cong et al..Source: §1 · p. 3

2 Methodology

The measurement model represents two-material attenuation with energy-dependent basis functions and line integrals. Acquisition differences enter through each measurement’s spectral response, allowing the same model to accommodate different detector and spectrum configurations.

2.1 Two-material basis measurement model

  • The model represents attenuation using two energy-dependent basis functions, with iodine and water line integrals as the two basis quantities in this work.The basis functions are the selected materials’ known mass-attenuation coefficients, and the corresponding density maps describe their spatial distributions.Source: §2.1 · p. 3
  • For K ≥ 2 measurements, each measurement has a spectral response that includes incident spectrum, filtration, exposure, acquisition geometry, and detector response.This response encodes the acquisition mechanism without changing the forward-model form.Source: §2.1 · p. 4
  • The discretized model uses energy samples and per-ray response contributions, so detector-specific spectral variation such as a bowtie filter requires no estimator change.Each q_ke is the expected unattenuated contribution of an energy sample to measurement k, including the energy-sampling interval.Source: §2.1 · p. 4
  • The expected measurement vector collects the expected counts for all K measurements at the same basis line integrals.Its components are λ_1(A,C) through λ_K(A,C).Source: §2.1 · p. 4

2.2 Statistical model and objective function

The statistical model treats measurements as conditionally independent Poisson variables given the two basis line integrals. The decomposition estimate is defined within a bounded feasible region so its one-dimensional searches remain finite.

2.2 Statistical model and objective function

  • After scatter correction, measured counts are modeled as conditionally independent Poisson variables with means given by the forward model.For ideal nonoverlapping photon-counting windows, this follows from Poisson splitting; residual detector correlations are treated as model mismatch.Source: §2.2 · p. 4
  • The estimator defines the joint material-decomposition estimate using the negative log-likelihood, omitting terms independent of the basis line integrals.The passage introduces this objective after defining observed measurements and their Poisson model.Source: §2.2 · p. 4
  • The optimization is restricted to a bounded feasible box whose limits keep both one-dimensional searches finite.The box is Θ = [A_min, A_max] × [C_min, C_max].Source: §2.2 · p. 4

2.3 Gradient and Fisher information

The method uses derivatives of the expected counts to form the objective’s gradient and Fisher information. The inverse Fisher information at the true basis line integrals gives the CRLB used to assess estimator efficiency.

2.3 Gradient and Fisher information

  • The update procedure uses the negative log-likelihood gradient and Fisher information as its expected curvature.Both are formed using derivatives of the expected counts.Source: §2.3 · p. 4
  • For the independent-Poisson model, the Fisher information matrix is computed from the expected-count derivatives.The passage presents the derivative-based gradient and Fisher information components.Source: §2.3 · pp. 4–5
  • At the true basis line integrals, the inverse Fisher information is the CRLB, a lower bound on covariance for unbiased estimators under usual regularity conditions.The reported CRLB values are standard deviations obtained from the square roots of its diagonal elements.Source: §2.3 · p. 5

2.4 Initialization

A stabilized log-linear estimate initializes the nonlinear decomposition, using all K measurements. Total-count matching can then refine the initial water-basis value before likelihood optimization.

2.4 Initialization

  • An approximate log-linear solution initializes the nonlinear estimator using stabilized log measurements and effective attenuation coefficients.A positive ε prevents taking the logarithm of zero.Source: §2.4 · p. 5
  • The initial basis estimate comes from a linear least-squares problem whose 2 × 2 normal equations use all K measurements.If the equations are singular, the estimate is set to (A, C) = (0, 20) g cm−2; it is projected onto Θ before nonlinear optimization.Source: §2.4 · p. 5
  • At fixed initial A, total-count matching refines C by bisection when a valid bracket exists; otherwise, the least-squares estimate is retained.This refinement initializes the likelihood optimization and does not constrain the final estimate.Source: §2.4 · p. 5

2.5 Nested optimization

Nested profile minimization rewrites the two-parameter estimator as searches over A and its conditional minimizer C⋆(A). Projected Fisher-scoring updates provide the inner and outer numerical implementation, with safeguards and diagnostics.

2.5 Nested optimization

  • The profile-minimization formulation reduces the joint inversion to one-dimensional searches while retaining equivalence to the original optimization when each inner problem is solved globally.It does not guarantee convexity of either objective.Source: §2.5 · p. 6

2.5.1 Inner optimization

  • At fixed A, the inner problem uses projected Fisher-scoring updates for C, stopping at the tolerance or after Ninner updates.The update is projected onto the feasible interval and limited in magnitude.Source: §2.5.1 · p. 6

2.5.2 Outer optimization

  • The outer optimization approximates profile-objective curvature with the Fisher-information Schur complement and reprofiles C after the final A update.Finite iteration limits, parameter bounds, step limits, curvature floors, and solver diagnostics provide numerical safeguards.Source: §2.5.2 · pp. 6–7

2.6 Bounded numerical reference

A bounded reference solver independently evaluates the nested profile minimization for comparison with the accelerated implementation.

2.6 Bounded numerical reference

  • The reference solver samples 65 C values over [Cmin, Cmax] at each tested A, refines detected basins with 80 golden-section iterations, and evaluates both endpoints.It applies the same sampling-and-refinement approach to the outer profile over [Amin, Amax], independently of the accelerated solver's initialization and updates.Source: §2.6 · p. 7

2.7 Denoising and VMI generation

Photon-counting basis sinograms receive likelihood-based bilateral denoising before reconstruction, while energy-integrating sinograms do not. Reconstructed basis maps are denoised jointly before VMI generation and CT-number calculation.

2.7 Denoising and VMI generation

  • Photon-counting basis sinograms are denoised with a total-likelihood bilateral filter, which is not applied to energy-integrating acquisitions.The filter weights neighbors by the Poisson likelihood of measured photon counts, which energy-integrating detectors do not record.Source: §2.7 · p. 7
  • The two basis sinograms are reconstructed independently with FBP, then anticorrelated noise reduction is applied jointly to the basis images.VMIs are generated from the denoised basis images at energy E0.Source: §2.7 · p. 7
  • CT numbers are computed relative to the water attenuation coefficient at E0, using the VMI and denoised iodine and water density maps.The reference CT values use the same definition.Source: §2.7 · p. 7

2.8 Experiments

The experiments compare decomposition across energy-integrating and photon-counting CT configurations using simulated acquisitions and physical scans of the Gammex phantom. Evaluation covers decomposition, VMI accuracy and noise, statistical efficiency, numerical sensitivity, and runtime.

2.8.1 Simulated acquisitions and phantoms

  • Three simulated configurations represented two energy-integrating systems and one photon-counting system, with the Gammex phantom also scanned physically on all three systems.A second tissue-characterization phantom was evaluated only in simulation.Source: §2.8.1 · pp. 8–9, Figure 1 · p. 8

2.8.2 Benchmarking

  • The proposed method was compared with Cong et al.’s two-measurement method on dual-energy data and on photon-counting data using either four windows or two merged measurements.The comparisons separate the effect of measurement count from that of the algorithm.Source: §2.8.2 · pp. 9–10

2.8.3 Photon-counting corrections

  • Photon-counting corrections addressed pile-up and scatter in the counts, while charge sharing and fluorescence escape were modeled in the spectral response.An additional no-correction acquisition left pile-up and scatter in the counts; decomposition used either an ideal or pile-up-adjusted response and four or two merged windows.Source: §2.8.3 · p. 10

2.8.4 Evaluation metrics

  • Evaluation assessed line-integral accuracy against known values and a bounded reference, and compared variability over noise realizations with the CRLB.Noiseless tests covered 145 known water-and-iodine path-length pairs; statistical efficiency used 20,000 Poisson realizations for a representative ray.Source: §2.8.4 · pp. 10, 12
  • VMI accuracy was measured against calculated reference CT numbers and physical-scan VMIs, using MAE, maximum error, regression slope, and residual error.VMIs were treated as a readout of the decomposition; physical-scan noise was shown only for reference because vendor processing did not allow controlled noise comparison.Source: §2.8 · p. 8, §2.8.4 · p. 10
  • Image noise, numerical sensitivity, and decomposition time were also assessed, with runtime measured per million rays on an A100 GPU.Noise was measured as HU standard deviation in solid water and mean within-insert standard deviation; sensitivity controls were varied from their default settings.Source: §2.8.4 · pp. 10, 12

3 Results

The nested method performed accurately and efficiently across the evaluated configurations; retaining four photon-counting measurements improved several noise and accuracy metrics, while low-energy calcium errors and uncorrected scatter remained important qualifications.

3.1 Decomposition performance

  • Noiseless iodine line-integral error stayed within 2.7 × 10−6 g cm−2, empirical variability was within 6% of the CRLB, and runtime was 0.34–0.67 s per million rays on the GPU.The statistical result covered six acquisition–method combinations and 20,000 Poisson realizations.Source: §3.1 · p. 13
  • For photon counting, using K = 4 instead of K = 2 reduced the iodine CRLB standard deviation by 19%, and image noise was 8–12% lower at higher VMI energies.At 40 keV, K = 2 increased mean within-insert SD by approximately 19% relative to K = 4.Source: §3.1 · p. 13, §3.2 · p. 15

3.2 VMI evaluation

  • At 70 keV and above, dual-energy regression slopes ranged from 0.96 to 0.99; four photon-counting measurements improved CT-number accuracy at 70, 100, and 140 keV, but calcium error rose at 40 keV.At 40 keV, Gammex 472 MAE increased from 35.7 to 38.3 HU; iodine MAE fell from 62.9 to 33.5 HU while calcium MAE rose from 17.7 to 70.7 HU. The calcium mismatch reflects the water/iodine basis representation, and every configuration under-read high-density calcium inserts at 40 keV.Source: §3.2 · p. 15

3.3 Comparison with physical scan measurements

  • For photon counting, K = 4 agreed more closely with physical scans than K = 2 across evaluated VMI energies, with 40-keV MAE falling from 70.1 to 41.7 HU.For the dual-energy systems, simulated-to-physical MAE was approximately 6–11 HU at 70 keV and above for rapid kVp switching, and approximately 11 HU at 70 keV rising to 22 HU at 140 keV for dual-source acquisition.Source: §3.3 · p. 18

3.4 Numerical sensitivity

  • Numerical sensitivity was limited over the evaluated ranges: remaining estimator parameters changed reconstructed CT numbers by less than 0.04 HU, while substantially restrictive bounds or update sizes caused appreciable degradation.Reducing the upper water bound below the approximately 33 g cm−2 maximum path length changed insert HU by approximately 1450 HU.Source: §3.4 · p. 19

3.5 Effect of photon-counting corrections

  • With uncorrected counts, an ideal response produced severe fitting failure, whereas including pile-up in the response recovered much of the performance; residual scatter left the largest error at 40 keV.With pile-up-adjusted response, K = 4 had lower MAE than K = 2 at every energy; at 40 keV MAE was 242.4 versus 281.3 HU, while all-corrections MAE was 38.3 HU.Source: §3.5 · p. 20, Table 6 · p. 22

4 Discussion

The method generalizes projection-domain decomposition to different spectral CT systems and measurement counts, with estimates near the CRLB and vendor VMI agreement at 70 keV and above. Four photon-counting measurements improve noise and quantitative accuracy over two, while calcium basis mismatch and simulation-only evaluation bound the conclusions.

Discussion

  • Vendor VMIs of physical scans agreed with the proposed method’s VMIs to within 6–22 HU MAE at 70 keV and above across three simulated acquisitions.For a representative ray, estimates approached the CRLB with small iodine bias (at most 0.04 standard deviations).Source: §4 · p. 23
  • For K = 2, the proposed method closely matched the univariate method in VMIs, CT-number accuracy, regression performance, and image noise, while reducing noiseless iodine error 10–15-fold.It required 1.8–2.2 times the decomposition time but remained below one second per sinogram on the GPU.Source: §4 · p. 23
  • Four photon-counting measurements lowered image noise at every VMI energy and improved quantitative performance at 70 keV and above compared with two merged measurements.The four-window configuration also reduced iodine CRLB standard deviation by 23% relative to combining the four windows into two.Source: §4 · p. 23
  • At 40 keV, calcium error increased with K = 4 and was attributed to the water/iodine basis rather than the solver, which recovered noiseless iodine line integrals within 2.7 × 10−6 g cm−2.Physical 40-keV images read the densest calcium inserts below calculated values, so the deviation was not specific to the proposed estimator.Source: §4 · pp. 23–24
  • Across evaluated settings, iteration limits, Fisher condition limits, and convergence tolerance changed CT numbers by less than 0.04 HU unless bounds or step limits were too restrictive.Stable performance required a water bound above the longest object path and iodine bounds that allowed per-ray noise.Source: §4 · p. 24
  • The method tolerates modeled photon-counting detector nonidealities, but scatter still required correction in this simulation.With pile-up, charge sharing, and escape included in the response, results were within a few HU of the fully corrected result at 100 and 140 keV.Source: §4 · p. 24
  • The evaluation used simulated data, not measured projections, and each system was modeled from reported parameters without vendor spectra, calibrations, or processing.Applying the estimator to measured projections with calibrated spectral models remains future work; the method was not benchmarked against a joint two-dimensional likelihood search.Source: §4 · p. 24

5 Conclusion

The paper presents a unified projection-domain decomposition method for spectral CT with arbitrary measurement counts. Four photon-counting measurements improved noise and quantitative accuracy over two merged measurements, supporting one consistent approach for simulations of multiple spectral CT systems.

5 Conclusion

  • Using all four photon-counting measurements reduced noise at every VMI energy and improved quantitative accuracy at 70 keV and above compared with two merged measurements.The method was evaluated across rapid kVp-switching, dual-source, and photon-counting acquisitions.Source: §5 · p. 24

Data and code availability

The proposed estimator and the univariate method are included in the open-source BasisSimulator.jl package from version 0.15.0; simulations in the paper used version 0.14.0.

Data and code availability

  • The estimator and univariate method are available in BasisSimulator.jl from version 0.15.0, while the paper’s simulations used version 0.14.0.The simulations can be regenerated with the package, and the Gammex 472 phantom is included in its documentation.Source: Data and code availability · p. 25
Loading 2610.30440v1…