From 86a7e5b6327dbfbcc60d82d2a9ed77ea8fd93a46 Mon Sep 17 00:00:00 2001
From: EricSpencer00 <espencer2@luc.edu>
Date: Tue, 22 Sep 2026 05:55:00 -0500
Subject: [PATCH] perf: avoid full-frame promotion during alignment

---
 src/faciometry/pipeline/align.py | 14 ++++++++---
 tests/unit/test_align.py         | 43 ++++++++++++++++++++++++++++++++
 2 files changed, 54 insertions(+), 3 deletions(-)
 create mode 100644 tests/unit/test_align.py

diff --git a/src/faciometry/pipeline/align.py b/src/faciometry/pipeline/align.py
index 6e5d7f8..3087e56 100644
--- a/src/faciometry/pipeline/align.py
+++ b/src/faciometry/pipeline/align.py
@@ -190,9 +190,17 @@ def _bilinear_sample(image: Image, xs: NDArray[np.float64], ys: NDArray[np.float
     x1c = np.clip(x0 + 1, 0, w - 1)
     y0c = np.clip(y0, 0, h - 1)
     y1c = np.clip(y0 + 1, 0, h - 1)
-    img = image.astype(np.float64)
-    top = img[y0c, x0c] * (1.0 - fx) + img[y0c, x1c] * fx
-    bot = img[y1c, x0c] * (1.0 - fx) + img[y1c, x1c] * fx
+    # Keep the source image in its compact dtype until after indexing. The
+    # crop reads only four source pixels per output pixel; promoting the whole
+    # photograph to float64 here copies tens of megabytes that are never read.
+    # Each gathered neighbourhood is promoted before the same float64
+    # arithmetic, so interpolation and rounding remain byte-for-byte stable.
+    top = image[y0c, x0c].astype(np.float64) * (1.0 - fx) + image[
+        y0c, x1c
+    ].astype(np.float64) * fx
+    bot = image[y1c, x0c].astype(np.float64) * (1.0 - fx) + image[
+        y1c, x1c
+    ].astype(np.float64) * fx
     out = top * (1.0 - fy) + bot * fy
     return np.clip(np.rint(out), 0, 255).astype(np.uint8)
 
diff --git a/tests/unit/test_align.py b/tests/unit/test_align.py
new file mode 100644
index 0000000..09d5834
--- /dev/null
+++ b/tests/unit/test_align.py
@@ -0,0 +1,43 @@
+"""The canonical crop's interpolation must retain its pixel semantics."""
+
+from __future__ import annotations
+
+import numpy as np
+
+from faciometry.pipeline.align import _bilinear_sample
+
+
+def _full_frame_reference(
+    image: np.ndarray, xs: np.ndarray, ys: np.ndarray
+) -> np.ndarray:
+    """The pre-optimisation implementation, kept as an exact-output oracle."""
+    h, w = image.shape[:2]
+    x0 = np.floor(xs).astype(np.int64)
+    y0 = np.floor(ys).astype(np.int64)
+    fx = (xs - x0)[..., None]
+    fy = (ys - y0)[..., None]
+    x0c = np.clip(x0, 0, w - 1)
+    x1c = np.clip(x0 + 1, 0, w - 1)
+    y0c = np.clip(y0, 0, h - 1)
+    y1c = np.clip(y0 + 1, 0, h - 1)
+    full = image.astype(np.float64)
+    top = full[y0c, x0c] * (1.0 - fx) + full[y0c, x1c] * fx
+    bot = full[y1c, x0c] * (1.0 - fx) + full[y1c, x1c] * fx
+    return np.clip(np.rint(top * (1.0 - fy) + bot * fy), 0, 255).astype(np.uint8)
+
+
+def test_sampled_promotion_matches_full_frame_reference_at_image_edges():
+    rng = np.random.default_rng(20260922)
+    image = rng.integers(0, 256, size=(7, 9, 3), dtype=np.uint8)
+    xs = np.array(
+        [[-2.25, 0.0, 0.125, 7.5, 8.0, 8.75, 11.0]],
+        dtype=float,
+    )
+    ys = np.array(
+        [[-1.5, 0.0, 2.75, 5.125, 6.0, 6.5, 9.0]],
+        dtype=float,
+    )
+
+    np.testing.assert_array_equal(
+        _bilinear_sample(image, xs, ys), _full_frame_reference(image, xs, ys)
+    )