Recovering an object’s pose from its projection

A projected bone edge is three pixels to the left of where the renderer puts it. Moving the object might repair that disagreement. So might moving the source, changing the detector calibration, or deciding that the bright line was an instrument rather than bone. The optimiser receives none of this context. It receives numbers and permission to change six of them.

Pose recovery gives those six numbers a physical meaning. We need an image model, a rule for comparing its prediction with the observation, and a way of deciding whether the resulting movement recovered the object. These decisions are coupled: an objective that ignores intensity scale can tolerate an exposure error, but it also discards any geometric information carried only by that scale.

Even with anatomy and acquisition calibration fixed, some motions barely change the projection. Those weak directions determine how precisely the image can locate the object.

Build the pose solver around one operation: given a candidate rigid transform, return an image loss and its six pose derivatives. Connect the projector, transmission model and image objective in that order, then pass the objective’s image derivatives backwards through them. PrimaryPoseEvaluator collects that calculation, and recover_pose() chooses the trial transforms. We will first fix acquisition calibration and fit counts or log transmission. Chapter 7 will let calibration change too.

6.1 Pose recovery as an inverse problem

Let the attenuation field be specified in object coordinates, and let the source and detector geometry be known in world coordinates. The unknown rigid transform is $\mathbf{G}_{WO}\in SE(3)$, mapping object coordinates into world coordinates. At each candidate pose, transform the physical source-to-detector segments into the object frame, integrate attenuation, and predict the quantity represented by the observed image.

Write $\mathbf{f}(\mathbf{G}_{WO})$ for that prediction and $\mathbf{y}$ for the observation. The components of $\mathbf{f}$ might be expected photon counts, calibrated detector signal, or a specified processed image. They must describe the same measurement domain as $\mathbf{y}$. A count prediction cannot be compared directly with a display image whose greyscale has been inverted and windowed.

A pose estimate minimises a data discrepancy, optionally with prior information:

In equation (6.1), the regulariser $\mathcal{R}$ expresses information beyond the image: a tracked initial pose with stated uncertainty, for example, or a mechanically restricted motion. Setting it to zero leaves a data-only fit. A numerical damping term used to control an optimisation step is a different object. It need not represent any belief about the final pose.

The known quantities deserve as much attention as the six unknowns. An error in voxel spacing changes the physical object. An error in principal point moves every predicted ray. An attenuation field from another anatomical state changes what those rays cross. Holding these quantities fixed does not make their errors disappear: it gives pose the job of compensating for them.

The rigid model is appropriate when the relative positions of the represented material structures remain fixed over the acquisition. It can describe the pose of a single bone even when neighbouring structures move independently, provided the image comparison isolates the relevant anatomy. A rigid fit to a whole region containing several moving bones asks one transform to satisfy incompatible geometry.

6.2 What should the images agree on?

For independent ideal photon counts $N_p$ with strictly positive expected counts $\lambda_p(\mathbf{G})$, the Poisson negative log-likelihood, after dropping terms independent of pose, is

The fixed pixel set $\mathcal{M}$ in equation (6.2) identifies the observations included in the fit. Zero observed counts are allowed: the product $0\log\lambda_p$ contributes zero. The expected count must remain positive. Detector energy integration, electronic noise and processing require a measurement model appropriate to those operations.

The count likelihood developed in Chapter 2 supplies the statistical comparison. For a local pose coordinate $\xi_k$, differentiating this objective gives

Each pixel in equation (6.3) contributes its count discrepancy through a pose sensitivity. A bright residual on an edge can have a large geometric effect because the expected count changes rapidly with movement there. A discrepancy in a locally uniform region may have little pose sensitivity, however conspicuous it looks.

For additive Gaussian measurement errors with a known, pose-independent positive-definite covariance $\mathbf{C}$, use a whitened residual:

Whitening in equation (6.4) measures discrepancies relative to their uncertainty and correlation. Here $\mathbf{C}^{-1/2}$ denotes a whitening factor whose transpose times itself is $\mathbf{C}^{-1}$. With independent errors it divides each residual by its standard deviation. If covariance depends on pose, the full Gaussian likelihood also contains a log-determinant term and derivatives of covariance. Freezing the weights defines an approximation to that likelihood.

Sometimes the available image has uncertain scale and offset. Normalised cross-correlation compares centred image patterns. With $\widetilde{\mathbf{f}}=\mathbf{f}-\bar f\mathbf{1}$ and $\widetilde{\mathbf{y}}=\mathbf{y}-\bar y\mathbf{1}$ on the same fixed mask, define

For nonconstant images, equation (6.5) is unchanged by a positive scalar gain and a uniform offset applied to either image. That is useful when those quantities are unknown. It also makes them unavailable as evidence for pose. Constant images have zero centred norm and no defined correlation. Adding a denominator floor changes the objective and its derivative near that case.

A gradient-image objective can suppress slowly varying backgrounds while emphasising anatomical boundaries. Its benefit depends on which boundaries remain visible in both images. Tools, collimator edges and saturation boundaries can supply stronger gradients than the anatomy. Restricting the fit to a justified region may help, but the region should not quietly change as the candidate pose becomes less convenient. Table 6.1 compares what each objective assumes and which image information it discards.

If a motion changes only overall positive intensity scale, NCC cannot detect it. The invariances of an objective therefore limit which movements its gradient can recover.

Start with squared error in the chosen measurement domain or the Poisson count objective. ObjectiveSpec selects between them and carries fixed nonnegative pixel weights. NCC and gradient comparisons need their own image derivatives before they can enter the same loop. Diagonal weights represent diagonal precision, so use a separate whitening operation for correlated noise. For mean reduction, decide which denominator defines the objective: reduction="mean" divides by all entries, including those with zero weight.

For Poisson counts, the code subtracts the observation-only value of the negative log-likelihood at perfect agreement. The resulting half-deviance is nonnegative in exact arithmetic and has the same parameter derivative as equation (6.2). This avoids carrying a large irrelevant constant through a line search. The observed == 0 branch gives the mean as its loss and one as its count-space derivative, without evaluating a logarithm of zero. It includes a zero predicted mean through the continuous extension of this zero-observation loss, while a positive observation still requires a positive mean.

@cache
def evaluation_kernel(
    poisson: bool,
    weighted: bool,
    gradient: bool,
    masked: bool = False,
    precision: str = "float32",
):
    value_dtype = wp.float64 if precision == "float64" else wp.float32

    @wp.kernel(module="unique", module_options=OPTIONS)
    def evaluate(
        prediction: wp.array(dtype=value_dtype),
        observation: wp.array(dtype=wp.float32),
        weights: wp.array(dtype=wp.float32),
        valid: wp.array(dtype=wp.uint8),
        size: int,
        normalisation: wp.float64,
        seed: wp.array(dtype=value_dtype),
        partials: wp.array(dtype=wp.float64),
        status: wp.array(dtype=wp.int32),
    ):
        block, lane = wp.tid()
        p = block * TILE + lane
        loss = wp.float64(0.0)
        if p < size:
            active = bool(True)  # noqa: UP018 - Warp runtime Boolean.
            if wp.static(masked):
                if valid[p] > wp.uint8(1):
                    wp.atomic_or(status, 0, 1)
                active = valid[p] == wp.uint8(1)
            if active:
                predicted = wp.float64(prediction[p])
                observed = wp.float64(observation[p])
                weight = normalisation
                if wp.static(weighted):
                    weight *= wp.float64(weights[p])
                difference = predicted - observed
                derivative = difference
                loss = wp.float64(0.5) * difference * difference
                if wp.static(poisson):
                    if observed == wp.float64(0.0):
                        loss = predicted
                        derivative = wp.float64(1.0)
                    else:
                        derivative = difference / predicted
                        loss = poisson_half_deviance(predicted, observed)
                if wp.static(precision == "float64" and not poisson):
                    # A small fixed weight can keep the final squared loss
                    # representable even when the unweighted square overflows.
                    loss = (wp.float64(0.5) * weight * difference) * difference
                else:
                    loss *= weight
                if not wp.isfinite(loss):
                    wp.atomic_or(status, 0, 2)
                if wp.static(gradient):
                    weighted_derivative = weight * derivative
                    if wp.static(precision == "float64" and poisson):
                        if not wp.isfinite(derivative):
                            # FP32 observations/weights bound this numerator;
                            # do not overflow a removable unweighted quotient.
                            weighted_derivative = (weight * difference) / predicted
                    result = value_dtype(weighted_derivative)
                    seed[p] = result
                    if not wp.isfinite(result):
                        wp.atomic_or(status, 0, 2)
            elif wp.static(gradient):
                seed[p] = value_dtype(0.0)
        total = wp.tile_sum(wp.tile(loss))
        wp.tile_store(partials, total, offset=block)

    return evaluate

Near agreement, the Poisson residual is much smaller than either count. The kernel expresses the loss through a relative residual and log1p, and when even that subtraction would cancel, it uses the corresponding short series. Away from agreement it subtracts the two logarithms instead of forming a ratio that could lose range. Binary64 operands matter before the division, not only after it. The detector seed is rounded to binary32 for the existing transmission interface, while loss terms enter a fixed binary64 reduction tree. Loss and seed therefore share the chosen objective without requiring an image-sized loss buffer.

evaluate_objective() exposes these outputs through caller-owned out_loss and optional out_seed buffers. A Poisson observation must be a finite nonnegative integer-valued count in the stored representation. A zero predicted mean with a positive observed count is rejected: adding a floor would silently replace the likelihood. Mean reduction on an empty image is also rejected. The composed count-space path retains the range limits identified in Chapter 5, including an unrepresentable detector seed.

6.3 Parameter scales and pose updates

Use a local rigid increment $\boldsymbol{\xi}=[\boldsymbol{\rho};\boldsymbol{\phi}]$, with translation coordinates $\boldsymbol{\rho}$ in millimetres and rotation coordinates $\boldsymbol{\phi}$ in radians. For a left update,

The increment in equation (6.6) is expressed in world coordinates. A right update instead expresses its increment in the current object frame. Gradients, Jacobians and updates must use the same choice. For a finite increment, the exponential produces a translation $\mathbf{J}(\boldsymbol{\phi})\boldsymbol{\rho}$, with $\mathbf{J}$ the $SO(3)$ left Jacobian derived in Chapter 3. This translation is not generally equal to the coordinate $\boldsymbol{\rho}$ placed in the wedge matrix.

A raw Euclidean step in these six coordinates treats one millimetre and one radian as equally large. At a point 100 mm from the rotation origin, a rotation of 0.01 rad produces about 1 mm of transverse movement. That relationship supplies a useful initial scale.

Choose a reference displacement $\ell>0$ in millimetres and a characteristic radius $r_*>0$ in millimetres. A dimensionless local variable $\mathbf{u}$ can be mapped to the rigid increment by

With equation (6.7), unit translation and rotation coordinates correspond to movements of roughly the same size near radius $r_*$. The radius is measured from the origin about which the chosen increment rotates. A remote world origin can make rotation appear to contain a large translation of the anatomy. Choosing a useful origin and converting the pose consistently is often more effective than asking a diagonal scale matrix to repair that coupling.

For the whitened residual, let $\mathbf{J}_{\xi}$ be its Jacobian with respect to a local increment at the current pose. The chain rule gives

The rescaling in equation (6.8) changes both the derivative coordinates and the interpretation of a step norm. It does not add information to the image. A nearly invisible motion remains nearly invisible after rescaling, although a sensible scale makes its numerical severity easier to assess.

PoseChart selects the right-composed alternative and gives each solve one fixed anchor. Its default scales are one millimetre for translation and 0.01 radian for rotation, corresponding to the 100 mm lever-arm example. They are explicit choices the caller can replace. The chart rejects trial rotation vectors at or beyond its declared radius, which cannot exceed $\pi$, rather than wrapping them into coordinates that would invalidate the optimiser’s history.

Changing a fixed chart coordinate changes the vector inside the exponential from the anchor. A fresh right increment instead multiplies the current pose by another exponential. Because rigid increments do not generally commute, these operations need not have the same local derivative away from the anchor. PoseChart.gradient() therefore first pulls the projector’s local right cotangent back through right_jacobian_se3(), then applies the physical scales. The Jacobian includes translation–rotation coupling. Merely multiplying the six local cotangents by the scales would be correct at the anchor and wrong at a general trial pose, an awkward place for a registration routine to become less accurate.

Linearising the residual as $\mathbf{r}+\mathbf{J}_u\Delta\mathbf{u}$ yields a damped Gauss–Newton step for the data-only objective:

The damping in equation (6.9) limits movement when the linear approximation is poor or a direction is weakly constrained. Solve the linear system rather than forming a matrix inverse. If the objective includes a differentiable prior, its gradient and curvature approximation also belong in the step.

Evaluate the actual objective at the proposed rigid pose before accepting it. A large predicted reduction can fail because the ray intersections or visible structures changed substantially. Rejecting the proposal, increasing damping and trying a smaller step gives the current local model a chance to remain relevant.

6.4 Initialisation and capture range

A derivative describes nearby changes. It has no obligation to point towards the correct bone from the far side of the detector.

Consider a one-dimensional image feature $f(x)=\exp[-x^2/(2\sigma^2)]$, with position $x$ and width $\sigma$ in pixels. Compare it with a copy shifted by $\delta$ using normalised overlap over the real line:

The overlap in equation (6.10) follows by completing the square in the product of the two Gaussians. It isolates a common registration problem: corresponding features cease to interact when their separation greatly exceeds their width. This is an uncentred overlap on an infinite domain, distinct from the finite-mask NCC objective above.

For the discrepancy $1-c_{\mathrm{ov}}$, the derivative is

At large separation, equation (6.11) becomes very small even though the shift error is large (Figure 6.1). A stopping rule based only on gradient magnitude can therefore report convergence while the features barely overlap.

Smoothing broadens features and can extend the range over which their overlap changes appreciably. It also removes distinctions between nearby structures. Repeated vertebrae, ribs or approximately symmetric outlines can still produce several plausible alignments. A coarse image may make an individual minimum wider while making the choice between minima less informative.

Initialisation supplies information on that larger scale: calibrated positioning, tracking, a previous frame, a coarse search, or a learned predictor. A multi-start method spends more computation on competing hypotheses, while a learned initialiser spends training effort on predicting a useful hypothesis. Both are followed by an image fit whose local behaviour still needs examination.

DiffPose trains a patient-specific pose predictor on X-rays rendered from the preoperative CT, then refines its prediction with differentiable rendering and a sparse multiscale correlation objective. Otake and colleagues use normalised gradient information with a multi-start covariance-matrix adaptation search. Both devote substantial machinery to finding a useful starting region. Our loss-and-gradient evaluator supplies the local refinement calculation. Pair it with an initialisation strategy appropriate to the available pose information, and test that combination over the starting errors it will encounter. [23, 18]

Capture range is consequently a property of a specified anatomy, view, objective, initialisation distribution and optimisation procedure. Define the starting translations, rotation axes and angles, stopping rule and success criterion before summarising it. A cube of translation perturbations and a ball of rotations are different sampling distributions from independent perturbations of Euler angles. They should not acquire the same label simply because all were described as “random”.

6.5 Single-view ambiguity and identifiability

One view supplies many pixels, so a single projection does not automatically imply that six rigid parameters are unrecoverable. What matters is how independently those pixels respond to the six motions.

Perspective already reveals one source of coupling. In a source-centred coordinate system with detector plane at $Z=f_d>0$, a point $(X,Y,Z)$ with $Z>0$ projects to $(u,v)=(f_dX/Z,f_dY/Z)$. All coordinates here are in millimetres, and the geometric Jacobian is

The matrix in equation (6.12) annihilates the direction $(X,Y,Z)^{\mathsf T}$: moving a point along its source ray leaves its projected location unchanged (Figure 6.2). Several points on a known rigid object cannot generally move independently along their own rays. Their common distances and orientation can recover information that one point lacks. Unknown object scale, however, introduces a further coupling with depth.

For the fixed-covariance Gaussian model, the local data information in scaled coordinates is

The quadratic form in equation (6.13) measures the squared whitened image change generated by a small motion in direction $\mathbf{v}$. A zero value identifies a first-order invisible direction. A small value identifies a direction whose effect is weak relative to the specified noise.

The singular value decomposition makes these directions explicit:

In equation (6.14), the right singular vectors describe combinations of translations and rotations, and the squared singular values give their local information. Full column rank supports first-order local identifiability. It neither excludes a second distant solution nor proves that a local optimiser will reach this one.

A perfectly spherical, radially symmetric attenuation field is unchanged by rotations about its centre. No transmission image of that field can identify those rotations. An almost symmetric object may replace the exact ambiguity with a very weak direction. The optimiser then has enough information to select a number, but perhaps too little to select it reproducibly.

For independent views with known relative acquisition geometry and a shared pose parameterisation, their information matrices add:

A second view helps in equation (6.15) when it sees a motion that the first view constrains poorly. Repeating nearly the same view improves precision in already visible directions but cannot resolve an exact shared null direction. The relative geometry must be known or estimated with its uncertainty: allowing every view its own unconstrained camera pose can reintroduce the freedom we intended to remove.

Priors also add curvature. They can stabilise an estimate along a weak direction, but the resulting precision belongs partly to the prior. Report that distinction when interpreting an uncertainty estimate. An inverse information matrix is a local Gaussian approximation near a suitably regular solution, but it does not represent distant alternative alignments.

For an offline view of the weak directions, assemble a small diagnostic Jacobian from selected observations and call sensitivity_spectrum() with the parameter scales. Retain its singular values, directions and rank threshold. Keep this matrix out of the repeated pose evaluation, which only needs the loss gradient.

def sensitivity_spectrum(
    jacobian_rows: Sequence[Sequence[float]],
    parameter_scales: Sequence[float],
    *,
    precision_weights: Sequence[float] | None = None,
    relative_threshold: float = 1e-8,
    absolute_threshold: float = 0.0,
) -> SensitivitySpectrum:
    """SVD of sqrt(W) J S; right directions use dimensionless chart coordinates.

    Singular values and an explicitly chosen threshold describe local sensitivity,
    not global identifiability or clinical recovery. Duplicate views add no new
    independent directions, although weighting can change threshold-defined rank.
    Noise correlations require an explicitly prewhitened Jacobian; diagonal
    precision weights must not stand in for an unknown covariance model.
    """
    rows, columns = len(jacobian_rows), len(parameter_scales)
    if not 1 <= rows <= 4096 or not 1 <= columns <= 32:
        raise ContractError("offline diagnostic is bounded to 4096 observations and 32 parameters")
    if any(len(row) != columns for row in jacobian_rows):
        raise ContractError("every Jacobian row must contain each active parameter")
    units = tuple(finite_scalar(value, "parameter scale", minimum=0) for value in parameter_scales)
    if min(units) == 0:
        raise ContractError("parameter scales must be positive")
    relative = finite_scalar(relative_threshold, "relative threshold", minimum=0)
    absolute = finite_scalar(absolute_threshold, "absolute threshold", minimum=0)
    if relative >= 1:
        raise ContractError("relative rank threshold must be less than one")
    weights = tuple(precision_weights) if precision_weights is not None else (1.0,) * rows
    if len(weights) != rows:
        raise ContractError("one fixed precision weight is required per observation")
    scaled = [
        [_scaled_entry(value, units[column], weights[index]) for column, value in enumerate(row)]
        for index, row in enumerate(jacobian_rows)
    ]
    if not all(math.isfinite(value) for row in scaled for value in row):
        raise NumericalError("scaled sensitivity exceeds binary64 range")
    # Zero rows preserve the missing directions when observations < parameters,
    # without requesting the large square left-singular-vector matrix.
    scaled.extend([[0.0] * columns for _ in range(max(0, columns - rows))])
    np: Any = importlib.import_module("numpy")
    _, values, right = np.linalg.svd(np.asarray(scaled, dtype=np.float64), full_matrices=False)
    singular = tuple(float(value) for value in values)
    if not all(map(math.isfinite, singular)):
        raise NumericalError("sensitivity decomposition produced nonfinite singular values")
    threshold = max(absolute, relative * singular[0])
    rank = sum(value > threshold for value in singular)
    condition = singular[0] / singular[-1] if rank == columns else None
    if condition is not None and not math.isfinite(condition):
        condition = None
    return SensitivitySpectrum(
        singular,
        tuple(tuple(float(value) for value in row) for row in right),
        rank,
        threshold,
        condition,
        rows,
        columns,
    )

The optional precision weights scale rows by their square roots. Correlated noise requires a Jacobian already transformed by the appropriate whitening operator: diagonal weights cannot manufacture that operator. Parameter scales act on columns, so the returned right singular vectors describe dimensionless combinations of the declared inputs. Rank depends on the chosen threshold and noise scale. The production solver below requests VJPs and does not retain this diagnostic matrix.

6.6 Assemble the optimisation loop

The evaluator supplies an objective and six derivatives for each pose. L-BFGS can use that interface directly: it estimates curvature from accepted steps and changes in their gradients. Build its history in the fixed, scaled chart from section 6.3. Compose each candidate from the original anchor and pull its derivative back into the same coordinates. If you choose a new anchor, restart the curvature history. The damped Gauss–Newton step above is an alternative when you also compute the required residual-Jacobian products.

The predicted image, residual and derivative must use the same mask, interpolation and acquisition model. A renderer that changes sampling density between the objective and its derivative can supply a direction for a different function. A line search that changes the mask between candidates can reward dropping difficult pixels. Table 6.2 follows the calculations and retained state through one iteration.

Image-sized arrays should remain on the device through rendering, comparison and differentiation. Copying a full image to the CPU merely to reduce it to one scalar introduces avoidable transfers. The host may need that scalar for step acceptance, but it rarely needs all its contributing pixels. For six pose parameters, a small dense curvature system is inexpensive, but storing a full image Jacobian or an autodiff tape can still dominate memory. The Jacobian products derived in Chapter 5 can support a solve without retaining every derivative entry.

Allocate the volume, observation and image workspaces before entering the optimiser. The evaluation below then projects optical depth, forms the detector expectation, evaluates the loss and pulls its derivative back to pose. For Poisson counts, composing these steps in depth space gives the weighted cotangent $N-\lambda$. The default binary64 path retains that composition through the objective, avoiding an intervening binary32 count cotangent. The explicit binary32 path uses the separate transmission and objective interfaces. In either case, only the scalar loss and six chart derivatives need to reach the optimiser.

    def __call__(self, parameters: Vector) -> Evaluation:
        require_no_tape()
        try:
            return self._evaluate(parameters)
        finally:
            # Even a failed launch can follow an enqueued asynchronous H2D copy.
            # Release its read of pinned staging before a retry may overwrite it.
            self.context.wp.synchronize_stream(self.context.stream)

    def _evaluate(self, parameters: Vector) -> Evaluation:
        pose = self.chart.pose(parameters)
        ctx, wp, problem = self.context, self.context.wp, self.problem
        validate = not self._inputs_checked
        self.projection.clear_status()
        if self.transmission is not None:
            self.transmission.clear_status()
        self.objective.clear_status()
        with ctx.scope():
            self._pose_view[:] = pose.packed()
            wp.copy(self.pose_device, self._pose_host, stream=ctx.stream)
            project_optical_depth(
                problem.attenuation,
                self.pose_device,
                out_L=self.optical_depth,
                workspace=self.projection,
                stream=ctx.stream,
                validate=validate,
            )
            if problem.precision == "float64":
                evaluate_primary_objective(
                    self.optical_depth,
                    problem.observation,
                    open_beam=problem.open_beam,
                    weights=problem.weights,
                    out_prediction=self.prediction,
                    out_loss=self.loss,
                    out_depth_seed=self.depth_seed,
                    workspace=self.objective,
                    stream=ctx.stream,
                    validate=validate,
                )
            else:
                assert self.transmission is not None
                counts = problem.objective.domain == "counts"
                transmit(
                    self.optical_depth,
                    problem.open_beam if counts else None,
                    out_counts=self.prediction if counts else None,
                    out_log_T=None if counts else self.prediction,
                    workspace=self.transmission,
                    stream=ctx.stream,
                    validate=validate,
                )
                evaluate_objective(
                    self.prediction,
                    problem.observation,
                    weights=problem.weights,
                    out_loss=self.loss,
                    out_seed=self.image_seed,
                    workspace=self.objective,
                    stream=ctx.stream,
                    validate=validate,
                )
                transmission_vjp(
                    self.optical_depth,
                    problem.open_beam if counts else None,
                    seed_counts=self.image_seed if counts else None,
                    seed_log_T=None if counts else self.image_seed,
                    out_grad_L=self.depth_seed,
                    workspace=self.transmission,
                    stream=ctx.stream,
                    validate=validate,
                )
            projection_vjp(
                problem.attenuation,
                self.pose_device,
                adj_L=self.depth_seed,
                out_pose=self.pose_seed,
                workspace=self.projection,
                stream=ctx.stream,
                validate=validate,
            )
            self.projection.check_status()
            if self.transmission is not None:
                self.transmission.check_status()
            self.objective.check_status()
            self._inputs_checked = True
            wp.copy(self._loss_host, self.loss, stream=ctx.stream)
            wp.copy(self._gradient_host, self.pose_seed, stream=ctx.stream)
            wp.synchronize_stream(ctx.stream)
        gradient = tuple(float(value) for value in self._gradient_view)
        return Evaluation(float(self._loss_view[0]), self.chart.gradient(parameters, gradient))

Validate fixed device inputs when preparing the evaluator and keep them immutable during the solve. Each trial then supplies a pose and checks the resulting numerical status. Reuse the pinned host staging buffer only after its transfer has finished. The finally synchronisation in the listing enforces that lifetime even when an evaluation fails.

The L-BFGS two-loop recursion applies an inverse-curvature approximation using a bounded history of parameter steps and gradient differences. This is the limited-memory construction of Liu and Nocedal. In our renderer, the useful saving is also in the interface: each evaluation returns a gradient without retaining a full image Jacobian. The six-coordinate history and dot products are small host operations. [34]

def _direction(gradient: Vector, pairs: list[tuple[Vector, Vector, float]]) -> Vector:
    """Two-loop recursion; curvature vectors all belong to the same scaled chart."""
    if not pairs:
        return _steepest_direction(gradient)
    q = gradient
    coefficients: list[float] = []
    for displacement, change, reciprocal in reversed(pairs):
        coefficient = reciprocal * _dot(displacement, q)
        coefficients.append(coefficient)
        q = tuple(a - coefficient * b for a, b in zip(q, change, strict=True))
    scale = 1.0
    if pairs:
        displacement, change, _ = pairs[-1]
        denominator = _dot(change, change)
        if not math.isfinite(denominator) or denominator <= 0.0:
            return _steepest_direction(gradient)
        scale = _dot(displacement, change) / denominator
        if not math.isfinite(scale) or scale <= 0.0:
            return _steepest_direction(gradient)
    result = tuple(scale * value for value in q)
    for (displacement, change, reciprocal), coefficient in zip(
        pairs, reversed(coefficients), strict=True
    ):
        beta = reciprocal * _dot(change, result)
        result = tuple(
            a + (coefficient - beta) * b for a, b in zip(result, displacement, strict=True)
        )
    return tuple(-value for value in result)

Curvature pairs are useful only when the step and gradient difference share coordinates and provide usable positive curvature. The driver checks their inner product and reciprocal before retaining them. If the proposed direction is nonfinite or fails the descent test, it discards the history and tries steepest descent.

def recover_parameters(
    evaluator: Evaluator,
    initial: Vector,
    *,
    policy: RecoveryPolicy | None = None,
    observe: Callable[[IterationRecord], None] | None = None,
    cancelled: Callable[[], bool] | None = None,
) -> RecoveryResult:
    """Minimise a deterministic objective in an immutable, dimensionless chart.

    The evaluator must use the same observation, forward model and parameter
    coordinates for every call. It may reuse device scratch, but must complete
    numerical-status checks before returning. Rejected trials never replace the
    accepted parameters in the result. Observer exceptions are not swallowed.
    """
    selected = policy or RecoveryPolicy()
    x: Vector = tuple(finite_scalar(value, "initial parameter") for value in initial)
    if not x:
        raise ContractError("recovery requires at least one active parameter")
    calls = 0

    def evaluate(point: Vector) -> Evaluation:
        nonlocal calls
        if calls >= selected.max_evaluations:
            raise _BudgetExhaustedError
        calls += 1
        return _checked(evaluator(point), len(x))

    current = evaluate(x)
    history: list[IterationRecord] = []
    pairs: list[tuple[Vector, Vector, float]] = []
    reason: StopReason = "iteration_budget"
    history.append(IterationRecord(0, current.loss, _norm(current.gradient), 0.0, calls, 0))
    for iteration in range(1, selected.max_iterations + 1):
        if cancelled is not None and cancelled():
            reason = "cancelled"
            break
        if _norm(current.gradient) <= selected.gradient_tolerance:
            reason = "gradient_tolerance"
            break
        direction = _direction(current.gradient, pairs)
        slope = _dot(direction, current.gradient)
        if not all(map(math.isfinite, direction)) or not math.isfinite(slope) or slope >= 0:
            pairs.clear()
            direction = _steepest_direction(current.gradient)
        try:
            accepted = _line_search(evaluate, x, current, direction, selected)
            if accepted is None:
                # A bad inverse-Hessian estimate need not terminate a sound
                # forward model: try steepest descent once with empty memory.
                pairs.clear()
                direction = _steepest_direction(current.gradient)
                accepted = _line_search(
                    evaluate,
                    x,
                    current,
                    direction,
                    selected,
                    initial_step=max(
                        selected.minimum_step, 2.0 ** -min(16, selected.line_search_evaluations)
                    ),
                )
        except _BudgetExhaustedError:
            reason = "evaluation_budget"
            break
        if accepted is None:
            reason = "line_search_failed"
            break
        alpha, new_value = accepted
        step: Vector = tuple(alpha * value for value in direction)
        new_x: Vector = tuple(a + b for a, b in zip(x, step, strict=True))
        change = tuple(a - b for a, b in zip(new_value.gradient, current.gradient, strict=True))
        curvature = _dot(step, change)
        threshold = 1e-12 * math.hypot(*step) * math.hypot(*change)
        if (
            math.isfinite(curvature)
            and curvature > max(0.0, threshold)
            and math.isfinite(1.0 / curvature)
        ):
            pairs.append((step, change, 1.0 / curvature))
            if len(pairs) > selected.memory:
                del pairs[0]
        old_loss = current.loss
        x, current = new_x, new_value
        record = IterationRecord(
            iteration, current.loss, _norm(current.gradient), alpha, calls, len(pairs)
        )
        history.append(record)
        if observe is not None:
            observe(record)
        if _norm(current.gradient) <= selected.gradient_tolerance:
            reason = "gradient_tolerance"
            break
        if _norm(step) <= selected.step_tolerance * max(1.0, _norm(x)):
            reason = "step_tolerance"
            break
        if abs(old_loss - current.loss) <= selected.relative_loss_tolerance * max(
            1.0, abs(old_loss), abs(current.loss)
        ):
            reason = "loss_stagnation"
            break
    return RecoveryResult(x, current, reason, calls, tuple(history))

Loss decrease and curvature serve different purposes in this search. Armijo requires at least a fixed fraction of the decrease predicted by the initial directional derivative, so a merely lower loss is insufficient. The strong-Wolfe condition also requires the magnitude of the directional derivative to fall sufficiently. Along a smooth descent line, satisfying that condition gives a positive inner product between the step and its gradient change in exact arithmetic, which is what the L-BFGS update needs. The maximum-step exception below does not guarantee that curvature condition, so the driver’s separate check of each stored pair still matters.

The line search scales the first steepest-descent direction, then brackets and bisects candidate step lengths to seek sufficient decrease and the strong-Wolfe curvature condition. At the configured maximum step, strict decrease satisfying Armijo is accepted even if the curvature condition is still unmet. A failed search may retry conservative steepest descent even before curvature memory exists, and exhausting the declared budget still leaves the incumbent intact. Out-of-chart and nonfinite trials cannot replace the incumbent. Input-contract errors still propagate as errors, rather than being mistaken for evidence that a smaller pose step would fix malformed data. The returned parameters and objective identify the last accepted state even if the most recent evaluation failed or was rejected.

Keep the accepted pose separate from the evaluator’s scratch image, which may contain a rejected trial. Re-render the returned pose when inspecting the final prediction. Record whether the solver stopped on a gradient, step, stagnation or budget criterion.

Log the objective at accepted poses, the norm of the gradient in scaled coordinates, the accepted step norm, and the number of forward and backward evaluations. Iteration count alone hides rejected proposals and line-search work. Record the data term separately from a prior so that apparent progress cannot be explained solely by a stronger preference for the initial pose.

Use scaled gradients, accepted step norms and actual objective changes to decide when further optimisation is useful. Evaluate the recovered geometry separately, as described next.

6.7 Evaluate recovery independently

Image agreement measures what the algorithm was asked to improve. A pose reference measures whether that improvement recovered the intended geometry.

Suppose the reference and estimated transforms use the same object origin, world frame and physical calibration, with rotations $\mathbf{R}_{\mathrm{ref}},\widehat{\mathbf{R}}$ and translations $\mathbf{t}_{\mathrm{ref}},\widehat{\mathbf{t}}$. Their translation and rotation discrepancies are

The errors in equation (6.16) have units of millimetres and radians. Clipping protects the inverse cosine from tiny floating-point excursions outside its domain. Translation error depends on the chosen object origin: changing that origin changes the translation accompanying a rotation error. Report the origin, particularly when comparing methods.

For application-relevant object points $\mathbf{x}_m^O$, a root-mean-square target registration error is

Unlike separate parameter errors, equation (6.17) measures displacement of specified physical targets. It exposes the consequence of a small angular error far from the rotation origin. Do not rigidly realign the estimated targets to the reference before evaluating this quantity, because that would remove the error the registration was meant to measure.

dpt.validation.recovery.pose_error() reports these quantities independently of the fitted image objective. It keeps origin displacement, rotation angle and target errors separate instead of adding millimetres to radians through an unexplained weight. The target set comes from the evaluation protocol, not from whichever points happened to align best.

def pose_error(
    recovered: RigidTransform,
    reference: RigidTransform,
    landmarks_object_mm: Sequence[Vector3],
) -> PoseError:
    """Report origin displacement, geodesic rotation and held-out target errors.

    A translation norm depends on the chosen object origin. Landmark errors also
    expose rotation about that origin, so they are reported rather than hidden
    inside a single weighted pose norm with arbitrary mixed units.
    """
    if not landmarks_object_mm:
        raise ContractError("recovery metrics need at least one declared landmark")
    # Relative rotation R_reference.T @ R_recovered, independent scalar indexing.
    relative = tuple(
        math.fsum(reference.rotation[3 * k + i] * recovered.rotation[3 * k + j] for k in range(3))
        for i in range(3)
        for j in range(3)
    )
    cosine = max(-1.0, min(1.0, (relative[0] + relative[4] + relative[8] - 1.0) / 2.0))
    sine = 0.5 * math.hypot(
        relative[7] - relative[5], relative[2] - relative[6], relative[3] - relative[1]
    )
    distances: list[float] = []
    for landmark in landmarks_object_mm:
        point = vector3(landmark, "landmark")
        a, b = recovered.point(point), reference.point(point)
        distances.append(math.hypot(*(x - y for x, y in zip(a, b, strict=True))))
    translation = math.hypot(
        *(a - b for a, b in zip(recovered.translation_mm, reference.translation_mm, strict=True))
    )
    rms = math.hypot(*(value / math.sqrt(len(distances)) for value in distances))
    if not math.isfinite(translation) or not math.isfinite(rms):
        raise NumericalError("recovery metric exceeds the finite coordinate range")
    return PoseError(translation, math.atan2(sine, cosine), rms, max(distances), len(distances))

The rotation calculation uses atan2 of sine and cosine information from the relative rotation. It measures the same principal angle as the clipped inverse-cosine expression in equation (6.16), while retaining more useful information near zero from the antisymmetric entries. Scaled hypot evaluations avoid squaring large coordinates unnecessarily. These small CPU calculations consume the returned rigid transform and a declared landmark set, and they do not derive a reference pose from the objective being judged.

A local example shows the role of target layout. Let $\mathbf{x}_m$ be reference target offsets from their centroid, expressed in world axes, so their mean is zero. Let $\delta\mathbf{t}$ be the displacement error at that centroid and $\delta\boldsymbol{\phi}$ a small world-frame rotation error. Then

The cross term vanishes in equation (6.18) because the targets are centred. Rotation about an axis barely moves targets near that axis and strongly moves distant ones. A target set concentrated in a small region can therefore report a small error while clinically relevant locations elsewhere move substantially.

Projected landmark error is useful for checking image alignment, but it measures a different quantity. Motion along a source ray may produce very little projected displacement. A low two-dimensional error cannot by itself bound the three-dimensional error of those landmarks.

Reference uncertainty should accompany the reported result. A tracked reference inherits tracking and calibration errors, while manually identified landmarks inherit localisation errors. Their agreement with a reconstruction is informative only to the precision those procedures support.

Keep the reference out of initialisation, hyperparameter choice and stopping when it is intended for final evaluation. Separate cases used to select those decisions from cases used to estimate recovery performance. For a capture study, report the starting-error distribution, success thresholds, number of trials and failures, alongside the distribution of final pose errors. Repeated noise realisations at one geometry answer a precision question, while independent geometries answer a different question about where recovery succeeds.

To test the assembled solver, generate a noiseless expected image from a prescribed pose of an analytic quadratic field, then recover it from several starting coordinates. Use a fresh optimiser history for each start and measure the returned transform with the geometric metrics above. Record both the numerical stopping measure and geometric error.

Using the same projector and discretisation for target generation isolates whether the solver’s pieces agree. After that check, vary acquisition assumptions and geometry to examine recovery under model mismatch.

6.8 When a lower loss is the wrong answer

Suppose the true pose is $\mathbf{G}_*$, but the observed mean includes an image component absent from the renderer. In the fixed whitening coordinates of the Gaussian objective, write

The mismatch $\mathbf{d}$ in equation (6.19) might arise from scatter, an exposure error, deformation or miscalibration. Here it denotes a mean discrepancy, and measurement noise can be considered separately. A small scaled pose displacement changes the prediction by approximately $\mathbf{J}_u\Delta\mathbf{u}$.

For a data-only least-squares fit with full-column-rank $\mathbf{J}_u$, the resulting first-order displacement is

In equation (6.20), the part of the mismatch that resembles a pose-induced image change drives the estimate away from the true pose. Weakly constrained directions can amplify that displacement. The formula is a local linearisation. Large shifts can leave the region in which its prediction is accurate.

The mismatch left after this fit is

The projection in equation (6.21) separates the mismatch into components that do and do not resemble pose changes (Figure 6.3), explaining why residual inspection alone can miss a systematic pose error. If $\mathbf{d}=\mathbf{J}_u\mathbf{b}$ for a small displacement $\mathbf{b}$, the linear model fits the entire mismatch by moving the object through $\mathbf{b}$. The residual vanishes. The pose error remains.

A slightly displaced detector origin provides a concrete example. Over a restricted field of view, its image shift may closely resemble an object translation. Letting the object move can then produce excellent image agreement while assigning the movement to the wrong physical component. A calibrated detector position or an additional informative view can distinguish which component moved.

When a fit improves the objective while worsening independent pose error, inspect which residual patterns the movement removed and which assumptions created them. Compare geometric calibration, represented anatomy and acquisition response before choosing a different optimisation algorithm. Chapter 7 develops the acquisition terms that can share an image direction with pose, and what remains identifiable when we estimate them together.