Call BigStitcher by command line and read the results (#7)

* Add macro that calls big stitcher

* Add function to call BigStitcher via the command line

* Add functions to read the output of BigStitcher

* Add stitch button to napari widget

* Added tests for big_stitcher_bridge

* Added tests for file_utils

* Fixed tests to account for new defaults in stitch function

* Added tests for new functions in image_mosaic

* Added tests for new ImageMosaic and StitchingWidget functions

* Updated docstring from big_stitcher_bridge

* Updated docs for file_utils

* Added docstrings for ImageMosaic and StitchingWidget

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* WIP apply suggestions from code review

* More code review changes

* Move test constants out of test_image_mosaic to conftest

* Moved constants out of test_file_utils

* Moved constants out of test_big_stitcher_bridge

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* Make sure all ImageMosaic objects are explicitly cleaned up

* Update error for wrong imagej path

* WIP adding docstrings to test_stitching_widgets

* Added docstrings to test_stitching_widget

* Added comments for conftest.py

* Added docstrings for test_big_stitcher_bridge

* Added tests for test_image_mosaic

* Add docstring to test_file_utils

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

MANIFEST.in

@@ -1,6 +1,7 @@
 include LICENSE
 include README.md
 include brainglobe_stitch/napari.yaml
+include brainglobe_stitch/bigstitcher_macro.ijm
 exclude .pre-commit-config.yaml
 
 recursive-include brainglobe_stitch *.py

brainglobe_stitch/big_stitcher_bridge.py

@@ -0,0 +1,67 @@
+import subprocess
+from pathlib import Path
+from platform import system
+
+
+def run_big_stitcher(
+    imagej_path: Path,
+    xml_path: Path,
+    tile_config_path: Path,
+    all_channels: bool = False,
+    selected_channel: int = 488,
+    downsample_x: int = 4,
+    downsample_y: int = 4,
+    downsample_z: int = 1,
+) -> subprocess.CompletedProcess:
+    """
+    Run the BigStitcher ImageJ macro. Output is captured and returned as part
+    of the subprocess.CompletedProcess.
+
+    Parameters
+    ----------
+    imagej_path : Path
+        The path to the ImageJ executable.
+    xml_path : Path
+        The path to the BigDataViewer XML file.
+    tile_config_path : Path
+        The path to the BigStitcher tile configuration file.
+    all_channels : bool, optional
+        Whether to stitch based on all channels (default False).
+    selected_channel : int, optional
+        The channel on which to base the stitching (default 488).
+    downsample_x : int, optional
+        The downsample factor in the x-dimension for the stitching (default 4).
+    downsample_y : int, optional
+        The downsample factor in the y-dimension for the stitching (default 4).
+    downsample_z : int, optional
+        The downsample factor in the z-dimension for the stitching (default 1).
+
+    Returns
+    -------
+    subprocess.CompletedProcess
+        The result of the subprocess run.
+
+    Raises
+    ------
+    subprocess.CalledProcessError
+        If the subprocess returns a non-zero exit status.
+    """
+    stitch_macro_path = (
+        Path(__file__).resolve().parent / "bigstitcher_macro.ijm"
+    )
+
+    if system().startswith("Darwin"):
+        imagej_path = imagej_path / "Contents/MacOS/ImageJ-macosx"
+
+    command = (
+        f"{imagej_path} --ij2 "
+        f"--headless -macro {stitch_macro_path} "
+        f'"{xml_path} {tile_config_path} {int(all_channels)} '
+        f'{selected_channel} {downsample_x} {downsample_y} {downsample_z}"'
+    )
+
+    result = subprocess.run(
+        command, capture_output=True, text=True, check=True, shell=True
+    )
+
+    return result

brainglobe_stitch/bigstitcher_macro.ijm

@@ -0,0 +1,57 @@
+args = getArgument();
+args = split(args, " ");
+xmlPath = args[0];
+tilePath = args[1];
+allChannels = args[2];
+selectedChannel = args[3];
+downSampleX = args[4];
+downSampleY = args[5];
+downSampleZ = args[6];
+
+print("Stitching " + xmlPath);
+print("Loading TileConfiguration from " + tilePath)
+run(
+    "Load TileConfiguration from File...",
+    "browse=" + xmlPath +
+    " select=" + xmlPath +
+    " tileconfiguration=" + tilePath +
+    " use_pixel_units keep_metadata_rotation"
+);
+
+print("Calculating pairwise shifts");
+if (allChannels == 1) {
+    run(
+        "Calculate pairwise shifts ...",
+        "select=" + xmlPath +
+        " process_angle=[All angles] process_channel=[All channels] process_illumination=[All illuminations]" +
+        " process_tile=[All tiles] process_timepoint=[All Timepoints]" +
+        " method=[Phase Correlation] channels=[Average Channels]" +
+        " downsample_in_x=" + downSampleX +
+        " downsample_in_y=" + downSampleY +
+        " downsample_in_z=" + downSampleZ
+    );
+} else {
+    run(
+        "Calculate pairwise shifts ...",
+        "select=" + xmlPath +
+        " process_angle=[All angles] process_channel=[All channels] process_illumination=[All illuminations]" +
+        " process_tile=[All tiles] process_timepoint=[All Timepoints]" +
+        " method=[Phase Correlation] channels=[use Channel " + selectedChannel + " nm]" +
+        " downsample_in_x=" + downSampleX +
+        " downsample_in_y=" + downSampleY +
+        " downsample_in_z=" + downSampleZ
+    );
+}
+
+print("Optimizing globally and applying shifts");
+run(
+    "Optimize globally and apply shifts ...",
+    "select=" + xmlPath +
+    " process_angle=[All angles] process_channel=[All channels] process_illumination=[All illuminations]" +
+    " process_tile=[All tiles] process_timepoint=[All Timepoints] relative=2.500 absolute=3.500" +
+    " global_optimization_strategy=[Two-Round using Metadata to align unconnected Tiles and iterative dropping of bad links]" +
+    " fix_group_0-0,"
+);
+
+print("Done");
+eval("script", "System.exit(0);");

brainglobe_stitch/file_utils.py

@@ -34,16 +34,16 @@ def create_pyramid_bdv_h5(
 
     Parameters
     ----------
-    input_file: Path
+    input_file : Path
         The path to the input HDF5 file.
-    yield_progress: bool, optional
+    yield_progress : bool, optional
         Whether to yield progress. If True, the function will yield the
         progress as a percentage.
-    resolutions_array: npt.NDArray, optional
+    resolutions_array : npt.NDArray, optional
         The downsampling factors to use for each resolution level.
         This is a 2D array where each row represents a resolution level and the
         columns represent the downsampling factors for x, y, and z.
-    subdivisions_array: npt.NDArray, optional
+    subdivisions_array : npt.NDArray, optional
         The size of the blocks at each resolution level.
         This is a 2D array where each row represents a resolution level and the
         columns represent the size of the blocks for x, y, and z.
@@ -113,7 +113,7 @@ def parse_mesospim_metadata(
 
     Parameters
     ----------
-    meta_file_name: Path
+    meta_file_name : Path
         The path to the h5_meta.txt file.
 
     Returns
@@ -163,7 +163,7 @@ def check_mesospim_directory(
 
     Parameters
     ----------
-    mesospim_directory: Path
+    mesospim_directory : Path
         The path to the mesoSPIM directory.
 
     Returns
@@ -203,9 +203,9 @@ def get_slice_attributes(
 
     Parameters
     ----------
-    xml_path: Path
+    xml_path : Path
         The path to the XML file.
-    tile_names: List[str]
+    tile_names : List[str]
         The names of the tiles.
 
     Returns
@@ -226,3 +226,140 @@ def get_slice_attributes(
         slice_attributes[name] = sub_dict
 
     return slice_attributes
+
+
+def get_big_stitcher_transforms(xml_path: Path) -> npt.NDArray:
+    """
+    Get the translations for each tile from a Big Data Viewer XML file.
+    The translations are calculated by BigStitcher.
+
+    Parameters
+    ----------
+    xml_path : Path
+        The path to the Big Data Viewer XML file.
+
+    Returns
+    -------
+    npt.NDArray
+        A numpy array of shape (num_tiles, num_dim) with the translations.
+        Each row corresponds to a tile and each column to a dimension.
+    """
+    tree = ET.parse(xml_path)
+    root = tree.getroot()
+
+    stitch_transforms = safe_find_all(
+        root, ".//ViewTransform/[Name='Stitching Transform']/affine"
+    )
+
+    # Stitching Transforms are there if aligning to grid is done manually
+    # Translation from Tile Configurations are there if aligned automatically
+    grid_transforms = safe_find_all(
+        root,
+        ".//ViewTransform/[Name='Translation from Tile Configuration']/affine",
+    )
+    if len(grid_transforms) == 0:
+        grid_transforms = safe_find_all(
+            root,
+            ".//ViewTransform/[Name='Translation to Regular Grid']/affine",
+        )
+
+    z_scale_str = safe_find(
+        root, ".//ViewTransform/[Name='calibration']/affine"
+    )
+
+    if not z_scale_str.text:
+        raise ValueError("No z scale found in XML")
+
+    z_scale = float(z_scale_str.text.split()[-2])
+
+    deltas = np.ones((len(stitch_transforms), 3))
+    grids = np.ones((len(grid_transforms), 3))
+    for i in range(len(stitch_transforms)):
+        delta_nums_text = stitch_transforms[i].text
+        grid_nums_text = grid_transforms[i].text
+
+        if not delta_nums_text or not grid_nums_text:
+            raise ValueError("No translation values found in XML")
+
+        delta_nums = delta_nums_text.split()
+        grid_nums = grid_nums_text.split()
+
+        # Extract the translation values from the transform.
+        # Swap the order of the axis (x,y,z) to (z,y,x).
+        # The input values are a flattened 4x4 matrix where
+        # the translation values in the last column.
+        deltas[i] = np.array(delta_nums[11:2:-4])
+        grids[i] = np.array(grid_nums[11:2:-4])
+
+    # Divide the z translation by the z scale
+    deltas[:, 0] /= z_scale
+    grids[:, 0] /= z_scale
+
+    # Round the translations to the nearest integer
+    grids = grids.round().astype(np.int32)
+    deltas = deltas.round().astype(np.int32)
+
+    # Normalise the grid transforms by subtracting the minimum value
+    norm_grids = grids - grids.min(axis=0)
+    # Calculate the maximum delta (from BigStitcher) for each dimension
+    max_delta = np.absolute(deltas).max(axis=0)
+
+    # Calculate the start and end coordinates for each tile such that the
+    # first tile is at 0,0,0 and provide enough padding to account for the
+    # transforms from BigStitcher
+    translations = norm_grids + deltas + max_delta
+
+    return translations
+
+
+def safe_find_all(root: ET.Element, query: str) -> List[ET.Element]:
+    """
+    Find all elements matching a query in an ElementTree root. If no
+    elements are found, return an empty list.
+
+    Parameters
+    ----------
+    root : ET.Element
+        The root of the ElementTree.
+    query : str
+        The query to search for.
+
+    Returns
+    -------
+    List[ET.Element]
+        A list of elements matching the query.
+    """
+    elements = root.findall(query)
+    if elements is None:
+        return []
+
+    return elements
+
+
+def safe_find(root: ET.Element, query: str) -> ET.Element:
+    """
+    Find the first element matching a query in an ElementTree root.
+    Raise a ValueError if no element found.
+
+    Parameters
+    ----------
+    root : ET.Element
+        The root of the ElementTree.
+    query : str
+        The query to search for.
+
+    Returns
+    -------
+    ET.Element
+        The element matching the query or None.
+
+    Raises
+    ------
+    ValueError
+        If no element is found.
+    """
+    element = root.find(query)
+    if element is None or element.text is None:
+        raise ValueError(f"No element found for query {query}")
+
+    return element