diff --git a/fractal/tasks/create_zarr_structure.py b/fractal/tasks/create_zarr_structure.py index 49c268a1..60d769c0 100644 --- a/fractal/tasks/create_zarr_structure.py +++ b/fractal/tasks/create_zarr_structure.py @@ -186,9 +186,6 @@ def create_zarr_structure( site_metadata, total_files = parse_yokogawa_metadata( mrf_path, mlf_path ) - # FIXME: hardcoded - image_size = {"x": 2560, "y": 2160} - well_size_z = 10 # Extract pixel sizes pixel_size_z = site_metadata["pixel_size_z"][0] @@ -348,8 +345,6 @@ def create_zarr_structure( # Prepare and write anndata table of FOV ROIs FOV_ROIs_table = prepare_FOV_ROI_table( site_metadata.loc[f"{row+column}"], - image_size=image_size, - well_size_z=well_size_z, ) group_tables = group_FOV.create_group("tables/") # noqa: F841 write_elem(group_tables, "FOV_ROI_table", FOV_ROIs_table) diff --git a/fractal/tasks/illumination_correction.py b/fractal/tasks/illumination_correction.py index 2f5fdb01..fdd4ea3b 100644 --- a/fractal/tasks/illumination_correction.py +++ b/fractal/tasks/illumination_correction.py @@ -22,6 +22,9 @@ from fractal.tasks.lib_pyramid_creation import write_pyramid from fractal.tasks.lib_regions_of_interest import convert_ROI_table_to_indices +from fractal.tasks.lib_regions_of_interest import ( + extract_zyx_pixel_sizes_from_zattrs, +) from fractal.tasks.lib_regions_of_interest import ( split_3D_indices_into_z_layers, ) @@ -135,6 +138,7 @@ def illumination_correction( newzarrurl += "/" # Hard-coded values for the image size + # FIXME can we parse this from somewhere? img_size_y = 2160 img_size_x = 2560 @@ -192,9 +196,15 @@ def illumination_correction( # Read FOV ROIs FOV_ROI_table = ad.read_zarr(f"{zarrurl}tables/FOV_ROI_table") + # Read pixel sizes from zattrs file + pixel_sizes_zyx = extract_zyx_pixel_sizes_from_zattrs(zarrurl + ".zattrs") + # Create list of indices for 3D FOVs spanning the entire Z direction list_indices = convert_ROI_table_to_indices( - FOV_ROI_table, level=0, coarsening_xy=coarsening_xy + FOV_ROI_table, + level=0, + coarsening_xy=coarsening_xy, + pixel_sizes_zyx=pixel_sizes_zyx, ) # Create the final list of single-Z-layer FOVs @@ -216,8 +226,11 @@ def illumination_correction( data_zyx_new = da.empty_like(data_zyx) for indices in list_indices: + s_z, e_z, s_y, e_y, s_x, e_x = indices[:] shape = [e_z - s_z, e_y - s_y, e_x - s_x] + if min(shape) == 0: + raise Exception(f"ERROR: ROI indices lead to shape {shape}") new_img = delayed_correct( data_zyx[s_z:e_z, s_y:e_y, s_x:e_x], illum_img, diff --git a/fractal/tasks/image_labeling.py b/fractal/tasks/image_labeling.py index dc47d2d7..265113e8 100644 --- a/fractal/tasks/image_labeling.py +++ b/fractal/tasks/image_labeling.py @@ -17,7 +17,6 @@ import time from concurrent.futures import ThreadPoolExecutor -import anndata as ad import dask import dask.array as da import numpy as np @@ -26,6 +25,9 @@ from cellpose import models from fractal.tasks.lib_pyramid_creation import write_pyramid +from fractal.tasks.lib_regions_of_interest import ( + extract_zyx_pixel_sizes_from_zattrs, +) def segment_FOV( @@ -127,11 +129,11 @@ def image_labeling( do_3D = data_zyx.shape[0] > 1 if do_3D: if anisotropy is None: - - adata = ad.read_zarr(f"{zarrurl}tables/FOV_ROI_table") - pixel_size_x = adata[:, "pixel_size_x"].X[0, 0] - pixel_size_y = adata[:, "pixel_size_y"].X[0, 0] - pixel_size_z = adata[:, "pixel_size_z"].X[0, 0] + # Read pixel sizes from zattrs file + pxl_zyx = extract_zyx_pixel_sizes_from_zattrs( + zarrurl + ".zattrs", level=labeling_level + ) + pixel_size_z, pixel_size_y, pixel_size_x = pxl_zyx[:] if not np.allclose(pixel_size_x, pixel_size_y): raise Exception( "ERROR: XY anisotropy detected\n" @@ -139,11 +141,6 @@ def image_labeling( f"pixel_size_y={pixel_size_y}" ) anisotropy = pixel_size_z / pixel_size_x - if labeling_level > 0: - raise NotImplementedError( - "TODO: fix automatic anisotropy " - "detection for higher levels" - ) # Check model_type if model_type not in ["nuclei", "cyto2", "cyto"]: diff --git a/fractal/tasks/lib_regions_of_interest.py b/fractal/tasks/lib_regions_of_interest.py index bb5e16b1..89de4b96 100644 --- a/fractal/tasks/lib_regions_of_interest.py +++ b/fractal/tasks/lib_regions_of_interest.py @@ -1,6 +1,7 @@ +import json import math -from typing import Dict from typing import List +from typing import Tuple from typing import Union import anndata as ad @@ -10,24 +11,21 @@ def prepare_FOV_ROI_table( df: pd.DataFrame, - image_size: Union[Dict, None] = None, - well_size_z: int = None, ) -> ad.AnnData: - if image_size is None: - raise Exception("Missing image_size arg in prepare_ROIs_table") - if well_size_z is None: - raise Exception("Missing well_size_z arg in prepare_ROIs_table") + for mu in ["x", "y", "z"]: - # Reset reference values for coordinates - df["x_micrometer"] -= df["x_micrometer"].min() - df["y_micrometer"] -= df["y_micrometer"].min() - df["z_micrometer"] -= df["z_micrometer"].min() + # Reset reference values for coordinates + df[f"{mu}_micrometer"] -= df[f"{mu}_micrometer"].min() - # Obtain box size in physical units - df["len_x_micrometer"] = image_size["x"] * df["pixel_size_x"] - df["len_y_micrometer"] = image_size["y"] * df["pixel_size_y"] - df["len_z_micrometer"] = well_size_z * df["pixel_size_z"] + # Obtain box size in physical units + df[f"len_{mu}_micrometer"] = df[f"{mu}_pixel"] * df[f"pixel_size_{mu}"] + + # Remove information about pixel sizes in physical units + df.drop(f"pixel_size_{mu}", inplace=True, axis=1) + + # Remove information about array size in pixels + df.drop(f"{mu}_pixel", inplace=True, axis=1) # Remove unused column df.drop("bit_depth", inplace=True, axis=1) @@ -52,22 +50,17 @@ def prepare_FOV_ROI_table( return adata -def convert_FOV_ROIs_3D_to_2D(adata: ad.AnnData = None) -> ad.AnnData: +def convert_FOV_ROIs_3D_to_2D( + adata: ad.AnnData = None, pixel_size_z: float = None +) -> ad.AnnData: - # FIXME: use this piece of code (or similar) in tests - """ - adata = ad.AnnData(X=np.zeros((2, 3)), dtype=int) - adata.obs_names = ["FOV1", "FOV2"] - adata.var_names = ["z", "len_z", "pixel_size_z"] - adata[:, "z"].X = [[2], [4]] - adata[:, "len_z"].X = [[5], [7]] - adata[:, "pixel_size_z"].X = [[2], [2]] - """ + if pixel_size_z is None: + raise Exception("Missing pixel_size_z in convert_FOV_ROIs_3D_to_2D") # Compress a 3D stack of images to a single Z plane, # with thickness equal to pixel_size_z df = adata.to_df() - df["len_z_micrometer"] = df["pixel_size_z"] + df["len_z_micrometer"] = pixel_size_z # Assign dtype explicitly, to avoid # >> UserWarning: X converted to numpy array with dtype float64 @@ -88,10 +81,13 @@ def convert_ROI_table_to_indices( ROI: ad.AnnData, level: int = 0, coarsening_xy: int = 2, -) -> List[List]: + pixel_sizes_zyx: Union[List[float], Tuple[float]] = None, +) -> List[List[int]]: list_indices = [] + pixel_size_z, pixel_size_y, pixel_size_x = pixel_sizes_zyx + for FOV in sorted(ROI.obs_names): # Extract data from anndata table @@ -101,9 +97,6 @@ def convert_ROI_table_to_indices( len_x_micrometer = ROI[FOV, "len_x_micrometer"].X[0, 0] len_y_micrometer = ROI[FOV, "len_y_micrometer"].X[0, 0] len_z_micrometer = ROI[FOV, "len_z_micrometer"].X[0, 0] - pixel_size_x = ROI[FOV, "pixel_size_x"].X[0, 0] - pixel_size_y = ROI[FOV, "pixel_size_y"].X[0, 0] - pixel_size_z = ROI[FOV, "pixel_size_z"].X[0, 0] # Set pyramid-level pixel sizes prefactor = coarsening_xy**level @@ -151,7 +144,10 @@ def split_3D_indices_into_z_layers( def _inspect_ROI_table( - path: str = None, level: int = 0, coarsening_xy: int = 2 + path: str = None, + level: int = 0, + coarsening_xy: int = 2, + pixel_sizes_zyx=[1.0, 0.1625, 0.1625], ) -> None: adata = ad.read_zarr(path) @@ -161,7 +157,10 @@ def _inspect_ROI_table( print() list_indices = convert_ROI_table_to_indices( - adata, level=level, coarsening_xy=coarsening_xy + adata, + level=level, + coarsening_xy=coarsening_xy, + pixel_sizes_zyx=pixel_sizes_zyx, ) list_indices = split_3D_indices_into_z_layers(list_indices) @@ -174,10 +173,120 @@ def _inspect_ROI_table( print() +def temporary_test(): + + pixel_size_z = 1.0 + pixel_size_y = 0.1625 + pixel_size_x = 0.1625 + + df = pd.DataFrame(np.zeros((2, 10)), dtype=int) + df.index = ["FOV1", "FOV2"] + df.columns = [ + "x_micrometer", + "y_micrometer", + "z_micrometer", + "x_pixel", + "y_pixel", + "z_pixel", + "pixel_size_x", + "pixel_size_y", + "pixel_size_z", + "bit_depth", + ] + df["x_micrometer"] = [0.0, 416.0] + df["y_micrometer"] = [0.0, 0.0] + df["z_micrometer"] = [0.0, 0.0] + df["x_pixel"] = [2560, 2560] + df["y_pixel"] = [2160, 2160] + df["z_pixel"] = [5, 5] + df["pixel_size_x"] = [pixel_size_x] * 2 + df["pixel_size_y"] = [pixel_size_y] * 2 + df["pixel_size_z"] = [pixel_size_z] * 2 + df["bit_depth"] = [16.0, 16.0] + + print("DataFrame") + print(df) + print() + + adata = prepare_FOV_ROI_table(df) + + print("AnnData table") + print(adata.var_names) + print(adata.obs_names) + print(adata.X) + print() + + print("Indices 3D") + pixel_sizes_zyx = [pixel_size_z, pixel_size_y, pixel_size_x] + list_indices = convert_ROI_table_to_indices( + adata, level=0, coarsening_xy=2, pixel_sizes_zyx=pixel_sizes_zyx + ) + for indices in list_indices: + print(indices) + print() + + print("Indices 3D / split") + list_indices = split_3D_indices_into_z_layers(list_indices) + for indices in list_indices: + print(indices) + print() + + print("Indices 2D") + adata = convert_FOV_ROIs_3D_to_2D(adata, pixel_size_z) + list_indices = convert_ROI_table_to_indices( + adata, level=0, coarsening_xy=2, pixel_sizes_zyx=pixel_sizes_zyx + ) + for indices in list_indices: + print(indices) + print() + + +def extract_zyx_pixel_sizes_from_zattrs(zattrs_path: str, level: int = 0): + with open(zattrs_path, "r") as jsonfile: + zattrs = json.load(jsonfile) + + try: + + # Identify multiscales + multiscales = zattrs["multiscales"] + + # Check that there is a single multiscale + if len(multiscales) > 1: + raise Exception(f"ERROR: There are {len(multiscales)} multiscales") + + # Check that there are no datasets-global transformations + if "coordinateTransformations" in multiscales[0].keys(): + raise Exception( + "ERROR: coordinateTransformations at the multiscales " + "level are not currently supported" + ) + + # Identify all datasets (AKA pyramid levels) + datasets = multiscales[0]["datasets"] + + # Select highest-resolution dataset + transformations = datasets[level]["coordinateTransformations"] + for t in transformations: + if t["type"] == "scale": + return t["scale"] + raise Exception( + "ERROR:" + f" no scale transformation found for level {level}" + f" in {zattrs_path}" + ) + + except KeyError as e: + raise KeyError( + "extract_zyx_pixel_sizes_from_zattrs failed, for {zattrs_path}\n", + e, + ) + + if __name__ == "__main__": - import sys + # import sys + # args = sys.argv[1:] + # types = [str, int, int] + # args = [types[ix](x) for ix, x in enumerate(args)] + # _inspect_ROI_table(*args) - args = sys.argv[1:] - types = [str, int, int] - args = [types[ix](x) for ix, x in enumerate(args)] - _inspect_ROI_table(*args) + temporary_test() diff --git a/fractal/tasks/replicate_zarr_structure_mip.py b/fractal/tasks/replicate_zarr_structure_mip.py index 2217d8b8..0a8ab2df 100644 --- a/fractal/tasks/replicate_zarr_structure_mip.py +++ b/fractal/tasks/replicate_zarr_structure_mip.py @@ -19,6 +19,9 @@ from anndata.experimental import write_elem from fractal.tasks.lib_regions_of_interest import convert_FOV_ROIs_3D_to_2D +from fractal.tasks.lib_regions_of_interest import ( + extract_zyx_pixel_sizes_from_zattrs, +) def replicate_zarr_structure_mip(zarrurl): @@ -113,7 +116,7 @@ def replicate_zarr_structure_mip(zarrurl): group_FOV = group_well.create_group(f"{FOV}/") # Copy .zattrs file at the COL/ROW/FOV level - path_FOV_zattrs = zarrurl + f"{row}/{column}/{FOV}/.zattrs" + path_FOV_zattrs = f"{zarrurl}{row}/{column}/{FOV}/.zattrs" with open(path_FOV_zattrs) as FOV_zattrs_file: FOV_zattrs = json.load(FOV_zattrs_file) for key in FOV_zattrs.keys(): @@ -121,11 +124,17 @@ def replicate_zarr_structure_mip(zarrurl): # Read FOV ROI table FOV_ROI_table = ad.read_zarr( - zarrurl + f"{row}/{column}/0/tables/FOV_ROI_table" + f"{zarrurl}{row}/{column}/0/tables/FOV_ROI_table" ) + # Read pixel sizes from zattrs file + pixel_sizes_zyx = extract_zyx_pixel_sizes_from_zattrs(path_FOV_zattrs) + pixel_size_z = pixel_sizes_zyx[0] + # Convert 3D FOVs to 2D - new_FOV_ROI_table = convert_FOV_ROIs_3D_to_2D(FOV_ROI_table) + new_FOV_ROI_table = convert_FOV_ROIs_3D_to_2D( + FOV_ROI_table, pixel_size_z + ) # Create table group and write new table group_tables = group_FOV.create_group("tables/")