From bbb4510effdae2834eb50042655044ee445fbfda Mon Sep 17 00:00:00 2001
From: Constantinos Eleftheriou <constantinos.eleftheriou@ed.ac.uk>
Date: Thu, 13 Aug 2026 15:02:32 +0100
Subject: [PATCH 01/11] fix inconsistent reading of landmark x,y coordinates

---
 src/mesoscopy/io.py                           |   8 +-
 src/mesoscopy/register/__init__.py            |   3 +-
 src/mesoscopy/register/landmarks_gui.py       |  36 ++++--
 src/mesoscopy/register/qa.py                  |  26 +++-
 src/mesoscopy/register/transform.py           |   7 +-
 ...eso_date-2025-05-15T16_29_37_landmarks.csv |  10 ++
 tests/test_registration.py                    | 121 ++++++++++++++++++
 7 files changed, 187 insertions(+), 24 deletions(-)
 create mode 100644 sub-TT712_exp-10hrb25ce_meso_date-2025-05-15T16_29_37_landmarks.csv

diff --git a/src/mesoscopy/io.py b/src/mesoscopy/io.py
index 9319f37..592ffc0 100644
--- a/src/mesoscopy/io.py
+++ b/src/mesoscopy/io.py
@@ -188,11 +188,13 @@ def load_deltaf(path: str, nwb: bool = False) -> tuple[str, np.ndarray, np.ndarr
 def read_points(path: str) -> dict[str, tuple[float, float]]:
     """Read a landmark points file.
 
+    Coordinates are read as (x, y) tuples, i.e. (column, row).
+
     Args:
         path (str): Path to the points file.
 
     Returns:
-        dict[str, tuple[float, float]]: Dictionary with the landmark names as keys and their x-y coordinates
+        dict[str, tuple[float, float]]: Dictionary with the landmark names as keys and their (x, y) coordinates
 
     Raises:
         ValueError: If the file format is unsupported.
@@ -293,9 +295,11 @@ def _read_csv_points(path: str) -> dict[str, tuple[float, float]]:
 def write_points(path: str, points: dict[str, tuple[float, float]]) -> None:
     """Write a dictionary of landmark points to a CSV file.
 
+    Coordinates are written as (x, y) tuples, i.e. (column, row).
+
     Args:
         path (str): Path to output CSV file.
-        points (dict[str, tuple[float, float]]): Dictionary with the landmark names as keys and their x-y coordinates
+        points (dict[str, tuple[float, float]]): Dictionary with the landmark names as keys and their (x, y) coordinates
     """
     if not path.endswith(".csv"):
         path += ".csv"
diff --git a/src/mesoscopy/register/__init__.py b/src/mesoscopy/register/__init__.py
index e68cbc8..31b4dab 100644
--- a/src/mesoscopy/register/__init__.py
+++ b/src/mesoscopy/register/__init__.py
@@ -77,7 +77,8 @@ def label_cmd(path, out_dir, template_points, session_id) -> dict:
 
     Returns:
         dict: Dictionary with the landmarks and their x-y coordinates.
-              Dictionary keys are landmark names, while x-y coordinates are stored as an (y, x)
+              Dictionary keys are landmark names, while x-y coordinates are stored as an (x, y)
+              tuple, i.e. (column, row).
     """
     click.echo("Loading imaging data...")
     nwb = bool(path.endswith(".nwb"))
diff --git a/src/mesoscopy/register/landmarks_gui.py b/src/mesoscopy/register/landmarks_gui.py
index 61972d6..10ec9c3 100644
--- a/src/mesoscopy/register/landmarks_gui.py
+++ b/src/mesoscopy/register/landmarks_gui.py
@@ -55,13 +55,17 @@ def mark_landmarks(
     - rpRSP: Right posterior aspect of the retrosplenial cortex.
     - aIPB: Anterior aspect of the interparietal bone.
 
+    Landmark coordinates are stored as (x, y) tuples, i.e. (column, row), matching the convention used
+    by the template landmark files and by ``skimage.transform``. Napari works in (row, column) order, so
+    coordinates are transposed on the way into and out of the viewer.
+
     Args:
         maxip_image (npt.NDArray): Maximum intensity projection image. Could be either channel.
         alt_image (npt.NDArray): Alternative image to be displayed alongside the maximum intensity projection image. Usually a second channel.
-        template_landmarks (dict, optional): Dictionary with the landmarks and their x-y coordinates. Dictionary keys are landmark names, while x-y coordinates are stored as an (y, x) tuple. Defaults to {}.
+        template_landmarks (dict, optional): Dictionary with the landmarks and their x-y coordinates, used to seed the initial point positions. Dictionary keys are landmark names, while x-y coordinates are stored as an (x, y) tuple. Defaults to {}.
 
     Returns:
-        dict: Dictionary with the landmarks and their x-y coordinates. Dictionary keys are landmark names, while x-y coordinates are stored as an (y, x) tuple.
+        dict: Dictionary with the landmarks and their x-y coordinates. Dictionary keys are landmark names, while x-y coordinates are stored as an (x, y) tuple.
 
     Raises:
         ValueError: If the maximum intensity projection and alternative image do not have the same dimensions.
@@ -85,21 +89,24 @@ def mark_landmarks(
         viewer.add_image(alt_image, name="alt_maxip")
 
     default_landmark_locations = template_landmarks or {
-        "bregma": (maxip_height / 2, maxip_width / 2),
-        "cFP": (maxip_height / 7, maxip_width / 2),
-        "rFP": (maxip_height / 7, maxip_width / 1.5),
-        "lFP": (maxip_height / 7, maxip_width / 3),
-        "rPB": (maxip_height / 4, maxip_width / 1.25),
-        "lPB": (maxip_height / 4, maxip_width / 5),
-        "lpRSP": (maxip_height / 1.25, maxip_width / 2.25),
-        "rpRSP": (maxip_height / 1.25, maxip_width / 1.75),
-        "aIPB": (maxip_height / 1.4, maxip_width / 2),
+        "bregma": (maxip_width / 2, maxip_height / 2),
+        "cFP": (maxip_width / 2, maxip_height / 7),
+        "rFP": (maxip_width / 1.5, maxip_height / 7),
+        "lFP": (maxip_width / 3, maxip_height / 7),
+        "rPB": (maxip_width / 1.25, maxip_height / 4),
+        "lPB": (maxip_width / 5, maxip_height / 4),
+        "lpRSP": (maxip_width / 2.25, maxip_height / 1.25),
+        "rpRSP": (maxip_width / 1.75, maxip_height / 1.25),
+        "aIPB": (maxip_width / 2, maxip_height / 1.4),
     }
 
     landmarks = list(default_landmark_locations.keys())
 
+    # Landmarks are stored as (x, y), napari points are (row, column) - transpose on the way in.
+    seed_points = np.array([(y, x) for x, y in default_landmark_locations.values()], dtype=float)
+
     points_layer = viewer.add_points(
-        data=np.array(list(default_landmark_locations.values())),
+        data=seed_points,
         name="landmarks",
         ndim=2,
         properties={"label": landmarks},
@@ -120,7 +127,10 @@ def mark_landmarks(
 
     napari.run()
 
-    return OrderedDict(zip(landmarks, points_layer.data, strict=True))
+    # Transpose back from napari's (row, column) to the stored (x, y) convention.
+    return OrderedDict(
+        (landmark, (float(x), float(y))) for landmark, (y, x) in zip(landmarks, points_layer.data, strict=True)
+    )
 
 
 def _create_label_menu(points_layer, labels):
diff --git a/src/mesoscopy/register/qa.py b/src/mesoscopy/register/qa.py
index 891d190..d24d248 100644
--- a/src/mesoscopy/register/qa.py
+++ b/src/mesoscopy/register/qa.py
@@ -20,6 +20,7 @@
 #  SOFTWARE.
 """Module for quality assurance (QA) functions in the mesoscopy registration pipeline."""
 
+import numpy as np
 import numpy.typing as npt
 import plotly.express as px
 import plotly.graph_objects as go
@@ -30,6 +31,16 @@ def plot_landmarks(
     target_landmarks: npt.NDArray,
     as_html: bool = False,
 ) -> str | go.Figure:
+    """Plots two sets of landmark points against each other using Plotly.
+
+    Args:
+        source_landmarks (npt.NDArray): Source landmarks as an (n, 2) array of (x, y) coordinates.
+        target_landmarks (npt.NDArray): Target landmarks as an (n, 2) array of (x, y) coordinates.
+        as_html (bool, optional): If True, returns the plot as an HTML string. If False, returns a Plotly Figure object.
+
+    Returns:
+        str | go.Figure: The plot as an HTML string if `as_html` is True, otherwise a Plotly Figure object.
+    """
     fig = go.Figure(
         layout=go.Layout(
             margin={"l": 20, "r": 20, "t": 20, "b": 20},
@@ -40,10 +51,10 @@ def plot_landmarks(
     )
 
     fig.add_trace(
-        go.Scatter(x=source_landmarks[:, 1], y=source_landmarks[:, 0], mode="markers", name="Source landmarks")
+        go.Scatter(x=source_landmarks[:, 0], y=source_landmarks[:, 1], mode="markers", name="Source landmarks")
     )
     fig.add_trace(
-        go.Scatter(x=target_landmarks[:, 1], y=target_landmarks[:, 0], mode="markers", name="Target landmarks")
+        go.Scatter(x=target_landmarks[:, 0], y=target_landmarks[:, 1], mode="markers", name="Target landmarks")
     )
 
     if as_html:
@@ -51,11 +62,13 @@ def plot_landmarks(
     return fig
 
 
-def plot_frame(data: npt.NDArray, landmarks=npt.NDArray | None, as_html: bool = False) -> str | go.Figure:
+def plot_frame(data: npt.NDArray, landmarks: npt.NDArray | None = None, as_html: bool = False) -> str | go.Figure:
     """Plots a single frame of data using Plotly.
 
     Args:
         data (npt.NDArray): The frame data to plot, as a NumPy array.
+        landmarks (npt.NDArray, optional): Landmarks to overlay on the frame, as an (n, 2) array of (x, y)
+            coordinates. Defaults to None.
         as_html (bool, optional): If True, returns the plot as an HTML string. If False, returns a Plotly Figure object.
 
     Returns:
@@ -68,11 +81,12 @@ def plot_frame(data: npt.NDArray, landmarks=npt.NDArray | None, as_html: bool =
         }
     )
 
-    if landmarks:
+    if landmarks is not None:
+        landmarks = np.asarray(landmarks)
         fig.add_trace(
             go.Scatter(
-                x=landmarks[:, 1],
-                y=landmarks[:, 0],
+                x=landmarks[:, 0],
+                y=landmarks[:, 1],
                 mode="markers",
                 name="Template landmarks",
                 marker={"color": "Green"},
diff --git a/src/mesoscopy/register/transform.py b/src/mesoscopy/register/transform.py
index 4b853a5..ab76c88 100644
--- a/src/mesoscopy/register/transform.py
+++ b/src/mesoscopy/register/transform.py
@@ -37,10 +37,13 @@ def landmarks_affine(
 ) -> tuple[np.ndarray, trf.ProjectiveTransform]:
     """Warp a DeltaF/F series to match a template using anatomical landmarks.
 
+    Both landmark sets must use the same (x, y) — i.e. (column, row) — coordinate convention as
+    ``skimage.transform``.
+
     Args:
         deltaf_series (np.ndarray): DeltaF/F series.
-        recording_landmarks (dict): Recording landmarks.
-        template_landmarks (dict): Template landmarks.
+        recording_landmarks (dict): Recording landmarks, as {name: (x, y)}.
+        template_landmarks (dict): Template landmarks, as {name: (x, y)}.
         crop_x (int, optional): Crop x-axis. Defaults to 0.
         crop_y (int, optional): Crop y-axis. Defaults to 0.
 
diff --git a/sub-TT712_exp-10hrb25ce_meso_date-2025-05-15T16_29_37_landmarks.csv b/sub-TT712_exp-10hrb25ce_meso_date-2025-05-15T16_29_37_landmarks.csv
new file mode 100644
index 0000000..14c2354
--- /dev/null
+++ b/sub-TT712_exp-10hrb25ce_meso_date-2025-05-15T16_29_37_landmarks.csv
@@ -0,0 +1,10 @@
+landmark,x,y
+bregma,71.0,63.24942780179567
+cFP,72.08314260059856,17.91685739940145
+lFP,47.028477131138644,17.19476233233575
+rFP,96.13780807005845,16.833714798802895
+lPB,30.0,31.75057219820434
+rPB,113.1662852011971,36.444190134131404
+lpRSP,60.722095067065695,116.44419013413142
+rpRSP,81.0,117.1662852011971
+aIPB,69.91685739940145,103.08314260059856
diff --git a/tests/test_registration.py b/tests/test_registration.py
index dfd243a..1ca48ed 100644
--- a/tests/test_registration.py
+++ b/tests/test_registration.py
@@ -7,7 +7,10 @@
 
 import mesoscopy
 import mesoscopy.register as reg
+import mesoscopy.register.landmarks_gui as reg_gui
+import mesoscopy.register.qa as reg_qa
 import mesoscopy.resources as res
+from mesoscopy import io
 from mesoscopy.register.transform import landmarks_affine
 
 
@@ -134,3 +137,121 @@ def test_register_landmarks_cli_default_points_h5(preproc_h5, output_dir):
 
 
 def test_mark_landmarks_gui(): ...
+
+
+# ---------------------------------------------------------------------------
+# Landmark coordinate convention
+#
+# Landmarks are stored as (x, y) == (column, row) everywhere on disk, matching
+# skimage.transform. Napari is the only (row, column) surface, and mark_landmarks
+# is responsible for transposing on the way in and out.
+# ---------------------------------------------------------------------------
+
+
+def test_template_landmarks_are_xy_not_rowcol():
+    """Paired left/right template landmarks must hit the same atlas region in each hemisphere.
+
+    This only holds when the shipped landmark file is read as (x, y); reading it as
+    (row, column) puts every pair in the wrong place.
+    """
+    landmarks = res.get_default_landmarks()
+    left_aba, right_aba = res.get_atlas()
+    midline = left_aba.shape[1] / 2
+
+    for left_name, right_name in (("lFP", "rFP"), ("lPB", "rPB"), ("lpRSP", "rpRSP")):
+        left_x, left_y = landmarks[left_name]
+        right_x, right_y = landmarks[right_name]
+
+        # Lateral landmarks sit on their own side of the midline...
+        assert left_x < midline
+        assert right_x > midline
+
+        # ...and pick out the same labelled region in the matching hemisphere map.
+        left_region = left_aba[int(left_y), int(left_x)]
+        right_region = right_aba[int(right_y), int(right_x)]
+        assert left_region != 0
+        assert left_region == right_region
+
+    # Midline landmarks sit on the midline.
+    for name in ("bregma", "cFP", "aIPB"):
+        assert landmarks[name][0] == pytest.approx(midline, abs=1.0)
+
+
+def test_points_roundtrip_preserves_xy(tmp_path):
+    """write_points -> read_points must not transpose the coordinates."""
+    points = {"bregma": (71.0, 60.0), "lFP": (51.0, 19.0)}
+    path = str(tmp_path / "points.csv")
+
+    io.write_points(path, points)
+
+    assert io.read_points(path) == points
+
+
+def test_mark_landmarks_transposes_between_napari_and_storage(monkeypatch):
+    """Seed points enter napari as (row, col); marked points come back out as (x, y)."""
+    captured = {}
+
+    class FakePointsLayer:
+        def __init__(self, data):
+            self.data = data
+            self.mode = None
+            self.face_color_mode = None
+
+    class FakeWindow:
+        def add_dock_widget(self, widget):
+            pass
+
+    class FakeViewer:
+        def __init__(self):
+            self.window = FakeWindow()
+
+        def add_image(self, *args, **kwargs):
+            pass
+
+        def add_points(self, data, **kwargs):
+            captured["seed"] = np.asarray(data, dtype=float)
+            # Simulate the user dragging the first point by +3 rows / +5 columns.
+            marked = np.asarray(data, dtype=float).copy()
+            marked[0] += (3.0, 5.0)
+            return FakePointsLayer(marked)
+
+    monkeypatch.setattr(reg_gui.napari, "view_image", lambda *args, **kwargs: FakeViewer())
+    monkeypatch.setattr(reg_gui.napari, "run", lambda *args, **kwargs: None)
+    monkeypatch.setattr(reg_gui, "_create_label_menu", lambda points_layer, labels: None)
+
+    template = {"a": (10.0, 20.0), "b": (30.0, 40.0)}
+    marked = reg_gui.mark_landmarks(np.zeros((100, 120)), None, template)
+
+    # (x, y) seeds are handed to napari as (row, col).
+    np.testing.assert_array_equal(captured["seed"], [[20.0, 10.0], [40.0, 30.0]])
+
+    # Marked points come back as (x, y), so the +3 row / +5 col drag reads as +5 x / +3 y.
+    assert marked["a"] == (15.0, 23.0)
+    assert marked["b"] == (30.0, 40.0)
+
+
+def test_qa_plot_landmarks_uses_xy():
+    """The QA scatter must plot x on the x-axis, not the row index."""
+    source = np.array([[10.0, 20.0], [30.0, 40.0]])
+    target = np.array([[11.0, 21.0], [31.0, 41.0]])
+
+    fig = reg_qa.plot_landmarks(source, target)
+
+    np.testing.assert_array_equal(fig.data[0].x, source[:, 0])
+    np.testing.assert_array_equal(fig.data[0].y, source[:, 1])
+    np.testing.assert_array_equal(fig.data[1].x, target[:, 0])
+    np.testing.assert_array_equal(fig.data[1].y, target[:, 1])
+
+
+def test_qa_plot_frame_without_landmarks():
+    """plot_frame must be callable without landmarks."""
+    fig = reg_qa.plot_frame(np.random.default_rng(0).random((10, 12)))
+    assert len(fig.data) == 1
+
+
+def test_qa_plot_frame_landmarks_use_xy():
+    landmarks = np.array([[10.0, 20.0], [30.0, 40.0]])
+    fig = reg_qa.plot_frame(np.random.default_rng(0).random((50, 50)), landmarks=landmarks)
+
+    np.testing.assert_array_equal(fig.data[1].x, landmarks[:, 0])
+    np.testing.assert_array_equal(fig.data[1].y, landmarks[:, 1])

From 0fe48ee35b94a8b32c4c552f13d6c88e58497b3e Mon Sep 17 00:00:00 2001
From: Constantinos Eleftheriou <constantinos.eleftheriou@ed.ac.uk>
Date: Thu, 13 Aug 2026 15:02:55 +0100
Subject: [PATCH 02/11] remove accidentally committed csv

---
 ...hrb25ce_meso_date-2025-05-15T16_29_37_landmarks.csv | 10 ----------
 1 file changed, 10 deletions(-)
 delete mode 100644 sub-TT712_exp-10hrb25ce_meso_date-2025-05-15T16_29_37_landmarks.csv

diff --git a/sub-TT712_exp-10hrb25ce_meso_date-2025-05-15T16_29_37_landmarks.csv b/sub-TT712_exp-10hrb25ce_meso_date-2025-05-15T16_29_37_landmarks.csv
deleted file mode 100644
index 14c2354..0000000
--- a/sub-TT712_exp-10hrb25ce_meso_date-2025-05-15T16_29_37_landmarks.csv
+++ /dev/null
@@ -1,10 +0,0 @@
-landmark,x,y
-bregma,71.0,63.24942780179567
-cFP,72.08314260059856,17.91685739940145
-lFP,47.028477131138644,17.19476233233575
-rFP,96.13780807005845,16.833714798802895
-lPB,30.0,31.75057219820434
-rPB,113.1662852011971,36.444190134131404
-lpRSP,60.722095067065695,116.44419013413142
-rpRSP,81.0,117.1662852011971
-aIPB,69.91685739940145,103.08314260059856

From 1cd9d9bf92080e35d320bc0a2db03207c05064f1 Mon Sep 17 00:00:00 2001
From: Constantinos Eleftheriou <constantinos.eleftheriou@ed.ac.uk>
Date: Thu, 13 Aug 2026 15:19:59 +0100
Subject: [PATCH 03/11] replace crop-x/y with output-width/height options,
 ensure registered output matches ABA or predefined shape

---
 src/mesoscopy/register/__init__.py  | 31 +++++++---
 src/mesoscopy/register/transform.py | 23 +++++---
 tests/test_registration.py          | 88 +++++++++++++++++++++--------
 3 files changed, 104 insertions(+), 38 deletions(-)

diff --git a/src/mesoscopy/register/__init__.py b/src/mesoscopy/register/__init__.py
index 31b4dab..2c0365d 100644
--- a/src/mesoscopy/register/__init__.py
+++ b/src/mesoscopy/register/__init__.py
@@ -142,15 +142,25 @@ def label_cmd(path, out_dir, template_points, session_id) -> dict:
     type=click.Path(dir_okay=False),
     help="Path to template landmark points in Fiji XML points format",
 )
-@click.option("--crop-x", default=0, help="Crop recording along the x-axis.")
-@click.option("--crop-y", default=0, help="Crop recording along the y-axis.")
+@click.option(
+    "--output-width",
+    type=int,
+    default=None,
+    help="Width of the registered frames. Defaults to the width of the Allen CCF template.",
+)
+@click.option(
+    "--output-height",
+    type=int,
+    default=None,
+    help="Height of the registered frames. Defaults to the height of the Allen CCF template.",
+)
 def landmarks_cmd(
     path: str,
     out_dir: str,
     recording_points: str,
     template_points: str,
-    crop_x: int = 0,
-    crop_y: int = 0,
+    output_width: int | None = None,
+    output_height: int | None = None,
 ) -> str:
     """Register a recording to a template based on defined landmarks.
 
@@ -159,8 +169,8 @@ def landmarks_cmd(
         out_dir (str): Output directory for registered recording.
         recording_points (str, optional): Path to recording landmark points in CSV or Fiji XML points format.
         template_points (str, optional): Path to template landmark points in CSV or Fiji XML points format.
-        crop_x (int, optional): Number of pixels to crop from the x-axis of the recording. Defaults to 0.
-        crop_y (int, optional): Number of pixels to crop from the y-axis of the recording. Defaults to 0.
+        output_width (int, optional): Width of the registered frames. Defaults to the Allen CCF template width.
+        output_height (int, optional): Height of the registered frames. Defaults to the Allen CCF template height.
 
     Returns:
         str: Path to the registered recording file.
@@ -194,12 +204,17 @@ def landmarks_cmd(
             raise ValueError(msg)
     recording_landmarks = io.read_points(recording_points)
 
+    # Registered frames land in template space, so default their shape to that of the CCF atlas.
+    output_shape = None
+    if output_width or output_height:
+        atlas_height, atlas_width = res.get_atlas()[0].shape
+        output_shape = (output_height or atlas_height, output_width or atlas_width)
+
     warped, tform = trf.landmarks_affine(
         deltaf_series,
         recording_landmarks,
         template_landmarks,
-        crop_x=crop_x,
-        crop_y=crop_y,
+        output_shape=output_shape,
     )
 
     # Save warped frames and timestamps
diff --git a/src/mesoscopy/register/transform.py b/src/mesoscopy/register/transform.py
index ab76c88..66318e2 100644
--- a/src/mesoscopy/register/transform.py
+++ b/src/mesoscopy/register/transform.py
@@ -27,30 +27,36 @@
 import numpy as np
 from skimage import transform as trf
 
+import mesoscopy.resources as res
+
 
 def landmarks_affine(
     deltaf_series: np.ndarray,
     recording_landmarks: dict,
     template_landmarks: dict,
-    crop_x: int = 0,
-    crop_y: int = 0,
+    output_shape: tuple[int, int] | None = None,
 ) -> tuple[np.ndarray, trf.ProjectiveTransform]:
     """Warp a DeltaF/F series to match a template using anatomical landmarks.
 
     Both landmark sets must use the same (x, y) — i.e. (column, row) — coordinate convention as
     ``skimage.transform``.
 
+    The registered frames are in template space, so their shape is that of the template rather than
+    that of the recording. It defaults to the shape of the Allen CCF atlas, but can be overridden with ``output_shape``.
+
     Args:
         deltaf_series (np.ndarray): DeltaF/F series.
         recording_landmarks (dict): Recording landmarks, as {name: (x, y)}.
         template_landmarks (dict): Template landmarks, as {name: (x, y)}.
-        crop_x (int, optional): Crop x-axis. Defaults to 0.
-        crop_y (int, optional): Crop y-axis. Defaults to 0.
+        output_shape (tuple[int, int], optional): Shape of the registered frames, as (height, width).
+            Defaults to the shape of the Allen CCF atlas template.
 
     Returns:
         tuple[np.ndarray, trf.ProjectiveTransform]: Registered DeltaF/F series and affine
             transformation matrix.
     """
+    if output_shape is None:
+        output_shape = res.get_atlas()[0].shape
     if not isinstance(deltaf_series, np.ndarray):
         click.echo("Loading imaging data into memory...")
         deltaf_series = np.asarray(deltaf_series)
@@ -67,15 +73,16 @@ def landmarks_affine(
     n_frames = deltaf_series.shape[0]
 
     def _warp_frame(idx: int) -> np.ndarray:
-        if crop_x > 0 or crop_y > 0:
-            return trf.warp(deltaf_series[idx, :crop_y, :crop_x], tform, order=3)
-        return trf.warp(deltaf_series[idx], tform, order=3)
+        return trf.warp(deltaf_series[idx], tform, order=3, output_shape=output_shape)
 
     start = time.time()
     results: list[np.ndarray | None] = [None] * n_frames
     n_workers = os.cpu_count() or 1
     with (
-        click.progressbar(length=n_frames, label="Registering recording to template...") as bar,
+        click.progressbar(
+            length=n_frames,
+            label=f"Registering recording to template ({output_shape[1]}x{output_shape[0]} frames)...",
+        ) as bar,
         ThreadPoolExecutor(max_workers=n_workers) as executor,
     ):
         futures = {executor.submit(_warp_frame, i): i for i in range(n_frames)}
diff --git a/tests/test_registration.py b/tests/test_registration.py
index 1ca48ed..acac0f2 100644
--- a/tests/test_registration.py
+++ b/tests/test_registration.py
@@ -37,11 +37,17 @@ def landmark_pair():
 # Tests
 # ---------------------------------------------------------------------------
 
-def test_landmarks_affine_output_shape(small_series, landmark_pair):
-    """Warped series must have the same shape as the input."""
+def test_landmarks_affine_defaults_to_atlas_shape(small_series, landmark_pair):
+    """Registered frames are in template space, so they default to the atlas shape.
+
+    This is the shape mesoscopy.process.region requires; the shape of the recording is irrelevant.
+    """
     recording_lm, template_lm = landmark_pair
     warped, _ = landmarks_affine(small_series, recording_lm, template_lm)
-    assert warped.shape == small_series.shape
+
+    atlas_shape = res.get_atlas()[0].shape
+    assert atlas_shape != small_series.shape[1:]  # guard: the fixture must not already be atlas-shaped
+    assert warped.shape == (small_series.shape[0], *atlas_shape)
 
 
 def test_landmarks_affine_returns_projective_transform(small_series, landmark_pair):
@@ -70,46 +76,50 @@ def test_landmarks_affine_matches_sequential_reference(small_series, landmark_pa
     template  = np.array(list(template_lm.values()),  dtype=np.float32)
     recording = np.array(list(recording_lm.values()), dtype=np.float32)
     ref_tform = trf.estimate_transform("affine", template, recording)
-    reference = np.array([trf.warp(small_series[i], ref_tform, order=3) for i in range(small_series.shape[0])])
+    reference = np.array(
+        [trf.warp(small_series[i], ref_tform, order=3, output_shape=(40, 40)) for i in range(small_series.shape[0])]
+    )
 
-    warped, _ = landmarks_affine(small_series, recording_lm, template_lm)
+    warped, _ = landmarks_affine(small_series, recording_lm, template_lm, output_shape=(40, 40))
 
     np.testing.assert_array_equal(warped, reference)
 
 
-def test_landmarks_affine_crop(landmark_pair):
-    """With crop_x / crop_y the warped frames should be cropped to those dimensions."""
+def test_landmarks_affine_explicit_output_shape(landmark_pair):
+    """An explicit output_shape overrides the atlas default."""
     rng = np.random.default_rng(1)
     series = rng.random((5, 40, 40), dtype=np.float32)
     recording_lm, template_lm = landmark_pair
-    crop_x, crop_y = 30, 25
 
-    warped, _ = landmarks_affine(series, recording_lm, template_lm, crop_x=crop_x, crop_y=crop_y)
+    warped, _ = landmarks_affine(series, recording_lm, template_lm, output_shape=(25, 30))
 
-    assert warped.shape == (5, crop_y, crop_x)
+    assert warped.shape == (5, 25, 30)
 
 
-def test_landmarks_affine_crop_matches_sequential_reference(landmark_pair):
-    """Cropped parallel output must be identical to the sequential cropped loop."""
-    rng = np.random.default_rng(2)
-    series = rng.random((5, 40, 40), dtype=np.float32)
+def test_landmarks_affine_output_shape_does_not_crop_the_source(landmark_pair):
+    """Sizing the output must not discard source pixels.
+
+    The previous implementation cropped the *input* to (crop_y, crop_x), which threw away exactly
+    the source pixels the warp needs whenever the brain was not in the top-left corner.
+    """
+    series = np.zeros((1, 200, 200), dtype=np.float32)
+    series[0, 120:160, 130:170] = 1.0  # signal well outside the top-left 40x40 corner
     recording_lm, template_lm = landmark_pair
-    crop_x, crop_y = 30, 25
 
-    template  = np.array(list(template_lm.values()),  dtype=np.float32)
-    recording = np.array(list(recording_lm.values()), dtype=np.float32)
-    ref_tform = trf.estimate_transform("affine", template, recording)
-    reference = np.array([trf.warp(series[i, :crop_y, :crop_x], ref_tform, order=3) for i in range(series.shape[0])])
+    # Landmarks describing a pure translation that brings the blob into a 60x60 output window.
+    recording = {"a": (130.0, 120.0), "b": (170.0, 120.0), "c": (130.0, 160.0)}
+    template = {"a": (10.0, 10.0), "b": (50.0, 10.0), "c": (10.0, 50.0)}
 
-    warped, _ = landmarks_affine(series, recording_lm, template_lm, crop_x=crop_x, crop_y=crop_y)
+    warped, _ = landmarks_affine(series, recording, template, output_shape=(60, 60))
 
-    np.testing.assert_array_equal(warped, reference)
+    assert warped.shape == (1, 60, 60)
+    assert warped[0].sum() > 0.9 * series[0].sum()
 
 
 def test_landmarks_affine_identity_landmarks(small_series):
     """When recording and template landmarks are identical the warp is a near-identity."""
     landmarks = {"A": (5.0, 5.0), "B": (35.0, 5.0), "C": (5.0, 35.0)}
-    warped, _ = landmarks_affine(small_series, landmarks, landmarks)
+    warped, _ = landmarks_affine(small_series, landmarks, landmarks, output_shape=(40, 40))
 
     # Interior pixels should be reproduced almost exactly (boundary pixels may differ
     # due to the constant-zero padding convention of skimage warp).
@@ -118,6 +128,40 @@ def test_landmarks_affine_identity_landmarks(small_series):
     np.testing.assert_allclose(warped_interior, interior, atol=1e-5)
 
 
+def test_landmarks_affine_registers_a_known_transform_to_the_atlas():
+    """Ground truth: warp the atlas by a known affine, then register it back.
+
+    Unlike the reference tests above, the expected result here is not derived from the
+    implementation - the recording is synthesised from the atlas itself.
+    """
+    atlas = res.get_atlas()[0].astype(np.float32)
+    template_lm = res.get_default_landmarks()
+
+    # Synthesise a "recording": the atlas seen through a known rotation, scale and shift.
+    true_tform = trf.AffineTransform(scale=(1.8, 1.7), rotation=np.deg2rad(8), translation=(40, 25))
+    recording = trf.warp(atlas, true_tform.inverse, output_shape=(300, 320), order=1)
+    recording_lm = {
+        name: tuple(true_tform(np.array([point]))[0]) for name, point in template_lm.items()
+    }
+
+    warped, tform = landmarks_affine(recording[None], recording_lm, template_lm)
+
+    # Output is in atlas space...
+    assert warped.shape == (1, *atlas.shape)
+
+    # ...the recovered transform is the one we applied...
+    np.testing.assert_allclose(tform.params, true_tform.params, atol=1e-3)
+
+    # ...and every landmark lands within a pixel of its template position.
+    landed = tform.inverse(np.array(list(recording_lm.values())))
+    expected = np.array(list(template_lm.values()))
+    assert np.linalg.norm(landed - expected, axis=1).max() < 1.0
+
+    # The registered image reproduces the atlas it was synthesised from.
+    labelled = atlas > 0
+    assert np.abs(warped[0][labelled] - atlas[labelled]).mean() < 0.02 * atlas.max()
+
+
 def test_update_nwb(nwbfile, preproc_h5):
     tform_mock = np.array([[(1, 2, 3), (1, 2, 4)]])
     nwb = reg.update_nwb(nwbfile, preproc_h5, tform_mock)

From 6635160a73332f2c493eca17bd27ae4c42de19d0 Mon Sep 17 00:00:00 2001
From: Constantinos Eleftheriou <constantinos.eleftheriou@ed.ac.uk>
Date: Thu, 13 Aug 2026 15:33:38 +0100
Subject: [PATCH 04/11] fix landmark identification gui issue where points
 could be silently dropped if deleted. add warnings if landmarks are missing.
 fix ABA scaling issue

---
 src/mesoscopy/register/__init__.py      |  5 +-
 src/mesoscopy/register/landmarks_gui.py | 75 +++++++++++++++++++++++--
 2 files changed, 73 insertions(+), 7 deletions(-)

diff --git a/src/mesoscopy/register/__init__.py b/src/mesoscopy/register/__init__.py
index 2c0365d..2f7b9cd 100644
--- a/src/mesoscopy/register/__init__.py
+++ b/src/mesoscopy/register/__init__.py
@@ -103,11 +103,14 @@ def label_cmd(path, out_dir, template_points, session_id) -> dict:
 
     click.echo("Loading template landmarks...")
     template_landmarks = res.get_default_landmarks()
+    template_shape = res.get_atlas()[0].shape
     if template_points:
         template_landmarks = io.read_points(template_points)
+        # The image a user-supplied template was marked on is unknown, so seed points can't be scaled.
+        template_shape = None
 
     click.echo("Launching landmark identification GUI...")
-    recording_landmarks = reg_gui.mark_landmarks(maxip, isosb_maxip, template_landmarks)
+    recording_landmarks = reg_gui.mark_landmarks(maxip, isosb_maxip, template_landmarks, template_shape=template_shape)
 
     click.echo("Saving recording landmarks...")
     outpath = out_dir + os.sep + session_id + "_landmarks.csv"
diff --git a/src/mesoscopy/register/landmarks_gui.py b/src/mesoscopy/register/landmarks_gui.py
index 10ec9c3..b65a196 100644
--- a/src/mesoscopy/register/landmarks_gui.py
+++ b/src/mesoscopy/register/landmarks_gui.py
@@ -20,6 +20,7 @@
 #  SOFTWARE.
 from collections import OrderedDict
 
+import click
 import magicgui.widgets as mgw
 import napari
 import numpy as np
@@ -40,7 +41,10 @@
 
 
 def mark_landmarks(
-    maxip_image: npt.NDArray | da.Array, alt_image: npt.NDArray | da.Array | None, template_landmarks: dict = {}
+    maxip_image: npt.NDArray | da.Array,
+    alt_image: npt.NDArray | da.Array | None,
+    template_landmarks: dict = {},
+    template_shape: tuple[int, int] | None = None,
 ) -> dict:
     """Launch the napari viewer to identify anatomical landmarks on a maximum intensity projection image.
 
@@ -59,13 +63,19 @@ def mark_landmarks(
     by the template landmark files and by ``skimage.transform``. Napari works in (row, column) order, so
     coordinates are transposed on the way into and out of the viewer.
 
+    Marked points are identified by their ``label`` property rather than by their position in the
+    points layer, since deleting and re-marking a point moves it to the end of the layer.
+
     Args:
         maxip_image (npt.NDArray): Maximum intensity projection image. Could be either channel.
         alt_image (npt.NDArray): Alternative image to be displayed alongside the maximum intensity projection image. Usually a second channel.
         template_landmarks (dict, optional): Dictionary with the landmarks and their x-y coordinates, used to seed the initial point positions. Dictionary keys are landmark names, while x-y coordinates are stored as an (x, y) tuple. Defaults to {}.
+        template_shape (tuple[int, int], optional): Shape of the image the template landmarks were
+            marked on, as (height, width). When given, the seed points are scaled to the size of the
+            recording so they do not bunch up in a corner. Defaults to None (no scaling).
 
     Returns:
-        dict: Dictionary with the landmarks and their x-y coordinates. Dictionary keys are landmark names, while x-y coordinates are stored as an (x, y) tuple.
+        dict: Dictionary with the landmarks and their x-y coordinates, in template order. Dictionary keys are landmark names, while x-y coordinates are stored as an (x, y) tuple.
 
     Raises:
         ValueError: If the maximum intensity projection and alternative image do not have the same dimensions.
@@ -102,6 +112,18 @@ def mark_landmarks(
 
     landmarks = list(default_landmark_locations.keys())
 
+    # Template landmarks are in template pixels, which may be a very different size to the recording.
+    # Scale them so the seed points land roughly on the anatomy instead of bunching up in a corner.
+    if template_landmarks and template_shape is not None:
+        template_height, template_width = template_shape
+        scale_x = maxip_width / template_width
+        scale_y = maxip_height / template_height
+        if not np.isclose(scale_x, 1.0) or not np.isclose(scale_y, 1.0):
+            click.echo(f"Scaling template seed points by {scale_x:.2f} (x) and {scale_y:.2f} (y) to fit the recording.")
+            default_landmark_locations = {
+                landmark: (x * scale_x, y * scale_y) for landmark, (x, y) in default_landmark_locations.items()
+            }
+
     # Landmarks are stored as (x, y), napari points are (row, column) - transpose on the way in.
     seed_points = np.array([(y, x) for x, y in default_landmark_locations.values()], dtype=float)
 
@@ -127,10 +149,51 @@ def mark_landmarks(
 
     napari.run()
 
-    # Transpose back from napari's (row, column) to the stored (x, y) convention.
-    return OrderedDict(
-        (landmark, (float(x), float(y))) for landmark, (y, x) in zip(landmarks, points_layer.data, strict=True)
-    )
+    return _collect_marked_points(points_layer, landmarks)
+
+
+def _collect_marked_points(points_layer: napari.layers.Points, landmarks: list[str]) -> OrderedDict:
+    """Map each landmark name to the position of the point labelled with it.
+
+    Points are matched by their ``label`` property, not by their index in the layer: deleting a point
+    and re-marking it appends it to the end of the layer, so point order does not track landmark
+    order. Results are returned in landmark (template) order, which is the order the registration
+    transform pairs the two point sets in.
+
+    Args:
+        points_layer (napari.layers.Points): The points layer the user annotated.
+        landmarks (list[str]): Landmark names, in template order.
+
+    Returns:
+        OrderedDict: Landmark names mapped to their (x, y) coordinates, in template order.
+
+    Raises:
+        ValueError: If the points layer has no usable label property.
+    """
+    points = np.asarray(points_layer.data, dtype=float)
+    labels = np.asarray(points_layer.properties.get("label", []), dtype=str)
+
+    if len(labels) != len(points):
+        msg = "Points layer labels do not correspond to the marked points, cannot identify landmarks."
+        raise ValueError(msg)
+
+    marked = OrderedDict()
+    for landmark in landmarks:
+        matches = np.flatnonzero(labels == landmark)
+        if len(matches) == 0:
+            click.echo(f"⚠️ No point is labelled '{landmark}' - this landmark will be missing from the output.")
+            continue
+        if len(matches) > 1:
+            click.echo(f"⚠️ {len(matches)} points are labelled '{landmark}' - using the most recently marked one.")
+        # napari points are (row, column), landmarks are stored as (x, y).
+        row, column = points[matches[-1]]
+        marked[landmark] = (float(column), float(row))
+
+    unknown = sorted(set(labels) - set(landmarks))
+    if unknown:
+        click.echo(f"⚠️ Ignoring points with unrecognised labels: {', '.join(unknown)}")
+
+    return marked
 
 
 def _create_label_menu(points_layer, labels):

From e4c929c62387a026ba83e46f88d3e5899b5c6bab Mon Sep 17 00:00:00 2001
From: Constantinos Eleftheriou <constantinos.eleftheriou@ed.ac.uk>
Date: Thu, 13 Aug 2026 15:33:48 +0100
Subject: [PATCH 05/11] add tests for registration gui

---
 tests/test_registration.py | 131 ++++++++++++++++++++++++++++++++-----
 1 file changed, 114 insertions(+), 17 deletions(-)

diff --git a/tests/test_registration.py b/tests/test_registration.py
index acac0f2..6cb499c 100644
--- a/tests/test_registration.py
+++ b/tests/test_registration.py
@@ -1,5 +1,7 @@
 from importlib import resources
+from types import SimpleNamespace
 
+import h5py
 import numpy as np
 import pytest
 from click.testing import CliRunner
@@ -231,38 +233,50 @@ def test_points_roundtrip_preserves_xy(tmp_path):
     assert io.read_points(path) == points
 
 
-def test_mark_landmarks_transposes_between_napari_and_storage(monkeypatch):
-    """Seed points enter napari as (row, col); marked points come back out as (x, y)."""
-    captured = {}
+class _FakePointsLayer:
+    """Stands in for a napari Points layer: (row, col) data plus an aligned 'label' property."""
 
-    class FakePointsLayer:
-        def __init__(self, data):
-            self.data = data
-            self.mode = None
-            self.face_color_mode = None
+    def __init__(self, data, labels):
+        self.data = np.asarray(data, dtype=float)
+        self.properties = {"label": np.asarray(labels, dtype=str)}
+        self.mode = None
+        self.face_color_mode = None
 
-    class FakeWindow:
-        def add_dock_widget(self, widget):
-            pass
+
+def _patch_napari(monkeypatch, final_state=None, captured=None):
+    """Replace the napari viewer with a fake that records seeds and returns a scripted end state.
+
+    `final_state` is the (data, labels) pair the points layer holds when the user closes the viewer,
+    i.e. the result of whatever marking, dragging and deleting they did. Defaults to the seed points
+    untouched.
+    """
 
     class FakeViewer:
         def __init__(self):
-            self.window = FakeWindow()
+            self.window = SimpleNamespace(add_dock_widget=lambda widget: None)
 
         def add_image(self, *args, **kwargs):
             pass
 
         def add_points(self, data, **kwargs):
-            captured["seed"] = np.asarray(data, dtype=float)
-            # Simulate the user dragging the first point by +3 rows / +5 columns.
-            marked = np.asarray(data, dtype=float).copy()
-            marked[0] += (3.0, 5.0)
-            return FakePointsLayer(marked)
+            if captured is not None:
+                captured["seed"] = np.asarray(data, dtype=float)
+                captured["labels"] = list(kwargs["properties"]["label"])
+            if final_state is None:
+                return _FakePointsLayer(data, kwargs["properties"]["label"])
+            return _FakePointsLayer(*final_state)
 
     monkeypatch.setattr(reg_gui.napari, "view_image", lambda *args, **kwargs: FakeViewer())
     monkeypatch.setattr(reg_gui.napari, "run", lambda *args, **kwargs: None)
     monkeypatch.setattr(reg_gui, "_create_label_menu", lambda points_layer, labels: None)
 
+
+def test_mark_landmarks_transposes_between_napari_and_storage(monkeypatch):
+    """Seed points enter napari as (row, col); marked points come back out as (x, y)."""
+    captured = {}
+    # The user drags the first point by +3 rows / +5 columns.
+    _patch_napari(monkeypatch, final_state=([[23.0, 15.0], [40.0, 30.0]], ["a", "b"]), captured=captured)
+
     template = {"a": (10.0, 20.0), "b": (30.0, 40.0)}
     marked = reg_gui.mark_landmarks(np.zeros((100, 120)), None, template)
 
@@ -274,6 +288,89 @@ def add_points(self, data, **kwargs):
     assert marked["b"] == (30.0, 40.0)
 
 
+def test_mark_landmarks_identifies_points_by_label_not_position(monkeypatch):
+    """Deleting and re-marking a point appends it to the layer; labels must still hold.
+
+    Zipping names against point order - as the previous implementation did - silently attaches every
+    name to the wrong coordinates as soon as the user re-marks anything.
+    """
+    template = {"bregma": (71.0, 60.0), "cFP": (71.0, 19.0), "lPB": (30.0, 35.0)}
+    # bregma was deleted and re-marked, so it sits last in the layer and the others shifted up.
+    _patch_napari(
+        monkeypatch,
+        final_state=([[19.0, 71.0], [35.0, 30.0], [63.0, 71.0]], ["cFP", "lPB", "bregma"]),
+    )
+
+    marked = reg_gui.mark_landmarks(np.zeros((140, 142)), None, template)
+
+    assert marked == {"bregma": (71.0, 63.0), "cFP": (71.0, 19.0), "lPB": (30.0, 35.0)}
+    # Returned in template order, which is the order the transform pairs the point sets in.
+    assert list(marked) == list(template)
+
+
+def test_mark_landmarks_reports_missing_duplicate_and_unknown_points(monkeypatch, capsys):
+    """Edited point sets are reported rather than silently mismapped."""
+    template = {"bregma": (71.0, 60.0), "cFP": (71.0, 19.0), "lPB": (30.0, 35.0)}
+    _patch_napari(
+        monkeypatch,
+        final_state=(
+            [[19.0, 71.0], [10.0, 10.0], [40.0, 40.0], [35.0, 30.0]],
+            ["cFP", "bregma", "bregma", "notALandmark"],
+        ),
+    )
+
+    marked = reg_gui.mark_landmarks(np.zeros((140, 142)), None, template)
+
+    assert marked["cFP"] == (71.0, 19.0)
+    assert marked["bregma"] == (40.0, 40.0)  # the most recently marked of the two
+    assert "lPB" not in marked  # never marked
+
+    output = capsys.readouterr().out
+    assert "lPB" in output
+    assert "bregma" in output
+    assert "notALandmark" in output
+
+
+def test_mark_landmarks_scales_seed_points_to_the_recording(monkeypatch):
+    """Template seeds are scaled to the recording so they don't bunch up in a corner."""
+    captured = {}
+    _patch_napari(monkeypatch, captured=captured)
+
+    # Recording is twice the size of the template in both axes.
+    reg_gui.mark_landmarks(np.zeros((280, 284)), None, {"bregma": (71.0, 60.0)}, template_shape=(140, 142))
+
+    np.testing.assert_allclose(captured["seed"], [[120.0, 142.0]])  # (row, col) of (x=142, y=120)
+
+
+def test_mark_landmarks_does_not_scale_seeds_without_a_template_shape(monkeypatch):
+    """Without a known template shape the seeds are used as-is."""
+    captured = {}
+    _patch_napari(monkeypatch, captured=captured)
+
+    reg_gui.mark_landmarks(np.zeros((280, 284)), None, {"bregma": (71.0, 60.0)})
+
+    np.testing.assert_allclose(captured["seed"], [[60.0, 71.0]])
+
+
+def test_register_label_cli_writes_landmarks(monkeypatch, tmp_path):
+    """The label command round-trips GUI points to a landmarks CSV in the (x, y) convention."""
+    preproc = tmp_path / "sub-test_preprocessed.h5"
+    with h5py.File(preproc, "w") as f:
+        f.create_dataset("/qa/gcamp_maxip_projection", data=np.random.default_rng(0).random((140, 142)))
+        f.create_dataset("/qa/isosb_maxip_projection", data=np.random.default_rng(1).random((140, 142)))
+
+    _patch_napari(monkeypatch)  # the user closes the viewer without moving anything
+
+    runner = CliRunner()
+    result = runner.invoke(mesoscopy.cli, args=f"register label {preproc} -o {tmp_path}")
+
+    assert result.exit_code == 0
+
+    # Recording is atlas-sized, so unmoved seeds must round-trip to the template landmarks exactly.
+    written = io.read_points(str(tmp_path / "sub-test_landmarks.csv"))
+    assert written == res.get_default_landmarks()
+
+
 def test_qa_plot_landmarks_uses_xy():
     """The QA scatter must plot x on the x-axis, not the row index."""
     source = np.array([[10.0, 20.0], [30.0, 40.0]])

From c9acc74aa9e0dd868cda34b328292cabd1c0fd65 Mon Sep 17 00:00:00 2001
From: Constantinos Eleftheriou <constantinos.eleftheriou@ed.ac.uk>
Date: Thu, 13 Aug 2026 15:42:20 +0100
Subject: [PATCH 06/11] fix nwb update

---
 src/mesoscopy/register/__init__.py | 84 +++++++++++++++++-------------
 tests/test_registration.py         | 57 +++++++++++++++++---
 2 files changed, 98 insertions(+), 43 deletions(-)

diff --git a/src/mesoscopy/register/__init__.py b/src/mesoscopy/register/__init__.py
index 2f7b9cd..d291341 100644
--- a/src/mesoscopy/register/__init__.py
+++ b/src/mesoscopy/register/__init__.py
@@ -25,6 +25,7 @@
 import click
 import h5py
 import numpy as np
+from pynwb import NWBFile
 from pynwb import TimeSeries
 from pynwb.image import ImageSeries
 from pynwb.ophys import CorrectedImageStack
@@ -237,7 +238,7 @@ def landmarks_cmd(
 
     if nwb:
         click.echo("Updating NWB file...")
-        update_nwb(path, outpath, tform)
+        update_nwb(path, outpath, tform.params)
         click.echo(f"Updated NWB file at {path}")
 
     return outpath
@@ -259,51 +260,60 @@ def load_maxips(path: str) -> tuple[np.ndarray, np.ndarray]:
     return gcamp_maxip_projection, isosb_maxip_projection
 
 
-def update_nwb(nwb_path: str, h5_path: str, tform_params: np.ndarray) -> None:
+def update_nwb(nwb_path: str, h5_path: str, tform_params: np.ndarray) -> NWBFile:
     """Update an NWB file with registered imaging data stored in an HDF5 file.
 
     Creates a link between the NWB file and the HDF5 file. See https://pynwb.readthedocs.io/en/stable/tutorials/advanced_io/linking_data.html.
 
+    The registration is a single global affine, so the same 3x3 matrix is stored for every frame of
+    the xy_translation series, giving it a shape of (n_timestamps, 3, 3).
+
     Args:
         nwb_path (str): Path to the NWB file.
         h5_path (str): Path to the HDF5 file containing the registered images.
-        tform_params (np.ndarray): Affine transformation parameters.
+        tform_params (np.ndarray): Affine transformation matrix, as a 3x3 array.
+
+    Returns:
+        NWBFile: The updated NWB file object. Note that its link to the HDF5 file is closed on
+            return, so the registered image data is only readable by re-opening nwb_path.
     """
     nwbfile, nwbio = io.read_nwb(nwb_path, return_io=True)
-    f = h5py.File(h5_path, "r")
-
-    try:
-        ophys_module = nwbfile.create_processing_module(name="ophys", description="optical physiology processed data")
-    except ValueError:
-        click.echo("Processing module already exists...")
-        ophys_module = nwbfile.processing["ophys"]
-
-    registered_series = ImageSeries(
-        name="corrected",
-        data=f["/F"],
-        timestamps=f["/timestamps"],
-        unit="df/f",
-        description="dF/F widefield cortical imaging series.",
-        comments="This is the haemodynamic corrected series registered to the Allen Brain Atlas CCFv3.",
-    )
-
-    xy_translation = TimeSeries(
-        name="xy_translation",
-        data=np.repeat(tform_params, len(f["/timestamps"]), axis=0),
-        unit="pixels",
-        timestamps=f["/timestamps"],
-        description="Affine transformation parameters for image registration to the ABA CCFv3.",
-    )
-
-    corrected_image_stack = CorrectedImageStack(
-        name="CCFRegisteredSeries",
-        corrected=registered_series,
-        original=nwbfile.acquisition["DualChannelImagingSeries"],
-        xy_translation=xy_translation,
-    )
-
-    ophys_module.add(corrected_image_stack)
 
-    io.write_nwb(nwb_path, nwbfile, io=nwbio)
+    with h5py.File(h5_path, "r") as f:
+        try:
+            ophys_module = nwbfile.create_processing_module(
+                name="ophys", description="optical physiology processed data"
+            )
+        except ValueError:
+            click.echo("Processing module already exists...")
+            ophys_module = nwbfile.processing["ophys"]
+
+        registered_series = ImageSeries(
+            name="corrected",
+            data=f["/F"],
+            timestamps=f["/timestamps"],
+            unit="df/f",
+            description="dF/F widefield cortical imaging series.",
+            comments="This is the haemodynamic corrected series registered to the Allen Brain Atlas CCFv3.",
+        )
+
+        xy_translation = TimeSeries(
+            name="xy_translation",
+            data=np.tile(tform_params, (len(f["/timestamps"]), 1, 1)),
+            unit="pixels",
+            timestamps=f["/timestamps"],
+            description="Affine transformation parameters for image registration to the ABA CCFv3.",
+        )
+
+        corrected_image_stack = CorrectedImageStack(
+            name="CCFRegisteredSeries",
+            corrected=registered_series,
+            original=nwbfile.acquisition["DualChannelImagingSeries"],
+            xy_translation=xy_translation,
+        )
+
+        ophys_module.add(corrected_image_stack)
+
+        io.write_nwb(nwb_path, nwbfile, io=nwbio)
 
     return nwbfile
diff --git a/tests/test_registration.py b/tests/test_registration.py
index 6cb499c..363c04a 100644
--- a/tests/test_registration.py
+++ b/tests/test_registration.py
@@ -142,9 +142,7 @@ def test_landmarks_affine_registers_a_known_transform_to_the_atlas():
     # Synthesise a "recording": the atlas seen through a known rotation, scale and shift.
     true_tform = trf.AffineTransform(scale=(1.8, 1.7), rotation=np.deg2rad(8), translation=(40, 25))
     recording = trf.warp(atlas, true_tform.inverse, output_shape=(300, 320), order=1)
-    recording_lm = {
-        name: tuple(true_tform(np.array([point]))[0]) for name, point in template_lm.items()
-    }
+    recording_lm = {name: tuple(true_tform(np.array([point]))[0]) for name, point in template_lm.items()}
 
     warped, tform = landmarks_affine(recording[None], recording_lm, template_lm)
 
@@ -165,13 +163,60 @@ def test_landmarks_affine_registers_a_known_transform_to_the_atlas():
 
 
 def test_update_nwb(nwbfile, preproc_h5):
-    tform_mock = np.array([[(1, 2, 3), (1, 2, 4)]])
-    nwb = reg.update_nwb(nwbfile, preproc_h5, tform_mock)
-    assert nwb.processing["ophys"]["CCFRegisteredSeries"].corrected.data
+    tform_params = trf.AffineTransform(scale=(1.2, 1.1), rotation=0.05, translation=(3, -4)).params
+
+    nwb = reg.update_nwb(nwbfile, preproc_h5, tform_params)
+
     assert nwb.processing["ophys"]["CCFRegisteredSeries"].original.data
     assert nwb.processing["ophys"]["CCFRegisteredSeries"].xy_translation.data.any()
 
 
+def test_update_nwb_stores_one_transform_matrix_per_timestamp(nwbfile, preproc_h5):
+    """xy_translation must be (n_timestamps, 3, 3), not the flattened (3 * n_timestamps, 3).
+
+    pynwb only warns when data and timestamps disagree in length, so a mis-shaped array is written
+    without complaint.
+    """
+    tform_params = trf.AffineTransform(scale=(1.2, 1.1), rotation=0.05, translation=(3, -4)).params
+
+    reg.update_nwb(nwbfile, preproc_h5, tform_params)
+
+    # Read the file back rather than trusting the in-memory object.
+    with h5py.File(preproc_h5, "r") as f:
+        n_timestamps = len(f["/timestamps"])
+
+    written = io.read_nwb(nwbfile, mode="r")
+    stack = written.processing["ophys"]["CCFRegisteredSeries"]
+    xy_translation = np.asarray(stack.xy_translation.data)
+
+    assert xy_translation.shape == (n_timestamps, 3, 3)
+    assert len(xy_translation) == len(stack.xy_translation.timestamps)
+    # Every frame carries the same global transform.
+    for frame_params in xy_translation:
+        np.testing.assert_allclose(frame_params, tform_params)
+
+
+def test_register_landmarks_cli_nwb_stores_the_transform_matrix(preproc_nwb, output_dir):
+    """End-to-end NWB path: the stored transform must be the matrix, correctly shaped per frame."""
+    points = str(resources.files(res).joinpath("ccf_template_landmarks_140x142.csv"))
+    runner = CliRunner()
+
+    result = runner.invoke(
+        mesoscopy.cli,
+        args=f"register landmarks {preproc_nwb} -r {points} -o {output_dir}",
+    )
+
+    assert result.exit_code == 0, result.output
+
+    written = io.read_nwb(preproc_nwb, mode="r")
+    stack = written.processing["ophys"]["CCFRegisteredSeries"]
+    xy_translation = np.asarray(stack.xy_translation.data)
+
+    assert xy_translation.shape == (len(stack.xy_translation.timestamps), 3, 3)
+    # Recording and template landmarks are the same file here, so the transform is the identity.
+    np.testing.assert_allclose(xy_translation[0], np.eye(3), atol=1e-6)
+
+
 def test_register_landmarks_cli_default_points_h5(preproc_h5, output_dir):
     mock_recording_landmarks = str(resources.files(res).joinpath("ccf_template_landmarks_140x142.csv"))
     runner = CliRunner()

From f6564e4a9af97fa5fb2294263bf5e92bdaf8352a Mon Sep 17 00:00:00 2001
From: Constantinos Eleftheriou <constantinos.eleftheriou@ed.ac.uk>
Date: Thu, 13 Aug 2026 15:51:33 +0100
Subject: [PATCH 07/11] fix automagic landmark file discovery, fix maxip
 discovery for nwb files

---
 src/mesoscopy/register/__init__.py | 118 +++++++++++++++++++++++------
 tests/test_registration.py         |  95 ++++++++++++++++++++++-
 2 files changed, 189 insertions(+), 24 deletions(-)

diff --git a/src/mesoscopy/register/__init__.py b/src/mesoscopy/register/__init__.py
index d291341..d807bc0 100644
--- a/src/mesoscopy/register/__init__.py
+++ b/src/mesoscopy/register/__init__.py
@@ -25,12 +25,12 @@
 import click
 import h5py
 import numpy as np
+from dask import array as da
 from pynwb import NWBFile
 from pynwb import TimeSeries
 from pynwb.image import ImageSeries
 from pynwb.ophys import CorrectedImageStack
 
-import mesoscopy.preprocess as preproc
 import mesoscopy.preprocess.compute as preproc_compute
 import mesoscopy.register.landmarks_gui as reg_gui
 import mesoscopy.register.transform as trf
@@ -85,22 +85,27 @@ def label_cmd(path, out_dir, template_points, session_id) -> dict:
     nwb = bool(path.endswith(".nwb"))
 
     if not session_id:
-        session_id = path.split("/")[-1].replace(".nwb", "") if nwb else path.split("/")[-1].replace(".h5", "")
-        session_id = session_id.replace("_preprocessed", "")
+        session_id = session_id_from_path(path)
 
     maxip = None
     isosb_maxip = None
 
-    # Load maxip from preprocessed file.
-    # if it does not exist, generate maxip from raw data.
-    if path.endswith("_preprocessed.h5"):
-        maxip, isosb_maxip = load_maxips(path)
-    else:
-        # generate maxip from raw data
-        click.echo("⚠️ No maximum intensity projection found. Generating from raw data, this might take some time...")
+    # Prefer the projections written by preprocessing: they are in the same pixel space as the dF/F
+    # series, which is the space the landmarks have to be marked in. An NWB file links its dF/F
+    # series to the preprocessed HDF5 file, so the projections can be read from there.
+    source = linked_preprocessed_path(path) if nwb else path
+    if source and pathlib.Path(source).exists():
+        maxip, isosb_maxip = load_maxips(source)
+
+    if maxip is None:
+        # Fall back to projecting the dF/F series itself. It is a poorer anatomical image than the
+        # gcamp projection, but it is guaranteed to be in the same pixel space as the data being
+        # registered - projecting the raw frames instead would be off by the preprocessing crop and
+        # binning factor, silently scaling the transform.
+        click.echo("⚠️ No preprocessing maximum intensity projection found, projecting the ∆F/F series instead.")
         with timer.Timer("Generating maximum intensity projection"):
-            _, raw_data, _ = preproc.load_raw(path, nwb=nwb)
-            maxip = preproc_compute.projections(raw_data)["maxip"]
+            _, deltaf_series, _ = io.load_deltaf(path, nwb=nwb)
+            maxip = preproc_compute.projections(da.from_array(deltaf_series))["maxip"]
 
     click.echo("Loading template landmarks...")
     template_landmarks = res.get_default_landmarks()
@@ -199,13 +204,17 @@ def landmarks_cmd(
         template_landmarks = io.read_points(template_points)
 
     if not recording_points:
-        if nwb and pathlib.Path(path.replace(".nwb", "_landmarks.csv")).exists():
-            recording_points = path.replace(".nwb", "_landmarks.csv")
-        elif pathlib.Path(path.replace(".h5", "_landmarks.csv")).exists():
-            recording_points = path.replace(".h5", "_landmarks.csv")
-        else:
-            msg = "Path to recording landmarks could not be inferred. Please supply a recording landmarks file."
+        candidates = landmarks_path_candidates(path, out_dir)
+        found = next((candidate for candidate in candidates if candidate.exists()), None)
+        if found is None:
+            searched = "\n  ".join(str(candidate) for candidate in candidates)
+            msg = (
+                "Path to recording landmarks could not be inferred. Searched:\n  "
+                f"{searched}\nPlease supply a recording landmarks file with -r/--recording-points."
+            )
             raise ValueError(msg)
+        recording_points = str(found)
+        click.echo(f"Using recording landmarks at {recording_points}")
     recording_landmarks = io.read_points(recording_points)
 
     # Registered frames land in template space, so default their shape to that of the CCF atlas.
@@ -244,22 +253,85 @@ def landmarks_cmd(
     return outpath
 
 
-def load_maxips(path: str) -> tuple[np.ndarray, np.ndarray]:
+def session_id_from_path(path: str) -> str:
+    """Derive a session identifier from a recording path.
+
+    Args:
+        path (str): Path to a recording file.
+
+    Returns:
+        str: The file name without its extension or the "_preprocessed" suffix.
+    """
+    return pathlib.Path(path).stem.replace("_preprocessed", "")
+
+
+def landmarks_path_candidates(path: str, out_dir: str | None = None) -> list[pathlib.Path]:
+    """List the paths a recording's landmarks file may have been written to.
+
+    ``register label`` names its output after the session ID, which drops the "_preprocessed"
+    suffix, so the landmarks file rarely sits at ``<recording>_landmarks.csv``.
+
+    Args:
+        path (str): Path to the recording file.
+        out_dir (str, optional): Output directory the landmarks may have been written to.
+
+    Returns:
+        list[pathlib.Path]: Candidate landmark file paths, in search order, without duplicates.
+    """
+    recording = pathlib.Path(path)
+    names = [f"{session_id_from_path(path)}_landmarks.csv", f"{recording.stem}_landmarks.csv"]
+    directories = [recording.parent]
+    if out_dir:
+        directories.append(pathlib.Path(out_dir))
+
+    candidates = [directory / name for directory in directories for name in names]
+    return list(dict.fromkeys(candidates))
+
+
+def load_maxips(path: str) -> tuple[np.ndarray | None, np.ndarray | None]:
     """Load maximum intensity projections from a preprocessed HDF5 file.
 
     Args:
         path (str): Path to the preprocessed HDF5 file.
 
     Returns:
-        tuple[np.ndarray, np.ndarray]: Maximum intensity projection for gcamp and isosb channels.
+        tuple[np.ndarray | None, np.ndarray | None]: Maximum intensity projection for the gcamp and
+            isosb channels, or (None, None) if the file holds no gcamp projection.
     """
-    f_preproc = h5py.File(path, "r")
-    gcamp_maxip_projection = np.array(f_preproc["/qa/gcamp_maxip_projection"])
-    isosb_maxip_projection = np.array(f_preproc["/qa/isosb_maxip_projection"])
+    if not h5py.is_hdf5(path):
+        return None, None
+
+    with h5py.File(path, "r") as f_preproc:
+        if "/qa/gcamp_maxip_projection" not in f_preproc:
+            return None, None
+
+        gcamp_maxip_projection = np.array(f_preproc["/qa/gcamp_maxip_projection"])
+        isosb_maxip_projection = None
+        if "/qa/isosb_maxip_projection" in f_preproc:
+            isosb_maxip_projection = np.array(f_preproc["/qa/isosb_maxip_projection"])
 
     return gcamp_maxip_projection, isosb_maxip_projection
 
 
+def linked_preprocessed_path(nwb_path: str) -> str | None:
+    """Resolve the preprocessed HDF5 file that an NWB file links its dF/F series to.
+
+    Args:
+        nwb_path (str): Path to the NWB file.
+
+    Returns:
+        str | None: Path to the linked HDF5 file, or None if the dF/F series is not an external link.
+    """
+    with h5py.File(nwb_path, "r") as f:
+        link = f.get("processing/ophys/DeltaFSeries/data", getlink=True)
+
+    if not isinstance(link, h5py.ExternalLink):
+        return None
+
+    # External link targets are stored relative to the NWB file (absolute paths pass through).
+    return str(pathlib.Path(nwb_path).parent / link.filename)
+
+
 def update_nwb(nwb_path: str, h5_path: str, tform_params: np.ndarray) -> NWBFile:
     """Update an NWB file with registered imaging data stored in an HDF5 file.
 
diff --git a/tests/test_registration.py b/tests/test_registration.py
index 363c04a..c7890cf 100644
--- a/tests/test_registration.py
+++ b/tests/test_registration.py
@@ -1,3 +1,5 @@
+import pathlib
+import shutil
 from importlib import resources
 from types import SimpleNamespace
 
@@ -8,6 +10,7 @@
 from skimage import transform as trf
 
 import mesoscopy
+import mesoscopy.preprocess as preproc
 import mesoscopy.register as reg
 import mesoscopy.register.landmarks_gui as reg_gui
 import mesoscopy.register.qa as reg_qa
@@ -217,6 +220,91 @@ def test_register_landmarks_cli_nwb_stores_the_transform_matrix(preproc_nwb, out
     np.testing.assert_allclose(xy_translation[0], np.eye(3), atol=1e-6)
 
 
+# ---------------------------------------------------------------------------
+# Landmark file discovery and maxip sourcing
+# ---------------------------------------------------------------------------
+
+
+def test_session_id_drops_extension_and_preprocessed_suffix():
+    assert reg.session_id_from_path("/data/sub-01_preprocessed.h5") == "sub-01"
+    assert reg.session_id_from_path("/data/sub-01.nwb") == "sub-01"
+    # A directory containing ".h5" must not be mangled - the old str.replace did exactly that.
+    assert reg.session_id_from_path("/data.h5archive/sub-01_preprocessed.h5") == "sub-01"
+
+
+def test_landmarks_inference_finds_the_file_label_actually_writes(preproc_h5, output_dir, tmp_path):
+    """`label` names its output after the session ID, so it drops the "_preprocessed" suffix.
+
+    The previous inference looked for "<recording>_landmarks.csv", which for a preprocessed HDF5
+    could never match what `label` had written.
+    """
+    recording = tmp_path / "sub-01_preprocessed.h5"
+    shutil.copy(preproc_h5, recording)
+    # Exactly what `register label` would have written for this recording.
+    io.write_points(str(tmp_path / "sub-01_landmarks.csv"), res.get_default_landmarks())
+
+    runner = CliRunner()
+    result = runner.invoke(mesoscopy.cli, args=f"register landmarks {recording} -o {output_dir}")
+
+    assert result.exit_code == 0, result.output
+    assert "sub-01_landmarks.csv" in result.output
+
+
+def test_landmarks_inference_error_lists_where_it_looked(preproc_h5, output_dir):
+    runner = CliRunner()
+    result = runner.invoke(mesoscopy.cli, args=f"register landmarks {preproc_h5} -o {output_dir}")
+
+    assert result.exit_code != 0
+    assert isinstance(result.exception, ValueError)
+    assert "Searched" in str(result.exception)
+    assert "_landmarks.csv" in str(result.exception)
+
+
+def test_load_maxips_returns_none_when_absent(preproc_h5):
+    """preproc_h5 holds only /F and /timestamps, so there is nothing to load."""
+    assert reg.load_maxips(preproc_h5) == (None, None)
+
+
+def test_linked_preprocessed_path_resolves_the_external_link(nwbfile, preproc_h5):
+    """Preprocessing links the NWB dF/F series to the HDF5 file holding the QA projections."""
+    preproc.update_nwb(nwbfile, preproc_h5)
+
+    resolved = reg.linked_preprocessed_path(nwbfile)
+
+    assert resolved is not None
+    assert pathlib.Path(resolved).resolve() == pathlib.Path(preproc_h5).resolve()
+
+
+def test_linked_preprocessed_path_is_none_without_an_external_link(preproc_nwb):
+    """A dF/F series stored inline rather than linked resolves to nothing."""
+    assert reg.linked_preprocessed_path(preproc_nwb) is None
+
+
+def test_label_falls_back_to_a_deltaf_projection(monkeypatch, tmp_path):
+    """Without QA projections the maxip must come from the dF/F series, not the raw frames.
+
+    The raw frames are neither cropped nor binned, so projecting them yields landmarks in a
+    different pixel space to the data being registered.
+    """
+    recording = tmp_path / "sub-01_preprocessed.h5"
+    rng = np.random.default_rng(0)
+    deltaf = rng.random((10, 40, 40), dtype=np.float32)
+    with h5py.File(recording, "w") as f:
+        f.create_dataset("/F", data=deltaf)
+        f.create_dataset("/timestamps", data=np.arange(10.0))
+
+    captured = {}
+    _patch_napari(monkeypatch, captured=captured)
+
+    runner = CliRunner()
+    result = runner.invoke(mesoscopy.cli, args=f"register label {recording} -o {tmp_path}")
+
+    assert result.exit_code == 0, result.output
+    assert "∆F/F" in result.output
+    # The projection is of the dF/F series, in the dF/F pixel space.
+    np.testing.assert_allclose(captured["image"], deltaf.max(axis=0))
+
+
 def test_register_landmarks_cli_default_points_h5(preproc_h5, output_dir):
     mock_recording_landmarks = str(resources.files(res).joinpath("ccf_template_landmarks_140x142.csv"))
     runner = CliRunner()
@@ -311,7 +399,12 @@ def add_points(self, data, **kwargs):
                 return _FakePointsLayer(data, kwargs["properties"]["label"])
             return _FakePointsLayer(*final_state)
 
-    monkeypatch.setattr(reg_gui.napari, "view_image", lambda *args, **kwargs: FakeViewer())
+    def view_image(image, *args, **kwargs):
+        if captured is not None:
+            captured["image"] = np.asarray(image)
+        return FakeViewer()
+
+    monkeypatch.setattr(reg_gui.napari, "view_image", view_image)
     monkeypatch.setattr(reg_gui.napari, "run", lambda *args, **kwargs: None)
     monkeypatch.setattr(reg_gui, "_create_label_menu", lambda points_layer, labels: None)
 

From e02f37c2aecf0fab40ff942fd9f5914e7e30692e Mon Sep 17 00:00:00 2001
From: Constantinos Eleftheriou <constantinos.eleftheriou@ed.ac.uk>
Date: Thu, 13 Aug 2026 16:02:45 +0100
Subject: [PATCH 08/11] fix pair name matching for points, add qa check for fit
 quality

---
 src/mesoscopy/register/__init__.py  |  11 ++-
 src/mesoscopy/register/transform.py | 125 +++++++++++++++++++++++-
 tests/test_registration.py          | 141 ++++++++++++++++++++++++++++
 3 files changed, 271 insertions(+), 6 deletions(-)

diff --git a/src/mesoscopy/register/__init__.py b/src/mesoscopy/register/__init__.py
index d807bc0..ae4824e 100644
--- a/src/mesoscopy/register/__init__.py
+++ b/src/mesoscopy/register/__init__.py
@@ -230,6 +230,9 @@ def landmarks_cmd(
         output_shape=output_shape,
     )
 
+    # Store the name-matched point pairs, so every QA dataset corresponds row for row.
+    landmark_names, aligned_template, aligned_recording = trf.align_landmarks(recording_landmarks, template_landmarks)
+
     # Save warped frames and timestamps
     outpath = out_dir + os.sep + session_id + "_registered.h5"
     outpath = io.write_h5(
@@ -238,9 +241,11 @@ def landmarks_cmd(
             "/F": warped,
             "/timestamps": timestamps,
             "/tform": tform.params,
-            "/qa/recording_landmarks": np.array(list(recording_landmarks.values())),
-            "/qa/template_landmarks": np.array(list(template_landmarks.values())),
-            "/qa/registered_landmarks": tform.inverse(np.array(list(recording_landmarks.values()))),
+            "/qa/landmark_names": np.array(landmark_names, dtype="S"),
+            "/qa/recording_landmarks": aligned_recording,
+            "/qa/template_landmarks": aligned_template,
+            "/qa/registered_landmarks": tform.inverse(aligned_recording),
+            "/qa/landmark_residuals": trf.landmark_residuals(tform, aligned_template, aligned_recording),
         },
     )
     click.echo(f"Saved registered frames at {outpath}")
diff --git a/src/mesoscopy/register/transform.py b/src/mesoscopy/register/transform.py
index 66318e2..0fd539a 100644
--- a/src/mesoscopy/register/transform.py
+++ b/src/mesoscopy/register/transform.py
@@ -29,6 +29,93 @@
 
 import mesoscopy.resources as res
 
+# An affine fit needs at least three non-collinear point pairs.
+MIN_LANDMARKS = 3
+
+# Warn when the fit leaves landmarks further than this fraction of the template's extent from their
+# target positions. 0.05 is roughly 6 px on the Allen CCF template.
+RESIDUAL_WARN_FRACTION = 0.05
+
+
+def align_landmarks(recording_landmarks: dict, template_landmarks: dict) -> tuple[list[str], np.ndarray, np.ndarray]:
+    """Pair two landmark sets by name.
+
+    The transform is fitted from two arrays of points, so the pairing between them is positional.
+    Building those arrays from ``dict.values()`` silently mispairs the points whenever the two sets
+    are ordered differently or one of them is missing a landmark, so they are matched by name here.
+
+    Args:
+        recording_landmarks (dict): Recording landmarks, as {name: (x, y)}.
+        template_landmarks (dict): Template landmarks, as {name: (x, y)}.
+
+    Returns:
+        tuple[list[str], np.ndarray, np.ndarray]: The shared landmark names in template order, and
+            the matching template and recording points as (n, 2) arrays of (x, y) coordinates.
+
+    Raises:
+        ValueError: If the two sets share fewer than MIN_LANDMARKS landmarks.
+    """
+    shared = [name for name in template_landmarks if name in recording_landmarks]
+
+    unmarked = [name for name in template_landmarks if name not in recording_landmarks]
+    if unmarked:
+        click.echo(f"⚠️ Template landmarks with no matching recording landmark: {', '.join(unmarked)}")
+
+    unknown = [name for name in recording_landmarks if name not in template_landmarks]
+    if unknown:
+        click.echo(f"⚠️ Recording landmarks with no matching template landmark: {', '.join(unknown)}")
+
+    if len(shared) < MIN_LANDMARKS:
+        msg = (
+            f"Only {len(shared)} landmark(s) are common to the recording and the template, "
+            f"at least {MIN_LANDMARKS} are needed to fit an affine transform."
+        )
+        raise ValueError(msg)
+
+    template = np.array([template_landmarks[name] for name in shared], dtype=np.float64)
+    recording = np.array([recording_landmarks[name] for name in shared], dtype=np.float64)
+
+    return shared, template, recording
+
+
+def _is_collinear(points: np.ndarray) -> bool:
+    """Check whether a set of points lies on a single line (or is a single repeated point).
+
+    Args:
+        points (np.ndarray): Points as an (n, 2) array.
+
+    Returns:
+        bool: True if the points span fewer than two dimensions.
+    """
+    centred = points - points.mean(axis=0)
+    singular_values = np.linalg.svd(centred, compute_uv=False)
+    return bool(singular_values[1] <= 1e-8 * singular_values[0])
+
+
+def landmark_residuals(
+    tform: trf.ProjectiveTransform,
+    template_points: np.ndarray,
+    recording_points: np.ndarray,
+) -> np.ndarray:
+    """Measure how far each marked landmark lands from its template position once registered.
+
+    Residuals are reported in template pixels, i.e. in the space of the registered frames, so they
+    are comparable across recordings of different sizes.
+
+    Note that this does not detect a transposed coordinate convention: an affine least-squares fit
+    solves each output axis independently, so swapping the axes of one point set permutes the rows
+    of the fitted matrix and leaves the residuals unchanged.
+
+    Args:
+        tform (trf.ProjectiveTransform): The fitted transform, mapping template to recording space.
+        template_points (np.ndarray): Template points as an (n, 2) array of (x, y) coordinates.
+        recording_points (np.ndarray): Matching recording points, in the same order.
+
+    Returns:
+        np.ndarray: Residual distance per landmark, in template pixels.
+    """
+    return np.linalg.norm(tform.inverse(recording_points) - template_points, axis=1)
+
 
 def landmarks_affine(
     deltaf_series: np.ndarray,
@@ -54,6 +141,9 @@ def landmarks_affine(
     Returns:
         tuple[np.ndarray, trf.ProjectiveTransform]: Registered DeltaF/F series and affine
             transformation matrix.
+
+    Raises:
+        ValueError: If the landmarks do not define a usable affine transform.
     """
     if output_shape is None:
         output_shape = res.get_atlas()[0].shape
@@ -61,13 +151,42 @@ def landmarks_affine(
         click.echo("Loading imaging data into memory...")
         deltaf_series = np.asarray(deltaf_series)
 
-    template = np.array(list(template_landmarks.values()), dtype=np.float32)
-    recording = np.array(list(recording_landmarks.values()), dtype=np.float32)
+    names, template, recording = align_landmarks(recording_landmarks, template_landmarks)
 
-    click.echo("Estimating transform...")
+    click.echo(f"Estimating transform from {len(names)} landmarks...")
     start = time.time()
     tform = trf.estimate_transform("affine", template, recording)
     end = time.time()
+
+    # skimage returns a least-squares solution without complaint for degenerate point sets, so the
+    # points have to be checked rather than the fit: collinear landmarks leave the fit
+    # underdetermined but still produce a plausible-looking, non-singular matrix.
+    for name, points in (("template", template), ("recording", recording)):
+        if _is_collinear(points):
+            msg = (
+                f"The {name} landmarks are collinear or coincident, so they do not define an affine "
+                "transform. Check the marked points."
+            )
+            raise ValueError(msg)
+
+    if not np.isfinite(tform.params).all():
+        msg = "Could not estimate a transform from these landmarks, the fit did not converge."
+        raise ValueError(msg)
+
+    residuals = landmark_residuals(tform, template, recording)
+    worst = int(np.argmax(residuals))
+    rmse = float(np.sqrt((residuals**2).mean()))
+    click.echo(
+        f"Landmark fit: RMSE {rmse:.2f} px, worst is '{names[worst]}' at {residuals[worst]:.2f} px "
+        "(in template pixels)."
+    )
+
+    template_extent = float(np.linalg.norm(template.max(axis=0) - template.min(axis=0)))
+    if rmse > RESIDUAL_WARN_FRACTION * template_extent:
+        click.echo(
+            f"⚠️ Landmark fit RMSE is over {RESIDUAL_WARN_FRACTION:.0%} of the template extent - "
+            "the registration may be poor. Check the landmarks and the registration QA report."
+        )
     click.echo(f"Transform estimated in {end - start} s")
 
     n_frames = deltaf_series.shape[0]
diff --git a/tests/test_registration.py b/tests/test_registration.py
index c7890cf..35be3c3 100644
--- a/tests/test_registration.py
+++ b/tests/test_registration.py
@@ -16,6 +16,7 @@
 import mesoscopy.register.qa as reg_qa
 import mesoscopy.resources as res
 from mesoscopy import io
+import mesoscopy.register.transform as trf_mod
 from mesoscopy.register.transform import landmarks_affine
 
 
@@ -165,6 +166,146 @@ def test_landmarks_affine_registers_a_known_transform_to_the_atlas():
     assert np.abs(warped[0][labelled] - atlas[labelled]).mean() < 0.02 * atlas.max()
 
 
+# ---------------------------------------------------------------------------
+# Landmark pairing and fit quality
+# ---------------------------------------------------------------------------
+
+
+def test_align_landmarks_pairs_by_name_not_order():
+    """Point pairing is positional once the arrays are built, so names must drive the ordering."""
+    template = {"bregma": (71.0, 60.0), "cFP": (71.0, 19.0), "lPB": (30.0, 35.0)}
+    recording = {"lPB": (32.0, 36.0), "bregma": (70.0, 61.0), "cFP": (72.0, 20.0)}
+
+    names, template_points, recording_points = trf_mod.align_landmarks(recording, template)
+
+    assert names == ["bregma", "cFP", "lPB"]
+    np.testing.assert_allclose(template_points, [[71, 60], [71, 19], [30, 35]])
+    np.testing.assert_allclose(recording_points, [[70, 61], [72, 20], [32, 36]])
+
+
+def test_landmarks_affine_is_unaffected_by_landmark_ordering(small_series):
+    """Reordering a landmark file must not change the fitted transform."""
+    template = {"a": (10.0, 10.0), "b": (30.0, 10.0), "c": (10.0, 25.0)}
+    recording = {"a": (12.0, 13.0), "b": (32.0, 13.0), "c": (12.0, 28.0)}
+    shuffled = {"c": recording["c"], "a": recording["a"], "b": recording["b"]}
+
+    _, tform = landmarks_affine(small_series, recording, template)
+    _, shuffled_tform = landmarks_affine(small_series, shuffled, template)
+
+    np.testing.assert_allclose(tform.params, shuffled_tform.params, atol=1e-9)
+
+
+def test_align_landmarks_reports_unmatched_names(capsys):
+    template = {"a": (10.0, 10.0), "b": (30.0, 10.0), "c": (10.0, 25.0), "onlyInTemplate": (1.0, 1.0)}
+    recording = {"a": (12.0, 13.0), "b": (32.0, 13.0), "c": (12.0, 28.0), "onlyInRecording": (2.0, 2.0)}
+
+    names, _, _ = trf_mod.align_landmarks(recording, template)
+
+    assert names == ["a", "b", "c"]
+    output = capsys.readouterr().out
+    assert "onlyInTemplate" in output
+    assert "onlyInRecording" in output
+
+
+def test_align_landmarks_rejects_too_few_shared_landmarks():
+    template = {"a": (10.0, 10.0), "b": (30.0, 10.0), "c": (10.0, 25.0)}
+    recording = {"a": (12.0, 13.0), "b": (32.0, 13.0)}
+
+    with pytest.raises(ValueError, match="at least 3"):
+        trf_mod.align_landmarks(recording, template)
+
+
+def test_landmarks_affine_rejects_collinear_landmarks(small_series):
+    """Collinear points leave the fit underdetermined.
+
+    skimage reports no error for them - it returns a wild but non-singular matrix - so the check has
+    to be on the geometry of the points, not on the determinant of the fitted transform.
+    """
+    template = {"a": (10.0, 10.0), "b": (20.0, 20.0), "c": (30.0, 30.0), "d": (40.0, 40.0)}
+    recording = {name: (x + 1.0, y + 1.0) for name, (x, y) in template.items()}
+
+    with pytest.raises(ValueError, match="collinear"):
+        landmarks_affine(small_series, recording, template)
+
+
+def test_landmarks_affine_rejects_coincident_landmarks(small_series):
+    template = {"a": (10.0, 10.0), "b": (30.0, 10.0), "c": (10.0, 25.0)}
+    recording = dict.fromkeys(template, (5.0, 5.0))
+
+    with pytest.raises(ValueError, match="collinear or coincident"):
+        landmarks_affine(small_series, recording, template)
+
+
+def test_landmarks_affine_reports_fit_residuals(small_series, capsys):
+    """A perfect fit reports ~0 px; a displaced landmark raises the RMSE."""
+    template = {"a": (10.0, 10.0), "b": (30.0, 10.0), "c": (10.0, 25.0), "d": (30.0, 25.0)}
+    exact = {name: (x + 2.0, y + 3.0) for name, (x, y) in template.items()}
+
+    landmarks_affine(small_series, exact, template)
+    assert "RMSE 0.00 px" in capsys.readouterr().out
+
+    displaced = dict(exact)
+    displaced["d"] = (exact["d"][0] + 8.0, exact["d"][1])
+    landmarks_affine(small_series, displaced, template)
+
+    output = capsys.readouterr().out
+    assert "RMSE 0.00 px" not in output
+    assert "worst is" in output
+
+
+def test_landmarks_affine_warns_about_a_poor_fit(small_series, capsys):
+    """A fit that leaves landmarks far from their targets is called out."""
+    template = {"a": (10.0, 10.0), "b": (30.0, 10.0), "c": (10.0, 25.0), "d": (30.0, 25.0)}
+    good = {name: (x + 2.0, y + 3.0) for name, (x, y) in template.items()}
+
+    landmarks_affine(small_series, good, template)
+    assert "⚠️" not in capsys.readouterr().out
+
+    # Scramble one point far enough that no affine can fit the set.
+    poor = dict(good)
+    poor["d"] = (good["d"][0] - 40.0, good["d"][1] + 30.0)
+    landmarks_affine(small_series, poor, template)
+
+    assert "registration may be poor" in capsys.readouterr().out
+
+
+def test_landmark_residuals_measured_in_template_pixels():
+    """Residuals are distances in the registered frame, independent of the recording's scale."""
+    template_points = np.array([[10.0, 10.0], [30.0, 10.0], [10.0, 25.0]])
+    # Recording is 4x the template, so a landmark off by 4 recording px is off by 1 template px.
+    tform = trf.AffineTransform(scale=(4.0, 4.0))
+    recording_points = tform(template_points)
+    recording_points[2] += (4.0, 0.0)
+
+    residuals = trf_mod.landmark_residuals(tform, template_points, recording_points)
+
+    np.testing.assert_allclose(residuals, [0.0, 0.0, 1.0], atol=1e-9)
+
+
+def test_registered_h5_stores_aligned_landmark_qa(preproc_h5, output_dir):
+    """Every landmark QA dataset must correspond row for row."""
+    points = str(resources.files(res).joinpath("ccf_template_landmarks_140x142.csv"))
+    runner = CliRunner()
+
+    result = runner.invoke(mesoscopy.cli, args=f"register landmarks {preproc_h5} -r {points} -o {output_dir}")
+    assert result.exit_code == 0, result.output
+
+    with h5py.File(pathlib.Path(output_dir) / "preproc_registered.h5", "r") as f:
+        names = [name.decode() for name in f["/qa/landmark_names"][:]]
+        template_points = f["/qa/template_landmarks"][:]
+        recording_points = f["/qa/recording_landmarks"][:]
+        registered_points = f["/qa/registered_landmarks"][:]
+        residuals = f["/qa/landmark_residuals"][:]
+
+    assert names == list(res.get_default_landmarks())
+    for array in (template_points, recording_points, registered_points, residuals):
+        assert len(array) == len(names)
+
+    # Recording and template points are the same file here, so the fit is exact.
+    np.testing.assert_allclose(residuals, 0.0, atol=1e-9)
+    np.testing.assert_allclose(registered_points, template_points, atol=1e-9)
+
+
 def test_update_nwb(nwbfile, preproc_h5):
     tform_params = trf.AffineTransform(scale=(1.2, 1.1), rotation=0.05, translation=(3, -4)).params
 

From 956f2efa54750ee77abcd6bee35a3cc4d69735f2 Mon Sep 17 00:00:00 2001
From: Constantinos Eleftheriou <constantinos.eleftheriou@ed.ac.uk>
Date: Thu, 13 Aug 2026 16:10:36 +0100
Subject: [PATCH 09/11] performance improvements, preallocate output array
 instead of building a li st and stacking

---
 src/mesoscopy/register/transform.py | 19 ++++++++++++---
 tests/test_registration.py          | 38 +++++++++++++++++++++++++++++
 2 files changed, 53 insertions(+), 4 deletions(-)

diff --git a/src/mesoscopy/register/transform.py b/src/mesoscopy/register/transform.py
index 0fd539a..e282bdb 100644
--- a/src/mesoscopy/register/transform.py
+++ b/src/mesoscopy/register/transform.py
@@ -190,12 +190,22 @@ def landmarks_affine(
     click.echo(f"Transform estimated in {end - start} s")
 
     n_frames = deltaf_series.shape[0]
+    if n_frames == 0:
+        msg = "The DeltaF/F series contains no frames."
+        raise ValueError(msg)
 
     def _warp_frame(idx: int) -> np.ndarray:
         return trf.warp(deltaf_series[idx], tform, order=3, output_shape=output_shape)
 
     start = time.time()
-    results: list[np.ndarray | None] = [None] * n_frames
+
+    first_frame = _warp_frame(0)
+    registered = np.empty((n_frames, *first_frame.shape), dtype=first_frame.dtype)
+    registered[0] = first_frame
+
+    def _warp_frame_into_output(idx: int) -> None:
+        registered[idx] = _warp_frame(idx)
+
     n_workers = os.cpu_count() or 1
     with (
         click.progressbar(
@@ -204,11 +214,12 @@ def _warp_frame(idx: int) -> np.ndarray:
         ) as bar,
         ThreadPoolExecutor(max_workers=n_workers) as executor,
     ):
-        futures = {executor.submit(_warp_frame, i): i for i in range(n_frames)}
+        bar.update(1)  # the first frame is already warped
+        futures = [executor.submit(_warp_frame_into_output, idx) for idx in range(1, n_frames)]
         for future in as_completed(futures):
-            results[futures[future]] = future.result()
+            future.result()  # re-raise anything a worker hit
             bar.update(1)
     end = time.time()
     click.echo(f"Recording registered in {end - start} s")
 
-    return np.stack(results), tform
+    return registered, tform
diff --git a/tests/test_registration.py b/tests/test_registration.py
index 35be3c3..2aac8fa 100644
--- a/tests/test_registration.py
+++ b/tests/test_registration.py
@@ -122,6 +122,44 @@ def test_landmarks_affine_output_shape_does_not_crop_the_source(landmark_pair):
     assert warped[0].sum() > 0.9 * series[0].sum()
 
 
+def test_landmarks_affine_preserves_frame_dtype(landmark_pair):
+    """Frames are written into a preallocated array, so its dtype must match what warp produces."""
+    recording_lm, template_lm = landmark_pair
+
+    for dtype in (np.float32, np.float64):
+        series = np.random.default_rng(0).random((4, 40, 40)).astype(dtype)
+        warped, _ = landmarks_affine(series, recording_lm, template_lm, output_shape=(40, 40))
+
+        expected_tform = trf.estimate_transform(
+            "affine",
+            np.array(list(template_lm.values())),
+            np.array(list(recording_lm.values())),
+        )
+        reference = trf.warp(series[0], expected_tform, order=3, output_shape=(40, 40))
+
+        assert warped.dtype == reference.dtype
+        np.testing.assert_array_equal(warped[0], reference)
+
+
+def test_landmarks_affine_writes_every_frame(landmark_pair):
+    """Each frame must land at its own index - workers complete out of order."""
+    recording_lm, template_lm = landmark_pair
+    # Each frame is a distinct constant, so a misplaced frame is immediately visible.
+    series = np.stack([np.full((40, 40), value, dtype=np.float32) for value in range(1, 21)])
+
+    warped, _ = landmarks_affine(series, recording_lm, template_lm, output_shape=(40, 40))
+
+    # The landmarks are a pure translation, so an interior pixel keeps its frame's value.
+    np.testing.assert_allclose(warped[:, 20, 20], np.arange(1, 21), atol=1e-4)
+
+
+def test_landmarks_affine_rejects_an_empty_series(landmark_pair):
+    recording_lm, template_lm = landmark_pair
+
+    with pytest.raises(ValueError, match="no frames"):
+        landmarks_affine(np.zeros((0, 40, 40), dtype=np.float32), recording_lm, template_lm)
+
+
 def test_landmarks_affine_identity_landmarks(small_series):
     """When recording and template landmarks are identical the warp is a near-identity."""
     landmarks = {"A": (5.0, 5.0), "B": (35.0, 5.0), "C": (5.0, 35.0)}

From d4467e80486165ca7177216d2583a93b8b325fd4 Mon Sep 17 00:00:00 2001
From: Constantinos Eleftheriou <constantinos.eleftheriou@ed.ac.uk>
Date: Thu, 13 Aug 2026 16:15:29 +0100
Subject: [PATCH 10/11] update docs

---
 docs/typical-workflow.md | 36 ++++++++++++++++++++++++++++++++++--
 1 file changed, 34 insertions(+), 2 deletions(-)

diff --git a/docs/typical-workflow.md b/docs/typical-workflow.md
index 412c5bb..e16a719 100644
--- a/docs/typical-workflow.md
+++ b/docs/typical-workflow.md
@@ -23,12 +23,44 @@ mesoscopy preprocess /path/to/example-recording.nwb
 
 ## Register to the Allen Brain Atlas
 
+First mark the anatomical landmarks on the recording. This opens the napari landmark GUI, seeded
+with a point per landmark: drag each one onto its anatomical location, then press **Save and Close**.
+
 ```bash
-mesoscopy register mark-landmarks /path/to/example-recording.nwb
+mesoscopy register label /path/to/example-recording.nwb
+```
+
+This writes `example-recording_landmarks.csv` to the output directory (`-o`, the current directory
+by default), holding each landmark's `(x, y)` position in the pixel space of the ∆F/F series.
+
+Then warp the recording onto the atlas. The landmarks file is found automatically if it sits next to
+the recording or in the output directory; pass `-r/--recording-points` to point at it explicitly.
+
+```bash
+mesoscopy register landmarks /path/to/example-recording.nwb
+```
+
+The registered frames are written in Allen CCF template space at the atlas's own dimensions, which
+is what `mesoscopy process area-responses` expects. Use `--output-width` / `--output-height` only if
+you are registering onto a different template.
+
+!!! note
+    `-t/--template-points` supplies the *template* landmarks being registered onto, not your
+    recording's landmarks. Leave it unset to use the Allen CCF landmarks that ship with mesoscopy.
+
+Registration reports how far each landmark ends up from its template position, and warns if the fit
+is poor:
+
+```
+Estimating transform from 9 landmarks...
+Landmark fit: RMSE 2.28 px, worst is 'rFP' at 3.94 px (in template pixels).
 ```
 
+To check the alignment visually, generate the QA report for the registered HDF5 file written by the
+step above — its path is echoed as `Saved registered frames at ...`:
+
 ```bash
-mesoscopy register landmarks --template-points example-recording_landmarks.csv /path/to/example-recording.nwb
+mesoscopy report /path/to/example-recording_registered.h5
 ```
 
 ## Extract area responses

From cdc4b48ebdeed43575f2d57917cf8a2b913897e6 Mon Sep 17 00:00:00 2001
From: Constantinos Eleftheriou <constantinos.eleftheriou@ed.ac.uk>
Date: Thu, 13 Aug 2026 16:19:11 +0100
Subject: [PATCH 11/11] fix docs generation

---
 docs/index.md                        | 2 +-
 src/mesoscopy/io.py                  | 6 ++++--
 src/mesoscopy/preprocess/__init__.py | 8 ++++----
 3 files changed, 9 insertions(+), 7 deletions(-)

diff --git a/docs/index.md b/docs/index.md
index bc6d6a5..2174801 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -1,4 +1,4 @@
-![https://www.mesoscopy.org](../assets/mesoscopy-logo-banner.png)
+![https://www.mesoscopy.org](assets/mesoscopy-logo-banner.png)
 
 ---
 
diff --git a/src/mesoscopy/io.py b/src/mesoscopy/io.py
index 592ffc0..6daf35c 100644
--- a/src/mesoscopy/io.py
+++ b/src/mesoscopy/io.py
@@ -59,14 +59,16 @@ def read_nwb(path: str, mode: str = "a", return_io: bool = False) -> NWBFile | t
     return nwbfile
 
 
-def write_nwb(path: str, nwbfile: NWBFile, mode: str = "w", io: NWBHDF5IO = None, **kwargs) -> None:
+def write_nwb(path: str, nwbfile: NWBFile, mode: str = "w", io: NWBHDF5IO = None, **kwargs: typing.Any) -> None:
     """Write an NWB file.
 
     Args:
         path (str): Path to the NWB file.
         nwbfile (NWBFile): NWB file object.
         mode (str, optional): File write mode (i.e. write/append). Defaults to "w".
-        **kwargs: Parameters passed to NWBHDF5IO.write.
+        io (NWBHDF5IO, optional): An already open IO object to write through. When given, the file
+            is not reopened and `path` and `mode` are ignored. Defaults to None.
+        **kwargs (typing.Any): Parameters passed to NWBHDF5IO.write.
     """
     if io:
         return io.write(nwbfile, **kwargs)
diff --git a/src/mesoscopy/preprocess/__init__.py b/src/mesoscopy/preprocess/__init__.py
index dcb15d3..d86334d 100644
--- a/src/mesoscopy/preprocess/__init__.py
+++ b/src/mesoscopy/preprocess/__init__.py
@@ -127,10 +127,10 @@ def run_preprocessing(
         chunks (int, optional): Number of chunks to load in memory. Defaults to 100.
         crop (int, optional): Number of pixels to crop from the edges of the recording. Defaults to 0.
         bins (int, optional): Recording pixel binning factor. Defaults to 2.
-        channel_means_only (bool, optional): Extract the channel means and exit without extracting a delta F series.
-        Defaults to False.
-        use_means (bool, optional): Use means histogram instead of standard deviation to separate channels.
-        Defaults to False.
+        channel_means_only (bool, optional): Extract the channel means and exit without extracting a
+            delta F series. Defaults to False.
+        use_means (bool, optional): Use means histogram instead of standard deviation to separate
+            channels. Defaults to False.
         flip_channels (bool, optional): Flip extracted channel order. Defaults to False.
         interim_dir (str, optional): Path to the interim directory. Defaults to "interim/".
         skip_start (int, optional): Number of frames to skip at the start of the recording. Defaults to None.