From an X-ray spectrum to a detector measurement

A thick path removes photons unevenly across the spectrum. The photons that survive are therefore a different mixture from those that entered. Even if the object and detector remain perfectly still, changing the source spectrum changes which material differences the image emphasises. A single attenuation coefficient cannot follow that change merely by being given a more impressive variable name.

We will retain straight, uncollided paths and let photon energy vary. The geometry still supplies physical material paths. Energy-dependent attenuation acts along those paths, and the detector assigns a response to the surviving photons. This order matters: the energy integral belongs outside the exponential, and the detector response belongs inside the energy integral.

The result is an expected primary signal. Scatter requires the additional paths in Chapter 9.

Reuse the geometric projector to integrate one field per material, then pass those material paths to spectral_signal. That second operation combines them with attenuation tables and response-weighted photon populations. Its output enters the detector model and objective. In reverse, it returns material-path seeds to the projector. This is the chain needed by Chapter 7’s joint evaluator, and it lets us change the energy model without changing ray geometry.

8.1 The source spectrum and its units

Let $n_{0,p}(E)$ be the expected number of open-beam photons per unit energy associated with pixel $p$ over one exposure. We measure $E$ in keV, so this density has units photons per keV. It is defined at the detector, before detection losses. Source output, geometric acceptance and source-to-detector distance have already entered this quantity, as they did for $n_{0,p}$ in Chapter 2.

For a nonzero open-beam population, separate the total expected count $n_{0,p}$ from a normalised spectral shape $f_p(E)$:

In equation (8.1), $f_p$ has units $\mathrm{keV}^{-1}$. Multiplying exposure while preserving shape changes $n_{0,p}$, while changing filtration can alter both factors. If $n_{0,p}=0$, no normalised shape is determined by the empty spectrum. Use the zero expected signal directly rather than dividing by zero.

A spectrum supplied as photons per area, per energy and per exposure must be integrated over the appropriate detector acceptance before it becomes $n_{0,p}(E)$. A source angular distribution requires the corresponding geometric transformation. Conversely, a spectrum already reported as photons in each energy bin is a bin integral, not a density to multiply by the bin width again.

Energy support is part of the input contract. Record the source model, tube potential, filtration, energy grid and normalisation convention. An energy outside a table’s support does not acquire a physical coefficient by clipping its index to the last row. If the source has discrete lines, treat their expected populations separately or integrate them consistently into bins. A delta-function line cannot be sampled faithfully by hoping that a quadrature node lands on it.

The source shape may depend on pixel through filtration or angular emission. A single shared shape is a useful reduction when justified by the acquisition, not a property of the notation. On the GPU, a justified shared spectrum lets pixels reuse one set of spectral weights.

Spectrum records the distinction between density and integrated populations before either reaches the GPU. Its bin_fluence property supplies the expected photons associated with each discrete term. It multiplies density values by the caller’s quadrature weights, but returns bin-integrated populations unchanged. A density with no weights, or a bin population with an extra set of weights, is rejected at ingestion.

@dataclass(frozen=True, slots=True)
class Spectrum:
    """A supplied open-beam spectrum at the detector, before response losses.

    Density values are photons/keV and require positive quadrature weights in
    keV. Bin values are already integrated expected photon populations and must
    not receive a second energy weight. No normalisation or inverse-square
    correction is guessed. Discrete lines may be included as integrated bins.
    """

    energies_kev: tuple[float, ...]
    values: tuple[float, ...]
    representation: Literal["density", "bin-fluence"]
    provenance: Provenance
    quadrature_weights_kev: tuple[float, ...] | None = None

    def __post_init__(self) -> None:
        finite_tuple(self.energies_kev, "energies_kev", positive=True)
        finite_tuple(self.values, "spectrum")
        if len(self.values) != len(self.energies_kev):
            raise ContractError("spectrum values and energies must have equal length")
        if any(a >= b for a, b in zip(self.energies_kev, self.energies_kev[1:], strict=False)):
            raise ContractError("spectrum energy nodes must be strictly increasing")
        if not isinstance(self.provenance, Provenance):
            raise ContractError("spectrum provenance is required")
        if self.representation == "density":
            if self.quadrature_weights_kev is None:
                raise ContractError("a density needs explicit energy quadrature weights")
            finite_tuple(self.quadrature_weights_kev, "quadrature weights", positive=True)
            if len(self.quadrature_weights_kev) != len(self.values):
                raise ContractError("one quadrature weight is required for each density value")
        elif self.representation == "bin-fluence":
            if self.quadrature_weights_kev is not None:
                raise ContractError("bin-integrated fluence must not be weighted twice")
        else:
            raise ContractError("spectrum representation must be density or bin-fluence")

    @property
    def bin_fluence(self) -> tuple[float, ...]:
        """Materialise energy-integrated host values at the explicit ingestion boundary."""
        if self.quadrature_weights_kev is None:
            return self.values
        result = tuple(
            value * weight
            for value, weight in zip(self.values, self.quadrature_weights_kev, strict=True)
        )
        if any(not math.isfinite(value) for value in result):
            raise ContractError("energy integration overflowed")
        if any(
            source > 0 and target == 0 for source, target in zip(self.values, result, strict=True)
        ):
            raise ContractError("energy integration underflowed a positive bin population")
        return result

Prepare the spectrum once, then keep its photon populations in CUDA arrays across image evaluations. Attach the source model, units, energy grid and calibration information to the input so that a later change of spectrum has a defined physical meaning.

8.2 Energy-dependent material attenuation

Keep the physical attenuation field in the object frame. For a mixture described by elemental mass fractions $w_m(\mathbf{x})$, density $\rho(\mathbf{x})$ and elemental mass attenuation coefficients $\kappa_m(E)$, the independent-constituent mixture model is

With $\rho$ in $\mathrm{g\,mm}^{-3}$ and $\kappa_m$ in $\mathrm{mm}^2\mathrm{g}^{-1}$, equation (8.2) produces $\mu$ in $\mathrm{mm}^{-1}$. Tabulations often use $\mathrm{g\,cm}^{-3}$ and $\mathrm{cm}^2\mathrm{g}^{-1}$ instead. Their product is inverse centimetres, so divide that product by ten to obtain inverse millimetres. Convert once at a declared boundary.

Mass fractions and volume fractions are different mixture descriptions. A volume-fraction combination of reference materials uses their linear coefficients and its own density assumptions. Mixing one scheme’s fractions with another scheme’s coefficients can yield plausible positive numbers without representing the intended material. NIST’s attenuation tables use mass-fraction additivity for compounds and distinguish total attenuation from energy-absorption coefficients. That distinction fixes both the mixture weights and the quantity needed in our exponent. [8]

The coefficient here removes photons from the primary beam. It includes scattering out of the uncollided population. The mass energy-absorption coefficient answers a different question about energy deposited through secondary charged particles. It cannot replace the attenuation coefficient in a transmission exponent.

For computation, a separable basis is often convenient. Let $a_m(\mathbf{x})$ be non-negative, dimensionless material fields and let $\mu_m(E)$ be linear attenuation coefficients at declared reference densities. Define the material paths $A_{pm}$ along the finite ray $\ell_p$:

The basis in equation (8.3) separates spatial integration from energy evaluation. $A_{pm}$ has units of millimetres, and it is a literal material length when $a_m$ is a material indicator. A fitted basis field can instead be a relative concentration, so its interpretation must follow its definition. The basis is a modelling assumption: a single CT number does not uniquely specify all its coefficients, as discussed in Appendix A.11.

This separation also underlies Alvarez–Macovski spectral reconstruction: energy-dependent basis functions multiply spatial coefficients, whose line integrals determine the measurements under each spectrum. That reconstruction solves for the line integrals from spectral measurements. Here we project supplied material fields to obtain them, then evaluate the same form of spectral forward map. Choosing material reference curves as our basis makes those paths directly interpretable in the renderer. [11, 12]

Interpolation needs the same care as the field representation. Preserve non-negative coefficients and distinguish the two sides of an absorption edge. Smoothing a discontinuity into a continuous curve changes the spectral operator near that edge. Refining the energy quadrature cannot repair a coefficient table that has already removed the feature. XCOM supplies values immediately above and below constituent absorption edges, reflecting the discontinuities that a material mixture inherits. Its isolated-atom data also have limits near edges where molecular or solid-state effects matter. [9]

Construct MaterialTable from the coefficient values and their units. Supply reference density for mass attenuation coefficients, then convert to inverse millimetres before combining them with millimetre paths. A linear attenuation table already includes density. Check that unit conversion preserves finite positive values.

@dataclass(frozen=True, slots=True)
class MaterialTable:
    """Total primary attenuation at one declared material composition/density.

    Arrays are immutable host metadata, never an alternative production backend.
    Energies are keV. Duplicate adjacent energies encode below/above edge values,
    in that order; at most two entries may share an energy. The density attached
    to a mass table is in g/cm³. Linear coefficients already include density.
    """

    name: str
    energies_kev: tuple[float, ...]
    coefficients: tuple[float, ...]
    units: Literal["mm^-1", "cm^-1", "cm^2/g"]
    provenance: Provenance
    density_g_cm3: float | None = None
    interpolation: Literal["linear", "log-log"] = "log-log"

    def __post_init__(self) -> None:
        if not self.name.strip():
            raise ContractError("a material table needs an unambiguous material/composition name")
        if not isinstance(self.provenance, Provenance):
            raise ContractError("material-table provenance is required")
        finite_tuple(self.energies_kev, "energies_kev", positive=True)
        finite_tuple(self.coefficients, "coefficients")
        if len(self.energies_kev) != len(self.coefficients):
            raise ContractError("energies and coefficients must have equal lengths")
        if self.units not in ("mm^-1", "cm^-1", "cm^2/g"):
            raise ContractError("unsupported coefficient units")
        if self.interpolation not in ("linear", "log-log"):
            raise ContractError("interpolation must be linear or log-log")
        if self.interpolation == "log-log" and any(c <= 0 for c in self.coefficients):
            raise ContractError("log-log interpolation requires strictly positive coefficients")
        for index, energy in enumerate(self.energies_kev[1:], 1):
            if energy < self.energies_kev[index - 1]:
                raise ContractError("energy support must be non-decreasing")
            if index > 1 and energy == self.energies_kev[index - 2]:
                raise ContractError("an absorption edge has exactly two one-sided entries")
        if self.units == "cm^2/g":
            density = self.density_g_cm3
            density = finite_scalar(density, "reference density in g/cm³", minimum=0)
            if density <= 0:
                raise ContractError("reference density must be strictly positive")
        elif self.density_g_cm3 is not None:
            raise ContractError(
                "linear coefficients already include density; do not apply it twice"
            )

    @property
    def linear_mm_inverse(self) -> tuple[float, ...]:
        """Convert source values once, without inventing a different mixture model."""
        scale = 1.0
        if self.units == "cm^-1":
            scale = 0.1
        elif self.units == "cm^2/g":
            assert self.density_g_cm3 is not None
            scale = self.density_g_cm3 / 10.0
        converted = tuple(c * scale for c in self.coefficients)
        if any(not math.isfinite(c) for c in converted):
            raise ContractError("unit conversion overflowed; coefficient table is unusable")
        if any(
            source > 0 and target == 0
            for source, target in zip(self.coefficients, converted, strict=True)
        ):
            raise ContractError("unit conversion underflowed a positive attenuation coefficient")
        return converted

    def at_energies(
        self,
        energies_kev: tuple[float, ...],
        *,
        edge_side: Literal["below", "above"],
    ) -> tuple[float, ...]:
        """Interpolate in mm⁻¹, requiring the one-sided convention at exact edges."""
        if edge_side not in ("below", "above"):
            raise ContractError("edge_side must be explicitly below or above")
        finite_tuple(energies_kev, "query energies", positive=True)
        grid = self.energies_kev
        coefficients = self.linear_mm_inverse
        result: list[float] = []
        for energy in energies_kev:
            if energy < grid[0] or energy > grid[-1]:
                raise ContractError(f"{self.name}: {energy} keV lies outside supplied support")
            left, right = bisect_left(grid, energy), bisect_right(grid, energy)
            if left != right:
                result.append(coefficients[left if edge_side == "below" else right - 1])
                continue
            low, high = left - 1, left
            if self.interpolation == "log-log":
                fraction = _log_ratio(energy, grid[low]) / _log_ratio(grid[high], grid[low])
                value = math.exp(
                    (1.0 - fraction) * math.log(coefficients[low])
                    + fraction * math.log(coefficients[high])
                )
            else:
                fraction = (energy - grid[low]) / (grid[high] - grid[low])
                value = (1.0 - fraction) * coefficients[low] + fraction * coefficients[high]
            result.append(value)
        return tuple(result)

Represent an absorption edge with two adjacent entries at the same energy: below-edge value first, above-edge value second. at_energies then requires a side for an exact-edge query and interpolates between distinct neighbouring nodes elsewhere. Use log–log interpolation for positive coefficients or the linear option when zero is allowed, staying within the supplied energy support.

8.3 Polychromatic transmission

At each energy, primary survival is still exponential. The expected spectral population reaching the detector is $n_p(E)=n_{0,p}(E)e^{-L_p(E)}$. Let $R_p(E)$ be the expected detector output per incident photon of energy $E$, under a spatially local, linear response model. The expected primary output before an additive electronic offset is

Equation (8.4) has units of detector output: photons per keV times output per photon times keV. The ideal counter has $R_p(E)=1$. An ideal full-energy integrator expressed in keV has $R_p(E)=E$. Actual response may include detection efficiency and incomplete energy deposition. Multiplying by $E$ alone does not calibrate a real detector.

For non-negative response and a positive open-beam output, define the normalised response-weighted spectrum $\alpha_p(E)$. The measured-domain transmission and its apparent optical depth are

Equation (8.5) exposes two averages that must not be exchanged. The detector averages transmitted contributions, and it does not average optical depths and then exponentiate. For a finite positive signal, Jensen’s inequality gives $\bar L_p\leq\int\alpha_p(E)L_p(E)\,\mathrm dE$. Equality holds when optical depth is constant on the contributing energy support.

Consider a homogeneous slab of thickness $d$ with coefficient $\mu(E)$ and fixed $\alpha(E)$. Define the transmitted, normalised energy weights $\beta_d(E)$. Differentiating the log transmission gives

As the slab thickens, components with smaller attenuation coefficients contribute more of the surviving signal. The slope of its apparent optical depth is their weighted mean coefficient, and its derivative with respect to thickness is the negative coefficient variance (equation (8.6)). When attenuation decreases with energy over the relevant band, this preferential survival of higher energies is beam hardening. Absorption edges can complicate the energy ordering without invalidating the variance identity for fixed non-negative coefficients and weights.

For an analytic two-energy beam, assign equal incident counting weights and illustrative coefficients $0.02$ and $0.04\,\mathrm{mm}^{-1}$. At $d=50\,\mathrm{mm}$,

Equation (8.7) compares the transmitted mixture with the mean optical depth. As thickness increases, the lower coefficient contributes more of the surviving signal. One fixed effective coefficient can therefore match one thickness and miss the next.

Figure 8.1 shows how the transmitted mixture changes with thickness.

8.4 Energy quadrature and material paths

Assume first a continuous spectrum on a fixed interval. Choose energy nodes $E_k$ and non-negative quadrature weights $w_k$ in keV. After the geometric operator computes the material paths, the discrete spectral calculation is

The coefficient $c_{pk}$ in equation (8.8) already includes an energy weight. For a bin-integrated photon count, replace $w_k n_{0,p}(E_k)$ by that bin population and state the approximation used for response and attenuation within the bin. These alternatives must not be applied together.

Each material path is reused across energy nodes. This saves repeated volume traversal when the basis is exact for the chosen field model. It also makes the sources of approximation separable: ray quadrature determines $A_{pm}$, material representation determines which fields exist, and energy quadrature determines the final sum. Doubling all three resolutions at once can reduce error without revealing its source.

Apply Chapter 4’s sampler and finite-ray quadrature to each material field. Bind the field and output views once with prepare_material_projection, then let project_material_paths schedule a projection per material. CUDA threads perform the ray integrals. Store the results as (material, pixel) with pixels contiguous, so adjacent threads read adjacent paths for a given material.

def project_material_paths(
    pose: Any,
    *,
    workspace: MaterialProjectionWorkspace,
    stream: Any = None,
    validate: bool = True,
) -> None:
    """Overwrite material-major path lengths in mm, using one canonical projector.

    All material data stay on CUDA. The host loop schedules M operations and
    never visits a voxel, pixel or quadrature sample. Explicit VJPs compose
    with spectral_signal; ambient tape recording is rejected.
    """
    require_no_tape()
    ctx, projection = workspace.context, workspace.projection
    ctx.assert_stream(stream)
    ctx.array(pose, "pose", dtype=ctx.wp.float64, shape=(12,))
    ctx.disjoint([("pose", pose)], [("paths", workspace.paths)])
    if validate:
        workspace.validate_inputs(stream=ctx.stream)
    for material in range(workspace.materials):
        _project(
            workspace.field_views[material],
            pose,
            workspace=projection,
            out_L=workspace.path_views[material],
            stream=ctx.stream,
            validate=validate and material == 0,
        )
    if validate:
        workspace.check_status()

Choose the meaning of each input field before projecting it. With field_domain="fractions", the fields are fixed-density material fractions in [0, 1], their voxelwise sum is at most one, and any remainder is vacuum. Use the oriented sample-centre grid and constant extension to its half-cell support faces defined for the scalar projector.

Check the fraction sum in binary64. Independently rounding a nonnegative partition to binary32 can leave an excess of at most $2^{-24}$ above one, which the validator allows. Each fraction must still lie in [0, 1]. Reject larger violations instead of renormalising the fields: changing their fractions would change the material being projected.

The explicit field_domain="nonnegative" option instead accepts finite, nonnegative dimensionless equivalent-basis coefficients without upper or sum caps. Those coefficients are not volume fractions: their normalisation belongs in the supplied basis attenuation, in $\mathrm{mm}^{-1}$.

With paths available, spectral_signal launches the following pixel-owned sum. Coefficients are material-major with energy contiguous. Weights and response are either shared energy vectors or energy-major pixel arrays, and wp.static removes the unused indexing branch for the prepared layout. The number of materials and energies belongs to that specialisation, rather than varying within a warp.

@cache
def get_forward(
    materials: int, energies: int, shared_weights: bool, shared_response: bool, wide: bool = False
):
    dtype = wp.float64 if wide else wp.float32
    optical_depth = get_optical_depth(wide)
    attenuate = get_attenuation(wide)
    product = get_product(wide)
    store_signal = get_signal_store(wide)

    @wp.kernel(module="unique", module_options=STRICT)
    def forward(
        paths: wp.array(dtype=dtype),
        coefficients: wp.array(dtype=wp.float32),
        weights: wp.array(dtype=wp.float32),
        response: wp.array(dtype=wp.float32),
        pixels: int,
        mean: wp.array(dtype=dtype),
        status: wp.array(dtype=wp.int32),
    ):
        pixel = wp.tid()
        total = wp.float64(0.0)
        compensation = wp.float64(0.0)
        for energy in range(energies):
            wi = energy * pixels + pixel
            ri = wi
            if wp.static(shared_weights):
                wi = energy
            if wp.static(shared_response):
                ri = energy
            depth = optical_depth(paths, coefficients, pixel, energy, pixels, materials, energies)
            contribution = product(
                wp.float64(weights[wi]) * wp.float64(response[ri]), attenuate(depth, status), status
            )
            # All terms are non-negative, but compensation retains small bins
            # when a broad response places many decades in the same sum.
            corrected = contribution - compensation
            updated = total + corrected
            compensation = (updated - total) - corrected
            total = updated
        mean[pixel] = store_signal(total, status)

    return forward

Accumulate material paths times coefficients in binary64 inside optical_depth, and retain that precision through the exponential and weight–response product. Round only at the destination. This postpones rounding until each spectral contribution to the exponential survival law has been formed.

By default, paths and image outputs use binary32 storage. Setting SpectralSpec(precision="float64") retains binary64 paths, image means, image cotangents and path cotangents; the supplied coefficients, weights and response remain binary32. Chapter 12 uses this option because rounding a predicted count can obscure the small loss decrease near a solution. Retaining precision through the loss and its reverse products lets the line search see the improvement that the gradient predicts.

Compensated accumulation retains small energy-bin contributions when a broad spectrum places very different magnitudes in one pixel’s sum. No pixel-by-energy image is written. This saves both its storage and the later traffic needed to read it back, at the cost of doing the energy loop within each pixel thread. Register demand consequently depends on the active material count, especially in reverse execution, so a timing for one material basis cannot stand in for all bases.

Figure 8.2 applies that separation to CT-derived anatomy. The material paths and view stay fixed while three detector channels weight the transmitted spectrum differently. The material assignments and response are declared simulation inputs, not measured patient composition or scanner calibration.

20–55 keV

Negative log display: 0–8

55–80 keV

Negative log display: 0–5.5

80–120 keV

Negative log display: 0–4

120 kVp tungsten source · nominal 2.5 mm Al filtration · common view and detector samples. The recorded model omits scatter, blur, pileup and charge sharing. Download the numerical record. Data sources and display rules.

Partition the energy interval at absorption edges and response thresholds. Integrate each smooth piece with its own nodes, respecting the one-sided coefficient values. Discrete source lines contribute explicit weighted terms. A high-order rule spanning a jump does not retain its smooth-integrand error guarantee merely because it has more decimal places in its weights.

For a fixed material basis, record convergence of the signal and of the derivatives that optimisation will use. Very small absolute changes in a nearly extinguished signal can coexist with large relative error. Choose an absolute signal tolerance, a relative tolerance where the reference is sufficiently nonzero, and a separate derivative tolerance.

Strong attenuation can underflow an exponential even when the desired log signal remains representable. For strictly positive terms, put $z_{pk}=\log c_{pk}-L_{pk}$ and $m_p=\max_k z_{pk}$. Then

Equation (8.9) evaluates the same positive sum in log space. Zero-weight terms are omitted, and if all weights vanish, the expected signal is zero and its logarithm is not finite. This numerical rearrangement does not authorise inserting a positive floor into the physical measurement model.

Choose the output domain before choosing its storage. spectral_signal returns a linear mean, with binary64 intermediates and a binary32 or binary64 destination. A log-mean output would instead need to evaluate equation (8.9) before rounding to linear storage. In the retained-binary64 mode, an intermediate that underflows is rejected because a later large factor could make the full contribution representable.

8.5 What the detector measures

A detector does more than decide whether a photon arrived. Let $Q_p$ be the random output caused by one photon incident at energy $E$. Its conditional mean is $R_p(E)$. A non-detection contributes zero, so efficiency can be included without a separate correction if it is already part of this response distribution.

For a photon-counting detector with energy bins, let $P_{pb}(E)$ be the probability that an incident photon produces one accepted count in bin $b$. In an ideal exclusive-bin model,

The inequality in equation (8.10) allows missed events. Perfect energy classification is a limiting response, not a consequence of calling the instrument photon-counting. Energy redistribution broadens the response. Charge sharing or multiple threshold counts can violate the exclusive-event model and require a different joint response law.

An energy integrator accumulates deposited energy or a related charge. If every detected photon deposits its full energy and detection is Bernoulli with efficiency $\eta_p(E)$, the mean response is $\eta_p(E)E$ in keV per photon. The second moment is $\eta_p(E)E^2$, not the square of the mean. This distinction determines noise.

Saturation and pulse pile-up couple events. A response that depends on arrival rate cannot be represented by one fixed $R_p(E)$ for every exposure. Either restrict the model to a regime where independent linear response is justified, or introduce a rate-dependent detector model and validate it with suitable calibration data. Exposure scaling alone cannot identify an arbitrary spectrum and an arbitrary detector response because their product enters the mean.

Different calibration observations constrain different factors. An open-beam image constrains the integrated response-weighted spectrum. Energy-resolved measurements constrain its shape, while independent detector characterisation constrains redistribution and efficiency. An image fit should not be expected to manufacture all three from one projection.

8.6 Spatial response and measurement noise

A spatial response can distribute a photon’s output across several pixels. Let $\bar z_r$ be the local expected output before that distribution and let $H_{pr}$ map it to readout pixel $p$. With additive mean electronic offsets $b_p$,

Equation (8.11) assumes a linear, specified spatial response. A shift-invariant blur is a special case. For a lossless redistribution on a complete readout, column sums equal one, while cropped boundaries or detection losses require different accounting. An energy-dependent blur belongs inside the spectral integral, so a single post-integral matrix is justified only when that separation is valid.

Represent a fixed spatial response with BlurSpec: a non-negative finite stencil of total weight at most one, with zero outside the detector. Construct its transpose by reversing the offset direction while retaining the weights and boundary support. The shared kernel factory below generates both operations from that rule.

@cache
def get_blur(height: int, width: int, kernel_height: int, kernel_width: int, transpose: bool):
    @wp.kernel(module="unique", module_options=STRICT)
    def blur(
        source: wp.array(dtype=wp.float32),
        weights: wp.array(dtype=wp.float64),
        output: wp.array(dtype=wp.float32),
        status: wp.array(dtype=wp.int32),
    ):
        p = wp.tid()
        row = p // width
        column = p % width
        total = wp.float64(0.0)
        compensation = wp.float64(0.0)
        for kr in range(kernel_height):
            for kc in range(kernel_width):
                dr = kr - kernel_height // 2
                dc = kc - kernel_width // 2
                if wp.static(transpose):
                    dr = -dr
                    dc = -dc
                sr = row + dr
                sc = column + dc
                # Zero extension loses signal crossing the finite detector edge.
                # Do not renormalise a boundary row: that changes both B and Bᵀ.
                if sr >= 0 and sr < height and sc >= 0 and sc < width:
                    term = wp.float64(weights[kr * kernel_width + kc]) * wp.float64(
                        source[sr * width + sc]
                    )
                    corrected = term - compensation
                    updated = total + corrected
                    compensation = (updated - total) - corrected
                    total = updated
        output[p] = checked_store(total, status)

    return blur

Use an asymmetric stencil when checking the forward–transpose pairing. A symmetric example can hide an accidental reuse of the forward operation. Keep the finite detector boundaries in both directions, including signal lost outside the image. Binary64 stencil storage and compensated accumulation retain weak tails whose products with bright pixels remain significant.

Mean blur alone does not determine noise correlations. Model a Poisson population of independent photon events, and let $Q_p(E)$ be one event’s output in pixel $p$. For a single incident stream with spectral intensity $n(E)$, compound-Poisson statistics give

In equation (8.12), $n(E)$ is the spectral population reaching this detector response, already attenuated. Independent incident streams contribute sums of these moments. The covariance uses a raw second moment because the number of events is itself Poisson. Independent electronic noise adds its covariance, and it does not change that photon term.

For a local ideal counter, $Q=1$ for every arrival, and mean equals variance. For a full-energy integrator, $Q=E$, giving variance $\int n(E)E^2\,\mathrm dE$ in squared energy units. With two arrival energies, the mean weights them by energy and the variance by energy squared. A Poisson likelihood written directly for deposited keV would confuse these two quantities.

Exclusive assignment of independent Poisson events to different pixels can yield independent thinned counts. Sharing one event’s charge between pixels instead produces joint contributions and generally nonzero covariance. Both mechanisms may produce a blurred expected image, but their noise is different. A processed image can add further dependencies through interpolation, subtraction or temporal filtering.

Figure 8.3 shows why a mean blur does not determine the noise covariance.

Only after choosing this observation model should we choose a discrepancy. Independent Poisson bins justify a Poisson count likelihood under the corresponding assumptions. A Gaussian approximation needs an appropriate count regime and covariance model. Applying a negative logarithm changes both the noise distribution and its sensitivity to calibration, especially near zero signal.

Generate observations after evaluating the expected image. For counts, sample_poisson_counts draws integers from count-domain means. For fixed per-event scores, sample_compound_poisson draws arrivals in energy bins and sums their contributions. Add electronic noise with add_gaussian_read_noise when required. A charge-sharing detector needs its own joint event response, as derived above.

ObservationIdentity attaches randomness to an observation and global pixel, with separate domains for energy bins and read noise. Keeping that identity reproduces a draw when an image is processed in different chunks, while a fresh observation requires a fresh identity. energy_offset selects the first energy-bin counter domain, so chunks can share the same observation without reusing a bin.

The Poisson sampler uses inversion at small means and rejection at larger means, with binary64 acceptance arithmetic. draw_budget limits random rejection proposals, not the sampled photon count. Small-mean inversion uses one proposal and a separate bounded recurrence.

Exhausting either numerical bound marks the whole realisation invalid. The zero stored on that failure path is a placeholder for a failed calculation, never an observed absence of photons.

8.7 Differentiate the spectral image

Suppose a scalar parameter $\theta$ changes the material paths, spectrum or response. For fixed energy limits and sufficient regularity to exchange derivative and integral, differentiation of the primary signal gives

Equation (8.13) separates changes in the incident/response weighting from changes in attenuation. For pose alone, with a fixed acquisition model, the first term vanishes. If an energy endpoint or threshold moves with the parameter, boundary contributions may be required, and the fixed-limit expression does not silently include them.

With fixed basis spectra, the derivative with respect to a material path has a particularly useful sign:

For non-negative inputs, equation (8.14) is non-positive. A positive derivative under those assumptions is a sign or modelling error. Its units are output per millimetre. A pose derivative then contracts these path sensitivities with the derivatives of $A_{pm}$ with respect to the chosen translation or rotation coordinates.

For the discrete operator, reverse accumulation can avoid storing a Jacobian indexed by every pixel and energy. If $\gamma_p=\partial\mathcal J/\partial\bar y_p^{(K)}$ is the incoming objective adjoint, the material-path adjoint is

Equation (8.15) is the reverse rule for the implemented spectral sum. It still needs the geometric operator’s derivative to reach pose. Recomputing the exponentials in reverse trades arithmetic for saved intermediates, while retaining them trades memory traffic for arithmetic. The derivative is unchanged only if recomputation uses the same inputs and numerical conventions.

spectral_vjp consumes the objective’s pixel seed and writes whichever declared input derivatives the caller requests. The pixel kernel below recomputes each energy depth from the original paths and coefficients. Material-path derivatives accumulate in compensated binary64 registers and are stored only after the energy loop. Per-pixel weight or response derivatives can be written directly by their owning pixel thread.

@cache
def get_pixel_vjp(
    materials: int,
    energies: int,
    shared_weights: bool,
    shared_response: bool,
    write_paths: bool,
    write_weights: bool,
    write_response: bool,
    wide: bool = False,
):
    dtype = wp.float64 if wide else wp.float32
    optical_depth = get_optical_depth(wide)
    attenuate = get_attenuation(wide)
    product = get_product(wide)
    store_signal = get_signal_store(wide)
    gradient_vector = wp.types.vector(length=materials, dtype=wp.float64)

    @wp.kernel(module="unique", module_options=STRICT)
    def vjp(
        paths: wp.array(dtype=dtype),
        coefficients: wp.array(dtype=wp.float32),
        weights: wp.array(dtype=wp.float32),
        response: wp.array(dtype=wp.float32),
        seed: wp.array(dtype=dtype),
        pixels: int,
        grad_paths: wp.array(dtype=dtype),
        grad_weights: wp.array(dtype=wp.float32),
        grad_response: wp.array(dtype=wp.float32),
        status: wp.array(dtype=wp.int32),
    ):
        pixel = wp.tid()
        path_gradient = gradient_vector()
        path_compensation = gradient_vector()
        for energy in range(energies):
            wi = energy * pixels + pixel
            ri = wi
            if wp.static(shared_weights):
                wi = energy
            if wp.static(shared_response):
                ri = energy
            depth = optical_depth(paths, coefficients, pixel, energy, pixels, materials, energies)
            weighted_seed = product(wp.float64(seed[pixel]), attenuate(depth, status), status)
            if wp.static(write_weights):
                grad_weights[wi] = checked_store(
                    product(weighted_seed, wp.float64(response[ri]), status), status
                )
            if wp.static(write_response):
                grad_response[ri] = checked_store(
                    product(weighted_seed, wp.float64(weights[wi]), status), status
                )
            if wp.static(write_paths):
                common = product(
                    product(weighted_seed, wp.float64(weights[wi]), status),
                    wp.float64(response[ri]),
                    status,
                )
                for material in range(materials):
                    term = product(
                        -common, wp.float64(coefficients[material * energies + energy]), status
                    )
                    corrected = term - path_compensation[material]
                    updated = path_gradient[material] + corrected
                    path_compensation[material] = (updated - path_gradient[material]) - corrected
                    path_gradient[material] = updated
        if wp.static(write_paths):
            for material in range(materials):
                grad_paths[material * pixels + pixel] = store_signal(
                    path_gradient[material], status
                )

    return vjp

Shared parameters need a further sum over pixels. Separate bounded reductions handle shared spectrum, response and coefficient derivatives, and they do not let every pixel atomically update one scalar. The coefficient reduction reuses the computed energy depth across material parameters while preserving each parameter’s accumulation order. wp.static removes inactive output branches, so fitting pose alone does not also compute an unused calibration Jacobian.

Keep the original arrays unchanged until their reverse consumers finish. Wire the spectral material-path derivatives into the projector VJP, and handle gain, exposure and offset in the calibration operation. The manual spectral VJP computes first derivatives at fixed energy nodes and fixed table interpolation. It runs outside ambient tape recording. Spectrum and response derivatives can use the same interface when the inverse model includes those parameters.

Normalisation can carry its own derivative. For a parameter-dependent open beam and positive signals,

The second term in equation (8.16) disappears only when the reference is held fixed. Jointly fitting exposure while differentiating a normalised image as though the open beam were constant creates the wrong objective derivative.

Parameter constraints also affect identifiability. One scalar spectral measurement per ray generally cannot determine several independently unknown material paths: its local Jacobian has one row. Multiple independent spectra or energy bins can provide additional rows, but full rank and useful conditioning still depend on their responses and the object. Merely creating more parameters in the optimiser does not create another measurement.

8.8 Implement and validate the extended operator

Assemble and check the spectral renderer in the same order as its dependencies. First compare the material paths from Listing 8.3 with known slab lengths. Feed those paths to Listing 8.4 and check the one- and two-energy signals. Then use Listing 8.6 to return an image seed to the paths, and compare its derivative with the analytic slab derivative. Only then attach the pose derivative and a numerical optimiser.

A ray-owned kernel can keep a small set of material paths in registers and traverse energy nodes locally. If the material count becomes large, those registers can spill. An alternative stages paths in device memory and evaluates energies in a separate kernel. That exposes parallelism across energy but can increase intermediate storage and require reductions. Choose after measuring the intended material and energy dimensions.

Keep tables on the device across optimisation iterations. Record coefficient layout, sharing of source weights, launch count, register use, spills, occupancy and actual memory traffic. Avoid host copies in the objective loop. If reverse execution uses retained optical depths, record the pixel-by-energy allocation explicitly, and if it recomputes them, measure both passes. Table 8.1 gives independent limiting cases for checking the extended operator before using it in a fit.

Use analytic cases before comparison with independent physical software, then match material composition, source spectrum, interaction definitions and detector response for that comparison. Convergence towards the same discrete answer is weaker than agreement with a separately specified reference. Keep physical-table uncertainty separate from numerical integration error.

For a monochromatic slab the derivative reduces exactly to $-n_0R\mu e^{-\mu d}$. For the two-energy example it is the weighted sum of the two such terms. These references can be computed without calling the spectral quadrature being tested. Run step-size sweeps in material paths and pose, respecting non-negativity at boundaries and the declared rigid-update convention.

Exercise the assembled calculation with supplied field, coefficient, spectrum and observation arrays. Inspect the expected image and its derivatives before attempting a pose fit.

Photons that scatter into the detector follow paths absent from the spectral primary calculation. Including those contributions requires the event-by-event path model in Chapter 9.