From c0da896d141c4b10e1db100d16b03403b406a7c2 Mon Sep 17 00:00:00 2001 From: HarjotDhanota <152345739+HarjotDhanota@users.noreply.github.com> Date: Fri, 28 Aug 2026 11:54:29 -0700 Subject: [PATCH 1/2] feat(camera): clubhead-pose geometry, mesh tooling and teed-ball fit Club mesh with provenance-checked normalization, calibrated camera model with silhouette projection, and sub-pixel teed-ball location that anchors the world frame. First stage of the clubface impact-location work; not yet wired into the shot pipeline. --- docs/CHANGELOG.md | 7 + scripts/analysis/download_club_mesh.py | 185 ++++++ src/openflight/camera/clubpose/__init__.py | 9 + src/openflight/camera/clubpose/mesh.py | 563 ++++++++++++++++++ .../camera/clubpose/meshes/SOURCES.md | 82 +++ .../camera/clubpose/meshes/assets/.gitignore | 2 + src/openflight/camera/clubpose/projection.py | 330 ++++++++++ src/openflight/camera/clubpose/teed_ball.py | 321 ++++++++++ tests/test_clubpose_ball_detect.py | 144 +++++ 9 files changed, 1643 insertions(+) create mode 100644 scripts/analysis/download_club_mesh.py create mode 100644 src/openflight/camera/clubpose/__init__.py create mode 100644 src/openflight/camera/clubpose/mesh.py create mode 100644 src/openflight/camera/clubpose/meshes/SOURCES.md create mode 100644 src/openflight/camera/clubpose/meshes/assets/.gitignore create mode 100644 src/openflight/camera/clubpose/projection.py create mode 100644 src/openflight/camera/clubpose/teed_ball.py create mode 100644 tests/test_clubpose_ball_detect.py diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index aa528d58..fd1a42f4 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -7,6 +7,13 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +### Added +- **Clubhead-pose geometry (`openflight.camera.clubpose`).** The club mesh and + its normalization tooling, the calibrated camera model with silhouette + projection, and sub-pixel location of the teed ball, which anchors the world + frame and the camera range. First stage of the clubface impact-location work; + not yet wired into the shot pipeline. + ### Added - **Automatic OV9281 exposure control.** High-speed camera capture now measures the impact area every five seconds, restores the last known-good setting at diff --git a/scripts/analysis/download_club_mesh.py b/scripts/analysis/download_club_mesh.py new file mode 100644 index 00000000..c306ca39 --- /dev/null +++ b/scripts/analysis/download_club_mesh.py @@ -0,0 +1,185 @@ +"""Acquire and normalize the club meshes without vendoring them. + +The 7-iron used by this research is a GrabCAD community model. It is not +redistributed, so you fetch your own copy under GrabCAD's terms and point this +script at it; see SOURCES.md for the link, the expected SHA-256 and the licence +position. Run from the repository root: + + uv run python scripts/analysis/download_club_mesh.py \n --local-iron "/path/to/690CB 7-iron.STL" + +""" + +from __future__ import annotations + +import argparse +import json +import sys +from pathlib import Path +from typing import Any + +ROOT = Path(__file__).resolve().parents[2] +sys.path.insert(0, str(ROOT / "src")) + +from openflight.camera.clubpose.mesh import ( # noqa: E402 + ACTIVE_MESH_SOURCES, + CATEGORY_DIMENSIONS_MM, + MESH_SOURCES, + MeshSource, + admit_mesh, + default_mesh_asset_root, + detect_face_plane, + face_detection_record, + load_binary_stl, + load_normalized_mesh, + normalize_clubhead, + save_normalized_mesh, +) + +_NORMALIZATION_VERSION = "geometric-face-anchor-v2" + + +def validate_source_metadata(source: MeshSource, metadata: dict[str, Any]) -> None: + """Fail closed if identity, author, download status, or license drifted.""" + if str(metadata.get("uid")) != source.uid or str(metadata.get("name")) != source.name: + raise ValueError(f"source identity changed for {source.club}") + if not bool(metadata.get("isDownloadable")): + raise ValueError(f"source is no longer downloadable for {source.club}") + author = str(metadata.get("user", {}).get("displayName", "")) + if author != source.author: + raise ValueError(f"source author changed for {source.club}: {author!r}") + license_payload = metadata.get("license", {}) + license_url = str(license_payload.get("url", "")).replace("http://", "https://") + if str(license_payload.get("label")) != "CC Attribution" or license_url.rstrip( + "/" + ) != source.license_url.rstrip("/"): + raise ValueError(f"source license changed for {source.club}") + + +def import_local_stl( + source_path: Path | str, + output_root: Path, + *, + expected_sha256: str | None = None, +) -> dict[str, Any]: + """Import the registered maintainer-local 7-iron without copying its STL.""" + source = MESH_SOURCES["poc_7iron"] + registered_hash = source.expected_source_sha256 + if expected_sha256 is not None and registered_hash is not None: + if expected_sha256.lower() != registered_hash.lower(): + raise ValueError("caller SHA-256 does not match the frozen local-source registration") + required_hash = expected_sha256 or registered_hash + loaded = load_binary_stl(source_path, source_uid=source.uid, expected_sha256=required_hash) + admission = admit_mesh( + loaded, + category_dimensions_mm=CATEGORY_DIMENSIONS_MM[source.club], + source_units_mm=True, + ) + if not admission.accepted: + raise ValueError(f"mesh admission failed for {source.club}: {admission.reasons}") + assert admission.face is not None + normalized = normalize_clubhead( + loaded, + CATEGORY_DIMENSIONS_MM[source.club], + source_units_mm=True, + ) + normalized_face = detect_face_plane(normalized) + asset_metadata = { + "source_uid": source.uid, + "source_name": source.name, + "author": source.author, + "page_url": source.page_url, + "license_spdx": source.license_spdx, + "license_url": source.license_url, + "source_file_sha256": loaded.source_sha256, + "download_format": "binary_stl_maintainer_local", + "redistribution": "prohibited; local research use only", + "normalization": _NORMALIZATION_VERSION, + "source_units_mm": True, + "category_dimensions_mm": CATEGORY_DIMENSIONS_MM[source.club], + "geometry_sha256": admission.geometry_sha256, + "component_count_after_weld": admission.component_count, + "boundary_edge_count_after_weld": admission.boundary_edge_count, + "boundary_edge_fraction_after_weld": admission.boundary_edge_fraction, + "dimensions_before_normalization_mm": admission.dimensions_mm, + "face_detection_source": face_detection_record(admission.face), + "face_detection_normalized": face_detection_record(normalized_face), + "source_vertex_count": int(len(loaded.vertices_local_mm)), + "source_triangle_count": int(len(loaded.faces)), + "clubhead_vertex_count": int(len(normalized.vertices_local_mm)), + "clubhead_triangle_count": int(len(normalized.faces)), + "trademark_note": "synthetic truth only; no Titleist endorsement implied", + } + asset_path = output_root / f"{source.club}.npz" + asset_sha256 = save_normalized_mesh(asset_path, normalized, asset_metadata) + record = {**asset_metadata, "asset_path": asset_path.name, "asset_sha256": asset_sha256} + print(json.dumps(record, indent=2, sort_keys=True)) + return record + + +def _existing_record(source: MeshSource, output_root: Path) -> dict[str, Any] | None: + asset_path = output_root / f"{source.club}.npz" + if not asset_path.is_file(): + return None + mesh, metadata, asset_sha256 = load_normalized_mesh(str(asset_path.resolve())) + if mesh.source_uid != source.uid: + raise ValueError(f"cached source identity mismatch for {source.club}") + if source.expected_source_sha256 is not None and ( + mesh.source_sha256 != source.expected_source_sha256 + ): + raise ValueError(f"cached source SHA-256 mismatch for {source.club}") + if source.expected_asset_sha256 is not None and (asset_sha256 != source.expected_asset_sha256): + raise ValueError(f"cached asset SHA-256 mismatch for {source.club}") + if metadata.get("normalization") != _NORMALIZATION_VERSION: + return None + return {**metadata, "asset_path": asset_path.name, "asset_sha256": asset_sha256} + + +def main() -> int: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--output", type=Path, default=default_mesh_asset_root()) + parser.add_argument("--local-iron", type=Path) + args = parser.parse_args() + args.output.mkdir(parents=True, exist_ok=True) + iron = ACTIVE_MESH_SOURCES["poc_7iron"] + iron_record = _existing_record(iron, args.output) + if iron_record is None: + if args.local_iron is None: + parser.error("--local-iron is required to import the missing maintainer-local 690CB") + stl = Path(args.local_iron).expanduser() + if not stl.is_file(): + parser.error( + f"no STL at {stl}. Fetch the 690CB 7-iron from the GrabCAD page in " + "src/openflight/camera/clubpose/meshes/SOURCES.md (free account, " + "their terms) and point --local-iron at the downloaded file." + ) + iron_record = import_local_stl(stl, args.output) + records = [iron_record] + (args.output / "manifest.json").write_text( + json.dumps( + { + "sources": records, + "retired_sources": [ + { + "club": source.club, + "source_uid": source.uid, + "status": source.status, + "reason": source.status_reason, + } + for source in MESH_SOURCES.values() + if source.status != "active" + ], + }, + indent=2, + sort_keys=True, + ) + + "\n", + encoding="utf-8", + ) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) + (admit_mesh,) + (detect_face_plane,) + (face_detection_record,) diff --git a/src/openflight/camera/clubpose/__init__.py b/src/openflight/camera/clubpose/__init__.py new file mode 100644 index 00000000..41bbea13 --- /dev/null +++ b/src/openflight/camera/clubpose/__init__.py @@ -0,0 +1,9 @@ +"""Clubhead pose from the behind-ball camera and the radars. + +This package lands in stages. This stage carries the geometry: the club mesh +(mesh), the camera model and silhouette projection (projection), and sub-pixel +location of the teed ball (teed_ball), which anchors the world frame. The pose +fit that consumes them follows, and nothing here is wired into the shot +pipeline yet. See docs/clubface-impact-location-report.md for what is +validated and what is not. +""" diff --git a/src/openflight/camera/clubpose/mesh.py b/src/openflight/camera/clubpose/mesh.py new file mode 100644 index 00000000..e6d188c6 --- /dev/null +++ b/src/openflight/camera/clubpose/mesh.py @@ -0,0 +1,563 @@ +"""License-pinned mesh acquisition helpers and a NumPy silhouette renderer.""" + +from __future__ import annotations + +import base64 +import hashlib +import json +import math +import struct +import zipfile +from dataclasses import dataclass +from functools import lru_cache +from pathlib import Path +from typing import Any + +import numpy as np + + +@dataclass(frozen=True) +class MeshSource: + club: str + uid: str + name: str + author: str + page_url: str + license_spdx: str + license_url: str + downloadable: bool + published_triangles: int + source_kind: str = "maintainer_local_binary_stl" + expected_source_sha256: str | None = None + expected_asset_sha256: str | None = None + status: str = "active" + status_reason: str = "" + + +MESH_SOURCES = { + "poc_7iron": MeshSource( + club="poc_7iron", + uid="grabcad:titleist-7-iron-golf-club-1:690cb-right-handed", + name="Titleist 690CB 7-iron golf club", + author="GrabCAD Community contributor", + page_url="https://grabcad.com/library/titleist-7-iron-golf-club-1", + license_spdx="LicenseRef-GrabCAD-Local-Research-Only", + license_url=("https://help.grabcad.com/article/246-how-can-models-be-used-and-shared"), + downloadable=False, + published_triangles=26_238, + source_kind="maintainer_local_binary_stl", + expected_source_sha256=("f35936799295e6ce344279e557f0265ccbb8acef69c4508daff80d219d03cb85"), + ), +} + +ACTIVE_MESH_SOURCES = { + club: source for club, source in MESH_SOURCES.items() if source.status == "active" +} +CATEGORY_DIMENSIONS_MM = { + "poc_7iron": {"width": 80.0, "height": 50.0, "depth": 38.0}, +} + + +@dataclass(frozen=True) +class TriangleMesh: + """Triangle mesh in club-local coordinates: +x depth, +y width, +z height.""" + + vertices_local_mm: np.ndarray + faces: np.ndarray + source_uid: str + source_sha256: str + + def __post_init__(self) -> None: + vertices = np.asarray(self.vertices_local_mm, dtype=float) + faces = np.asarray(self.faces, dtype=np.int32) + if vertices.ndim != 2 or vertices.shape[1] != 3: + raise ValueError("mesh vertices must have shape [N,3]") + if faces.ndim != 2 or faces.shape[1] != 3: + raise ValueError("mesh faces must have shape [M,3]") + if faces.size and (int(faces.min()) < 0 or int(faces.max()) >= len(vertices)): + raise ValueError("mesh face index is outside the vertex array") + if not np.all(np.isfinite(vertices)): + raise ValueError("mesh vertices must be finite") + object.__setattr__(self, "vertices_local_mm", vertices) + object.__setattr__(self, "faces", faces) + + +@dataclass(frozen=True) +class FacePlaneDetection: + normal_source: np.ndarray + width_axis_source: np.ndarray + height_axis_source: np.ndarray + centroid_source: np.ndarray + coherent_area_mm2: float + face_span_mm: np.ndarray + triangle_indices: np.ndarray + + +@dataclass(frozen=True) +class MeshAdmission: + accepted: bool + reasons: tuple[str, ...] + component_count: int + boundary_edge_count: int + boundary_edge_fraction: float + geometry_sha256: str + dimensions_mm: dict[str, float] + face: FacePlaneDetection | None + + +def face_detection_record(face: FacePlaneDetection) -> dict[str, Any]: + """Serialize the required face-plane provenance without embedding geometry.""" + return { + "normal": face.normal_source.tolist(), + "coherent_area_mm2": face.coherent_area_mm2, + "face_span_mm": face.face_span_mm.tolist(), + "triangle_count": int(len(face.triangle_indices)), + } + + +_COMPONENT_DTYPES = { + 5120: np.dtype("i1"), + 5121: np.dtype("u1"), + 5122: np.dtype(" np.ndarray: + accessor = document["accessors"][index] + if "sparse" in accessor: + raise ValueError("sparse glTF accessors are not supported") + view = document["bufferViews"][accessor["bufferView"]] + dtype = _COMPONENT_DTYPES[int(accessor["componentType"])] + width = _TYPE_WIDTHS[str(accessor["type"])] + count = int(accessor["count"]) + offset = int(view.get("byteOffset", 0)) + int(accessor.get("byteOffset", 0)) + stride = int(view.get("byteStride", dtype.itemsize * width)) + raw = buffers[int(view["buffer"])] + if stride == dtype.itemsize * width: + values = np.frombuffer(raw, dtype=dtype, count=count * width, offset=offset) + return values.reshape(count, width).copy() + output = np.empty((count, width), dtype=dtype) + for row in range(count): + output[row] = np.frombuffer(raw, dtype=dtype, count=width, offset=offset + row * stride) + return output + + +def _quaternion_matrix(value: list[float]) -> np.ndarray: + x, y, z, w = (float(item) for item in value) + norm = math.sqrt(x * x + y * y + z * z + w * w) + if norm == 0.0: + return np.eye(4) + x, y, z, w = x / norm, y / norm, z / norm, w / norm + return np.array( + [ + [1 - 2 * (y * y + z * z), 2 * (x * y - z * w), 2 * (x * z + y * w), 0], + [2 * (x * y + z * w), 1 - 2 * (x * x + z * z), 2 * (y * z - x * w), 0], + [2 * (x * z - y * w), 2 * (y * z + x * w), 1 - 2 * (x * x + y * y), 0], + [0, 0, 0, 1], + ], + dtype=float, + ) + + +def _node_matrix(node: dict[str, Any]) -> np.ndarray: + if "matrix" in node: + return np.asarray(node["matrix"], dtype=float).reshape(4, 4, order="F") + translation = np.eye(4) + translation[:3, 3] = np.asarray(node.get("translation", [0, 0, 0]), dtype=float) + scale = np.eye(4) + scale[np.arange(3), np.arange(3)] = np.asarray(node.get("scale", [1, 1, 1]), dtype=float) + return translation @ _quaternion_matrix(node.get("rotation", [0, 0, 0, 1])) @ scale + + +def _archive_buffer(bundle: zipfile.ZipFile, gltf_dir: str, uri: str) -> bytes: + if uri.startswith("data:"): + return base64.b64decode(uri.split(",", 1)[1]) + return bundle.read(str(Path(gltf_dir, uri)).replace("\\", "/")) + + +def load_gltf_archive( + archive_path: Path | str, *, source_uid: str, source_sha256: str +) -> TriangleMesh: + """Load triangle primitives and scene-node transforms from a Sketchfab glTF ZIP.""" + with zipfile.ZipFile(archive_path) as bundle: + gltf_names = sorted(name for name in bundle.namelist() if name.lower().endswith(".gltf")) + if len(gltf_names) != 1: + raise ValueError("download archive must contain exactly one .gltf scene") + gltf_name = gltf_names[0] + document = json.loads(bundle.read(gltf_name)) + gltf_dir = str(Path(gltf_name).parent) + buffers = [ + _archive_buffer(bundle, gltf_dir, str(item["uri"])) for item in document["buffers"] + ] + + vertices: list[np.ndarray] = [] + faces: list[np.ndarray] = [] + + def visit(node_index: int, parent: np.ndarray) -> None: + node = document["nodes"][node_index] + world = parent @ _node_matrix(node) + if "mesh" in node: + mesh = document["meshes"][int(node["mesh"])] + for primitive in mesh["primitives"]: + if int(primitive.get("mode", 4)) != 4: + raise ValueError("only glTF TRIANGLES primitives are supported") + position = _read_accessor( + document, buffers, int(primitive["attributes"]["POSITION"]) + ).astype(float) + transformed = np.column_stack([position, np.ones(len(position))]) @ world.T + transformed = transformed[:, :3] / transformed[:, 3, None] + if "indices" in primitive: + triangle = _read_accessor(document, buffers, int(primitive["indices"])).reshape( + -1 + ) + else: + triangle = np.arange(len(position), dtype=np.int32) + if len(triangle) % 3: + raise ValueError("triangle index count is not divisible by three") + offset = sum(len(item) for item in vertices) + vertices.append(transformed) + faces.append(triangle.reshape(-1, 3).astype(np.int32) + offset) + for child in node.get("children", []): + visit(int(child), world) + + scene_index = int(document.get("scene", 0)) + for root in document["scenes"][scene_index].get("nodes", []): + visit(int(root), np.eye(4)) + if not vertices or not faces: + raise ValueError("glTF scene contains no triangle geometry") + return TriangleMesh(np.vstack(vertices), np.vstack(faces), source_uid, source_sha256) + + +def load_binary_stl( + path: Path | str, *, source_uid: str, expected_sha256: str | None +) -> TriangleMesh: + """Decode a binary STL after an optional fail-closed source-hash check.""" + payload = Path(path).read_bytes() + digest = hashlib.sha256(payload).hexdigest() + if expected_sha256 is not None and digest.lower() != expected_sha256.lower(): + raise ValueError( + f"binary STL SHA-256 mismatch: expected {expected_sha256.lower()}, got {digest}" + ) + if len(payload) < 84: + raise ValueError("binary STL is shorter than its 84-byte header") + triangle_count = struct.unpack_from(" list[np.ndarray]: + # Weld identical seam vertices before topology traversal. glTF material + # primitives commonly duplicate vertices along an otherwise connected shell. + _, welded = np.unique(np.round(vertices, decimals=9), axis=0, return_inverse=True) + welded_faces = welded[faces] + parent = np.arange(int(welded.max()) + 1, dtype=np.int32) + + def root(index: int) -> int: + while parent[index] != index: + parent[index] = parent[parent[index]] + index = int(parent[index]) + return index + + for a, b, c in welded_faces: + anchor = root(int(a)) + for other in (int(b), int(c)): + other_root = root(other) + if anchor != other_root: + parent[other_root] = anchor + labels = np.array([root(int(face[0])) for face in welded_faces], dtype=np.int32) + return [np.flatnonzero(labels == label) for label in np.unique(labels)] + + +def _welded_faces(vertices: np.ndarray, faces: np.ndarray) -> tuple[np.ndarray, np.ndarray]: + welded_vertices, welded = np.unique( + np.round(np.asarray(vertices, dtype=float), decimals=8), axis=0, return_inverse=True + ) + return welded_vertices, welded[np.asarray(faces, dtype=np.int32)] + + +def geometry_hash(mesh: TriangleMesh) -> str: + """Hash triangle geometry independent of source vertex/face ordering.""" + triangles = np.round(mesh.vertices_local_mm[mesh.faces], decimals=8) + canonical_triangles = [] + for triangle in triangles: + order = np.lexsort((triangle[:, 2], triangle[:, 1], triangle[:, 0])) + canonical_triangles.append(triangle[order].reshape(-1)) + canonical = np.asarray(canonical_triangles) + order = np.lexsort(tuple(canonical[:, column] for column in reversed(range(9)))) + return hashlib.sha256(np.ascontiguousarray(canonical[order]).tobytes()).hexdigest() + + +def _triangle_adjacency(welded_faces: np.ndarray) -> list[set[int]]: + edge_owners: dict[tuple[int, int], list[int]] = {} + for face_index, face in enumerate(welded_faces): + for first, second in ((face[0], face[1]), (face[1], face[2]), (face[2], face[0])): + edge = tuple(sorted((int(first), int(second)))) + edge_owners.setdefault(edge, []).append(face_index) + adjacency = [set() for _ in welded_faces] + for owners in edge_owners.values(): + for owner in owners: + adjacency[owner].update(other for other in owners if other != owner) + return adjacency + + +def detect_face_plane( + mesh: TriangleMesh, + *, + normal_tolerance_deg: float = 15.0, + aspect_bounds: tuple[float, float] = (0.35, 0.65), +) -> FacePlaneDetection: + """Find the largest coherent extremity plane with a clubface-like lens aspect.""" + vertices = mesh.vertices_local_mm + triangles = vertices[mesh.faces] + cross = np.cross(triangles[:, 1] - triangles[:, 0], triangles[:, 2] - triangles[:, 0]) + double_area = np.linalg.norm(cross, axis=1) + valid = double_area > 1e-10 + normals = np.zeros_like(cross) + normals[valid] = cross[valid] / double_area[valid, None] + areas = double_area / 2.0 + _, welded_faces = _welded_faces(vertices, mesh.faces) + adjacency = _triangle_adjacency(welded_faces) + cosine_limit = math.cos(math.radians(normal_tolerance_deg)) + unassigned = set(np.flatnonzero(valid).tolist()) + candidates: list[FacePlaneDetection] = [] + mesh_center = np.mean(vertices, axis=0) + while unassigned: + seed = max(unassigned, key=lambda index: float(areas[index])) + seed_normal = normals[seed] + region = set() + pending = [seed] + while pending: + current = pending.pop() + if current not in unassigned: + continue + if abs(float(normals[current] @ seed_normal)) < cosine_limit: + continue + unassigned.remove(current) + region.add(current) + pending.extend(adjacency[current] & unassigned) + if not region: + continue + indices = np.asarray(sorted(region), dtype=np.int32) + aligned = normals[indices] * np.sign(normals[indices] @ seed_normal)[:, None] + normal = np.sum(aligned * areas[indices, None], axis=0) + normal /= np.linalg.norm(normal) + vertex_indices = np.unique(mesh.faces[indices].reshape(-1)) + points = vertices[vertex_indices] + centroid = np.average(np.mean(triangles[indices], axis=1), weights=areas[indices], axis=0) + if float(normal @ (centroid - mesh_center)) < 0.0: + normal *= -1.0 + centered = points - centroid + planar = centered - np.outer(centered @ normal, normal) + _, _, axes = np.linalg.svd(planar, full_matrices=False) + width_axis = axes[0] - normal * float(axes[0] @ normal) + width_axis /= np.linalg.norm(width_axis) + if width_axis[int(np.argmax(np.abs(width_axis)))] < 0.0: + width_axis *= -1.0 + height_axis = np.cross(normal, width_axis) + spans = np.array([np.ptp(points @ width_axis), np.ptp(points @ height_axis)], dtype=float) + if spans[1] > spans[0]: + spans = spans[::-1] + width_axis, height_axis = height_axis, -width_axis + aspect = float(spans[1] / max(spans[0], 1e-12)) + projection = vertices @ normal + extremity_distance = min( + abs(float(centroid @ normal) - float(np.min(projection))), + abs(float(np.max(projection)) - float(centroid @ normal)), + ) + depth = float(np.ptp(projection)) + flatness = float(np.max(np.abs(centered @ normal))) + if ( + aspect_bounds[0] <= aspect <= aspect_bounds[1] + and extremity_distance <= max(1.0, 0.10 * depth) + and flatness <= max(1.0, 0.04 * spans[0]) + ): + candidates.append( + FacePlaneDetection( + normal_source=normal, + width_axis_source=width_axis, + height_axis_source=height_axis, + centroid_source=centroid, + coherent_area_mm2=float(np.sum(areas[indices])), + face_span_mm=spans, + triangle_indices=indices, + ) + ) + if not candidates: + raise ValueError("no coherent extremity plane has the registered clubface lens aspect") + return max(candidates, key=lambda item: item.coherent_area_mm2) + + +def _boundary_edge_count(mesh: TriangleMesh) -> int: + _, faces = _welded_faces(mesh.vertices_local_mm, mesh.faces) + counts: dict[tuple[int, int], int] = {} + for face in faces: + for first, second in ((face[0], face[1]), (face[1], face[2]), (face[2], face[0])): + edge = tuple(sorted((int(first), int(second)))) + counts[edge] = counts.get(edge, 0) + 1 + return sum(count != 2 for count in counts.values()) + + +def admit_mesh( + mesh: TriangleMesh, + *, + category_dimensions_mm: dict[str, float], + source_units_mm: bool, + tolerance_fraction: float = 0.15, +) -> MeshAdmission: + """Apply the frozen CAD-corpus admission checks before normalization.""" + components = _connected_face_components(mesh.vertices_local_mm, mesh.faces) + boundary_edges = _boundary_edge_count(mesh) + edge_denominator = max(1.0, 1.5 * len(mesh.faces)) + boundary_fraction = boundary_edges / edge_denominator + reasons = [] + if len(components) != 1: + reasons.append("component_count") + if boundary_fraction > 0.001: + reasons.append("open_boundary") + try: + face = detect_face_plane(mesh) + except ValueError: + face = None + reasons.append("face_plane_missing") + dimensions: dict[str, float] = {} + if face is not None: + dimensions = { + "width": float(face.face_span_mm[0]), + "height": float(face.face_span_mm[1]), + "depth": float(np.ptp(mesh.vertices_local_mm @ face.normal_source)), + } + if source_units_mm: + for name, nominal in category_dimensions_mm.items(): + relative = abs(dimensions[name] - float(nominal)) / float(nominal) + if relative > tolerance_fraction + 0.001: + reasons.append(f"dimension_{name}") + return MeshAdmission( + accepted=not reasons, + reasons=tuple(reasons), + component_count=len(components), + boundary_edge_count=boundary_edges, + boundary_edge_fraction=boundary_fraction, + geometry_sha256=geometry_hash(mesh), + dimensions_mm=dimensions, + face=face, + ) + + +def normalize_clubhead( + mesh: TriangleMesh, + dimensions_mm: dict[str, float], + *, + source_units_mm: bool = False, +) -> TriangleMesh: + """Anchor axes to the detected face plane; preserve trusted metric CAD scale.""" + components = _connected_face_components(mesh.vertices_local_mm, mesh.faces) + if len(components) != 1: + raise ValueError("mesh admission requires one welded connected component") + face = detect_face_plane(mesh) + axes = np.stack([face.normal_source, face.width_axis_source, face.height_axis_source]) + local = mesh.vertices_local_mm @ axes.T + local -= (np.min(local, axis=0) + np.max(local, axis=0)) / 2.0 + if not source_units_mm: + target = np.array( + [dimensions_mm["depth"], dimensions_mm["width"], dimensions_mm["height"]], + dtype=float, + ) + local *= target / np.ptp(local, axis=0) + normalized = TriangleMesh(local, mesh.faces, mesh.source_uid, mesh.source_sha256) + normalized_face = detect_face_plane(normalized) + angle = math.degrees( + math.acos(float(np.clip(normalized_face.normal_source @ np.array([1.0, 0.0, 0.0]), -1, 1))) + ) + if angle > 3.0: + raise ValueError(f"normalized face normal invariant failed: {angle:.3f} degrees") + return normalized + + +def rasterize_projected_triangles( + vertices_uv: np.ndarray, faces: np.ndarray, *, width: int, height: int +) -> np.ndarray: + """Rasterize the union of projected triangles with NumPy scanline intervals.""" + vertices = np.asarray(vertices_uv, dtype=float) + triangles = vertices[np.asarray(faces, dtype=np.int32)] + finite = np.all(np.isfinite(triangles), axis=(1, 2)) + triangles = triangles[finite] + mask = np.zeros((height, width), dtype=bool) + if not len(triangles): + return mask + min_y = np.min(triangles[:, :, 1], axis=1) + max_y = np.max(triangles[:, :, 1], axis=1) + edges_a = triangles[:, [0, 1, 2]] + edges_b = triangles[:, [1, 2, 0]] + for row in range(height): + y = row + 0.5 + active = (min_y <= y) & (max_y >= y) + if not np.any(active): + continue + a = edges_a[active] + b = edges_b[active] + dy = b[:, :, 1] - a[:, :, 1] + crosses = (np.abs(dy) > 1e-12) & ( + (y >= np.minimum(a[:, :, 1], b[:, :, 1])) & (y <= np.maximum(a[:, :, 1], b[:, :, 1])) + ) + safe_dy = np.where(crosses, dy, 1.0) + intersections = a[:, :, 0] + (y - a[:, :, 1]) * (b[:, :, 0] - a[:, :, 0]) / safe_dy + intersections = np.where(crosses, intersections, np.nan) + with np.errstate(all="ignore"): + left = np.nanmin(intersections, axis=1) + right = np.nanmax(intersections, axis=1) + valid = np.isfinite(left) & np.isfinite(right) + starts = np.maximum(0, np.ceil(left[valid] - 0.5).astype(int)) + ends = np.minimum(width - 1, np.floor(right[valid] - 0.5).astype(int)) + visible = starts <= ends + difference = np.zeros(width + 1, dtype=np.int32) + np.add.at(difference, starts[visible], 1) + np.add.at(difference, ends[visible] + 1, -1) + mask[row] = np.cumsum(difference[:-1]) > 0 + return mask + + +def save_normalized_mesh(path: Path | str, mesh: TriangleMesh, metadata: dict[str, Any]) -> str: + """Write a deterministic local cache and return its content SHA-256.""" + path = Path(path) + path.parent.mkdir(parents=True, exist_ok=True) + np.savez( + path, + vertices_local_mm=mesh.vertices_local_mm, + faces=mesh.faces, + source_uid=np.asarray(mesh.source_uid), + source_sha256=np.asarray(mesh.source_sha256), + metadata_json=np.asarray(json.dumps(metadata, sort_keys=True, separators=(",", ":"))), + ) + return hashlib.sha256(path.read_bytes()).hexdigest() + + +@lru_cache(maxsize=8) +def load_normalized_mesh(path: str) -> tuple[TriangleMesh, dict[str, Any], str]: + payload = np.load(path, allow_pickle=False) + mesh = TriangleMesh( + payload["vertices_local_mm"], + payload["faces"], + str(payload["source_uid"]), + str(payload["source_sha256"]), + ) + metadata = json.loads(str(payload["metadata_json"])) + return mesh, metadata, hashlib.sha256(Path(path).read_bytes()).hexdigest() + + +def default_mesh_asset_root() -> Path: + return Path(__file__).resolve().parent / "meshes" / "assets" diff --git a/src/openflight/camera/clubpose/meshes/SOURCES.md b/src/openflight/camera/clubpose/meshes/SOURCES.md new file mode 100644 index 00000000..b553a218 --- /dev/null +++ b/src/openflight/camera/clubpose/meshes/SOURCES.md @@ -0,0 +1,82 @@ +# Phase F1 club-mesh sources and licenses + +No third-party mesh is committed in this directory. Acquisition uses the +publisher's authenticated download endpoint, records the returned archive hash, +and leaves downloaded files ignored by Git. + +## Selected and retired sources + +| Club | Source | Model ID | License | Published geometry | +|---|---|---|---|---:| +| Driver (RETIRED) | [Callaway Maverik Golf Driver](https://sketchfab.com/3d-models/callaway-maverik-golf-driver-978d0740dc514c8695bbb02f4083f0e3), Paul Ekins | `978d0740dc514c8695bbb02f4083f0e3` | [CC BY 4.0](https://creativecommons.org/licenses/by/4.0/) | 41,855 triangles | +| 7-iron | [Titleist 7-iron golf club](https://grabcad.com/library/titleist-7-iron-golf-club-1), GrabCAD Community contributor; maintainer-downloaded 690CB right-handed STL | `grabcad:titleist-7-iron-golf-club-1:690cb-right-handed` | Local research use only; no redistribution | 26,238 triangles | + +The Sketchfab v3 metadata reports the driver as downloadable and labels its +license `CC Attribution`, with author credit required and commercial use allowed. +CC BY 4.0 permits sharing and adaptation with attribution. We still use a local +cache rather than vendoring the archive: it preserves source provenance, does not +add a third-party binary payload to OpenFlight, and follows the design +specification's existing no-redistribution boundary. + +The right-handed Titleist 690CB source is a maintainer-supplied, millimetre-scaled +binary STL from the named GrabCAD model. Its source SHA-256 is +`f35936799295e6ce344279e557f0265ccbb8acef69c4508daff80d219d03cb85`. +It is **local use only**: neither the STL nor its normalized NPZ may be committed. +Only its provenance, hash, attribution, and aggregate evaluation results enter +the repository. The local importer checks this hash before parsing the STL. + +The downloaded models are used only as synthetic truth. Their names or geometry +do not imply endorsement by Callaway Golf or Titleist. + +### Post-F1 source-quality correction + +The Maverik is retired and excluded from active manifests and evaluations. It is +a posed art scene containing grass, a ball, and tens of thousands of disconnected +shell components. Its face/sole geometry is ambiguous, and the former PCA +extent-order normalizer assigned its face normal to the height axis and then +anisotropically distorted the head using an incorrect 55 mm driver depth. It is +not salvageable as canonical driver truth. Driver arms are `HOLD_CAD_MESH` until +the maintainer supplies a locally admitted CAD driver; the corrected driver +category references are 118 mm width, 60 mm height, and 112 mm depth. + +The 690CB was re-imported from the same pinned STL using geometric face anchoring +and trusted source millimetres. Corrected normalized asset SHA-256 is +`d63bf7cf1224eb9ce0c7480967057201a4843f3cc2612e4f779ec48fd0839c8a`; +geometry hash is +`87cfacdf639f8c7203ffdfb7da2c9e7ba60a63ed302fc6dbb5db2aed2b9047e3`. +It has one welded component and 23 boundary edges out of approximately 39,357 +edges (0.058%, retained as a provenance diagnostic). Before normalization, the +detected coherent face patch is 79.739 x 42.497 mm, 863.296 mm2, with source +normal `(0.112969, -0.253749, -0.960650)`. After the rigid axis transform the +normal is `(1, 0, 0)` to numerical precision; no dimension scaling is applied. + +### Pre-outcome iron-source substitution + +The originally selected Sketchfab iron (`dc748ddd268c4acab25c54c4048b3912`) +failed the deliberately strict identity validator before any F1 outcome ran. Its +uploader display name changed from the pinned ASCII `real_slimshady` to +`β„œπ”’π”žπ”© 𝔖𝔩𝔦π”ͺ 𝔖π”₯π”žπ”‘π”Ά`. The validator was not relaxed or Unicode-normalized. The +maintainer substituted the higher-resolution local 690CB source above. This is +an acquisition/provenance amendment, not a grid, solver, criterion, or gate +change. + +## Rejected candidates + +- GrabCAD's [library-use guidance](https://help.grabcad.com/article/246-how-can-models-be-used-and-shared) + permits public non-commercial rendering with attribution but does not clearly + grant redistribution of a raw CAD file in an AGPL repository. The maintainer's + Titleist 7-iron is therefore accepted only as an uncommitted local input. +- CGTrader's [Terms and Conditions](https://www.cgtrader.com/pages/terms-and-conditions) + prohibit making a purchased product available as a separate file. Its + Royalty Free License permits an incorporated product, not raw-mesh + redistribution. The candidate `Golf Club 7 Iron` therefore was not purchased, + downloaded, or committed. + +## Required attribution in generated results + +Every F1 result bundle records the model page, attribution, model ID, license or +use boundary, source SHA-256, and normalized mesh SHA-256. Any redistributed +driver render or derived dataset must retain the CC BY credit and indicate that +the geometry was normalized into OpenFlight's calibrated club-local coordinate +frame. The local-use-only iron mesh and normalized asset must never be +redistributed. diff --git a/src/openflight/camera/clubpose/meshes/assets/.gitignore b/src/openflight/camera/clubpose/meshes/assets/.gitignore new file mode 100644 index 00000000..d6b7ef32 --- /dev/null +++ b/src/openflight/camera/clubpose/meshes/assets/.gitignore @@ -0,0 +1,2 @@ +* +!.gitignore diff --git a/src/openflight/camera/clubpose/projection.py b/src/openflight/camera/clubpose/projection.py new file mode 100644 index 00000000..c22ad977 --- /dev/null +++ b/src/openflight/camera/clubpose/projection.py @@ -0,0 +1,330 @@ +"""Frozen Phase 1b club-state solver, promoted to the Phase 3 fusion package.""" + +from __future__ import annotations + +import math +from dataclasses import dataclass + +import cv2 +import numpy as np + +NOMINAL_RANGE_MM = 1_575.0 +CAMERA_HEIGHT_MM = 209.55 +RADAR_STATIC_BIAS_MM = 66.0069821 +BALL_RADIUS_MM = 42.67 / 2.0 +FRAME_TO_IMPACT_S = 1.0e-3 +FRAME_PERIOD_S = 2.137e-3 +MAX_EXTRAPOLATION_S = 2.5e-3 +CENTROID_NOISE_PX = 0.5 +MOMENT_EDGE_NOISE_PX = 0.5 +RANGE_NOISE_MM = 3.0 +FIT_RESIDUAL_LIMIT_PX = 8.0 +AMBIGUITY_RATIO_MIN = 1.10 +MODEL_VERSION = "phase1b-v1-frozen" +RADAR_HEIGHT_MM = 152.4 + +_CAMERA_X_MM = math.sqrt(NOMINAL_RANGE_MM**2 - CAMERA_HEIGHT_MM**2) +CAMERA_CENTER_WORLD = np.array([-_CAMERA_X_MM, 0.0, CAMERA_HEIGHT_MM]) +TARGET_WORLD = np.zeros(3) +WORLD_RIGHT = np.array([0.0, 1.0, 0.0]) +WORLD_UP = np.array([0.0, 0.0, 1.0]) +FACE_NORMAL = np.array([1.0, 0.0, 0.0]) +_FORWARD = (TARGET_WORLD - CAMERA_CENTER_WORLD) / np.linalg.norm(TARGET_WORLD - CAMERA_CENTER_WORLD) +_DOWN = np.cross(WORLD_RIGHT, _FORWARD) +_DOWN /= np.linalg.norm(_DOWN) +_R_WC = np.stack([WORLD_RIGHT, _DOWN, _FORWARD]) +RADAR_CENTER_WORLD = np.array( + [ + -math.sqrt(NOMINAL_RANGE_MM**2 - RADAR_HEIGHT_MM**2), + 0.0, + RADAR_HEIGHT_MM, + ] +) + + +@dataclass(frozen=True) +class CameraPreset: + """One explicit camera/crop configuration from approved spec section 5.""" + + name: str + width: int + height: int + fx: float + fy: float + cx: float + cy: float + plate_scale_px_per_mm: float + sensor_crop: tuple[int, int, int, int] + sampling_increment: tuple[int, int] + isp_offset: tuple[int, int] + orientation: str + gate_b1_passed: bool + physical_status: str + + @property + def horizontal_fov_deg(self) -> float: + return math.degrees(2.0 * math.atan(self.width / (2.0 * self.fx))) + + +@dataclass(frozen=True) +class ClubTemplate: + """Named analytic rear-view silhouette and speed distribution.""" + + name: str + radius_u_mm: float + radius_v_mm: float + speed_mean_mm_s: float + speed_sd_mm_s: float + impact_u_limit_mm: float + impact_v_limit_mm: float + velocity_direction: tuple[float, float, float] + + +@dataclass(frozen=True) +class SilhouetteObservation: + centroid_uv: np.ndarray + covariance_px2: np.ndarray + + +@dataclass(frozen=True) +class ClubState: + ok: bool + reason: str | None + frame_center_world: np.ndarray | None + roll_rad: float | None + fit_residual_px: float | None + calibrated_range_mm: float | None + predicted_covariance_px2: np.ndarray | None + + +def camera_presets() -> dict[str, CameraPreset]: + """Return independent intrinsics; no fixed-FOV scaling is permitted.""" + return { + "A0": CameraPreset( + name="A0", + width=320, + height=200, + fx=1033.0, + fy=1033.0, + cx=160.0, + cy=100.0, + plate_scale_px_per_mm=0.656, + sensor_crop=(336, 150, 816, 516), + sampling_increment=(2, 2), + isp_offset=(4, 4), + orientation="landscape_register_window", + gate_b1_passed=False, + physical_status="existing_320x200_plus_10us_strobe", + ), + "A1": CameraPreset( + name="A1", + width=320, + height=200, + fx=2063.0, + fy=2063.0, + cx=160.0, + cy=100.0, + plate_scale_px_per_mm=1.31, + sensor_crop=(480, 150, 320, 200), + sampling_increment=(1, 1), + isp_offset=(0, 0), + orientation="landscape_crop_metadata_sensitivity", + gate_b1_passed=False, + physical_status="plate_scale_sensitivity_only", + ), + "B": CameraPreset( + name="B", + width=1280, + height=200, + fx=2095.0, + fy=2095.0, + cx=640.0, + cy=100.0, + plate_scale_px_per_mm=1.33, + sensor_crop=(0, 300, 1280, 200), + sampling_increment=(1, 1), + isp_offset=(0, 0), + orientation="portrait_experimental", + gate_b1_passed=False, + physical_status="experimental_gate_b1_not_run", + ), + } + + +def _project(points_world: np.ndarray, camera: CameraPreset) -> tuple[np.ndarray, np.ndarray]: + points = np.asarray(points_world, dtype=float).reshape(-1, 3) + cam = (points - CAMERA_CENTER_WORLD) @ _R_WC.T + in_front = cam[:, 2] > 1e-9 + safe_z = np.where(in_front, cam[:, 2], 1.0) + uv = np.column_stack( + [ + camera.fx * cam[:, 0] / safe_z + camera.cx, + camera.fy * cam[:, 1] / safe_z + camera.cy, + ] + ) + return uv, in_front + + +def _ray_world(uv: np.ndarray, camera: CameraPreset) -> np.ndarray: + xy = np.array([(uv[0] - camera.cx) / camera.fx, (uv[1] - camera.cy) / camera.fy, 1.0]) + ray = xy @ _R_WC + return ray / np.linalg.norm(ray) + + +def _backproject_range( + uv: np.ndarray, + range_mm: float, + camera: CameraPreset, + range_origin_world: np.ndarray = CAMERA_CENTER_WORLD, +) -> np.ndarray: + """Intersect a camera ray with a range sphere around the supplied sensor.""" + ray = _ray_world(uv, camera) + offset = CAMERA_CENTER_WORLD - np.asarray(range_origin_world, dtype=float) + projection = float(offset @ ray) + discriminant = projection**2 - float(offset @ offset) + float(range_mm) ** 2 + if discriminant < 0.0: + return np.full(3, np.nan) + distance = -projection + math.sqrt(discriminant) + return CAMERA_CENTER_WORLD + ray * distance + + +def _range_mm( + point_world: np.ndarray, range_origin_world: np.ndarray = CAMERA_CENTER_WORLD +) -> float: + return float(np.linalg.norm(np.asarray(point_world) - range_origin_world)) + + +def _face_axes(roll_rad: float) -> tuple[np.ndarray, np.ndarray]: + c = math.cos(float(roll_rad)) + s = math.sin(float(roll_rad)) + return c * WORLD_RIGHT + s * WORLD_UP, -s * WORLD_RIGHT + c * WORLD_UP + + +def _velocity(template: ClubTemplate, speed_mm_s: float, reverse: bool = False) -> np.ndarray: + direction = np.asarray(template.velocity_direction, dtype=float) + direction /= np.linalg.norm(direction) + return direction * float(speed_mm_s) * (-1.0 if reverse else 1.0) + + +def _projected_velocity( + center_world: np.ndarray, velocity_world: np.ndarray, camera: CameraPreset +) -> np.ndarray: + dt = 1.0e-5 + uv_pair, front = _project( + np.stack([center_world - velocity_world * dt / 2, center_world + velocity_world * dt / 2]), + camera, + ) + if not bool(np.all(front)): + return np.zeros(2) + return (uv_pair[1] - uv_pair[0]) / dt + + +def _projection_jacobian(center_world: np.ndarray, camera: CameraPreset) -> np.ndarray: + """Pixel derivative for one millimetre along world-right/world-up.""" + epsilon = 1.0e-3 + points = np.stack( + [ + center_world - WORLD_RIGHT * epsilon, + center_world + WORLD_RIGHT * epsilon, + center_world - WORLD_UP * epsilon, + center_world + WORLD_UP * epsilon, + ] + ) + uv, front = _project(points, camera) + if not bool(np.all(front)): + return np.full((2, 2), np.nan) + return np.column_stack([(uv[1] - uv[0]) / (2.0 * epsilon), (uv[3] - uv[2]) / (2.0 * epsilon)]) + + +def _silhouette_moments( + center_world: np.ndarray, + roll_rad: float, + velocity_world: np.ndarray, + exposure_us: float, + camera: CameraPreset, + template: ClubTemplate, +) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray]: + center_uv, front = _project(center_world[None, :], camera) + if not bool(front[0]): + nan = np.full(2, np.nan) + return nan, np.full((2, 2), np.nan), nan, nan, nan + center_uv = center_uv[0] + axis_u, axis_v = _face_axes(roll_rad) + jacobian = _projection_jacobian(center_world, camera) + if not bool(np.all(np.isfinite(jacobian))): + nan = np.full(2, np.nan) + return nan, np.full((2, 2), np.nan), nan, nan, nan + body_u = np.array([float(axis_u @ WORLD_RIGHT), float(axis_u @ WORLD_UP)]) + body_v = np.array([float(axis_v @ WORLD_RIGHT), float(axis_v @ WORLD_UP)]) + vector_u = jacobian @ body_u * template.radius_u_mm + vector_v = jacobian @ body_v * template.radius_v_mm + blur_vector = _projected_velocity(center_world, velocity_world, camera) * ( + float(exposure_us) * 1e-6 + ) + covariance = ( + np.outer(vector_u, vector_u) / 4.0 + + np.outer(vector_v, vector_v) / 4.0 + + np.outer(blur_vector, blur_vector) / 12.0 + ) + extents = np.sqrt(vector_u**2 + vector_v**2) + np.abs(blur_vector) / 2.0 + return center_uv, covariance, extents, vector_u, vector_v + + +def _visible(center_uv: np.ndarray, extents: np.ndarray, camera: CameraPreset) -> bool: + if not bool(np.all(np.isfinite(center_uv))) or not bool(np.all(np.isfinite(extents))): + return False + return bool( + center_uv[0] - extents[0] >= 0.0 + and center_uv[0] + extents[0] < camera.width + and center_uv[1] - extents[1] >= 0.0 + and center_uv[1] + extents[1] < camera.height + ) + + +def _ball_geometry( + ball_center_world: np.ndarray, camera: CameraPreset +) -> tuple[np.ndarray, np.ndarray]: + center_uv, front = _project(ball_center_world[None, :], camera) + if not bool(front[0]): + return np.full(2, np.nan), np.full(2, np.nan) + center_uv = center_uv[0] + endpoints, endpoint_front = _project( + np.stack( + [ + ball_center_world + WORLD_RIGHT * BALL_RADIUS_MM, + ball_center_world + WORLD_UP * BALL_RADIUS_MM, + ] + ), + camera, + ) + if not bool(np.all(endpoint_front)): + return np.full(2, np.nan), np.full(2, np.nan) + extents = np.abs(endpoints - center_uv).max(axis=0) + return center_uv, extents + + +def _silhouette_polygon( + center_uv: np.ndarray, vector_u: np.ndarray, vector_v: np.ndarray, blur_vector: np.ndarray +) -> np.ndarray: + theta = np.linspace(0.0, 2.0 * np.pi, 24, endpoint=False) + ellipse = ( + center_uv[None, :] + + np.cos(theta)[:, None] * vector_u[None, :] + + np.sin(theta)[:, None] * vector_v[None, :] + ) + points = np.vstack([ellipse - blur_vector / 2.0, ellipse + blur_vector / 2.0]) + return cv2.convexHull(points.astype(np.float32)).reshape(-1, 2) + + +def _polygon_iou(a: np.ndarray, b: np.ndarray) -> float: + area_a = abs(float(cv2.contourArea(a))) + area_b = abs(float(cv2.contourArea(b))) + if area_a <= 0.0 or area_b <= 0.0: + return 0.0 + intersection, _ = cv2.intersectConvexConvex(a.astype(np.float32), b.astype(np.float32)) + union = area_a + area_b - float(intersection) + return float(intersection / union) if union > 0.0 else 0.0 + + +def _normalize_roll(angle: float) -> float: + return (float(angle) + math.pi / 2.0) % math.pi - math.pi / 2.0 diff --git a/src/openflight/camera/clubpose/teed_ball.py b/src/openflight/camera/clubpose/teed_ball.py new file mode 100644 index 00000000..396bc6fb --- /dev/null +++ b/src/openflight/camera/clubpose/teed_ball.py @@ -0,0 +1,321 @@ +"""Polarity-agnostic ball detection with known-radius arc fitting. + +Replaces the fixed-polarity, fixed-magnitude rule used by the superseded +synthetic pipeline (on the fork's feat/silhouette-poc branch): + + ball_mask = frame >= percentile(frame, 10) + 210.0 + +That rule fails on real captures in both directions at once. Measured on the real +OpenFlight archive (`frames.npz`, 2026-08-24): + + at address ball 192 DN on a clipped 255 DN mat -> -62 DN WRONG SIGN + in flight ball 224 DN on a 102 DN dark wall -> +122 DN right sign, + still under +210 + +No constant fixes it, because the background swings 102..255 DN inside one shot. + +What this does instead: + 1. estimate the LOCAL background, so a clipped mat and a dark wall are both handled + 2. threshold the SIGNED residual in BOTH directions, scaled to measured noise + rather than an absolute DN offset + 3. confirm candidates by SHAPE, since the ball is the one object in frame with a + known fixed physical diameter + 4. fit a circle to the reliable part of the boundary, so a ball whose lit side has + blended into the background still yields an honest radius + +Rejections are named, matching the project's fail-closed convention. +""" + +from __future__ import annotations + +from dataclasses import dataclass, field + +import cv2 +import numpy as np + +SENSOR_NOISE_DN = 5.6 # measured: 50 static frames of the real archive at gain 15.94 +SATURATION_DN = 253 # within noise of the 8-bit ceiling +BALL_DIAMETER_MM = 42.67 + + +@dataclass +class BallFit: + center: np.ndarray # (x, y) subpixel + radius_px: float + polarity: str # "dark_on_light" | "light_on_dark" + contrast_dn: float # signed, ball minus local background + arc_coverage_deg: float # angular span of boundary used for the fit + visible_area_px: int + circularity: float + mask: np.ndarray = field(repr=False) + + @property + def px_per_mm(self) -> float: + return 2.0 * self.radius_px / BALL_DIAMETER_MM + + +class BallNotFound(ValueError): + """Named rejection, so a failure says which gate it failed.""" + + +def _local_background(frame: np.ndarray, ball_px: int) -> np.ndarray: + """Median over a window several ball-widths across. + + A ball is small relative to the window, so it is smoothed away while real scene + structure (the mat, the wall, the boundary between them) survives. + """ + k = min(int(max(15, ball_px * 3)), 99) | 1 + return cv2.medianBlur(frame, k) + + +def _robust_sigma(residual: np.ndarray) -> float: + """MAD-based noise estimate, floored at the sensor's measured read noise.""" + mad = float(np.median(np.abs(residual - np.median(residual)))) + return max(1.4826 * mad, SENSOR_NOISE_DN) + + +def _fit_circle(points: np.ndarray) -> tuple[np.ndarray, float]: + """Circle fit: algebraic (Kasa) seed, then geometric Gauss-Newton refinement. + + The refinement minimises true perpendicular distance, which is the statistically + correct objective and β€” unlike a purely algebraic fit β€” stays unbiased on the + short arcs this detector actually sees when the ball's lit side has blended into + a saturated background. + """ + x, y = points[:, 0].astype(float), points[:, 1].astype(float) + if len(x) < 4: + raise BallNotFound("too_few_boundary_points") + # Kasa seed: x^2 + y^2 = 2a.x + 2b.y + c, linear in (a, b, c) + design = np.column_stack([x, y, np.ones_like(x)]) + try: + sol, *_ = np.linalg.lstsq(design, x * x + y * y, rcond=None) + except np.linalg.LinAlgError: + raise BallNotFound("circle_fit_degenerate") from None + cx, cy = sol[0] / 2.0, sol[1] / 2.0 + seed = sol[2] + cx * cx + cy * cy + if not np.isfinite(seed) or seed <= 0.0: + raise BallNotFound("circle_fit_degenerate") + r = float(np.sqrt(seed)) + + for _ in range(50): # geometric refinement + dx, dy = x - cx, y - cy + dist = np.hypot(dx, dy) + if np.any(dist < 1e-9): + break + residual = dist - r + jac = np.column_stack([-dx / dist, -dy / dist, -np.ones_like(dist)]) + try: + step, *_ = np.linalg.lstsq(jac, -residual, rcond=None) + except np.linalg.LinAlgError: + break + cx, cy, r = cx + step[0], cy + step[1], r + step[2] + if not np.isfinite([cx, cy, r]).all() or r <= 0.0: + raise BallNotFound("circle_fit_degenerate") + if np.linalg.norm(step) < 1e-9: + break + return np.array([cx, cy]), float(r) + + +def _reliable_boundary(mask: np.ndarray, frame: np.ndarray) -> np.ndarray: + """Boundary pixels whose edge is trustworthy. + + Discards boundary where the BALL'S OWN pixel is saturated: there its lit surface + has merged into a clipped background, so the apparent edge marks where clipping + began, not where the ball ends, and including it biases the radius low. + + It deliberately does NOT discard boundary merely because the background beside it + is saturated. A dark ball against a blown-out mat has a real, sharp edge all the + way round; rejecting it there threw away 17 of 28 boundary points and left a + 45-degree arc that no circle fit can use. + """ + contours, _ = cv2.findContours(mask, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE) + if not contours: + raise BallNotFound("no_boundary") + pts = max(contours, key=cv2.contourArea).reshape(-1, 2) + h, w = mask.shape + keep = [ + (px, py) + for px, py in pts + if 0 < px < w - 1 and 0 < py < h - 1 and frame[py, px] < SATURATION_DN + ] + return np.array(keep if len(keep) >= 8 else pts, dtype=float) + + +def _arc_coverage(points: np.ndarray, center: np.ndarray) -> float: + """Angular span covered by the fitted boundary, in degrees. + + A radius fitted to a 90 degree arc is far less trustworthy than one fitted to + 300 degrees, and the caller should be able to tell the difference. + """ + ang = np.degrees(np.arctan2(points[:, 1] - center[1], points[:, 0] - center[0])) + occupied = np.zeros(72, dtype=bool) # 5 degree bins + occupied[((ang + 180.0) / 5.0).astype(int) % 72] = True + return float(occupied.sum() * 5) + + +def candidates( + frame: np.ndarray, + *, + expected_radius_px: float | None = None, + radius_tolerance: float = 0.6, + sigma_k: float = 3.5, + min_area_px: int = 20, +) -> list[BallFit]: + """Every region that could be a ball, whether darker or brighter than background.""" + if frame.ndim != 2: + raise BallNotFound("frame_not_2d") + f = frame.astype(np.float32) + nominal = int(expected_radius_px * 2) if expected_radius_px else 14 + background = _local_background(frame, nominal).astype(np.float32) + residual = f - background + sigma = _robust_sigma(residual) + + found: list[BallFit] = [] + for polarity, signed in (("light_on_dark", residual), ("dark_on_light", -residual)): + mask = (signed > sigma_k * sigma).astype(np.uint8) + mask = cv2.morphologyEx(mask, cv2.MORPH_OPEN, np.ones((2, 2), np.uint8)) + mask = cv2.morphologyEx(mask, cv2.MORPH_CLOSE, np.ones((3, 3), np.uint8)) + count, labels, stats, _ = cv2.connectedComponentsWithStats(mask, 8) + for idx in range(1, count): + w, h, area = stats[idx, 2], stats[idx, 3], stats[idx, 4] + if area < min_area_px: + continue + if max(w, h) / max(min(w, h), 1) > 3.0: # a long smear is not a ball + continue + comp = (labels == idx).astype(np.uint8) + try: + pts = _reliable_boundary(comp, frame) + if len(pts) < 8: + continue + center, radius = _fit_circle(pts) + except BallNotFound: + continue + if not np.isfinite(radius) or radius <= 1.0 or radius > 0.25 * max(frame.shape): + continue + if expected_radius_px is not None: + lo = expected_radius_px * (1.0 - radius_tolerance) + hi = expected_radius_px * (1.0 + radius_tolerance) + if not lo <= radius <= hi: + continue + circularity = float(area) / (np.pi * radius * radius) + if not 0.30 <= circularity <= 1.60: + continue + contrast = float(np.median(f[comp > 0]) - np.median(background[comp > 0])) + found.append( + BallFit( + center=center, + radius_px=float(radius), + polarity=polarity, + contrast_dn=contrast, + arc_coverage_deg=_arc_coverage(pts, center), + visible_area_px=int(area), + circularity=circularity, + mask=comp, + ) + ) + + return found + + +def fit_teed_ball(frame: np.ndarray, **kwargs) -> BallFit: + """Single-frame best guess. Shape only β€” see find_ball_at_address for the real one.""" + found = candidates(frame, **kwargs) + if not found: + raise BallNotFound("no_candidate_passed_shape_gate") + return max(found, key=lambda c: (c.arc_coverage_deg / 360.0) * min(c.circularity, 1.0)) + + +def fit_teed_ball_sequence( + frames: np.ndarray, + impact_index: int, + tee_region: tuple[int, int, int, int], + *, + stride: int = 2, + lead: int = 60, + settle: int = 2, + look_after: int = 4, + cluster_px: float = 5.0, + min_persistence: float = 0.25, + **kwargs, +) -> BallFit: + """Locate the teed ball inside a known tee region. + + `tee_region` is (x0, y0, x1, y1). This is not a shortcut: the camera is fixed to + the unit and aimed at the ball, so where the ball can be is a property of the rig, + established once at install. Searching the whole frame is a harder problem than + the system actually poses, and on the real capture it is not reliably solvable β€” + a bay contains signage, shoes and knees that fit a circle better than a ball whose + lit side has blended into the mat. + + Within the region, two physical facts do the work: + * the ball is stationary for the whole address period + * the ball is GONE once the club has sent it away + + The second is the strong one. Static clutter scores 0.00 on departure; the ball + scores 1.00. + + Ambient captures can change illumination mid-sequence (the real archive drops + 58% six frames after the trigger), so the "after" window is checked to be in the + same illumination regime and truncated if not β€” otherwise everything looks like + it departed. + """ + x0, y0, x1, y1 = tee_region + inside = lambda p: x0 <= p[0] <= x1 and y0 <= p[1] <= y1 # noqa: E731 + + lo, hi = max(1, impact_index - lead), max(2, impact_index - settle) + clusters: list[dict] = [] + for index in range(lo, hi, stride): + for cand in candidates(frames[index], **kwargs): + if not inside(cand.center): + continue + for cluster in clusters: + if np.linalg.norm(cand.center - cluster["center"]) <= cluster_px: + cluster["fits"].append(cand) + cluster["frames"].add(index) + cluster["center"] = np.mean([f.center for f in cluster["fits"]], axis=0) + break + else: + clusters.append({"center": cand.center.copy(), "fits": [cand], "frames": {index}}) + if not clusters: + raise BallNotFound("no_candidate_in_tee_region") + + # keep the post-impact window inside the same illumination regime + level = float(frames[lo:hi].reshape(hi - lo, -1).mean()) + after = [ + j + for j in range(impact_index + settle, min(impact_index + settle + look_after, len(frames))) + if abs(float(frames[j].mean()) - level) < 12.0 + ] + if not after: + raise BallNotFound("no_comparable_post_impact_frame") + + sampled = len(range(lo, hi, stride)) + scored = [] + for cluster in clusters: + persistence = len(cluster["frames"]) / max(sampled, 1) + if persistence < min_persistence: + continue + survives = sum( + any( + np.linalg.norm(c.center - cluster["center"]) <= cluster_px * 2 + for c in candidates(frames[j], **kwargs) + if inside(c.center) + ) + for j in after + ) + departed = 1.0 - survives / len(after) + if departed < 0.75: # a ball that is still there is not the ball + continue + scored.append((departed, persistence, cluster)) + if not scored: + raise BallNotFound("no_departing_candidate_in_tee_region") + + # Departure ranks above persistence. A bay can hold several balls, and a spare + # that the club sweeps across is briefly hidden β€” that scrapes past the gate as a + # partial departure. The struck ball is gone completely, so ranking on departure + # first separates them instead of leaving it to a persistence tie-break. + _, _, best = max(scored, key=lambda s: (round(s[0], 3), s[1])) + fit = max(best["fits"], key=lambda f: f.arc_coverage_deg) + fit.center = best["center"] + fit.radius_px = float(np.median([f.radius_px for f in best["fits"]])) + return fit diff --git a/tests/test_clubpose_ball_detect.py b/tests/test_clubpose_ball_detect.py new file mode 100644 index 00000000..654d6338 --- /dev/null +++ b/tests/test_clubpose_ball_detect.py @@ -0,0 +1,144 @@ +"""Properties of the polarity-agnostic ball detector, pinned by test. + +Each test here corresponds to a defect found by pointing the existing detector at +the first real OpenFlight capture (`frames.npz`, 2026-08-24). They are written so a +regression re-breaks the specific thing that was broken. +""" + +from __future__ import annotations + +import numpy as np +import pytest + +from openflight.camera.clubpose.teed_ball import ( + BallNotFound, + _fit_circle, + candidates, + fit_teed_ball_sequence, +) + + +def _disc(shape, center, radius, fg, bg, *, cap=None): + """Render a disc. `cap` clips everything above that row to the background, + reproducing a ball whose lit upper surface has blended into a saturated mat.""" + image = np.full(shape, bg, dtype=np.uint8) + ys, xs = np.mgrid[0 : shape[0], 0 : shape[1]] + inside = (xs - center[0]) ** 2 + (ys - center[1]) ** 2 <= radius**2 + if cap is not None: + inside &= ys >= cap + image[inside] = fg + return image + + +@pytest.mark.parametrize("arc_deg", [360, 300, 240, 180, 120, 90]) +def test_circle_fit_is_exact_on_partial_arcs(arc_deg): + """The previous algebraic fit was biased -44% at 120 degrees of arc. + + Partial arcs are the normal case, not the exception: at address the ball's lit + side blends into the background and only part of its outline survives. + """ + angles = np.radians(np.linspace(0.0, arc_deg, max(12, arc_deg // 6), endpoint=False)) + center, radius = np.array([127.0, 154.0]), 7.3 + points = np.column_stack( + [center[0] + radius * np.cos(angles), center[1] + radius * np.sin(angles)] + ) + fitted_center, fitted_radius = _fit_circle(points) + assert fitted_radius == pytest.approx(radius, abs=1e-6) + assert fitted_center == pytest.approx(center, abs=1e-6) + + +def test_finds_ball_darker_than_its_background(): + """The real teed ball is 62 DN DARKER than the clipped mat it sits on. + + The rule it replaces asked for `>= background + 210 DN`, which cannot express a + negative contrast at any threshold value. + """ + frame = _disc((200, 320), (127, 154), 7, fg=192, bg=255) + found = candidates(frame) + assert found, "a dark ball on a light background must be detectable" + best = min(found, key=lambda c: np.linalg.norm(c.center - np.array([127.0, 154.0]))) + assert best.polarity == "dark_on_light" + assert best.contrast_dn < 0 + assert best.center == pytest.approx(np.array([127.0, 154.0]), abs=1.5) + + +def test_finds_ball_brighter_than_its_background_below_the_old_threshold(): + """The real airborne ball is +122 DN, which the old `+210 DN` bar also failed.""" + frame = _disc((200, 320), (128, 84), 7, fg=224, bg=102) + found = candidates(frame) + best = min(found, key=lambda c: np.linalg.norm(c.center - np.array([128.0, 84.0]))) + assert best.polarity == "light_on_dark" + assert 0 < best.contrast_dn < 210 + assert best.center == pytest.approx(np.array([128.0, 84.0]), abs=1.5) + + +def test_detector_is_not_fooled_by_a_static_round_distractor(): + """A bay contains signage and shoes that fit a circle better than a blended ball. + + Shape alone picked the golfer's leg in 10 real frames out of 10. Departure is + what separates them: the ball leaves, the clutter does not. + """ + frames = [] + for index in range(40): + frame = _disc((200, 320), (40, 30), 9, fg=230, bg=100) # static distractor + if index < 30: # ball, gone after 30 + ys, xs = np.mgrid[0:200, 0:320] + frame[(xs - 127) ** 2 + (ys - 154) ** 2 <= 49] = 60 + frames.append(frame) + result = fit_teed_ball_sequence( + np.array(frames), 30, (80, 120, 200, 190), lead=30, look_after=4 + ) + assert result.center == pytest.approx(np.array([127.0, 154.0]), abs=3.0) + + +def test_rejection_is_named_not_silent(): + frames = np.full((20, 200, 320), 128, dtype=np.uint8) + with pytest.raises(BallNotFound) as excinfo: + fit_teed_ball_sequence(frames, 12, (80, 120, 200, 190), lead=10, look_after=3) + assert str(excinfo.value) in { + "no_candidate_in_tee_region", + "no_departing_candidate_in_tee_region", + "no_comparable_post_impact_frame", + } + + +def _mat_scene(n_frames, balls, occluder=None): + """A blown-out mat carrying several balls. `balls` is (x, y, departs_at|None).""" + frames = [] + for index in range(n_frames): + frame = np.full((200, 320), 255, dtype=np.uint8) + ys, xs = np.mgrid[0:200, 0:320] + for bx, by, departs in balls: + if departs is not None and index >= departs: + continue + frame[(xs - bx) ** 2 + (ys - by) ** 2 <= 49] = 192 + if occluder is not None: + ox, oy, active = occluder + if index in active: + frame[max(0, oy - 14) : oy + 14, max(0, ox - 16) : ox + 16] = 120 + frames.append(frame) + return np.array(frames) + + +TEE = (80, 120, 220, 190) + + +def test_picks_the_struck_ball_when_several_are_on_the_mat(): + """A real bay holds spares. Only the struck ball leaves.""" + scene = _mat_scene(40, [(127, 154, 30), (160, 168, None), (100, 140, None)]) + result = fit_teed_ball_sequence(scene, 30, TEE, lead=30, look_after=4) + assert result.center == pytest.approx(np.array([127.0, 154.0]), abs=4.0) + + +def test_a_spare_ball_swept_over_by_the_club_is_not_mistaken_for_the_struck_one(): + """Brief occlusion looks like a partial departure; full departure must outrank it.""" + scene = _mat_scene(40, [(127, 154, 30), (175, 170, None)], occluder=(175, 170, {32, 33, 34})) + result = fit_teed_ball_sequence(scene, 30, TEE, lead=30, look_after=4) + assert result.center == pytest.approx(np.array([127.0, 154.0]), abs=4.0) + + +def test_selection_is_by_departure_not_by_position_in_the_tee_region(): + """The struck ball need not be the one nearest the middle of the region.""" + scene = _mat_scene(40, [(200, 180, None), (95, 130, 30)]) + result = fit_teed_ball_sequence(scene, 30, TEE, lead=30, look_after=4) + assert result.center == pytest.approx(np.array([95.0, 130.0]), abs=4.0) From e18f6d01ffca8b6d930eebffa8879ebf3f8e9d4b Mon Sep 17 00:00:00 2001 From: HarjotDhanota <152345739+HarjotDhanota@users.noreply.github.com> Date: Fri, 28 Aug 2026 11:55:04 -0700 Subject: [PATCH 2/2] feat(camera): clubhead pose fit, scores and delivered angles Sequence fitting with physical bounds and temporal smoothness, boundary-distance and rotation-consistency scores (silhouette IoU measured inversely related to pose correctness on real masks), clubhead/shaft separation, and delivered loft/face/lie with plausibility envelopes. Includes the technical report: clubhead orientation has no accuracy figure against truth, and this change does not claim one. --- docs/CHANGELOG.md | 6 + docs/clubface-impact-location-report.md | 379 +++++++++++++++++ docs/clubface-impact-location.md | 96 +++++ src/openflight/camera/clubpose/angles.py | 134 ++++++ src/openflight/camera/clubpose/fit.py | 400 ++++++++++++++++++ src/openflight/camera/clubpose/head_split.py | 71 ++++ src/openflight/camera/clubpose/motion.py | 39 ++ src/openflight/camera/clubpose/scores.py | 193 +++++++++ tests/test_clubpose_club_angles.py | 98 +++++ tests/test_clubpose_fit_real_range_grid.py | 56 +++ ...test_clubpose_fit_sequence_pinned_range.py | 33 ++ tests/test_clubpose_head_split.py | 113 +++++ tests/test_clubpose_pose_scores.py | 197 +++++++++ tests/test_clubpose_rigid_motion.py | 46 ++ 14 files changed, 1861 insertions(+) create mode 100644 docs/clubface-impact-location-report.md create mode 100644 docs/clubface-impact-location.md create mode 100644 src/openflight/camera/clubpose/angles.py create mode 100644 src/openflight/camera/clubpose/fit.py create mode 100644 src/openflight/camera/clubpose/head_split.py create mode 100644 src/openflight/camera/clubpose/motion.py create mode 100644 src/openflight/camera/clubpose/scores.py create mode 100644 tests/test_clubpose_club_angles.py create mode 100644 tests/test_clubpose_fit_real_range_grid.py create mode 100644 tests/test_clubpose_fit_sequence_pinned_range.py create mode 100644 tests/test_clubpose_head_split.py create mode 100644 tests/test_clubpose_pose_scores.py create mode 100644 tests/test_clubpose_rigid_motion.py diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index fd1a42f4..92d3f6a2 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -13,6 +13,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 projection, and sub-pixel location of the teed ball, which anchors the world frame and the camera range. First stage of the clubface impact-location work; not yet wired into the shot pipeline. +- **Clubhead-pose fitting and delivered angles.** Sequence pose fit with + physical bounds, boundary-distance and rotation-consistency scores, clubhead + and shaft separation, and delivered loft/face/lie with plausibility + envelopes. Includes the technical report on what is validated (fused + radar+camera clubhead velocity, ball detection, impact timing) and what is + not (clubhead orientation has no accuracy figure against truth). ### Added - **Automatic OV9281 exposure control.** High-speed camera capture now measures diff --git a/docs/clubface-impact-location-report.md b/docs/clubface-impact-location-report.md new file mode 100644 index 00000000..7637cdaf --- /dev/null +++ b/docs/clubface-impact-location-report.md @@ -0,0 +1,379 @@ + + +# Markerless clubface impact location + +> This is the markdown rendering of the [technical report](https://claude.ai/code/artifact/c8817c34-c3ea-4455-9700-cf5a4e238b75). The web version carries 11 figures β€” real frames with the model's own projections overlaid β€” which are not reproduced here. + + +*OpenFlight Β· technical report Β· session 20260825_181734, 21 shots* + +An assessment of whether clubface impact location and face angle can be measured from a single behind-ball camera under ambient light, with no markers on the ball or club. This report states what has been measured, what has not, and what would resolve the remaining questions. + +## 01. Summary + +Ball flight measurement is working and validated against the radar. Club measurement is not, and the reason is now understood well enough to act on. + +> **βœ… Established** +> **Ball detection** succeeds on 21 of 22 captures with no false detections. **Impact timing** from camera and radar agree to **0.66 frames**. **Camera attitude** is measured rather than assumed: boresight pitch **βˆ’0.185Β° Β± 0.111Β°**, with the ball centre **163 mm** below the lens. The **7-iron/9-iron launch angle difference is real** β€” an independent reconstruction gives +2.91Β° against the radar estimator's +2.59Β°, and it survives a 13Β° change in assumed camera pitch. Camera and radar independently agree on the tee position to **1632 Β± 53 mm** against a taped **1581 mm**. + +> **⚠️ Not established** +> **Face angle, dynamic loft and impact location** are model-dependent inferences, not measurements. **No accuracy figure exists for any of them**, because the comparison harness does not yet cover club metrics. Radar club path and attack angle are **rejected on 21 of 21 shots**. Camera and radar disagree by approximately **5Β° in both axes**, and no current data arbitrates between them. + +> **⚑ Three findings that determine what to do next** +> **1. The pose fit is noise-limited, not model-limited.** One pixel of segmentation error is worth approximately **ten degrees** of face angle. The first 5Β° of face angle change the projected silhouette by **zero pixels**. +> **2. The clubhead range is modelled incorrectly.** Every fit renders the club at the range of the *ball*. The radar shows the club traversing **529 mm** of range during the frames being fitted. Constant range is rejected at **p β‰ˆ 0.04**. +> **3. Changing the fit metric does not help.** Two objectives with independent failure modes recover poses differing by a median of **12.7Β°** β€” the same magnitude as the segmentation noise. + +> **⚠️ Correction affecting several figures in this report** +> The acoustic trigger lag was previously stated as **6.0 frames**. It is **2.11 frames** on this rig β€” the sound's travel time over 1.575 m β€” and it varies with how far the unit sits from the ball. The wrong constant put contact at frame 68 on every shot instead of a per-shot value near 71.9, and it anchored the clubhead range model and every pose fit below. Those results are **superseded**; see Β§03. Production measurement paths are unaffected. + +The highest-value next step is not another estimator. The club-metric comparison harness has now been **built and shipped**; the step is a **session alongside a Trackman** to give it truth to compare against. Full priorities, including what was completed on 2026-08-27 and what each completion showed, in [Recommendations](#recommendations). + +## 02. Measurement setup + +A single monochrome global-shutter camera behind the ball, with two radars. All figures in this report come from one session of 21 correctly exposed shots (7-iron and 9-iron), captured 2026-08-25. + +| Parameter | Value | Source | +|---|---|---| +| Sensor / mode | OV9281, 320Γ—200 | 2Γ— subsampled readout | +| Frame rate | 467.6 fps | measured, 0 dropped frames | +| Exposure | 247–298 Β΅s | measured | +| Lens / focal length | 2.8 mm, fx = 466.7 px | datasheet optics, not a calibrated matrix | +| Plate scale at the ball | 0.2952 px/mm | derived; 1 px = 3.39 mm | +| Camera height | 203.2 mm | tape | +| Camera-to-ball range | 1581 mm | tape chain; radar agrees to 1632 Β± 53 mm | +| Boresight pitch | βˆ’0.185Β° Β± 0.111Β° | recovered from footage, 21 shots | +| Radars | OPS243-A + IWR6843 | 24 GHz Doppler; 62 GHz FMCW | + +**Two properties of this configuration constrain everything that follows.** The plate scale means the clubhead spans roughly 32 pixels. And there is no distortion model, no independently estimated principal point, and no separate fx/fy β€” the intrinsics are nominal, derived from the datasheet lens over the effective pixel pitch. + +> **Known gap in the setup record** +> The enclosure used for this session is not dimensioned in any drawing. Camera pitch was recovered from the footage rather than specified, so it is not reproducible on a second unit. Establishing the mounting geometry, and ideally deriving camera attitude from the existing LIS3DH inclinometer as the radar already does, would remove that dependency. + +## 03. Validated measurements + +These measurements are supported by cross-sensor agreement or by an independent reconstruction, and are the parts of the system suitable to build on. + +### Ball detection + +Detection succeeds on 21 of 22 captures with no false positives. The one failure is a capture taken at 495 Β΅s and gain 15, which saturated 99.8 % of the frame; it is excluded from all analysis in this report. + +Detection cannot rely on brightness alone. At address the ball sits against a mat driven past the sensor's ceiling and registers as a *dark* object; in flight it is the brightest thing in frame. The detector must accommodate both polarities. + +> **Figure 1** β€” The same ball, 120 ms apart, at 9Γ—. At address it sits on a mat driven past the sensor's ceiling, so it registers as a **dark** object β€” 192 DN against 255. Once airborne against the dark backdrop it is a **bright** object at +141 DN. The contrast **inverts sign within a single shot**. +> *(image in the [web version](https://claude.ai/code/artifact/c8817c34-c3ea-4455-9700-cf5a4e238b75))* + +### Ball geometry at the tee + +The teed ball's image is not circular. It is measurably flattened across the top, which biased an earlier radius estimate and, through it, the assumed camera-to-ball range. Fitting the boundary rather than thresholding the bright region removes the bias. + +> **Figure 2** β€” The teed ball at 24Γ—, four different shots. **Olive** is the detected mask, **violet** the fitted circle. The bright ball region is visibly flattened across the top, and the circle's upper arc passes through grey mat rather than white ball. The blob is wider than it is tall β€” which a sphere cannot be. +> *(image in the [web version](https://claude.ai/code/artifact/c8817c34-c3ea-4455-9700-cf5a4e238b75))* + +### Impact timing + +Contact precedes the acoustic trigger by the time the sound takes to reach the unit: + +``` +impact_time = trigger_time - distance_to_microphone / speed_of_sound(T) +``` + +On this rig the ball sits **1.575 m** from the unit, giving **4.59 ms** or **2.11 frames** at 468 fps. Two independent routes agree: + +| Method | Impact frame | +|---|---| +| Acoustic model β€” distance Γ· speed of sound | 71.85 | +| Ball departure, measured per shot (n=20) | 71.89 Β± 0.77 | + +They agree to **0.04 frames**. The SEN-14262 hardware path is about 10 Β΅s and is negligible; there is no unexplained detector latency. + +> **⚠️ This corrects a figure published earlier in this report** +> Earlier versions stated a **6.0-frame** lag and warned readers off the ball-track estimate. **Both were wrong.** The 6.0 constant put contact at frame 68 on every shot β€” off by **3.89 frames**, about 8 ms β€” and the ball track, which was dismissed, was correct. +> The error came from misreading a render. The clubhead is **27 px wide** and sits adjacent to the ball for two or three frames before it strikes, so "the head reaches the ball" is not contact. Contact is when the **ball starts moving**. Per-shot values run **70.65 to 73.64**, not a constant. +> **Everything anchored on that constant is superseded** β€” the clubhead range model in Β§04, the pose fits, and a claimed "post-impact frames were fitted" defect that was an artefact of the wrong anchor. + +> **⚑ The lag belongs to the installation, not to the software** +> It scales with how far the unit sits from the ball, so a fixed frame offset is only ever right for the rig it was measured on. At 468 fps, with the trigger at frame 74: +> +> | Ball to unit | Lag | Impact frame | +> |---|---|---| +> | 1.0 m | 2.91 ms | 72.64 | +> | **1.575 m β€” this rig** | 4.59 ms | 71.85 | +> | 2.5 m | 7.28 ms | 70.59 | +> | 3.5 m | 10.20 ms | 69.23 | + +Distance dominates: doubling it doubles the lag, while the whole 0–40 Β°C range moves the speed of sound about 7 %. Shipped as `src/openflight/acoustic.py` with 21 tests. `tee_range_m` in `iwr6843/calibration.py` is the distance source. + +> **βœ… Production already solved this; the research code did not use it** +> `iwr6843/shot.py:impact_time_s` back-extrapolates the *ball's own range walk* to the tee rather than trusting the trigger, and its docstring records why: assuming the trigger's ring position "is why the club-path estimator was fitting the follow-through". `camera/club_delivery.py` likewise *detects* the impact frame and uses the trigger only as a Β±8/+10 frame plausibility gate. +> **So the production measurement paths are not affected.** The defect was confined to the research scripts, which reinvented a solved problem and got it wrong. An earlier draft of this panel claimed a production-wide 4.6 ms bias; that claim was itself incorrect and is withdrawn. + +### Launch angle + +The 7-iron/9-iron launch angle difference is real and not an artefact of the estimator. An independent reconstruction from camera rays plus radar range walk gives **+2.91Β°** and **+4.22Β°** against the shipped estimator's +2.59Β° and +3.60Β°, and the result is invariant to a 13Β° change in assumed camera pitch. A separate hypothesis, that the estimator carries a club-dependent bias, was tested and refuted (βˆ’0.003Β°). + +## 04. Clubhead pose estimation + +The clubhead is located reliably. Its orientation is not, and the limiting factor is not the fitting method. + +### The club model and its reference frame + +The model is a triangle mesh of a Titleist 690CB 7-iron. Pose is solved as a 3D centre plus orientation, projected through the camera model and rasterised; every overlay in this report is the model's own projection, unpadded. + +> **Figure 3** β€” **Left** renders the surfaces facing `+x` β€” the cavity back, with its recessed oval and perimeter rim. **Right** renders `βˆ’x` β€” the striking face, flat with fine parallel grooves. Both are z-buffered per pixel. A plain silhouette render cannot tell front from back at all: it is the same outline, mirrored. +> *(image in the [web version](https://claude.ai/code/artifact/c8817c34-c3ea-4455-9700-cf5a4e238b75))* + +> **⚠️ Defect: the model's reference frame is anchored to the back of the club** +> `detect_face_plane` selects the plane by an extremity criterion. On a cavity-back iron the hosel protrudes past the striking face, so the criterion selects the **cavity rim on the reverse side**. Measured loft was reported as 17.5Β°; the true values are **33.10Β° loft and 61.19Β° lie**. +> **Every dynamic-loft figure derived from the mesh inherits this.** The two images below are the two candidate surfaces; the scorelines identify the correct one unambiguously. + +> **Figure 4** β€” The surface the code anchors the model’s frame to, seen face-on. It is the **back** of a cavity-back iron β€” a perimeter rim (orange) around a recessed cavity floor (blue), with the sole in green. Every pixel here is the model’s own geometry, coloured by which candidate surface each triangle belongs to. +> *(image in the [web version](https://claude.ai/code/artifact/c8817c34-c3ea-4455-9700-cf5a4e238b75))* + +> **Figure 5** β€” The other side, which both earlier passes had dismissed as β€œthe recessed cavity floor”. The scorelines settle it. **This is the striking face.** +> *(image in the [web version](https://claude.ai/code/artifact/c8817c34-c3ea-4455-9700-cf5a4e238b75))* + +### Fitting the model to real frames + +Position and scale behave. Orientation does not. The fit locates a clubhead-shaped object in approximately the right place on every pre-impact frame, and the printed angles are the part that should not yet be relied on. + +> **Figure 6** β€” Every second frame from F58 to F80 of shot 014, raw sensor pixels at 3Γ— nearest-neighbour. The clubhead now has a readable outline β€” sole, topline, hosel, face β€” for the whole pass, and the shaft is a bright unbroken line. On the old capture the impact zone was 83–94 % clipped and none of this existed. +> *(image in the [web version](https://claude.ai/code/artifact/c8817c34-c3ea-4455-9700-cf5a4e238b75))* + +> **Figure 7** β€” Shot 014, F63–F77. **Olive** is the observed silhouette the fit consumes. **Pink** is the projected 3D mesh at its fitted pose β€” the model’s own output, unpadded. The olive outline hugs the head in all 15 frames. The pink one wanders: on adjacent frames the fit reports pitch βˆ’40Β° then +70Β°, roll 0Β° then 135Β°, range 1340 mm then 1645 mm. Teal marks the teed ball until it departs. +> *(image in the [web version](https://claude.ai/code/artifact/c8817c34-c3ea-4455-9700-cf5a4e238b75))* + +> **Figure 8** β€” Left: real pixels with the observed silhouette (olive) and the projected mesh (violet). Middle: the same fitted pose **shaded, from the camera's own viewpoint**. Right: the same pose from a **fixed** viewpoint with world axes, so frames can be compared against each other. Both renders are the model's output β€” they show what the fit believes, which is the point. +> *(image in the [web version](https://claude.ai/code/artifact/c8817c34-c3ea-4455-9700-cf5a4e238b75))* + +### Information available in the silhouette + +The silhouette carries much less orientation information than it appears to. Measured directly on the mesh, projected width against face angle: + +| Face angle | 0Β° | 5Β° | 10Β° | 15Β° | 20Β° | 30Β° | +|---|---|---|---|---|---|---| +| Projected width | 32 px | 32 px | 30 px | 30 px | 28 px | 26 px | + +**The first five degrees of face angle change the silhouette by no pixels at all.** Ten degrees changes it by two. This is a property of the projection, not of the fitting method β€” a clubhead rotating about the vertical axis presents an almost stationary width, because the thickness rotating into view compensates for the face length rotating out of it. + +Sweeping each axis around a fitted pose on real frames and measuring how far it can move before either metric registers a change: + +| Axis | IoU (0.01 threshold) | Chamfer (0.1 px threshold) | +|---|---|---| +| Yaw β€” face angle | Β±7Β° | Β±10Β° | +| Pitch β€” dynamic loft | Β±8.5Β° | Β±14Β° | +| Roll β€” lie | Β±7Β° | Β±9.5Β° | + +### Segmentation as the binding constraint + +Perturbing the observed mask by **one pixel** β€” a dilation or erosion, which is what a change in segmentation threshold produces β€” moves the fit score by **0.40 to 1.45 times** as much as a **Β±30Β° pose error**. On four of six shots, the boundary metric moved further for one pixel of mask than across the entire 60Β° sweep. + +**One pixel of segmentation error is worth approximately ten degrees of face angle.** + +This also explains why silhouette overlap ran inversely to pose correctness in earlier work. Mask area changes by about **3 %** for 10Β° of face angle, but by roughly **25 %** for a one-pixel dilation β€” so an area-based metric is around eight times more responsive to segmentation quality than to twenty degrees of pose. It was measuring mask quality, and the arms with cleaner masks scored better while recovering worse poses. + +Consequence: the effective levers are sub-pixel edge extraction and plate scale, not the fitting algorithm. + +### The clubhead range model + +The fitter places the mesh at a fixed 1581 mm, the camera-to-ball distance. The clubhead is not at that range during the frames being fitted β€” it arrives from behind the ball and passes through it. + +Running the production club tracker over the raw radar cube places the clubhead at **1.042 m** five frames before impact and **1.571 m** at contact, a **529 mm** traverse, extrapolating to the tee with 29 mm of error. Across all 21 shots the radar's summary fields agree and independently confirm the taped ball position: + +``` +track start 1238 +- 24 mm ++ range_rate x span += track end 1632 +- 53 mm tee ball, by tape: 1581 mm +``` + +> **βœ… Test independent of orientation** +> If the club recedes, its projected area must fall as 1/rΒ². Observed clubhead mask area, last pre-impact frame divided by first, n=10 shots: +> **Observed 0.829 Β± 0.222** Β· radar-derived range model predicts **0.813** Β· constant range predicts **1.000**. +> Constant range is **rejected at p β‰ˆ 0.04**; the radar-derived model is consistent at p β‰ˆ 0.8. This test uses no orientation parameter, so refitting angles cannot account for it. + +This is also a mechanism for the metric behaviour above rather than a separate symptom. With the model systematically under-scaled on the early frames, the pose that projects largest matches best regardless of its orientation β€” and the orientation angles were the only free parameters available to absorb a scale error. + +> **Figure 9** β€” **The mismatch, at the range the tape says.** The model covers only **43–55 %** of the observed pixels. A shorter range makes the projection larger, so the fit was pulling the club nearer to close that gap. Note the thin cyan tail running up and left out of each silhouette: that is the shaft, and **the model has none** β€” its 62 mm protrusion measures out as hosel and ferrule (Β§11h). The model cannot cover that tail at any range. +> *(image in the [web version](https://claude.ai/code/artifact/c8817c34-c3ea-4455-9700-cf5a4e238b75))* + +**Correcting the render alone is not sufficient** (IoU βˆ’0.0035, chamfer +0.021 px). Orientation was fitted under the constant-range assumption and must be re-solved with the corrected range model. That work is outstanding. + +### Objective function selection + +Because overlap was suspect, a boundary-distance metric was implemented and evaluated as a replacement. It fails in a different way from overlap β€” area versus shape β€” so agreement between the two is informative. + +Refitting from identical seeds under each metric, the recovered poses differ by a **median of 12.7Β°** (n=18, range 5.8–28.0Β°). + +> **Figure 10** β€” **One shot, three frames before impact, fitted three ways.** Cyan is the observed silhouette the fit consumed; orange is the model’s own projection at the fitted pose. Rows: the shipped depth grid, the corrected grid, and range pinned at the measured 1581 mm. Reading across a row shows how coherent the pose sequence is; reading down a column shows what the depth treatment changed. The ball is the pale disc at bottom centre. +> *(image in the [web version](https://claude.ai/code/artifact/c8817c34-c3ea-4455-9700-cf5a4e238b75))* + +**When the choice of objective function moves the recovered orientation by more than ten degrees, the data is not determining the pose.** Two metrics with independent failure modes reaching the same limit is a stronger result than either alone: the fit is noise-limited, not objective-limited. + +Recommendation: no further effort on the objective function until segmentation quality or plate scale improves. + +## 05. Radar contribution + +The radar already measures quantities the pose fit does not use, and its own club-angle output is rejected on every shot for a reason worth diagnosing. + +### What the radar measures today + +Impact timing is available to approximately 33 Β΅s from the OPS243 30 kHz I/Q buffer. The IWR6843 tracks the clubhead's range and range rate through the approach. The range–time map separates the clubhead, the ball and static clutter cleanly: + +> **Figure 11** β€” Range–time map for shot 014, static returns removed. Three things separate cleanly. A **diagonal from 1.05 to 1.35 m in the 12 ms before impact** β€” the clubhead approaching. A **bright vertical band at 1.85–1.95 m that never walks in range** for the entire 72 ms β€” the golfer’s body and arms, which move but do not translate. And a second diagonal departing outward after impact β€” the ball. Time runs down; each block of twelve rows is one 3 ms frame. +> *(image in the [web version](https://claude.ai/code/artifact/c8817c34-c3ea-4455-9700-cf5a4e238b75))* + +Per shot, the radar currently reports and the pose fit currently discards: clubhead range at track start (**1238 Β± 24 mm**), range rate (approximately **33 m/s**), azimuth rate, and track span. The clubhead travels within about **25Β°** of the radar boresight, so 91 % of its speed is measured directly. + +The 22 raw capture files are included in the session export and are decoded by `src/openflight/iwr6843/dump.py`, which is production code. Prior to this report the silhouette work had never read them. + +### Club path and attack angle + +Club path and attack angle are rejected on **21 of 21 shots**, always with status `rejected_phase_span`. Two observations narrow the cause: + +| Observation | Value | Expected | +|---|---|---| +| Azimuth phase span | 2.18–3.91 rad | β‰ˆ1.3 rad; ceiling Ο€/2 | +| Attack angle, all shots | βˆ’25.3Β° to βˆ’37.3Β° (sd 2.8Β°) | β‰ˆβˆ’4Β° for a 7-iron | +| Club path, all shots | βˆ’8.6Β° to +37.1Β° | a few degrees | + +Attack angle returning a tightly clustered value near βˆ’31Β° on every shot regardless of the swing is a systematic artefact, not a measurement. The apparent azimuth swing is **three to ten times larger than the clubhead can physically produce**. + +A plausible mechanism is scatterer migration across an extended, rotating target: the clubhead subtends about 4Β° at 1.25 m while rotating at roughly 1300Β°/s, so the dominant scattering point moves between frames. **This is a hypothesis and has not been confirmed.** Note that the phase-span check deliberately does not unwrap, for documented reasons β€” unwrapping fabricated angles in earlier work. + +### Cross-range resolution and Doppler + +The radar cannot image the clubhead directly. Angular resolution is set by aperture, and the array is 19.3 mm wide: + +| Quantity | Value | +|---|---| +| Wavelength (62 GHz) | 4.835 mm | +| Range resolution | 46.9 mm | +| Aperture (8 virtual elements at Ξ»/2) | 19.3 mm | +| Beamwidth | 12.7Β° | +| **Cross-range cell at 1.25 m** | **277 mm** | +| Clubhead, for comparison | 90 mm β€” 0.33 cells | + +Clubhead, shaft and hands fall within a single angular cell. Additional pulses improve signal-to-noise and velocity resolution; they do not improve angular resolution. + +**A rotating target does synthesise an effective aperture** (inverse synthetic aperture radar). At 1300Β°/s the head rotates 15.1Β° across the 11.7 ms tracked, giving `Ξ»/(2Δθ)` = **9.2 mm** of cross-range resolution β€” a 30-fold improvement, approximately ten cells across the head. Micro-Doppler analysis gives the same cell count, as it must. + +Measured on the raw cube with range walk corrected, across **21 shots and 112 frames**: + +| Measurement | Doppler bins | +|---|---| +| Point-target floor, same estimator and window | 1.27 | +| Predicted from rotation alone (838 Hz) | 1.36 | +| Expected in quadrature | 1.86 | +| **Measured clubhead median** | **1.95** | + +The implied toe-to-heel spread is **893 Hz** against **838 Hz** predicted. The clubhead return is measurably broader than a point target, by close to the amount its rotation should produce. + +> **⚠️ This is not yet evidence of rotation** +> The discriminating test is whether the spread scales with club speed. It returned a negative correlation, but the test has **no statistical power on this data**: the predicted effect across the full speed range is **0.154 bins** against an observed scatter of **0.490 bins**, and club type is perfectly confounded with speed (7-iron 37.4–38.8 m/s, 9-iron 34.6–36.4 m/s, no overlap). Additionally, **23 % of frames return a width below the point-target floor**, which is unphysical and indicates a noisy estimator at 12 samples. +> Two measurements would resolve it: a session spanning a wide club-speed range, and a positive control on the ball, whose spin predicts roughly 8 bins of spread. Neither has been run. + +If the rotation signal is real, the useful output is not an image. **Focusing an ISAR image requires estimating the target's rotation rate and axis** β€” the two parameters the pose fit currently leaves free. + +## 06. Comparison with commercial systems + +Two commercial systems solve this problem behind the ball, both with less capable cameras than ours. + +| System | Camera | Illumination | Markerless impact location | +|---|---|---|---| +| **Trackman 4** | 720p @ 60 fps | ambient, 700–800 lux | yes | +| **Mevo Gen 2** | single phone-class module | ambient, 300 lux minimum | yes | +| **OpenFlight** | 468 fps | ambient | not yet | + +The relevant difference is architectural rather than optical. Trackman's approach fuses the camera with radar that supplies kinematics and timing at 40 kHz; at 60 fps the clubhead travels roughly 0.7 m between frames, so the camera cannot track impact independently and is not required to. Impact location is a product of the fusion, not of frame rate. + +OpenFlight has the ingredients: **8Γ— Trackman 4's frame rate**, impact timing to approximately 33 Β΅s, and radar kinematics from two devices. What is missing is the fusion model β€” and, as Β§04 shows, a camera term whose noise floor is currently around ten degrees per pixel. + +Comparator set is Trackman 4, Full Swing KIT and Mevo Gen 2. Trackman iO is excluded: it is ceiling-mounted, and its frame rate is the price of markerless spin from overhead rather than a behind-ball figure. + +## 07. Recommendations + +Ordered by dependency. Each item carries the measurement that justifies it. Items completed since this report was first issued are kept, with their outcomes, because several outcomes changed what the remaining items are worth. + +### Completed, 2026-08-27 β€” and what each one showed + +| Item | Outcome | +|---|---| +| **Extend the Trackman comparison to club metrics** | Written; in a separate PR. | +| `compare_trackman.py` pairs OpenFlight and Trackman shots but compares ball data only, so no club-side figure can be scored against truth. An extension covering twelve club-delivery metrics is written and tested, but OpenFlight can currently supply only two of them (attack angle and club path) and there is no Trackman session to run it against, so it is offered as its own change rather than bundled here. | +| **Impact timing from the installation** | Model validated; module on the fork. | +| Contact = trigger βˆ’ distance Γ· speed of sound: the model and measured ball departure agree to **0.04 frames**, correcting a 3.89-frame error that had anchored earlier fits (Β§03). The module ships when the capture path consumes it; nothing in the runtime calls it yet. | +| **Surface computed trajectory metrics** | Open; not in this branch. | +| The ballistics simulator computes apex, lateral deviation, flight time, landing speed and landing angle on every shot, and `server.py` reads only `carry_yards`. Five Trackman-parity outputs are discarded at no extra cost. Wiring them through is unrelated to camera vision, so it is left for a separate change. | +| **Radar range model + rotation-axis constraint** | Implemented; orientation still fails. | +| Range now comes from the radar's own range rate and the rotation axis from the fused velocity, dropping the fit from five free parameters to four. The velocity half **validates**: fused \|v\| matches the OPS243's independent club speed with a mean ratio of **0.970** (sd 0.029, spread 0.941–1.015 across 6 shots; worst shot 5.9% off). The OPS243 takes no part in the fit, so this is a genuine cross-sensor check. The orientation half does not: **0 of 6** shots land inside the physical envelope, and refitting with the corrected impact anchor moved the recovered angles substantially β€” a fit that sensitive to its time anchor is not extracting orientation from the pixels. **This is why the remaining items are ordered as they are.** | +| **Analyse the raw radar captures** | Opened; rotation unconfirmed. | +| All 22 `.l3dump` files decoded. The clubhead's Doppler width (1.95 bins median) exceeds the point-target floor (1.27) by close to the rotation-predicted amount, but the discriminating test has no statistical power on a 7-iron/9-iron-only session. Needs the wide-speed-range capture below. | + +### Tier 1 β€” prerequisite for any accuracy claim + +| Item | Justification | Cost | +|---|---|---| +| **A session alongside a Trackman, using the new club-metric comparison** | The harness exists; no club figure can be validated until it has truth to compare against. | 1 session | +| **Target at a taped position, visible to both sensors** | Resolves the **5Β°** camera/radar disagreement, which nothing in the current data can arbitrate. | 1 session | + +### Tier 2 β€” actionable with existing data + +| Item | Justification | Cost | +|---|---|---| +| **Extract more of the frames the club already appears in** | The club is visible for roughly **ten frames** before contact (about f62–f72 at this framing); the current extractor keeps **3–5**, losing the early frames against the dark netting. A better segmenter therefore roughly **doubles** the observations per shot β€” a bounded gain, not an open-ended one, since nothing recovers more frames than the club is in view for. Against four fit parameters, 5β†’10 observations changes the conditioning materially. | days | +| **Correct `detect_face_plane`** | Still anchors the model frame to the cavity rim on the reverse of the club (true loft **33.10Β°**, reported 17.5Β°). Interim workaround exists: `replay/club_angles.py` carries the measured axes, and its `square_pose()` is now the required seed for any fit β€” the mesh frame's origin is a *backwards* club. | hours | + +### Tier 3 β€” requires a new capture or hardware change + +| Item | Justification | Cost | +|---|---|---| +| **Capture at 1280Γ—800, 1:1** | Plate scale doubles to **0.655 px/mm**, so 10Β° of face angle becomes 4 px rather than 2, at the same frame rate and field of view. Multiplies with the segmentation item above β€” same frames, more pixels each. The optical half is certain; whether segmentation error stays near 1 px is not β€” earlier testing on real segmented edges gave a 0.78Γ— improvement, not 2Γ—. | 1 session | +| **Capture across a wide club-speed range** | Every shot here is a 7-iron or 9-iron with no speed overlap, which left the radar-rotation test unable to discriminate. A driver and a wedge in one session removes the confound and would settle whether the Doppler broadening is rotation. | 1 session | +| **Lux and exposure ladder in the bay** | Determines whether the optical route is a one-degree or four-degree instrument. Comparator anchors: Trackman 4 at 700–800 lux, Mevo Gen 2 at 300 lux, both continuous. | 1 session | +| **Dimension the enclosure; self-level the camera** | Camera pitch is currently recovered from footage rather than specified, so it is not reproducible on a second unit. `inclinometer.py` already tilt-compensates the radar. The acoustic timing fix also depends on a per-installation ball-to-unit distance, which belongs in the same calibration record. | days | + +### Candidate approaches, not yet evaluated + +| Approach | Rationale | Principal risk | +|---|---|---| +| **Sub-pixel edge extraction** | Masks come from a hard threshold. Given that one pixel is worth ~10Β° of face angle, boundary precision is worth more than any change to the fit. | Motion blur may already exceed sub-pixel scale β€” the club moves about 3 px during exposure. | +| **Mark the club on the measuring rig only** | The shipped product must be markerless; a calibration rig need not be. Provides per-frame truth to score the markerless estimator against. | None technical. Requires the rig to be built. | +| **Second camera** | Stereo resolves depth directly and would settle the range question outright. | Cost and synchronisation; does not address the segmentation limit, which is currently binding. | + +## 08. Reproducing this work + +Every figure here was produced by a script. Those scripts and their recorded results are kept on the fork ([`falsification/`](https://github.com/HarjotDhanota/openflight/tree/feat/silhouette-poc/research/silhouette_poc/falsification)) rather than in this repository: each answered one question once, and the answers are stated above. The library they exercise is `src/openflight/camera/clubpose/`, with its tests in `tests/`. + +| Result | Script | +|---|---| +| Silhouette information limits, per-axis | `pose_landscape.py` | +| Metric comparison and refit | `test_fusion_chamfer.py` | +| Range model test | `test_radar_range_ramp.py` | +| Doppler width / ISAR assessment | `test_isar_doppler_width.py` | +| Acoustic trigger timing | `src/openflight/acoustic.py`, `tests/test_acoustic.py` | +| Delivered loft / face angle / lie | `src/openflight/camera/clubpose/angles.py`, `tests/test_clubpose_club_angles.py` | +| Scoring primitives and their unit tests | `src/openflight/camera/clubpose/scores.py`, `tests/test_clubpose_pose_scores.py` | + +Scripts resolve the capture export via `OPENFLIGHT_SESSION`, a `--session` argument, or conventional paths, and fail with an actionable message if none is found. The club mesh is not redistributed; it is available from GrabCAD and is used for research only. + +**Contributing a capture is the most useful help.** A session with a driver and a wedge alongside the irons, at 1280Γ—800 1:1, with a lux reading at the ball, would unblock three separate items in Tier 3 at once. + +## 09. Appendix: corrections to earlier claims + +Figures published in earlier versions of this work that were subsequently found to be wrong. They are listed because some were circulated, and because the failure modes recur. + +| Claim as published | Correction | +|---|---| +| Measured loft 17.5Β° | **33.10Β°.** The detector anchored to the cavity rim on the reverse of the club. | +| The mesh has a 62 mm shaft stub | It has **no shaft**. The feature measures 63.8 mm Γ— 12.9–17.5 mm and is the hosel and ferrule. | +| Free-depth fits at 1180–1336 mm are errors of βˆ’245 to βˆ’401 mm against the tape | Those ranges lie **inside** the clubhead's physical range during the fitted frames. The fit was tracking the club; pinning it to 1581 mm moved it onto the ball. See Β§04. | +| Representative fit quality, IoU 0.636 | **Not reproducible** by any code in the repository β€” the committed tracker returns 0.292 on that shot, a careful rebuild 0.452. | +| The trigger lags impact by 2.11 frames β€” then retracted in favour of 6.0 frames | **The retraction was the error.** 2.11 frames is correct and equals the acoustic time of flight over 1.575 m. Confirmed by the ball track (71.89 Β± 0.77, n=20) and by the physics, agreeing to 0.04 frames. The 6.0 figure came from misreading a render. | +| Production carries a ~4.6 ms impact-timing bias | **Withdrawn.** Both production measurement paths derive impact from the data, not the trigger. | +| Field of view 2.17 m in the current mode | **1.08 m.** The calculation used the full sensor width where the capture reads half. | +| A 6 mm lens yields about one degree of face angle | Did not survive testing on real segmented edges, which improved by **0.78Γ—** rather than the projected 2Γ—. | +| `iwr_club_path_club_range_m` is the clubhead's range | It is the range at the **start of the radar track**, roughly 5.5 frames before impact. | +| Radar club path and attack angle are exactly equal and opposite β€” a degeneracy | True on one shot only. Across 22 shots the sum ranges +7.8Β° to βˆ’42.1Β°. | + +> **Common cause** +> All but one of the above came from generalising a single shot, a single measurement, or a geometric assumption that was never checked against the mesh or the data. The corrections came from cross-set checks and from rendering the thing in question and looking at it. Both are cheap; neither was applied first. diff --git a/docs/clubface-impact-location.md b/docs/clubface-impact-location.md new file mode 100644 index 00000000..0d2539f9 --- /dev/null +++ b/docs/clubface-impact-location.md @@ -0,0 +1,96 @@ +# Clubface impact location: status and how to help + +An investigation into measuring clubface impact location and face angle from +the existing hardware β€” the single behind-ball OV9281 camera plus the OPS243 +and IWR6843 radars, ambient light, no markers on the ball or club. + +## Read this first + +**[Technical report](clubface-impact-location-report.md)** β€” the full assessment: +what is measured, what is not, and what would resolve the rest. States every +retracted claim alongside what replaced it. The +[web version](https://claude.ai/code/artifact/c8817c34-c3ea-4455-9700-cf5a4e238b75) +carries eleven figures: real frames with the model’s own projections overlaid. + +**[Fusion status, frame by frame](https://claude.ai/code/artifact/ab9f69dd-de06-4335-83b1-29e3e29ee6b9)** +β€” two real shots with the model's own projections overlaid, nothing padded. +Thirty seconds of stepping through it conveys the state faster than any prose. + +**[Full working log](https://claude.ai/code/artifact/42a6f3f4-0b9b-4faf-bf9c-1ff45b4e94dd)** +β€” the chronological record, corrections applied in place, for tracing how any +conclusion was reached. + +## Where it stands, in three lines + +- **Validated:** ball detection (21/22), impact timing (camera and radar agree + to 0.66 frames), camera attitude (measured, not assumed), and the fused + radar+camera clubhead velocity, which matches the OPS243's independent club + speed with a mean ratio of 0.970 (sd 0.029, spread 0.941–1.015). +- **Not yet working:** clubhead orientation. Face angle, dynamic loft and + impact location remain model-dependent inferences with no accuracy figure + against truth. +- **Why:** the first 5Β° of face angle change the projected silhouette by zero + pixels; one pixel of segmentation error is worth about 10Β° of face angle; + and the club is segmentable for only ~10 pre-impact frames, of which the + current extractor keeps 3–5, against a four-parameter fit. + +## Where help is most valuable + +**Contributing a capture is the most useful thing you can do**, and you no +longer need anyone else's data to do it β€” see *Running it on your own device* +below. The current session is 21 shots of 7-iron and 9-iron from a single rig, +thin enough that several tests cannot discriminate. + +1. **A session recorded alongside a Trackman.** Nothing here has been scored + against a reference instrument, so no accuracy figure exists for any club + metric. This is the single measurement that would change that. +2. **Clubhead segmentation.** Extracting more of the ~10 frames the club + appears in roughly doubles the observations per shot. The masks currently + come from a hard background-difference threshold. +3. **A capture at 1280Γ—800 1:1** (doubles plate scale at the same frame rate) + **and across a wide club-speed range** (a driver and a wedge; the existing + session is 7-iron/9-iron with no speed overlap, which starves several + discriminating tests of power). + +## Running the code + +The library lives in `src/openflight/camera/clubpose/`, with its tests in `tests/`, +whose `README.md` maps every script to the question it answers. Per-shot +result JSONs are committed so conclusions can be re-analysed without repeating +fits that cost ~25 minutes per arm. + +### Running it on your own device + +Two inputs are not in git, and both fail closed with instructions when absent. + +**One-time setup β€” the club mesh.** Every analysis run projects the 7-iron +model, so this is needed whichever capture you use. It is a GrabCAD community +model used as local research truth and is **not redistributed**; +`src/openflight/camera/clubpose/meshes/SOURCES.md` records the source link, expected +SHA-256, and licence position, and you fetch your own copy under GrabCAD's +terms: + +```bash +uv run python \ + scripts/analysis/download_club_mesh.py --local-iron +``` + +**Then your own captures.** The library takes frames and a mesh; it has no +opinion about where your data lives. A session recorded by `start-kiosk.sh` +already contains everything needed β€” the camera `frames.npz` and the IWR6843 +`.l3dump` per shot β€” and `openflight.iwr6843.replay.inputs_from_session` +resolves those paths straight out of the session JSONL. + +Both the camera and the IWR6843 must be enabled while capturing; a shot +missing either one cannot be fitted. + +The reference **capture session** used throughout the report is +available from the maintainer if you want to reproduce its exact numbers. Your +own export works for everything else. The **7-iron mesh** is fetched from +GrabCAD (local research use only, no redistribution); +`src/openflight/camera/clubpose/meshes/SOURCES.md` has the provenance, hashes, and +download script. + +Deliberately excluded from this branch: the superseded synthetic-phase +evaluation, the old web studio, and the June–July simulation studies. They +remain on the fork's `feat/silhouette-poc` branch for archaeology. diff --git a/src/openflight/camera/clubpose/angles.py b/src/openflight/camera/clubpose/angles.py new file mode 100644 index 00000000..e752fcdd --- /dev/null +++ b/src/openflight/camera/clubpose/angles.py @@ -0,0 +1,134 @@ +"""Clubface orientation in the terms a golfer uses, and the poses that are possible. + +The fitter parameterises orientation as yaw, pitch and roll applied to the +mesh's own frame. Those numbers are not checkable by eye, are not what the +product reports, and are measured from a misleading origin: the mesh's local ++x axis points out the BACK of the club, so ``triad(0, 0, 0)`` is a clubface +aimed at the camera -- **face angle +178.7 deg, loft -19.7 deg**. + +That origin caused a real error. A fit seeded on a grid around zero searches +the neighbourhood of a backwards club and never reaches the square one, which +sits near ``pitch = -194 deg``. Seeds must come from :func:`square_pose`. + +The three angles here have known physical envelopes, so they validate a fit +without a reference instrument: a 7-iron delivered with negative loft, or a +shaft twenty degrees off its lie, is wrong whatever it scores. + +Club axes are measured off the mesh (see +the fork's feat/silhouette-poc branch), not assumed, and give +loft 33.10 deg / lie 61.19 deg -- consistent with a 690CB catalogue 34/62. +""" + +from __future__ import annotations + +import math + +import numpy as np + +from .fit import triad + +FACE_NORMAL_LOCAL = np.array([-0.941, 0.021, -0.337]) +FACE_NORMAL_LOCAL = FACE_NORMAL_LOCAL / np.linalg.norm(FACE_NORMAL_LOCAL) +SHAFT_LOCAL = np.array([-0.245, 0.295, -0.924]) +SHAFT_LOCAL = SHAFT_LOCAL / np.linalg.norm(SHAFT_LOCAL) + +STATIC_LOFT_DEG = 33.10 +STATIC_LIE_DEG = 61.19 + +# Wider than any real delivery: outside these a pose is wrong, not unusual. +ENVELOPE = { + "dynamic_loft_deg": (15.0, 50.0), + "face_angle_deg": (-25.0, 25.0), + "lie_deg": (45.0, 78.0), +} + + +def basis_from_angles(yaw_deg: float, pitch_deg: float, roll_deg: float) -> np.ndarray: + """Local-to-world matrix for a pose, columns being the mesh's own axes.""" + normal, width, height = triad(yaw_deg, pitch_deg, roll_deg) + return np.column_stack((normal, width, height)) + + +def delivered_angles(basis_world: np.ndarray) -> dict[str, float]: + """Dynamic loft, face angle and lie in degrees, from a local-to-world basis. + + World axes are x downrange, y right, z up. Face angle is the azimuth of the + face normal about vertical, so zero is square and positive is open for a + right-handed player. + """ + basis = np.asarray(basis_world, dtype=float) + face = basis @ FACE_NORMAL_LOCAL + shaft = basis @ SHAFT_LOCAL + face_norm, shaft_norm = np.linalg.norm(face), np.linalg.norm(shaft) + # Fail closed: NaN angles would read as "implausible pose", not "broken input". + if ( + not (np.isfinite(face_norm) and np.isfinite(shaft_norm)) + or min(face_norm, shaft_norm) < 1e-9 + ): + raise ValueError("degenerate basis: club axes collapse to zero length") + face = face / face_norm + shaft = shaft / shaft_norm + return { + "dynamic_loft_deg": math.degrees(math.asin(float(np.clip(face[2], -1.0, 1.0)))), + "face_angle_deg": math.degrees(math.atan2(float(face[1]), float(face[0]))), + "lie_deg": math.degrees(math.asin(float(np.clip(abs(shaft[2]), -1.0, 1.0)))), + } + + +def angles_from_pose(yaw_deg: float, pitch_deg: float, roll_deg: float) -> dict[str, float]: + """Convenience wrapper: pose angles straight to delivered angles.""" + return delivered_angles(basis_from_angles(yaw_deg, pitch_deg, roll_deg)) + + +def in_envelope(angles: dict[str, float]) -> bool: + """Could a real club have been delivered in this orientation?""" + return all(low <= angles[key] <= high for key, (low, high) in ENVELOPE.items()) + + +def square_pose( + dynamic_loft_deg: float = STATIC_LOFT_DEG, + face_angle_deg: float = 0.0, + lie_deg: float = STATIC_LIE_DEG, +) -> tuple[float, float, float]: + """The (yaw, pitch, roll) that delivers the requested angles. + + Solved rather than hard-coded, so it stays correct if the mesh, its + measured axes, or ``triad``'s convention change. Used to seed fits: a grid + around the origin searches the neighbourhood of a backwards club. + """ + from scipy.optimize import minimize # noqa: PLC0415 + + target = { + "dynamic_loft_deg": float(dynamic_loft_deg), + "face_angle_deg": float(face_angle_deg), + "lie_deg": float(lie_deg), + } + + def cost(params: np.ndarray) -> float: + angles = angles_from_pose(*params) + # Circular difference: near +-180 is not 360 degrees from its target. + total = 0.0 + for key, want in target.items(): + delta = angles[key] - want + if key == "face_angle_deg": + delta = math.remainder(delta, 360.0) + total += delta * delta + return total + + best = None + for yaw in (-90.0, 0.0, 90.0, 180.0): + for pitch in (-180.0, -90.0, 0.0, 90.0): + for roll in (-90.0, 0.0, 90.0, 180.0): + result = minimize( + cost, + np.array([yaw, pitch, roll], dtype=float), + method="Nelder-Mead", + options={"maxiter": 2000, "xatol": 1e-5, "fatol": 1e-10}, + ) + if best is None or result.fun < best.fun: + best = result + if best is None or best.fun > 1e-3: + raise RuntimeError( + f"no pose delivers {target}; residual {None if best is None else best.fun}" + ) + return tuple(float(v) for v in best.x) diff --git a/src/openflight/camera/clubpose/fit.py b/src/openflight/camera/clubpose/fit.py new file mode 100644 index 00000000..f5f2c6ad --- /dev/null +++ b/src/openflight/camera/clubpose/fit.py @@ -0,0 +1,400 @@ +"""Fit the real 3D club mesh to a real camera frame. + +Everything the POC has measured so far came from fitting the mesh to SYNTHETIC +silhouettes produced by the same mesh. This module points it at real pixels for +the first time, which is the actual thing under validation. + +Three corrections are required before the existing machinery can be used at all, +because the shipped `A0` preset describes a camera we do not have: + + A0 says fx = 1033 px, plate scale 0.656 px/mm, range 1575 mm + measured fx = 466.7 px, plate scale 0.295 px/mm, range ~1581 mm + +`fx` follows from the NOMINAL datasheet lens (2.8 mm) over the effective pixel +pitch of the shipped 320x200 mode (3.0 um at 2x subsample = 6.0 um). It is not +a calibrated camera matrix: there is no distortion model, no separately +estimated principal point, and no independent fx/fy. + +The RANGE was previously 1425 mm, from a 13.97 px ball. That is wrong. Across +the 21 correctly exposed shots of session 20260825_181734 the teed ball +measures 12.77 px, and the tape chain -- camera lens 203.2 mm above the floor +(kiosk log `mount_height_m`), ball centre 40 mm, radar slant tee range 1575 mm, +camera lateral offset -60.325 mm -- gives 1581 mm. The 13.97 px figure traces +to the capture that turned out 99.8 % clipped, where the ball bloomed. + +That matters more than a 10 % scale error. Both `range_grid_mm` defaults below +spanned 1300-1550 and 1325-1525, so neither CONTAINED the true range. The local +refinement below is not bounded by the grid, so it could in principle climb out +-- but it hill-climbs greedily from the best COARSE pose, and that pose was +selected at the wrong depth. Reaching the truth was therefore left to chance +rather than to the search. See +the technical report's clubhead-range section. + +The pose model here is the one the POC already uses: 3D centre plus roll about +the face normal. That is **4 degrees of freedom, not 6** - the face normal is +fixed by `FACE_NORMAL` rather than solved. Loft and lie are therefore baked into +the mesh's own frame and are not recovered. Anything reported here is a fit of +position and roll only. +""" + +from __future__ import annotations + +import math +from dataclasses import dataclass + +import numpy as np + +from openflight.camera.clubpose.mesh import rasterize_projected_triangles +from openflight.camera.clubpose.projection import ( + CAMERA_CENTER_WORLD, + FACE_NORMAL, + CameraPreset, + _face_axes, + _project, + _ray_world, +) + +# Measured configuration of the shipped camera. NOT the A0 preset. +LENS_MM = 2.8 +PITCH_UM = 3.0 +SUBSAMPLE = 2 +FOCAL_PX = LENS_MM / (PITCH_UM * SUBSAMPLE * 1e-3) # 466.7 +# Tape chain, not a ball measurement: see the module docstring. +CAMERA_BALL_RANGE_MM = 1581.0 + + +def measured_camera(width: int = 320, height: int = 200) -> CameraPreset: + """The camera we actually have, from datasheet optics and measured range.""" + return CameraPreset( + name="MEASURED", + width=width, + height=height, + fx=FOCAL_PX, + fy=FOCAL_PX, + cx=width / 2.0, + cy=height / 2.0, + plate_scale_px_per_mm=FOCAL_PX / CAMERA_BALL_RANGE_MM, + sensor_crop=(336, 150, 816, 516), + sampling_increment=(SUBSAMPLE, SUBSAMPLE), + isp_offset=(4, 4), + orientation="rot180", + gate_b1_passed=False, + physical_status="measured_from_real_capture", + ) + + +def render_mask(mesh, center_world, roll_rad, camera) -> tuple[np.ndarray, np.ndarray] | None: + """Project and rasterise the mesh, against an explicitly supplied camera. + + The evaluation copy of this resolves the camera by preset NAME, which hard- + codes the wrong intrinsics for real data. This one takes the camera object. + """ + axis_u, axis_v = _face_axes(float(roll_rad)) + local = mesh.vertices_local_mm + world = ( + np.asarray(center_world, dtype=float)[None, :] + + local[:, 0, None] * FACE_NORMAL[None, :] + + local[:, 1, None] * axis_u[None, :] + + local[:, 2, None] * axis_v[None, :] + ) + uv, front = _project(world, camera) + center_uv, center_front = _project(np.asarray(center_world, dtype=float)[None, :], camera) + if not bool(center_front[0]) or not front.any(): + return None + faces = mesh.faces[np.all(front[mesh.faces], axis=1)] + if not len(faces): + return None + mask = rasterize_projected_triangles(uv, faces, width=camera.width, height=camera.height) + return mask, center_uv[0] + + +def iou(a: np.ndarray, b: np.ndarray) -> float: + a = a.astype(bool) + b = b.astype(bool) + union = np.count_nonzero(a | b) + return float(np.count_nonzero(a & b) / union) if union else 0.0 + + +@dataclass +class RealFit: + ok: bool + reason: str + iou: float + center_world: np.ndarray | None + roll_deg: float | None + range_mm: float | None + observed_px: int + rendered_px: int + + +def fit_frame( + mesh, + observed_mask: np.ndarray, + camera: CameraPreset, + *, + range_grid_mm=np.arange(1250.0, 1651.0, 50.0), + roll_grid_deg=np.arange(-90.0, 90.0, 7.5), + refine: bool = True, +) -> RealFit: + """Best 3D centre + roll for one observed club silhouette, by direct IoU. + + The centre is constrained to the pixel ray through the observed centroid, so + the search is over range and roll rather than free 3D translation. That is + the same constraint the radar range sphere would impose, standing in for a + radar measurement this camera-only capture does not have. + """ + observed = observed_mask.astype(bool) + n_obs = int(observed.sum()) + if n_obs < 40: + return RealFit(False, "observed_mask_too_small", 0.0, None, None, None, n_obs, 0) + + ys, xs = np.nonzero(observed) + centroid = np.array([xs.mean(), ys.mean()], dtype=float) + ray = _ray_world(centroid, camera) + + def point_at(range_mm: float) -> np.ndarray: + # Range is from the CAMERA, not the world origin at the impact point. + return CAMERA_CENTER_WORLD + ray * float(range_mm) + + best = (0.0, None, None, None, 0) + for range_mm in range_grid_mm: + center = point_at(range_mm) + for roll_deg in roll_grid_deg: + out = render_mask(mesh, center, math.radians(float(roll_deg)), camera) + if out is None: + continue + score = iou(out[0], observed) + if score > best[0]: + best = (score, center, float(roll_deg), float(range_mm), int(out[0].sum())) + + if best[1] is None: + return RealFit(False, "no_pose_projected", 0.0, None, None, None, n_obs, 0) + + if refine: + score, center, roll_deg, range_mm, n_ren = best + for _ in range(3): + improved = False + for d_range in (-25.0, -10.0, 10.0, 25.0): + for d_roll in (-4.0, -1.5, 1.5, 4.0): + cand_center = point_at(range_mm + d_range) + out = render_mask(mesh, cand_center, math.radians(roll_deg + d_roll), camera) + if out is None: + continue + cand = iou(out[0], observed) + if cand > score: + score, center, roll_deg, range_mm = ( + cand, + cand_center, + roll_deg + d_roll, + range_mm + d_range, + ) + n_ren, improved = int(out[0].sum()), True + if not improved: + break + best = (score, center, roll_deg, range_mm, n_ren) + + score, center, roll_deg, range_mm, n_ren = best + return RealFit(True, "ok", score, center, roll_deg, range_mm, n_obs, n_ren) + + +# 6-DOF fitting; the earlier 4-DOF model could not represent loft, lie or face. + + +def _rot(axis: np.ndarray, angle_rad: float) -> np.ndarray: + """Rodrigues rotation about an arbitrary axis.""" + a = np.asarray(axis, dtype=float) + a = a / np.linalg.norm(a) + K = np.array([[0.0, -a[2], a[1]], [a[2], 0.0, -a[0]], [-a[1], a[0], 0.0]]) + return np.eye(3) + math.sin(angle_rad) * K + (1.0 - math.cos(angle_rad)) * (K @ K) + + +def triad(yaw_deg: float, pitch_deg: float, roll_deg: float) -> tuple[np.ndarray, ...]: + """Full orientation as an orthonormal (normal, u, v) triad. + + yaw - about world up, i.e. FACE ANGLE (open/closed) + pitch - about world right, i.e. DYNAMIC LOFT + roll - about the face normal, i.e. LIE / toe-up rotation + """ + from openflight.camera.clubpose.projection import WORLD_RIGHT, WORLD_UP + + R = _rot(WORLD_UP, math.radians(yaw_deg)) @ _rot(WORLD_RIGHT, math.radians(pitch_deg)) + n = R @ FACE_NORMAL + u = R @ WORLD_RIGHT + v = R @ WORLD_UP + Rr = _rot(n, math.radians(roll_deg)) + return n, Rr @ u, Rr @ v + + +def render_mask_6dof(mesh, center_world, yaw_deg, pitch_deg, roll_deg, camera): + """Project and rasterise with a FULL orientation rather than roll alone.""" + n, u, v = triad(yaw_deg, pitch_deg, roll_deg) + local = mesh.vertices_local_mm + world = ( + np.asarray(center_world, dtype=float)[None, :] + + local[:, 0, None] * n[None, :] + + local[:, 1, None] * u[None, :] + + local[:, 2, None] * v[None, :] + ) + uv, front = _project(world, camera) + _, center_front = _project(np.asarray(center_world, dtype=float)[None, :], camera) + if not bool(center_front[0]) or not front.any(): + return None + faces = mesh.faces[np.all(front[mesh.faces], axis=1)] + if not len(faces): + return None + return rasterize_projected_triangles(uv, faces, width=camera.width, height=camera.height) + + +def fit_frame_6dof( + mesh, + observed_mask, + camera, + *, + range_grid_mm=(1456.0, 1581.0, 1706.0), + yaw_grid=(-40.0, -20.0, 0.0, 20.0, 40.0), + pitch_grid=(-40.0, -20.0, 0.0, 20.0, 40.0), + roll_grid=(-60.0, -30.0, 0.0, 30.0, 60.0, 90.0), +): + """Best 6-DOF pose by direct IoU. Coarse grid, then local refinement.""" + observed = observed_mask.astype(bool) + if int(observed.sum()) < 40: + return {"ok": False, "reason": "observed_mask_too_small", "iou": 0.0} + ys, xs = np.nonzero(observed) + ray = _ray_world(np.array([xs.mean(), ys.mean()], dtype=float), camera) + + def score(rng, yaw, pitch, roll): + m = render_mask_6dof(mesh, CAMERA_CENTER_WORLD + ray * rng, yaw, pitch, roll, camera) + return (0.0, None) if m is None else (iou(m, observed), m) + + best = (0.0, None) + for rng in range_grid_mm: + for yaw in yaw_grid: + for pitch in pitch_grid: + for roll in roll_grid: + s, _ = score(rng, yaw, pitch, roll) + if s > best[0]: + best = (s, (rng, yaw, pitch, roll)) + if best[1] is None: + return {"ok": False, "reason": "no_pose_projected", "iou": 0.0} + + rng, yaw, pitch, roll = best[1] + step = [60.0, 10.0, 10.0, 15.0] + for _ in range(4): + improved = False + for k, deltas in enumerate(step): + for d in (-deltas, deltas): + cand = [rng, yaw, pitch, roll] + cand[k] += d + s, _ = score(*cand) + if s > best[0]: + best, (rng, yaw, pitch, roll), improved = (s, tuple(cand)), cand, True + if not improved: + step = [x / 2.0 for x in step] + return { + "ok": True, + "reason": "ok", + "iou": best[0], + "range_mm": rng, + "yaw_deg": yaw, + "pitch_deg": pitch, + "roll_deg": roll, + } + + +# Independent per-frame fits jumped >100 deg between frames: a 20-40 px +# silhouette under-determines six DOF, so frames share bounds and smoothness. + +# Loose sanity bounds only: yaw/pitch/roll are offsets from the mesh's own +# normalised frame, not face angle / loft / lie. +YAW_BOUND_DEG = 60.0 +PITCH_RANGE_DEG = (-40.0, 90.0) +ROLL_BOUND_DEG = 70.0 + + +def fit_sequence( + mesh, + masks: dict[int, np.ndarray], + camera: CameraPreset, + *, + smooth_deg: float = 70.0, + smooth_mm: float = 300.0, + range_grid_mm=(1481.0, 1581.0, 1681.0), + yaw_grid=(-40.0, -20.0, 0.0, 20.0, 40.0), + pitch_grid=(-30.0, 0.0, 30.0, 60.0, 85.0), + roll_grid=(-60.0, -30.0, 0.0, 30.0, 60.0), + refine_range: bool = True, +) -> dict[int, dict]: + """Fit an ordered run of frames, penalising jumps between consecutive poses. + + Score is IoU minus a smoothness penalty against the previous accepted pose. + `smooth_deg` and `smooth_mm` set how much orientation and range change costs + one unit of IoU, so a large IoU gain can still justify real motion while noise + cannot. Set ``refine_range=False`` with a singleton ``range_grid_mm`` to keep + an externally measured range hard-pinned during local refinement. + """ + out: dict[int, dict] = {} + prev = None + for i in sorted(masks): + observed = masks[i].astype(bool) + if int(observed.sum()) < 40: + continue + ys, xs = np.nonzero(observed) + ray = _ray_world(np.array([xs.mean(), ys.mean()], dtype=float), camera) + + def score(rng, yaw, pitch, roll): + if abs(yaw) > YAW_BOUND_DEG or not PITCH_RANGE_DEG[0] <= pitch <= PITCH_RANGE_DEG[1]: + return -1.0 + if abs(roll) > ROLL_BOUND_DEG: + return -1.0 + m = render_mask_6dof(mesh, CAMERA_CENTER_WORLD + ray * rng, yaw, pitch, roll, camera) + if m is None: + return -1.0 + value = iou(m, observed) + if prev is not None: + d_ang = ( + abs(yaw - prev["yaw_deg"]) + + abs(pitch - prev["pitch_deg"]) + + abs(roll - prev["roll_deg"]) + ) + d_rng = abs(rng - prev["range_mm"]) + value -= d_ang / (3.0 * smooth_deg) + d_rng / smooth_mm + return value + + best = (-1.0, None) + for rng in range_grid_mm: + for yaw in yaw_grid: + for pitch in pitch_grid: + for roll in roll_grid: + s = score(rng, yaw, pitch, roll) + if s > best[0]: + best = (s, (rng, yaw, pitch, roll)) + if best[1] is None: + continue + rng, yaw, pitch, roll = best[1] + step = [50.0, 8.0, 8.0, 8.0] + for _ in range(4): + improved = False + for k, delta in enumerate(step): + if k == 0 and not refine_range: + continue + for d in (-delta, delta): + cand = [rng, yaw, pitch, roll] + cand[k] += d + s = score(*cand) + if s > best[0]: + best = (s, tuple(cand)) + rng, yaw, pitch, roll = cand + improved = True + if not improved: + step = [x / 2.0 for x in step] + + m = render_mask_6dof(mesh, CAMERA_CENTER_WORLD + ray * rng, yaw, pitch, roll, camera) + prev = { + "range_mm": rng, + "yaw_deg": yaw, + "pitch_deg": pitch, + "roll_deg": roll, + "iou": iou(m, observed) if m is not None else 0.0, + } + out[i] = dict(prev, mask=m) + return out diff --git a/src/openflight/camera/clubpose/head_split.py b/src/openflight/camera/clubpose/head_split.py new file mode 100644 index 00000000..acff4dd5 --- /dev/null +++ b/src/openflight/camera/clubpose/head_split.py @@ -0,0 +1,71 @@ +"""Separate the clubhead from the shaft inside one moving component. + +MEASURED basis (session 20260825_181734, shots 002/005/014/021/029, +frames 58-88, 205 components): + + isolated shaft max inscribed radius 2.0 - 3.0 px, 0.0% of pixels >= 4 px + isolated head max inscribed radius 5.0 - 11.6 px, 12 - 49% of pixels >= 4 px + +The two never overlap, so "thicker than any shaft" is a measured discriminator, +not a tuned one. The shaft is a thin line; the head is a compact body. + +The head core (distance transform >= HEAD_CORE_PX) is used only as a SEED. The +returned mask is the original component's own pixels, partitioned by a watershed +on the distance transform, so the head's outline is the observed silhouette +boundary - never eroded, never padded. Only the cut across the hosel neck is +synthetic, and that is where the head genuinely ends. +""" + +from __future__ import annotations + +import cv2 +import numpy as np + +HEAD_CORE_PX = 4.0 # exceeds every measured shaft inscribed radius (max 3.0) +SHAFT_MAX_PX = 3.5 # below every measured head inscribed radius (min 5.0) +# Measured over 205 components: head alone reaches 4-41 px, head+shaft +# 145-187 px. 60 sits inside the empty gap. +SHAFT_REACH_PX = 60.0 + + +def split_head(component: np.ndarray) -> tuple[np.ndarray, np.ndarray] | None: + """(head_mask, shaft_mask) for one connected moving component. + + Returns None when the component contains no body thick enough to be a + clubhead - fail closed rather than hand back a piece of shaft. + """ + c = (component > 0).astype(np.uint8) + dt = cv2.distanceTransform(c, cv2.DIST_L2, 5) + core = (dt >= HEAD_CORE_PX).astype(np.uint8) + n, lab, st, _ = cv2.connectedComponentsWithStats(core, 8) + if n <= 1: + return None + head_label = 1 + int(np.argmax(st[1:, 4])) + head_core = (lab == head_label).astype(np.uint8) + + reach = cv2.distanceTransform(1 - head_core, cv2.DIST_L2, 5) + far = c.copy() + far[reach <= SHAFT_REACH_PX] = 0 # farther than any head extends + far[dt >= SHAFT_MAX_PX] = 0 # and genuinely thin + + markers = np.zeros(c.shape, np.int32) + markers[c == 0] = 1 # background + markers[head_core > 0] = 2 + markers[far > 0] = 3 + if not (markers == 3).any(): + return c, np.zeros_like(c) # nothing thin and remote: all head + + relief = (255.0 * (1.0 - dt / max(dt.max(), 1e-6))).astype(np.uint8) + cv2.watershed(cv2.cvtColor(relief, cv2.COLOR_GRAY2BGR), markers) + head = ((markers == 2) & (c > 0)).astype(np.uint8) + shaft = ((markers == 3) & (c > 0)).astype(np.uint8) + # watershed marks its ridge -1; award those pixels to whichever side is nearer + ridge = (markers == -1) & (c > 0) + if ridge.any(): + dh = cv2.distanceTransform(1 - head, cv2.DIST_L2, 5) + ds = cv2.distanceTransform(1 - np.clip(shaft, 0, 1), cv2.DIST_L2, 5) + head[ridge & (dh <= ds)] = 1 + shaft[ridge & (dh > ds)] = 1 + if int(head.sum()) < 60: + return None + return head, shaft diff --git a/src/openflight/camera/clubpose/motion.py b/src/openflight/camera/clubpose/motion.py new file mode 100644 index 00000000..1380e2bf --- /dev/null +++ b/src/openflight/camera/clubpose/motion.py @@ -0,0 +1,39 @@ +"""Sequence-level rigid-rotation primitives for clubhead pose experiments.""" + +from __future__ import annotations + +import math + +import numpy as np + + +def axis_from_angles(azimuth_deg: float, elevation_deg: float) -> np.ndarray: + """Return a unit world axis from azimuth and elevation in degrees.""" + azimuth = math.radians(float(azimuth_deg)) + elevation = math.radians(float(elevation_deg)) + horizontal = math.cos(elevation) + return np.asarray( + ( + horizontal * math.cos(azimuth), + horizontal * math.sin(azimuth), + math.sin(elevation), + ), + dtype=float, + ) + + +def constrained_omega_deg_s( + speed_mps: float, + swing_radius_m: float, + axis_azimuth_deg: float, + axis_elevation_deg: float, +) -> np.ndarray: + """Return angular velocity whose magnitude is fixed by ``omega = v / r``.""" + speed = float(speed_mps) + radius = float(swing_radius_m) + if not math.isfinite(speed) or speed <= 0.0: + raise ValueError("club speed must be finite and positive") + if not math.isfinite(radius) or radius <= 0.0: + raise ValueError("swing radius must be finite and positive") + magnitude_deg_s = math.degrees(speed / radius) + return magnitude_deg_s * axis_from_angles(axis_azimuth_deg, axis_elevation_deg) diff --git a/src/openflight/camera/clubpose/scores.py b/src/openflight/camera/clubpose/scores.py new file mode 100644 index 00000000..01975b6a --- /dev/null +++ b/src/openflight/camera/clubpose/scores.py @@ -0,0 +1,193 @@ +"""Scores for a clubhead pose that are not silhouette IoU. + +Section 11f measured IoU running *inversely* to pose correctness on real +segmented masks: the arms that recovered a worse pose scored a better overlap. +That makes every IoU-scored conclusion in this project unreadable rather than +wrong -- including the radar-constrained rotation experiment, which reports a +0.041 IoU penalty for imposing ``|omega| = v / r`` and cannot say whether that +penalty means the constraint is wrong or means the constraint is right and IoU +is punishing it for being right. + +Two replacements, chosen because they fail differently from IoU and from each +other: + +``chamfer_px`` + Symmetric mean distance between mask BOUNDARIES. IoU on a 200-pixel blob is + dominated by area, so a pose that is the right size in the right place + scores well while pointing the wrong way. Edge distance is dominated by + shape, which is the part we cannot currently read. + +``omega_residual_deg_s`` + A pose PAIR implies an angular velocity. The radar independently fixes its + magnitude at ``v / r``. This is the only score here that uses information + from outside the image, and it is the one Trackman's OERT leans on -- their + 720p 60 fps camera cannot track impact alone either. + +Neither is validated against truth on real pixels. Nothing in this module +should be quoted as an accuracy figure. ``openflight.camera.clubpose/tests/test_pose_scores.py`` shows +only that they rank a known pose first on synthetic masks, which establishes +they are not broken, not that they work. +""" + +from __future__ import annotations + +import json +import math +from dataclasses import dataclass +from pathlib import Path + +import cv2 +import numpy as np + +from .fit import iou as _iou, render_mask_6dof, triad + + +@dataclass(frozen=True) +class Pose: + """A full clubhead pose: where the centre is, and how it is oriented. + + Angles follow ``fit.triad``: yaw about world up is FACE ANGLE, pitch + about world right is DYNAMIC LOFT, roll about the face normal is LIE. + """ + + center_world: tuple[float, float, float] + yaw_deg: float + pitch_deg: float + roll_deg: float + + def as_dict(self) -> dict: + """Plain-JSON form, for annotation files.""" + return { + "center_world_mm": [float(v) for v in self.center_world], + "yaw_deg": float(self.yaw_deg), + "pitch_deg": float(self.pitch_deg), + "roll_deg": float(self.roll_deg), + } + + @classmethod + def from_dict(cls, d: dict) -> Pose: + """Rebuild a pose from `as_dict` output.""" + return cls( + tuple(float(v) for v in d["center_world_mm"]), + float(d["yaw_deg"]), + float(d["pitch_deg"]), + float(d["roll_deg"]), + ) + + def basis(self) -> np.ndarray: + """Orientation as a rotation matrix whose columns are the triad.""" + n, u, v = triad(self.yaw_deg, self.pitch_deg, self.roll_deg) + return np.column_stack((n, u, v)) + + def render(self, mesh, camera) -> np.ndarray | None: + """Rasterise this pose through the SAME renderer the fitter uses.""" + return render_mask_6dof( + mesh, + np.asarray(self.center_world, float), + self.yaw_deg, + self.pitch_deg, + self.roll_deg, + camera, + ) + + +@dataclass(frozen=True) +class PoseScores: + """Three readings of one pose. Directions differ, so they are named.""" + + iou: float # higher is better + chamfer_px: float # LOWER is better + omega_residual_deg_s: float | None # LOWER is better; None without motion + + +def mask_edge(mask: np.ndarray) -> np.ndarray: + """Boundary pixels of a mask: inside it, but touching the outside.""" + m = np.ascontiguousarray(mask.astype(np.uint8)) + eroded = cv2.erode(m, np.ones((3, 3), np.uint8), borderValue=0) + return (m.astype(bool)) & (~eroded.astype(bool)) + + +def chamfer_px(rendered: np.ndarray, observed: np.ndarray) -> float: + """Symmetric mean boundary distance in pixels. Lower is better. + + Returns ``inf`` when either mask has no boundary. That is a fail-closed + choice: an empty render against a real mask is a total failure, and + reporting 0.0 for it would make the worst possible pose the best scoring. + """ + ea, eb = mask_edge(rendered), mask_edge(observed) + if not ea.any() or not eb.any(): + return math.inf + dt_a = cv2.distanceTransform((~ea).astype(np.uint8), cv2.DIST_L2, 3) + dt_b = cv2.distanceTransform((~eb).astype(np.uint8), cv2.DIST_L2, 3) + return 0.5 * (float(dt_b[ea].mean()) + float(dt_a[eb].mean())) + + +def observed_omega_deg_s(pose_a: Pose, pose_b: Pose, dt_s: float) -> float: + """Angular speed implied by two poses, in degrees per second.""" + dt = float(dt_s) + if not math.isfinite(dt) or dt <= 0.0: + raise ValueError("dt_s must be finite and positive") + relative = pose_b.basis() @ pose_a.basis().T + cos_angle = (float(np.trace(relative)) - 1.0) / 2.0 + return math.degrees(math.acos(max(-1.0, min(1.0, cos_angle)))) / dt + + +def omega_residual_deg_s( + pose_a: Pose, pose_b: Pose, dt_s: float, speed_mps: float, swing_radius_m: float +) -> float: + """How far the pose pair's rotation rate sits from the radar's ``v / r``. + + The radar fixes the MAGNITUDE only; the axis stays free. So this constrains + one number, not three, and a zero residual does not mean the pose is right. + """ + speed, radius = float(speed_mps), float(swing_radius_m) + if not math.isfinite(speed) or speed <= 0.0: + raise ValueError("club speed must be finite and positive") + if not math.isfinite(radius) or radius <= 0.0: + raise ValueError("swing radius must be finite and positive") + required = math.degrees(speed / radius) + return abs(observed_omega_deg_s(pose_a, pose_b, dt_s) - required) + + +def score_pose(mesh, pose: Pose, observed_mask, camera, motion=None) -> PoseScores: + """Score one pose against one observed mask. + + ``motion``, when given, is ``(previous_pose, dt_s, speed_mps, radius_m)`` + and enables the only score that uses non-image information. + """ + observed = np.asarray(observed_mask).astype(bool) + rendered = pose.render(mesh, camera) + if rendered is None: + return PoseScores(0.0, math.inf, None) + residual = None + if motion is not None: + prev, dt_s, speed_mps, radius_m = motion + residual = omega_residual_deg_s(prev, pose, dt_s, speed_mps, radius_m) + return PoseScores(_iou(rendered, observed), chamfer_px(rendered, observed), residual) + + +def save_annotation( + path, *, shot: int, frame: int, pose: Pose, annotator: str, pass_index: int +) -> None: + """Append one hand-labelled pose. Never overwrites: the experiment IS the + comparison between repeated passes, so losing an earlier pass loses it.""" + record = { + "shot": int(shot), + "frame": int(frame), + "annotator": str(annotator), + "pass_index": int(pass_index), + "pose": pose.as_dict(), + } + p = Path(path) + p.parent.mkdir(parents=True, exist_ok=True) + with p.open("a", encoding="utf-8") as handle: + handle.write(json.dumps(record) + "\n") + + +def load_annotations(path) -> list[dict]: + """Every labelled pose from one file; empty when it does not exist yet.""" + p = Path(path) + if not p.exists(): + return [] + with p.open(encoding="utf-8") as handle: + return [json.loads(line) for line in handle if line.strip()] diff --git a/tests/test_clubpose_club_angles.py b/tests/test_clubpose_club_angles.py new file mode 100644 index 00000000..9f15e95d --- /dev/null +++ b/tests/test_clubpose_club_angles.py @@ -0,0 +1,98 @@ +"""Is the orientation measuring stick itself correct? + +Everything that judges a pose physically possible depends on this conversion, +so it is checked against facts known independently of the fitter: the mesh's +own catalogue geometry, and the fact that a square club must be expressible. +""" + +from __future__ import annotations + +import math + +import numpy as np +import pytest + +from openflight.camera.clubpose.angles import ( + ENVELOPE, + STATIC_LIE_DEG, + STATIC_LOFT_DEG, + angles_from_pose, + basis_from_angles, + delivered_angles, + in_envelope, + square_pose, +) + + +class TestSquarePose: + def test_a_square_club_is_expressible(self): + """If no pose delivers the club's own static geometry, the conversion + or the parameterisation is wrong and every verdict built on it is void.""" + got = angles_from_pose(*square_pose()) + assert got["dynamic_loft_deg"] == pytest.approx(STATIC_LOFT_DEG, abs=0.05) + assert got["face_angle_deg"] == pytest.approx(0.0, abs=0.05) + assert got["lie_deg"] == pytest.approx(STATIC_LIE_DEG, abs=0.05) + + def test_the_square_pose_is_nowhere_near_the_origin(self): + """This is the whole reason the module exists: a fit seeded on a grid + around zero never reaches the square club.""" + yaw, pitch, roll = square_pose() + assert abs(math.remainder(pitch, 360.0)) > 90.0, ( + f"pitch {pitch} is near zero, so the seeding hazard has gone away " + "and the warning in the docstring should be revisited" + ) + + def test_the_origin_is_a_backwards_club(self): + origin = angles_from_pose(0.0, 0.0, 0.0) + assert abs(origin["face_angle_deg"]) > 150.0 + assert not in_envelope(origin) + + def test_solves_for_a_requested_non_square_delivery(self): + got = angles_from_pose(*square_pose(dynamic_loft_deg=28.0, face_angle_deg=-4.0)) + assert got["dynamic_loft_deg"] == pytest.approx(28.0, abs=0.05) + assert got["face_angle_deg"] == pytest.approx(-4.0, abs=0.05) + + +class TestDeliveredAngles: + def test_basis_is_orthonormal(self): + basis = basis_from_angles(11.0, -37.0, 64.0) + assert np.allclose(basis.T @ basis, np.eye(3), atol=1e-9) + + def test_face_angle_sign_is_open_to_the_right(self): + """A right-handed player's open face points right, which is +y.""" + yaw, pitch, roll = square_pose(face_angle_deg=6.0) + assert angles_from_pose(yaw, pitch, roll)["face_angle_deg"] > 0 + + def test_angles_are_continuous_under_small_perturbation(self): + base = square_pose() + before = angles_from_pose(*base) + after = angles_from_pose(base[0] + 1.0, base[1], base[2]) + for key in ("dynamic_loft_deg", "lie_deg"): + assert abs(after[key] - before[key]) < 5.0 + assert abs(math.remainder(after["face_angle_deg"] - before["face_angle_deg"], 360.0)) < 5.0 + + def test_rejects_a_degenerate_basis(self): + with pytest.raises(ValueError): + delivered_angles([[0, 0, 0], [0, 0, 0], [0, 0, 0]]) + + +class TestEnvelope: + def test_a_square_club_is_inside(self): + assert in_envelope(angles_from_pose(*square_pose())) + + def test_negative_loft_is_rejected(self): + """A lofted iron cannot present a downward-facing face at contact. + An earlier fit produced exactly this and had to be caught.""" + assert not in_envelope({"dynamic_loft_deg": -7.1, "face_angle_deg": 16.7, "lie_deg": 28.2}) + + def test_a_backwards_face_is_rejected(self): + assert not in_envelope({"dynamic_loft_deg": 27.3, "face_angle_deg": -163.0, "lie_deg": 5.2}) + + def test_a_realistic_delivery_is_accepted(self): + assert in_envelope({"dynamic_loft_deg": 27.5, "face_angle_deg": -1.8, "lie_deg": 62.4}) + + def test_envelope_covers_the_static_geometry_with_margin(self): + low, high = ENVELOPE["dynamic_loft_deg"] + assert low < STATIC_LOFT_DEG < high + low, high = ENVELOPE["lie_deg"] + assert low < STATIC_LIE_DEG < high diff --git a/tests/test_clubpose_fit_real_range_grid.py b/tests/test_clubpose_fit_real_range_grid.py new file mode 100644 index 00000000..3b0dd361 --- /dev/null +++ b/tests/test_clubpose_fit_real_range_grid.py @@ -0,0 +1,56 @@ +"""The mesh fitter's depth search must contain the real camera-to-ball range. + +`fit.py` inherited a 1425 mm camera-to-ball range from a 13.97 px ball +measurement. Across the 21 correctly exposed shots of session 20260825_181734 +the teed ball measures 12.77 px, giving 1560 mm, and the tape gives 1581 mm +(camera lens 203.2 mm, ball centre 40 mm, radar slant tee range 1575 mm). The +13.97 px figure traces to the capture that was 99.8 % clipped, where the ball +bloomed. + +That is not an inaccuracy, it is a fail-closed violation: both shipped +`range_grid_mm` grids spanned 1300-1550 and 1325-1525, so the true range sat +at or outside the edge and the search could not reach it. +""" + +from __future__ import annotations + +import inspect +import math + +from openflight.camera.clubpose import fit + +# Tape chain, all measured: kiosk log mount_height_m, cal json radar_height_m, +# server.py --iwr6843-tee-m / --iwr6843-ball-height-m defaults. +CAMERA_HEIGHT_M = 0.2032 +CAMERA_LATERAL_M = -0.060325 +RADAR_HEIGHT_M = 0.1524 +TEE_RANGE_M = 1.575 +BALL_HEIGHT_M = 0.040 + + +def tape_camera_ball_range_mm() -> float: + forward = math.sqrt(TEE_RANGE_M**2 - (BALL_HEIGHT_M - RADAR_HEIGHT_M) ** 2) + return 1000.0 * math.hypot( + math.hypot(forward, CAMERA_LATERAL_M), CAMERA_HEIGHT_M - BALL_HEIGHT_M + ) + + +def test_tape_chain_gives_about_1580_mm(): + """Guard the reference value itself, so the grids below have an anchor.""" + assert tape_camera_ball_range_mm() == __import__("pytest").approx(1580.0, abs=8.0) + + +def test_measured_camera_preset_uses_the_tape_range(): + camera = fit.measured_camera() + implied = fit.FOCAL_PX / camera.plate_scale_px_per_mm + assert implied == __import__("pytest").approx(tape_camera_ball_range_mm(), abs=25.0) + + +def test_range_grids_bracket_the_true_range(): + """Every default depth grid must contain the measured range, with margin.""" + truth = tape_camera_ball_range_mm() + for func in (fit.fit_frame, fit.fit_frame_6dof, fit.fit_sequence): + grid = inspect.signature(func).parameters["range_grid_mm"].default + assert min(grid) < truth < max(grid), ( + f"{func.__name__} range_grid_mm={grid} does not contain {truth:.0f} mm" + ) diff --git a/tests/test_clubpose_fit_sequence_pinned_range.py b/tests/test_clubpose_fit_sequence_pinned_range.py new file mode 100644 index 00000000..df1b50c0 --- /dev/null +++ b/tests/test_clubpose_fit_sequence_pinned_range.py @@ -0,0 +1,33 @@ +from __future__ import annotations + +import numpy as np + +from openflight.camera.clubpose import fit + + +def test_fit_sequence_can_keep_a_singleton_range_hard_pinned(monkeypatch): + monkeypatch.setattr(fit, "CAMERA_CENTER_WORLD", np.zeros(3)) + monkeypatch.setattr(fit, "_ray_world", lambda _uv, _camera: np.asarray([1.0, 0.0, 0.0])) + monkeypatch.setattr( + fit, + "render_mask_6dof", + lambda _mesh, centre, _yaw, _pitch, _roll, _camera: np.asarray([[centre[0]]]), + ) + monkeypatch.setattr( + fit, + "iou", + lambda rendered, _observed: 1.0 - abs(float(rendered[0, 0]) - 1531.0) / 1000.0, + ) + + result = fit.fit_sequence( + object(), + {0: np.ones((10, 10), dtype=np.uint8)}, + object(), + range_grid_mm=(1581.0,), + yaw_grid=(0.0,), + pitch_grid=(0.0,), + roll_grid=(0.0,), + refine_range=False, + ) + + assert result[0]["range_mm"] == 1581.0 diff --git a/tests/test_clubpose_head_split.py b/tests/test_clubpose_head_split.py new file mode 100644 index 00000000..7c72d4e5 --- /dev/null +++ b/tests/test_clubpose_head_split.py @@ -0,0 +1,113 @@ +"""Tests for separating the clubhead from the shaft. + +The shipped tracker treats one moving connected component as the clubhead. That +held only because the reference capture was so overexposed that the shaft had no +contrast. On a properly exposed capture (session 20260825_181734) the shaft is a +strong moving object and merges with the head for the frames closest to impact, +so the merged component's centroid sits halfway up the shaft. + +These use synthetic masks with the geometry measured from that session, so they +run without the capture archive. +""" + +from __future__ import annotations + +import cv2 +import numpy as np +import pytest + +from openflight.camera.clubpose.head_split import ( + HEAD_CORE_PX, + SHAFT_MAX_PX, + SHAFT_REACH_PX, + split_head, +) + + +def _head(shape=(200, 320), centre=(150, 148), half=(22, 15)) -> np.ndarray: + """A compact body, the size a real clubhead measures at this plate scale.""" + m = np.zeros(shape, np.uint8) + cv2.ellipse(m, centre, half, 0, 0, 360, 1, -1) + return m + + +def _shaft(shape=(200, 320), start=(160, 138), end=(250, 20), width=5) -> np.ndarray: + """A thin line. Measured shaft width on real frames is 4-6 px.""" + m = np.zeros(shape, np.uint8) + cv2.line(m, start, end, 1, width) + return m + + +def test_isolated_head_is_returned_whole(): + """A component that is entirely clubhead must come back unchanged. + + An earlier version seeded the shaft marker by distance from the head core + alone, which cut the thin toe and sole extremities off 54 of 119 real heads. + """ + head = _head() + out = split_head(head) + assert out is not None + got, shaft = out + assert int(shaft.sum()) == 0 + assert int(got.sum()) == int(head.sum()) + + +def test_pure_shaft_is_refused(): + """Fail closed: a thin line is not a clubhead, so there is nothing to return.""" + assert split_head(_shaft()) is None + + +def test_merged_component_splits_at_the_neck(): + head, shaft = _head(), _shaft() + merged = np.clip(head + shaft, 0, 1).astype(np.uint8) + n, _ = cv2.connectedComponents(merged, 8) + assert n == 2, "fixture must actually be one connected component" + + out = split_head(merged) + assert out is not None + got_head, got_shaft = out + + assert not np.any(got_head & (merged == 0)) + assert not np.any(got_shaft & (merged == 0)) + assert int((got_head | got_shaft).sum()) == int(merged.sum()) + + overlap = int((got_head & head).sum()) / int(head.sum()) + assert overlap > 0.9, f"recovered only {overlap:.1%} of the head" + stolen = int((got_head & (shaft & ~head)).sum()) + assert stolen < 0.25 * int(shaft.sum()), "kept too much shaft" + + +def test_merged_centroid_is_wrong_which_is_why_this_exists(): + """The regression this guards: the merged centroid is not the clubhead.""" + head, shaft = _head(), _shaft() + merged = np.clip(head + shaft, 0, 1).astype(np.uint8) + + def centroid(m): + ys, xs = np.nonzero(m) + return np.array([xs.mean(), ys.mean()]) + + drift = float(np.linalg.norm(centroid(merged) - centroid(head))) + assert drift > 20.0, "fixture should reproduce the centroid drift" + + got_head, _ = split_head(merged) + fixed = float(np.linalg.norm(centroid(got_head) - centroid(head))) + # A shaft stub stays attached in this fixture (no sharp constriction); + # what matters is the centroid is the head's, not the merged blob's. + assert fixed < drift / 3.0, f"split centroid {fixed:.1f} px off vs merged {drift:.1f}" + + +def test_thresholds_sit_in_the_measured_gaps(): + """The constants are read off measurements, not tuned. Keep them there.""" + assert 3.0 < HEAD_CORE_PX < 5.0, "must exceed every shaft, sit below every head" + assert 3.0 < SHAFT_MAX_PX < 5.0 + assert 40.8 < SHAFT_REACH_PX < 145.4, "must sit in the measured reach gap" + + +@pytest.mark.parametrize("width", [4, 5, 6, 7]) +def test_split_survives_the_measured_shaft_width_range(width): + head = _head() + merged = np.clip(head + _shaft(width=width), 0, 1).astype(np.uint8) + out = split_head(merged) + assert out is not None + got_head, _ = out + assert int((got_head & head).sum()) / int(head.sum()) > 0.85 diff --git a/tests/test_clubpose_pose_scores.py b/tests/test_clubpose_pose_scores.py new file mode 100644 index 00000000..eca19f85 --- /dev/null +++ b/tests/test_clubpose_pose_scores.py @@ -0,0 +1,197 @@ +"""Do the pose scores rank a KNOWN-correct pose first? + +Every fit quality number this project has published came from silhouette IoU, +and section 11f measured IoU running *inversely* to pose correctness on real +segmented masks. Before replacing it we have to establish that the replacement +is implemented correctly, and the only place correctness is knowable is +synthetic data: render the mesh at a pose we chose, then ask each score to +recover it. + +Passing here does NOT mean a score works on real pixels. It means the score is +not broken. Those are different claims and this file only supports the second. +""" + +from __future__ import annotations + +import math + +import numpy as np +import pytest + +from openflight.camera.clubpose.fit import measured_camera, render_mask_6dof +from openflight.camera.clubpose.mesh import TriangleMesh +from openflight.camera.clubpose.projection import CAMERA_CENTER_WORLD, _ray_world +from openflight.camera.clubpose.scores import ( + Pose, + chamfer_px, + load_annotations, + mask_edge, + observed_omega_deg_s, + omega_residual_deg_s, + save_annotation, + score_pose, +) + +# World X is DOWNRANGE, Z is up -- so a centre must be built from a camera ray, +# exactly as fit_frame_6dof does, not written down as a bare triple. +_CAM = measured_camera(320, 200) +CENTER = tuple(CAMERA_CENTER_WORLD + _ray_world(np.array([_CAM.cx, _CAM.cy]), _CAM) * 1581.0) + + +def box_mesh() -> TriangleMesh: + """A clubhead-proportioned box: thin along the face normal, wide, medium tall. + + Deliberately not square in u/v so that yaw and pitch produce DIFFERENT + silhouette changes. A symmetric shape would hide a score that cannot tell + the two axes apart. + + Deliberately 4x a real clubhead. At true scale this box spans ~27 px, where + a 24 deg yaw changes its width by about one pixel -- which is the actual + measurement problem, and makes a poor unit test of arithmetic. Whether the + scores work at 27 px is what the landscape sweep measures; this file only + establishes that they are not broken. + """ + half = np.array([20.0, 180.0, 90.0]) + corners = np.array([[i, j, k] for i in (-1, 1) for j in (-1, 1) for k in (-1, 1)], float) + verts = corners * half + faces = np.array( + [ + [0, 1, 3], + [0, 3, 2], + [4, 6, 7], + [4, 7, 5], + [0, 4, 5], + [0, 5, 1], + [2, 3, 7], + [2, 7, 6], + [0, 2, 6], + [0, 6, 4], + [1, 5, 7], + [1, 7, 3], + ], + dtype=np.int32, + ) + return TriangleMesh(verts, faces, "unit_box", "b" * 64) + + +@pytest.fixture(name="cam") +def _cam(): + return measured_camera(320, 200) + + +@pytest.fixture(name="mesh") +def _mesh(): + return box_mesh() + + +def render(mesh, cam, yaw=0.0, pitch=0.0, roll=0.0): + m = render_mask_6dof(mesh, np.asarray(CENTER), yaw, pitch, roll, cam) + assert m is not None and m.sum() > 0, "fixture pose must project on-sensor" + return m + + +class TestChamfer: + def test_zero_for_identical_masks(self, mesh, cam): + m = render(mesh, cam) + assert chamfer_px(m, m) == pytest.approx(0.0, abs=1e-9) + + def test_symmetric_in_its_arguments(self, mesh, cam): + a, b = render(mesh, cam), render(mesh, cam, yaw=12.0) + assert chamfer_px(a, b) == pytest.approx(chamfer_px(b, a), abs=1e-9) + + def test_edge_is_a_boundary_not_the_body(self, mesh, cam): + m = render(mesh, cam) + e = mask_edge(m) + assert e.sum() < m.sum(), "an edge that is the whole mask is not an edge" + assert np.all(m[e]), "edge pixels must lie inside the mask" + + def test_empty_mask_is_reported_not_silently_zero(self, mesh, cam): + m = render(mesh, cam) + blank = np.zeros_like(m) + assert math.isinf(chamfer_px(m, blank)), "no overlap must fail closed, not score 0" + + +class TestRanksTruthOnSyntheticData: + """The mask IS the mesh at a known pose, so the true pose is recoverable.""" + + @pytest.mark.parametrize("err", [4.0, 8.0, 16.0]) + def test_chamfer_penalises_yaw_error(self, mesh, cam, err): + truth = render(mesh, cam) + assert chamfer_px(render(mesh, cam, yaw=err), truth) > chamfer_px(truth, truth) + + def test_chamfer_increases_monotonically_with_yaw_error(self, mesh, cam): + truth = render(mesh, cam) + d = [chamfer_px(render(mesh, cam, yaw=e), truth) for e in (0.0, 4.0, 8.0, 16.0, 24.0)] + assert d == sorted(d), f"chamfer must worsen as pose worsens, got {d}" + + def test_iou_also_ranks_truth_first_when_the_mask_is_exact(self, mesh, cam): + """IoU is not broken -- it is misleading on REAL masks (section 11f). + + Pinning that down here matters: it stops anyone reading the anti- + correlation result as 'IoU is buggy'. On an exact mask it behaves. + """ + truth = render(mesh, cam) + best = score_pose(mesh, Pose(CENTER, 0.0, 0.0, 0.0), truth, cam).iou + for err in (4.0, 8.0, 16.0): + assert score_pose(mesh, Pose(CENTER, err, 0.0, 0.0), truth, cam).iou < best + + def test_yaw_and_pitch_are_distinguishable(self, mesh, cam): + truth = render(mesh, cam) + assert chamfer_px(render(mesh, cam, yaw=10.0), truth) != pytest.approx( + chamfer_px(render(mesh, cam, pitch=10.0), truth), abs=1e-6 + ) + + +class TestOmegaResidual: + def test_observed_omega_recovers_a_known_rate(self): + dt, rate = 0.001, 2000.0 + a = Pose(CENTER, 0.0, 0.0, 0.0) + b = Pose(CENTER, rate * dt, 0.0, 0.0) + assert observed_omega_deg_s(a, b, dt) == pytest.approx(rate, rel=1e-6) + + def test_residual_is_zero_when_motion_matches_v_over_r(self): + speed, radius = 35.0, 1.6 + expected = math.degrees(speed / radius) + dt = 0.001 + a = Pose(CENTER, 0.0, 0.0, 0.0) + b = Pose(CENTER, expected * dt, 0.0, 0.0) + assert omega_residual_deg_s(a, b, dt, speed, radius) == pytest.approx(0.0, abs=1e-6) + + def test_residual_grows_when_the_pose_pair_rotates_too_fast(self): + speed, radius, dt = 35.0, 1.6, 0.001 + expected = math.degrees(speed / radius) + a = Pose(CENTER, 0.0, 0.0, 0.0) + slow = Pose(CENTER, expected * dt, 0.0, 0.0) + fast = Pose(CENTER, expected * dt * 3.0, 0.0, 0.0) + assert omega_residual_deg_s(a, fast, dt, speed, radius) > omega_residual_deg_s( + a, slow, dt, speed, radius + ) + + def test_rejects_nonphysical_inputs(self): + a = Pose(CENTER, 0.0, 0.0, 0.0) + b = Pose(CENTER, 1.0, 0.0, 0.0) + with pytest.raises(ValueError): + omega_residual_deg_s(a, b, 0.0, 35.0, 1.6) + with pytest.raises(ValueError): + omega_residual_deg_s(a, b, 0.001, -1.0, 1.6) + + +class TestAnnotationStorage: + def test_roundtrips_a_labelled_pose(self, tmp_path): + path = tmp_path / "labels.jsonl" + pose = Pose(CENTER, 3.0, -7.5, 12.0) + save_annotation(path, shot=2, frame=64, pose=pose, annotator="reviewer", pass_index=1) + got = load_annotations(path) + assert len(got) == 1 + assert got[0]["shot"] == 2 and got[0]["frame"] == 64 + assert got[0]["pose"]["yaw_deg"] == pytest.approx(3.0) + assert Pose.from_dict(got[0]["pose"]).roll_deg == pytest.approx(12.0) + + def test_appends_rather_than_overwrites(self, tmp_path): + """A second labelling pass must never destroy the first -- the whole + experiment is the COMPARISON between passes.""" + path = tmp_path / "labels.jsonl" + p = Pose(CENTER, 0.0, 0.0, 0.0) + save_annotation(path, shot=2, frame=64, pose=p, annotator="a", pass_index=1) + save_annotation(path, shot=2, frame=64, pose=p, annotator="a", pass_index=2) + assert len(load_annotations(path)) == 2 diff --git a/tests/test_clubpose_rigid_motion.py b/tests/test_clubpose_rigid_motion.py new file mode 100644 index 00000000..81cfed61 --- /dev/null +++ b/tests/test_clubpose_rigid_motion.py @@ -0,0 +1,46 @@ +"""Physical constraints for sequence-level clubhead rotation.""" + +from __future__ import annotations + +import math + +import numpy as np +import pytest + +from openflight.camera.clubpose.motion import axis_from_angles, constrained_omega_deg_s + + +@pytest.mark.parametrize( + ("azimuth_deg", "elevation_deg", "expected"), + ( + (0.0, 0.0, (1.0, 0.0, 0.0)), + (90.0, 0.0, (0.0, 1.0, 0.0)), + (0.0, 90.0, (0.0, 0.0, 1.0)), + ), +) +def test_axis_angles_form_a_unit_world_axis(azimuth_deg, elevation_deg, expected): + axis = axis_from_angles(azimuth_deg, elevation_deg) + np.testing.assert_allclose(axis, expected, atol=1e-12) + assert np.linalg.norm(axis) == pytest.approx(1.0) + + +def test_constrained_omega_has_v_over_r_magnitude(): + omega = constrained_omega_deg_s( + speed_mps=36.6, + swing_radius_m=1.6, + axis_azimuth_deg=30.0, + axis_elevation_deg=-20.0, + ) + + assert np.linalg.norm(omega) == pytest.approx(math.degrees(36.6 / 1.6)) + np.testing.assert_allclose( + omega / np.linalg.norm(omega), + axis_from_angles(30.0, -20.0), + atol=1e-12, + ) + + +@pytest.mark.parametrize("speed_mps,swing_radius_m", ((0.0, 1.6), (36.6, 0.0), (-1.0, 1.6))) +def test_constrained_omega_rejects_non_physical_inputs(speed_mps, swing_radius_m): + with pytest.raises(ValueError): + constrained_omega_deg_s(speed_mps, swing_radius_m, 0.0, 0.0)