Chapter 11Rev. 1.0.0

Registering a CT volume to calibrated radiographs

A known CT volume and calibrated X-rays give us a concrete registration problem. We put pose recovery to work across views and assess the recovered geometry independently of the image fit.

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 established the pose chart, image derivatives and optimiser. Here we assemble an application around them: supplied data with documented units, several views that must agree on one pose, and evaluation information withheld from the fit. The executable example is python/dpt/examples/registration.py. It fits calibrated primary photon counts with a fixed volume and fixed acquisition geometry. More general image comparisons and mismatch models have a place in the study, but require the corresponding measurement operators.

Our worked example starts with the supplied CT-ORG case 2 volume and generates two calibrated, orthogonal radiographs. From the identity transform, it recovers the shared pose in 16 accepted updates and reduces RMS target error from 3.77 mm to 0.010 mm. It also reaches the numerical stopping tolerance and predicts two reserved views accurately. The complete recipe includes the prepared input, its provenance and one driver connecting generation, fitting and independent evaluation. The anatomy comes from acquired CT; its 80 keV attenuation assignment and the photon-count radiographs are simulated.

11.1 Define the registration case

The case begins with an attenuation field μO(x)\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 TWO\mathbf{T}_{WO} map object points into world coordinates. For a physical target point xO\mathbf{x}_O and calibrated projection map πv\pi_v of view vv,

xW=TWOxO,uv=πv(xW).\mathbf{x}_W=\mathbf{T}_{WO}\mathbf{x}_O, \qquad \mathbf{u}_v=\pi_v(\mathbf{x}_W).
(11.1)

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 (nz,ny,nx)(n_z,n_y,n_x), with x varying fastest. GridSpec supplies spacing in (x,y,z)(x,y,z) order, the first sample-centre position, and a right-handed orientation matrix. Attenuation is in inverse millimetres. CT numbers are not that attenuation field without a declared conversion, and a polychromatic CT reconstruction does not automatically provide a monochromatic material map. Preserve the conversion recipe and its calibration source alongside the original volume. Chapter 4 describes the interpolation and finite-support assumptions that the projector will actually use.

Each detector specifies shape (H,W)(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.

One representative calibrated view

A fixed acquisition and an object-frame attenuation field A planar section shows the source S at left, the detector D at right, and a volume attached to O. W is the common fixed world frame. Circles with central dots denote axes pointing out of the page. The detector's first axis points upwards, second out of the page and third to the right, forming a right-handed triad.xyxyxzxySO · volumeDW

The field and sampling grid move together with O.

Supplied volumeμO(x), grid geometry\mu_O(\mathbf x),\ \text{grid geometry}

Attenuation in mm⁻¹, with spacing and origin in mm.

Supplied acquisition for each view vsvW, qvpW, N0v\mathbf s_v^W,\ \mathbf q_{vp}^W,\ N_{0v}

Source position, detector pixel centres and open-beam counts, with calibrated detector axes and pitches.

Estimate one shared rigid transformTWO:xOxW\mathbf T_{WO}:\mathbf x_O\longmapsto\mathbf x_W
Project a transformed target into each viewuv=πv(TWOxO)\mathbf u_v=\pi_v(\mathbf T_{WO}\mathbf x_O)
Figure 11.1The registration data and coordinate framesThe attenuation field stays fixed in O and the calibrated apparatus in W. One object-to-world transform must explain every view.

11.2 Compare compatible image quantities

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

λvp(T)=N0vexp[Lvp(T)],Lvp(T)=rvpμO(T1x)ds.\lambda_{vp}(\mathbf{T}) = N_{0v}\exp[-L_{vp}(\mathbf{T})], \qquad L_{vp}(\mathbf{T}) = \int_{r_{vp}}\mu_O(\mathbf{T}^{-1}\mathbf{x})\,ds.
(11.2)

The integration follows the finite source-to-pixel segment rvpr_{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 NvpN_{vp} and a fixed binary inclusion mask mvpm_{vp}, the implemented objective is a sum of Poisson half-deviances:

D(T)=vpmvpd(λvp(T),Nvp),d(λ,N)=λN+Nlog ⁣(Nλ).\begin{aligned} \mathcal{D}(\mathbf{T}) &=\sum_v\sum_p m_{vp}\,d(\lambda_{vp}(\mathbf{T}),N_{vp}),\\ d(\lambda,N) &=\lambda-N+N\log\!\left(\frac{N}{\lambda}\right). \end{aligned}
(11.3)

For N=0N=0, the limiting expression is d(λ,0)=λd(\lambda,0)=\lambda. For positive NN, 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 qkq_k gives

Dqk=v,pmvp(1Nvpλvp)λvpqk.\frac{\partial\mathcal{D}}{\partial q_k} = \sum_{v,p}m_{vp} \left(1-\frac{N_{vp}}{\lambda_{vp}}\right) \frac{\partial\lambda_{vp}}{\partial q_k}.
(11.4)

The FP64 primary composition forms the weighted depth cotangent mvp(Nvpλvp)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 λvp\lambda_{vp}. Assigning each image the same arbitrary scale would erase part of that model.

A binary mask excludes specified pixels from the sum. Fix it before the solve and retain it with the input record. Excluding bad detector elements is defensible when their identity is independently known. Excluding whatever disagrees with the current pose changes the objective and can conceal the anatomy that would disprove the fit. Saturation also needs a declared treatment: the stored ceiling is not an ordinary unsaturated Poisson count. The example accepts finite integer-valued FP32 observations through 2242^{24}. Every excluded pixel must contain zero in the supplied preprocessed observation, so an underflowed prediction there cannot invalidate the count calculation before the mask is applied. Those zeros are excluded placeholders. Preserve the raw source and record the preprocessing recipe with the supplied array’s provenance.

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 logT=L\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 supplied 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.

Loading recorded arrays…

Display and array data

Both views retain all 128 × 128 detector samples. The same logarithmic count window is used for each panel; displayed pixels are not smoothed.

Reaching the stopping criterion

Scaled gradient infinity norm against Accepted update. The vertical scale is logarithmic. Accepted pose and Required ≤ 0.001.Scaled gradient infinity norm0.0010.1101000100,0000481216Accepted update
  • Accepted pose
  • Required ≤ 0.001

Recovering the generating pose

RMS target error (mm) against Accepted update. The vertical scale is logarithmic. Independent evaluation.RMS target error (mm)0.010.110481216Accepted update
  • Independent evaluation

The optimiser sees the two count images, the known volume and fixed calibration. After fixing the result, the evaluator compares its action on the eight target points with their generating positions. The final gradient and target error answer different questions, and both pass.

Recorded result and acceptance criteria
QuantityResultCriterion
Solver terminationGradient tolerance; 16 accepted updatesScaled gradient ≤ 0.001
Final scaled gradient0.000119134≤ 0.001
RMS target error3.769 → 0.00998 mm≤ 1 mm
Worst target error0.01220 mm≤ 2 mm
Reserved-view mean error0.00448 Poisson SD RMS≤ 1 Poisson SD RMS in both views
Reproduction and recorded data
Figure 11.2From two radiographs to a recovered poseRecorded CT-ORG case 2 teaching example: assigned 80 keV attenuation, simulated counts, identity start and independently evaluated recovery.

11.3 Initialise within a defensible search region

An initial pose can come from tracked geometry, independently acquired landmarks, or an explicitly trained proposal method. Record what information produced it. A pose obtained by perturbing the reference transform is a controlled capture-range experiment, and it is not an initialisation procedure available on a new acquisition without that reference.

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

T(q)=T0exp ⁣(Sq^),S=diag(sx,sy,sz,sωx,sωy,sωz).\mathbf{T}(\mathbf{q}) = \mathbf{T}_0\exp\!\left(\widehat{\mathbf{S}\mathbf{q}}\right), \qquad \mathbf{S}=\operatorname{diag}(s_x,s_y,s_z,s_{\omega x},s_{\omega y},s_{\omega z}).
(11.5)

The first three scales have units of millimetres and the last three radians. The hat maps the scaled six-vector into se(3)\mathfrak{se}(3). A rotation scale should reflect the lever arm of the anatomy being fitted: for a point at distance rr from the object origin, a small angle δθ\delta\theta produces displacement of order rδθ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 q\mathbf{q}. PoseChart.gradient() performs that calculation. Simply rescaling six local derivatives is insufficient at a general trial pose. The declared rotation radius also bounds the chart, and crossing it requires a new solve with cleared optimiser history.

The driver starts at q=0\mathbf{q}=\mathbf{0} and performs one local solve. For multistart registration, give each starting pose its own chart and reuse the same fixed observations. Choose the result without consulting the evaluation reference. Charge every proposal its full evaluation cost. Coarse-to-fine registration similarly needs a documented physical image pyramid, with geometry and comparison updated at every level. Blurring counts and retaining the independent Poisson interpretation is generally inconsistent.

DiffPose provides a concrete alternative initialisation strategy: a patient-specific network learns pose from rendered projections of the known CT, then differentiable image optimisation refines its prediction. [23] Its relevance here is the separation between a proposal mechanism and local refinement. Adopting that approach would require its own training record and evaluation of acquisition mismatch, as well as the local solver’s tests.

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 T(q)\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.

The implementation gives each view its own prepared evaluator, all holding the same PoseChart object. Each evaluator returns six derivatives already expressed in that chart. The host sums those small results:

Combining view losses and pose gradientspython/dpt/examples/registration.pyL55–70
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 jvp=qλvp\mathbf{j}_{vp}=\nabla_{\mathbf{q}}\lambda_{vp} at a declared pose. Under the positive-mean independent Poisson model, the expected information is

F(q)=v,pmvpλvpjvpjvpT=vFv(q).\mathbf{F}(\mathbf{q}) = \sum_{v,p}\frac{m_{vp}}{\lambda_{vp}} \mathbf{j}_{vp}\mathbf{j}_{vp}^{\mathsf T} = \sum_v\mathbf{F}_v(\mathbf{q}).
(11.6)

For any direction a\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. The restricted-view example at the end of the chapter then shows what happens when that information is removed.

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 (q,η)(\mathbf{q},\boldsymbol{\eta}), including any justified prior, and assume its nuisance block is positive definite. Eliminating those locally identifiable nuisance directions gives

Hpose,eff=HqqHqηHηη1Hηq.\mathbf{H}_{\mathrm{pose,eff}} = \mathbf{H}_{qq} - \mathbf{H}_{q\eta}\mathbf{H}_{\eta\eta}^{-1}\mathbf{H}_{\eta q}.
(11.7)

The inverse requires a nonsingular nuisance block, possibly after a justified prior has been included. 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.

The supplied-count driver therefore keeps geometry and open-beam counts fixed. It does not estimate scatter, detector offsets or spectral parameters. Where calibrated spectral data and a material representation are available, SpectralPoseEvaluator provides the canonical multiview and shared-nuisance composition discussed in §7.6. Its identified gain/exposure convention still matters: two unrestricted multiplicative scales cannot be recovered separately from their product.

A mismatch study should begin with a separately checked matched-model reference, then change one declared factor at a time. Keep the supplied anatomy and image evidence traceable through every variant. For real acquisitions, residual patterns can motivate a model revision, but cannot by themselves establish which physical mechanism caused the error. A rigid optimiser cannot repair an incorrect anatomical model by becoming more enthusiastic.

11.6 Assemble the complete application

The supplied-input JSON case declares its input arrays by name, with file paths, SHA-256 digests, source descriptions, redistribution rights and units. The application section supplies the grid and chart, a positive ray-sample count, the recovery policy and calibrated views. The accompanying python/dpt/examples/README.md describes the complete input contract and the command for running the example. The worked recipe supplies a prepared, attributed input and writes new observations and execution records to a separate output directory. Its orchestration driver generates those observations before calling the supplied-count fitter; the fitter itself receives only the declared arrays and calibration.

Preparation validates the supplied arrays and uploads the attenuation field once. Each view retains its own observations and binary mask, while sharing that attenuation allocation and a CUDA stream. The source below connects those persistent inputs to the existing projection, transmission and objective operators:

Preparing a calibrated, masked count viewpython/dpt/examples/registration.pyL150–164
                    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 needs one callback and one initial six-vector. Its policy records evaluation and iteration budgets separately, together with stopping tolerances. The returned state is the last accepted state:

Fitting and recording the shared posepython/dpt/examples/registration.pyL232–246
        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,
            },
        )

A gradient-tolerance stop records numerical stationarity in this chart. A small step or loss change records stagnation. Neither establishes geometric success. Preserve line-search failures and exhausted budgets in the study record rather than excluding inconvenient runs after looking at their images.

One final evaluation per view is deliberately outside the solver budget. Reused prediction scratch may contain a rejected trial, so the example recomputes every view at the accepted pose before exporting its count prediction. Those export evaluations are counted separately in fit-report.json. Multiplying solver calls by the number of views gives a reported upper bound on view evaluations: a rejected numerical trial can stop before all views are evaluated. The report also retains each view’s half-deviance and included-pixel count.

Near the solution, rounding intermediate counts can obscure the small loss decrease predicted by the gradient. The worked configuration therefore retains optical depth, predicted counts and their cotangents in FP64. 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 VJP includes the moving endpoint terms at the field’s external support. Support-topology changes and rays lying on interpolation planes still require care: the implemented branch derivative does not establish a unique derivative there. Observation generation uses a separately checked path with 4,096 midpoint samples per ray. Both projection calculations are checked before any Poisson observations are drawn.

The gradient threshold remains 0.001 in the declared dimensionless chart. Early exits for small steps or small relative loss changes are disabled for this worked example, so neither can be mistaken for convergence. This matters near the solution: a substantial noise residual can remain, making the relative loss change tiny while the pose still has a useful descent direction.

A profile of two warmed two-view gradient evaluations on the GB10 records 28 kernel launches, 384 bytes sent to the device and 256 bytes returned, with no CUDA allocation inside the measured range. The persistent buffer addresses also remain unchanged. Those observations support the intended division of work: image calculations stay on the GPU while the host handles the small pose and optimisation state. Setup and final image export are outside that range.

11.7 Evaluate geometric recovery independently

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

et=t^t2,eR=cos1 ⁣(tr(RTR^)12).e_t=\|\widehat{\mathbf{t}}-\mathbf{t}_\star\|_2, \qquad e_R=\cos^{-1}\!\left( \frac{\operatorname{tr}(\mathbf{R}_\star^{\mathsf T}\widehat{\mathbf{R}})-1}{2} \right).
(11.8)

The rotation formula is the geodesic angle in [0,π][0,\pi], and numerical evaluation needs a stable implementation near zero and a half turn. The canonical metric uses an atan2 construction. Translation error measures displacement of the selected object origin, so retain that origin in the report. Moving it changes the interpretation when rotational error is present.

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 xi\mathbf{x}_i:

TRERMS=1Ki=1KT^xiTxi22.\operatorname{TRE}_{\mathrm{RMS}} = \sqrt{\frac{1}{K}\sum_{i=1}^{K} \left\| \widehat{\mathbf{T}}\mathbf{x}_i-\mathbf{T}_\star\mathbf{x}_i \right\|_2^2}.
(11.9)

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.

The example reads optional evaluation landmarks and the reference pose only after fixing the accepted result. They cannot affect its initialisation, line search or stopping rule. Source and uncertainty descriptions are mandatory when a reference is supplied. Without one, the report states that geometric evaluation is unavailable, and it does not replace the missing reference with the fitted transform.

Exporting predictions and evaluating the recovered geometrypython/dpt/examples/registration.pyL250–268
        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),
            },
        )

A capture-range study varies initial error under a predeclared sampling protocol and records both successes and failures. Specify success thresholds before examining the evaluation cases. Use case-level summaries when several perturbations belong to the same acquisition, because hundreds of restarts on one image do not supply hundreds of independent patients. The present eight object-frame probes are the corners of the cube with coordinates ±40 mm. Their locations are fixed before fitting; the known generating transform provides their simulation reference.

11.8 Complete the worked example

Begin with the maintained CUDA environment and the prepared input in public/generated/worked-examples/inputs/registration/. Its known.json records the acquired source, crop, voxel spacing and attenuation assignment; the driver checks the volume hash before using it. Follow the linked recipe with a new output directory. The driver qualifies every required view against an independent CPU integral, then draws the two 128 × 128 count images at an open-beam expectation of 10,000 photons per pixel per view.

The unknown is one six-parameter rigid pose. Calibration and the attenuation field remain fixed, and the initial transform is the identity. The following call passes the generated observations to the canonical fitter without supplying their generating transform:

Running the supplied-count registration from the declared initial poseexperiments/worked-applications/run.pyL393–402
            registration = _fit(
                known_root,
                output,
                public,
                protocol,
                "registration",
                reg_ids,
                RigidTransform(),
                device,
            )

Read fits/registration/optimisation.json next. In the recorded execution, the half-deviance falls from 364,067.02 to 16,359.81. After 16 accepted updates, the largest absolute scaled gradient is 0.000119, below the fixed 0.001 threshold. A nonzero deviance is expected: the observations contain Poisson noise, and one rigid transform cannot reproduce every random pixel fluctuation. Driving that residual to zero would be a rather suspicious achievement.

Only after fixing the result does the evaluator use the generating transform to measure the eight probe errors. RMS error falls from 3.7686 mm to 0.009976 mm, and the worst probe error is 0.01220 mm. The 45° and 135° views were absent from fitting. Their predicted means differ from the generating means by 0.00448 and 0.00412 Poisson standard deviations RMS, respectively. These checks establish successful recovery for this controlled case. They do not turn assigned CT attenuation into patient composition measurements or certify a clinical imaging system.

The same command continues into the view-selection exercise in Chapter 13. Its summary.json passes only if every prescribed solve converges, meets the geometric criteria and predicts the reserved views accurately, and the acquisition comparison passes its declared criterion. The complete record retains source and input hashes, accepted states, rejected trials and independently evaluated outcomes.

Why the earlier registration runs stopped, and what the restricted-view case teaches

The earlier study is retained below because it exposes two different problems. Its full-view fits recovered the controlled geometry, but all 40 fits failed the numerical stationarity test. Replaying a prescribed full-view case isolated rounding of depth and count cotangents in FP32. Retaining FP64 intermediates repaired the inconsistency between the objective and its derivative, but fixed midpoint quadrature still introduced a gradient jump when one sample crossed an interpolation plane. Cell-wise integration removes that sampling artefact. With the original 0.001 gradient threshold and the explicitly revised stopping policy, the same historical case reaches a gradient of 0.0000495. The original stopped record remains unchanged.

Figure 11.3 shows the deliberately restricted case. It combines one view, reduced detector support and fewer photons, so it cannot establish the separate effect of any one restriction. Its failed recovery leaves numerical completion and available image information as separate questions; the stopped fit alone cannot establish that the pose is unrecoverable. Restore the qualified full views to reproduce the successful example above before experimenting with these restrictions.

Restricted, low-count

Probe RMS: 4.9403 → 1.6679 mm · geometric criterion failed

0° view

Observed counts
Observed counts at 0 degrees. Only the central eight detector columns contain observations. The transparent remainder is excluded.

Negative-log display: 0–8

Initial prediction
Initial prediction at 0 degrees. CT-derived simulated radiograph, displayed at the recorded common window.

Negative-log display: 0–8

Accepted prediction
Accepted prediction at 0 degrees. CT-derived simulated radiograph, displayed at the recorded common window.

Negative-log display: 0–8

Hatching: no observations. Rust boundaries: fitted columns 60–67. Predictions outside the strip are excluded from the objective.

128 × 128 detector samples · 8 mm pitch · 1,024 × 1,024 mm detector extent · image top is +v. Display: −log((counts + 0.5) / open beam), 0–8. The pseudocount is absent from the Poisson objective. Termination: line search failed. The accepted pose is not stationary.

Accepted iterates and the unchanged stationarity test

Scaled gradient infinity norm against Accepted iteration. The vertical scale is logarithmic. Orthogonal views, Near-parallel views, Restricted, low-count and Original stationarity threshold.Scaled gradient infinity norm0.0010.1101000100,00001020Accepted iteration
  • Orthogonal views
  • Near-parallel views
  • Restricted, low-count
  • Original stationarity threshold

Geometric error at the recorded accepted poses

125-probe RMS displacement / mm against Accepted iteration (0 is the initial anchor). The vertical scale is logarithmic. Orthogonal views, Near-parallel views, Restricted, low-count and RMS gate (1 mm).125-probe RMS displacement / mm0.010.111001020Accepted iteration (0 is the initial anchor)
  • Orthogonal views
  • Near-parallel views
  • Restricted, low-count
  • RMS gate (1 mm)

Both plots retain the same three examples. Error uses 125 fixed mathematical probes at every recorded pose. Reference poses were disclosed after fitting, without guiding the solver or selecting an iterate.

All 136 application outcomes and numerical records · Sources, rights and display rules. CT anatomy is acquired. Attenuation assignments and radiographs are simulated.

Figure 11.3When a close projection leaves the pose wrongA close match within eight detector columns can coexist with a wrong pose. All eight restricted, low-count fits miss the geometric criterion when coverage, view count and photon population are reduced together.

Figure 11.4 retains all 40 historical outcomes, including the eight restricted cases that fail the geometric criterion. These are the original stopped fits, not additional converged runs of the revised example.

Prescribed starts and final geometric error

Final 125-probe RMS displacement / mm against Initial 125-probe RMS displacement / mm. The vertical scale is logarithmic. Orthogonal, Near-parallel, Restricted, low-count and RMS success threshold.Final 125-probe RMS displacement / mm0.010.1105101520Initial 125-probe RMS displacement / mmOrthogonal: 4.9403, 0.017949Orthogonal: 8.6206, 0.018001Orthogonal: 7.625, 0.017984Orthogonal: 17.834, 0.017984Orthogonal: 4.9403, 0.013281Orthogonal: 8.6206, 0.013283Orthogonal: 7.625, 0.013221Orthogonal: 17.834, 0.013288Orthogonal: 5.0912, 0.0066511Orthogonal: 8.7995, 0.0066779Orthogonal: 7.6751, 0.0066771Orthogonal: 17.871, 0.0066779Orthogonal: 5.0912, 0.010672Orthogonal: 8.7995, 0.010723Orthogonal: 7.6751, 0.01067Orthogonal: 17.871, 0.010686Near-parallel: 4.9403, 0.014455Near-parallel: 8.6206, 0.014451Near-parallel: 7.625, 0.014468Near-parallel: 17.834, 0.01446Near-parallel: 4.9403, 0.017639Near-parallel: 8.6206, 0.017678Near-parallel: 7.625, 0.017751Near-parallel: 17.834, 0.017698Near-parallel: 5.0912, 0.0069554Near-parallel: 8.7995, 0.0069764Near-parallel: 7.6751, 0.0069844Near-parallel: 17.871, 0.0069595Near-parallel: 5.0912, 0.02184Near-parallel: 8.7995, 0.021812Near-parallel: 7.6751, 0.021818Near-parallel: 17.871, 0.021866Restricted, low-count: 4.9403, 1.6679Restricted, low-count: 8.6206, 2.757Restricted, low-count: 7.625, 1.6686Restricted, low-count: 17.834, 1.6703Restricted, low-count: 5.0912, 1.5336Restricted, low-count: 8.7995, 1.5654Restricted, low-count: 7.6751, 1.5341Restricted, low-count: 17.871, 1.5655
  • Orthogonal
  • Near-parallel
  • Restricted, low-count
  • RMS success threshold

Paired full-view registration errors

Near-parallel 125-probe RMS displacement / mm against Orthogonal 125-probe RMS displacement / mm. Case 0 and Case 1.Near-parallel 125-probe RMS displacement / mm00.010.0200.010.02Orthogonal 125-probe RMS displacement / mmCase 0: 0.017949, 0.014455Case 0: 0.018001, 0.014451Case 0: 0.017984, 0.014468Case 0: 0.017984, 0.01446Case 0: 0.013281, 0.017639Case 0: 0.013283, 0.017678Case 0: 0.013221, 0.017751Case 0: 0.013288, 0.017698Case 1: 0.0066511, 0.0069554Case 1: 0.0066779, 0.0069764Case 1: 0.0066771, 0.0069844Case 1: 0.0066779, 0.0069595Case 1: 0.010672, 0.02184Case 1: 0.010723, 0.021812Case 1: 0.01067, 0.021818Case 1: 0.010686, 0.021866
  • Case 0
  • Case 1

All 32 full-view fits passed: probe RMS ≤1 mm, maximum error ≤2 mm, rotation error ≤0.5°. All eight restricted fits failed. The 125 probes are mathematical points, not anatomical landmarks. Overlapping markers are listed separately below.

All 40 fits: pose errors and stopping tests
Case / arrangement / noise / startInitial RMS / mmFinal RMS / mmMaximum / mmRotation / °Gradient norm
0 / Orthogonal views / 0 / 04.94030.01790.02980.01381.866
0 / Near-parallel views / 0 / 04.94030.01450.02460.00990.1013
0 / Orthogonal views / 0 / 18.62060.01800.02990.01380.4486
0 / Near-parallel views / 0 / 18.62060.01450.02460.00990.3637
0 / Orthogonal views / 0 / 27.62500.01800.02980.01380.09378
0 / Near-parallel views / 0 / 27.62500.01450.02460.00990.1409
0 / Orthogonal views / 0 / 317.83410.01800.02980.01380.07825
0 / Near-parallel views / 0 / 317.83410.01450.02460.00990.02584
0 / Orthogonal views / 1 / 04.94030.01330.02060.00530.4154
0 / Near-parallel views / 1 / 04.94030.01760.03090.00941.411
0 / Orthogonal views / 1 / 18.62060.01330.02060.00530.1103
0 / Near-parallel views / 1 / 18.62060.01770.03100.00940.2947
0 / Orthogonal views / 1 / 27.62500.01320.02050.00521.566
0 / Near-parallel views / 1 / 27.62500.01780.03120.00940.3468
0 / Orthogonal views / 1 / 317.83410.01330.02060.00530.5080
0 / Near-parallel views / 1 / 317.83410.01770.03110.00940.4537
0 / Restricted, low-count / 0 / 04.94031.66792.64450.65200.01123
0 / Restricted, low-count / 0 / 18.62062.75704.25321.11910.01632
0 / Restricted, low-count / 0 / 27.62501.66862.64560.65210.01158
0 / Restricted, low-count / 0 / 317.83411.67032.64780.65240.005984
1 / Orthogonal views / 0 / 05.09120.00670.01080.00432.281
1 / Near-parallel views / 0 / 05.09120.00700.01310.00640.8148
1 / Orthogonal views / 0 / 18.79950.00670.01090.00430.3454
1 / Near-parallel views / 0 / 18.79950.00700.01310.00642.310
1 / Orthogonal views / 0 / 27.67510.00670.01090.00430.02352
1 / Near-parallel views / 0 / 27.67510.00700.01310.00641.325
1 / Orthogonal views / 0 / 317.87100.00670.01090.00430.08475
1 / Near-parallel views / 0 / 317.87100.00700.01310.00640.06669
1 / Orthogonal views / 1 / 05.09120.01070.01750.00950.1847
1 / Near-parallel views / 1 / 05.09120.02180.03440.00982.122
1 / Orthogonal views / 1 / 18.79950.01070.01760.00950.2892
1 / Near-parallel views / 1 / 18.79950.02180.03430.00980.2143
1 / Orthogonal views / 1 / 27.67510.01070.01750.00950.06192
1 / Near-parallel views / 1 / 27.67510.02180.03430.00980.09484
1 / Orthogonal views / 1 / 317.87100.01070.01750.00950.05459
1 / Near-parallel views / 1 / 317.87100.02190.03450.00980.4019
1 / Restricted, low-count / 0 / 05.09121.53362.39531.09410.01145
1 / Restricted, low-count / 0 / 18.79951.56542.51511.05590.03269
1 / Restricted, low-count / 0 / 27.67511.53412.39661.09400.01033
1 / Restricted, low-count / 0 / 317.87101.56552.51471.05600.002869

Every solve stopped with line_search_failed. Every final gradient norm exceeds the fixed 0.001 threshold. Reference transforms and two complete reserved views per fit were assessed after the accepted results were fixed.

All 136 application outcomes and numerical records · Sources, rights and display rules. CT anatomy is acquired. Attenuation assignments and radiographs are simulated.

Figure 11.4Every prescribed registration outcomeAcross all 40 fits, the plots compare initial and final pose errors and the two full-view arrangements. All 32 full-view fits meet the geometric gates, while the eight restricted, low-count fits miss them.

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.

References

  1. Otake, Yoshito, Wang, Adam S., Stayman, J. Webster, Uneri, Ali, Kleinszig, Gerhard, Vogt, Sebastian, Khanna, A. Jay, Gokaslan, Ziya L. and Siewerdsen, Jeffrey H. (2013). Robust 3D–2D image registration: application to spine interventions and vertebral labeling in the presence of anatomical deformation. Physics in Medicine and Biology, 58(23), 8535-8553. https://doi.org/10.1088/0031-9155/58/23/8535
  2. Gopalakrishnan, Vivek, Dey, Neel and Golland, Polina (2024). Intraoperative 2D/3D Image Registration via Differentiable X-ray Rendering. Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 11662-11672. https://doi.org/10.1109/CVPR52733.2024.01108