Choosing the next view and where to spend the photons

A registration can fit its first X-ray closely while leaving the object free to move along a poorly observed direction. Another image might resolve that uncertainty. Taking almost the same view again mostly provides more photons along directions already measured. Moving the C-arm incurs its own costs, and the angle with the prettiest derivative image may be physically inaccessible. We need to predict which feasible measurement will reduce uncertainty in the task we care about.

We will build a next-view selector around the renderer from Chapter 11. At each candidate geometry it predicts how photon counts respond to pose, converts those derivatives into information, and scores the resulting target uncertainty. The worked example then generates the chosen CT-derived simulated radiograph and refits the pose.

The forward model and pose convention come from Chapter 11. Chapter 12 supplies the harder extension in which the unknown is a volume. We will keep track of the information available before each decision, because a policy that has already seen tomorrow’s X-ray tends to perform suspiciously well.

13.1 Define the task and the information available

Suppose the current data are $\mathcal{Y}_0$ and the next acquisition is described by $a$: its source and detector geometry together with calibrated exposure settings. A decision rule can use the known volume, $\mathcal{Y}_0$, the current pose estimate and a declared uncertainty model. It may render hypothetical measurements at candidate acquisitions. It cannot use their future measured pixel values.

For registration, the quantity of interest need not be an equally weighted combination of translation and rotation. A small rotation produces a larger displacement at a distant surgical target than near the rotation centre. We therefore supply object-frame target points and score their predicted world-coordinate variance. These points describe the task using anatomy already available to the policy. Withheld reference transformations are used later to measure error.

A fixed design chooses the entire sequence before seeing any new measurements. A sequential design chooses the next acquisition, receives its observation, updates the estimate and chooses again. The executable example implements one decision in the latter sequence. Its current precision must describe uncertainty at that decision, including previous measurements exactly once. Adding their information again would reward the policy for counting its homework twice.

Figure 13.1 separates the current estimate and candidate predictions from future observations. The selector evaluates expected information at the current pose; the next measured image enters only after a view has been chosen.

13.2 Choose a task-specific design objective

Keep the pose chart fixed while comparing candidates. At the current transform $\widehat{\mathbf G}_{WO}$, use dimensionless coordinates $\mathbf z$ with positive diagonal scale matrix $\mathbf S$:

The first three entries of $\mathbf S$ set translation scales in millimetres, and the remaining entries set rotation scales in radians. All derivatives in this section are evaluated at $\mathbf z=\mathbf 0$. Equation (13.1) is a local right update. A precision matrix expressed in another frame, at another chart centre or with different scales must be transformed before use.

For candidate $a$, assume conditionally independent ideal counts $N_{a,p}\sim\operatorname{Poisson}(\lambda_{a,p}(\mathbf z))$. At positive means, the score contribution is $(N_{a,p}/\lambda_{a,p}-1)\nabla_{\mathbf z}\lambda_{a,p}$. Its expectation is zero and its variance uses $\operatorname{Var}(N_{a,p})=\lambda_{a,p}$. Independence removes cross-pixel score covariances, giving the expected Fisher information

A large image derivative contributes little useful information if its noise is large or if it repeats an already well-observed parameter direction. For the primary model, let $L_{a,p}$ be the optical depth and $n_{0,a,p}$ the fixed open-beam expectation. Since $\lambda_{a,p}=n_{0,a,p}\exp(-L_{a,p})$, equation (13.2) becomes

where $\mathbf J_{L,a}$ has one row $\partial L_{a,p}/\partial\mathbf z$ per pixel. This form avoids dividing by a tiny predicted count. A zero incident beam contributes zero information. For correlated read noise or energy-integrating measurements, derive the weighting from the corresponding likelihood.

Let the current pose distribution be approximated locally by a Gaussian with positive-definite precision $\boldsymbol\Lambda_0$ in this same chart. Adding expected information gives

Equation (13.4) is a local Gaussian approximation built from expected likelihood curvature. It is useful when uncertainty is concentrated around one pose. Include earlier measurements once in $\boldsymbol\Lambda_0$; an optimiser’s damping matrix has a different role and should not be treated as measured information.

For an object-frame target $\mathbf x_q$, differentiating its world position under equation (13.1) gives

Here $[\mathbf x_q]_\times\mathbf v=\mathbf x_q\times\mathbf v$. The translation part of $\mathbf H_q$ moves every target alike, while its rotation part retains the distance and direction from the object origin. Averaging target-position variance gives the scalar criterion

Its units are square millimetres. The criterion measures variance, so bias needs a separate assessment. Listing 13.1 evaluates it with Cholesky solves rather than an explicit inverse. Positive-definite current precision makes those solves well defined.

def target_jacobians(np: Any, pose: RigidTransform, targets: Any, scales: Any) -> Any:
    """Return d(world target)/dz in mm for T_WO exp((S z)^) at z=0."""
    rotation = np.asarray(pose.rotation, dtype=np.float64).reshape(3, 3)
    result = np.empty((len(targets), 3, 6), dtype=np.float64)
    for index, (x, y, z) in enumerate(targets):
        skew = np.asarray([[0.0, -z, y], [z, 0.0, -x], [-y, x, 0.0]])
        result[index, :, :3] = rotation
        result[index, :, 3:] = -rotation @ skew
    result *= scales[None, None, :]
    if not np.isfinite(result).all():
        raise NumericalError("target derivatives exceed FP64; check target coordinates and scales")
    return result


def target_variance(
    np: Any, precision: Any, target_jacobian: Any, maximum_condition: float
) -> tuple[float, Any, float]:
    """Mean target-position variance in mm², without forming an inverse to score."""
    if not np.isfinite(precision).all():
        raise NumericalError("pose precision is non-finite; discard this candidate score")
    values = np.linalg.eigvalsh(precision)
    if values[0] <= 0 or values[-1] / values[0] > maximum_condition:
        raise NumericalError(
            "pose precision must be positive definite and within maximum_precision_condition; "
            "check the stated uncertainty and parameter scales, without adding numerical jitter"
        )
    factor = np.linalg.cholesky(precision)
    right = target_jacobian.reshape(-1, 6).T
    whitened = np.linalg.solve(factor, right)
    if not np.isfinite(whitened).all():
        raise NumericalError(
            "target covariance solve exceeds FP64; check the precision and chart scales"
        )
    # hypot scales its sum of squares; individual products must not underflow
    # before they contribute to an otherwise representable total variance.
    normalisation = math.sqrt(len(target_jacobian))
    root_variance = math.hypot(*(float(value) / normalisation for value in whitened.flat))
    variance = root_variance * root_variance
    if not math.isfinite(variance) or variance <= 0:
        raise NumericalError(
            "positive target-position variance is outside FP64 range; check the stated "
            "uncertainty and physical units, or use a range-preserving score representation"
        )
    inverse_factor = np.linalg.solve(factor, np.eye(6))
    covariance = inverse_factor.T @ inverse_factor
    if not np.isfinite(covariance).all() or bool((np.diag(covariance) <= 0).any()):
        raise NumericalError(
            "local covariance is outside FP64 range; check the precision and chart scales "
            "before exporting this candidate"
        )
    return variance, covariance, float(values[-1] / values[0])

Nuisance parameters can substantially change the ranking. If $\boldsymbol\eta$ describes uncertain calibration and the joint precision is partitioned into pose and nuisance blocks, the marginal pose precision is the Schur complement $\boldsymbol\Lambda_{zz}-\boldsymbol\Lambda_{z\eta}\boldsymbol\Lambda_{\eta\eta}^{-1}\boldsymbol\Lambda_{\eta z}$, assuming the required positive definiteness. Holding nuisance parameters fixed instead uses the pose block directly and can understate uncertainty. Our example makes that restriction explicit: acquisition calibration and attenuation are fixed, and only pose is uncertain.

13.3 Separate candidate selection from continuous optimisation

For a finite eligible set $\mathcal A$, the decision is simply $a^\star=\operatorname{argmin}_{a\in\mathcal A}\Phi(a)$. Score every candidate against the same current precision, pose and targets. Do not update the current precision after scoring a hypothetical candidate, because that would make the answer depend on enumeration order. Exact score ties in the example prefer lower cost, then the candidate name.

Enumerating candidate views uses the pose derivatives we have already built. Continuous optimisation of source or detector position would also require derivatives of the mean and pose Jacobian with respect to acquisition geometry. Choosing a finite set therefore gives us a complete selector using first-order pose sensitivities; extending it to continuous motion would require mixed pose–acquisition derivatives.

A criterion that includes a reconstruction solve adds another dependency. Differentiating a fixed number of solver iterations gives a derivative of that algorithm. Implicit differentiation instead uses an equation characterising the solution and requires appropriate local uniqueness, regularity and convergence. These quantities need not agree for an unfinished solve or at an active-set change. Holding the current estimate fixed does not differentiate its dependence on the acquisition. Record that approximation explicitly.

13.4 Restrict choices to feasible acquisitions

The candidate set should come from an imaging-system specification. Source and detector positions must respect working distance, angular travel and detector coverage. A reachable endpoint can still require a motion path that collides with equipment. Exposure settings may be discrete, and moving to a candidate may consume more time than taking its image. Geometry and movement costs therefore belong to each candidate record.

Apply mechanical constraints before scoring information: retain only reachable source–detector configurations with acceptable motion paths. Figure 13.2 illustrates an endpoint constraint by marking excluded configurations. The selector then compares the remaining candidates under the acquisition budget.

Calibrate the open-beam expectation for each candidate’s source setting and geometry. Moving the source changes distance and coverage, so exposure alone may not determine the photon population at each pixel. For a budget in expected detector photons, sum the open-beam array over the detector.

Drag to orbit · Scroll to zoom

• Allowed detector endpoint
• Source-clearance exclusion
• Outside candidate domain

At +60°: source clearance 353.72 mm, detector configuration allowed.

$\mathbf s=R_s\mathbf u,\quad\mathbf d=-R_d\mathbf u$relative to O, with $\mathbf u=\cos\alpha\,\mathbf e_0+\sin\alpha\,\mathbf e_1$.

$\|\mathbf s-\mathbf c\|\geq R_s/2,\quad\mathbf c=R_s\mathbf e_0$Allowed: $\arccos(7/8)\leq|\alpha|\leq\pi/2$.

Model-world distances: $R_s=353.72\,\mathrm{mm}$,$R_d=348.28\,\mathrm{mm}$. $\mathbf e_0=(0,-1,0),\ \mathbf e_1=(1,0,0)$. The ±90° candidate domain and point-clearance disk are illustrative. They do not specify device travel or a collision-free path.

Use the angle slider to inspect the endpoint constraint. Camera orbiting changes only the view.

13.5 Choose a second view for registration

Begin with the attenuation field, current pose, local precision and object-frame target points. For each candidate, call the pose-sensitivity operation from Chapter 5 at its calibrated geometry. The same volume allocation can be reused throughout the search.

The pose derivatives use physical millimetres and radians. Multiply their columns by the chart scales to form $\mathbf J_{L,a}$, then weight each row by the square root of the predicted count. Forming $\log\lambda=\log n_0-L$ before exponentiation preserves the optical-depth calculation’s range. A GPU matrix product reduces the $P\times6$ weighted Jacobian to the $6\times6$ information matrix.

Listing 13.2 implements this calculation. Retaining the Jacobian and a weighted copy costs $96P$ bytes in double precision. Process one candidate at a time to reuse that storage; the candidate comparison itself needs only the small information matrices.

def candidate_information(
    *,
    np: Any,
    torch: Any,
    ctx: Any,
    grid: GridSpec,
    geometry: DetectorGeometry,
    attenuation: Any,
    pose: Any,
    beam_host: Any,
    scales: Any,
    samples_per_ray: int,
    precision: Literal["float32", "float64"] = "float32",
    integration: Literal["midpoint", "cell_gauss"] = "midpoint",
) -> tuple[Any, dict[str, Any]]:
    """Use canonical depth derivatives and an FP64 CUDA Fisher reduction.

    Preparation allocates one candidate's O(6P) Jacobian and weighted copy.
    Torch views alias Warp storage on the same CUDA stream. There is no
    per-pixel Python loop and no image/Jacobian download. This is an offline
    candidate evaluation, not a captured or profiled optimisation hot path.
    """
    wp, pixels = ctx.wp, geometry.pixels
    try:
        projection = prepare_projection(
            grid,
            geometry,
            ProjectionSpec(samples_per_ray, precision=precision, integration=integration),
            device=str(ctx.device),
            stream=ctx.stream,
        )
        transmission = (
            prepare_transmission(
                TransmissionSpec(beam="none"),
                max_pixels=pixels,
                device=str(ctx.device),
                stream=ctx.stream,
            )
            if precision == "float32"
            else None
        )
        with ctx.scope():
            beam = wp.array(beam_host.reshape(-1), dtype=wp.float32, device=ctx.device)
            depth = wp.empty(pixels, dtype=projection.dtype, device=ctx.device)
            log_transmission = wp.empty(
                pixels if precision == "float32" else 0, dtype=wp.float32, device=ctx.device
            )
            ones = wp.ones(pixels, dtype=projection.dtype, device=ctx.device)
            jacobian = wp.empty(6 * pixels, dtype=wp.float64, device=ctx.device)
        project_optical_depth(
            attenuation,
            pose,
            workspace=projection,
            out_L=depth,
            stream=ctx.stream,
        )
        if transmission is not None:
            transmit(
                depth,
                out_log_T=log_transmission,
                workspace=transmission,
                stream=ctx.stream,
            )
        projection_pose_sensitivities(
            attenuation,
            pose,
            ones,
            out_jacobian=jacobian,
            workspace=projection,
            stream=ctx.stream,
        )
        beam_tensor = wp.to_torch(beam)
        log_mean = beam_tensor.to(dtype=torch.float64).log_()
        if precision == "float64":
            log_mean.sub_(wp.to_torch(depth))
        else:
            log_mean.add_(wp.to_torch(log_transmission))
        illuminated = beam_tensor > 0
        # A zero supplied beam carries no information. Positive means outside
        # the normal FP64 range are rejected; no tail, floor or pixel is hidden.
        log_tiny = math.log(float(np.finfo(np.float64).tiny))
        if bool((illuminated & (log_mean < log_tiny)).any().item()):
            raise NumericalError(
                "a positive candidate count mean is below normal FP64 range; "
                "use a range-preserving design implementation before ranking this case"
            )
        root_mean = log_mean.mul_(0.5).exp_()
        depth_jacobian = wp.to_torch(jacobian).reshape(pixels, 6)
        weighted = depth_jacobian.clone()
        weighted.mul_(scales).mul_(root_mean[:, None])
        if not bool(torch.isfinite(weighted).all().item()):
            raise NumericalError("weighted pose derivatives exceed FP64; discard this candidate")
        if bool(((depth_jacobian != 0) & illuminated[:, None] & (weighted == 0)).any().item()):
            raise NumericalError(
                "a nonzero weighted derivative underflowed; discard this candidate"
            )
        magnitude = weighted.abs()
        minimum = math.sqrt(float(np.finfo(np.float64).tiny))
        maximum = math.sqrt(float(np.finfo(np.float64).max) / (6.0 * pixels))
        if bool((((magnitude > 0) & (magnitude < minimum)) | (magnitude > maximum)).any().item()):
            raise NumericalError(
                "candidate Fisher products exceed the declared FP64 accumulation range; "
                "rescale the pose chart or use a range-preserving reduction"
            )
        # Weighted depth derivatives give J_L^T diag(lambda) J_L directly.
        # Avoid a division by rounded, potentially zero count predictions.
        information = weighted.T @ weighted
        result = information.cpu().numpy().copy()
        if not np.isfinite(result).all():
            raise NumericalError("candidate information is non-finite; discard this score")
        diagnostics = {
            "pixels": pixels,
            "precision": precision,
            "integration": integration,
            "zero_beam_pixels": int((~illuminated).sum().item()),
            "omitted_positive_beam_pixels": 0,
            "jacobian_and_weighted_copy_bytes": 2 * 6 * pixels * 8,
            "fisher_reduction": "FP64 CUDA; only the 6 by 6 information matrix is downloaded",
        }
        return 0.5 * (result + result.T), diagnostics
    finally:
        # Warp owners and the host beam must survive every pending shared-stream operation.
        wp.synchronize_stream(ctx.stream)

Use adequate precision for the exponential weights and matrix reduction. Strongly attenuated rays contribute little information, while poor parameter scaling can make some weighted derivatives much larger than others. Check their range before interpreting the smallest information directions.

Listing 13.3 adds each candidate’s information to the same current precision, computes the target criterion and records threshold-defined local Fisher rank. Rank depends on the stated tolerance and parameter scaling. A view can reduce variance along an existing direction without increasing rank, so rank alone is not the selection rule.

                for candidate in eligible:
                    information, diagnostics = candidate_information(
                        np=np,
                        torch=torch,
                        ctx=ctx,
                        grid=grid,
                        geometry=candidate["geometry"],
                        attenuation=attenuation,
                        pose=pose_device,
                        beam_host=candidate["beam"],
                        scales=scales_device,
                        samples_per_ray=samples,
                        precision=cfg.get("precision", "float32"),
                        integration=cfg.get("integration", "midpoint"),
                    )
                    variance, covariance, condition = target_variance(
                        np, prior + information, target_derivatives, maximum_condition
                    )
                    eigenvalues = np.linalg.eigvalsh(information)
                    largest = float(eigenvalues[-1])
                    if eigenvalues[0] < -rank_tolerance * max(largest, np.finfo(np.float64).tiny):
                        raise NumericalError(
                            "candidate Fisher matrix is indefinite beyond rank tolerance"
                        )
                    rank = int(np.count_nonzero(eigenvalues > rank_tolerance * largest))
                    rows.append(
                        {
                            "name": candidate["name"],
                            "cost": candidate["cost"],
                            "mean_target_variance_mm2": variance,
                            "local_fisher_rank": rank,
                            "posterior_precision_condition": condition,
                            "diagnostics": diagnostics,
                        }
                    )
                    matrices.append(information)
                    covariances.append(covariance)
                # Exact score ties prefer the smaller stated cost, then the name.
                selected = min(
                    rows,
                    key=lambda row: (row["mean_target_variance_mm2"], row["cost"], row["name"]),
                )

Start with a 0° image at 1,000 expected open-beam photons per pixel and fit its pose. Compare 5°, 30°, 60° and 90° candidates at the same exposure. Use the base view’s Fisher matrix as $\boldsymbol\Lambda_0$ and score variance at the eight object-frame targets with coordinates ±40 mm.

In the first replicate, the selector chooses 60° even though 90° is farther from the base view. The criterion depends on the anatomy and target locations, not angular separation alone. Generate the chosen view’s count image, refit the shared pose, and compare the predicted uncertainty with the realised target errors in Figure 13.3.

The 5° baseline spends the second exposure on nearly the same direction. It illustrates the benefit of constraining a weakly observed motion. A fixed orthogonal view would be a more demanding baseline. If the first image admits several distant plausible poses, replace the local Gaussian approximation with a pose ensemble or a multimodal uncertainty model.

13.6 Design views for reconstruction

An attenuation volume replaces the six pose coordinates with perhaps millions of unknown coefficients. The design question still needs a task: variance of a regional average, visibility of a specified feature or reconstruction error under a declared prior ensemble. Summing uncertainty equally over all voxels can spend acquisitions on parts of the field that do not matter to that task.

A local approximation has a precision operator of the form

where $\mathbf H_0$ describes current local precision and $\mathbf W_a$ follows the candidate’s measurement model. A dense volume covariance would require quadratic storage. Instead, an implementation could apply precision to vectors and solve for the covariance action needed by a small set of task functionals. For a linear task $q(\boldsymbol\mu)=\mathbf c^{\mathsf T}\boldsymbol\mu$, the local variance is $\mathbf c^{\mathsf T}\mathbf H_a^{-1}\mathbf c$, and a linear solve can obtain it without constructing the inverse.

To implement equation (13.7) without a dense matrix, combine Jacobian-vector and transpose products with the regulariser Hessian, then solve for the task’s covariance action. Reuse the volume reverse pass for the transpose operation and add a matching forward directional derivative for the other half. At active nonnegativity constraints, define the local approximation on feasible directions.

Burger and colleagues use sequential Bayesian design to choose projection angle and lateral position in two-dimensional parallel-beam tomography. Their Gaussian prior and additive Gaussian noise permit posterior-covariance criteria such as A- and D-optimality. [38] We use the same sequential principle—update uncertainty, score candidates, acquire one view—but derive local information from Poisson pose sensitivities and score variance at selected targets.

Compare reconstruction policies at the same measurement budget and include the time spent deciding which view to take. An expensive design search can be appropriate for a costly acquisition. It can also become a remarkably elaborate method of avoiding the next reconstruction. Spectral-channel choices add the calibrated response and material-identifiability conditions developed in §12.6. They cannot be ranked from geometric coverage alone.

13.7 Allocate exposure and state the budget correctly

Under a fixed spectrum and linear independent-count model, suppose a relative exposure $e$ scales the expected counts as $\lambda_{a,p}(e)=e\lambda_{a,p}(1)$. The depth Jacobian stays fixed, so

Information grows linearly in this model, but the resulting target variance generally does not, because it depends on the inverse of the sum of prior and measurement information. Detector saturation, spectral changes or exposure-dependent processing invalidate the simple scaling. Each requires a forward and noise model that follows the actual control setting.

The example permits several candidate records with the same geometry and different calibrated beam/cost pairs. This evaluates discrete exposure alternatives for one next view. A multi-view allocation would instead choose exposures jointly under a total budget and per-view bounds, and would need to retain the dependence of its criterion on every allocation. Movement time can be included only when the cost is defined consistently with the current system state.

Expected detector photons, tube-current–time product in mAs and acquisition time in seconds are different budgets. None is absorbed dose. A dose-constrained extension must define whose absorbed energy per unit mass is being constrained and use a separately validated deposition model that accounts for spectrum, geometry and material. The number of Monte Carlo histories used to estimate a prediction is a computational setting, not the exposure given to a subject.

The worked comparison holds expected open-beam photons fixed and changes the second angle. It therefore demonstrates an angular decision at equal exposure. An exposure-allocation curve would need several controlled photon populations and must retain that axis label, even if “dose reduction” would fit rather nicely in the abstract.

13.8 Complete and check the acquisition sequence

The complete sequence is: fit the first image, rank feasible views, acquire the selected image, then refit using both. The acquisition-design driver supplies the ranking step for specified candidates; compose it with the registration solver for the two fits. In the simulated comparison, the selected policy and the 5° baseline share the same first image. Generate a future observation only after choosing its geometry; this preserves the information available at the decision.

Repeat the sequence with eight independent noise pairs. Keep the anatomy, exposure, targets and candidate angles fixed so the differences reflect the noisy acquisitions. Assess the final poses at the same targets and use the reserved views to inspect prediction.

In this example, selection reduces mean squared target error by about 71% relative to the near-parallel baseline at equal exposure. That result summarises the eight noise pairs on this anatomy.

In two of the eight pairs, the selected view produces a larger target error than the 5° baseline. The Fisher criterion predicts variance averaged over possible future count images, while each refit receives one new noisy image alongside the shared base image. A view with lower predicted variance can therefore give the worse result in an individual pair. The mean improvement includes both reversals, without requiring every noisy observation to be unusually obliging.

Try adding a fixed orthogonal baseline, restricting the available angles, or moving the target points. Each change asks a concrete question: how much does selection add, how does accessibility limit the choice, and which anatomy should the acquisition prioritise?

13.9 A closing word

Volume I ends with a calculation that can help choose its next measurement. We began with photon survival through a slab. Reaching this decision has required a working renderer and derivatives that can guide an inverse calculation. The complete examples now recover unknown parameters from simulated observations and check the returned solutions against independent references. It is quite a lot of machinery for choosing another X-ray, but we can now account for what the machinery does.

The registration and reconstruction examples reuse the forward model while changing the unknowns from six pose coordinates to a material field. A new application starts with the same question: which parts of that measurement model still describe the data?

Volume II is planned to cover neural representations and learned methods, and Volume III to develop volume reconstruction further. Those subjects have a place to begin in the forward model we have built here. Thank you for working through the physics and code with me. I hope that, when you next encounter a difficult imaging problem, a working calculation feels within reach.