Photon-Counting CT Virtual Monoenergetic Imaging

A photon-counting CT scan of a Gammex 472 multi-energy phantom, taken all the way to virtual monoenergetic images (VMIs). The scanner is the Siemens NAEOTOM Alpha model of the basis-spectral-denoising study: a CdTe detector that sorts every photon into four energy windows (thresholds 20 / 35 / 55 / 70 keV), a 140 kVp, 0.5 s axial acquisition at the dose of the physical Gammex scan, and every view integrated over the arc the gantry turns while it is read.

The notebook shows three things:

1. The acquisition. Each simulate! call returns the four corrected energy-bin sinograms and the air count of every ray in every bin, I0[col, row, bin]. The bowtie makes that air response vary across the fan, so the spectral model is resolved per ray. The phantom is projected once; the noise-free reference, a calibration draw and the measured draw all reuse that projection.

2. The VMI chain. spectral_basis(ws; I0) builds the model the simulation applied, and one vmi_pipeline call turns the four bins into VMIs at 40, 70, 100 and 140 keV, with the denoiser, the reconstruction windows and the spectral pair of basis-spectral-denoising.

3. The check. Every calcium and iodine rod is compared with its theoretical HU from first-principles attenuation, and with the same chain applied to the noise-free draw. The solid-water background gives the HU accuracy and the noise at each energy.

Pipeline

simulate!(…; keep_projection = true)   noise-free draw; keeps the projection
simulate!(…; projection)               calibration and measured draws: only the noise is new
  →  4 corrected bins h_k = −log(y_k / I0_k)  +  I0[col, row, bin]
  →  spectral_basis(ws; I0)            per-ray absolute response Φ[col, row, E, k]

calibration draw  →  spectral_pair     →  pair_basis (E*, β), composite_energy
measured draw     →  vmi_pipeline(; denoiser = SpectralHYPR(), fbp_filter = PairFilter(…),
                                     pair_basis, composite_energy)
                       projection HYPR-LR on the counts of each detector row
                       K-channel maximum-likelihood decomposition → iodine + water sinograms
                       FDK of the composite and its complement, each with its own window
                       ACNR on the complement
                       image HYPR-LR on the complement
                       VMI synthesis at 40 / 70 / 100 / 140 keV
noise-free draw   →  the same reconstruction without denoising → the reference

All four bins enter the decomposition as separate measurements; nothing is merged into "low" and "high". Every setting is the one the basis-spectral-denoising study settled on for the NAEOTOM Alpha.

Notebook Setup

Backend detected: CuArray

The Acquisition

Phantom

The Gammex 472: a 33 cm solid-water body with seven calcium rods (50 to 600 mg/mL) on an inner ring and seven iodine rods (2 to 20 mg/mL) on an outer ring, 28 mm each. The label mask (512 × 512 at 35 cm, 16 slices over 1 cm) goes to the GPU; its material table stays on the host.

Scanner

The NAEOTOM Alpha geometry: 0.275 × 0.322 mm native CdTe dexels read out 2 × 2 binned, an arc detector with the quarter-detector offset, a 50 cm scan field behind the large-body bowtie, and a 1.6 mm CdTe sensor with 10 keV energy resolution, charge sharing and 5 ns dead time. Pile-up is simulated and corrected, and so is scatter; the counts are exact Poisson. These are the values of the basis-spectral-denoising and semmd-bayesian model. See the scanners page for where each value comes from.

Protocol and Dose

140 kVp, 0.5 s rotation, 1200 views, 5 mm collimation, 0.9 mm titanium added: the physical Gammex scan, which ran at 174 mA and reported a CTDIvol of 10.12 mGy. A tube's output per mA is its own and the simulator's source is generic, so the tube current is scaled by one factor, DOSE_SCALE, chosen so that the simulated CTDIvol (BasisSimulator's Monte Carlo dose in the 32 cm body phantom) equals the physical scan's. The noise is then that of the dose the physical scan delivered.

Simulation and Reconstruction Options

The physics switches of the basis-spectral-denoising Alpha model. Quantum noise, Compton scatter and the focal-spot blur are on. The switches turned off model a scintillator (fill factor, optical crosstalk, the scintillator efficiency table, lag): the CdTe detector's efficiency, fill and charge sharing are already in its Monte Carlo response. The heel effect is off.

View integration. A detector integrates while the gantry turns, so each view reads the transmitted intensity averaged over the arc it sweeps: 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 per view (the midpoint rule; at the edge of a 35 cm field they are 0.2 mm apart, a third of a pixel), and view_arc = 1.0 because a photon-counting detector reads for the whole view period. Without it the simulated signal is sharper, relative to its noise, than a physical scanner's.

Three draws, one projection. Everything before the noise (the projection through the phantom, sub-views included) is the same for every draw, so it is computed once, on the noise-free draw (keep_projection = true), and every other draw reuses it (projection) with only its own noise drawn. This is bit-identical to simulating each draw afresh.

The reconstruction grid is the physical scan's clinical series: 512 × 512 over 35 cm, twelve 0.4 mm slices.

Simulate: simulate! and spectral_basis

simulate! runs the whole photon-counting detector model on the GPU: polychromatic forward projection of every sub-view through the Monte Carlo detector response, scatter injection, Poisson counts in every bin, pile-up migration, then pile-up and scatter correction. It returns

• pcct_sino.bins: four corrected sinograms $h_k = -\log(y_k / I_{0,k})$, [col, row, view];

• I0_bins: the air count of every ray in every bin, [col, row, bin];

• projection, with keep_projection = true: everything before the noise, for the next draw;

• dose: the DoseReport of the beam.

spectral_basis(ws; I0) is read from the same workspace: the response the simulation applied, passed through each ray's own bowtie transmission, so the decomposition inverts exactly what was simulated and needs no calibration scan. It checks that the response sums to I0 in every ray and bin and refuses a basis that does not. The model does not depend on the noise, so one basis serves every draw.

The Four Energy Windows

Left: what each bin detects on the central ray, the absolute response the decomposition uses, as a function of the energy the photon arrived with. The windows overlap because the 10 keV energy resolution, charge sharing and K-fluorescence escape let a photon register below its true energy; the spikes are tungsten's characteristic lines. Right: the air count across the fan. The bowtie cuts the flux toward the edge and hardens what is left, so the high bins fall off less than the low ones. One scalar I0 per bin could not describe this, which is why the basis is resolved per ray.

The VMI Chain

The Denoiser

Generalized HYPR-LR in both domains, BS.SpectralHYPR():

• projection domain (ProjectionHYPR): within each detector row, each ray's split of its total count across the four bins is pooled over a 3 × 3 (column × view) window whose view neighbours are two views apart (view_stride = 2), so the odd and even views stay independent. Neighbours are weighted by how likely their total count is under the ray's own, and each ray keeps its own total. The count dispersion of every bin is measured from the air rays.

• image domain (ImageHYPR): the reconstructed basis pair is rewritten as the minimum-noise VMI (the composite) and a complement whose noise is uncorrelated with it. The composite is kept as reconstructed, so its resolution and noise texture are FDK's; the complement, which carries only spectral information, is pooled within its slice with likelihood weights on the composite, over the window width (3 to 41 pixels) with the least estimated risk. Every noise level comes from the difference of the odd- and even-view reconstructions.

Nothing in the chain averages detector rows or slices. These are SpectralHYPR()'s defaults, the chain of basis-spectral-denoising.

The Reconstruction Windows

The water and iodine images are the wrong pair to filter differently: their noise is strongly anti-correlated, and a VMI near the minimum-noise energy is quiet only where both were filtered alike. BS.PairFilter therefore filters the uncorrelated pair: the composite, whose window sets the resolution and noise texture at the minimum-noise energy, and the complement, whose window sets how that changes away from it. Each window is an apodized ramp, $W(f) = \exp(-(f/f_c)^p)$ on the grid-Nyquist axis, sampled at 11 knots as a CustomFilter.

basis-spectral-denoising fitted the Alpha's pair once, to the physical Gammex scan's MTF at every energy and its NPS shape (VMI, Br36f, QIR 0, 0.4 mm, 350 mm): composite $f_c = 0.5$, $p = 2$; complement $f_c = 0.45$, $p = 6$. A window is defined on the frequency axis of the grid it is applied on, so on a grid other than the fitted one (512² over 35 cm) its argument is rescaled by the ratio r of the grids' bandlimits (BS.grid_bandlimit). This notebook uses the fitted grid, so r = 1.

The Spectral Pair, Fixed from the Calibration Draw

spectral_pair measures, from the odd- and even-view halves of a reconstruction, the minimum-noise energy $E^*$ and the complement's coefficient $\beta$. 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; it is a property of the scanner and protocol, so it is measured once, on the calibration draw, and fixed for every other draw, the noise-free one included:

• pair_basis = (Estar, β): from the calibration draw decomposed without denoising and reconstructed with a standard soft-tissue window (SoftFilter); the pair the windows act on;

• composite_energy: the minimum-noise energy of the same draw after the projection-domain HYPR, which moves it; the composite ACNR and the image HYPR act on.

vmi_pipeline

One call runs the whole chain on the measured draw. Every detector row is kept and reconstructed onto the twelve-slice grid; the decomposition is the K-channel maximum-likelihood estimator (method = :nchannel) with its published controls; ACNR (on by default) acts on the complement only, before the image HYPR. The noise-free draw goes through the same call without the denoiser and ACNR: the reference, reconstructed exactly as its noisy counterpart.

What the chain fixed and measured:

Basis Images and VMIs

The chain's basis pair: iodine as a mass density (mg/mL) and water as a density (g/mL), the central slice. Every VMI is synthesized from this one pair with the monoenergetic attenuation of water and iodine; no energy-dependent filtering is applied.

Results

The phantom's labels are resampled onto the reconstruction grid (resample_to_recon). Each rod is measured in a disk of 60 % of its radius, and the background in the solid water eroded 12 pixels (8 mm) away from every edge. Means pool all twelve slices; the noise is the standard deviation within one slice, averaged (RMS) over the slices. The theoretical HU of every material is $1000\,(\mu_E - \mu_{E,\mathrm{water}})/\mu_{E,\mathrm{water}}$ from its composition. Every quantity is also measured on the noise-free reference, the same decomposition and reconstruction without noise (and so without the denoiser): its error against theory is the bias of the basis model itself, and the measured draw's departure from it is what the noise and the denoising add. The noise of the measured draw is the standard deviation of its difference from the reference, in which the phantom's structure and any deterministic artefact cancel.

Water ROI

The eroded solid-water region (red) on the central slice of the 70 keV VMI, and its mean HU at each energy, measured (bars) and on the noise-free reference (black diamonds). The phantom's solid water is modelled with the composition of water, so its theoretical HU is 0 at every energy (dashed line).

Noise versus Energy

The noise in the same solid-water ROI, per 0.4 mm slice: the standard deviation of the difference between the measured VMI and the noise-free reference, in which the phantom's structure and any deterministic artefact cancel. The standard deviation of the measured VMI itself differs from it by at most 0.06 HU at any energy, so the ROI holds noise only. The noise is highest at 40 keV, where the synthesis weights the iodine image most heavily and so amplifies its noise, and changes little from 70 keV up, where the VMI is mostly the composite and the water image.

Rods against Theory

Solid lines: measured HU of each calcium and iodine rod. Dashed lines: its theoretical HU.

The decomposition represents every voxel as iodine plus water. The solid water and the iodine rods are close to that basis; calcium is not exactly a combination of iodine and water at every energy, so its VMIs carry a small, energy-dependent model error on top of the noise. The per-rod table below separates the two with the noise-free reference.

Linear Regression

Measured rod HU against theoretical HU at each energy, with the least-squares line. The dashed line is identity.

Per-Rod Error

Measured minus theoretical HU for every rod at every energy, and in brackets the same for the noise-free reference: the bias of the iodine–water basis model itself, which no amount of dose would remove. The difference between the two is what the noise and the denoiser add. The largest relative error of the measured draw is 66.2 % over the iodine rods and 6.1 % over the calcium rods; of the noise-free reference, 8.7 % and 4.8 %.

Verification: 9 of 9 checks pass

Summary

ws    = BS.create_workspace(scanner, protocol, opts, recon_opts, phantom)
sim   = BS.simulate!(ws, phantom, protocol, opts; keep_projection = true)  # then (…; projection)
basis = BS.spectral_basis(ws; I0 = sim.I0_bins)          # ray-resolved, no calibration scan
vmi   = BS.vmi_pipeline(; channels = sim.pcct_sino.bins, basis, geom = ws.geom,
                          to_backend, matrix_size, vmi_energies = (40, 70, 100, 140),
                          denoiser = BS.SpectralHYPR(), fbp_filter = PAIR_FILTER.filter,
                          pair_basis, composite_energy)

Four photon-counting energy windows go into the decomposition as four measurements, and the spectral model is the one the simulation applied, ray by ray through the bowtie, so no calibration scan is involved. The acquisition integrates every view over five sub-views, and the noise-free reference, the calibration draw and the measured draw share one projection.

At 40 / 70 / 100 / 140 keV the solid water reads 5.9 / 1.7 / 0.6 / 0.1 HU (noise-free: 4.6 / 0.9 / -0.1 / -0.5 HU; theory 0 HU), with a noise of 42.6 / 35.3 / 33.7 / 33.1 HU per 0.4 mm slice at the physical scan's 10.12 mGy. Measured rod HU track theory at every energy (the regression slopes and R² in the verification table). The per-rod table separates each rod's error into the basis model's own bias, largest for calcium, which the iodine–water basis represents only approximately, and what the noise and the denoiser add to it.