Repository: brainglobe/brainglobe-template-builder (language: python)

Issue:
implement dask-downsampling
## Description

**What is this PR**

- [ ] Bug fix
- [x] Addition of a new feature
- [ ] Other

**Why is this PR needed?**

It's useful for at least the tadpole data (and hopefully many others) preprocessing to effectively downsample data, both for initial low-res atlases (for testing), and later higher-res production-ready atlases.

**What does this PR do?**

Adds a function to downsample an image in xy, and then z, with Dask.
This could be moved to utils, maybe, in a later step.

## References

\ 

## How has this PR been tested?

Simple tests added, trying not to over-engineer at this point in the package's history.

## Is this a breaking change?

No

## Does this PR require an update to the documentation?

Docstrings added. Documenting this package on a higher level is out-of-scope for this PR

## Checklist:

- [x] The code has been tested locally
- [x] Tests have been added to cover all new functionality (unit & integration)
- [x] The documentation has been updated to reflect any changes
- [x] The code has been formatted with [pre-commit](https://pre-commit.com/)



Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
index 618521d..31b66ad 100644
--- a/.gitignore
+++ b/.gitignore
@@ -77,7 +77,7 @@ venv/
 .python-version
 
 # OS
-.DS_Store
+**/.DS_Store
 
 # written by setuptools_scm
 **/_version.py
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 9f673d0..1a04567 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -1,5 +1,5 @@
-exclude: 'conf.py'
-# Configuring https://pre-commit.ci/ bot
+# ignore everything in scripts/optimized_antsMultiVariateTemplateConstruction
+exclude: ^$|^(scripts/optimized_antsMultiVariateTemplateConstruction)
 ci:
     autoupdate_schedule: monthly
 repos:
diff --git a/MANIFEST.in b/MANIFEST.in
index 1038ec6..5637c2b 100644
--- a/MANIFEST.in
+++ b/MANIFEST.in
@@ -1,9 +1,14 @@
 include LICENSE
 include README.md
+include brainglobe_template_builder/napari/napari.yaml
 exclude *.yaml
 exclude cruft.json
 
+
 recursive-exclude * __pycache__
 recursive-exclude * *.py[co]
 recursive-exclude docs *
 recursive-exclude tests *
+recursive-exclude scripts *
+recursive-exclude examples *
+recursive-exclude data *
diff --git a/README.md b/README.md
index 3e47bee..ac0ee5c 100644
--- a/README.md
+++ b/README.md
@@ -10,6 +10,9 @@ Build unbiased anatomical templates from individual images
 
 - [Overview](#overview)
 - [Installation](#installation)
+  - [Pre-requisites](#pre-requisites)
+  - [Create a `conda` environment](#create-a-conda-environment)
+  - [Install `brainglobe-template-builder` with `pip`](#install-brainglobe-template-builder-with-pip)
 - [Background](#background)
   - [On templates and atlases](#on-templates-and-atlases)
   - [Single-subject vs population templates](#single-subject-vs-population-templates)
@@ -32,32 +35,53 @@ Build unbiased anatomical templates from individual images
 
 ## Installation
 
-> **Warning**
-> - [ANTs](http://stnava.github.io/ANTs/), which [we depend on](#template-construction-with-ants), is a large package. The installation may take a while.
-> - ANTs is only available for Linux and macOS. If you are on Windows, you can use [WSL](https://learn.microsoft.com/en-us/windows/wsl/install)
+### Pre-requisites
+
+- A Unix-like operating system (Linux or MacOS)
+- A `bash` shell (if you are using MacOS, where `zsh` is the default, you may need to install `bash` via Homebrew).
+- A working installation of [Advanced Normalisation Tools (ANTs)](http://stnava.github.io/ANTs/). We recommend following [these instructions](https://github.com/ANTsX/ANTs/wiki/Compiling-ANTs-on-Linux-and-Mac-OS) to compile ANTs from source and to follow all recommended post-installation steps.
+- The scripts from the [optimised ANTs template construction pipeline](https://github.com/CoBrALab/optimized_antsMultivariateTemplateConstruction/tree/master), which should be available in your `PATH`. According to the instructions in that repository's README, you should do the following:
+
+  ```sh
+  git clone --recursive https://github.com/CoBrALab/optimized_antsMultivariateTemplateConstruction.git
+  cd optimized_antsMultivariateTemplateConstruction
+  echo "export PATH=$(pwd):\$PATH" >> $HOME/.bashrc
+  echo "export PATH=$(pwd)/minc-toolkit-extras:\$PATH" >> $HOME/.bashrc
+  source $HOME/.bashrc  # or restart your terminal
+  ```
+
+  If you are on MacOS, you may need to replace `$HOME/.bashrc` with `$HOME/.bash_profile`.
 
+
+
+### Create a `conda` environment
 We recommend installing `brainglobe-template-builder` within a [conda](https://docs.conda.io/en/latest/) or [mamba](https://mamba.readthedocs.io/en/latest/index.html) environment. Instructions assume `conda` usage, but `mamba`/`micromamba` are interchangeable.
 
-Use the provided `environment.yaml` file to create a new environment.
 
 ```sh
-git clone https://github.com/brainglobe/brainglobe-template-builder
-cd brainglobe-template-builder
-conda env create -f environment.yaml -n template-builder
+conda env create -n template-builder -f environment.yaml
 conda activate template-builder
 ```
 
-We have called the environment `template-builder`, but you can choose any name you like.
+We have called the environment "template-builder", but you can choose any name you like.
 
-To install the latest development version of `brainglobe-template-builder`,
+This environment contains all dependencies for running the optimised ANTs template construction pipeline, but if you want to use the pre- and
+post-processing funcitonalities of `brainglobe-template-builder`, you will need to also pip install the package in editable mode (see below).
+
+### Install `brainglobe-template-builder` with `pip`
+
+
+To install the latest development version of `brainglobe-template-builder`, first clone the repository:
 
 ```sh
-pip install -e .[dev]
+git clone https://github.com/brainglobe/brainglobe-template-builder
+cd brainglobe-template-builder
 ```
-For zsh users (default shell on macOS):
+
+Then, install the package in editable mode with the following command:
 
 ```sh
-pip install -e '.[dev]'
+pip install -e .[dev]
 ```
 
 ## Background
diff --git a/brainglobe_template_builder/io.py b/brainglobe_template_builder/io.py
new file mode 100644
index 0000000..d88b1fc
--- /dev/null
+++ b/brainglobe_template_builder/io.py
@@ -0,0 +1,206 @@
+from pathlib import Path
+from typing import Literal
+
+import nibabel as nib
+import numpy as np
+import pandas as pd
+from brainglobe_utils.IO.image.load import load_any
+from brainglobe_utils.IO.image.save import to_tiff
+
+
+def get_unique_folder_in_dir(search_dir: Path, search_str: str) -> Path:
+    """
+    Find a folder in a directory that contains a unique string.
+
+    Parameters
+    ----------
+    search_dir : Path
+        Directory to search in
+    search_str : str
+        String to search for in folder names
+
+    Returns
+    -------
+    folder : Path
+        Path to the folder that contains the search string
+    """
+    all_folders = [x for x in search_dir.iterdir() if x.is_dir()]
+    folders_with_str = [x for x in all_folders if search_str in x.name]
+    if len(folders_with_str) == 0:
+        raise ValueError(f"No folders with {search_str} found")
+    if len(folders_with_str) > 1:
+        raise ValueError(f"Multiple folders with {search_str} found")
+    return folders_with_str[0]
+
+
+def load_tiff(tiff_path: Path):
+    """
+    Load an image to napari
+
+    Parameters
+    ----------
+    tiff_path : pathlib.Path
+        path to the tiff image
+
+    Returns
+    -------
+    image : np.ndarray
+    """
+    valid_extensions = [".tif", ".tiff"]
+    if tiff_path.suffix not in valid_extensions:
+        raise ValueError(
+            f"File extension {tiff_path.suffix} is not valid. "
+            f"Expected one of {valid_extensions}"
+        )
+    image = load_any(tiff_path.as_posix())
+    return image
+
+
+def save_3d_points_to_csv(points: np.ndarray, file_path: Path):
+    """
+    Save 3D points to a csv file
+    """
+
+    if points.shape[1] != 3:
+        raise ValueError(
+            f"Points must be of shape (n, 3). Got shape {points.shape}"
+        )
+    if file_path.suffix != ".csv":
+        raise ValueError(
+            f"File extension {file_path.suffix} is not valid. "
+            f"Expected file path to end in .csv"
+        )
+
+    points_df = pd.DataFrame(points, columns=["z", "y", "x"])
+    points_df.to_csv(file_path, index=False)
+
+
+def save_nii(
+    stack: np.ndarray,
+    vox_sizes: list,
+    dest_path: Path,
+    kind: Literal["image", "label"] = "image",
+):
+    """
+    Save 3D image stack to dest_path as a nifti image.
+
+    This function assumes that the image is in the ASR orientation
+    and sets the qform and sform of the nifti header accordingly
+    (so that the image is displayed correctly in nifti viewers like ITK-SNAP).
+
+    Parameters
+    ----------
+    stack : np.ndarray
+        3D image stack
+    vox_sizes : list
+        list of voxel dimensions in mm. The order is 'x', 'y', 'z'
+    dest_path : pathlib.Path
+        path to save the nifti image
+    kind : Literal["image", "label"]
+        Whether the stack is an image or a label. If label, the dtype
+        is converted to uint8.
+    """
+    if kind == "label":
+        stack = stack.astype(np.uint8)
+
+    affine = _get_transf_matrix_from_res(vox_sizes)
+    nii_img = nib.Nifti1Image(stack, affine, dtype=stack.dtype)
+    # Set qform and sform to match axes orientation, assuming ASR
+    reorient = np.array(
+        [
+            [0, 0, -1, 0],
+            [-1, 0, 0, 0],
+            [0, -1, 0, 0],
+            [0, 0, 0, 1],
+        ]
+    )
+    new_form = reorient @ affine
+    nii_img.set_qform(new_form, code=3)
+    nii_img.set_sform(new_form, code=3)
+    # save the nifti image
+    nib.save(nii_img, dest_path.as_posix())
+
+
+def file_path_with_suffix(path: Path, suffix: str, new_ext=None) -> Path:
+    """
+    Return a new path with the given suffix added before the extension.
+
+    "suffix" the string to add before the left-most period, while
+    extension is the string after that. For example, if ``suffix=="_new"``,
+    and ``new_ext=".nii.gz"``, the output path will end with "_new.nii.gz".
+
+    Parameters
+    ----------
+    path : pathlib.Path
+        The file path to modify.
+    suffix : str
+        The suffix to add before the extension.
+    new_ext : pathlib.Path, optional
+        If given, replace the current extension with this one.
+        Should include the leading period.
+
+    Returns
+    -------
+    Path
+        The new path to the file with the given suffix.
+
+    """
+    suffixes = "".join(path.suffixes)
+    pure_stem = str(path.stem).rstrip(suffixes)
+    if new_ext is not None:
+        new_name = f"{pure_stem}{suffix}{new_ext}"
+    else:
+        new_name = f"{pure_stem}{suffix}{suffixes}"
+    return path.with_name(new_name)
+
+
+def _get_transf_matrix_from_res(vox_sizes: list) -> np.ndarray:
+    """Create transformation matrix from a dictionary of voxel dimensions.
+
+    Parameters
+    ----------
+    vox_sizes : list
+        list of voxel dimensions in mm. The order is 'x', 'y', 'z'
+
+    Returns
+    -------
+    np.ndarray
+        A (4, 4) transformation matrix with the voxel dimensions
+        on the first 3 diagonal entries.
+    """
+    transformation_matrix = np.eye(4)
+    for i in range(3):
+        transformation_matrix[i, i] = vox_sizes[i]
+    return transformation_matrix
+
+
+def tiff_to_nifti(tiff_path: Path, nifti_path: Path, vox_sizes: list):
+    """
+    Convert a tiff image to a nifti image
+
+    Parameters
+    ----------
+    tiff_path : pathlib.Path
+        path to the tiff image
+    nifti_path : pathlib.Path
+        path to save the nifti image
+    vox_sizes : list
+        list of voxel dimensions in mm. The order is 'x', 'y', 'z'
+    """
+    stack = load_any(tiff_path.as_posix())
+    save_nii(stack, vox_sizes, nifti_path)
+
+
+def nifti_to_tiff(nifti_path: Path, tiff_path: Path):
+    """
+    Convert a nifti image to a tiff image
+
+    Parameters
+    ----------
+    nifti_path : pathlib.Path
+        path to the nifti image
+    tiff_path : pathlib.Path
+        path to save the tiff image
+    """
+    stack = load_any(nifti_path.as_posix())
+    to_tiff(stack, tiff_path.as_posix())
diff --git a/brainglobe_template_builder/napari/__init__.py b/brainglobe_template_builder/napari/__init__.py
new file mode 100644
index 0000000..c6e4a8e
--- /dev/null
+++ b/brainglobe_template_builder/napari/__init__.py
@@ -0,0 +1,15 @@
+import setuptools_scm
+
+try:
+    release = setuptools_scm.get_version(root="../..", relative_to=__file__)
+    __version__ = release.split("+")[0]  # remove git hash
+except LookupError:
+    __version__ = "unknown"
+
+from brainglobe_template_builder.napari._reader import napari_get_reader
+from brainglobe_template_builder.napari._widget import PreprocWidgets
+
+__all__ = (
+    "napari_get_reader",
+    "PreprocWidgets",
+)
diff --git a/brainglobe_template_builder/napari/_reader.py b/brainglobe_template_builder/napari/_reader.py
new file mode 100644
index 0000000..19fb080
--- /dev/null
+++ b/brainglobe_template_builder/napari/_reader.py
@@ -0,0 +1,77 @@
+"""
+This module is an example of a barebones numpy reader plugin for napari.
+
+It implements the Reader specification, but your plugin may choose to
+implement multiple readers or even other plugin contributions. see:
+https://napari.org/stable/plugins/guides.html?#readers
+"""
+
+from pathlib import Path
+
+from brainglobe_utils.IO.image.load import load_any
+
+
+def napari_get_reader(path):
+    """A basic implementation of a Reader contribution.
+
+    Parameters
+    ----------
+    path : str or list of str
+        Path to file, or list of paths.
+
+    Returns
+    -------
+    function or None
+        If the path is a recognized format, return a function that accepts the
+        same path or list of paths, and returns a list of layer data tuples.
+    """
+    if isinstance(path, list):
+        # reader plugins may be handed single path, or a list of paths.
+        # if it is a list, it is assumed to be an image stack...
+        # so we are only going to look at the first file.
+        path = path[0]
+
+    # if we know we cannot read the file, we immediately return None.
+    valid_extensions = (".tif", ".tiff")
+    if Path(path).suffix not in valid_extensions:
+        print(f"Invalid file extension. Expected one of: {valid_extensions}")
+        return None
+
+    # otherwise we return the *function* that can read ``path``.
+    return reader_function
+
+
+def reader_function(path):
+    """Take a path or list of paths and return a list of LayerData tuples.
+
+    Readers are expected to return data as a list of tuples, where each tuple
+    is (data, [add_kwargs, [layer_type]]), "add_kwargs" and "layer_type" are
+    both optional.
+
+    Parameters
+    ----------
+    path : str or list of str
+        Path to file, or list of paths.
+
+    Returns
+    -------
+    layer_data : list of tuples
+        A list of LayerData tuples where each tuple in the list contains
+        (data, metadata, layer_type), where data is a numpy array, metadata is
+        a dict of keyword arguments for the corresponding viewer.add_* method
+        in napari, and layer_type is a lower-case string naming the type of
+        layer. Both "meta", and "layer_type" are optional. napari will
+        default to layer_type=="image" if not provided
+    """
+    layer_type = "image"
+    # handle both a string and a list of strings
+    paths = [path] if isinstance(path, str) else path
+    # load all files into arrays
+    arrays = [load_any(_path) for _path in paths]
+    # Specify names for the layers
+    names = [Path(_path).stem for _path in paths]
+    # return a list of tuples of the form (data, add_kwargs, layer_type)
+    return [
+        (array, {"name": name}, layer_type)
+        for array, name in zip(arrays, names)
+    ]
diff --git a/brainglobe_template_builder/napari/_widget.py b/brainglobe_template_builder/napari/_widget.py
new file mode 100644
index 0000000..c120772
--- /dev/null
+++ b/brainglobe_template_builder/napari/_widget.py
@@ -0,0 +1,49 @@
+from brainglobe_utils.qtpy.collapsible_widget import CollapsibleWidgetContainer
+from napari.viewer import Viewer
+
+from brainglobe_template_builder.napari.align_widget import AlignMidplane
+from brainglobe_template_builder.napari.mask_widget import CreateMask
+from brainglobe_template_builder.napari.reorient_widget import Reorient
+from brainglobe_template_builder.napari.save_widget import SaveFiles
+
+
+class PreprocWidgets(CollapsibleWidgetContainer):
+    def __init__(self, napari_viewer: Viewer, parent=None):
+        super().__init__()
+
+        self.add_widget(
+            Reorient(napari_viewer, parent=self),
+            collapsible=True,
+            widget_title="Reorient to standard space",
+        )
+
+        self.add_widget(
+            CreateMask(napari_viewer, parent=self),
+            collapsible=True,
+            widget_title="Create mask",
+        )
+
+        self.add_widget(
+            AlignMidplane(napari_viewer, parent=self),
+            collapsible=True,
+            widget_title="Align midplane",
+        )
+
+        self.add_widget(
+            SaveFiles(napari_viewer, parent=self),
+            collapsible=True,
+            widget_title="Save files",
+        )
+
+        (
+            self.reorient_widget,
+            self.mask_widget,
+            self.midplane_widget,
+            self.save_widget,
+        ) = self.collapsible_widgets
+        # expand first widget by default
+        self.reorient_widget.expand()
+        # refresh dropdowns when midline widget is toggled
+        self.midplane_widget.toggled.connect(
+            self.midplane_widget.content().refresh_dropdowns
+        )
diff --git a/brainglobe_template_builder/napari/align_widget.py b/brainglobe_template_builder/napari/align_widget.py
new file mode 100644
index 0000000..e1f3bb9
--- /dev/null
+++ b/brainglobe_template_builder/napari/align_widget.py
@@ -0,0 +1,228 @@
+from pathlib import Path
+
+from napari.layers import Image, Labels, Layer, Points
+from napari.utils.notifications import show_info
+from napari.viewer import Viewer
+from qtpy.QtWidgets import (
+    QComboBox,
+    QFileDialog,
+    QFormLayout,
+    QGroupBox,
+    QPushButton,
+    QWidget,
+)
+
+from brainglobe_template_builder.preproc import (
+    MidplaneAligner,
+    MidplaneEstimator,
+)
+
+
+class AlignMidplane(QWidget):
+    """Widget to align the plane of symmetry to the midplane of the image."""
+
+    def __init__(self, napari_viewer: Viewer, parent=None):
+        super().__init__(parent=parent)
+        self.viewer = napari_viewer
+        self.setLayout(QFormLayout())
+        self._create_estimate_group()
+        self._create_align_group()
+
+    def _create_estimate_group(self):
+        """Create the group of widgets concerned with estimating midplane
+        points."""
+        self.estimate_groupbox = QGroupBox("Estimate points along midplane")
+        self.estimate_groupbox.setLayout(QFormLayout())
+        self.layout().addRow(self.estimate_groupbox)
+
+        # Add dropdown to select labels layer (mask)
+        self.select_mask_dropdown = QComboBox(parent=self.estimate_groupbox)
+        self.select_mask_dropdown.addItems(self._get_layers_by_type(Labels))
+        self.select_mask_dropdown.currentTextChanged.connect(
+            self._on_dropdown_selection_change
+        )
+        self.estimate_groupbox.layout().addRow(
+            "mask:", self.select_mask_dropdown
+        )
+
+        # Add dropdown to select axis
+        self.select_axis_dropdown = QComboBox(parent=self.estimate_groupbox)
+        self.select_axis_dropdown.addItems(["x", "y", "z"])
+        self.estimate_groupbox.layout().addRow(
+            "symmetry axis:", self.select_axis_dropdown
+        )
+
+        # Initialise button to estimate midplane points
+        self.estimate_points_button = QPushButton(
+            "Estimate points", parent=self.estimate_groupbox
+        )
+        self.estimate_points_button.setEnabled(False)
+        self.estimate_points_button.clicked.connect(
+            self._on_estimate_button_click
+        )
+        self.estimate_groupbox.layout().addRow(self.estimate_points_button)
+
+    def _create_align_group(self):
+        """Create the group of widgets concerned with aligning the image to
+        the midplane."""
+
+        self.align_groupbox = QGroupBox("Align image to midplane")
+        self.align_groupbox.setLayout(QFormLayout())
+        self.layout().addRow(self.align_groupbox)
+
+        # Add dropdown to select image layer
+        self.select_image_dropdown = QComboBox(parent=self.align_groupbox)
+        self.select_image_dropdown.addItems(self._get_layers_by_type(Image))
+        self.select_image_dropdown.currentTextChanged.connect(
+            self._on_dropdown_selection_change
+        )
+        self.align_groupbox.layout().addRow(
+            "image:", self.select_image_dropdown
+        )
+
+        # Add dropdown to select points layer
+        self.select_points_dropdown = QComboBox(parent=self.align_groupbox)
+        self.select_points_dropdown.addItems(self._get_layers_by_type(Points))
+        self.select_points_dropdown.currentTextChanged.connect(
+            self._on_dropdown_selection_change
+        )
+        self.align_groupbox.layout().addRow(
+            "points:", self.select_points_dropdown
+        )
+
+        # Add button to align image to midplane
+        self.align_image_button = QPushButton(
+            "Align image", parent=self.align_groupbox
+        )
+        self.align_image_button.setEnabled(False)
+        self.align_image_button.clicked.connect(self._on_align_button_click)
+        self.align_groupbox.layout().addRow(self.align_image_button)
+
+        # Add button to save transformed image
+        self.save_transform_button = QPushButton(
+            "Save transform", parent=self.align_groupbox
+        )
+        self.save_transform_button.setEnabled(False)
+        self.save_transform_button.clicked.connect(
+            self._on_save_transform_click
+        )
+        self.align_groupbox.layout().addRow(self.save_transform_button)
+
+    def _get_layers_by_type(self, layer_type: Layer) -> list:
+        """Return a list of napari layers of a given type."""
+        return [
+            layer.name
+            for layer in self.viewer.layers
+            if isinstance(layer, layer_type)
+        ]
+
+    def refresh_dropdowns(self):
+        """Refresh the dropdowns to reflect the current layers."""
+        for layer_type, dropdown in zip(
+            [Labels, Image, Points],
+            [
+                self.select_mask_dropdown,
+                self.select_image_dropdown,
+                self.select_points_dropdown,
+            ],
+        ):
+            dropdown.clear()
+            dropdown.addItems(self._get_layers_by_type(layer_type))
+
+    def _on_estimate_button_click(self):
+        """Estimate midplane points and add them to the viewer."""
+        # Estimate 9 midplane points based on the selected mask
+        mask_name = self.select_mask_dropdown.currentText()
+        mask = self.viewer.layers[mask_name]
+        axis = self.select_axis_dropdown.currentText()
+        estimator = MidplaneEstimator(mask.data, symmetry_axis=axis)
+        points = estimator.get_points()
+        prefix = mask_name.split("_label-")[0]
+        points_name = f"{prefix}_points-midplane"
+
+        # Point layer attributes
+        point_attrs = {
+            "properties": {"label": list(range(9))},
+            "face_color": "green",
+            "symbol": "cross",
+            "edge_width": 0,
+            "opacity": 0.6,
+            "size": 6,
+            "ndim": mask.ndim,
+            "name": points_name,
+        }
+
+        self.viewer.add_points(points, **point_attrs)
+        self.refresh_dropdowns()
+        # Move viewer to show z-plane of first point
+        self.viewer.dims.set_point(0, points[0][0])
+        # Enable "Select points" mode
+        self.viewer.layers[points_name].mode = "select"
+        show_info("Please move all 9 estimated points exactly to the midplane")
+
+    def _on_align_button_click(self):
+        """Align image and mask to midplane and add them to the viewer."""
+        image_name = self.select_image_dropdown.currentText()
+        image_data = self.viewer.layers[image_name].data
+        mask_name = self.select_mask_dropdown.currentText()
+        mask_data = self.viewer.layers[mask_name].data
+        points_name = self.select_points_dropdown.currentText()
+        points_data = self.viewer.layers[points_name].data
+        axis = self.select_axis_dropdown.currentText()
+
+        aligner = MidplaneAligner(
+            image_data,
+            points_data,
+            symmetry_axis=axis,
+        )
+        aligned_image = aligner.transform_image(image_data)
+        aligned_image_name = f"{image_name}_aligned"
+        self.viewer.add_image(aligned_image, name=aligned_image_name)
+        aligned_mask = aligner.transform_image(mask_data)
+        aligned_mask_name = f"{mask_name}_aligned"
+        self.viewer.add_labels(
+            aligned_mask, name=aligned_mask_name, opacity=0.5
+        )
+        aligned_halves = aligner.label_halves(aligned_image)
+        aligned_halves_name = aligned_mask_name.replace("-brain", "-halves")
+        self.viewer.add_labels(
+            aligned_halves, name=aligned_halves_name, opacity=0.5
+        )
+        # Hide original image, mask, and points layers
+        self.viewer.layers[image_name].visible = False
+        self.viewer.layers[mask_name].visible = False
+        self.viewer.layers[points_name].visible = False
+        # Hide aligned mask layer
+        self.viewer.layers[aligned_mask_name].visible = False
+        # Make aligner object accessible to other methods
+        self.aligner = aligner
+        # Enable save transform button
+        self.save_transform_button.setEnabled(True)
+
+    def _on_dropdown_selection_change(self):
+        # Enable estimate button if mask dropdown has a value
+        if self.select_mask_dropdown.currentText() == "":
+            self.estimate_points_button.setEnabled(False)
+        else:
+            self.estimate_points_button.setEnabled(True)
+
+        # Enable align button if both image and points dropdowns have a value
+        if (
+            self.select_image_dropdown.currentText() == ""
+            or self.select_points_dropdown.currentText() == ""
+        ):
+            self.align_image_button.setEnabled(False)
+        else:
+            self.align_image_button.setEnabled(True)
+
+    def _on_save_transform_click(self):
+        """Save the midplane alignment tranform to a text file."""
+        if not hasattr(self, "aligner"):
+            show_info("Please align the image to the midplane first")
+            return
+        dlg = QFileDialog()
+        dlg.setFileMode(QFileDialog.AnyFile)
+        dlg.setAcceptMode(QFileDialog.AcceptSave)
+        if dlg.exec_():
+            path = dlg.selectedFiles()[0]
+            self.aligner.save_transform(Path(path))
diff --git a/brainglobe_template_builder/napari/mask_widget.py b/brainglobe_template_builder/napari/mask_widget.py
new file mode 100644
index 0000000..e463e38
--- /dev/null
+++ b/brainglobe_template_builder/napari/mask_widget.py
@@ -0,0 +1,75 @@
+from napari.layers import Image
+from napari.utils.notifications import show_info
+from napari.viewer import Viewer
+from qtpy.QtWidgets import (
+    QComboBox,
+    QFormLayout,
+    QGroupBox,
+    QPushButton,
+    QSpinBox,
+    QWidget,
+)
+
+from brainglobe_template_builder.preproc import create_mask
+
+
+class CreateMask(QWidget):
+    """Widget to create a mask from a selected image layer."""
+
+    def __init__(self, napari_viewer: Viewer, parent=None):
+        super().__init__(parent=parent)
+        self.viewer = napari_viewer
+        self.setLayout(QFormLayout())
+
+        self._create_mask_group()
+
+    def _create_mask_group(self):
+        """Create the group of widgets concerned with creating a mask."""
+        self.mask_groupbox = QGroupBox("Create mask to exclude background")
+        self.mask_groupbox.setLayout(QFormLayout())
+        self.layout().addRow(self.mask_groupbox)
+
+        self.gauss_sigma = QSpinBox(parent=self.mask_groupbox)
+        self.gauss_sigma.setRange(0, 20)
+        self.gauss_sigma.setValue(3)
+        self.mask_groupbox.layout().addRow("gaussian sigma:", self.gauss_sigma)
+
+        self.threshold_method = QComboBox(parent=self.mask_groupbox)
+        self.threshold_method.addItems(["triangle", "otsu", "isodata"])
+        self.mask_groupbox.layout().addRow(
+            "threshold method:", self.threshold_method
+        )
+
+        self.closing_size = QSpinBox(parent=self.mask_groupbox)
+        self.closing_size.setRange(0, 20)
+        self.closing_size.setValue(5)
+        self.mask_groupbox.layout().addRow("closing size:", self.closing_size)
+
+        self.create_mask_button = QPushButton("Create mask", parent=self)
+        self.mask_groupbox.layout().addRow(self.create_mask_button)
+        self.create_mask_button.clicked.connect(self._on_button_click)
+
+    def _on_button_click(self):
+        """Create a mask from the selected image layer, using the parameters
+        specified in the widget, and add it to the napari viewer.
+        """
+
+        if len(self.viewer.layers.selection) != 1:
+            show_info("Please select exactly one Image layer")
+            return None
+
+        image = list(self.viewer.layers.selection)[0]
+
+        if not isinstance(image, Image):
+            show_info("The selected layer is not an Image layer")
+            return None
+
+        mask_data = create_mask(
+            image.data,
+            gauss_sigma=self.gauss_sigma.value(),
+            threshold_method=self.threshold_method.currentText(),
+            closing_size=self.closing_size.value(),
+        )
+
+        mask_name = f"{image.name}_label-brain"
+        self.viewer.add_labels(mask_data, name=mask_name, opacity=0.5)
diff --git a/brainglobe_template_builder/napari/napari.yaml b/brainglobe_template_builder/napari/napari.yaml
new file mode 100644
index 0000000..5282ed8
--- /dev/null
+++ b/brainglobe_template_builder/napari/napari.yaml
@@ -0,0 +1,19 @@
+name: brainglobe-template-builder
+display_name: brainglobe-template-builder
+contributions:
+  commands:
+    - id: brainglobe-template-builder.get_reader
+      python_name: brainglobe_template_builder.napari._reader:napari_get_reader
+      title: Open data with brainglobe-template-builder
+    - id: brainglobe-template-builder.make_preproc_widget
+      python_name: brainglobe_template_builder.napari._widget:PreprocWidgets
+      title: Preprocess
+  readers:
+    - command: brainglobe-template-builder.get_reader
+      accepts_directories: false
+      filename_patterns:
+      - '*.tif'
+      - '*.tiff'
+  widgets:
+    - command: brainglobe-template-builder.make_preproc_widget
+      display_name: Preprocess
diff --git a/brainglobe_template_builder/napari/reorient_widget.py b/brainglobe_template_builder/napari/reorient_widget.py
new file mode 100644
index 0000000..ed090e7
--- /dev/null
+++ b/brainglobe_template_builder/napari/reorient_widget.py
@@ -0,0 +1,128 @@
+from brainglobe_space import AnatomicalSpace
+from napari.layers import Image, Labels, Points
+from napari.utils.notifications import show_info
+from napari.viewer import Viewer
+from qtpy.QtWidgets import (
+    QComboBox,
+    QFormLayout,
+    QGroupBox,
+    QHBoxLayout,
+    QPushButton,
+    QWidget,
+)
+
+
+class Reorient(QWidget):
+    """Widget to reorient images to a standard anatomical space."""
+
+    def __init__(self, napari_viewer: Viewer, parent=None):
+        super().__init__(parent=parent)
+        self.viewer = napari_viewer
+        self.setLayout(QFormLayout())
+
+        self._create_reorient_group()
+
+    def _create_reorient_group(self):
+        self.reorient_groupbox = QGroupBox("Reorient to standard space")
+        self.reorient_groupbox.setLayout(QFormLayout())
+        self.layout().addRow(self.reorient_groupbox)
+
+        avail_labels = AnatomicalSpace.lims_labels
+
+        # Enter labels for source space
+        source_origin_row = QHBoxLayout()
+        self.source_origin1 = QComboBox(parent=self.reorient_groupbox)
+        self.source_origin1.addItems(avail_labels)
+        self.source_origin1.setCurrentIndex(0)
+        source_origin_row.addWidget(self.source_origin1)
+        self.source_origin2 = QComboBox(parent=self.reorient_groupbox)
+        self.source_origin2.addItems(avail_labels)
+        self.source_origin2.setCurrentIndex(2)
+        source_origin_row.addWidget(self.source_origin2)
+        self.source_origin3 = QComboBox(parent=self.reorient_groupbox)
+        self.source_origin3.addItems(avail_labels)
+        self.source_origin3.setCurrentIndex(5)
+        source_origin_row.addWidget(self.source_origin3)
+        self.reorient_groupbox.layout().addRow(
+            "Source origin:", source_origin_row
+        )
+
+        # Enter labels for target space
+        target_origin_row = QHBoxLayout()
+        self.target_origin1 = QComboBox(parent=self.reorient_groupbox)
+        self.target_origin1.addItems(avail_labels)
+        self.target_origin1.setCurrentIndex(1)
+        target_origin_row.addWidget(self.target_origin1)
+        self.target_origin2 = QComboBox(parent=self.reorient_groupbox)
+        self.target_origin2.addItems(avail_labels)
+        self.target_origin2.setCurrentIndex(2)
+        target_origin_row.addWidget(self.target_origin2)
+        self.target_origin3 = QComboBox(parent=self.reorient_groupbox)
+        self.target_origin3.addItems(avail_labels)
+        self.target_origin3.setCurrentIndex(5)
+        target_origin_row.addWidget(self.target_origin3)
+        self.reorient_groupbox.layout().addRow(
+            "Target origin:", target_origin_row
+        )
+
+        self.reorient_button = QPushButton(
+            "Reorient selected layers", parent=self.reorient_groupbox
+        )
+        self.reorient_button.clicked.connect(self.reorient_layers)
+        self.reorient_groupbox.layout().addRow(self.reorient_button)
+
+    def reorient_layers(self):
+        """Reorient selected layers to the target space."""
+        selected_layers = self.viewer.layers.selection
+        selected_stacks = [
+            layer
+            for layer in selected_layers
+            if isinstance(layer, Image) or isinstance(layer, Labels)
+        ]
+        selected_points = [
+            layer for layer in selected_layers if isinstance(layer, Points)
+        ]
+
+        if not selected_layers:
+            show_info("No layers selected. Select layers to reorient.")
+            return
+
+        source_origin = [
+            self.source_origin1.currentText(),
+            self.source_origin2.currentText(),
+            self.source_origin3.currentText(),
+        ]
+        target_origin = [
+            self.target_origin1.currentText(),
+            self.target_origin2.currentText(),
+            self.target_origin3.currentText(),
+        ]
+        source_origin = "".join(source_origin)
+        target_origin = "".join(target_origin)
+
+        for layer in selected_points:
+            if not selected_stacks:
+                show_info(
+                    "Cannot reorient Points layer alone. Please select "
+                    "at least one Image or Labels layer along with the "
+                    "Points layer."
+                )
+                return
+            source_stack_shape = selected_stacks[0].data.shape
+            source_space = AnatomicalSpace(
+                source_origin, shape=source_stack_shape
+            )
+            layer.data = source_space.map_points_to(target_origin, layer.data)
+            layer.name = f"{layer.name}_space-{target_origin}"
+
+        for layer in selected_stacks:
+            interp_order = 0 if isinstance(layer, Labels) else 3
+            source_space = AnatomicalSpace(
+                source_origin, shape=layer.data.shape
+            )
+            layer.data = source_space.map_stack_to(
+                target_origin, layer.data, interp_order=interp_order
+            )
+            layer.name = f"{layer.name}_orig-{target_origin}"
+
+        show_info("Layers reoriented to target space.")
diff --git a/brainglobe_template_builder/napari/save_widget.py b/brainglobe_template_builder/napari/save_widget.py
new file mode 100644
index 0000000..15d2503
--- /dev/null
+++ b/brainglobe_template_builder/napari/save_widget.py
@@ -0,0 +1,126 @@
+from pathlib import Path
+
+from napari.layers import Image, Labels, Points
+from napari.utils.notifications import show_info
+from napari.viewer import Viewer
+from qtpy.QtWidgets import (
+    QFileDialog,
+    QFormLayout,
+    QGroupBox,
+    QHBoxLayout,
+    QLineEdit,
+    QPushButton,
+    QWidget,
+)
+
+from brainglobe_template_builder.io import save_3d_points_to_csv, save_nii
+
+
+class SaveFiles(QWidget):
+    def __init__(self, napari_viewer: Viewer, parent=None):
+        super().__init__(parent=parent)
+        self.viewer = napari_viewer
+        self.setLayout(QFormLayout())
+
+        self._create_save_group()
+
+    def _create_save_group(self):
+        self.save_groupbox = QGroupBox("Save to output directory")
+        self.save_groupbox.setLayout(QFormLayout())
+        self.layout().addRow(self.save_groupbox)
+
+        self._create_voxel_size_widget()
+        self._create_save_path_widget()
+
+        self.save_button = QPushButton(
+            "Save selected layers", parent=self.save_groupbox
+        )
+        self.save_button.clicked.connect(self.save_selected_layers)
+        self.save_groupbox.layout().addRow(self.save_button)
+
+    def _create_voxel_size_widget(self):
+        """Create 3 fields for entering the voxel size."""
+        self.voxel_size_layout = QHBoxLayout()
+        self.z_size = QLineEdit()
+        self.y_size = QLineEdit()
+        self.x_size = QLineEdit()
+        self.z_size.setText("1")
+        self.y_size.setText("1")
+        self.x_size.setText("1")
+        self.voxel_size_layout.addWidget(self.z_size)
+        self.voxel_size_layout.addWidget(self.y_size)
+        self.voxel_size_layout.addWidget(self.x_size)
+
+        self.save_groupbox.layout().addRow(
+            "Voxel size (zyx) in mm:", self.voxel_size_layout
+        )
+
+    def _create_save_path_widget(self):
+        """Create a line edit and browse button for selecting a save path.
+
+        The path has to be a valid directory path.
+        """
+        self.path_edit = QLineEdit()
+        self.browse_button = QPushButton("browse")
+        self.browse_button.clicked.connect(self._open_save_dialog)
+
+        self.path_layout = QHBoxLayout()
+        self.path_layout.addWidget(self.path_edit)
+        self.path_layout.addWidget(self.browse_button)
+        self.save_groupbox.layout().addRow(
+            "Output directory:", self.path_layout
+        )
+
+    def _open_save_dialog(self):
+        """Select an existing directory path to save files to."""
+        dlg = QFileDialog()
+        dlg.setFileMode(QFileDialog.Directory)
+        dlg.setOption(QFileDialog.ShowDirsOnly, True)
+        dlg.AcceptMode(QFileDialog.AcceptSave)
+        if dlg.exec_():
+            path = dlg.selectedFiles()[0]
+            self.path_edit.setText(path)
+
+    def save_selected_layers(self):
+        save_dir = self.path_edit.text()
+        if not save_dir:
+            return
+
+        # Get voxel sizes
+        try:
+            vox_sizes = [
+                float(self.z_size.text()),
+                float(self.y_size.text()),
+                float(self.x_size.text()),
+            ]
+        except ValueError:
+            show_info("Please enter valid voxel sizes in mm.")
+
+        selected_layers = self.viewer.layers.selection
+        saved_layer_names = []
+        for layer in selected_layers:
+            if isinstance(layer, Points):
+                csv_path = f"{save_dir}/{layer.name}.csv"
+                layer.save(csv_path)  # native napari save for points
+                df_path = Path(f"{save_dir}/{layer.name}_df.csv")
+                save_3d_points_to_csv(layer.data, df_path)  # with headers
+                saved_layer_names.append(layer.name)
+            elif isinstance(layer, Image) or isinstance(layer, Labels):
+                tif_path = f"{save_dir}/{layer.name}.tif"
+                if isinstance(layer, Image):
+                    layer.data = layer.data.astype("float32")
+                layer.save(tif_path)  # native napari save for images
+                nii_path = Path(f"{save_dir}/{layer.name}.nii.gz")
+                save_nii(  # custom save to nii
+                    layer.data,
+                    vox_sizes=vox_sizes,
+                    dest_path=nii_path,
+                )
+                saved_layer_names.append(layer.name)
+            else:
+                raise UserWarning(
+                    f"Layer {layer.name} is of type {type(layer)} and cannot "
+                    "be saved. Only Image, Labels and Points layers are saved."
+                )
+        info_msg = "Saved layers: " + ", ".join(saved_layer_names)
+        show_info(info_msg)
diff --git a/brainglobe_template_builder/preproc/__init__.py b/brainglobe_template_builder/preproc/__init__.py
new file mode 100644
index 0000000..5f6f61c
--- /dev/null
+++ b/brainglobe_template_builder/preproc/__init__.py
@@ -0,0 +1,5 @@
+from brainglobe_template_builder.preproc.masking import create_mask
+from brainglobe_template_builder.preproc.alignment import (
+    MidplaneAligner,
+    MidplaneEstimator,
+)
diff --git a/brainglobe_template_builder/preproc/alignment.py b/brainglobe_template_builder/preproc/alignment.py
new file mode 100644
index 0000000..44fbe07
--- /dev/null
+++ b/brainglobe_template_builder/preproc/alignment.py
@@ -0,0 +1,236 @@
+from itertools import product
+from pathlib import Path
+from typing import Literal
+
+import numpy as np
+from skimage import measure
+
+from brainglobe_template_builder.preproc.transform_utils import (
+    apply_transform,
+    get_rotation_from_vectors,
+)
+
+NAPARI_AXIS_ORDER = "zyx"
+
+
+class MidplaneEstimator:
+    """Class to estimates points along the midplane of a 3D object, based on a
+    binary mask of the object."""
+
+    def __init__(
+        self,
+        mask: np.ndarray,
+        symmetry_axis: Literal["x", "y", "z"] = "x",
+    ):
+        """Initialise the estimator.
+
+        Parameters
+        ----------
+        mask : np.ndarray
+            A 3D binary mask of the object.
+        symmetry_axis : str
+            Axis of symmetry, one of 'x', 'y', and 'z'.
+            Defaults to 'x'. For brains, this would be the left-right axis.
+            Keep in mind that the axis order is zyx in napari.
+        """
+        self.mask = mask
+        self.symmetry_axis = symmetry_axis
+
+        self._validate_inputs()
+        self.symmetry_axis_idx = NAPARI_AXIS_ORDER.index(symmetry_axis)
+        print("Axis index:", self.symmetry_axis_idx)
+
+    def _validate_inputs(self):
+        """Validate the inputs to the aligner."""
+
+        if self.mask.ndim != 3:
+            raise ValueError("Mask must be 3D")
+        if self.symmetry_axis not in ["x", "y", "z"]:
+            raise ValueError("Symmetry axis must be one of 'x', 'y', or 'z'")
+        try:
+            self.mask = self.mask.astype(bool)
+        except ValueError:
+            raise ValueError("Mask must be binary")
+
+    def _get_mask_properties(self):
+        """Get properties of the mask, specifically the centroid and the
+        dimensions of the mask's bounding box (in pixels)."""
+        props = measure.regionprops(measure.label(self.mask))[0]
+        self.centroid = np.array(props.centroid)
+        bbox = np.array(props.bbox).reshape(2, 3).T
+        self.mask_dims = bbox[:, 1] - bbox[:, 0]
+        print("Mask dims:", self.mask_dims)
+
+    def get_points(self):
+        """Estimate 9 points along the midplane of the mask.
+
+        The midplane is initialised as the plane perpendicular to the symmetry
+        axis that contains the centroid of the mask. The midplane is then
+        intersected with 3 slices along each of the other two axes, resulting
+        in 9 points.
+        """
+
+        self._get_mask_properties()
+        # Find middle of the symmetry axis
+        mid_plane = self.centroid[self.symmetry_axis_idx]
+        # Find slices at 1/4, 2/4, and 3/4 of the two other dimensions
+        a, b = [i for i in range(3) if i != self.symmetry_axis_idx]
+        print("Other axes:", a, b)
+        other_planes_a = [self.mask_dims[a] / 4 * i for i in [1, 2, 3]]
+        other_planes_b = [self.mask_dims[b] / 4 * i for i in [1, 2, 3]]
+        # Find points at the intersection the symmetry axis midplane
+        # with quarter-planes of the other two axes (produces 9 points)
+        points = list(product(other_planes_a, other_planes_b, [mid_plane]))
+        self.points = np.array(points)  # (9, 3)
+        return self.points
+
+
+class MidplaneAligner:
+    """Class for aligning a given plane (as defined by a set of 3D points)
+    to the midplane of an image along a given symmetry axis."""
+
+    def __init__(
+        self,
+        image: np.ndarray,
+        points: np.ndarray,
+        symmetry_axis: Literal["x", "y", "z"] = "x",
+    ):
+        """Initialise the aligner.
+
+        Parameters
+        ----------
+        image : np.ndarray
+            A 3D image to align.
+        points : np.ndarray
+            An array of shape (n_points, 3) containing 3D point coordinates.
+            At least 3 points are required, and they must not be colinear.
+        symmetry_axis : str
+            Axis of symmetry, one of 'x', 'y', and 'z'.
+            Defaults to 'x'. For brains, this would be the left-right axis.
+            Keep in mind that the axis order is zyx in napari.
+        """
+        self.image = image
+        self.points = points
+        self.symmetry_axis = symmetry_axis
+
+        self._validate_inputs()
+        self.symmetry_axis_idx = NAPARI_AXIS_ORDER.index(symmetry_axis)
+
+    def _validate_inputs(self):
+        """Validate the inputs to the aligner."""
+        if self.image.ndim != 3:
+            raise ValueError("Image must be 3D")
+        if self.points.ndim != 2 or self.points.shape[1] != 3:
+            raise ValueError("Points must be an array of shape (n_points, 3)")
+        if self.points.shape[0] < 3:
+            raise ValueError("At least 3 points are required")
+        if np.linalg.matrix_rank(self.points) < 3:
+            raise ValueError("Points must not be colinear")
+        if self.symmetry_axis not in ["x", "y", "z"]:
+            raise ValueError("Symmetry axis must be one of 'x', 'y', or 'z'")
+
+    def _fit_plane_to_points(self):
+        """Fit a plane to the points.
+
+        The plane is fitted using SVD, and the normal vector is computed as the
+        last row of the V matrix. The normal vector is then inverted if it
+        points in the opposite direction of the symmetry axis unit vector.
+        """
+        # Use SVD to get the normal vector to the plane
+        centroid = np.mean(self.points, axis=0)
+        _, _, vh = np.linalg.svd(self.points - centroid)
+        normal_vector = vh[-1]
+
+        # Unit vector along the symmetry axis
+        symmetry_axis_vector = np.zeros(3)
+        symmetry_axis_vector[self.symmetry_axis_idx] = 1
+
+        # invert the normal vector if it points in the opposite direction of
+        # symmetry axis unit vector
+        if np.dot(normal_vector, symmetry_axis_vector) < 0:
+            normal_vector = -normal_vector
+
+        self.centroid = centroid
+        self.normal_vector = normal_vector
+        self.symmetry_axis_vector = symmetry_axis_vector
+
+    def _compute_transform(self):
+        """Find the transformation matrix that aligns the plane defined by the
+        points to the midplane of the image along the symmetry axis.
+        """
+        self._fit_plane_to_points()
+
+        # Compute the necessary transforms
+        # 1. translate to origin (so that centroid is at origin)
+        translation_to_origin = np.eye(4)
+        translation_to_origin[:3, 3] = -self.centroid
+        # 2. rotate to align fitted plane with symmetry axis
+        rotation = np.eye(4)
+        rotation[:3, :3] = get_rotation_from_vectors(
+            self.normal_vector, self.symmetry_axis_vector
+        )
+        # 3. translate to mid-axis (so that centroid is at middle of axis)
+        translation_to_mid_axis = np.eye(4)
+        offset = (
+            self.image.shape[self.symmetry_axis_idx] / 2
+            - self.centroid[self.symmetry_axis_idx]
+        ) * self.symmetry_axis_vector
+        translation_to_mid_axis[:3, 3] = self.centroid + offset
+        # Combine the transforms
+        self.transform = (
+            translation_to_mid_axis @ rotation @ translation_to_origin
+        )
+
+    def transform_image(self, image: np.ndarray = None):
+        """Transform the image using the transformation matrix.
+
+        Parameters
+        ----------
+        image : np.ndarray
+            The image to transform. If None, the image passed to the
+            constructor is used.
+        """
+        if not hasattr(self, "transform"):
+            self._compute_transform()
+        if image is None:
+            image = self.image
+        self.transformed_image = apply_transform(image, self.transform)
+        return self.transformed_image
+
+    def save_transform(self, dest_path: Path):
+        """Save the midplane alignment transform to a text file."""
+        if not hasattr(self, "transform"):
+            raise ValueError("Please align the image to the midplane first")
+        np.savetxt(dest_path, self.transform)
+
+    def label_halves(self, image: np.ndarray) -> np.ndarray:
+        """Label each half of the image along the symmetry axis with different
+        integer values, to help diagnose issues with the splitting process.
+
+        Parameters
+        ----------
+        image : np.ndarray
+            The image to label. This should be aligned along the symmetry axis.
+
+        Returns
+        -------
+        np.ndarray
+            An array of the same shape as the input image, with each half
+            labelled with a different integer value (2 and 3).
+        """
+        axi = self.symmetry_axis_idx
+        axis_len = image.shape[axi]
+        half_len = axis_len // 2
+        labelled_halves = np.zeros_like(image, dtype=int)
+
+        # Create slicing objects for each half
+        slicer_half1 = [slice(None)] * image.ndim
+        slicer_half1[axi] = slice(0, half_len)
+        slicer_half2 = [slice(None)] * image.ndim
+        slicer_half2[axi] = slice(half_len, axis_len)
+
+        # Apply different integer values to each half
+        labelled_halves[tuple(slicer_half1)] = 2
+        labelled_halves[tuple(slicer_half2)] = 3
+
+        return labelled_halves
diff --git a/brainglobe_template_builder/preproc/masking.py b/brainglobe_template_builder/preproc/masking.py
new file mode 100644
index 0000000..f81c5e3
--- /dev/null
+++ b/brainglobe_template_builder/preproc/masking.py
@@ -0,0 +1,121 @@
+from typing import Literal, Union
+
+import numpy as np
+from skimage import filters, measure, morphology
+
+
+def _extract_largest_object(binary_image):
+    """Keep only the largest object in a binary image.
+
+    Parameters
+    ----------
+    binary_image : np.ndarray
+        A binary image.
+
+    Returns
+    -------
+    np.ndarray
+        A binary image containing only the largest object.
+    """
+    labeled_image = measure.label(binary_image)
+    regions = measure.regionprops(labeled_image)
+    largest_region = max(regions, key=lambda region: region.area)
+    return labeled_image == largest_region.label
+
+
+def _threshold_image(
+    image: np.ndarray,
+    method: Literal["triangle", "otsu", "isodata"] = "triangle",
+) -> Union[np.ndarray, None]:
+    """Threshold an image using the specified method to get a binary mask.
+
+    Parameters
+    ----------
+    image : np.ndarray
+        Image to threshold.
+    method : str
+        Thresholding method to use. One of 'triangle', 'otsu', and 'isodata'
+        (corresponding to methods from the skimage.filters module).
+        Defaults to 'triangle'.
+
+    Returns
+    -------
+    np.ndarray
+        A binary mask.
+    """
+
+    method_to_func = {
+        "triangle": filters.threshold_triangle,
+        "otsu": filters.threshold_otsu,
+        "isodata": filters.threshold_isodata,
+    }
+    if method in method_to_func.keys():
+        thresholded = method_to_func[method](image)
+        return image > thresholded
+    else:
+        raise ValueError(f"Unknown thresholding method {method}")
+
+
+def create_mask(
+    image: np.ndarray,
+    gauss_sigma: float = 3,
+    threshold_method: Literal["triangle", "otsu", "isodata"] = "triangle",
+    closing_size: int = 5,
+    erode_size: int = 0,
+) -> np.ndarray:
+    """Threshold image and create a mask for the largest object.
+
+    The mask is generated by applying a Gaussian filter to the image,
+    thresholding the smoothed image, keeping only the largest object, and
+    applying a binary closing operation to the mask, optionally followed by
+    binary erosion.
+
+    Parameters
+    ----------
+    image : np.ndarray
+        A 3D image to generate the mask from.
+    gauss_sigma : float
+        Standard deviation for Gaussian kernel (in pixels) to smooth image
+        before thresholding. Set to 0 to skip smoothing.
+    threshold_method : str
+        Thresholding method to use. One of 'triangle', 'otsu', and 'isodata'
+        (corresponding to methods from the skimage.filters module).
+        Defaults to 'triangle'.
+    closing_size : int
+        Size of the binary closing footprint (in pixels) to apply to the mask.
+        Set to 0 to skip closing. This is useful to remove small holes in the
+        mask. Larger values will remove larger holes.
+    erode_size : int
+        Size of the binary erosion footprint (in pixels) to apply to the mask.
+        Set to 0 to skip erosion. This is useful to remove small structures
+        at the mask boundary. Larger values will remove larger structures.
+
+    Returns
+    -------
+    mask : np.ndarray
+        A binary mask of the largest object in the image.
+    """
+
+    # Check input
+    if image.ndim != 3:
+        raise ValueError("Image must be 3D")
+
+    if gauss_sigma > 0:
+        data_smoothed = filters.gaussian(image, sigma=gauss_sigma)
+    else:
+        data_smoothed = image
+
+    binary = _threshold_image(data_smoothed, method=threshold_method)
+    mask = _extract_largest_object(binary)
+
+    if closing_size > 0:
+        mask = morphology.binary_closing(
+            mask, footprint=np.ones((closing_size,) * image.ndim)
+        )
+
+    if erode_size > 0:
+        mask = morphology.binary_erosion(
+            mask, footprint=np.ones((erode_size,) * image.ndim)
+        )
+
+    return mask
diff --git a/brainglobe_template_builder/preproc/splitting.py b/brainglobe_template_builder/preproc/splitting.py
new file mode 100644
index 0000000..0873b61
--- /dev/null
+++ b/brainglobe_template_builder/preproc/splitting.py
@@ -0,0 +1,136 @@
+"""Utilities for splitting arrays into two halves along a symmetry axis."""
+
+from pathlib import Path
+
+import numpy as np
+
+from brainglobe_template_builder.io import save_nii
+
+
+def get_right_and_left_slices(array: np.ndarray) -> tuple:
+    """
+    Get the slices for equally splitting an array into right and left halves.
+
+    Parameters
+    ----------
+    array : np.ndarray
+        Array to split. Must be a 3D array in ASR orientation.
+
+    Returns
+    -------
+    A tuple of slice objects for the right and left halves of the array.
+    """
+    slices = [slice(None)] * array.ndim
+    axis = -1  # right-left axis is always the last in ASR orientation
+
+    right_slices = slices.copy()
+    right_slices[axis] = slice(0, array.shape[axis] // 2)
+
+    left_slices = slices.copy()
+    left_slices[axis] = slice(array.shape[axis] // 2, array.shape[axis])
+    return tuple(right_slices), tuple(left_slices)
+
+
+def generate_arrays_4template(
+    subject: str, brain: np.ndarray, mask: np.ndarray, pad: int = 0
+) -> dict[str, np.ndarray]:
+    """Generate all needed arrays for the template building process.
+
+    Parameters
+    ----------
+    subject : str
+        Subject ID to use as a prefix for the keys in the output dictionary,
+        e.g. "sub-01".
+    brain : np.ndarray
+        The aligned brain image to split into hemispheres and symmetrise.
+    mask : np.ndarray
+        The aligned mask to split into hemispheres and symmetrise.
+        Must be a binary mask and have the same shape as the image.
+    pad : int, optional
+        Number of planes to zero-pad arrays with. Default is 0 (no padding).
+        The same number of planes will be added to all arrays, along each axis.
+        This is useful to ensure that registration algorithms can correctly
+        handle the array edges.
+
+    Returns
+    -------
+    A dictionary containing the multiple arrays needed for the template
+    building process. The keys start with the subject ID and end with the
+    following suffixes (e.g. "sub-01_asym-brain"):
+    - asym-brain: the input aligned image (asymmetric brain)
+    - asym-mask: the input aligned mask (asymmetric mask)
+    - right-hemi-brain: the right hemisphere of the image
+    - right-hemi-mask: the right hemisphere of the mask
+    - left-hemi-xflip-brain: the reflection of the left hemisphere
+    - left-hemi-xflip-mask: the reflection of the left hemisphere mask
+    - right-sym-brain: the right hemisphere merged with its reflection
+    - right-sym-mask: the right hemisphere mask merged with its reflection
+    - left-sym-brain: the left hemisphere merged with its reflection
+    - left-sym-mask: the left hemisphere mask merged with its reflection
+    """
+
+    # ensure mask in uint8
+    mask = mask.astype(np.uint8)
+
+    # ensure shapes match
+    assert brain.ndim == 3
+    assert brain.shape == mask.shape
+
+    # Put the input arrays into the dictionary
+    out_dict = {
+        f"{subject}_asym-brain": brain,
+        f"{subject}_asym-mask": mask,
+    }
+
+    right_half, left_half = get_right_and_left_slices(brain)
+
+    for label, arr in zip(("brain", "mask"), (brain, mask)):
+        right_half_arr = arr[right_half]
+        left_half_arr = arr[left_half]
+        right_half_arr_xflip = np.flip(right_half_arr, axis=-1)
+        left_half_arr_xflip = np.flip(left_half_arr, axis=-1)
+        right_sym_arr = np.dstack([right_half_arr, right_half_arr_xflip])
+        left_sym_arr = np.dstack([left_half_arr_xflip, left_half_arr])
+        out_dict.update(
+            {
+                f"{subject}_right-hemi-{label}": right_half_arr,
+                f"{subject}_left-hemi-xflip-{label}": left_half_arr_xflip,
+                f"{subject}_right-sym-{label}": right_sym_arr,
+                f"{subject}_left-sym-{label}": left_sym_arr,
+            }
+        )
+
+    if pad > 0:
+        for key, arr in out_dict.items():
+            out_dict[key] = np.pad(arr, pad_width=pad, mode="constant")
+
+    return out_dict
+
+
+def save_array_dict_to_nii(
+    array_dict: dict[str, np.ndarray],
+    save_dir: Path,
+    vox_sizes: list[float],
+):
+    """Save arrays in a dictionary to NIfTI files in the specified directory.
+
+    Parameters
+    ----------
+    array_dict : dict
+        A dictionary containing numpy arrays to save. The keys should be the
+        filenames (without extension) and the values should be the arrays.
+        If the arrays are masks, the keys should contain the string "mask".
+    save_dir : Path
+        Directory to save the files to.
+    vox_sizes : tuple
+        Tuple of voxel sizes in mm.
+    """
+
+    assert save_dir.is_dir(), f"Directory {save_dir} does not exist."
+
+    for key, data in array_dict.items():
+        file_path = save_dir / f"{key}.nii.gz"
+        if "mask" in key:
+            save_nii(data, vox_sizes, file_path, kind="label")
+        else:
+            save_nii(data, vox_sizes, file_path, kind="image")
diff --git a/brainglobe_template_builder/preproc/transform_utils.py b/brainglobe_template_builder/preproc/transform_utils.py
new file mode 100644
index 0000000..fcd0185
--- /dev/null
+++ b/brainglobe_template_builder/preproc/transform_utils.py
@@ -0,0 +1,135 @@
+import dask.array as da
+import numpy as np
+from scipy.ndimage import affine_transform
+from skimage import transform
+
+
+def get_rotation_from_vectors(vec1: np.ndarray, vec2: np.ndarray):
+    """Find the rotation matrix that aligns vec1 to vec2. Implementation
+    adapted from StackOverflow [1]_.
+
+    Parameters
+    ----------
+    vec1 : np.ndarray
+        The 3D "source" vector
+    vec2 : np.ndarray
+        The 3D "target" vector
+
+    Returns
+    -------
+    A rotation matrix (3x3) that, when applied to vec1, aligns it with vec2.
+
+    References
+    ----------
+    .. [1] https://stackoverflow.com/questions/45142959
+    """
+    a = (vec1 / np.linalg.norm(vec1)).reshape(3)
+    b = (vec2 / np.linalg.norm(vec2)).reshape(3)
+    v = np.cross(a, b)
+    c = np.dot(a, b)
+    s = np.linalg.norm(v)
+    kmat = np.array([[0, -v[2], v[1]], [v[2], 0, -v[0]], [-v[1], v[0], 0]])
+    rotation_matrix = np.eye(3) + kmat + kmat.dot(kmat) * ((1 - c) / (s**2))
+    return rotation_matrix
+
+
+def apply_transform(
+    data: np.ndarray,
+    transform: np.ndarray,
+) -> np.ndarray:
+    """Apply a rigid transformation to an image.
+
+    Parameters
+    ----------
+    data : np.ndarray
+        A 3D image to transform.
+    transform : np.ndarray
+        A 4x4 transformation matrix.
+
+    Returns
+    -------
+    np.ndarray
+        The transformed data.
+
+    Notes
+    -----
+    This function inverts the affine and flips the offset when passing the
+    data to `scipy.ndimage.affine_transform`. This is because the
+    transforms are given in the 'forward' direction, transforming input
+    to output, whereas `scipy.ndimage.affine_transform` does 'backward'
+    resampling, transforming the output space to the input.
+    """
+
+    if data.ndim != 3:
+        raise ValueError("Data must be 3D")
+    if transform.shape != (4, 4):
+        raise ValueError("Transform must be a 4x4 matrix")
+
+    transformed = affine_transform(
+        data,
+        np.linalg.inv(transform[:3, :3]),
+        offset=-transform[:3, 3],
+    )
+    # Preserve original data range and type
+    transformed = np.clip(transformed, data.min(), data.max())
+    return transformed.astype(data.dtype)
+
+
+def downsample_anisotropic_image_stack(
+    stack: da.Array, in_plane_factor: int, axial_factor: int
+) -> np.ndarray:
+    """
+
+    Lazily downsamples a dask array first along axes 1,2 (in-plane) and then
+    along axis 0 (axial), using a local mean of the pixels. The image is
+    zero-padded to allow for the correct dimensions of the averaging
+    neighbourhood, since it uses `skimage.transform.downscale_local_mean`
+    under the hood.
+
+    This setup is typical for certain types of microscopy,
+    where axial resolution is lower than in-plane resolution.
+
+    The input dask array must be chunked by plane. The (smaller) array
+    is returned in memory (numpy) form at the end.
+
+    Parameters:
+    ----------
+    stack : da.Array
+        The input dask array representing the image stack.
+    in_plane_factor : int
+        The in-plane downsampling factor (axes=1,2).
+    axial_factor : int
+        The downsampling factor in axial direction (axis=0).
+    Returns:
+    -------
+    np.ndarray
+        The computed downsampled (numpy) array.
+    Raises:
+    ------
+    AssertionError
+        If the array is not chunked by plane along axis 0.
+    """
+    # check we have expected slice chunks
+    assert np.all(
+        np.array(stack.chunks[0]) == 1
+    ), f"Array not chunked by plane! Chunks on axis 0 are {stack.chunks[0]}"
+
+    # we have xy slices as chunks, so apply downscaling in xy first
+    downsampled_inplane = stack.map_blocks(
+        transform.downscale_local_mean,
+        (1, in_plane_factor, in_plane_factor),
+        dtype=np.float64,
+    )
+
+    # rechunk so we can map_blocks along z
+    downsampled_inplane = downsampled_inplane.rechunk(
+        {0: downsampled_inplane.shape[0], 1: -1, 2: -1}
+    )
+
+    # downsample in z
+    downsampled_axial = downsampled_inplane.map_blocks(
+        transform.downscale_local_mean,
+        (axial_factor, 1, 1),
+        dtype=np.float64,
+    )
+    return downsampled_axial.compute()
diff --git a/data/SWC_BlackCap_use-for-template.csv b/data/SWC_BlackCap_use-for-template.csv
new file mode 100644
index 0000000..2c5fceb
--- /dev/null
+++ b/data/SWC_BlackCap_use-for-template.csv
@@ -0,0 +1,11 @@
+subject_id,channel,color,hemi,comment
+BC58,3,green,right,shadow artifacts in left hemi
+BC73,3,green,both,small damage in midbrain
+BC38,3,green,both,small damage in midbrain and cerebellum
+BC63,4,blue,both,blue channel has the best SNR
+BC15,4,blue,both,may need to mask injection and missing few slices at cerebellum
+BC08,3,green,both,need to mask injection and missing small part of midbrain
+BC71,3,green,both,need to mask injection
+BC46,3,green,left,some shadow artifacts mostly in rigth hemi and missing parts of cerebellum
+BC74,3,green,both,not great but acceptable SNR
+BC41o,3,green,both,excellent quality
diff --git a/data/SWC_brain-list_for-atlases_2024-04-15.csv b/data/SWC_brain-list_for-atlases_2024-04-15.csv
new file mode 100644
index 0000000..6bd00be
--- /dev/null
+++ b/data/SWC_brain-list_for-atlases_2024-04-15.csv
@@ -0,0 +1,64 @@
+Common name,Scientific name,Strain,Specimen ID,Age (Days),Gender (M / F),Injections,Perfusion (whole / drop fix),Imaging approach (serial 2p / lightsheet / other) ,Imaging Res X (microns),Imaging Res Y,Imaging Res Z,Imaging channel to use ,Imaging quality ,Data path (raw),comments
+Zebra finch ,Taeniopygia guttata,bird,SW_ZF7927 red,unknown,unknown,none,whole,serial 2p,2.273,2.273,5,ch3 (green),"some shadow artefacts, use ch3 for atlas",nfs/winstor/margrie/SimonWeiler/RawData/Tracing_Imaging/serial2p/Bird_brains/zebra_finch_atlas,use for atlas
+Zebra finch ,Taeniopygia guttata,bird,SW_ZF8003 red,unknown,unknown,none,whole,serial 2p,2.273,2.273,5,ch3 (green),perfect for atlas: usee chr3,nfs/winstor/margrie/SimonWeiler/RawData/Tracing_Imaging/serial2p/Bird_brains/try_out_brains,use for atlas
+Zebra finch ,Taeniopygia guttata,bird,ZF 3890 yellow,unknown,unknown, AAV-F(129L)-GFP,whole,serial 2p,2.278,2.278,5,ch2 (red),injection trace visible. Right side of brain very dimly imaged (no idea why),nfs/winstor/margrie/SimonWeiler/RawData/Tracing_Imaging/serial2p/Bird_brains/try_out_brains, unusable
+Zebra finch ,Taeniopygia guttata,bird,ZF 4006 yellow,unknown,unknown," left hemisphere AAV-1-GFP, right hemisphere AAV-7-GFP",whole,serial 2p,2.212,2.212,5,ch2 (red),"left hemisphere, supe rbright expression ",nfs/winstor/margrie/SimonWeiler/RawData/Tracing_Imaging/serial2p/Bird_brains/zebra_finch_atlas,"maybe only left hemisphere, because green injection is brigh even in red channel"
+Zebra finch ,Taeniopygia guttata,bird,ZF 65 grey,unknown,unknown, AAV-8-eGFP,whole,serial 2p,2.212,2.212,5,ch3 (green),no expression ,nfs/winstor/margrie/SimonWeiler/RawData/Tracing_Imaging/serial2p/Bird_brains/zebra_finch_atlas,use for atlas
+Zebra finch ,Taeniopygia guttata,bird,ZF3925,unknown,unknown, AAV-9-eGFP,whole,serial 2p,2.212,2.212,5,ch2 (red),expression,nfs/winstor/margrie/SimonWeiler/RawData/Tracing_Imaging/serial2p/Bird_brains/zebra_finch_atlas,use for atlas
+Zebra finch ,Taeniopygia guttata,bird,ZF_28,unknown,unknown,none,whole,serial 2p,2.278,2.278,5,ch3 (green),good,nfs/winstor/margrie/SimonWeiler/RawData/Tracing_Imaging/serial2p/Bird_brains/zebra_finch_atlas,use for atlas
+Zebra finch ,Taeniopygia guttata,bird,ZF_69,unknown,unknown,none,drop fix,serial 2p,2.278,2.278,5,ch3 (green),good,nfs/winstor/margrie/SimonWeiler/RawData/Tracing_Imaging/serial2p/Bird_brains/zebra_finch_atlas,use for atlas
+Zebra finch ,Taeniopygia guttata,bird,ZF1,unknown,unknown,none,drop fix,serial 2p,2.278,2.278,5,ch3 (green),good,nfs/winstor/margrie/SimonWeiler/RawData/Tracing_Imaging/serial2p/Bird_brains/zebra_finch_atlas,use for atlas
+Zebra finch ,Taeniopygia guttata,bird,ZF2,unknown,unknown,none,drop fix,serial 2p,2.278,2.278,5,ch3 (green),good,nfs/winstor/margrie/SimonWeiler/RawData/Tracing_Imaging/serial2p/Bird_brains/zebra_finch_atlas,use for atlas
+Zebra finch ,Taeniopygia guttata,bird,ZF3,unknown,unknown,none,drop fix,serial 2p,2.278,2.278,5,ch3 (green),good,nfs/winstor/margrie/SimonWeiler/RawData/Tracing_Imaging/serial2p/Bird_brains/zebra_finch_atlas,use for atlas
+Zebra finch ,Taeniopygia guttata,bird,ZF5,unknown,unknown,none,drop fix,serial 2p,2.278,2.278,5,ch3 (green),good,nfs/winstor/margrie/SimonWeiler/RawData/Tracing_Imaging/serial2p/Bird_brains/zebra_finch_atlas,use for atlas
+Zebra finch ,Taeniopygia guttata,bird,ZF6,unknown,unknown,none,drop fix,serial 2p,2.278,2.278,5,ch3 (green),good,nfs/winstor/margrie/SimonWeiler/RawData/Tracing_Imaging/serial2p/Bird_brains/zebra_finch_atlas,use for atlas
+Zebra finch ,Taeniopygia guttata,bird,ZF 3883 yellow,unknown,unknown, AAV-PHP.V1-eGFP,whole,serial 2p,2.278,2.278,5,ch2 (red),"good, little expression",nfs/winstor/margrie/SimonWeiler/RawData/Tracing_Imaging/serial2p/Bird_brains/zebra_finch_atlas,use for atlas
+Black cap,Sylvia atricapilla,bird,SW_BC58,unknown,unknown,none ,whole,serial 2p,2.278,2.278,5,ch3 (green),good,nfs/winstor/margrie/SimonWeiler/RawData/Tracing_Imaging/serial2p/Bird_brains/black_cap_atlas,use for atlas
+Black cap,Sylvia atricapilla,bird,SW_BC73,unknown,unknown,none ,whole,serial 2p,2.278,2.278,5,ch3 (green),good,nfs/winstor/margrie/SimonWeiler/RawData/Tracing_Imaging/serial2p/Bird_brains/black_cap_atlas,use for atlas
+Black cap,Sylvia atricapilla,bird,SW_BC38 violet,unknown,unknown,none,whole,serial 2p,2.273,2.273,5,ch3 (green),good,nfs/winstor/margrie/SimonWeiler/RawData/Tracing_Imaging/serial2p/Bird_brains/black_cap_atlas,use for atlas
+Black cap,Sylvia atricapilla,bird,SW_BC40 violet,unknown,unknown,"tracer injection CtB Alexa488, CtB Alexa594",whole,serial 2p,2.273,2.273,5,ch4 (blue),good,nfs/winstor/margrie/SimonWeiler/RawData/Tracing_Imaging/serial2p/Bird_brains/black_cap_atlas,use for atlas
+Black cap,Sylvia atricapilla,bird,SW_BC63 violet,unknown,unknown,"tracer injection CtB Alexa594, CtB Alexa647",whole,serial 2p,2.273,2.273,5,ch4 (blue),good,nfs/winstor/margrie/SimonWeiler/RawData/Tracing_Imaging/serial2p/Bird_brains/black_cap_atlas,use for atlas
+Black cap,Sylvia atricapilla,bird,SW_BC15 grey,unknown,unknown,"tracer injection CtBAlexa488, CtB Alex647",whole,serial 2p,2.273,2.273,5,ch4 (blue),good,nfs/winstor/margrie/SimonWeiler/RawData/Tracing_Imaging/serial2p/Bird_brains/black_cap_atlas,use for atlas
+Black cap,Sylvia atricapilla,bird,SW_BC08 violet,unknown,unknown,scAAV-DJ/9/2 EGFP + beads FluoSpheres,whole,serial 2p,2.278,2.278,5,ch4 (blue),good,nfs/winstor/margrie/SimonWeiler/RawData/Tracing_Imaging/serial2p/Bird_brains/black_cap_atlas,use for atlas
+Black cap,Sylvia atricapilla,bird,SW_BC71 violet,unknown,unknown,scAAV-DJ/9/2 EGFP + beads FluoSpheres,whole,serial 2p,2.278,2.278,5,ch4 (blue),good,nfs/winstor/margrie/SimonWeiler/RawData/Tracing_Imaging/serial2p/Bird_brains/black_cap_atlas,use for atlas
+Black cap,Sylvia atricapilla,bird,SW_BC78 yellow,unknown,unknown,scAAV-DJ/9/2 EGFP beak,whole,serial 2p,2.278,2.278,5,ch3 (green),some shadow artefacts,nfs/winstor/margrie/SimonWeiler/RawData/Tracing_Imaging/serial2p/Bird_brains/black_cap_atlas,use for atlas
+Black cap,Sylvia atricapilla,bird,SW_BC26 yellow,unknown,unknown,AAV mRuby beak,whole,serial 2p,2.278,2.278,5,ch3 (green),good,nfs/winstor/margrie/SimonWeiler/RawData/Tracing_Imaging/serial2p/Bird_brains/black_cap_atlas,use for atlas
+Black cap,Sylvia atricapilla,bird,SW_BC46 yellow,unknown,unknown,scAAV-DJ/9/2 EGFP + AAV mRuby beak,whole,serial 2p,1.667,1.667,5,ch3 (green),missing part of the brain ,nfs/winstor/margrie/SimonWeiler/RawData/Tracing_Imaging/serial2p/Bird_brains/black_cap_atlas,use for atlas
+Black cap,Sylvia atricapilla,bird,SW_BC33 (86 yellow),unknown,unknown,scAAV-DJ/9/2 EGFP,whole,serial 2p,1.667,1.667,5,ch4 (blue),"perfect for atlas (ch4), beautiful",nfs/winstor/margrie/SimonWeiler/RawData/Tracing_Imaging/serial2p/Bird_brains/black_cap_atlas,use for atlas
+Black cap,Sylvia atricapilla,bird,SW_BC74white,unknown,unknown,none,whole,serial 2p,2.278,2.278,5,ch3 (green),USED FOR ANNOTATION IN A SINGLE BRAIN WITH MOURITSEN LAB,nfs/winstor/margrie/SimonWeiler/RawData/Tracing_Imaging/serial2p/Bird_brains/black_cap_atlas,use for atlas
+Black cap,Sylvia atricapilla,bird,SW_BC41v,unknown,unknown,none ,whole,serial 2p,2.278,2.278,5,ch3 (green),part of cerebellum missing,nfs/ceph/margrie/sweiler/RawData/Tracing_Imaging/serial2p/Bird_brains/black_cap_atlas,use for atlas
+Black cap,Sylvia atricapilla,bird,SW_BC41o,unknown,unknown,none ,whole,serial 2p,2.278,2.278,5,ch3 (green),good,nfs/ceph/margrie/sweiler/RawData/Tracing_Imaging/serial2p/Bird_brains/black_cap_atlas,use for atlas
+rat,Rattus norvegicus domestica ,Lister Hooded,CAP01,60-80,F,"IV,  AAV (low or no expression)",,serial 2p,2.278,2.278,5,red,missing tiles & chipped cerebellum and cortex,/nfs/winstor/akrami/anatomy/virus_tests,
+rat,Rattus norvegicus domestica ,Lister Hooded,CAP02,60-80,F,"IV,  AAV (low or no expression)",,serial 2p,2.476,2.476,5.7,red,dark tiles & some labelled cells,/nfs/winstor/akrami/anatomy/virus_tests,
+rat,Rattus norvegicus domestica ,Lister Hooded,CAP03,60-80,F,"IV,  AAV (low or no expression)",,serial 2p,2.278,2.278,5.7,green,looks ok,/nfs/winstor/akrami/anatomy/virus_tests,
+rat,Rattus norvegicus domestica ,Lister Hooded,CAP04,60-80,F,"IV,  AAV (low or no expression)",,serial 2p,2.212,2.212,5.7,green,looks ok,/nfs/winstor/akrami/anatomy/virus_tests,
+rat,Rattus norvegicus domestica ,Lister Hooded,CAP05,60-80,F,"IV,  AAV (low or no expression)",,serial 2p,2.212,2.212,5.6,red,left cortical damage & some labelled cells,/nfs/winstor/akrami/anatomy/virus_tests,
+rat,Rattus norvegicus domestica ,Lister Hooded,CAP08,60-80,F,"IV,  AAV (low or no expression)",,serial 2p,2.212,2.212,6,green,missing tiles & some labelled neurons,/nfs/winstor/akrami/anatomy/virus_tests,
+rat,Rattus norvegicus domestica ,Lister Hooded,CAP09,60-80,F,"IV,  AAV (low or no expression)",,serial 2p,2.212,2.212,6.7,red,labelled cells in the cerebellum,/nfs/winstor/akrami/anatomy/virus_tests,
+rat,Rattus norvegicus domestica ,Lister Hooded,CAP11,60-80,F,"IV,  AAV (low or no expression)",,serial 2p,2.476,2.476,7.1,red,missing tiles,/nfs/winstor/akrami/anatomy/virus_tests,
+rat,Rattus norvegicus domestica ,Lister Hooded,CAP12,60-80,F,"IV,  AAV (low or no expression)",,serial 2p,2.212,2.212,6.3,red,missing tiles & some labelled cells,/nfs/winstor/akrami/anatomy/virus_tests,
+rat,Rattus norvegicus domestica ,Lister Hooded,CAP14,60-80,F,"IV,  AAV (low or no expression)",,serial 2p,2.212,2.212,6.3,green,few missing tiles,/nfs/winstor/akrami/anatomy/virus_tests,
+rat,Rattus norvegicus domestica ,Lister Hooded,CAP18,60-80,F,"IV,  AAV (low or no expression)",,serial 2p,2.476,2.476,7.1,green,missing tiles,/nfs/winstor/akrami/anatomy/virus_tests,
+rat,Rattus norvegicus domestica ,Lister Hooded,CAP22,60-80,F,"IV,  AAV (low or no expression)",,serial 2p,2.476,2.476,7.1,red,missing tiles,/nfs/winstor/akrami/anatomy/virus_tests,
+rat,Rattus norvegicus domestica ,Lister Hooded,CAP24,60-80,F,"IV,  AAV (low or no expression)",,serial 2p,2.278,2.278,7.1,red,chipped cerebellum & missing tiles & few labelled cells,/nfs/winstor/akrami/anatomy/virus_tests,
+rat,Rattus norvegicus domestica ,Lister Hooded,CAP25,60-80,F,"IV,  AAV (low or no expression)",,serial 2p,2.212,2.212,6.3,red,chipped cortex,/nfs/winstor/akrami/anatomy/virus_tests,
+rat,Rattus norvegicus domestica ,Lister Hooded,CAP26,60-80,F,"IV,  AAV (low or no expression)",,serial 2p,2.476,2.476,6.3,red,looks OK,/nfs/winstor/akrami/anatomy/virus_tests,
+rat,Rattus norvegicus domestica ,Lister Hooded,CAP27,60-80,F,"IV,  AAV (low or no expression)",,serial 2p,2.504,2.504,6.3,red,missing tiles & bleedthrough,/nfs/winstor/akrami/anatomy/virus_tests,
+rat,Rattus norvegicus domestica ,Lister Hooded,CAP28,60-80,F,"IV,  AAV (low or no expression)",,serial 2p,2.278,2.278,6.3,blue,looks good might not be bright enough,/nfs/winstor/akrami/anatomy/virus_tests,
+rat,Rattus norvegicus domestica ,Lister Hooded,CAP29,60-80,F,"IV,  AAV (low or no expression)",,serial 2p,2.212,2.212,6.3,blue,looks good might not be bright enough,/nfs/winstor/akrami/anatomy/virus_tests,
+rat,Rattus norvegicus domestica ,Lister Hooded,CAP30,60-80,F,"IV,  AAV (low or no expression)",,serial 2p,2.41,2.41,6.3,blue,looks good might not be bright enough,/nfs/winstor/akrami/anatomy/virus_tests,
+rat,Rattus norvegicus domestica ,Lister Hooded,CAP31,60-80,F,"IV,  AAV (low or no expression)",,serial 2p,2.504,2.504,7.1,red,damaged cerebellum,/ceph/akrami/capsid_testing/imaging/2p,
+rat,Rattus norvegicus domestica ,Lister Hooded,CAP32,60-80,F,"IV,  AAV (low or no expression)",,serial 2p,2.504,2.504,6.3,red,maybe a few missing slices otherwise OK,/ceph/akrami/capsid_testing/imaging/2p,
+rat,Rattus norvegicus domestica ,Lister Hooded,CAP33,60-80,F,"IV,  AAV (low or no expression)",,serial 2p,2.278,2.278,6.3,red,looks good,/ceph/akrami/capsid_testing/imaging/2p,
+rat,Rattus norvegicus domestica ,Lister Hooded,CAP34,60-80,F,"IV,  AAV (low or no expression)",,serial 2p,2.504,2.504,6.3,red,few missing tiles,/ceph/akrami/capsid_testing/imaging/2p,
+rat,Rattus norvegicus domestica ,Lister Hooded,CAP35,60-80,F,"IV,  AAV (low or no expression)",,serial 2p,2.278,2.278,6.3,red,damaged brainstem,/ceph/akrami/capsid_testing/imaging/2p,
+rat,Rattus norvegicus domestica ,Lister Hooded,CAP38,60-80,F,"IV,  AAV (low or no expression)",,serial 2p,1.205,1.205,6.3,green,damaged cortex,/ceph/akrami/capsid_testing/imaging/2p,
+rat,Rattus norvegicus domestica ,Lister Hooded,CAP40,60-80,F,"IV,  AAV (low or no expression)",,serial 2p,2.504,2.504,6.3,red,right cortex is damaged,/ceph/akrami/capsid_testing/imaging/2p,
+rat,Rattus norvegicus domestica ,Lister Hooded,CAP42,60-80,F,"IV,  AAV (low or no expression)",,serial 2p,,,,green,cortex is damaged from recording,/ceph/akrami/capsid_testing/imaging/2p,
+rat,Rattus norvegicus domestica ,Lister Hooded,CAP43,60-80,F,"IV,  AAV (low or no expression)",,serial 2p,2.278,2.278,6.3,red,some cortical damage,/ceph/akrami/capsid_testing/imaging/2p,
+rat,Rattus norvegicus domestica ,Lister Hooded,CAP44,60-80,F,"IV,  AAV (low or no expression)",,serial 2p,2.504,2.504,6.3,red,few dark tiles & some cortical damage,/ceph/akrami/capsid_testing/imaging/2p,
+rat,Rattus norvegicus domestica ,Lister Hooded,CAP45,60-80,F,"IV,  AAV (low or no expression)",,serial 2p,2.278,2.278,6.3,green,cortical damage,/ceph/akrami/capsid_testing/imaging/2p,
+rat,Rattus norvegicus domestica ,Lister Hooded,CAP48,60-80,F,"IV,  AAV (low or no expression)",,serial 2p,2.0,2.0,6.3,red,looks ok,/ceph/akrami/capsid_testing/imaging/2p,
+rat,Rattus norvegicus domestica ,Lister Hooded,CAP49,60-80,F,"IV,  AAV (low or no expression)",,serial 2p,2.0,2.0,6.3,green,little cortical damage,/ceph/akrami/capsid_testing/imaging/2p,
+rat,Rattus norvegicus domestica ,Lister Hooded,CAP50,60-80,F,"IV,  AAV (low or no expression)",,serial 2p,2.0,2.0,6.3,green,damaged cerebellum,/ceph/akrami/capsid_testing/imaging/2p,
+rat,Rattus norvegicus domestica ,Lister Hooded,CAP54,60-80,F,"IV,  AAV (low or no expression)",,serial 2p,2.504,2.504,6.3,green,little cerebellar damage,/ceph/akrami/capsid_testing/imaging/2p,
diff --git a/docs/source/conf.py b/docs/source/conf.py
index 06fe014..3623d07 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -13,7 +13,7 @@
 
 # Used when building API docs, put the dependencies
 # of any class you are documenting here
-autodoc_mock_imports = []
+autodoc_mock_imports: list[str] = []
 
 # Add the module path to sys.path here.
 # If the directory is relative to the documentation root,
@@ -97,7 +97,7 @@
             # Label for this link
             "name": "GitHub",
             # URL where the link will redirect
-            "url": "https://github.com/brainglobe/brainglobe-template-builder",  # required
+            "url": "https://github.com/brainglobe/brainglobe-template-builder",
             # Icon class (if "type": "fontawesome"),
             # or path to local image (if "type": "local")
             "icon": "fa-brands fa-github",
diff --git a/environment.yaml b/environment.yaml
index da5b9d2..a61ef32 100644
--- a/environment.yaml
+++ b/environment.yaml
@@ -1,11 +1,11 @@
 channels:
-  - aramislab
   - conda-forge
-  - defaults
 dependencies:
-  - python=3.10
-  - conda-forge::parallel
-  - conda-forge::napari
-  - conda-forge::pyqt
-  - aramislab::ants
+  - python=3.11
   - pip
+  - parallel
+  - qbatch
+  - scipy
+  - simpleitk
+  - vtk
+  - numpy<2.0
diff --git a/examples/BlackCap_build_template.sh b/examples/BlackCap_build_template.sh
new file mode 100755
index 0000000..f3139fe
--- /dev/null
+++ b/examples/BlackCap_build_template.sh
@@ -0,0 +1,112 @@
+#!/bin/bash
+
+# This script is used to build a BlackCap brain template on a
+# desktop computer (not HPC). It is assumed that the user has already
+# run the BlackCap_source_images.py and BlackCap_build_template.py scripts.
+
+# Start the timer
+start_time=$(date +%s)
+
+# Initialize CLI parameters with default values
+atlas_dir="/ceph/neuroinformatics/neuroinformatics/atlas-forge"
+species=""
+template_name=""
+average_type="mean"
+
+# Function to display help message
+usage() {
+  echo "Usage: $0 [--atlas-dir <path>] --species <string> --template-name <string> [--average-type <string>]"
+  echo ""
+  echo "Options:"
+  echo "  --atlas-dir <path>         Path to the atlas-forge directory (default: /ceph/neuroinformatics/neuroinformatics/atlas-forge)"
+  echo "  --species <string>         Species name (e.g., BlackCap) [REQUIRED]"
+  echo "  --template-name <string>   The name of the appropriate subfolder within templates (e.g., template_asym_res-50um_n-8) [REQUIRED]"
+  echo "  --average-type <string>    The type of average to use (default: mean)."
+  exit 1
+}
+
+# Check for help flag first
+if [[ "$1" == "--help" ]]; then
+  usage
+fi
+
+# Parse command line arguments
+while [[ $# -gt 0 ]]; do
+  case $1 in
+    --atlas-dir)
+      atlas_dir="$2"
+      shift 2
+      ;;
+    --species)
+      species="$2"
+      shift 2
+      ;;
+    --template-name)
+      template_name="$2"
+      shift 2
+      ;;
+    --average-type)
+      average_type="$2"
+      shift 2
+      ;;
+    --help)
+      usage
+      ;;
+    *)
+      echo "Unknown option $1"
+      usage
+      ;;
+  esac
+done
+
+# Check for required parameters
+if [ -z "$species" ] || [ -z "$template_name" ]; then
+  echo "Error: --species and --template-name are required."
+  usage
+fi
+
+echo "atlas-dir: ${atlas_dir}"
+echo "Building the template ${template_name} for the species ${species}..."
+
+# Set the path to the working directory where the template will be built
+species_dir="${atlas_dir}/${species}"
+templates_dir="${species_dir}/templates"
+working_dir="${templates_dir}/${template_name}"
+
+# If average type is trimmed_mean or efficient_trimean, we need python
+average_prog="ANTs"
+if [ "${average_type}" == "trimmed_mean" ] || [ "${average_type}" == "efficient_trimean" ]; then
+  average_prog="python"
+fi
+
+# Verify that the working directory exists before changing directory
+if [ ! -d "${working_dir}" ]; then
+  echo "The working directory does not exist: ${working_dir}"
+  exit 1
+fi
+
+cd "${working_dir}" || exit
+
+echo "Results will be written to ${working_dir}"
+
+# Execute the actual template building
+echo "Starting to build the template..."
+bash modelbuild.sh --output-dir "${working_dir}" \
+    --starting-target first \
+    --stages rigid,similarity,affine,nlin \
+    --masks "${working_dir}/mask_paths.txt" \
+    --average-type "${average_type}" \
+    --average-prog "${average_prog}" \
+    --reuse-affines \
+    --dry-run \
+    "${working_dir}/brain_paths.txt"
+
+echo "Finished building the template!"
+
+# Write execution time (in HH:MM format) to a file
+end_time=$(date +%s)
+execution_time=$((end_time - start_time))
+hours=$((execution_time / 3600))
+minutes=$(( (execution_time % 3600) / 60))
+formatted_time=$(printf "%02d:%02d" $hours $minutes)
+echo "Execution time: $formatted_time" > "${working_dir}/execution_time.txt"
diff --git a/examples/BlackCap_post_build.py b/examples/BlackCap_post_build.py
new file mode 100644
index 0000000..df39907
--- /dev/null
+++ b/examples/BlackCap_post_build.py
@@ -0,0 +1,298 @@
+"""Post-processing after building the BlackCap template.
+
+This script is run after the BlackCap template has been built. It performs
+the following steps:
+1. Select the appropriate template and pad to generate the atlas reference.
+2. Resample the annotated subject image and mask to the atlas resolution.
+3. Generate brain masks for the atlas reference and the annotated subject.
+4. Register the annotated subject image to the atlas reference image.
+5. Transform the annotations to the atlas reference image space.
+"""
+
+# %%
+# Imports
+# -------
+import os
+import subprocess
+from datetime import date
+from pathlib import Path
+
+import ants
+import numpy as np
+from loguru import logger
+
+from brainglobe_template_builder.io import (
+    save_nii,
+)
+from brainglobe_template_builder.preproc.masking import create_mask
+
+# %%
+# Setup
+# -----
+
+# Species-specific atlas-forge directory
+atlas_dir = Path("/ceph/neuroinformatics/neuroinformatics/atlas-forge")
+species_id = "BlackCap"
+species_dir = atlas_dir / species_id
+res_um = 25  # resolution in microns
+
+# Define the path to the final template file to use as the atlas reference
+res_str = f"res-{res_um}um"
+template_name = f"template_sym_{res_str}_n-18"
+# Voxel size in mm (as in nifti)
+res_mm = res_um * 1e-3
+vox_sizes = [res_mm, res_mm, res_mm]
+
+template_build_dir = species_dir / "templates" / template_name
+template_final_dir = template_build_dir / "final" / "average"
+template_file = template_final_dir / "template_sharpen_shapeupdate.nii.gz"
+assert (
+    template_final_dir.is_dir()
+), f"Directory {template_final_dir} does not exist."
+assert template_file.is_file(), f"File {template_file} does not exist."
+
+# Define paths to the annotation file and the correspondind subject's image
+annot_subj = "sub-BC74"
+annot_dir = species_dir / "rawdata" / annot_subj / "annotations"
+annot_img_file = (
+    annot_dir
+    / "1001010_ds_SW_BC74white_220217_120749_10_10_ch03_chan_3_green_raw_oriented.nii"  # noqa
+)
+assert annot_img_file.is_file(), f"File {annot_img_file} does not exist."
+annot_labels_file = (
+    annot_dir / "BC74white_100um_annoations_IM_200522_16_Feb_2024.nii"
+)
+assert annot_labels_file.is_file(), f"File {annot_labels_file} does not exist."
+
+# Set up logging
+today = date.today()
+current_script_name = os.path.basename(__file__).replace(".py", "")
+logger.add(species_dir / "logs" / f"{today}_{current_script_name}.log")
+
+# Create the output directory
+output_dir = species_dir / "templates" / template_name / "for_atlas"
+output_dir.mkdir(parents=True, exist_ok=True)
+logger.info(f"Post-processing {species_id} {template_name}.")
+logger.info(f"Outputs will be saved to {output_dir}.")
+
+
+# %%
+# Create atlas reference image
+# ----------------------------
+# The atlas reference image is the final template with some zero-padding
+# (to avoid masks/labels being cut off at the edges).
+
+n_pad = 6  # Number of planes to pad on each side
+
+# Load the template file and add n_pad planes of zero-padding to each side
+template = ants.image_read(template_file.as_posix())
+template_padded = np.pad(template.numpy(), pad_width=n_pad, mode="constant")
+logger.debug(f"Template zero-padded to shape: {template_padded.shape}.")
+
+# Save the padded template to the output directory
+# Henceforth we will refer to it as the atlas reference image.
+atlas_reference_file = output_dir / f"reference_{res_str}_image.nii.gz"
+save_nii(template_padded, vox_sizes, atlas_reference_file, kind="image")
+logger.debug(f"Saved atlas reference image as {atlas_reference_file.name}.")
+
+# Plot the atlas reference image to check
+atlas_reference = ants.image_read(atlas_reference_file.as_posix())
+ants.plot(
+    atlas_reference,
+    axis=1,
+    title="Atlas reference image",
+    filename=output_dir / f"reference_{res_str}_image.png",
+)
+
+# %%
+# Resample annotated subject image and mask to the atlas reference resolution
+# ---------------------------------------------------------------------------
+
+annot_img = ants.image_read(annot_img_file.as_posix())
+annot_labels = ants.image_read(annot_labels_file.as_posix())
+
+annot_img_res = ants.resample_image(
+    image=annot_img,
+    resample_params=vox_sizes,  # Target voxel sizes
+    use_voxels=False,
+    interp_type=4,  # B-spline interpolation
+)
+annot_labels_res = ants.resample_image(
+    image=annot_labels,
+    resample_params=vox_sizes,  # Target voxel sizes
+    use_voxels=False,
+    interp_type=1,  # Nearest neighbor interpolation
+)
+
+annot_img_res_file = output_dir / f"{annot_subj}_{res_str}_image.nii.gz"
+ants.image_write(annot_img_res, annot_img_res_file.as_posix())
+logger.debug(f"Resampled annotated {annot_subj} image to {res_str} ")
+logger.debug(f"and saved as {annot_img_res_file.name}.")
+
+annot_labels_res_file = output_dir / f"{annot_subj}_{res_str}_labels.nii.gz"
+ants.image_write(annot_labels_res, annot_labels_res_file.as_posix())
+logger.debug(f"Resampled annotated {annot_subj} labels to {res_str} ")
+logger.debug(f"and saved as {annot_labels_res_file.name}.")
+
+
+# %%
+# Generate brain masks
+# --------------------
+# Two masks in the atlas reference image space:
+# - A loose mask for registration, e.g. "reference_res-50um_mask-4reg"
+# - A tight mask tos serve as root label, e.g. "reference_res-50um_mask-brain"
+# One mask in the resampled annotated subject image space:
+# - A loose mask for registration, e.g. "sub-BC74_res-50um_mask-4reg"
+
+# Loose mask in atlas reference space (for registration)
+reg_mask_arr = create_mask(
+    atlas_reference.numpy(),
+    gauss_sigma=3,
+    threshold_method="triangle",
+    closing_size=5,
+    erode_size=0,
+)
+reg_mask = atlas_reference.new_image_like(reg_mask_arr.astype(np.uint8))
+reg_mask_file = output_dir / f"reference_{res_str}_mask-4reg.nii.gz"
+ants.image_write(reg_mask, reg_mask_file.as_posix())
+
+# Tight mask in atlas reference space (for root label)
+brain_mask_arr = create_mask(
+    atlas_reference.numpy(),
+    gauss_sigma=2,
+    threshold_method="triangle",
+    closing_size=4,
+    erode_size=4,
+)
+brain_mask = atlas_reference.new_image_like(brain_mask_arr.astype(np.uint8))
+brain_mask_file = output_dir / f"reference_{res_str}_mask-brain.nii.gz"
+ants.image_write(brain_mask, brain_mask_file.as_posix())
+
+# Loose mask in annotated subject image space (for registration)
+subj_reg_mask_arr = create_mask(
+    annot_img_res.numpy(),
+    gauss_sigma=3,
+    threshold_method="triangle",
+    closing_size=5,
+    erode_size=0,
+)
+subj_reg_mask = annot_img_res.new_image_like(
+    subj_reg_mask_arr.astype(np.uint8)
+)
+subj_reg_mask_file = output_dir / f"{annot_subj}_{res_str}_mask-4reg.nii.gz"
+ants.image_write(subj_reg_mask, subj_reg_mask_file.as_posix())
+
+logger.debug("Generated brain masks.")
+
+# Plot the masks over the reference image to check
+for bg_name, bg_img, mask_name, mask_img in zip(
+    ("reference", "reference", annot_subj),
+    (atlas_reference, atlas_reference, annot_img_res),
+    ("4reg", "brain", "4reg"),
+    (reg_mask, brain_mask, subj_reg_mask),
+):
+    mask_plot_path = (
+        output_dir / f"{bg_name}_{res_str}_mask-{mask_name}_overlay.png"
+    )
+    ants.plot(
+        bg_img,
+        mask_img,
+        overlay_alpha=0.5,
+        axis=1,
+        title=f"mask-{mask_name} over {bg_name} image",
+        filename=mask_plot_path.as_posix(),
+    )
+logger.debug("Plotted overlays to visually check masks.")
+
+
+# %%
+# Register annotated subject image to the atlas reference image
+# --------------------------------------------------------------
+
+output_prefix = output_dir / f"{annot_subj}_to_reference_{res_str}_"
+
+# Call the antsRegistration_affine_SyN.sh script with the provided parameters
+
+cmd = "antsRegistration_affine_SyN.sh "
+cmd += f"--moving-mask {subj_reg_mask_file.as_posix()} "
+cmd += f"--fixed-mask {reg_mask_file.as_posix()} "
+cmd += "--clobber "  # Overwrite existing files
+cmd += f"{annot_img_res_file.as_posix()} "  # moving image
+cmd += f"{atlas_reference_file.as_posix()} "  # fixed image
+cmd += f"{output_prefix.as_posix()}"  # output prefix
+
+logger.debug(f"Running the following ANTs registration script: {cmd}.")
+subprocess.run(cmd, shell=True, check=True)
+logger.debug(
+    f"Registered {annot_subj} {res_str} image to atlas reference image."
+)
+
+# %%
+# Transform the annotations to the atlas reference image space
+# -------------------------------------------------------------
+
+# Order of transforms is as in linear algebra (last applied first)
+transforms = [
+    output_prefix.as_posix() + "1Warp.nii.gz",
+    output_prefix.as_posix() + "0GenericAffine.mat",
+]
+
+# Apply the transforms to the resampled annotated subject image
+aligned_img = ants.apply_transforms(
+    fixed=atlas_reference,
+    moving=annot_img_res,
+    transformlist=transforms,
+    interpolator="bSpline",
+)
+# Save the transformed image to the output directory
+aligned_img_file = (
+    output_dir / f"{annot_subj}_{res_str}_image_aligned-to-reference.nii.gz"
+)
+ants.image_write(aligned_img, aligned_img_file.as_posix())
+logger.debug(
+    f"Transformed {annot_subj} {res_str} image to atlas reference image space."
+)
+logger.debug(f"Saved transformed image as {aligned_img_file.name}.")
+
+
+aligned_labels = ants.apply_transforms(
+    fixed=atlas_reference,
+    moving=annot_labels_res,
+    transformlist=transforms,
+    interpolator="genericLabel",
+)
+# Save the transformed annotations to the output directory
+aligned_labels_file = (
+    output_dir / f"{annot_subj}_{res_str}_labels_aligned-to-reference.nii.gz"
+)
+ants.image_write(aligned_labels, aligned_labels_file.as_posix())
+logger.debug(
+    f"Transformed {annot_subj} {res_str} annotations to atlas reference image."
+)
+logger.debug(f"Saved transformed annotations as {aligned_labels_file.name}.")
+
+
+# Plot the aligned image over the reference image to check registration
+ants.plot(
+    atlas_reference,
+    aligned_img,
+    overlay_alpha=0.5,
+    overlay_cmap="plasma",
+    vminol=-250,
+    vmaxol=7000,
+    axis=1,
+    title=f"{annot_subj} {res_str} image aligned to the atlas reference",
+    filename=output_prefix.as_posix() + "image-overlay.png",
+)
+
+# Plot the aligned annotation over the reference image to check registration
+ants.plot(
+    atlas_reference,
+    aligned_labels,
+    overlay_alpha=0.5,
+    overlay_cmap="tab20",
+    axis=1,
+    title=f"{annot_subj} {res_str} annotations aligned to the atlas reference",
+    filename=output_prefix.as_posix() + "labels-overlay.png",
+)
+logger.debug("Plotted overlays to visually check registration.")
diff --git a/examples/BlackCap_prep_highres.py b/examples/BlackCap_prep_highres.py
new file mode 100644
index 0000000..546dd6b
--- /dev/null
+++ b/examples/BlackCap_prep_highres.py
@@ -0,0 +1,334 @@
+"""
+Prepare high-resolution BlackCap images for template construction
+================================================================
+The following operations are performed on the lowest-resolution images to
+prepare them for template construction:
+- Upsample aligned space into high resolution
+- Import each image as tiff, re-orient to ASR and save as nifti
+- Perform N4 Bias field correction using ANTs
+- Upsample from brain mask from lowres to highres
+- Transform the image and mask into the aligned space using existing transforms
+- Split the image and mask into hemispheres and reflect each hemisphere
+- Generate symmetric brains using either the left or right hemisphere
+- Save all resulting images as nifti files to be used for template construction
+"""
+
+# %%
+# Imports
+# -------
+import os
+from datetime import date
+from pathlib import Path
+
+import ants
+import numpy as np
+import pandas as pd
+from brainglobe_space import AnatomicalSpace
+from loguru import logger
+from tqdm import tqdm
+
+from brainglobe_template_builder.io import (
+    file_path_with_suffix,
+    load_tiff,
+    save_nii,
+)
+from brainglobe_template_builder.preproc.splitting import (
+    generate_arrays_4template,
+    save_array_dict_to_nii,
+)
+
+# %%
+# Setup
+# -----
+# Define some global variables, including paths to input and output directories
+# and set up logging.
+
+# Define voxel size(in microns) of the lowest resolution image
+lowres = 50
+highres = 25
+# String to identify the resolution in filenames
+lowres_str = f"res-{lowres}um"
+highres_str = f"res-{highres}um"
+# Define voxel sizes in mm (for Nifti saving)
+lowre_vox_sizes = [lowres * 1e-3] * 3  # in mm
+highres_vox_sizes = [highres * 1e-3] * 3  # in mm
+
+# Prepare directory structure
+atlas_dir = Path("/media/ceph-niu/neuroinformatics/atlas-forge")
+species_id = "BlackCap"
+species_dir = atlas_dir / species_id
+raw_dir = species_dir / "rawdata"
+deriv_dir = species_dir / "derivatives"
+assert deriv_dir.exists(), f"Could not find derivatives directory {deriv_dir}."
+
+# Set up logging
+today = date.today()
+current_script_name = os.path.basename(__file__).replace(".py", "")
+logger.add(species_dir / "logs" / f"{today}_{current_script_name}.log")
+
+# %%
+# Define registration target
+# --------------------------
+# We need three images, all in ASR orientation and nifti format:
+# 1. The reference image aligned to the mid-sagittal plane
+# 2. The brain mask of the reference image
+# 3. The mask of the halves of the reference image (left and right)
+
+# Here we used the previously generated template as a target
+# We have manually aligned the template to the mid-sagittal plane with napari
+# and prepared the masks accordingly.
+
+target_dir = species_dir / "templates" / "template_asym_res-50um_n-3"
+target_image_lowres = ants.image_read(
+    (target_dir / "template_orig-asr_aligned.nii.gz").as_posix()
+)
+target_mask_lowres = ants.image_read(
+    (target_dir / "template_orig-asr_label-brain_aligned.nii.gz").as_posix()
+)
+target_halves_mask_lowres = ants.image_read(
+    (target_dir / "template_orig-asr_label-halves_aligned.nii.gz").as_posix()
+)
+
+# Upsample the targets to high resolution
+lowres_shape = target_image_lowres.shape
+factor = lowres // highres
+highres_shape = tuple([int(dim * factor) for dim in lowres_shape])
+# For the image use B-spline interpolation (interp_type=4)
+target_image = ants.resample_image(
+    target_image_lowres, highres_shape, use_voxels=True, interp_type=4
+)
+# For masks use Nearest Neighbor interpolation (interp_type=1)
+target_mask = ants.resample_image(
+    target_mask_lowres, highres_shape, use_voxels=True, interp_type=1
+)
+target_halves_mask = ants.resample_image(
+    target_halves_mask_lowres, highres_shape, use_voxels=True, interp_type=1
+)
+
+# %%
+# Load a dataframe with image paths to use for the template
+# ---------------------------------------------------------
+
+source_csv_dir = species_dir / "templates" / "use_for_template.csv"
+df = pd.read_csv(source_csv_dir)
+n_subjects = len(df)
+
+# take the hemi column and split it into two boolean columns
+# named use_left and use_right (for easier filtering)
+df["use_left"] = df["hemi"].isin(["left", "both"])
+df["use_right"] = df["hemi"].isin(["right", "both"])
+
+df.head(n_subjects)
+
+# %%
+# Run the pipeline for each subject
+# ---------------------------------
+
+# Create a dictionary to store the paths to the use4template directories
+# per subject. These will contain all necessary images for template building.
+use4template_dirs = {}
+
+for idx, row in tqdm(df.iterrows(), total=n_subjects):
+    # Figure out input-output paths
+    row = df.iloc[idx]
+    subject = "sub-" + row["subject_id"]
+    channel_str = "channel-" + row["color"]
+    file_prefix = f"{subject}_{highres_str}_{channel_str}"
+    file_prefix_lowres = f"{subject}_{lowres_str}_{channel_str}"
+    deriv_subj_dir = deriv_dir / subject
+    raw_subj_dir = raw_dir / subject
+    tiff_path = raw_subj_dir / f"{file_prefix}.tif"
+
+    logger.info(f"Starting to process {file_prefix}...")
+    logger.info(f"Will save outputs to {deriv_subj_dir}/")
+
+    # Load the image
+    image = load_tiff(tiff_path)
+    logger.debug(f"Loaded image {tiff_path.name} with shape: {image.shape}.")
+
+    # Reorient the image to ASR
+    source_origin = ["P", "S", "R"]
+    target_origin = ["A", "S", "R"]
+    source_space = AnatomicalSpace(source_origin, shape=image.shape)
+    image_asr = source_space.map_stack_to(target_origin, image)
+    logger.debug(f"Reoriented image from {source_origin} to {target_origin}.")
+    logger.debug(f"Reoriented image shape: {image_asr.shape}.")
+
+    # Save the reoriented image as nifti
+    nii_path = deriv_subj_dir / f"{file_prefix}_orig-asr.nii.gz"
+    save_nii(image_asr, highres_vox_sizes, nii_path, kind="image")
+    logger.debug(f"Saved reoriented image as {nii_path.name}.")
+
+    # Bias field correction (to homogenise intensities)
+    image_ants = ants.image_read(nii_path.as_posix())
+    image_n4 = ants.n4_bias_field_correction(image_ants)
+    image_n4_path = file_path_with_suffix(nii_path, "_N4")
+    ants.image_write(image_n4, image_n4_path.as_posix())
+    logger.debug(
+        f"Created N4 bias field corrected image as {image_n4_path.name}."
+    )
+
+    # Upsample the brain mask from lowres to highres
+    mask_lowres_path = (
+        deriv_subj_dir / f"{file_prefix_lowres}_orig-asr_N4_mask.nii.gz"
+    )
+    mask_lowres = ants.image_read(mask_lowres_path.as_posix())
+    highres_shape = image_n4.shape
+    mask = ants.resample_image(
+        mask_lowres, highres_shape, use_voxels=True, interp_type=1
+    )
+    mask_path = file_path_with_suffix(nii_path, "_N4_mask")
+    ants.image_write(mask, mask_path.as_posix())
+    logger.debug(
+        f"Upsampled brain mask from {mask_lowres_path.name} "
+        f"to {mask_path.name}."
+    )
+
+    # Plot the mask over the image to check
+    mask_plot_path = (
+        deriv_subj_dir / f"{file_prefix}_orig-asr_N4_mask-overlay.png"
+    )
+    ants.plot(
+        image_n4,
+        mask,
+        overlay_alpha=0.5,
+        axis=1,
+        title="Brain mask over image",
+        filename=mask_plot_path.as_posix(),
+    )
+    logger.debug("Plotted overlay to visually check mask.")
+
+    # Transform the reoriented image to an already aligned target
+    # find a file in deriv_subj_dir that ends with _0GenericAffine.mat
+    # this is the transform we want to use
+    transforms = list(deriv_subj_dir.glob("*_0GenericAffine.mat"))
+    rigid_transform = [transforms[0].as_posix()]
+    assert len(rigid_transform) == 1, "Expected one affine transform."
+
+    # Transform the image to the aligned space
+    aligned_image = ants.apply_transforms(
+        fixed=target_image,
+        moving=image_n4,
+        transformlist=rigid_transform,
+        interpolator="bSpline",
+    )
+    logger.debug("Transformed image to aligned space.")
+
+    # Transform the brain mask to the aligned image space
+    aligned_mask = ants.apply_transforms(
+        fixed=target_image,
+        moving=mask,
+        transformlist=rigid_transform,
+        interpolator="nearestNeighbor",
+    )
+    logger.debug("Transformed brain mask to aligned space.")
+
+    output_prefix = file_path_with_suffix(nii_path, "_N4_aligned_", new_ext="")
+    # Plot the aligned image over the target to check registration
+    ants.plot(
+        target_image,
+        aligned_image,
+        overlay_alpha=0.5,
+        overlay_cmap="plasma",
+        axis=1,
+        title="Aligned image over target (rigid registration)",
+        filename=output_prefix.as_posix() + "target-overlay.png",
+    )
+    # Plot the halves mask over the aligned image to check the split
+    ants.plot(
+        aligned_image,
+        target_halves_mask,
+        overlay_alpha=0.5,
+        axis=1,
+        title="Aligned image split into right and left halves",
+        filename=output_prefix.as_posix() + "halves-overlay.png",
+    )
+    logger.debug("Plotted overlays to visually check alignment.")
+
+    # Generate arrays for template construction and save as niftis
+    use4template_dir = Path(output_prefix.as_posix() + "padded_use4template")
+    # if it exists, delete existing files in it
+    if use4template_dir.exists():
+        logger.warning(f"Removing existing files in {use4template_dir}.")
+        for file in use4template_dir.glob("*"):
+            file.unlink()
+    use4template_dir.mkdir(exist_ok=True)
+
+    array_dict = generate_arrays_4template(
+        subject, aligned_image.numpy(), aligned_mask.numpy(), pad=2
+    )
+    save_array_dict_to_nii(array_dict, use4template_dir, highres_vox_sizes)
+    use4template_dirs[subject] = use4template_dir
+    logger.info(
+        f"Saved images for template construction in {use4template_dir}."
+    )
+    logger.info(f"Finished processing {file_prefix}.")
+
+
+# %%
+# Generate lists of file paths for template construction
+# -----------------------------------------------------
+# Use the paths to the use4template directories to generate lists of file paths
+# for the template construction pipeline. Three kinds of template will be
+# generated, and each needs the corresponging brain image and mask files:
+# 1. All asym* files for subjects where hemi=both. These will be used to
+#    generate an asymmetric brain template.
+# 2. All right-sym* files for subjects where use_right is True, and
+#    all left-sym* files for subjects where use_left is True.
+#    These will be used to generate a symmetric brain template.
+# 3. All right-hemi* files for subjects where use_right is True,
+#    and all left-hemi-xflip* files for subjects where use_left is True.
+#    These will be used to generate a symmetric hemisphere template.
+
+filepath_lists: dict[str, list] = {
+    "asym-brain": [],
+    "asym-mask": [],
+    "sym-brain": [],
+    "sym-mask": [],
+    "hemi-brain": [],
+    "hemi-mask": [],
+}
+
+for _, row in df.iterrows():
+    subject = "sub-" + row["subject_id"]
+    use4template_dir = use4template_dirs[subject]
+
+    if row["hemi"] == "both":
+        # Add paths for the asymmetric brain template
+        for label in ["brain", "mask"]:
+            filepath_lists[f"asym-{label}"].append(
+                use4template_dir / f"{subject}_asym-{label}.nii.gz"
+            )
+
+    if row["use_right"]:
+        for label in ["brain", "mask"]:
+            # Add paths for the symmetric brain template
+            filepath_lists[f"sym-{label}"].append(
+                use4template_dir / f"{subject}_right-sym-{label}.nii.gz"
+            )
+            # Add paths for the hemispheric template
+            filepath_lists[f"hemi-{label}"].append(
+                use4template_dir / f"{subject}_right-hemi-{label}.nii.gz"
+            )
+
+    if row["use_left"]:
+        for label in ["brain", "mask"]:
+            # Add paths for the symmetric brain template
+            filepath_lists[f"sym-{label}"].append(
+                use4template_dir / f"{subject}_left-sym-{label}.nii.gz"
+            )
+            # Add paths for the hemispheric template
+            filepath_lists[f"hemi-{label}"].append(
+                use4template_dir / f"{subject}_left-hemi-xflip-{label}.nii.gz"
+            )
+
+# %%
+# Save the file paths to text files, each in a separate directory
+
+for key, paths in filepath_lists.items():
+    kind, label = key.split("-")  # e.g. "asym" and "brain"
+    n_images = len(paths)
+    template_name = f"template_{kind}_{highres_str}_n-{n_images}"
+    template_dir = species_dir / "templates" / template_name
+    template_dir.mkdir(exist_ok=True)
+    np.savetxt(template_dir / f"{label}_paths.txt", paths, fmt="%s")
diff --git a/examples/BlackCap_prep_lowres.py b/examples/BlackCap_prep_lowres.py
new file mode 100644
index 0000000..d1d9473
--- /dev/null
+++ b/examples/BlackCap_prep_lowres.py
@@ -0,0 +1,321 @@
+"""
+Prepare low-resolution BlackCap images for template construction
+================================================================
+The following operations are performed on the lowest-resolution images to
+prepare them for template construction:
+- Import each image as tiff, re-orient to ASR and save as nifti
+- Perform N4 Bias field correction using ANTs
+- Generate brain mask based on N4-corrected image
+- Rigid-register the re-oriented image to an already aligned target (with ANTs)
+- Transform the image and mask into the aligned space
+- Split the image and mask into hemispheres and reflect each hemisphere
+- Generate symmetric brains using either the left or right hemisphere
+- Save all resulting images as nifti files to be used for template construction
+"""
+
+# %%
+# Imports
+# -------
+import os
+from datetime import date
+from pathlib import Path
+
+import ants
+import numpy as np
+import pandas as pd
+from brainglobe_space import AnatomicalSpace
+from loguru import logger
+from tqdm import tqdm
+
+from brainglobe_template_builder.io import (
+    file_path_with_suffix,
+    load_tiff,
+    save_nii,
+)
+from brainglobe_template_builder.preproc.masking import create_mask
+from brainglobe_template_builder.preproc.splitting import (
+    generate_arrays_4template,
+    save_array_dict_to_nii,
+)
+
+# %%
+# Setup
+# -----
+# Define some global variables, including paths to input and output directories
+# and set up logging.
+
+# Define voxel size(in microns) of the lowest resolution image
+lowres = 50
+# String to identify the resolution in filenames
+res_str = f"res-{lowres}um"
+# Define voxel sizes in mm (for Nifti saving)
+vox_sizes = [lowres * 1e-3] * 3  # in mm
+
+# Prepare directory structure
+atlas_dir = Path("/ceph/neuroinformatics/neuroinformatics/atlas-forge")
+species_id = "BlackCap"
+species_dir = atlas_dir / species_id
+deriv_dir = species_dir / "derivatives"
+assert deriv_dir.exists(), f"Could not find derivatives directory {deriv_dir}."
+
+# Set up logging
+today = date.today()
+current_script_name = os.path.basename(__file__).replace(".py", "")
+logger.add(species_dir / "logs" / f"{today}_{current_script_name}.log")
+
+# %%
+# Define registration target
+# --------------------------
+# We need three images, all in ASR orientation and nifti format:
+# 1. The reference image aligned to the mid-sagittal plane
+# 2. The brain mask of the reference image
+# 3. The mask of the halves of the reference image (left and right)
+
+# Here we used the previously generated template as a target
+# We have manually aligned the template to the mid-sagittal plane with napari
+# and prepared the masks accordingly.
+
+target_dir = species_dir / "templates" / "template_asym_res-50um_n-3"
+target_image = ants.image_read(
+    (target_dir / "template_orig-asr_aligned.nii.gz").as_posix()
+)
+target_mask = ants.image_read(
+    (target_dir / "template_orig-asr_label-brain_aligned.nii.gz").as_posix()
+)
+target_halves_mask = ants.image_read(
+    (target_dir / "template_orig-asr_label-halves_aligned.nii.gz").as_posix()
+)
+
+# %%
+# Load a dataframe with image paths to use for the template
+# ---------------------------------------------------------
+
+source_csv_dir = species_dir / "templates" / "use_for_template.csv"
+df = pd.read_csv(source_csv_dir)
+n_subjects = len(df)
+
+# take the hemi column and split it into two boolean columns
+# named use_left and use_right (for easier filtering)
+df["use_left"] = df["hemi"].isin(["left", "both"])
+df["use_right"] = df["hemi"].isin(["right", "both"])
+
+df.head(n_subjects)
+
+# %%
+# Run the pipeline for each subject
+# ---------------------------------
+
+# Create a dictionary to store the paths to the use4template directories
+# per subject. These will contain all necessary images for template building.
+use4template_dirs = {}
+
+for idx, row in tqdm(df.iterrows(), total=n_subjects):
+    # Figure out input-output paths
+    row = df.iloc[idx]
+    subject = "sub-" + row["subject_id"]
+    channel_str = "channel-" + row["color"]
+    file_prefix = f"{subject}_{res_str}_{channel_str}"
+    deriv_subj_dir = deriv_dir / subject
+    tiff_path = deriv_subj_dir / f"{file_prefix}.tif"
+
+    logger.info(f"Starting to process {file_prefix}...")
+    logger.info(f"Will save outputs to {deriv_subj_dir}/")
+
+    # Load the image
+    image = load_tiff(tiff_path)
+    logger.debug(f"Loaded image {tiff_path.name} with shape: {image.shape}.")
+
+    # Reorient the image to ASR
+    source_origin = ["P", "S", "R"]
+    target_origin = ["A", "S", "R"]
+    source_space = AnatomicalSpace(source_origin, shape=image.shape)
+    image_asr = source_space.map_stack_to(target_origin, image)
+    logger.debug(f"Reoriented image from {source_origin} to {target_origin}.")
+    logger.debug(f"Reoriented image shape: {image_asr.shape}.")
+
+    # Save the reoriented image as nifti
+    nii_path = file_path_with_suffix(tiff_path, "_orig-asr", new_ext=".nii.gz")
+    save_nii(image_asr, vox_sizes, nii_path, kind="image")
+    logger.debug(f"Saved reoriented image as {nii_path.name}.")
+
+    # Bias field correction (to homogenise intensities)
+    image_ants = ants.image_read(nii_path.as_posix())
+    image_n4 = ants.n4_bias_field_correction(image_ants)
+    image_n4_path = file_path_with_suffix(nii_path, "_N4")
+    ants.image_write(image_n4, image_n4_path.as_posix())
+    logger.debug(
+        f"Created N4 bias field corrected image as {image_n4_path.name}."
+    )
+
+    # Generate a brain mask based on the N4-corrected image
+    mask_data = create_mask(
+        image_n4.numpy(),
+        gauss_sigma=3,
+        threshold_method="triangle",
+        closing_size=5,
+    )
+    mask_path = file_path_with_suffix(nii_path, "_N4_mask")
+    mask = image_n4.new_image_like(mask_data.astype(np.uint8))
+    ants.image_write(mask, mask_path.as_posix())
+    logger.debug(
+        f"Generated brain mask with shape: {mask.shape} "
+        f"and saved as {mask_path.name}."
+    )
+
+    # Plot the mask over the image to check
+    mask_plot_path = (
+        deriv_subj_dir / f"{file_prefix}_orig-asr_N4_mask-overlay.png"
+    )
+    ants.plot(
+        image_n4,
+        mask,
+        overlay_alpha=0.5,
+        axis=1,
+        title="Brain mask over image",
+        filename=mask_plot_path.as_posix(),
+    )
+    logger.debug("Plotted overlay to visually check mask.")
+
+    # Rigid-register the reoriented image to an already aligned target
+    output_prefix = file_path_with_suffix(nii_path, "_N4_aligned_", new_ext="")
+
+    xfm = ants.registration(
+        fixed=target_image,
+        moving=image_n4,
+        mask=target_mask,
+        moving_mask=mask,
+        type_of_transform="Rigid",
+        initial_transform=None,
+        verbose=False,
+        outprefix=output_prefix.as_posix(),
+    )
+    logger.debug(
+        "Aligned the reoriented image to the target via rigid registration."
+    )
+    aligned_image = xfm["warpedmovout"]
+
+    # Ensure we only grab the correct transform (forward rigid)
+    rigid_transform = [
+        file_path
+        for file_path in xfm["fwdtransforms"]
+        if file_path.endswith("0GenericAffine.mat")
+    ]
+    assert len(rigid_transform) == 1, "Expected one affine transform."
+
+    # Transform the brain mask to the aligned image space
+    aligned_mask = ants.apply_transforms(
+        fixed=target_image,
+        moving=mask,
+        transformlist=rigid_transform,
+        interpolator="nearestNeighbor",
+    )
+    logger.debug("Transformed brain mask to aligned space.")
+
+    # Plot the aligned image over the target to check registration
+    ants.plot(
+        target_image,
+        aligned_image,
+        overlay_alpha=0.5,
+        overlay_cmap="plasma",
+        axis=1,
+        title="Aligned image over target (rigid registration)",
+        filename=output_prefix.as_posix() + "target-overlay.png",
+    )
+    # Plot the halves mask over the aligned image to check the split
+    ants.plot(
+        aligned_image,
+        target_halves_mask,
+        overlay_alpha=0.5,
+        axis=1,
+        title="Aligned image split into right and left halves",
+        filename=output_prefix.as_posix() + "halves-overlay.png",
+    )
+    logger.debug("Plotted overlays to visually check alignment.")
+
+    # Generate arrays for template construction and save as niftis
+    use4template_dir = Path(output_prefix.as_posix() + "padded_use4template")
+    # if it exists, delete existing files in it
+    if use4template_dir.exists():
+        logger.warning(f"Removing existing files in {use4template_dir}.")
+        for file in use4template_dir.glob("*"):
+            file.unlink()
+    use4template_dir.mkdir(exist_ok=True)
+
+    array_dict = generate_arrays_4template(
+        subject, aligned_image.numpy(), aligned_mask.numpy(), pad=2
+    )
+    save_array_dict_to_nii(array_dict, use4template_dir, vox_sizes)
+    use4template_dirs[subject] = use4template_dir
+    logger.info(
+        f"Saved images for template construction in {use4template_dir}."
+    )
+    logger.info(f"Finished processing {file_prefix}.")
+
+
+# %%
+# Generate lists of file paths for template construction
+# -----------------------------------------------------
+# Use the paths to the use4template directories to generate lists of file paths
+# for the template construction pipeline. Three kinds of template will be
+# generated, and each needs the corresponging brain image and mask files:
+# 1. All asym* files for subjects where hemi=both. These will be used to
+#    generate an asymmetric brain template.
+# 2. All right-sym* files for subjects where use_right is True, and
+#    all left-sym* files for subjects where use_left is True.
+#    These will be used to generate a symmetric brain template.
+# 3. All right-hemi* files for subjects where use_right is True,
+#    and all left-hemi-xflip* files for subjects where use_left is True.
+#    These will be used to generate a symmetric hemisphere template.
+
+filepath_lists: dict[str, list] = {
+    "asym-brain": [],
+    "asym-mask": [],
+    "sym-brain": [],
+    "sym-mask": [],
+    "hemi-brain": [],
+    "hemi-mask": [],
+}
+
+for _, row in df.iterrows():
+    subject = "sub-" + row["subject_id"]
+    use4template_dir = use4template_dirs[subject]
+
+    if row["hemi"] == "both":
+        # Add paths for the asymmetric brain template
+        for label in ["brain", "mask"]:
+            filepath_lists[f"asym-{label}"].append(
+                use4template_dir / f"{subject}_asym-{label}.nii.gz"
+            )
+
+    if row["use_right"]:
+        for label in ["brain", "mask"]:
+            # Add paths for the symmetric brain template
+            filepath_lists[f"sym-{label}"].append(
+                use4template_dir / f"{subject}_right-sym-{label}.nii.gz"
+            )
+            # Add paths for the hemispheric template
+            filepath_lists[f"hemi-{label}"].append(
+                use4template_dir / f"{subject}_right-hemi-{label}.nii.gz"
+            )
+
+    if row["use_left"]:
+        for label in ["brain", "mask"]:
+            # Add paths for the symmetric brain template
+            filepath_lists[f"sym-{label}"].append(
+                use4template_dir / f"{subject}_left-sym-{label}.nii.gz"
+            )
+            # Add paths for the hemispheric template
+            filepath_lists[f"hemi-{label}"].append(
+                use4template_dir / f"{subject}_left-hemi-xflip-{label}.nii.gz"
+            )
+
+# %%
+# Save the file paths to text files, each in a separate directory
+
+for key, paths in filepath_lists.items():
+    kind, label = key.split("-")  # e.g. "asym" and "brain"
+    n_images = len(paths)
+    template_name = f"template_{kind}_{res_str}_n-{n_images}"
+    template_dir = species_dir / "templates" / template_name
+    template_dir.mkdir(exist_ok=True)
+    np.savetxt(template_dir / f"{label}_paths.txt", paths, fmt="%s")
diff --git a/examples/BlackCap_source_images.py b/examples/BlackCap_source_images.py
new file mode 100644
index 0000000..e48f4e8
--- /dev/null
+++ b/examples/BlackCap_source_images.py
@@ -0,0 +1,250 @@
+"""
+Identify source images for BlackCap template
+============================================
+Set up project directory with all relevant data and identify source images
+to be used for building the BlackCap template.
+"""
+
+# %%
+# Imports
+# -------
+import os
+import shutil
+from datetime import date
+from pathlib import Path
+
+import pandas as pd
+from loguru import logger
+from tqdm import tqdm
+
+from brainglobe_template_builder.io import get_unique_folder_in_dir
+
+# %%
+# Set up directories and logging
+# ------------------------------
+
+# Prepare directory structure
+atlas_dir = Path("/ceph/neuroinformatics/neuroinformatics/atlas-forge")
+species_id = "BlackCap"
+species_common_name = "Black cap"  # as in the source csv table
+species_dir = atlas_dir / species_id
+species_dir.mkdir(parents=True, exist_ok=True)
+# Make "rawdata", "derivatives", "templates", and "logs" directories
+for folder in ["rawdata", "derivatives", "templates", "logs"]:
+    (species_dir / folder).mkdir(exist_ok=True)
+
+# Set up logging
+today = date.today()
+current_script_name = os.path.basename(__file__).replace(".py", "")
+logger.add(species_dir / "logs" / f"{today}_{current_script_name}.log")
+logger.info(f"Will save outputs to {species_dir}.")
+
+# %%
+# Load a dataframe with all SWC brains used for atlases
+# ------------------------------------------------------
+
+source_csv_dir = Path.cwd() / "data"  # relative to repo root
+source_csv_path = source_csv_dir / "SWC_brain-list_for-atlases_2024-04-15.csv"
+df = pd.read_csv(source_csv_path)
+
+# Strip trailing space from str columns and from column names
+df = df.apply(lambda x: x.str.strip() if x.dtype == "object" else x)
+df.columns = df.columns.str.strip()
+# Select only data for the species of interest
+df = df[df["Common name"] == species_common_name]
+logger.info(f"Found {len(df)} {species_id} subjects.")
+
+# %%
+# Exract subject IDs and paths
+# ----------------------------
+
+df["subject_id"] = (
+    df["Specimen ID"]
+    .apply(lambda x: x[3:8].strip())
+    .apply(lambda x: x[:-1] if x.endswith("w") else x)
+)
+# Assert that subject IDs are unique
+assert len(df["subject_id"].unique()) == len(df), "Non-unique subject IDs"
+
+data_path_col = "Data path (raw)"
+df[data_path_col] = df[data_path_col].apply(
+    lambda x: x.replace("nfs/ceph", "ceph")
+)
+df["subject_path"] = df.apply(
+    lambda x: get_unique_folder_in_dir(
+        Path("/" + x[data_path_col]), x["subject_id"]
+    ),
+    axis=1,
+)
+
+sub_ids = df["subject_id"].values
+logger.info(f"Identified {sub_ids.size} unique subject IDs. {sub_ids}")
+sub_paths_msg = "Subject paths:\n"
+for idx, row in df.iterrows():
+    sub_id = row["subject_id"]
+    sub_path = row["subject_path"].as_posix()
+    sub_paths_msg += f"  {sub_id}: {sub_path}\n"
+logger.debug(sub_paths_msg)
+
+# %%
+# Identify all downsampled images
+# -------------------------------
+# We will identify all downsampled images for each subject, by iterating over
+# the downsampled_stacks folders and the various channels.
+# We will aggregate the image information and paths in a dataframe.
+
+valid_colors = ["far_red", "red", "green", "blue"]  # order matters
+
+images_list_of_dicts = []
+for _, row in df.iterrows():
+    sub = row["subject_id"]
+    sub_path_down = row["subject_path"] / "downsampled_stacks"
+    stacks = [f for f in os.listdir(sub_path_down) if f.endswith("micron")]
+
+    for stack in sorted(stacks):
+        microns = int(stack.split("_")[0])
+        stack_path = sub_path_down / stack
+
+        images = [
+            f
+            for f in os.listdir(stack_path)
+            if f.endswith(".tif") and "_ch0" in f
+        ]
+        for img in sorted(images):
+            ch_num = int(img.split("_ch0")[1].split("_")[0])
+            ch_color = [
+                c for c in valid_colors if img.split(".tif")[0].endswith(c)
+            ][0]
+            ch_color = ch_color if ch_color != "far_red" else "farred"
+            image_id = f"sub-{sub}_res-{microns}um_channel-{ch_color}"
+            images_list_of_dicts.append(
+                {
+                    "subject_id": sub,
+                    "microns": microns,
+                    "channel": ch_num,
+                    "color": ch_color,
+                    "image_id": image_id,
+                    "image_path": stack_path / img,
+                }
+            )
+
+images_df = pd.DataFrame(images_list_of_dicts)
+n_img = len(images_df)
+logger.info(f"Found {n_img} images across subjects, resolutions and channels.")
+
+# %%
+# Check for missing downsampled stacks
+# -------------------------------------
+# Logs a warning if any of the expected resolutions are missing
+
+expected_microns = [10, 25, 50]
+for sub in df["subject_id"]:
+    for microns in expected_microns:
+        if not (
+            images_df[
+                (images_df["subject_id"] == sub)
+                & (images_df["microns"] == microns)
+            ].shape[0]
+        ):
+            logger.warning(f"Subject {sub} lacks {microns} micron stack")
+
+# %%
+# Save dataframes and images to rawdata
+# -------------------------------------
+# Save the dataframes to the species rawdata directory.
+
+rawdata_dir = species_dir / "rawdata"
+subjects_csv = rawdata_dir / f"{today}_subjects.csv"
+df.to_csv(subjects_csv, index=False)
+logger.info(f"Saved subject information to csv: {subjects_csv}")
+
+images_csv = rawdata_dir / f"{today}_images.csv"
+images_df.to_csv(images_csv, index=False)
+logger.info(f"Saved image information to csv: {images_csv}")
+
+# %%
+# Save images to the species rawdata directory
+# (doesn't overwrite existing images).
+# High-resolution images (10um) are just symlinked to avoid duplication
+# of large files.
+
+n_copied, n_symlinked = 0, 0
+for idx in tqdm(images_df.index):
+    row = images_df.loc[idx, :]
+    sub = row["subject_id"]
+    sub_dir = rawdata_dir / f"sub-{sub}"
+    sub_dir.mkdir(exist_ok=True)
+    microns = row["microns"]
+    image_id = row["image_id"]
+    img_source_path = row["image_path"]
+    img_dest_path = sub_dir / f"{image_id}.tif"
+
+    # if the destination path exists, skip
+    if img_dest_path.exists():
+        logger.debug(f"Skipping {img_dest_path} as it already exists.")
+        continue
+
+    # if image is 10 microns, symlink it (they are large)
+    if microns == 10:
+        img_dest_path.symlink_to(img_source_path)
+        logger.debug(f"Symlinked {img_dest_path} to {img_source_path}")
+        n_symlinked += 1
+    # else copy the image
+    else:
+        shutil.copyfile(img_source_path, img_dest_path)
+        logger.debug(f"Copied {img_source_path} to {img_dest_path}")
+        n_copied += 1
+
+logger.info(
+    f"Copied {n_copied} and symlinked {n_symlinked} "
+    f"images to {rawdata_dir}."
+)
+
+# %%
+# Select data to use for template creation
+# ----------------------------------------
+# Here we select the subset of subjects to use for template creation.
+# The selection is based on visual inspection of the lowest-resolution
+# images (here 50um).
+#
+# The selected subjects are supplied as a .csv file named
+# "use_for_template.csv" which should be placed in the derivatives folder.
+
+# The resolution to use for visual inspection + preprocessing
+lowres = "50um"
+
+template_dir = species_dir / "templates"
+deriv_dir = species_dir / "derivatives"
+
+# Check for the .csv file
+use_for_template_csv = template_dir / "use_for_template.csv"
+if use_for_template_csv.exists():
+    use_for_template = pd.read_csv(use_for_template_csv)
+
+    assert (
+        use_for_template.subject_id.is_unique
+    ), "Subject IDs must be unique in use_for_template.csv"
+    n_subjects_to_use = len(use_for_template)
+    logger.info(
+        f"Will use {n_subjects_to_use} subjects for template creation."
+    )
+
+    for _, row in use_for_template.iterrows():
+        sub = row["subject_id"]
+        sub_str = f"sub-{sub}"
+        channel = row["color"]
+        hemi = row["hemi"]
+        filename = f"{sub_str}_res-{lowres}_channel-{channel}.tif"
+        # Create derivatives folder for each selected subject
+        (deriv_dir / sub_str).mkdir(parents=True, exist_ok=True)
+        logger.info(f"Created derivatives folder for {sub_str}")
+        # copy the low-res image from rawdata to derivatives for convenience
+        raw_file = rawdata_dir / sub_str / filename
+        derivative_file = deriv_dir / sub_str / filename
+        if derivative_file.is_file():
+            logger.info(f"Skipping {derivative_file} as it already exists.")
+        else:
+            shutil.copyfile(raw_file, derivative_file)
+            logger.info(
+                f"Copied {filename} to the {sub_str} derivatives folder."
+            )
diff --git a/examples/slurm_batch_across_average_types.sh b/examples/slurm_batch_across_average_types.sh
new file mode 100644
index 0000000..b697e1c
--- /dev/null
+++ b/examples/slurm_batch_across_average_types.sh
@@ -0,0 +1,47 @@
+#!/bin/bash
+
+#SBATCH -J ave_type # job name
+#SBATCH -p gpu # partition
+#SBATCH -N 1   # number of nodes
+#SBATCH --mem 64G # memory pool for all cores
+#SBATCH -n 20 # number of cores
+#SBATCH -t 3-00:00 # time (D-HH:MM)
+#SBATCH -o slurm.%x.%N.%j.out # write STDOUT
+#SBATCH -e slurm.%x.%N.%j.err # write STDERR
+#SBATCH --mail-type=ALL
+#SBATCH --mail-user=ucqfnsi@ucl.ac.uk
+#SBATCH --array=0-2
+
+# Load the required modules
+module load miniconda
+module load ants
+
+# Activate the conda environment
+conda activate template-builder
+
+# Define the atlas-forge directory and species name
+ATLAS_DIR="/ceph/neuroinformatics/neuroinformatics/atlas-forge"
+SPECIES="BlackCap"
+
+# Define the average types as an array
+AVE_TYPES=(median trimmed_mean efficient_trimean)
+AVE_TYPE=${AVE_TYPES[$SLURM_ARRAY_TASK_ID]}
+
+# Original template created via mean averaging
+TEMP_NAME_PREFIX="template_hemi_res-50um_n-18"
+# Get only the first 4 characters of the average type
+AVE_TYPE_SHORT=$(echo "$AVE_TYPE" | cut -c1-4)
+# Define the new template name
+TEMP_NAME="${TEMP_NAME_PREFIX}_average-${AVE_TYPE_SHORT}"
+
+# Create the working directory where the template will be built
+mkdir -p "${ATLAS_DIR}/${SPECIES}/templates/${TEMP_NAME}"
+# Copy brain_paths.txt and mask_paths.txt from the original template
+cp "${ATLAS_DIR}/${SPECIES}/templates/${TEMP_NAME_PREFIX}/brain_paths.txt" "${ATLAS_DIR}/${SPECIES}/templates/${TEMP_NAME}"
+cp "${ATLAS_DIR}/${SPECIES}/templates/${TEMP_NAME_PREFIX}/mask_paths.txt" "${ATLAS_DIR}/${SPECIES}/templates/${TEMP_NAME}"
+
+# Path to the bash script that builds the template
+BUILD_SCRIPT="${ATLAS_DIR}/${SPECIES}/scripts/BlackCap_build_template.sh"
+
+# Run the script to build the template
+bash "$BUILD_SCRIPT" --atlas-dir "$ATLAS_DIR" --species "$SPECIES" --template-name "$TEMP_NAME" --average-type "$AVE_TYPE"
diff --git a/examples/slurm_batch_across_templates.sh b/examples/slurm_batch_across_templates.sh
new file mode 100644
index 0000000..0cf4189
--- /dev/null
+++ b/examples/slurm_batch_across_templates.sh
@@ -0,0 +1,37 @@
+#!/bin/bash
+
+#SBATCH -J temp_build # job name
+#SBATCH -p cpu # partition
+#SBATCH -N 1   # number of nodes
+#SBATCH --mem 128G # memory pool for all cores
+#SBATCH -n 20 # number of cores
+#SBATCH -t 5-00:00 # time (D-HH:MM)
+#SBATCH -o slurm.%x.%N.%j.out # write STDOUT
+#SBATCH -e slurm.%x.%N.%j.err # write STDERR
+#SBATCH --mail-type=ALL
+#SBATCH --mail-user=n.sirmpilatze@ucl.ac.uk
+#SBATCH --array=0-2
+
+# Load the required modules
+module load miniconda
+module load ants
+# Activate the conda environment
+conda activate template-builder
+
+# Define the atlas-forge directory and species name
+ATLAS_DIR="/ceph/neuroinformatics/neuroinformatics/atlas-forge"
+SPECIES="BlackCap"
+
+# Define the template names as an array
+TEMP_NAME1="template_asym_res-50um_n-8"
+TEMP_NAME2="template_sym_res-50um_n-18"
+TEMP_NAME3="template_hemi_res-50um_n-18"
+TEMP_NAMES=($TEMP_NAME1 $TEMP_NAME2 $TEMP_NAME3)
+
+CURRENT_TEMP_NAME=${TEMP_NAMES[$SLURM_ARRAY_TASK_ID]}
+
+# Path to the bash script that builds the template
+BUILD_SCRIPT="${ATLAS_DIR}/${SPECIES}/scripts/BlackCap_build_template.sh"
+
+# Run the script to build the template
+bash $BUILD_SCRIPT --atlas-dir $ATLAS_DIR --species $SPECIES --template-name $CURRENT_TEMP_NAME
diff --git a/examples/slurm_build_template.sh b/examples/slurm_build_template.sh
new file mode 100644
index 0000000..fb557bb
--- /dev/null
+++ b/examples/slurm_build_template.sh
@@ -0,0 +1,30 @@
+#!/bin/bash
+
+#SBATCH -J temp_build # job name
+#SBATCH -p cpu # partition
+#SBATCH -N 1   # number of nodes
+#SBATCH --mem 64G # memory pool for all cores
+#SBATCH -n 20 # number of cores
+#SBATCH -t 5-00:00 # time (D-HH:MM)
+#SBATCH -o slurm.%x.%N.%j.out # write STDOUT
+#SBATCH -e slurm.%x.%N.%j.err # write STDERR
+#SBATCH --mail-type=ALL
+#SBATCH --mail-user=ucqfnsi@ucl.ac.uk
+
+# Load the required modules
+module load miniconda
+module load ants
+# Activate the conda environment
+conda activate template-builder
+
+# Define atlas-forge directory, species and template names, and average type
+ATLAS_DIR="/ceph/neuroinformatics/neuroinformatics/atlas-forge"
+SPECIES="BlackCap"
+TEMP_NAME="template_sym_res-50um_n-18"
+AVE_TYPE="mean"
+
+# Path to the bash script that builds the template
+BUILD_SCRIPT="${ATLAS_DIR}/${SPECIES}/scripts/BlackCap_build_template.sh"
+
+# Run the script to build the template
+bash $BUILD_SCRIPT --atlas-dir $ATLAS_DIR --species $SPECIES --template-name $TEMP_NAME --average-type $AVE_TYPE
diff --git a/examples/slurm_post_build.sh b/examples/slurm_post_build.sh
new file mode 100644
index 0000000..3abb626
--- /dev/null
+++ b/examples/slurm_post_build.sh
@@ -0,0 +1,21 @@
+#!/bin/bash
+
+#SBATCH -J postbuild # job name
+#SBATCH -p cpu # partition
+#SBATCH -N 1   # number of nodes
+#SBATCH --mem 64G # memory pool for all cores
+#SBATCH -n 16 # number of cores
+#SBATCH -t 1-00:00 # time (D-HH:MM)
+#SBATCH -o slurm.%x.%N.%j.out # write STDOUT
+#SBATCH -e slurm.%x.%N.%j.err # write STDERR
+#SBATCH --mail-type=ALL
+#SBATCH --mail-user=niko.sirbiladze@gmail.com
+
+# Load the required modules
+module load miniconda
+module load ants
+# Activate the conda environment
+conda activate atlas-forge
+
+# Run the post-build Python script
+python BlackCap_post_build.py
diff --git a/pyproject.toml b/pyproject.toml
index cc8a916..71def27 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,32 +1,34 @@
 [project]
 name = "brainglobe-template-builder"
-authors = [{name = "Niko Sirmpilatze", email= "niko.sirbiladze@gmail.com"}]
+authors = [{ name = "Niko Sirmpilatze", email = "niko.sirbiladze@gmail.com" }]
 description = "Build unbiased anatomical templates from individual images"
 readme = "README.md"
-requires-python = ">=3.9.0"
+requires-python = ">=3.10.0"
 dynamic = ["version"]
 
-license = {text = "BSD-3-Clause"}
+license = { text = "BSD-3-Clause" }
 
 classifiers = [
     "Development Status :: 2 - Pre-Alpha",
     "Programming Language :: Python",
     "Programming Language :: Python :: 3",
-    "Programming Language :: Python :: 3.9",
     "Programming Language :: Python :: 3.10",
     "Programming Language :: Python :: 3.11",
+    "Programming Language :: Python :: 3.12",
     "Operating System :: OS Independent",
     "License :: OSI Approved :: BSD License",
 ]
 
 dependencies = [
-    "numpy",
-    "scipy",
-    "simpleitk",
-    "vtk",
-    "qbatch"
+    "napari[all]==0.4.19",
+    "brainglobe-utils>=0.5",
+    "brainglobe-atlasapi>=2.0.7",
+    "loguru",
+    "antspyx",
 ]
 
+entry-points."napari.manifest".brainglobe-template-builder = "brainglobe_template_builder.napari:napari.yaml"
+
 [project.urls]
 "Homepage" = "https://github.com/brainglobe/brainglobe-template-builder"
 "Bug Tracker" = "https://github.com/brainglobe/brainglobe-template-builder/issues"
@@ -36,23 +38,19 @@ dependencies = [
 
 [project.optional-dependencies]
 dev = [
-  "pytest",
-  "pytest-cov",
-  "coverage",
-  "tox",
-  "black",
-  "mypy",
-  "pre-commit",
-  "ruff",
-  "setuptools_scm",
+    "pytest",
+    "pytest-cov",
+    "coverage",
+    "tox",
+    "black",
+    "mypy",
+    "pre-commit",
+    "ruff",
+    "setuptools_scm",
 ]
 
 [build-system]
-requires = [
-    "setuptools>=45",
-    "wheel",
-    "setuptools_scm[toml]>=6.2",
-]
+requires = ["setuptools>=45", "wheel", "setuptools_scm[toml]>=6.2"]
 build-backend = "setuptools.build_meta"
 
 [tool.setuptools]
@@ -67,7 +65,7 @@ exclude = ["tests", "docs*"]
 addopts = "--cov=brainglobe_template_builder"
 
 [tool.black]
-target-version = ['py39', 'py310', 'py311']
+target-version = ['py310', 'py311', 'py312']
 skip-string-normalization = false
 line-length = 79
 
@@ -75,39 +73,34 @@ line-length = 79
 
 [tool.check-manifest]
 ignore = [
-  ".yaml",
-  "tox.ini",
-  "tests/",
-  "tests/test_unit/",
-  "tests/test_integration/",
-  "docs/",
-  "docs/source/",
+    ".yaml",
+    "tox.ini",
+    "tests/",
+    "tests/test_unit/",
+    "tests/test_integration/",
+    "docs/",
+    "docs/source/",
+    "scripts/",
 ]
 
 
 [tool.ruff]
 line-length = 79
-exclude = ["__init__.py","build",".eggs"]
+exclude = ["__init__.py", "build", ".eggs", "scripts"]
 select = ["I", "E", "F"]
 fix = true
 
-[tool.cibuildwheel]
-build = "cp39-* cp310-* cp311-*"
-
-[tool.cibuildwheel.macos]
-archs = ["x86_64", "arm64"]
-
 [tool.tox]
 legacy_tox_ini = """
 [tox]
-envlist = py{39,310,311}
+envlist = py{310,311,312}
 isolated_build = True
 
 [gh-actions]
 python =
-    3.9: py39
     3.10: py310
     3.11: py311
+    3.12: py312
 
 [testenv]
 extras =


Test patch:
diff --git a/.github/workflows/test_and_deploy.yml b/.github/workflows/test_and_deploy.yml
index 2a4ba81..223d73c 100644
--- a/.github/workflows/test_and_deploy.yml
+++ b/.github/workflows/test_and_deploy.yml
@@ -27,16 +27,16 @@ jobs:
     strategy:
       matrix:
         # Run all supported Python versions on linux
-        python-version: ["3.9", "3.10", "3.11"]
+        python-version: ["3.10", "3.11", "3.12"]
         os: [ubuntu-latest]
         # Include one windows and ARM/Intel macos run
         include:
         - os: macos-13 # Intel Mac
-          python-version: "3.11"
+          python-version: "3.12"
         - os: macos-latest # ARM Mac
-          python-version: "3.11"
+          python-version: "3.12"
         - os: windows-latest
-          python-version: "3.10"
+          python-version: "3.11"
 
     steps:
       # Run tests
diff --git a/tests/test_unit/test_transform_utils.py b/tests/test_unit/test_transform_utils.py
new file mode 100644
index 0000000..0785d76
--- /dev/null
+++ b/tests/test_unit/test_transform_utils.py
@@ -0,0 +1,49 @@
+import dask.array as da
+import numpy as np
+import pytest
+from skimage import transform
+
+from brainglobe_template_builder.preproc.transform_utils import (
+    downsample_anisotropic_image_stack,
+)
+
+
+@pytest.fixture()
+def stack():
+    """Create a dask array representing an image stack"""
+    rng = np.random.default_rng()
+    data = rng.random(size=(10, 100, 100))
+    return da.from_array(data, chunks=(1, 100, 100))
+
+
+@pytest.fixture()
+def not_slicewise_stack(stack):
+    """Create a dask array representing an image stack"""
+    return stack.rechunk({0: 2, 1: 100, 2: 100})
+
+
+def test_downsample_anisotropic_image_stack(stack):
+    """Test that downsampling with dask gives same as without."""
+    in_plane = 20
+    axial = 2
+
+    downsampled_stack = downsample_anisotropic_image_stack(
+        stack, in_plane, axial
+    )
+
+    assert downsampled_stack.shape == (5, 5, 5)
+
+    expected = transform.downscale_local_mean(
+        stack.compute(), (1, in_plane, in_plane)
+    )
+    expected = transform.downscale_local_mean(expected, (axial, 1, 1))
+    assert np.all(
+        downsampled_stack == expected
+    ), "dask downsampling does not match expected skimage result"
+
+
+def test_downsample_anisotropic_image_stack_raises(not_slicewise_stack):
+    with pytest.raises(AssertionError, match="not chunked by plane!"):
+        downsample_anisotropic_image_stack(
+            not_slicewise_stack, in_plane_factor=20, axial_factor=2
+        )