Multispectral forward projection

In computed tomography (CT), especially when modeling polychromatic X-ray sources, the forward projection step can become a major computational bottleneck — particularly in iterative reconstruction (IR). A naïve approach loops over energy bins inside the ray traversal, increasing the cost proporionally to the number of spectral bins. But there’s a faster way: factor the geometry (ray marching) from the spectral mixing. This lets you compute the expensive line integrals once, then collapse them into a spectrum in a cheap per-pixel step. This method is well known in spectral CT literature, but often buried in MBIR math. This post walks through the concept, math, and a GPU-friendly implementation outline.

#ct #spectral #multispectral #forward-projection #gpu #iterative-reconstruction

Why polychromatic modeling matters

If we do not take into account the polychromatic nature of X-ray sources, we face several issues:

Beam hardening artifacts: cupping, streaks, HU bias.
Quantitative accuracy: important for material decomposition, dual-energy, and spectral CT.
Robustness: less dependent on empirical beam-hardening corrections.

But accuracy comes at a cost: polychromatic forward models require integrating over energy, which in naïve form is expensive.

The naïve way

For each ray:

Loop over K energy bins.
Traverse the ray through the volume for each bin.
Accumulate attenuation for that bin.
Sum over bins for the detector reading.

Cost: O(K × Nray samples) operations.

The factorization trick

Instead of re-traversing the volume for each energy, use a material basis approach. We assume the object is represented as a set of $M$ basis materials with attenuation spectra $μ_{m} (E)$ . For a ray $r$ , the monochromatic projection at energy $E$ is:

p (E) = \sum_{m = 1}^{M} L_{m} μ_{m} (E)

where:

$L_{m}$ is the total path length weighted by the local material density through material $m$ for that ray.
$μ_{m} (E)$ is the attenuation coefficient of material $m$ at energy $E$ .

The detected intensity for a polychromatic source spectrum $S (E)$ is:

I = \int S (E) \exp (- p (E)) d E

where:

$S (E)$ is the source spectrum (e.g., X-ray tube output) premultiplied by the detector response.

As an extension, the ray march can be performed for all materials at the same time, allowing reusage of the interpolation metrics.

Decoupled algorithm

Instead of computing $p (E)$ for each voxel during traversal, we accumulate $L_{m}$ once, then evaluate $I$ over the spectrum in a short loop.

Pseudocode:

cpp

// Step 1: Geometry traversal (expensive part)
float L[M] = {0};
for (sample : ray_samples) {
    int mat = classify_material(sample.position);
    L[mat] += sample.length;
}

// Step 2: Spectral evaluation (cheap part)
float I = 0.0f;
for (e = 0; e < num_energy_bins; ++e) {
    float p = 0.0f;
    for (m = 0; m < M; ++m)
        p += L[m] * mu[m][e];
    I += S[e] * exp(-p);
}

Performance

Traversal cost: unchanged compared to a monochromatic projector, scales with the number of materials.
Spectral cost: $O (M \cdot N_{E})$ operations per ray, where $M$ is typically small (2–4) and $N_{E}$ is the number of energy bins (often 5–20).
On a GPU, this extra cost is negligible compared to memory fetch latency from the 3D texture.

Limitations

Requires a basis-material decomposition of the volume (not just HU values). Nevertheless, can also be used as part of an iterative algorithm without a full decomposition, by updating the basis weights during optimization.
Assumes negligible scatter and ignores detector spectral response — those need separate modeling.
Accuracy depends on energy bin resolution and the quality of $μ_{m} (E)$ data.

Conclusion

For applications moving toward dual-layer or spectral CT, this decoupled approach gives most of the benefit of full multispectral forward projection without the prohibitive cost of spectral sampling inside the traversal loop. It’s especially appealing for iterative reconstruction, where the forward projector is called multiple times.

Multispectral forward projection ​

Why polychromatic modeling matters ​

The naïve way ​

The factorization trick ​

Decoupled algorithm ​

Performance ​

Limitations ​

Conclusion ​

Literature ​