Registration API

Core API

register_cmd()

Register recordings to a template.

Source code in src/mesoscopy/register/__init__.py

@click.group("register")
def register_cmd() -> None:
    """Register recordings to a template."""

label_cmd(path, out_dir, template_points, session_id)

Mark landmarks on a recording for registration to a template using the landmarks GUI.

Parameters:

• path (str) –

  Path to preprocessed HDF5 file or NWB file.

• out_dir (str) –

  Output directory for registration landmarks file.

• template_points (str) –

  Path to template landmark points in CSV or Fiji XML points format.

• session_id (str) –

  Session ID for the recording.

Returns:

• dict ( 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, i.e. (column, row).

Source code in src/mesoscopy/register/__init__.py

@register_cmd.command("label")
@click.argument(
    "path",
    type=click.Path(exists=True),
)
@click.option(
    "-o",
    "--out_dir",
    type=click.Path(dir_okay=True),
    default="./",
    help="Output directory for registered recording.",
)
@click.option(
    "-t",
    "--template-points",
    type=click.Path(dir_okay=False),
    help="Path to template landmark points in CSV or Fiji XML points format",
)
@click.option(
    "--session-id",
    type=str,
    help="Session ID for the recording.",
)
def label_cmd(path, out_dir, template_points, session_id) -> dict:
    """Mark landmarks on a recording for registration to a template using the landmarks GUI.

    Args:
        path (str): Path to preprocessed HDF5 file or NWB file.
        out_dir (str): Output directory for registration landmarks file.
        template_points (str): Path to template landmark points in CSV or Fiji XML points format.
        session_id (str): Session ID for the recording.

    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, i.e. (column, row).
    """
    click.echo("Loading imaging data...")
    nwb = bool(path.endswith(".nwb"))

    if not session_id:
        session_id = session_id_from_path(path)

    maxip = None
    isosb_maxip = None

    # 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"):
            _, 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()
    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, template_shape=template_shape)

    click.echo("Saving recording landmarks...")
    outpath = out_dir + os.sep + session_id + "_landmarks.csv"
    io.write_points(outpath, recording_landmarks)

    click.echo(f"Recording landmarks saved at {outpath}.")

    return recording_landmarks

landmarks_cmd(path, out_dir, recording_points, template_points, output_width=None, output_height=None)

Register a recording to a template based on defined landmarks.

Parameters:

• path (str) –

  Path to preprocessed recording HDF5 or NWB file.

• out_dir (str) –

  Output directory for registered recording.

• recording_points (str) –

  Path to recording landmark points in CSV or Fiji XML points format.

• template_points (str) –

  Path to template landmark points in CSV or Fiji XML points format.

• output_width (int, default: None ) –

  Width of the registered frames. Defaults to the Allen CCF template width.

• output_height (int, default: None ) –

  Height of the registered frames. Defaults to the Allen CCF template height.

Returns:

• str ( str ) –

  Path to the registered recording file.

Raises:

• ValueError –

  If the path to recording landmarks cannot be inferred.

Source code in src/mesoscopy/register/__init__.py

@register_cmd.command("landmarks")
@click.argument(
    "path",
    type=click.Path(exists=True),
)
@click.option(
    "-o",
    "--out_dir",
    type=click.Path(dir_okay=True),
    default="./",
    help="Output directory for registered recording.",
)
@click.option(
    "-r",
    "--recording-points",
    type=click.Path(dir_okay=False),
    help="Path to recording landmark points in Fiji XML points format",
)
@click.option(
    "-t",
    "--template-points",
    type=click.Path(dir_okay=False),
    help="Path to template landmark points in Fiji XML points format",
)
@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,
    output_width: int | None = None,
    output_height: int | None = None,
) -> str:
    """Register a recording to a template based on defined landmarks.

    Args:
        path (str): Path to preprocessed recording HDF5 or NWB file.
        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.
        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.

    Raises:
        ValueError: If the path to recording landmarks cannot be inferred.
    """
    click.echo(f"Registering recording {path} to template.")

    os.makedirs(out_dir, exist_ok=True)

    click.echo("Loading imaging data...")

    # Determine whether we're working with an NWB file
    nwb = True if path.endswith(".nwb") else False

    session_id, deltaf_series, timestamps = io.load_deltaf(path, nwb)

    click.echo("Loading landmarks...")
    template_landmarks = res.get_default_landmarks()
    if template_points:
        template_landmarks = io.read_points(template_points)

    if not recording_points:
        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.
    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,
        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(
        path=outpath,
        data={
            "/F": warped,
            "/timestamps": timestamps,
            "/tform": tform.params,
            "/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}")

    if nwb:
        click.echo("Updating NWB file...")
        update_nwb(path, outpath, tform.params)
        click.echo(f"Updated NWB file at {path}")

    return outpath

session_id_from_path(path)

Derive a session identifier from a recording path.

Parameters:

• path (str) –

  Path to a recording file.

Returns:

• str ( str ) –

  The file name without its extension or the "_preprocessed" suffix.

Source code in src/mesoscopy/register/__init__.py

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", "")

landmarks_path_candidates(path, out_dir=None)

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.

Parameters:

• path (str) –

  Path to the recording file.

• out_dir (str, default: None ) –

  Output directory the landmarks may have been written to.

Returns:

• list[Path] –

  list[pathlib.Path]: Candidate landmark file paths, in search order, without duplicates.

Source code in src/mesoscopy/register/__init__.py

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))

load_maxips(path)

Load maximum intensity projections from a preprocessed HDF5 file.

Parameters:

• path (str) –

  Path to the preprocessed HDF5 file.

Returns:

• tuple[ndarray | None, ndarray | None] –

  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.

Source code in src/mesoscopy/register/__init__.py

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 | None, np.ndarray | None]: Maximum intensity projection for the gcamp and
            isosb channels, or (None, None) if the file holds no gcamp 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

linked_preprocessed_path(nwb_path)

Resolve the preprocessed HDF5 file that an NWB file links its dF/F series to.

Parameters:

• nwb_path (str) –

  Path to the NWB file.

Returns:

• str | None –

  str | None: Path to the linked HDF5 file, or None if the dF/F series is not an external link.

Source code in src/mesoscopy/register/__init__.py

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)

update_nwb(nwb_path, h5_path, tform_params)

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).

Parameters:

• nwb_path (str) –

  Path to the NWB file.

• h5_path (str) –

  Path to the HDF5 file containing the registered images.

• tform_params (ndarray) –

  Affine transformation matrix, as a 3x3 array.

Returns:

• NWBFile ( 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.

Source code in src/mesoscopy/register/__init__.py

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 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)

    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

Transform API

align_landmarks(recording_landmarks, template_landmarks)

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.

Parameters:

• recording_landmarks (dict) –

  Recording landmarks, as {name: (x, y)}.

• template_landmarks (dict) –

  Template landmarks, as {name: (x, y)}.

Returns:

• tuple[list[str], ndarray, ndarray] –

  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.

Source code in src/mesoscopy/register/transform.py

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

landmark_residuals(tform, template_points, recording_points)

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.

Parameters:

• tform (ProjectiveTransform) –

  The fitted transform, mapping template to recording space.

• template_points (ndarray) –

  Template points as an (n, 2) array of (x, y) coordinates.

• recording_points (ndarray) –

  Matching recording points, in the same order.

Returns:

• ndarray –

  np.ndarray: Residual distance per landmark, in template pixels.

Source code in src/mesoscopy/register/transform.py

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)

landmarks_affine(deltaf_series, recording_landmarks, template_landmarks, output_shape=None)

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.

Parameters:

• deltaf_series (ndarray) –

  DeltaF/F series.

• recording_landmarks (dict) –

  Recording landmarks, as {name: (x, y)}.

• template_landmarks (dict) –

  Template landmarks, as {name: (x, y)}.

• output_shape (tuple[int, int], default: None ) –

  Shape of the registered frames, as (height, width). Defaults to the shape of the Allen CCF atlas template.

Returns:

• tuple[ndarray, ProjectiveTransform] –

  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.

Source code in src/mesoscopy/register/transform.py

def landmarks_affine(
    deltaf_series: np.ndarray,
    recording_landmarks: dict,
    template_landmarks: dict,
    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)}.
        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.

    Raises:
        ValueError: If the landmarks do not define a usable affine transform.
    """
    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)

    names, template, recording = align_landmarks(recording_landmarks, template_landmarks)

    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]
    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()

    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(
            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,
    ):
        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):
            future.result()  # re-raise anything a worker hit
            bar.update(1)
    end = time.time()
    click.echo(f"Recording registered in {end - start} s")

    return registered, tform

Landmarks GUI

mark_landmarks(maxip_image, alt_image, template_landmarks={}, template_shape=None)

Launch the napari viewer to identify anatomical landmarks on a maximum intensity projection image.

The following landmarks are identified for registration: - bregma: Just plain ol' bregma. - cFP: Frontal pole center. - rFP: Rightmost aspect of the frontal pole. - lFP: Leftmost aspect of the frontal pole. - rPB: Right lateral edge of the parietal bone. - lPB: Left lateral edge of the parietal bone. - lpRSP: Left posterior aspect of the retrosplenial cortex. - 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.

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.

Parameters:

• maxip_image (NDArray) –

  Maximum intensity projection image. Could be either channel.

• alt_image (NDArray) –

  Alternative image to be displayed alongside the maximum intensity projection image. Usually a second channel.

• template_landmarks (dict, default: {} ) –

  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], default: None ) –

  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 ( 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.

Source code in src/mesoscopy/register/landmarks_gui.py

def mark_landmarks(
    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.

    The following landmarks are identified for registration:
    - bregma: Just plain ol' bregma.
    - cFP: Frontal pole center.
    - rFP: Rightmost aspect of the frontal pole.
    - lFP: Leftmost aspect of the frontal pole.
    - rPB: Right lateral edge of the parietal bone.
    - lPB: Left lateral edge of the parietal bone.
    - lpRSP: Left posterior aspect of the retrosplenial cortex.
    - 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.

    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, 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.
    """
    if len(maxip_image.shape) == 3:
        maxip_height, maxip_width, _ = maxip_image.shape
    else:
        maxip_height, maxip_width = maxip_image.shape

    viewer = napari.view_image(maxip_image)

    if alt_image is not None:
        if len(alt_image.shape) == 3:
            alt_height, alt_width, _ = alt_image.shape
        else:
            alt_height, alt_width = alt_image.shape
        if (maxip_height, maxip_width) != (alt_height, alt_width):
            msg = "Maximum intensity projection and alternative image must have the same dimensions."
            raise ValueError(msg)

        viewer.add_image(alt_image, name="alt_maxip")

    default_landmark_locations = template_landmarks or {
        "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())

    # 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)

    points_layer = viewer.add_points(
        data=seed_points,
        name="landmarks",
        ndim=2,
        properties={"label": landmarks},
        property_choices={"label": landmarks},
        symbol="o",
        face_color="label",
        face_color_cycle=COLOR_CYCLE,
        border_width=0,  # fraction of point size
        size=5,
    )
    points_layer.face_color_mode = "cycle"

    # add the label menu widget to the viewer
    label_widget = _create_label_menu(points_layer, landmarks)
    viewer.window.add_dock_widget(label_widget)

    points_layer.mode = "select"

    napari.run()

    return _collect_marked_points(points_layer, landmarks)