diff --git a/.github/workflows/test_and_deploy.yml b/.github/workflows/test_and_deploy.yml index 241341b..83da4e3 100644 --- a/.github/workflows/test_and_deploy.yml +++ b/.github/workflows/test_and_deploy.yml @@ -37,10 +37,16 @@ jobs: python-version: "3.10" steps: + # these libraries enable testing on Qt on linux + - uses: pyvista/setup-headless-display-action@v2 + with: + qt: true + # Run tests - uses: neuroinformatics-unit/actions/test@v2 with: python-version: ${{ matrix.python-version }} + secret-codecov-token: ${{ secrets.CODECOV_TOKEN }} build_sdist_wheels: name: Build source distribution diff --git a/MANIFEST.in b/MANIFEST.in index 2231169..9175ede 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -1,5 +1,6 @@ include LICENSE include README.md +include brainglobe_stitch/napari.yaml exclude .pre-commit-config.yaml recursive-include brainglobe_stitch *.py diff --git a/brainglobe_stitch/file_utils.py b/brainglobe_stitch/file_utils.py index 788825e..8aa95c0 100644 --- a/brainglobe_stitch/file_utils.py +++ b/brainglobe_stitch/file_utils.py @@ -1,4 +1,10 @@ +import xml.etree.ElementTree as ET from pathlib import Path +from typing import Dict, List, Tuple, Union + +import h5py +import numpy as np +import numpy.typing as npt HEADERS = [ "[POSITION]", @@ -8,70 +14,135 @@ ] -def write_big_stitcher_tile_config(meta_file_name: Path) -> list[dict]: - tile_metadata = parse_mesospim_metadata(meta_file_name) - - output_file = str(meta_file_name)[:-12] + "_tile_config.txt" - - first_channel = tile_metadata[0]["Laser"] - num_channels = 1 - - for metadata in tile_metadata[1:]: - if metadata["Laser"] == first_channel: - break - else: - num_channels += 1 - - num_tiles = len(tile_metadata) // num_channels - tile_xy_locations = [] - - for i in range(0, len(tile_metadata), num_channels): - curr_tile_dict = tile_metadata[i] - - x = round(curr_tile_dict["x_pos"] / curr_tile_dict["Pixelsize in um"]) - y = round(curr_tile_dict["y_pos"] / curr_tile_dict["Pixelsize in um"]) - - tile_xy_locations.append((x, y)) - - relative_locations = [(0, 0)] - - for abs_tuple in tile_xy_locations[1:]: - rel_tuple = ( - abs(abs_tuple[0] - tile_xy_locations[0][0]), - abs(abs_tuple[1] - tile_xy_locations[0][1]), - ) - relative_locations.append(rel_tuple) - - with open(output_file, "w") as f: - f.write("dim=3\n") - for i in range(len(tile_metadata)): - f.write( - f"{i};;" - f"({relative_locations[i%num_tiles][0]}," - f"{relative_locations[i%num_tiles][1]},0)\n" - ) - - return tile_metadata - - -def parse_mesospim_metadata(meta_file_name: Path): +def create_pyramid_bdv_h5( + input_file: Path, + yield_progress: bool = False, + resolutions_array: npt.NDArray = np.array( + [[1, 1, 1], [2, 2, 1], [4, 4, 1], [8, 8, 1], [16, 16, 1]] + ), + subdivisions_array: npt.NDArray = np.array( + [[32, 32, 16], [32, 32, 16], [32, 32, 16], [32, 32, 16], [32, 32, 16]] + ), +): + """ + Create a resolution pyramid for a Big Data Viewer HDF5 file. The function + assumes no pyramid exists and creates a new one in place. By default, + the function creates a 5 level pyramid with downsampling factors of 1, 2, + 4, 8, and 16 in x and y, with no downsampling in z. Deletes the old + resolutions and subdivisions datasets and creates new ones. + + + Parameters + ---------- + input_file: Path + The path to the input HDF5 file. + yield_progress: bool, optional + Whether to yield progress. If True, the function will yield the + progress as a percentage. + 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 + 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. + """ + with h5py.File(input_file, "r+") as f: + data_group = f["t00000"] + num_done = 0 + num_slices = len(data_group) + + for curr_slice in data_group: + # Delete the old resolutions and subdivisions datasets + del f[f"{curr_slice}/resolutions"] + del f[f"{curr_slice}/subdivisions"] + f[f"{curr_slice}/resolutions"] = resolutions_array + f[f"{curr_slice}/subdivisions"] = subdivisions_array + + grp: h5py.Group = f[f"t00000/{curr_slice}"] + for i in range(1, resolutions_array.shape[0]): + downsampling_factors = ( + resolutions_array[i] // resolutions_array[i - 1] + ) + prev_resolution = grp[f"{i - 1}/cells"] + # Add 1 to account for odd dimensions + + grp.require_dataset( + f"{i}/cells", + dtype=prev_resolution.dtype, + # Data is stored in z,y,x, but the downsampling + # factors are in x,y,z, so need to reverse + # Adding 1 allows to account for dimensions of odd size, + # Only add 1 if the downsampling factor is greater than 1 + shape=( + ( + prev_resolution.shape[0] + + (downsampling_factors[2] > 1) + ) + // downsampling_factors[2], + ( + prev_resolution.shape[1] + + (downsampling_factors[1] > 1) + ) + // downsampling_factors[1], + ( + prev_resolution.shape[2] + + (downsampling_factors[0] > 1) + ) + // downsampling_factors[0], + ), + ) + grp[f"{i}/cells"][...] = prev_resolution[ + :: downsampling_factors[2], + :: downsampling_factors[1], + :: downsampling_factors[0], + ] + + num_done += 1 + + if yield_progress: + yield int(100 * num_done / num_slices) + + +def parse_mesospim_metadata( + meta_file_name: Path, +) -> List[Dict]: + """ + Parse the metadata from a mesoSPIM .h5_meta.txt file. + + Parameters + ---------- + meta_file_name: Path + The path to the h5_meta.txt file. + + Returns + ------- + List[Dict] + A list of dictionaries containing the metadata for each tile. + """ tile_metadata = [] with open(meta_file_name, "r") as f: lines = f.readlines() - curr_tile_metadata: dict[str, str | int | float] = {} + curr_tile_metadata: Dict[str, Union[str, int, float]] = {} for line in lines[3:]: line = line.strip() + # Tile metadata is separated by a line starting with [CFG] if line.startswith("[CFG"): tile_metadata.append(curr_tile_metadata) curr_tile_metadata = {} + # Skip the headers elif line in HEADERS: continue + # Skip empty lines elif not line: continue else: split_line = line.split("]") value = split_line[1].strip() + # Check if the value is an int or a float + # If it is neither, store it as a string if value.isdigit(): curr_tile_metadata[split_line[0][1:]] = int(value) else: @@ -82,3 +153,76 @@ def parse_mesospim_metadata(meta_file_name: Path): tile_metadata.append(curr_tile_metadata) return tile_metadata + + +def check_mesospim_directory( + mesospim_directory: Path, +) -> Tuple[Path, Path, Path]: + """ + Check that the mesoSPIM directory contains the expected files. + + Parameters + ---------- + mesospim_directory: Path + The path to the mesoSPIM directory. + + Returns + ------- + Tuple[Path, Path, Path] + The paths to the bdv.xml, h5_meta.txt, and bdv.h5 files. + """ + # List all files in the directory that do not start with a period + # But end in the correct file extension + xml_path = list(mesospim_directory.glob("[!.]*bdv.xml")) + meta_path = list(mesospim_directory.glob("[!.]*h5_meta.txt")) + h5_path = list(mesospim_directory.glob("[!.]*bdv.h5")) + + # Check that there is exactly one file of each type + if len(xml_path) != 1: + raise FileNotFoundError( + f"Expected 1 bdv.xml file, found {len(xml_path)}" + ) + + if len(meta_path) != 1: + raise FileNotFoundError( + f"Expected 1 h5_meta.txt file, found {len(meta_path)}" + ) + + if len(h5_path) != 1: + raise FileNotFoundError(f"Expected 1 h5 file, found {len(h5_path)}") + + return xml_path[0], meta_path[0], h5_path[0] + + +def get_slice_attributes( + xml_path: Path, tile_names: List[str] +) -> Dict[str, Dict]: + """ + Get the slice attributes from a Big Data Viewer XML file. Attributes + include the illumination id, channel id, and tile id, and angle id. + + Parameters + ---------- + xml_path: Path + The path to the XML file. + tile_names: List[str] + The names of the tiles. + + Returns + ------- + Dict[str, Dict] + A dictionary containing the slice attributes for each tile. + """ + tree = ET.parse(xml_path) + root = tree.getroot() + view_setups = root.findall(".//ViewSetup//attributes") + + slice_attributes = {} + for child, name in zip(view_setups, tile_names): + sub_dict = {} + for sub_child in child: + sub_dict[sub_child.tag] = sub_child.text + + slice_attributes[name] = sub_dict + + return slice_attributes diff --git a/brainglobe_stitch/image_mosaic.py b/brainglobe_stitch/image_mosaic.py new file mode 100644 index 0000000..a3d470f --- /dev/null +++ b/brainglobe_stitch/image_mosaic.py @@ -0,0 +1,263 @@ +from pathlib import Path +from typing import Dict, List, Optional, Tuple + +import dask.array as da +import h5py +import numpy as np +import numpy.typing as npt +from rich.progress import Progress + +from brainglobe_stitch.file_utils import ( + check_mesospim_directory, + create_pyramid_bdv_h5, + get_slice_attributes, + parse_mesospim_metadata, +) +from brainglobe_stitch.tile import Tile + + +class ImageMosaic: + """ + Class to represent an image as a collection of tiles. + + Attributes + ---------- + directory : Path + The directory containing the image data. + xml_path : Path | None + The path to the Big Data Viewer XML file. + meta_path : Path | None + The path to the mesoSPIM metadata file. + h5_path : Path | None + The path to the Big Data Viewer h5 file containing the raw data. + tile_config_path : Path | None + The path to the BigStitcher tile configuration file. + h5_file : h5py.File | None + An open h5py file object for the raw data. + channel_names : List[str] + The names of the channels in the image as strings. + tiles : List[Tile] + The tiles in the image. + tile_names : List[str] + The names of the image tiles from BigDataViewer. + x_y_resolution : float + The resolution of the image in the x and y dimensions + in micrometers per pixel. + z_resolution : float + The resolution of the image in the z dimension + in micrometers per pixel. + num_channels : int + The number of channels in the image. + """ + + def __init__(self, directory: Path): + self.directory: Path = directory + self.xml_path: Optional[Path] = None + self.meta_path: Optional[Path] = None + self.h5_path: Optional[Path] = None + self.tile_config_path: Optional[Path] = None + self.h5_file: Optional[h5py.File] = None + self.channel_names: List[str] = [] + self.tiles: List[Tile] = [] + self.tile_names: List[str] = [] + self.tile_metadata: List[Dict] = [] + self.x_y_resolution: float = 4.0 # um per pixel + self.z_resolution: float = 5.0 # um per pixel + self.num_channels: int = 1 + + self.load_mesospim_directory() + + def __del__(self): + if self.h5_file is not None: + self.h5_file.close() + self.h5_file = None + + def data_for_napari( + self, resolution_level: int = 0 + ) -> List[Tuple[da.Array, npt.NDArray]]: + """ + Return data for visualisation in napari. + + Parameters + ---------- + resolution_level: int + The resolution level to get the data for. + + Returns + ------- + List[Tuple[da.Array, npt.NDArray]] + A list of tuples containing the data and the translation for each + tile scaled to the selected resolution. + """ + data = [] + for tile in self.tiles: + scaled_tile = tile.data_pyramid[resolution_level] + scaled_translation = ( + np.array(tile.position) + // tile.resolution_pyramid[resolution_level] + ) + data.append((scaled_tile, scaled_translation)) + + return data + + def load_mesospim_directory(self) -> None: + """ + Load the mesoSPIM directory and its data into the ImageMosaic. + """ + try: + ( + self.xml_path, + self.meta_path, + self.h5_path, + ) = check_mesospim_directory(self.directory) + except FileNotFoundError: + print("Invalid mesoSPIM directory") + + assert self.xml_path is not None + assert self.meta_path is not None + assert self.h5_path is not None + + self.h5_file = h5py.File(self.h5_path, "r") + + # Check if resolution pyramid exists + # s00 should have more than 1 key if the resolution pyramid exists + # Each key in ["t00000/s00"] corresponds to a resolution level + if len(self.h5_file["t00000/s00"].keys()) <= 1: + print("Resolution pyramid not found.") + # Close the file as it's open as read only + self.h5_file.close() + print("Creating resolution pyramid.") + + # Create resolution pyramid + with Progress() as progress: + task = progress.add_task("Downsampling...", total=100) + + for update in create_pyramid_bdv_h5( + self.h5_path, + yield_progress=True, + ): + progress.update(task, advance=update) + + # Reopen the file + self.h5_file = h5py.File(self.h5_path, "r") + + self.tile_config_path = Path( + str(self.xml_path)[:-4] + "_tile_config.txt" + ) + + self.tile_metadata = parse_mesospim_metadata(self.meta_path) + + self.channel_names = [] + idx = 0 + while self.tile_metadata[idx]["Laser"] not in self.channel_names: + self.channel_names.append(self.tile_metadata[idx]["Laser"]) + idx += 1 + + self.x_y_resolution = self.tile_metadata[0]["Pixelsize in um"] + self.z_resolution = self.tile_metadata[0]["z_stepsize"] + self.num_channels = len(self.channel_names) + + # Each tile is a group under "t00000" + # Ordered in increasing order based on acquisition + # Names aren't always contiguous, e.g. s00, s01, s04, s05 is valid + tile_group = self.h5_file["t00000"] + self.tile_names = list(tile_group.keys()) + slice_attributes = get_slice_attributes(self.xml_path, self.tile_names) + + self.tiles = [] + for idx, tile_name in enumerate(self.tile_names): + tile = Tile(tile_name, idx, slice_attributes[tile_name]) + tile.channel_name = self.channel_names[tile.channel_id] + self.tiles.append(tile) + tile_data = [] + + for pyramid_level in tile_group[tile_name].keys(): + tile_data.append( + da.from_array( + tile_group[f"{tile_name}/{pyramid_level}/cells"], + ) + ) + + tile.data_pyramid = tile_data + + # Add the scaling factor for each resolution level of the pyramid. + resolutions = self.h5_file[f"{tile_name}/resolutions"] + tile.resolution_pyramid = np.ones( + (len(resolutions), 3), dtype=np.int16 + ) + + # Switch to z,y,x order from x,y,z order + for i, resolution in enumerate(resolutions): + tile.resolution_pyramid[i] = resolution[-1::-1] + + # Don't rewrite the tile config file if it already exists + # These will be used as the initial tile positions + # This file will be passed to BigStitcher + if not self.tile_config_path.exists(): + print("Tile positions not found. Writing tile config file.") + self.write_big_stitcher_tile_config(self.meta_path) + + # Read the tile config file to get the initial tile positions + with open(self.tile_config_path, "r") as f: + # Skip header + f.readline() + + for line, tile in zip(f.readlines(), self.tiles): + split_line = line.split(";")[-1].strip("()\n").split(",") + # BigStitcher uses x,y,z order + # Switch to z,y,x order + translation = [ + int(split_line[2]), + int(split_line[1]), + int(split_line[0]), + ] + tile.position = translation + + def write_big_stitcher_tile_config(self, meta_file_name: Path) -> None: + """ + Write the BigStitcher tile configuration file + (placement for each tile based on stage coordinates). + + Parameters + ---------- + meta_file_name: Path + The path to the mesoSPIM metadata file. + """ + # Remove .h5_meta.txt from the file name + output_file = str(meta_file_name)[:-12] + "_tile_config.txt" + + tile_xy_locations = [] + for i in range(0, len(self.tiles), self.num_channels): + curr_tile_dict = self.tile_metadata[i] + + # Get the x and y positions in pixels + x = round( + curr_tile_dict["x_pos"] / curr_tile_dict["Pixelsize in um"] + ) + y = round( + curr_tile_dict["y_pos"] / curr_tile_dict["Pixelsize in um"] + ) + + tile_xy_locations.append((x, y)) + + # Calculate relative pixel positions for each tile + # The first tile is at (0,0) + relative_locations = [(0, 0)] + for abs_tuple in tile_xy_locations[1:]: + rel_tuple = ( + abs(abs_tuple[0] - tile_xy_locations[0][0]), + abs(abs_tuple[1] - tile_xy_locations[0][1]), + ) + relative_locations.append(rel_tuple) + + # Write the tile config file based on what BigStitcher expects + with open(output_file, "w") as f: + f.write("dim=3\n") + for tile, tile_name in zip(self.tiles, self.tile_names): + f.write( + f"{tile_name[1:]};;" + f"({relative_locations[tile.tile_id][0]}," + f"{relative_locations[tile.tile_id][1]},0)\n" + ) + + return diff --git a/brainglobe_stitch/napari.yaml b/brainglobe_stitch/napari.yaml new file mode 100644 index 0000000..39ee91f --- /dev/null +++ b/brainglobe_stitch/napari.yaml @@ -0,0 +1,10 @@ +name: brainglobe-stitch +display_name: BrainGlobe Stitch +contributions: + commands: + - id: brainglobe-stitch.make_stitching_widget + python_name: brainglobe_stitch.stitching_widget:StitchingWidget + title: Make stitching widget + widgets: + - command: brainglobe-stitch.make_stitching_widget + display_name: BrainGlobe Stitch diff --git a/brainglobe_stitch/stitching_widget.py b/brainglobe_stitch/stitching_widget.py new file mode 100644 index 0000000..3b86ce2 --- /dev/null +++ b/brainglobe_stitch/stitching_widget.py @@ -0,0 +1,241 @@ +from pathlib import Path +from typing import List, Optional, Tuple + +import dask.array as da +import h5py +import napari +import numpy.typing as npt +from brainglobe_utils.qtpy.logo import header_widget +from napari import Viewer +from napari.qt.threading import create_worker +from napari.utils.notifications import show_warning +from qtpy.QtWidgets import ( + QFileDialog, + QHBoxLayout, + QLabel, + QLineEdit, + QProgressBar, + QPushButton, + QVBoxLayout, + QWidget, +) + +from brainglobe_stitch.file_utils import ( + check_mesospim_directory, + create_pyramid_bdv_h5, +) +from brainglobe_stitch.image_mosaic import ImageMosaic + + +def add_tiles_from_mosaic( + napari_data: List[Tuple[da.Array, npt.NDArray]], tile_names: List[str] +): + """ + Add tiles to the napari viewer from the ImageMosaic. + + Parameters + ------------ + napari_data : List[Tuple[da.Array, npt.NDArray]] + The data and position for each tile in the mosaic. + tile_names : List[str] + The list of tile names. + """ + + for data, tile_name in zip(napari_data, tile_names): + tile_data, tile_position = data + tile_layer = napari.layers.Image( + tile_data.compute(), + name=tile_name, + blending="translucent", + contrast_limits=[0, 4000], + multiscale=False, + ) + tile_layer.translate = tile_position + + yield tile_layer + + +class StitchingWidget(QWidget): + """ + napari widget for stitching large tiled 3d images. + + Parameters + ------------ + napari_viewer : napari.Viewer + The napari viewer to add the widget to. + + Attributes + ---------- + progress_bar : QProgressBar + The progress bar for the widget, reused for multiple function. + image_mosaic : Optional[ImageMosaic] + The ImageMosaic object representing the data that will be stitched. + tile_layers : List[napari.layers.Image] + The list of napari layers containing the tiles. + resolution_to_display : int + The resolution level of the pyramid to display in napari. + header : QWidget + The header widget for the StitchingWidget. + default_directory : Path + The default directory for the widget (home directory by default). + working_directory : Path + The working directory for the widget. + select_mesospim_directory : QWidget + The widget for selecting the mesoSPIM directory. + mesospim_directory_text_field : QLineEdit + The text field for the mesoSPIM directory. + open_file_dialog : QPushButton + The button for opening the file dialog. + create_pyramid_button : QPushButton + The button for creating the resolution pyramid. + add_tiles_button : QPushButton + The button for adding the tiles to the viewer. + """ + + def __init__(self, napari_viewer: Viewer): + super().__init__() + self._viewer = napari_viewer + self.progress_bar = QProgressBar(self) + self.image_mosaic: Optional[ImageMosaic] = None + self.tile_layers: List[napari.layers.Image] = [] + self.resolution_to_display: int = 3 + + self.setLayout(QVBoxLayout()) + + self.header = header_widget( + package_name="brainglobe-stitch", + package_tagline="Stitching mesoSPIM data", + ) + + self.layout().addWidget(self.header) + + self.default_directory = Path.home() + self.working_directory = self.default_directory + + self.select_mesospim_directory = QWidget() + self.select_mesospim_directory.setLayout(QHBoxLayout()) + + self.mesospim_directory_text_field = QLineEdit() + self.mesospim_directory_text_field.setText(str(self.default_directory)) + self.mesospim_directory_text_field.editingFinished.connect( + self._on_mesospim_directory_text_edited + ) + self.select_mesospim_directory.layout().addWidget( + self.mesospim_directory_text_field + ) + + self.open_file_dialog = QPushButton("Browse") + self.open_file_dialog.clicked.connect( + self._on_open_file_dialog_clicked + ) + self.select_mesospim_directory.layout().addWidget( + self.open_file_dialog + ) + + self.layout().addWidget(QLabel("Select mesospim directory:")) + self.layout().addWidget(self.select_mesospim_directory) + + self.create_pyramid_button = QPushButton("Create resolution pyramid") + self.create_pyramid_button.clicked.connect( + self._on_create_pyramid_button_clicked + ) + self.create_pyramid_button.setEnabled(False) + + self.layout().addWidget(self.create_pyramid_button) + + self.add_tiles_button = QPushButton("Add tiles to viewer") + self.add_tiles_button.clicked.connect( + self._on_add_tiles_button_clicked + ) + self.add_tiles_button.setEnabled(False) + self.layout().addWidget(self.add_tiles_button) + + def _on_open_file_dialog_clicked(self): + """ + Open a file dialog to select the mesoSPIM directory. + """ + self.working_directory = Path( + QFileDialog.getExistingDirectory( + self, "Select mesoSPIM directory", str(self.default_directory) + ) + ) + # Add the text to the mesospim directory text field + self.mesospim_directory_text_field.setText(str(self.working_directory)) + self.check_and_load_mesospim_directory() + + def _on_mesospim_directory_text_edited(self): + """ + Update the working directory when the text field is edited. + """ + self.working_directory = Path( + self.mesospim_directory_text_field.text() + ) + self.check_and_load_mesospim_directory() + + def _on_create_pyramid_button_clicked(self): + """ + Create the resolution pyramid for the input mesoSPIM h5 data. + """ + self.progress_bar.setValue(0) + self.progress_bar.setRange(0, 100) + + worker = create_worker( + create_pyramid_bdv_h5, + self.h5_path, + yield_progress=True, + ) + worker.yielded.connect(self.progress_bar.setValue) + worker.finished.connect(self.progress_bar.reset) + worker.start() + + self.create_pyramid_button.setEnabled(False) + self.add_tiles_button.setEnabled(True) + + def _on_add_tiles_button_clicked(self): + """ + Add the tiles from the mesoSPIM h5 file to the viewer. + """ + self.image_mosaic = ImageMosaic(self.working_directory) + + napari_data = self.image_mosaic.data_for_napari( + self.resolution_to_display + ) + + worker = create_worker( + add_tiles_from_mosaic, napari_data, self.image_mosaic.tile_names + ) + worker.yielded.connect(self._set_tile_layers) + worker.start() + + def _set_tile_layers(self, tile_layer: napari.layers.Image): + """ + Add the tile layer to the viewer and store it in the tile_layers list. + + Parameters + ---------- + tile_layer : napari.layers.Image + """ + tile_layer = self._viewer.add_layer(tile_layer) + self.tile_layers.append(tile_layer) + + def check_and_load_mesospim_directory(self): + """ + Check if the selected directory is a valid mesoSPIM directory, + if valid load the h5 file and check if the resolution pyramid + is present. If not present, enable the create pyramid button. + Otherwise, enable the add tiles button. + """ + try: + ( + self.xml_path, + self.meta_path, + self.h5_path, + ) = check_mesospim_directory(self.working_directory) + with h5py.File(self.h5_path, "r") as f: + if len(f["t00000/s00"].keys()) <= 1: + show_warning("Resolution pyramid not found") + self.create_pyramid_button.setEnabled(True) + else: + self.add_tiles_button.setEnabled(True) + except FileNotFoundError: + show_warning("mesoSPIM directory not valid") diff --git a/brainglobe_stitch/tile.py b/brainglobe_stitch/tile.py new file mode 100644 index 0000000..73daa34 --- /dev/null +++ b/brainglobe_stitch/tile.py @@ -0,0 +1,57 @@ +from typing import Dict, List, Optional + +import dask.array as da +import numpy as np +import numpy.typing as npt + + +class Tile: + """ + Tile class to store information about a single tile of an image. + + Attributes + ---------- + name : str + The name of the tile. + id : int + The id of the tile (unique for each tile). + position : List[int] + The position of the tile in the fused image. Initialized to 0, 0, 0. + neighbours : List[int] + The ids of the neighbouring tiles. + data_pyramid : List[da.Array] + The pyramid of data arrays for the tile. The full resolution data + is at index 0. + resolution_pyramid : npt.NDArray + The downsample factors for each pyramid level. + channel_id : int + The id of the channel. + channel_name : str + The name of the channel. + tile_id : int + The id of the tile based on its position in the image. + Two tiles with the same tile_id are in the same position in + different channels. + illumination_id : int + The id of the illumination side. + angle : float + The angle of sample rotation. + """ + + def __init__( + self, + tile_name: str, + tile_id: int, + attributes: Dict[str, str], + ): + self.name: str = tile_name + self.id: int = tile_id + self.position: List[int] = [0, 0, 0] + self.neighbours: List[int] = [] + self.data_pyramid: List[da.Array] = [] + self.resolution_pyramid: npt.NDArray = np.array([]) + self.channel_name: Optional[str] = None + self.channel_id: int = int(attributes["channel"]) + self.tile_id: int = int(attributes["tile"]) + self.illumination_id: int = int(attributes["illumination"]) + self.angle: float = float(attributes["angle"]) diff --git a/pyproject.toml b/pyproject.toml index 00d127f..106015f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -21,6 +21,18 @@ classifiers = [ "Topic :: Scientific/Engineering :: Image Processing", ] +dependencies = [ + "napari>=0.4.18", + "brainglobe-utils>=0.3.4", + "h5py", + "napari-ome-zarr", + "ome-zarr", + "zarr", + "numpy", + "qtpy", + "tifffile" +] + [project.urls] "Homepage" = "https://brainglobe.info" "Bug Tracker" = "https://github.com/brainglobe/brainglobe-stitch/issues" @@ -29,12 +41,19 @@ classifiers = [ "User Support" = "https://forum.image.sc/tag/brainglobe" +[project.entry-points."napari.manifest"] +mesospim-stitcher = "brainglobe_stitch:napari.yaml" + [project.optional-dependencies] dev = [ "pytest", "pytest-cov", + "pytest-mock", + "pytest-qt", + "pyqt5", "coverage", "tox", + "pooch", "black", "mypy", "pre-commit", @@ -99,8 +118,19 @@ python = 3.11: py311 [testenv] +platform = + macos: darwin + linux: linux + windows: win32 +passenv = + CI + GITHUB_ACTIONS + DISPLAY + XAUTHORITY + NUMPY_EXPERIMENTAL_ARRAY_FUNCTION + PYVISTA_OFF_SCREEN extras = dev commands = - pytest -v --color=yes --cov=brainglobe_stitch --cov-report=xml + pytest -v --color=yes --cov=brainglobe-stitch --cov-report=xml """ diff --git a/tests/test_unit/conftest.py b/tests/test_unit/conftest.py new file mode 100644 index 0000000..cf87a31 --- /dev/null +++ b/tests/test_unit/conftest.py @@ -0,0 +1,63 @@ +import shutil +from pathlib import Path + +import pooch +import pytest + +TEMP_DIR = Path.home() / "temp_test_directory" +TEST_DATA_URL = "https://gin.g-node.org/IgorTatarnikov/brainglobe-stitch-test/raw/master/brainglobe-stitch/brainglobe-stitch-test-data.zip" + + +@pytest.fixture(scope="session", autouse=True) +def download_test_data(): + TEMP_DIR.mkdir(exist_ok=True) + pooch.retrieve( + TEST_DATA_URL, + known_hash="9437cb05566f03cd78cd905da1cd939dd1bb439837f96e32c6fb818b9b010a6e", + processor=pooch.Unzip(extract_dir=str(TEMP_DIR)), + ) + + yield TEMP_DIR + + shutil.rmtree(TEMP_DIR) + + +@pytest.fixture(scope="session") +def test_data_directory(): + yield TEMP_DIR + + +@pytest.fixture(scope="module") +def naive_bdv_directory(): + test_dir = Path.home() / "test_directory" + + shutil.copytree( + TEMP_DIR, + test_dir, + dirs_exist_ok=True, + ignore=shutil.ignore_patterns("*_interpolated_bdv.h5"), + ) + # Create UNIX style hidden files that should be ignored + (test_dir / ".test_data_bdv.h5").touch() + (test_dir / ".test_data_bdv.h5_meta.txt").touch() + (test_dir / ".test_data_bdv.h5_meta.txt").touch() + + yield test_dir + + shutil.rmtree(test_dir) + + +@pytest.fixture +def bdv_directory_function_level(): + test_dir = Path.home() / "quick_test_directory" + + shutil.copytree( + TEMP_DIR, + test_dir, + dirs_exist_ok=True, + ignore=shutil.ignore_patterns("*_interpolated_bdv.h5"), + ) + + yield test_dir + + shutil.rmtree(test_dir) diff --git a/tests/test_unit/test_image_mosaic.py b/tests/test_unit/test_image_mosaic.py new file mode 100644 index 0000000..fabef43 --- /dev/null +++ b/tests/test_unit/test_image_mosaic.py @@ -0,0 +1,99 @@ +import os + +import pytest + +from brainglobe_stitch.image_mosaic import ImageMosaic + +# The tiles lie in one z-plane and are arranged in a 2x2 grid with 2 channels. +# Each tile is 128x128x107 pixels (x, y, z). +# The tiles overlap by 10% in x and y (13 pixels). +# The tiles are arranged in the following pattern: +# channel 0 | channel 1 +# 00 10 | 04 05 +# 01 11 | 14 15 +NUM_TILES = 8 +NUM_RESOLUTIONS = 5 +NUM_CHANNELS = 2 +TILE_SIZE = (107, 128, 128) +# Expected tile config for the test data in test_data_bdv.h5 +# The tile positions are in pixels in x, y, z order +EXPECTED_TILE_CONFIG = [ + "dim=3", + "00;;(0,0,0)", + "01;;(0,115,0)", + "04;;(0,0,0)", + "05;;(0,115,0)", + "10;;(115,0,0)", + "11;;(115,115,0)", + "14;;(115,0,0)", + "15;;(115,115,0)", +] +# Expected tile positions for the test data in test_data_bdv.h5 +# The tile positions are in pixels in z, y, x order +EXPECTED_TILE_POSITIONS = [ + [0, 0, 0], + [0, 115, 0], + [0, 0, 0], + [0, 115, 0], + [0, 0, 115], + [0, 115, 115], + [0, 0, 115], + [0, 115, 115], +] + + +@pytest.fixture(scope="module") +def image_mosaic(naive_bdv_directory): + os.remove(naive_bdv_directory / "test_data_bdv_tile_config.txt") + image_mosaic = ImageMosaic(naive_bdv_directory) + + yield image_mosaic + + # Explicit call to clean up open h5 files + image_mosaic.__del__() + + +def test_image_mosaic_init(image_mosaic, naive_bdv_directory): + image_mosaic = image_mosaic + assert image_mosaic.xml_path == naive_bdv_directory / "test_data_bdv.xml" + assert ( + image_mosaic.meta_path + == naive_bdv_directory / "test_data_bdv.h5_meta.txt" + ) + assert image_mosaic.h5_path == naive_bdv_directory / "test_data_bdv.h5" + assert ( + image_mosaic.meta_path + == naive_bdv_directory / "test_data_bdv.h5_meta.txt" + ) + assert image_mosaic.h5_file is not None + assert len(image_mosaic.channel_names) == NUM_CHANNELS + assert len(image_mosaic.tiles) == NUM_TILES + assert len(image_mosaic.tile_names) == NUM_TILES + assert image_mosaic.x_y_resolution == 4.08 + assert image_mosaic.z_resolution == 5.0 + assert image_mosaic.num_channels == NUM_CHANNELS + + +def test_write_big_stitcher_tile_config(image_mosaic, naive_bdv_directory): + if (naive_bdv_directory / "test_data_bdv_tile_config.txt").exists(): + os.remove(naive_bdv_directory / "test_data_bdv_tile_config.txt") + + image_mosaic.write_big_stitcher_tile_config( + naive_bdv_directory / "test_data_bdv.h5_meta.txt" + ) + + assert (naive_bdv_directory / "test_data_bdv_tile_config.txt").exists() + + with open(naive_bdv_directory / "test_data_bdv_tile_config.txt", "r") as f: + for idx, line in enumerate(f): + assert line.strip() == EXPECTED_TILE_CONFIG[idx] + + +def test_data_for_napari(image_mosaic): + data = image_mosaic.data_for_napari(0) + + assert len(data) == NUM_TILES + + for i in range(NUM_TILES): + assert data[i][0].shape == TILE_SIZE + assert (data[i][1] == EXPECTED_TILE_POSITIONS[i]).all() diff --git a/tests/test_unit/test_placeholder.py b/tests/test_unit/test_placeholder.py deleted file mode 100644 index 3ada1ee..0000000 --- a/tests/test_unit/test_placeholder.py +++ /dev/null @@ -1,2 +0,0 @@ -def test_placeholder(): - assert True diff --git a/tests/test_unit/test_stitching_widget.py b/tests/test_unit/test_stitching_widget.py new file mode 100644 index 0000000..980df0e --- /dev/null +++ b/tests/test_unit/test_stitching_widget.py @@ -0,0 +1,191 @@ +from pathlib import Path + +import dask.array as da +import napari.layers +import numpy as np +import pytest + +import brainglobe_stitch +from brainglobe_stitch.stitching_widget import ( + StitchingWidget, + add_tiles_from_mosaic, +) + + +def test_add_tiles_from_mosaic(): + num_tiles = 4 + + test_data = [] + for i in range(num_tiles): + test_data.append((da.ones((10, 10, 10)), np.array([i, i, i]))) + + tile_names = [f"s{i:02}" for i in range(num_tiles)] + + for data, tile in zip( + test_data, add_tiles_from_mosaic(test_data, tile_names) + ): + assert isinstance(tile, napari.layers.Image) + assert (tile.data == data[0]).all() + assert (tile.translate == data[1]).all() + + +def test_stitching_widget_init(make_napari_viewer_proxy): + viewer = make_napari_viewer_proxy() + stitching_widget = StitchingWidget(viewer) + + assert stitching_widget._viewer == viewer + assert stitching_widget.image_mosaic is None + assert len(stitching_widget.tile_layers) == 0 + assert stitching_widget.resolution_to_display == 3 + + +def test_on_open_file_dialog_clicked(make_napari_viewer_proxy, mocker): + viewer = make_napari_viewer_proxy() + stitching_widget = StitchingWidget(viewer) + test_dir = str(Path.home() / "test_dir") + mocker.patch( + "brainglobe_stitch.stitching_widget.QFileDialog.getExistingDirectory", + return_value=test_dir, + ) + mocker.patch( + "brainglobe_stitch.stitching_widget.StitchingWidget.check_and_load_mesospim_directory", + ) + + stitching_widget._on_open_file_dialog_clicked() + + assert stitching_widget.mesospim_directory_text_field.text() == test_dir + assert stitching_widget.working_directory == Path(test_dir) + + +def test_on_mesospim_directory_text_edited(make_napari_viewer_proxy, mocker): + viewer = make_napari_viewer_proxy() + stitching_widget = StitchingWidget(viewer) + test_dir = str(Path.home() / "test_dir") + mocker.patch( + "brainglobe_stitch.stitching_widget.StitchingWidget.check_and_load_mesospim_directory", + ) + + stitching_widget.mesospim_directory_text_field.setText(test_dir) + + stitching_widget._on_mesospim_directory_text_edited() + + assert stitching_widget.working_directory == Path(test_dir) + + +def test_on_create_pyramid_button_clicked(make_napari_viewer_proxy, mocker): + viewer = make_napari_viewer_proxy() + stitching_widget = StitchingWidget(viewer) + stitching_widget.h5_path = Path.home() / "test_path" + mock_create_worker = mocker.patch( + "brainglobe_stitch.stitching_widget.create_worker", + autospec=True, + ) + + stitching_widget._on_create_pyramid_button_clicked() + + mock_create_worker.assert_called_once_with( + brainglobe_stitch.file_utils.create_pyramid_bdv_h5, + stitching_widget.h5_path, + yield_progress=True, + ) + + assert not stitching_widget.create_pyramid_button.isEnabled() + assert stitching_widget.add_tiles_button.isEnabled() + + +def test_on_add_tiles_button_clicked( + make_napari_viewer_proxy, naive_bdv_directory, mocker +): + viewer = make_napari_viewer_proxy() + stitching_widget = StitchingWidget(viewer) + stitching_widget.working_directory = naive_bdv_directory + + mock_create_worker = mocker.patch( + "brainglobe_stitch.stitching_widget.create_worker", + autospec=True, + ) + + stitching_widget._on_add_tiles_button_clicked() + + mock_create_worker.assert_called_once() + + +@pytest.mark.parametrize("num_layers", [1, 2, 5]) +def test_set_tile_layers_multiple(make_napari_viewer_proxy, num_layers): + viewer = make_napari_viewer_proxy() + stitching_widget = StitchingWidget(viewer) + test_data = da.ones((10, 10, 10)) + + test_layers = [] + + for i in range(num_layers): + test_layer = napari.layers.Image(data=test_data) + stitching_widget._set_tile_layers(test_layer) + test_layers.append(test_layer) + + assert len(stitching_widget.tile_layers) == num_layers + for i in range(num_layers): + assert stitching_widget.tile_layers[i] == test_layers[i] + + +def test_check_and_load_mesospim_directory( + make_napari_viewer_proxy, naive_bdv_directory +): + viewer = make_napari_viewer_proxy() + stitching_widget = StitchingWidget(viewer) + stitching_widget.working_directory = naive_bdv_directory + + stitching_widget.check_and_load_mesospim_directory() + + assert stitching_widget.h5_path == naive_bdv_directory / "test_data_bdv.h5" + assert ( + stitching_widget.xml_path == naive_bdv_directory / "test_data_bdv.xml" + ) + assert ( + stitching_widget.meta_path + == naive_bdv_directory / "test_data_bdv.h5_meta.txt" + ) + assert stitching_widget.add_tiles_button.isEnabled() + + +def test_check_and_load_mesospim_directory_no_pyramid( + make_napari_viewer_proxy, bdv_directory_function_level, mocker +): + viewer = make_napari_viewer_proxy() + stitching_widget = StitchingWidget(viewer) + stitching_widget.working_directory = bdv_directory_function_level + + mock_show_warning = mocker.patch( + "brainglobe_stitch.stitching_widget.show_warning" + ) + + stitching_widget.check_and_load_mesospim_directory() + + mock_show_warning.assert_called_once_with("Resolution pyramid not found") + assert not stitching_widget.add_tiles_button.isEnabled() + assert stitching_widget.create_pyramid_button.isEnabled() + + +@pytest.mark.parametrize( + "file_to_remove", + ["test_data_bdv.h5", "test_data_bdv.xml", "test_data_bdv.h5_meta.txt"], +) +def test_check_and_load_mesospim_directory_missing_files( + make_napari_viewer_proxy, + bdv_directory_function_level, + mocker, + file_to_remove, +): + viewer = make_napari_viewer_proxy() + stitching_widget = StitchingWidget(viewer) + stitching_widget.working_directory = bdv_directory_function_level + error_message = "mesoSPIM directory not valid" + + mock_show_warning = mocker.patch( + "brainglobe_stitch.stitching_widget.show_warning" + ) + (bdv_directory_function_level / file_to_remove).unlink() + + stitching_widget.check_and_load_mesospim_directory() + + mock_show_warning.assert_called_once_with(error_message) diff --git a/tests/test_unit/test_tile.py b/tests/test_unit/test_tile.py new file mode 100644 index 0000000..eb9ed2d --- /dev/null +++ b/tests/test_unit/test_tile.py @@ -0,0 +1,29 @@ +from brainglobe_stitch.tile import Tile + + +def test_tile_init(): + tile_id = 0 + attribute_tile_id = 0 + channel_id = 0 + illumination_id = 0 + angle = 234.5 + + attributes = { + "channel": channel_id, + "tile": attribute_tile_id, + "illumination": illumination_id, + "angle": angle, + } + + tile = Tile("test", tile_id, attributes) + + assert tile.name == "test" + assert tile.tile_id == tile_id + assert tile.position == [0, 0, 0] + assert len(tile.neighbours) == 0 + assert len(tile.data_pyramid) == 0 + assert len(tile.resolution_pyramid) == 0 + assert tile.channel_id == channel_id + assert tile.channel_name is None + assert tile.illumination_id == illumination_id + assert tile.angle == angle