Dual-kVp Rapid-Switching Virtual Monoenergetic Imaging

A GE Revolution Apex Elite rapid-kVp-switching acquisition (80 + 140 kVp) of the Gammex 472 phantom, turned into virtual monoenergetic images (VMIs) at 50, 70, 100 and 140 keV by one package call, BS.vmi_pipeline, and checked rod by rod against a noise-free acquisition of the same scan and against the attenuation of the phantom's own materials.

The scanner, the acquisition and the chain are those of the basis-spectral-denoising paper: its GE model with its physics switches, the two exposures of a rotation with their duty cycle, view integration, the fitted pair of FDK windows and a calibration draw that fixes the spectral pair.

Pipeline

80 kVp + 140 kVp exposures of one rotation   (duty cycle 0.65 / 0.35, each view integrated over its arc)
  → one projection per exposure, three noise draws from it:
      the noise-free reference, the measured scan and a calibration scan
  → BS.spectral_basis_from_acquisitions       (per-ray absolute response of each kVp)
  → calibration draw → pair_basis (E*, β) and composite_energy
  → BS.vmi_pipeline(; denoiser = SpectralHYPR(), fbp_filter = PairFilter(…), pair_basis, composite_energy)
       projection-domain HYPR-LR on the counts of each detector row
     → K = 2 n-channel Poisson maximum-likelihood decomposition (iodine + water, every ray)
     → FDK of the composite and of its complement, each with its own window
     → ACNR on the complement → image-domain HYPR on the complement
     → VMI synthesis at 50 / 70 / 100 / 140 keV
  → scored against the noise-free reference and against theory

Both kVp channels stay separate all the way into the likelihood; nothing is rebinned into a single "effective" spectrum.

Notebook Setup

Backend detected: CuArray

Scan Setup and Simulation

1. Phantom() Struct

The Gammex 472 multi-energy phantom: a solid-water body with seven calcium (50–600 mg/mL) and seven iodine (2–20 mg/mL) inserts, 1 cm thick and z-invariant.

2. EICTScanner() Struct

The Apex Elite geometry (source–isocenter 625.6 mm, source–detector 1100 mm, 834 columns × 0.6 mm on an arc with the quarter-detector offset, 256 rows × 0.625 mm), its Gemstone scintillator (:lumex) and the GE large-body bowtie: the basis-spectral-denoising model, value for value.

Rapid kVp switching alternates the tube between 80 and 140 kVp within every view period, on one tube and one detector, so the two channels see the same rays. It is modelled as two co-registered acquisitions on this one geometry, one per kVp, each with its share of the view period: its duty cycle, 0.65 at 80 kVp and 0.35 at 140 kVp. A channel's tube current is the instantaneous current times its duty cycle, and its views integrate over that fraction of the arc the gantry turns per view (view_arc).

3. CTProtocol() Struct

One protocol per exposure: 984 views over a 0.5 s rotation, 5 mm collimation, 4.5 mm of added aluminium on both. The physical scan's currents are 407 mA at 80 kVp and 405 mA at 140 kVp, each multiplied by its duty cycle.

The simulator's tube output per mA is generic, not GE's, so the tube currents alone do not fix the dose. Both are scaled by one factor, DOSE_SCALE, which keeps the split between the two exposures and makes the simulated CTDIvol (the package's Monte Carlo dose in the 32 cm body phantom, summed over the exposures) equal to the 10.07 mGy the physical scan reported. Noise is then simulated at the dose the physical scan delivered.

DOSE_SCALE = 1.0314. The simulated exposures:

4. SimOptions() & ReconOptions()

The physics switches of the basis-spectral-denoising GE model: every effect of the energy-integrating path on except optical crosstalk, and the detector efficiency from the Monte Carlo table of the Lumex scintillator (detector_efficiency_mode = :mc_lut).

View integration. A detector integrates while the gantry turns, so each view records the transmitted intensity averaged over the arc it sweeps, r · view_arc · Δθ at radius r: a blur of the object, largest far from the isocentre, that the noise, counted once per view, does not share. view_samples = 5 samples that arc with five sub-views (the midpoint rule), 0.2 mm apart at the edge of a 35 cm field. Without it the simulated signal is sharper, relative to its noise, than a physical scanner's, and no reconstruction window can match both the physical resolution and the physical noise texture.

Draws. Each exposure is simulated three times, differing only in the noise: the noise-free expectation (draw 0, use_noise = false), the measured scan (draw 1) and a calibration scan (draw 9) of its own seeds. The reconstruction grid is the clinical series', 8 slices of 0.625 mm over the central 5 mm, 512² over 35 cm.

5. Forward Project: simulate!

Each exposure is projected once. The noise-free draw runs first with keep_projection = true, which returns everything before the noise; the measured and calibration draws pass it back (projection) and only draw their noise. That is bit-identical to simulating each draw afresh, at a fraction of the cost.

Each acquisition keeps what the spectral basis needs: its corrected log sinograms, its per-ray air counts I0_ray (the detector's air count × the bowtie air profile) and its per-ray detected spectrum from resolve_source_spectrum_full (source × filtration × bowtie × detector response), the same model the forward projector applied.

The three draws of the two exposures took 51.8 s for 80 kVp and 37.8 s for 140 kVp, one projection each. Each sinogram is 834 × 12 × 984 (columns × rows × views).

VMI Pipeline

1. Spectral Basis from the Two Acquisitions

spectral_basis_from_acquisitions merges the two kVp energy grids onto their union and scales each acquisition's per-ray spectrum by its own air counts, so the likelihood sees the absolute response $\Phi_k(E)$ of every ray and channel: source, filtration, bowtie and detector, the model that generated the data. It does not depend on the noise, so every draw shares it.

The basis holds 2 channels on a 274-point energy grid, 4 to 140 keV, for 834 × 12 rays; the response sums to the air counts to within 1.1e-8 (relative).

2. The Chain's Settings

The denoiser is generalized HYPR-LR in both domains, BS.SpectralHYPR(): before the decomposition, a 3 × 3 (column × view) window on the counts of each detector row, pooling views of one parity so the odd- and even-view halves stay independent, with each channel's dispersion measured on the air rays; after the reconstruction, the complement of the basis pair pooled within each slice, over the window with the least estimated risk. Neither ever pools across detector rows or slices. The decomposition is the K-channel estimator with its published controls; T-LBF (a photon-counting filter) is off; the view-direction antialias is on; ACNR is on (vmi_pipeline's default).

A pair of FDK windows. A VMI is $A(E)\,M + B(E)\,I_\perp$, with $M$ the minimum-noise VMI (the composite, at $E^\star$) and $I_\perp$ its complement, the iodine component whose noise is uncorrelated with $M$'s. BS.PairFilter gives each its own window: the composite's sets the resolution and noise texture at $E^\star$, the complement's how they change away from it, and one pair serves every energy. A window per basis image instead would stop the strongly anti-correlated noise of the water and iodine images cancelling wherever the two windows differ, leaving uncancelled noise in every VMI. Each window is an apodized ramp $W(f) = \exp(-(f/f_c)^p)$ on the grid-Nyquist axis, the GE pair fitted by basis-spectral-denoising to the physical Gammex scan's MTF at every energy and NPS shape, on this scanner's clinical grid, which is this notebook's grid.

3. The Calibration Draw: pair_basis and composite_energy

The pair the windows act on, $(E^\star, \beta)$, and the energy of the composite that ACNR and the image-domain HYPR act on are properties of the scanner and protocol. They are measured once, on the calibration draw, and fixed for every other draw, exactly as basis-spectral-denoising's PAIR_BASIS does:

1. the calibration draw decomposed as measured and after the projection-domain HYPR alone (vmi_pipeline(; keep_sinograms = true), with SpectralHYPR(image = nothing) for the second);

2. BS.spectral_pair of the first with a standard soft-tissue window: its minimum-noise energy $E^\star$ and complement $\beta$, the pair_basis;

3. BS.spectral_pair of the second on that basis: the minimum-noise energy of the pair after the projection-domain HYPR, which moves it, the composite_energy.

Near its minimum the VMI noise hardly changes with energy, so the argmin of one acquisition is itself noise and would move from draw to draw. Measured on a draw of its own, it cannot tune the chain to the draw being scored, and the noise-free reference, which has no noise to measure it from, is reconstructed with the same pair as the measured scan.

4. vmi_pipeline

One call, from the two measured sinograms to the VMI stack, reconstructed on the notebook's grid (recon_opts.matrix_size) with every detector row kept: the published chain, with the GE's windows and the calibration draw's pair.

The noise-free reference is the noise-free draw through the same reconstruction, the same windows and the same pair, with nothing to denoise: no HYPR and no ACNR. The measured VMIs are scored against it, so the scores measure what the noise and the denoising do, separately from what the physics the decomposition does not invert exactly (scatter, and the blurs of the focal spot, the detector and the view integration) does to both.

The calibration draw fixed the windows' pair at E* = 78 keV, β = 0.0174, and the composite at 76 keV. On the measured draw the decomposition solved 9847872 rays in 1.8 s, 2.0 outer iterations on average; 0.0% did not converge and 0.006% / 0.0% touched the iodine / water bounds. The projection HYPR measured dispersions, the variance-to-mean ratio of the counts on the air rays, of 1.32 and 1.3 for 80 and 140 kVp. The image HYPR pooled the complement over a 31 × 31 window.

5. Basis Maps and VMIs

The water and iodine basis pair after ACNR and image HYPR (mid slice), and the VMIs synthesized from it. Every VMI comes from the same basis pair.

Results

Per-rod measured, noise-free and theoretical HU at the four VMI energies, the water-region mean and noise, and the linear regression of measured on theoretical HU.

• Measured HU: mean over an 8-px-radius circle at each rod's centroid, over every reconstructed slice, of the measured VMIs.

• Noise-free HU: the same mean on the noise-free reference.

• Theoretical HU: $1000\,(\mu_r(E) - \mu_w(E))/\mu_w(E)$ of the rod's material from BS.compute_μ_at_energy, with no fitting.

• Water region: the solid-water body eroded by 12 px, away from every insert and edge. Its noise σ is the standard deviation of the measured VMI minus the noise-free reference, the noise alone, without the phantom's own nonuniformity.

Water ROI

The eroded solid-water region, overlaid on the 70 keV VMI, and its mean HU at each energy, measured and noise-free.

Water-Region Noise

The standard deviation of the measured VMI about the noise-free reference, over the same eroded water region, at each energy.

Per-Rod Measured vs Theoretical

Solid lines show measured HU, crosses the noise-free reference and dashed lines theoretical HU, for the calcium and iodine inserts across the four energies.

Accuracy against the noise-free reference and against theory. Water: mean HU of the measured / noise-free VMI, and the noise σ. Rods: the mean absolute error against theory of the measured / noise-free VMI, and the mean and the largest absolute bias of the 14 measured rods against the noise-free reference.

The solid water reads -0.2 / -1.1 / -1.5 / -1.7 HU at 50 / 70 / 100 / 140 keV, within 0.7 HU of the noise-free reference: an offset the reference shares comes from physics the decomposition does not invert exactly (scatter, and the blurs of the focal spot, the detector and the view integration, taken before the log), not from the noise or the denoising. The noise falls from σ = 38.0 HU at 50 keV to 30.5 HU at 140 keV. The largest mean rod error against theory, over both rod groups and every energy, is 46.0 HU measured and 55.9 HU noise-free.

That largest error is the calcium rods at 50 keV, 46.0 HU measured and 55.9 HU noise-free — the noise-free reference has it too, so it is the basis, not the noise or the denoiser: the decomposition writes every material as water plus iodine, which spans calcium's energy dependence only approximately, and least well at low energy, where the photoelectric term dominates. From 70 keV up the calcium error is 2.3 to 15.0 HU.

Linear Regression

Measured rod HU regressed on theoretical HU at each energy. The dashed identity line is perfect agreement.

✅ Verification: PASS

Summary

Simulate 80 + 140 kVp (duty cycle 0.65 / 0.35, view_samples = 5, dose-matched to 10.07 mGy)
   one projection per exposure → noise-free, measured and calibration draws
   → BS.spectral_basis_from_acquisitions (absolute per-ray response of each kVp)
   → calibration draw → pair_basis, composite_energy
   → BS.vmi_pipeline(; denoiser = SpectralHYPR(), fbp_filter = PAIR_FILTER, pair_basis, composite_energy)
        projection HYPR → K = 2 n-channel decomposition on every ray
        → FDK of the composite and complement, each with its window
        → ACNR on the complement → image HYPR → VMI synthesis at 50 / 70 / 100 / 140 keV
   → water, noise, per-rod and regression checks against the noise-free reference and theory

Every VMI is synthesized from the same water/iodine pair; no energy-dependent filtering is applied to any VMI.