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.

This chapter supplies the physical operators used inside Chapter 7’s joint evaluator. The existing geometric projector will integrate one field per material. spectral_signal will consume those paths, supplied attenuation tables and response-weighted photon populations. Its image feeds the detector and objective, and its reverse product supplies material-path seeds to the same projector. Keeping those boundaries explicit lets us change the energy model without rewriting 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

This is a small host-side preparation operation. The repeated image calculation receives already prepared CUDA arrays and never asks Python to reconstruct a spectrum for each pixel. Provenance records the original source identity, checksum and rights alongside its description. Those fields identify the supplied input, but they do not certify its calibration. An analytic spectrum must name its defining construction rather than borrowing the name of a real X-ray tube.

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.

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]

MaterialTable performs that preparation without bundling a material database. A mass table requires its reference density, while a linear table rejects another density factor because it already contains one. The conversion property returns inverse millimetres, the unit needed by our millimetre paths. It also rejects a conversion that overflows or turns a positive coefficient into zero.

@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)

Two equal adjacent energies encode the below-edge and above-edge values in that order. at_energies requires an explicit side for an exact-edge query and interpolates between distinct neighbouring nodes elsewhere. Its default log–log interpolation requires positive coefficients, while the linear option permits zero. Neither extrapolates beyond the supplied support. These are declared numerical interpretations of a table, not additional physical measurements, and the coefficient uncertainty remains after interpolation has been performed perfectly.

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

The slope in equation (8.6) is the surviving population’s average coefficient. Its decline with thickness follows from the coefficient variance. When attenuation decreases with energy over the relevant band, this reweighting preferentially retains higher energies: beam hardening. Absorption edges can complicate the ordering by energy, while the variance identity still holds for fixed non-negative coefficients and weights.

An exact two-energy example makes the effect visible without inventing a tube spectrum. 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 an exact transmitted mixture with the mean optical depth. The numbers define an analytic test, not tissue properties or a measured beam. At increasing thickness the lower coefficient contributes an increasing fraction of the surviving signal. Replacing the mixture by one fixed effective coefficient can match one thickness without matching 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.

The geometric composition uses the sampler and finite-ray quadrature already developed in Chapter 4. prepare_material_projection binds material-major field and output views once. project_material_paths then schedules the canonical projector for each material. The Python loop runs over the small material count, while CUDA threads still perform every ray integral. Its output layout is (material, pixel), flattened with pixels contiguous, so adjacent pixel threads read adjacent path values 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()

The adapter’s default field_domain="fractions" adopts a narrower basis than the general fields in equation (8.3): fixed-density material fractions lie in [0, 1] and sum to at most one at each voxel, and the remainder is vacuum. They use the same oriented sample-centre grid and constant extension to its half-cell support faces as the scalar projector.

The validator accumulates stored fractions in binary64 and allows an excess of at most $2^{-24}$ above one: independently rounding a nonnegative partition to binary32 can create that much excess without changing its intended occupancy. Individual fractions must still lie in [0, 1]. Larger violations are rejected, never rescaled to sum to one. Silent renormalisation would change material amounts and their derivatives while leaving the array shape reassuringly intact.

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

optical_depth combines material paths and coefficients in binary64. The exponential and the weight–response product retain that precision until the final store. Calling the scalar transmit interface on a prematurely narrowed depth would change this numerical contract. The physics is still the same exponential survival law, but its spectral inputs need the wider intermediate sum.

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. Chapter 12 compares the observations with a reconstructed prediction. Recorded 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.

The spectral_signal API returns the linear expected mean in its declared storage precision; it does not expose a spectral log-mean output. Binary64 intermediates extend the useful arithmetic range before storage, but cannot preserve a linear output below the destination’s range. The retained-binary64 mode also rejects a nonzero exponential or intermediate product that rounds to zero: a later large factor could make the full contribution representable, so silently discarding it would be wrong. Taking a logarithm of the stored zero would not implement equation (8.9). A log-domain output would have to evaluate the log sum before rounding to the linear destination.

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.

The implemented BlurSpec declares a fixed, non-negative finite stencil with total weight at most one. Its array is a matrix rule on a finite detector: samples outside the image are zero. The same kernel factory produces blur and blur_transpose, changing the offset direction while retaining the weights and boundary support.

@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

Flipping the offsets is necessary for an asymmetric stencil. Applying the forward blur to a reverse seed would give the wrong product unless the actual finite matrix happened to be symmetric. Nor do we renormalise an edge row to conceal lost signal: that would alter the forward operator and its transpose. The small stencil uses binary64 storage and compensated binary64 accumulation, because a weak tail multiplied by a bright pixel can still contribute representable signal. Tests pair the forward and transpose products with an asymmetric stencil, so symmetry cannot let an incorrect reverse rule pass unnoticed.

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.

Observation generation is a separate API. sample_poisson_counts consumes count-domain means and writes integer draws. sample_compound_poisson consumes energy-bin arrival rates and prescribed per-event scores, sums the scored draws, and can be followed by add_gaussian_read_noise. This implements independent Poisson arrivals with the declared deterministic scores, and it does not infer a charge-sharing event law from a blur kernel or from the mean response alone.

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.

Original arrays remain unchanged until their reverse consumers finish. The manual spectral VJP rejects ambient tape recording and supplies first derivatives only, and it does not differentiate energy nodes or host table interpolation. Chapter 7’s evaluator currently activates material paths to recover pose, with named gain/exposure/offset groups outside this operator. The lower-level spectrum and response derivatives are available for a separately constrained inverse model, but their presence alone does not establish identifiability.

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

The spectral calculation takes material paths, physical coefficient tables and a detector response. Listing 8.3 supplies the paths, Listing 8.4 evaluates the spectral sum and Listing 8.6 sends image derivatives back to those paths. This separation gives each stage an independent numerical reference before the stages are composed in a pose fit.

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.

experiments/spectral-projection/run.py executes this composition from supplied field, coefficient, spectrum and observation arrays. It records source and input digests with the objective and derivatives, and an explicit prediction export downloads numerical arrays after evaluation.

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.