Derivatives of the projection, from voxels to pose

Now we move the estimated vertebra a fraction of a millimetre. Some rays acquire more bone, some lose it, and some barely change. A registration method needs to know how each input changes its mismatch with the measured image (this is the heart of differentiability). Rendering both poses gives a finite difference, and a derivative computes that local response without requiring a second render for every parameter (this is the heart of differentiability without going bankrupt from GPU bills).

The exponential is the easy part. Chapter 2 already gave its derivative. Most of the work lies upstream: the pose changes sample locations, the interpolant changes beneath them, and clipping can move the integration interval. A plausible gradient image can omit one of those dependencies just as a plausible projection can contain a half-voxel error. The forward model is still the primary, fixed-energy measurement from the preceding chapters.

The recorded pelvic example below shows the projection and its six local count sensitivities. The anatomy is synthetic, the beam is fixed, and each displayed pose selects a separately recorded calculation. Translation and rotation retain their different derivative units. The image alone is not a derivative check.

Pelvic anatomyx right · y anterior · z superior (RAS)

The 3D anatomy loads as this figure comes into view. The recorded images are available below.

TranslateRotateReset

t_x

0 mm

t_y

0 mm

t_z

0 mm

ω_x

0 rad

ω_y

0 rad

ω_z

0 rad

Drag an axis handle. Each move uses one recorded coordinate. The other coordinates return to zero.

Rendered projection

Δ

Logarithmic count display · 0–1000 counts

Reference pose

Translation

counts/mm

Translation x (right)

Translation y (anterior)

Translation z (superior)

−252.80+252.8

Symmetric asinh colour scale

Rotation

counts/rad

Rotation about x

Rotation about y

Rotation about z

−352500+35250

Symmetric asinh colour scale

The six derivatives are evaluated at the reference pose, including when a perturbed projection is selected.

Column Row zero-based detector indices

Projectioncounts

—

t_x counts/mm

—

t_y counts/mm

—

t_z counts/mm

—

ω_x counts/rad

—

ω_y counts/rad

—

ω_z counts/rad

—

Pixel values load from the recorded execution data.

Independent finite-difference check status loads alongside the recorded samples.

What is shown, and where the data came from

The CT was generated using MAISI [36] and the bone mask was enhanced for visual clarity.

The illustrative 80 keV water-equivalent model converts synthetic CT values to attenuation with μ = 0.01837 × max(0, 1 + HU/1000) mm⁻¹. A primary-only cone beam and ideal count response use 1000 expected open-beam counts per pixel. The source is 750 mm anterior to the sacral pivot and the detector is 450 mm posterior to it. Scatter and spectral effects are omitted. The pivot and RAS coordinates are fixed across all recorded poses. Rotational derivatives use radians.

Display transforms preserve the stored values used by the pixel inspector. The grey level is log1p(F) / log1p(1000). For a family limit L, the signed colour coordinate is asinh(J / (0.02 L)) / asinh(50). The translation and rotation limits are separate. No values are clipped.

Drag an arrow to translate the pelvis or a ring to rotate it. Movements snap to the recorded poses. Translations span -2 to 2 mm in 0.5 mm increments, and rotations span -0.05 to 0.05 rad in 0.01 rad increments. Each move uses one coordinate at a time, while the other coordinates return to zero. The camera stays fixed. Image right is model right and the top is superior. The full CT supplies the projection. The 3D view shows its pelvic bone surfaces.

The Δ toggle subtracts the previous displayed projection from the current one, in counts. Its blue–rust colour scale is symmetric about zero, with limits set by the largest absolute difference in that pair. Before the first movement, the difference is zero. The pixel inspector reports the displayed counts or difference. The six derivatives stay at the reference pose. Switching the toggle does not change the comparison frame.

Projection and sensitivity colours use the recorded display scales. The browser selects recorded projections, subtracts their stored counts in Δ mode and transforms the supplied surfaces. It does not calculate photon transport. Recorded CUDA calculation and validation. Bone-label source: TotalSegmentator v2, CC BY 4.0.

Downloads: static figure · figure metadata · recorded float32 values.

We now extend the existing projector with dpt.projection.projection_vjp. It receives the optical-depth cotangent produced by transmission_vjp and returns derivatives for the selected inputs. Chapter 6’s pose evaluator uses this operation to reduce the image’s contribution to six pose derivatives, without retaining a full image Jacobian.

5.1 Which projection derivative do we need?

Let $\boldsymbol{\theta}$ collect the active inputs: perhaps the six local pose coordinates, perhaps volume coefficients, perhaps both. Hold the other inputs fixed. The renderer produces a vector of expected primary counts, and a scalar objective compares that prediction with observed data $\mathbf{y}$:

Equation (5.1) separates the renderer’s output from the comparison applied to it. Differentiating an expected count does not require sampling photons. A noise model matters when we choose $\mathcal D$ and interpret the observations.

If there are $P$ detector samples and $K$ active scalar parameters, the full Jacobian has $P K$ entries. Two products are usually more useful:

In equation (5.2), a Jacobian-vector product, or JVP, predicts the image change along the parameter direction $\boldsymbol{\eta}$. A vector-Jacobian product, or VJP, carries the detector weights $\overline{\boldsymbol{\lambda}}$ back to the inputs. In this chapter’s reverse rules, an overbar denotes a cotangent: the sensitivity of a downstream scalar to that variable.

For the objective above, choose $\overline{\boldsymbol{\lambda}}=\nabla_{\boldsymbol{\lambda}}\mathcal D$. Then $\mathbf g$ is the parameter gradient. Reverse mode is attractive when one scalar depends on many inputs. With only six pose parameters, a small set of directional derivatives can also be practical. The mathematical requirement is the same product, regardless of the differentiation strategy.

The implemented API computes the VJP. adj_L is a caller-owned detector array whose precision matches the prepared optical-depth output: binary32 by default, or binary64 when explicitly requested. The requested destinations identify the derivative coordinates: out_pose holds six binary64 cotangents for a local right increment, while out_mu holds the flat binary32 coefficient cotangents. A separate out_pose_matrix destination contains twelve unconstrained storage partials for a tape chain rule. The six-coordinate result describes changes along valid local rigid motions. A tape differentiating the construction of the stored transform instead needs partials with respect to its entries, which it can pass through that construction to the upstream parameters. Those twelve entries are not six pose parameters with spare capacity. Requesting both pose representations in one call is rejected, and grid or acquisition calibration is fixed in this projector.

A derivative with respect to translation has count per millimetre units, while one with respect to rotation has count per radian units. To compare or optimise scaled coordinates, write a local increment as

Here $\mathbf z$ is dimensionless and the diagonal matrix $\mathbf D$ supplies declared physical step scales. Equation (5.3) changes the coordinates of the derivative, but it does not add information to the image. Pose increments still compose with the current transform on the side chosen in Chapter 3.

5.2 Attenuation and transmission derivatives

For one detector sample, the expected primary count is $\lambda_p=n_{0,p}\exp(-\widehat L_p)$. Let both optical depth and open-beam expectation be active. Their differential is

Equation (5.4) supplies two distinct paths for an image change. Increasing illumination raises expected counts, while increasing attenuation lowers them. Holding the open beam fixed removes the first term. If it is itself calibrated by active acquisition parameters, that dependency remains part of the derivative.

For a fixed ray and sample pattern, Chapter 4’s linear optical-depth operator gives

The coefficient $W_{pv}$ in equation (5.5) has units of millimetres and includes both interpolation weights and physical path weights. It is zero when the ray’s sampled basis support never reaches coefficient $v$. Increasing a nonnegative coefficient cannot increase the expected primary count for this field model.

The corresponding reverse rules are

Equation (5.6) turns an image-space cotangent into an attenuation-field gradient. If one open-beam scalar is shared across pixels, its cotangent is the sum of the per-pixel contributions. A per-pixel illumination array has a different derivative shape, even if all its values happen to be equal.

Strong attenuation makes $\lambda_p$ small. For a fixed, bounded count-space seed, its contribution to $\overline L_p$ then becomes small too. The objective can alter that seed substantially. With known positive $n_{0,p}$, the Poisson negative log-likelihood, after dropping terms independent of $L_p$, has

Equation (5.7) can be evaluated without forming $1-y_p/\lambda_p$, which becomes numerically troublesome when the predicted count is tiny. Its shortened objective assumes fixed illumination. Fitting $n_{0,p}$ requires restoring the terms that depend on it. This is the same statistical distinction used in Chapter 2.

Derivatives must also retain the forward numerical policy. Computing a weighted derivative from an already-underflowed stored transmission can lose a representable result when the cotangent is large. Recompute the required exponential product with the declared precision and range, as for expected counts. The derivative of the intended smooth formula is the target, with floating-point error measured against it.

This is why the composition calls the existing transmission VJP with the original optical depths and illumination. Its optical-depth destination becomes adj_L for the projector. In a pose-only fit, the volume remains fixed and the open beam is a declared scalar, so the only shared reduction still needed is the six-coordinate pose gradient. If a later problem activates coefficients, the same projection traversal can request their gradient too.

For a Poisson fit, dpt’s recovery adapter currently composes these operators through a stored count-space seed. It rejects a positive observed count when the stored predicted mean is zero, and rejects an unrepresentable seed. Equation (5.7) permits a separately composed log-domain Poisson path for more extreme ranges, but the current adapter does not evaluate that path. Its intermediate count and seed arrays limit its numerical domain.

5.3 Geometry through the sampled field

A pose update moves the ray in object coordinates while the volume coefficients stay attached to the anatomy. For a fixed world point $\mathbf{x}^W$, let $\mathbf{x}^O=\mathbf{R}^{\mathsf T}(\mathbf{x}^W-\mathbf{t})$. Inverting the left and right updates from Chapter 3 gives the first-order motions

The minus signs in equation (5.8) arise because we query the inverse pose. Translating the object to the right makes a fixed world ray visit object points farther to the left. These expressions differ from the active point motions in Chapter 3 because the point being held fixed is different.

Within an open trilinear interpolation cell, coefficient differences give the gradient with respect to grid coordinates. For example,

Equation (5.9) differentiates the interpolation weights while the selected cell stays fixed. The other components follow by exchanging directions. The zero-extended reference uses zero outside coefficients. The implemented half-cell sampler instead suppresses the slope along a clamped edge axis, and the interior formula is unchanged. Multiplying by $\mathbf{Q}\mathbf{S}^{-1}$ converts this grid gradient into $\nabla_{\mathbf{x}^O}\mu_A$ with units $\mathrm{mm}^{-2}$.

Now let one scalar parameter be $\theta$. For a fixed sample count in a region where the clipping and interpolation branches are differentiable, the discrete optical-depth derivative is

Equation (5.10) displays both ways geometry changes a quadrature: it changes where we sample and how much physical distance each sample represents. This expression holds coefficients fixed. Active coefficient changes add the linear contribution from equation (5.5).

For samples interpolated between object-space endpoints, their derivative is

The last term in equation (5.11) vanishes for a sample parameter fixed on the full source-to-detector segment. It generally remains when sample positions are redistributed over a pose-dependent clipped interval. Detaching the clipping result from automatic differentiation discards that motion.

In the dpt library, finite_interval() from Listing 4.3 returns the derivatives of whichever box faces currently determine entry and exit. A bound still fixed at the source or detector has zero derivative with respect to object motion. The VJP can therefore distinguish a sample moving because the inverse pose changed from a sample moving because its fraction of the clipped interval changed. It needs both effects even when every sampled coefficient happens to have the same value.

For equal midpoint quadrature with fixed $M_p$, write $c_r=(r+1/2)/M_p$. Differentiating the actual node and weight definitions gives

Equation (5.12) makes the clipping dependency explicit (Figure 5.2). Rigid object motion preserves the source-to-detector distance, so $\partial_\theta d_p=0$ for those parameters. Moving source or detector parameters can change it. The active box face determines the derivative of each clipped endpoint away from ties.

5.4 Discretise and differentiate

We have differentiated a finite weighted sum. The derivative of a continuous line integral provides another useful comparison, but it has its own hypotheses. Suppose $f(t,\theta)$ is differentiable on the changing interval, with enough regularity to differentiate under the integral. Then

For our ray, $f=d_p\mu_A(\mathbf{x}_p^O(t))$. The endpoint terms in equation (5.13) account for changing integration bounds. At an exterior support face of the continuous zero-extended hat field, the field value is zero, so that face’s continuous endpoint term vanishes. Differentiating a finite midpoint sum still requires its moving weights and sample positions: those terms combine to approximate the continuous derivative.

The implemented clamped field can have a nonzero trace at its outer half-cell face. Its continuous endpoint term then remains, as does the corresponding contribution from the changing finite quadrature. This is one reason to keep the zero-extended reference and the executable field distinct: importing the reference’s vanishing endpoint term would remove a real sensitivity from the code’s declared model.

A sharp boundary shows why the distinction matters. Consider a one-dimensional material of coefficient $\mu_0>0$ occupying $[a,b]$ inside a fixed ray segment, with $b$ fixed and $a$ moving. Its integral is

The boundary derivative in equation (5.14) is nonzero even though the coefficient is constant everywhere inside the material. Sampling its hard indicator at a fixed finite set of points produces a staircase as $a$ moves: automatic differentiation of each selected branch returns zero between sample crossings. That is a derivative of the sampled staircase, and it does not approximate the moving-boundary derivative pointwise by merely increasing the number of samples.

Figure 5.3 separates the derivative of the continuous interval from that of the sampled staircase.

Exact cell traversal represents boundary motion through changing segment lengths. A continuous interpolant represents a softened boundary through spatial gradients over a finite width. Each can support a useful derivative, but they differentiate different finite-resolution fields. The boundary convention changes the gradient seen by the inverse problem.

DiffDRR expresses ray-plane intersections, sorting and length-weighted sums as tensor operations, allowing PyTorch automatic differentiation to propagate geometric sensitivities through its projection. It provides a useful comparison with our sampled trilinear operator: one exposes lengths through intersected cells, the other exposes basis weights at quadrature nodes. [22]

The sample-count rule adds another discrete decision. If $M_p=\lceil\ell_p/h\rceil$, changing pose can change both the number and positions of samples. At the transition, two different midpoint rules need not produce the same value. Within a fixed-count branch, equation (5.10) differentiates the selected quadrature. It does not describe a derivative across that transition.

One option is to keep sample counts fixed during each local optimisation stage and refine the discretisation between stages. Another is to use fixed sample parameters on the complete finite ray, accepting extra samples in zero-valued regions. Either choice has computational consequences. What matters for gradient verification is that the perturbed forward evaluations and the derivative use the same declared discretisation.

5.5 Reverse-mode differentiation in Warp

First obtain the cotangent of each optical depth from the detector objective. Then revisit the contributing samples. Each sample sends a contribution to the volume coefficients it read and, when geometry is active, a contribution through its position and weight to the pose parameters.

For one trilinear sample, with basis value $B_v(\mathbf{x}_{pr}^O)$ for coefficient $v$, the local reverse contributions are

Equation (5.15) is an execution guide for a VJP. Many rays can contribute to one volume coefficient, so those writes require a correct accumulation strategy. Pose-only differentiation avoids the volume-gradient scatter, but still needs to reduce per-ray contributions into a small shared parameter vector.

With $P=1024^2$ rays and $M=512$ samples per ray, storing only three binary32 position coordinates per sample uses 6 GiB. Coefficient indices, interpolation weights and their adjoints add to that total. Reconstructing positions from the endpoints and sample index can replace that retained tensor with arithmetic in the backward pass.

A recomputed backward must reproduce the forward’s clipping, interpolation, sampling and arithmetic policies. If its branch selection differs because inputs were mutated or precision changed, it can scatter cotangents into different coefficients. Retain the original inputs and any branch decisions needed for exact replay, or prove that recomputation obtains them under the recorded configuration.

Warp 1.17 documents that dynamic loops are not automatically replayed or unrolled in the backward pass. A generated adjoint of a ray loop therefore needs an explicit account of its intermediate values. The loop counter, selected interpolation cell and any intermediate consumed after the loop must have the values required by the backward. Trace each of those reads to retained storage or a recomputation. A small differentiable sample function can expose the local calculation, and a custom VJP using equation (5.15) can organise the replay traversal explicitly.

The custom kernel reconstructs the same clipped midpoint sequence as the forward operator. Its accumulators separate changes to the ray origin and direction from changes to the two interval bounds. lower_seed and upper_seed combine the field sum with motion of the sample positions. Their opposite signs include the loss or gain of physical interval length. The transpose of object_to_grid then converts these grid-coordinate cotangents back to physical object coordinates before the final right-pose contraction.

@cache
def get_vjp(
    active_volume: bool,
    matrix_gradient: bool,
    active_pose: bool,
    *,
    per_ray: bool = False,
    double_seed: bool = False,
    cell_gauss: bool = False,
):
    if per_ray and (active_volume or matrix_gradient or not active_pose):
        raise ValueError("per-ray diagnostics require only the six local pose coordinates")
    components = 12 if matrix_gradient else 6
    seed_dtype = wp.float64 if double_seed else wp.float32

    @wp.kernel(module="unique", module_options=OPTIONS)
    def vjp(
        field: wp.array(dtype=wp.float32),
        pose: wp.array(dtype=wp.float64),
        config: Configuration,
        seeds: wp.array(dtype=seed_dtype),
        volume_gradient: wp.array(dtype=wp.float32),
        partials: wp.array(dtype=wp.float64),
        status: wp.array(dtype=wp.int32),
    ):
        block, lane = wp.tid()
        pixel = block * BLOCK + lane
        gradient = Vec12d()
        if pixel < config.pixels:
            ray = ray_for_pixel(config, pose, pixel)
            if not valid_ray(ray):
                wp.atomic_or(status, 0, 2)
            else:
                interval = finite_interval(ray.origin, ray.direction, config.shape)
                span = interval.upper - interval.lower
                if span > wp.float64(0.0):
                    origin_gradient = wp.vec3d()
                    direction_gradient = wp.vec3d()
                    factor = wp.float64(seeds[pixel]) * ray.length
                    if wp.static(cell_gauss):
                        traversal = begin_cells(ray, interval, config.shape)
                        start = interval.lower
                        segments = wp.int64(0)
                        limit = (
                            wp.int64(config.shape[0])
                            + wp.int64(config.shape[1])
                            + wp.int64(config.shape[2])
                            + wp.int64(1)
                        )
                        while start < interval.upper and segments < limit:
                            end = wp.min(
                                traversal.next[0], wp.min(traversal.next[1], traversal.next[2])
                            )
                            if end <= start:
                                wp.atomic_or(status, 0, 2)
                                break
                            radius = (end - start) * wp.float64(0.5)
                            centre = start + radius
                            offset = radius / wp.sqrt(wp.float64(3.0))
                            for node in range(2):
                                position = centre + wp.float64(2 * node - 1) * offset
                                point = ray.origin + position * ray.direction
                                if wp.static(active_pose):
                                    value = sample_field(field, point, config.shape)
                                    origin_gradient += radius * value.gradient
                                    direction_gradient += radius * position * value.gradient
                                if wp.static(active_volume):
                                    scatter_field(
                                        volume_gradient,
                                        point,
                                        config.shape,
                                        factor * radius,
                                        status,
                                    )
                            start = end
                            traversal = advance_cells(
                                traversal, ray, config.shape, end, interval.upper
                            )
                            segments += wp.int64(1)
                        if start < interval.upper:
                            wp.atomic_or(status, 0, 2)
                        if wp.static(active_pose):
                            # Field continuity cancels all internal moving-cell terms.
                            # The zero-extended support still has external endpoint terms.
                            lower_value = support_value(
                                field, ray.origin + interval.lower * ray.direction, config.shape
                            )
                            upper_value = support_value(
                                field, ray.origin + interval.upper * ray.direction, config.shape
                            )
                            origin_gradient += (
                                upper_value * interval.upper_origin_gradient
                                - lower_value * interval.lower_origin_gradient
                            )
                            direction_gradient += (
                                upper_value * interval.upper_direction_gradient
                                - lower_value * interval.lower_direction_gradient
                            )
                    else:
                        mean = wp.float64(0.0)
                        direct_origin = wp.vec3d()
                        direct_direction = wp.vec3d()
                        shift_lower = wp.float64(0.0)
                        shift_upper = wp.float64(0.0)
                        for sample in range(config.samples):
                            alpha = (wp.float64(sample) + wp.float64(0.5)) / wp.float64(
                                config.samples
                            )
                            position = interval.lower + alpha * span
                            point = ray.origin + position * ray.direction
                            if not finite_point(point):
                                wp.atomic_or(status, 0, 2)
                            if wp.static(active_pose):
                                value = sample_field(field, point, config.shape)
                                mean += value.value
                                direct_origin += value.gradient
                                direct_direction += position * value.gradient
                                along = wp.dot(value.gradient, ray.direction)
                                shift_lower += (wp.float64(1.0) - alpha) * along
                                shift_upper += alpha * along
                            if wp.static(active_volume):
                                scatter_field(
                                    volume_gradient,
                                    point,
                                    config.shape,
                                    wp.float64(seeds[pixel])
                                    * ray.length
                                    * span
                                    / wp.float64(config.samples),
                                    status,
                                )
                        if wp.static(active_pose):
                            lower_seed = -mean + span * shift_lower
                            upper_seed = mean + span * shift_upper
                            origin_gradient = (
                                span * direct_origin
                                + lower_seed * interval.lower_origin_gradient
                                + upper_seed * interval.upper_origin_gradient
                            )
                            direction_gradient = (
                                span * direct_direction
                                + lower_seed * interval.lower_direction_gradient
                                + upper_seed * interval.upper_direction_gradient
                            )
                        factor /= wp.float64(config.samples)
                    if wp.static(active_pose):
                        source_adjoint = factor * (
                            wp.transpose(config.object_to_grid) * origin_gradient
                        )
                        direction_adjoint = factor * (
                            wp.transpose(config.object_to_grid) * direction_gradient
                        )
                        if wp.static(matrix_gradient):
                            for row in range(3):
                                for column in range(3):
                                    gradient[3 * row + column] = (
                                        ray.world_source_relative[row] * source_adjoint[column]
                                        + ray.world_direction[row] * direction_adjoint[column]
                                    )
                            translation_adjoint = -(unpack_rotation(pose) * source_adjoint)
                            for axis in range(3):
                                gradient[9 + axis] = translation_adjoint[axis]
                        else:
                            rotation_adjoint = wp.cross(
                                source_adjoint, ray.object_source
                            ) + wp.cross(direction_adjoint, ray.object_direction)
                            for axis in range(3):
                                gradient[axis] = -source_adjoint[axis]
                                gradient[axis + 3] = rotation_adjoint[axis]
        # The diagnostic writes the same ray derivative before any reduction.
        # Its caller owns O(6 P) storage; the production specialisation retains
        # only reduction partials and performs exactly its existing reduction.
        if wp.static(per_ray):
            if pixel < config.pixels:
                for component in range(components):
                    component_value = gradient[component]
                    partials[wp.int64(pixel) * wp.int64(components) + wp.int64(component)] = (
                        component_value
                    )
                    if not wp.isfinite(component_value):
                        wp.atomic_or(status, 0, 2)
        # All lanes, including padded rays, participate in the fixed reduction.
        elif wp.static(active_pose):
            for component in range(components):
                values = wp.tile(gradient[component])
                total = wp.tile_sum(values)
                wp.tile_store(partials, total, offset=block * components + component)

    return vjp

The flags passed to get_vjp() specialise the kernel when it is prepared. A field-only VJP needs basis weights but no field-value reads or pose reduction. A pose-only VJP needs values and spatial slopes but no scatter into a volume gradient. wp.static removes the unrequested work from each specialisation, so the flags do not branch separately for every ray at runtime.

For pose, each 128-lane block forms binary64 partial sums and subsequent kernels reduce those partials to the requested six or twelve entries. Padded lanes contribute zero and still participate in the block reduction. This avoids a numerical global atomic update for every ray and pose component. The active-volume path does use binary32 scatter atomics, because rays share coefficients. Their addition order is not deterministic, and a post-scatter check detects overflow after individually finite contributions. A fixed volume avoids that storage and contention altogether.

Warp records launches with wp.Tape, and arrays participating in differentiation use requires_grad=True. A supplied output-cotangent array seeds the VJP through Tape.backward(grads=...). The FirstOrderPass boundary checks that the seed has the output’s exact shape and dtype, contiguous storage and the same CUDA device before Warp copies it into the cotangent buffer. The pinned differentiation guide documents these calls and the buffer-overwrite rules. The forward inputs needed by the reverse calculation must remain available. Allocate the output and cotangent buffers before the repeated calls, keep them on the same CUDA device and use a declared stream order.

Table 5.1 identifies the reverse work and shared writes for each choice of active inputs.

Measure both forward and backward traffic and occupancy, including contention where many rays update the same volume coefficients and memory retained between passes.

Clear or initialise gradient buffers according to the call’s accumulation contract. An isolated VJP can overwrite its output, while a node inside a larger reverse graph usually adds to an input’s existing cotangent. Repeated backward calls need explicit seed and reset rules.

FirstOrderPass makes a recorded evaluation’s lifetime explicit when the projector and transmission are used with Warp’s tape. Enter its context to record the forward work, pass its .tape to each participating operator, then call backward() after leaving the recording context. The caller supplies the detector output and its seed, and the workspaces perform their numerical checks before the pass is considered complete.

@dataclass(slots=True)
class FirstOrderPass:
    """Own one forward/reverse evaluation and its explicit diagnostic checkpoint.

    Use as a context manager; pass .tape into canonical forward operators.
    Call backward after exiting the context. Keep original inputs immutable
    until close(). A new iteration requires a fresh context entry after close.
    """

    workspaces: Sequence[CheckedWorkspace]
    tape: Any = field(default=None, init=False, repr=False)
    _state: str = field(default="new", init=False)

    def __enter__(self) -> FirstOrderPass:
        if self._state != "new":
            raise ContractError("a FirstOrderPass owns exactly one evaluation")
        require_no_tape()
        for workspace in self.workspaces:
            workspace.clear_status()
        self.tape = load_warp().Tape()
        self.tape.__enter__()
        self._state = "recording"
        return self

    def __exit__(self, exc_type: Any, exc_value: Any, traceback: Any) -> None:
        self.tape.__exit__(exc_type, exc_value, traceback)
        self._state = "recorded"
        if exc_type is not None:
            self.close()

    def backward(self, output: Any, seed: Any) -> None:
        """Seed a caller-owned output cotangent and complete numerical checks."""
        if self._state != "recorded":
            raise ContractError("backward requires one completed, unconsumed forward pass")
        try:
            wp = load_warp()
            if not isinstance(output, wp.array) or output.grad is None:
                raise ContractError("output must be a Warp CUDA array with a preallocated gradient")
            ctx = prepare_context(device=output.device)
            ctx.array(output, "output", dtype=output.dtype, ndim=output.ndim)
            ctx.array(seed, "seed", dtype=output.dtype, shape=tuple(output.shape), ndim=output.ndim)
            self.tape.backward(grads={output: seed})
            for workspace in self.workspaces:
                workspace.check_status()
            self._state = "complete"
        except BaseException:
            self.close()
            raise

    def close(self) -> None:
        """Clear accumulated cotangents and release every abandoned recording."""
        if self._state == "recording":
            raise ContractError("exit the recording context before closing it")
        if self.tape is not None and self._state != "closed":
            self.tape.zero()
            self.tape.reset()
            for workspace in self.workspaces:
                discard = getattr(workspace, "discard_recording", None)
                if discard is not None:
                    discard(self.tape)
        self._state = "closed"

The original pose and field must remain unchanged through the reverse calculation. close() clears accumulated cotangents and resets the recording before releasing the workspaces. It also handles an abandoned forward pass, whose scratch cannot be reused while a later backward might still depend on it. These checks catch an ownership mistake before it turns into a gradient for a different pose. They do not add higher-order differentiation: the custom projection rule promises first-order products only. The recovery adapter uses explicit VJP calls and prepares its buffers once, which is sufficient for its single forward/reverse chain.

5.6 Independent directional checks

Choose a base pose and a direction in dimensionless scaled coordinates $\mathbf z$. Build perturbed poses using the same left or right composition rule as the analytic derivative. For a vector output, a central difference estimates the JVP:

Equation (5.16) uses a dimensionless step $h$, and the physical perturbations come from the scaling matrix in equation (5.3). For pose parameters, the notation means composing the scaled increment with the base pose, not adding entries to a rotation matrix.

For the implemented VJP, hold a detector seed fixed and take its dot product with each perturbed forward output before forming the central difference. Compare that scalar difference with the dot product of the VJP cotangent, expressed in the same scaled parameter coordinates, and the chosen direction. Both calculations then predict the change in the same seeded scalar, so this check needs neither a stored image Jacobian nor a separately implemented JVP. The seed must correspond to the output being checked: optical depth for the projector alone, or expected counts for the composed rendering chain.

Check a range of steps. In a smooth region, central-difference truncation error decreases quadratically as $h$ shrinks, until subtraction and floating-point errors dominate. One favourable step size can conceal a wrong derivative. A useful discrepancy measure is

The scale $s_0$ in equation (5.17) has the same units as the image derivative. It prevents an almost-zero reference from turning a harmless absolute error into a large relative ratio. Report the absolute error too. The normalisation must not conceal a materially wrong signal change.

A scalar Taylor test avoids dividing two nearly equal output values by a very small step. For a differentiable scalar objective and direction, form

With sufficient local smoothness, equation (5.18) has a quadratic remainder. Leaving the linear term out gives a first-order change when the directional derivative is nonzero. Seeing those two slopes separates a meaningful derivative correction from an objective that is simply insensitive to the chosen direction.

Figure 5.4 compares the error patterns used by the central-difference and Taylor checks.

Analytic fields provide references that share no interpolation code with the implementation. For $\mu(\mathbf{x}^O)=\mu_0+\mathbf{b}^{\mathsf T}\mathbf{x}^O$ along a finite segment entirely inside an affine interpolation region,

The translation derivative in equation (5.19) holds the world ray and object rotation fixed and differentiates with respect to world translation. The vector $\mathbf b$ has units $\mathrm{mm}^{-2}$. Choose $\mu_0$ so attenuation remains nonnegative on the tested neighbourhood. This case checks inverse-pose signs, the physical spacing conversion and distance weights without a moving support boundary.

Test the volume derivative separately using coefficient directions and the fixed-geometry linear operator. Then test pure translations, pure rotations and mixed directions. Mixed directions matter: separate component checks can miss an error in the way shared intermediate contributions are accumulated.

experiments/projection-gradients/run.py carries out the six right-pose directional sweeps with a prescribed signed detector seed and an asymmetric quadratic field. Its gradients.json keeps two comparisons: differences of actual CUDA forward outputs, and differences of an independent binary64 implementation of the same fixed-count quadrature over the uploaded coefficients. The former includes rounding each optical depth to binary32, while the latter exposes the derivative of the unrounded discrete formula. Their small-step error floors should differ. Requiring both to reach the lower floor would confuse the stored interface with the arithmetic used to evaluate its derivative.

The GPU regression suite also includes exact small-slope cases beside a large orthogonal field background. These exercise the corner-difference arrangement in the sampler: an ordinary smooth phantom can pass while a representable transverse derivative disappears through cancellation. The acceptance target remains the independent derivative, not whichever answer looks most plausible after plotting it.

5.7 Boundaries, conditioning and failure cases

Between its exterior support faces, a trilinear field is continuous but not generally differentiable at its knot planes. The implemented field can additionally jump at a nonzero half-cell boundary. In one dimension, let adjacent slopes on either side of a sample centre be $g_-$ and $g_+$. At the knot,

If the two values in equation (5.20) differ, no ordinary derivative exists there. A symmetric finite difference tends to their average, while an executed interpolation branch may report either one-sided slope. Agreement with that average is therefore not the acceptance criterion for a branch derivative. Move the base point into a regular cell for a smooth test, and test the knot with its stated one-sided behaviour.

A physical tangency can also create a large derivative. For a homogeneous sphere of radius $R$ and a line whose perpendicular distance from the centre is $b<R$, assume the finite source-detector segment contains the whole chord. Then

Equation (5.21) diverges as the line approaches tangency from inside. This is a property of the sharp sphere model. A sampled smooth boundary can regularise it at the cost of changing that model. A large gradient near a boundary and a missing gradient caused by a detached clipping operation require different repairs.

Figure 5.5 shows why a small signal can coexist with a steep boundary derivative.

A uniform region may barely change when translated along directions that preserve the traversed material. Symmetry can make several poses indistinguishable. Strong attenuation can suppress count-space sensitivity. The derivative can be correct and small in all these cases. Adding a numerical floor to force a gradient invents a response that the forward model does not have.

Use spatially asymmetric fields, unequal spacings and off-centre rotation pivots to expose coordinate errors. Also include deliberately insensitive cases, where the expected derivative is zero. Tests consisting only of large nonzero gradients will miss accidental offsets and spurious boundary responses.

During debugging, inspect the first layer at which agreement is lost: point motion, interpolated field value, optical depth, detector expectation or scalar objective. A forward finite difference of the final loss combines all of them. Layered analytic cases make the sign or scale error much easier to locate without lowering the final acceptance standard.

5.8 Gradient validation as an executable contract

A derivative result is identified by more than its output array. Record the active parameters, their frame and units, the update side and centre, the field extension, interpolation, sampling policy and precision. Those choices define which mathematical function the derivative claims to approximate.

For any consistent JVP and VJP under Euclidean inner products, parameter direction $\boldsymbol\eta$ and detector seed $\overline{\boldsymbol{\lambda}}$ satisfy

Equation (5.22) checks adjoint consistency without storing the Jacobian. Use several independent, asymmetric directions and seeds. This check catches missing accumulation terms and mismatched interpolation weights, but two mutually consistent implementations of the wrong operator can still pass. The analytic forward and directional checks remain necessary.

If the desired inner products include quadrature or voxel-volume weights, the corresponding adjoint changes. With positive definite weighting matrices $\mathbf M_\theta$ and $\mathbf M_y$,

Equation (5.23) distinguishes that weighted adjoint from the ordinary array transpose computed by a Euclidean VJP. Introduce such weights through the objective and parameter convention deliberately. Multiplying by a voxel volume in one direction only breaks the identity being tested. Record the conventions and evidence in Table 5.2 alongside each gradient check.

Run fresh-buffer and reused-buffer cases, repeat the backward pass under its documented reset rule, and test shared-input accumulation inside a larger computation. For a CUDA operator, include real device execution, memory checks and both ordinary and exceptional numerical regimes in that record.

Choosing a data term, controlling step size and deciding whether a fitted pose explains the acquisition are the next tasks. Chapter 6 uses these derivatives to recover pose.