Registering a CT volume to calibrated radiographs

The CT tells us where a structure was when the volume was acquired. The radiograph tells us where its projected attenuation appears now. Registration estimates the rigid transform connecting those two descriptions. Once that transform is known, a point selected in the volume can be located in the imaging system and projected into each calibrated view. A convincing overlay is useful, but it will have to earn its keep as evidence of geometric recovery.

Chapter 6 gave us a pose chart, image derivatives and an optimiser. We now connect them to the volume projector and count model to build a multiview registration algorithm. The design choice is to fit one shared pose while holding attenuation and acquisition geometry fixed. That leaves six unknowns, so we can study initialisation, image weighting and geometric recovery without changing the renderer.

The worked example uses CT-ORG case 2 to define an attenuation field at 80 keV and generates two orthogonal radiographs. The anatomy comes from acquired CT; the attenuation assignment and photon counts are simulated. We will assemble the fitting algorithm around the registration driver, which accepts a supplied volume, calibrated views and count images.

11.1 Define the registration case

The case begins with an attenuation field $\mu_O(\mathbf{x})$, expressed in object coordinates, and a set of radiographs acquired with known source and detector geometry in a shared world frame. Let $\mathbf{T}_{WO}$ map object points into world coordinates. For a physical target point $\mathbf{x}_O$ and calibrated projection map $\pi_v$ of view $v$,

Homogeneous coordinates are implicit in the transform product. This convention fixes what the reported transform means: its translation is expressed in world millimetres, and its rotation acts on object-frame column vectors. A camera-pose estimate from another package may encode the inverse map. Compare the action on physical points before comparing twelve numbers printed by two libraries.

The volume array has shape $(n_z,n_y,n_x)$, with x varying fastest. Use spacing in $(x,y,z)$ order, the first sample-centre position, and a right-handed orientation matrix, as in Chapter 4. Convert CT numbers to attenuation in inverse millimetres before projection; the conversion specifies the energy or material model represented by the field.

Each detector specifies shape $(H,W)$, its first pixel centre, column and row basis vectors, and their physical pitches. Distortion correction, cropping or resampling changes how image indices map to rays. The geometry must describe the array passed to the objective, including any updated origin and spacing. Image orientation guesses are expensive when the optimiser takes them seriously.

A registration case also needs a physical claim about rigidity. A single bone can often be treated as one rigid object for a specified acquisition interval, but neighbouring anatomy may violate that assumption. Removing a neighbouring structure from the volume does not remove its attenuation from the radiograph. Define the region whose agreement will be tested and explain how unmodelled material along those same rays is handled. Figure 11.1 shows the coordinate frames and the shared unknown pose.

11.2 Compare compatible image quantities

The example assumes independent photon counts with a spatially uniform, calibrated open-beam expectation $N_{0v}>0$ in each view. The primary prediction is

The integration follows the finite source-to-pixel segment $r_{vp}$ in world coordinates. Rigid motion preserves length, so the object-frame integrator retains millimetres along that segment. Exposure belongs to the view, while the attenuation volume and the object pose are shared.

For observed counts $N_{vp}$ and a fixed binary inclusion mask $m_{vp}$, the implemented objective is a sum of Poisson half-deviances:

For $N=0$, the limiting expression is $d(\lambda,0)=\lambda$. For positive $N$, a zero predicted mean is invalid. The half-deviance differs from the negative log-likelihood by an observation-only constant, so both give the same pose gradient. With positive predicted means, differentiation in any smooth pose coordinate $q_k$ gives

The FP64 primary composition forms the weighted depth cotangent $m_{vp}(N_{vp}-\lambda_{vp})$ directly; the projection VJP then pulls it back to pose. The example uses reduction="sum", preserving the likelihood contribution of every included observation. Dividing each view by its own pixel count would define a different relative weighting. If exposures differ, their count statistics already enter through $\lambda_{vp}$. Assigning each image the same arbitrary scale would erase part of that model.

Choose the binary mask before fitting, for example to exclude known defective detector elements. Apply it consistently to the objective and gradient. Saturated measurements require a censoring model because the stored ceiling represents an interval of possible counts. The preprocessing required by the runnable example is described in its input guide.

Most energy-integrating or processed fluoroscopy images need a different comparison. Dark subtraction and flat-field division generally destroy the raw-integer-count interpretation. A log-transmission objective must use $\log T=-L$ consistently, while conventional positive log projections use the opposite sign. Where a display pipeline cannot be inverted, correlation or gradient-based similarity can still define a registration criterion, without acquiring a measurement likelihood by wishful notation. Otake and colleagues used normalised gradient information with multistart optimisation to address deformation and content mismatch in spine registration. [18] Figure 11.2 follows the two simulated count images from the initial pose to the recovered one. Its logarithmic window makes the attenuated anatomy visible; the optimiser still compares counts through the half-deviance above. The gradient history and independent target errors show why the fitted images deserve more confidence than an overlay alone.

11.3 Initialise within a defensible search region

An initial pose can come from tracked geometry, independently acquired landmarks or a learned predictor. It needs to lie within the local solver’s capture range. To measure that range in a simulation, perturb the known transform and examine which starting errors the solver can recover.

The example anchors one PoseChart at the initial transform $\mathbf{T}_0$. Its dimensionless coordinates $\mathbf{q}\in\mathbb{R}^6$ define

The first three scales have units of millimetres and the last three radians. The hat maps the scaled six-vector into $\mathfrak{se}(3)$. A rotation scale should reflect the lever arm of the anatomy being fitted: for a point at distance $r$ from the object origin, a small angle $\delta\theta$ produces displacement of order $r\delta\theta$. This relationship gives translation and rotation steps comparable physical meaning without pretending that millimetres and radians are interchangeable.

The anchor and scales remain fixed during a solve. Away from the anchor, the projector’s derivative for a fresh local right increment must be pulled back through the right Jacobian of the exponential before it becomes a derivative in $\mathbf{q}$. PoseChart.gradient() performs that calculation. Simply rescaling six local derivatives is insufficient at a general trial pose. The chosen rotation radius also bounds the chart, and crossing it requires a new solve with cleared optimiser history.

Begin at $\mathbf{q}=\mathbf{0}$ and perform one local solve. A multistart method repeats the solve from several anchors and chooses the best fitted objective. A coarse-to-fine method instead solves on progressively finer image grids, updating the detector geometry and image comparison at each level. Smoothing photon counts introduces correlations, so the coarse images need a compatible objective.

DiffPose tackles initialisation by training a patient-specific network on projections rendered from the known CT, then refining its pose prediction with differentiable image-similarity optimisation. [23] The same separation is useful here: a proposal method gets us near a solution, and the renderer’s derivatives refine it. We use a specified starting pose and calibrated count likelihood to make the local solve explicit. A learned initialiser could later replace that starting rule without replacing the forward model.

11.4 Add views to constrain the same pose

With simultaneous views, or sequential acquisitions during which the object remains still, every view evaluates the same $\mathbf{T}(\mathbf{q})$. Each source and detector can occupy a different position in the common world frame. One view’s almost invisible depth motion may then become a conspicuous lateral motion in another. That benefit depends on the object and field of view as well as angular separation.

Build one evaluator per view and give them the same PoseChart. Each evaluator applies the projection, transmission and objective chain, returning a loss and six chart derivatives. Sum those values to form the callback used by the optimiser:

class SharedPoseObjective:
    """Sum views in one immutable chart; each evaluator owns its image scratch."""

    def __init__(self, views: tuple[PreparedView, ...], chart: PoseChart) -> None:
        if not views or any(view.evaluator.chart is not chart for view in views):
            raise ContractError("all views must share the same PoseChart object")
        self.views, self.chart = views, chart

    def __call__(self, parameters: Vector) -> Evaluation:
        values = tuple(view.evaluator(parameters) for view in self.views)
        return Evaluation(
            math.fsum(value.loss for value in values),
            tuple(math.fsum(value.gradient[k] for value in values) for k in range(6)),
        )

Averaging transforms fitted independently to each view does not minimise this joint objective. Gradients must also use the same frame and scales before they can be added. Requiring one chart object keeps the anchor and scales identical in every view.

A local information calculation explains the geometric benefit. Write $\mathbf{j}_{vp}=\nabla_{\mathbf{q}}\lambda_{vp}$ at the pose being examined. Under the positive-mean independent Poisson model, the expected information is

For any direction $\mathbf{a}$, its quadratic form is a sum of nonnegative squared directional sensitivities. Adding an independent view cannot reduce this information matrix in positive-semidefinite order under the stated fixed model. It can still leave a common null direction unresolved. Nearly repeated views may improve precision where information already exists while contributing little to the weakest direction.

Local curvature says nothing about a distant pose producing another plausible projection. A repeated anatomical pattern can admit competing solutions even where the neighbourhood of each solution is well conditioned. Likewise, if the object moves between exposures, one shared pose no longer describes the acquisition. Fit time-specific poses or account for motion before interpreting a sharper local minimum as stronger evidence. The worked case uses two full orthogonal views so that the reader can first complete this shared-pose recovery.

11.5 Account for imperfect agreement

A primary renderer fitted to measured data will often leave structured residuals. Miscalibrated geometry can shift edges across the image. Scatter can add a smooth component. Anatomical change can create a disagreement confined to one region but impossible for a single rigid transform to repair. Inspecting residuals by view and location helps distinguish these explanations, while the total objective hides that structure.

Allowing a nuisance parameter $\boldsymbol{\eta}$ changes what the images can determine about pose. Take a positive-semidefinite joint information matrix for $(\mathbf{q},\boldsymbol{\eta})$, including any justified prior, and assume its nuisance block is positive definite. Eliminating those locally identifiable nuisance directions gives

Directions that can be explained by either pose or nuisance variation lose pose information. For example, unconstrained source motion can imitate changes in the object’s position relative to that source. A numerical optimiser cannot distinguish parameter names that produce the same measurements.

For this calculation, hold geometry and open-beam counts fixed. If calibration must be estimated as well, add only identifiable nuisance parameters using the joint composition from §7.6. For example, an exposure scale and a detector gain appear only through their product unless one is independently constrained.

To diagnose a structured residual, vary one plausible cause at a time: geometry, scatter or the attenuation field. Each change asks a different physical question. A rigid optimiser cannot repair an incorrect anatomical model by becoming more enthusiastic.

11.6 Assemble the complete application

Assemble the application in three steps. First prepare the fixed attenuation field and, for each view, the calibrated geometry, observed counts and mask. Then define a callback that renders every view at a trial pose and sums its loss and gradient. Finally pass that callback and the initial pose coordinates to the optimiser from Chapter 6. The listings below implement these steps with the operators we have already built.

Upload the attenuation field once and reuse it across views. Keep each view’s observations, mask and intermediate images on the GPU. Sharing a CUDA stream makes the forward and reverse operations execute in order while those arrays remain alive:

                    problem = PrimaryPoseProblem(
                        grid=grid,
                        geometry=geometry,
                        attenuation=attenuation,
                        observation=wp.array(observed.reshape(-1), dtype=wp.float32, device=device),
                        objective=ObjectiveSpec(
                            kind="poisson", domain="counts", reduction="sum", weighted=True
                        ),
                        open_beam=beam,
                        weights=wp.array(mask.reshape(-1), dtype=wp.float32, device=device),
                        samples_per_ray=samples,
                        precision=configuration.get("precision", "float64"),
                        integration=configuration.get("integration", "midpoint"),
                    )
                    evaluator = PrimaryPoseEvaluator(problem, chart, device=device, stream=stream)

Image-sized intermediate arrays remain on the device throughout optimisation. The host exchanges rigid-transform values, a scalar objective and six derivatives with the device, together with numerical-status checks. Summing those derivatives on the host does not require copying the detector images. The example evaluates views sequentially on its shared stream.

The optimiser receives the callback and an initial six-vector. A line search evaluates trial poses and keeps the last accepted state; iteration and evaluation limits bound the work:

        initial: Vector = (0.0,) * 6
        result = recover_parameters(objective, initial, policy=policy)
        recovered = chart.pose(result.parameters)
        run.write_json(
            "optimisation.json",
            {
                "pose_object_to_world": asdict(recovered),
                "chart": asdict(chart),
                "policy": asdict(policy),
                "precision": views[0].evaluator.problem.precision,
                "integration": views[0].evaluator.problem.integration,
                "result": asdict(result),
                "stationary": result.stationary,
            },
        )

Stop when the scaled gradient is small enough for the desired numerical accuracy. A small step or loss change can also signal stagnation, so retain the gradient when diagnosing a stopped solve. Geometric accuracy is assessed separately in §11.7.

After the solve, render each view once more at the accepted pose. A line search may leave a rejected trial in reused image buffers; refreshing them ensures that the exported predictions describe the returned transform.

Near the solution, rounding intermediate counts can obscure the small loss decrease predicted by the gradient. Retaining optical depth, predicted counts and their cotangents in FP64 preserves those small differences. Fixed midpoint samples cause a separate problem: when one crosses an interpolation plane, the sampled derivative can jump. With integration="cell_gauss", the ray is split at those planes and two Gauss nodes integrate each cubic segment exactly up to floating-point arithmetic. This preserves the trilinear field introduced in Chapter 4 while removing that sampling artefact.

The pose derivative includes the moving endpoints where a ray enters and leaves the field’s external support. Rays on interpolation planes or at changes in intersection topology need separate attention because the local branch derivative may not be unique.

Near a solution, the count noise can leave a substantial residual even when the remaining pose error is small. Relative loss changes may then become tiny before the gradient does. This is why we use the scaled gradient, rather than the appearance of a flat loss curve, to decide whether to continue.

11.7 Evaluate geometric recovery independently

Suppose an independently established reference transform is $\mathbf{T}_\star=(\mathbf{R}_\star,\mathbf{t}_\star)$ and the estimate is $\widehat{\mathbf{T}}=(\widehat{\mathbf{R}},\widehat{\mathbf{t}})$. Separate translation and rotation errors:

The rotation formula is the geodesic angle in $[0,\pi]$. Use a stable numerical evaluation near zero and a half turn, such as the atan2 construction in the accompanying implementation. Translation error refers to the selected object origin; rotational error makes a different origin move by a different amount.

Translation and rotation errors describe the transform, but the practical question is how far its errors displace a point in the anatomy. Rotation makes that displacement depend on the point’s location, so we combine both errors at independently selected object-frame target locations $\mathbf{x}_i$:

This is a three-dimensional target registration error in millimetres. A projected target error in detector pixels answers a different question: a large depth error can produce little displacement in one projection. Report the RMS target error and the worst target error, with the target distribution and reference uncertainty. A reference produced by registering the same radiograph with a related objective may be useful for comparison, but its errors are not independent merely because its output file has a different name.

Evaluate the returned transform against reference landmarks and pose only after fitting. Keeping this information out of initialisation and stopping lets target error assess the image-based calculation independently.

        for view in views:
            evaluator = view.evaluator
            # A rejected line-search trial may be the last occupant of scratch.
            final_value = evaluator(result.parameters)
            final_losses[view.identifier] = final_value.loss
            prediction = evaluator.prediction.numpy().reshape(evaluator.problem.geometry.shape)
            write_array(run, f"prediction-{view.identifier}.npy", prediction)
        run.write_json(
            "fit-report.json",
            {
                "final_half_deviance_by_view": final_losses,
                "included_pixels_by_view": {v.identifier: v.included_pixels for v in views},
                "solver_evaluations": result.evaluations,
                # A rejected trial can fail before some view evaluators are called.
                "solver_view_evaluation_upper_bound": result.evaluations * len(views),
                "additional_final_export_evaluations": len(views),
                "evaluation": evaluate_reference(case, configuration, recovered),
            },
        )

To examine capture range, repeat the fit from different initial errors and plot target error against the starting displacement. Choose target locations across the structure rather than near the rotation centre alone. In the worked case, eight probes occupy the corners of a cube at object-frame coordinates ±40 mm.

11.8 Complete the worked example

The worked case uses two 128 × 128 simulated count images with an open-beam expectation of 10,000 photons per pixel in each view. The generating pose is kept for evaluation; fitting receives the images and calibration.

Start the six-parameter pose at the identity and keep calibration and attenuation fixed. Use the shared-view objective and optimiser assembled above to fit the supplied observations.

The fit recovers the pose in 16 accepted updates, reducing RMS target error from about 3.77 mm to 0.01 mm. The scaled gradient meets the stopping criterion, and predictions agree closely with the simulated means in two views withheld from fitting. These checks test numerical completion, geometry and prediction separately.

The deviance remains nonzero because the observations contain Poisson noise: one rigid transform cannot reproduce every random pixel fluctuation. Driving that residual to zero would be a rather suspicious achievement.

Try changing the angular separation while holding exposure and initial pose fixed. The resulting errors show how a second view constrains motion that was weakly observed in the first. Chapter 13 turns that choice into a calculation.

The known volume has constrained every pose trial in this chapter. Chapter 12 removes that advantage: geometry is fixed, and the unknowns become the attenuation values themselves. The renderer stays recognisable. The space of plausible explanations becomes much larger.