python/dpt/validation/recovery.py

1"""Independent geometric recovery metrics; objective decrease is a separate fact.23These CPU reports consume explicit small landmark sets and rigid transforms.4They never estimate a reference pose from the image used by the optimiser.5"""67from __future__ import annotations89import math10from collections.abc import Sequence11from dataclasses import dataclass1213from dpt.contracts import ContractError, NumericalError14from dpt.geometry import DetectorGeometry, RigidTransform, Vector3, vector3151617@dataclass(frozen=True, slots=True)18class PoseError:19    translation_mm: float20    rotation_radians: float21    landmark_rms_mm: float22    landmark_max_mm: float23    landmarks: int242526# region book:recovery-geometric-metrics27def pose_error(28    recovered: RigidTransform,29    reference: RigidTransform,30    landmarks_object_mm: Sequence[Vector3],31) -> PoseError:32    """Report origin displacement, geodesic rotation and held-out target errors.3334    A translation norm depends on the chosen object origin. Landmark errors also35    expose rotation about that origin, so they are reported rather than hidden36    inside a single weighted pose norm with arbitrary mixed units.37    """38    if not landmarks_object_mm:39        raise ContractError("recovery metrics need at least one declared landmark")40    # Relative rotation R_reference.T @ R_recovered, independent scalar indexing.41    relative = tuple(42        math.fsum(reference.rotation[3 * k + i] * recovered.rotation[3 * k + j] for k in range(3))43        for i in range(3)44        for j in range(3)45    )46    cosine = max(-1.0, min(1.0, (relative[0] + relative[4] + relative[8] - 1.0) / 2.0))47    sine = 0.5 * math.hypot(48        relative[7] - relative[5], relative[2] - relative[6], relative[3] - relative[1]49    )50    distances: list[float] = []51    for landmark in landmarks_object_mm:52        point = vector3(landmark, "landmark")53        a, b = recovered.point(point), reference.point(point)54        distances.append(math.hypot(*(x - y for x, y in zip(a, b, strict=True))))55    translation = math.hypot(56        *(a - b for a, b in zip(recovered.translation_mm, reference.translation_mm, strict=True))57    )58    rms = math.hypot(*(value / math.sqrt(len(distances)) for value in distances))59    if not math.isfinite(translation) or not math.isfinite(rms):60        raise NumericalError("recovery metric exceeds the finite coordinate range")61    return PoseError(translation, math.atan2(sine, cosine), rms, max(distances), len(distances))626364# endregion book:recovery-geometric-metrics656667def detector_coordinate(geometry: DetectorGeometry, point_world_mm: Vector3) -> tuple[float, float]:68    """Continuous (column, row) of a point's perspective projection.6970    Coordinates outside the detector remain meaningful residuals. A point on71    the source-parallel plane or behind the source is rejected; no clipped pixel72    is substituted for an undefined projection.73    """74    point = vector3(point_world_mm, "world landmark")75    u, v = geometry.u, geometry.v76    normal = (u[1] * v[2] - u[2] * v[1], u[2] * v[0] - u[0] * v[2], u[0] * v[1] - u[1] * v[0])77    offset = tuple(o - s for o, s in zip(geometry.origin_mm, geometry.source_mm, strict=True))78    ray = tuple(p - s for p, s in zip(point, geometry.source_mm, strict=True))79    plane_distance = math.fsum(a * b for a, b in zip(offset, normal, strict=True))80    direction = math.fsum(a * b for a, b in zip(ray, normal, strict=True))81    if direction == 0 or not math.isfinite(direction):82        raise ContractError("landmark projection is undefined on a source-parallel plane")83    scale = plane_distance / direction84    if not math.isfinite(scale) or scale <= 0:85        raise ContractError("landmark must project onto the forward detector plane")86    hit = tuple(scale * r - o for r, o in zip(ray, offset, strict=True))87    coordinate = (88        math.fsum(a * b for a, b in zip(hit, u, strict=True)) / geometry.spacing_mm[0],89        math.fsum(a * b for a, b in zip(hit, v, strict=True)) / geometry.spacing_mm[1],90    )91    if not all(map(math.isfinite, coordinate)):92        raise NumericalError("projected landmark exceeds the finite pixel range")93    return coordinate949596def reprojection_rms_pixels(97    recovered: RigidTransform,98    reference: RigidTransform,99    geometry: DetectorGeometry,100    landmarks_object_mm: Sequence[Vector3],101) -> float:102    if not landmarks_object_mm:103        raise ContractError("reprojection metrics need at least one declared landmark")104    differences: list[float] = []105    for landmark in landmarks_object_mm:106        point = vector3(landmark, "landmark")107        a = detector_coordinate(geometry, recovered.point(point))108        b = detector_coordinate(geometry, reference.point(point))109        differences.extend(x - y for x, y in zip(a, b, strict=True))110    result = math.hypot(*(x / math.sqrt(len(landmarks_object_mm)) for x in differences))111    if not math.isfinite(result):112        raise NumericalError("reprojection residual exceeds the finite pixel range")113    return result114