add wingdisc functionality

brainglobe_template_builder/preproc/load_wingdisc.py

@@ -0,0 +1,62 @@
+from pathlib import Path
+
+import bioio_czi
+import bioio_sldy
+import numpy as np
+import tifffile as tiff
+from bioio import BioImage
+
+
+def load_images(path: Path, file_type: str):
+    if file_type == "czi":
+        array = BioImage(path, reader=bioio_czi.Reader)
+        channel_array = array.get_image_data(
+            "ZYX", T=0
+        )  # will only read channel 1 into memory
+    elif file_type == "sldy":
+        array = BioImage(path, reader=bioio_sldy.Reader)
+        channel_array = array.get_image_data("ZYX", T=0)
+    elif file_type == "tif":
+        channel_array = tiff.imread(path)
+    else:
+        print("File format not supported!")
+    return channel_array
+
+
+def load_channel_data(
+    image_path: Path, image_format: str, channel: int
+) -> np.ndarray:
+    """load the images for certain formats
+
+    Args:
+        image_path: path to Path
+        image_format: image format ('czi' or 'sldy', can add more simplely
+        channel: the image channel that need to be loaded
+
+    Returns:
+        np.ndarray:
+
+    Raises:
+        ValueError:
+    """
+    # make a dictionary
+    readers = {"czi": bioio_czi.Reader, "sldy": bioio_sldy.Reader}
+
+    # verify if the specified image format is supported
+    if image_format not in readers:
+        raise ValueError(
+            f"Unsupported format: {image_format}. Supported: {list(readers.keys())}"
+        )
+
+    # load images
+    bio_image = BioImage(image_path, reader=readers[image_format])
+    return bio_image.get_image_data("ZYX", T=0, C=channel)
+
+
+"""
+# Check the image dimensions
+original_file_path = Path('D:/UCL/Postgraduate_programme/wing_disc_imaging_data/wd2_overview-Airyscan-Processing-02/wd2_overview-Airyscan-Processing-02.czi')
+assert original_file_path.exists()
+czi_array = load_channel_data(original_file_path,image_format='czi',channel=1)
+print(czi_array.shape)
+"""

brainglobe_template_builder/preproc/mirroring_wingdisc.py

@@ -0,0 +1,20 @@
+import numpy as np
+
+
+def mirroring(
+    image_array: np.ndarray,
+) -> np.ndarray:
+    """Create a Mirror image of a 3D image along the (1,0) plane.
+
+    Parameters
+    ----------
+    image_array : np.ndarray
+        Array for image to mirror.
+    """
+    z, x, y = image_array.shape
+    mirrored_image = np.zeros((z, x, y))
+    for k in range(0, z):
+        for i in range(0, x):
+            for j in range(0, y):
+                mirrored_image[k, -(i + 1), j] = image_array[k, i, j]
+    return mirrored_image

examples/drosophila/preprocessing_wing-disc.py

@@ -0,0 +1,221 @@
+import argparse
+from pathlib import Path
+
+import ants
+import numpy as np
+import pandas as pd
+from brainglobe_utils.IO.image import load_any, save_any
+from dask import array as da
+from loguru import logger
+
+from brainglobe_template_builder.io import (
+    file_path_with_suffix,
+    save_as_asr_nii,
+)
+from brainglobe_template_builder.preproc.masking import create_mask
+from brainglobe_template_builder.preproc.transform_utils import (
+    downsample_anisotropic_image_stack,
+)
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(description="Download source image")
+    parser.add_argument(
+        "--source_data_root",
+        type=str,
+        help="Path to the source data folder. The source data should contain"
+        "a subfolder per subject, with picture files within it",
+        required=True,
+    )
+    parser.add_argument(
+        "--template_building_root",
+        type=str,
+        help="Path to the template-building root folder.Results will be "
+        "written to the rawdata folder.",
+        required=True,
+    )
+    parser.add_argument(
+        "--target_isotropic_resolution",
+        type=int,
+        help="Target isotropic resolution",
+        required=True,
+    )
+
+    parser.add_argument(
+        "--data_catalog",
+        type=str,
+        help="The full path to the data catalog file",
+        required=True,
+    )
+
+    parser.add_argument(
+        "--dataset",
+        type=str,
+        help="The dataset used to process",
+        required=True,
+    )
+
+    args = parser.parse_args()
+
+    source_data = Path(args.source_data_root)
+    template_building_root = Path(args.template_building_root)
+    target_isotropic_resolution = int(args.target_isotropic_resolution)
+
+    in_plane_resolution = 0.55
+    out_of_plane_resolution = 1
+
+    in_plane_factor = int(
+        np.ceil(target_isotropic_resolution / in_plane_resolution)
+    )
+    axial_factor = int(
+        np.ceil(target_isotropic_resolution / out_of_plane_resolution)
+    )
+
+    template_raw_data = template_building_root / "rawdata"
+    template_raw_data.mkdir(exist_ok=True, parents=True)
+
+    # Load the data catalog from argument
+    data_catalog_path = Path(args.data_catalog)
+    data_catalog = pd.read_csv(data_catalog_path)
+    dataset = args.dataset
+
+    # Specified the dataset to process
+    dataset_catalog = data_catalog[data_catalog["dataset"] == dataset]
+    logger.debug(f"Loaded {dataset} dataset catalog from {data_catalog_path}.")
+
+    # Check if there are right wing discs is in the dataset
+    right_wingdisc_catalog = dataset_catalog[dataset_catalog["is_left"] == "n"]
+    if right_wingdisc_catalog.empty:
+        logger.info(f"No right wing discs found in {dataset} dataset.")
+
+    for sample_folder in source_data.iterdir():
+        # Load the images and specify the filename of processed images
+        logger.info(f"Downsampling {sample_folder}...")
+        sample_id = str(sample_folder.name).split("_")[0].lower()
+        channel = "membrane"
+        downsampled_filename = (
+            f"{sample_id}_res-{target_isotropic_resolution}"
+            f"um_channel-{channel}.tif"
+        )
+        assert Path(sample_folder).exists(), f"{sample_folder} not found"
+        original_file_path = Path(sample_folder) / f"{sample_folder.name}.tif"
+        assert Path(
+            original_file_path
+        ).exists(), f"Filepath {original_file_path} not found"
+        image_array = load_any(original_file_path)
+
+        # Do mirroring if the sample is right wing disc
+        if (
+            str(sample_folder.name)
+            in right_wingdisc_catalog["filename"].astype(str).tolist()
+        ):
+            image_array = np.flip(image_array, axis=2)
+            logger.info(f"Mirrored {sample_folder.name}.")
+
+        # Downsample the image array
+        image_dask = da.from_array(image_array, chunks={0: 1, 1: -1, 2: -1})
+        down_sampled_image = downsample_anisotropic_image_stack(
+            image_dask, in_plane_factor, axial_factor
+        )
+        down_sampled_image = down_sampled_image.astype(np.uint16)
+
+        # Save the downsampled image as tif
+        saving_folder = (
+            template_raw_data
+            / f"{source_data.name}"
+            / downsampled_filename.split(".")[0]
+        )
+        saving_folder.mkdir(exist_ok=True, parents=True)
+        assert Path(
+            saving_folder
+        ).exists(), f"Filepath {saving_folder} not found"
+        saving_path = saving_folder / downsampled_filename
+        save_any(down_sampled_image, saving_path)
+        logger.info(
+            f"{sample_folder} downsampled, saved as {downsampled_filename}"
+        )
+
+        # Save the downsampled image as nifti
+        nii_path = file_path_with_suffix(
+            saving_path, "_downsampled", new_ext=".nii.gz"
+        )
+        vox_sizes = [
+            target_isotropic_resolution,
+        ] * 3
+        save_as_asr_nii(down_sampled_image, vox_sizes, nii_path)
+        logger.info(f"Saved downsampled image as {nii_path.name}.")
+
+        # Generate the wingdisc mask
+        image_ants = ants.image_read(nii_path.as_posix())
+        mask_data = create_mask(
+            image_ants.numpy(),
+            gauss_sigma=10,
+            threshold_method="triangle",
+            closing_size=5,
+        )
+        mask_path = file_path_with_suffix(nii_path, "_mask")
+        mask = image_ants.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 = (
+            saving_folder / f"{sample_id}_downsampled_mask-overlay.png"
+        )
+        ants.plot(
+            image_ants,
+            mask,
+            overlay_alpha=0.5,
+            axis=1,
+            title="Wingdisc mask over image",
+            filename=mask_plot_path.as_posix(),
+        )
+        logger.debug("Plotted overlay to visually check mask.")
+
+        # Write a text file to record the file path to the downsampled image and mask created
+        with open("brain_path.txt", "a") as f:
+            f.write(f"{nii_path} + \n")
+        with open("mask_path.txt", "a") as f:
+            f.write(f"{mask_path} + \n")
+        logger.debug(
+            "Recorded the file path to the brain_path.txt and mask_path.txt created."
+        )
+
+        # Write a text file to record the file path to the downsampled image and mask created
+        nii_path_str = str(nii_path)
+        downsampled_txt_file_name = "brain_paths.txt"
+        downsampled_txt_file_path = (
+            template_raw_data
+            / f"{source_data.name}"
+            / downsampled_txt_file_name
+        )
+        mask_path_str = str(mask_path)
+        mask_txt_file_name = "mask_paths.txt"
+        mask_txt_file_path = (
+            template_raw_data / f"{source_data.name}" / mask_txt_file_name
+        )
+        # Read the file(if exist first) to check if the path already exists
+        try:
+            with open(downsampled_txt_file_path, "r", encoding="utf-8") as f:
+                lines_downsampled = (
+                    f.read().splitlines()
+                )  # Read lines without newlines
+            with open(mask_txt_file_path, "r", encoding="utf-8") as f:
+                lines_mask = f.read().splitlines()
+        except FileNotFoundError:
+            lines_downsampled = []
+            lines_mask = []  # If file doesn't exist, start with an empty list
+        if nii_path_str not in lines_downsampled:
+            with open(downsampled_txt_file_path, "a", encoding="utf-8") as f:
+                f.write(f"{nii_path_str}\n")
+                logger.debug(
+                    f"Recorded {nii_path_str} in {downsampled_txt_file_path}."
+                )
+        if mask_path_str not in lines_mask:
+            with open(mask_txt_file_path, "a", encoding="utf-8") as f:
+                f.write(f"{mask_path_str}\n")
+                logger.debug(
+                    f"Recorded {mask_path_str} in {mask_txt_file_path}."
+                )

pyproject.toml

@@ -48,6 +48,11 @@ dev = [
     "ruff",
     "setuptools_scm",
 ]
+wingdisc = [
+    'bioio',
+    'bioio_czi',
+    'bioio_sldy'
+]
 
 [build-system]
 requires = ["setuptools>=45", "wheel", "setuptools_scm[toml]>=6.2"]

tests/test_unit/test_load_wingdisc.py

@@ -0,0 +1,19 @@
+from pathlib import Path
+
+import tifffile as tiff
+
+from brainglobe_template_builder.preproc.load_wingdisc import load_channel_data
+
+
+def test_load_images_czi():
+    path_czi = Path(
+        "D:/UCL/Postgraduate_programme/wing_disc_imaging_data/wd2_overview-Airyscan-Processing-02/wd2_overview-Airyscan-Processing-02.czi"
+    )
+    assert path_czi.exists()
+    path_tiff = Path(
+        "D:/UCL/Postgraduate_programme/wing_disc_imaging_data/wd2_overview-Airyscan-Processing-02/wd2_overview-Airyscan-Processing-02_C1.tif"
+    )
+    assert path_tiff.exists()
+    czi_array = load_channel_data(path_czi, image_format="czi", channel=1)
+    tiff_array = tiff.imread(path_tiff)
+    assert (czi_array == tiff_array).all()

tests/test_unit/test_mirroring_wingdisc.py

@@ -0,0 +1,10 @@
+import numpy as np
+
+from brainglobe_template_builder.preproc.mirroring_wingdisc import mirroring
+
+
+def test_mirroring():
+    shape = np.random.randint(1, 10, size=3)
+    random_image_array = np.random.randint(0, 100, size=shape)
+    mirrored_image = mirroring(random_image_array)
+    assert (mirrored_image == np.flip(random_image_array, 1)).all()