From 6711b2b9694308678cbc578510932694d22df4c3 Mon Sep 17 00:00:00 2001 From: IgorTatarnikov Date: Thu, 11 Jan 2024 18:16:52 +0000 Subject: [PATCH] Using Zarr writer to write direct to disk, write ome zarr metadata after --- mesospim_stitcher/fuse.py | 117 ++++++++++++++++++++++++++++++++-- mesospim_stitcher/workflow.py | 50 ++++++--------- 2 files changed, 131 insertions(+), 36 deletions(-) diff --git a/mesospim_stitcher/fuse.py b/mesospim_stitcher/fuse.py index e7e2c32..49c167c 100644 --- a/mesospim_stitcher/fuse.py +++ b/mesospim_stitcher/fuse.py @@ -8,8 +8,9 @@ import numpy as np import zarr from numcodecs import Blosc +from ome_zarr.dask_utils import downscale_nearest from ome_zarr.io import parse_url -from ome_zarr.writer import write_image +from ome_zarr.writer import write_image, write_multiscales_metadata from mesospim_stitcher.file_utils import write_bdv_xml @@ -106,6 +107,113 @@ def fuse_image( max(translation[1] for translation in translations), ) + fuse_to_bdv_h5( + fused_image_shape, + group, + output_path, + tile_names, + translations, + xml_path, + ) + fuse_to_zarr( + fused_image_shape, group, output_path, tile_names, translations + ) + + input_file.close() + + +def fuse_to_zarr( + fused_image_shape: tuple[int, int, int], + group: h5py.Group, + output_path: Path, + tile_names: list, + translations: list, +): + num_tiles = len(tile_names) + + tiles = [] + + for child in group: + curr_tile = group[f"{child}/0/cells"] + tiles.append(da.from_array(curr_tile)) + + # fused_image = da.zeros(fused_image_shape, dtype="int16") + store = zarr.NestedDirectoryStore(str(output_path)) + root = zarr.group(store=store) + fused_image_store = root.create( + "0", + shape=fused_image_shape, + chunks=(4, 2048, 2048), + dtype="i2", + compressor=Blosc(cname="zstd", clevel=4, shuffle=Blosc.SHUFFLE), + ) + + coordinate_transformations = [ + [{"type": "scale", "scale": [5.0, 2.6, 2.6]}], + [{"type": "scale", "scale": [5.0, 5.2, 5.2]}], + [{"type": "scale", "scale": [5.0, 10.4, 10.4]}], + [{"type": "scale", "scale": [5.0, 20.8, 20.8]}], + [{"type": "scale", "scale": [5.0, 41.6, 41.6]}], + ] + + for i in range(num_tiles - 1, -1, -1): + x_s, x_e, y_s, y_e, z_s, z_e = translations[i] + curr_tile = tiles[i] + fused_image_store[z_s:z_e, y_s:y_e, x_s:x_e] = curr_tile + + print("Done tile " + str(i)) + + for i in range(1, len(coordinate_transformations)): + prev_resolution = da.from_zarr(root[f"{i - 1}"]) + downsampled_image = downscale_nearest(prev_resolution, (1, 2, 2)) + downsampled_shape = downsampled_image.shape + downsampled_store = root.require_dataset( + f"{i}", + shape=downsampled_shape, + chunks=(4, (2048 // 2**i), (2048 // 2**i)), + dtype="i2", + compressor=Blosc(cname="zstd", clevel=4, shuffle=Blosc.SHUFFLE), + ) + downsampled_image.to_zarr(downsampled_store) + + datasets = [] + + for i, transform in enumerate(coordinate_transformations): + datasets.append( + {"path": f"{i}", "coordinateTransformations": transform} + ) + + write_multiscales_metadata( + group=root, + datasets=datasets, + axes=[ + {"name": "z", "type": "space", "unit": "micrometer"}, + {"name": "y", "type": "space", "unit": "micrometer"}, + {"name": "x", "type": "space", "unit": "micrometer"}, + ], + ) + + root.attrs["omero"] = { + "channels": [ + { + "window": {"start": 0, "end": 1200, "min": 0, "max": 65535}, + "label": "stitched", + "active": True, + } + ] + } + + # write_ome_zarr(output_path, fused_image, overwrite=True) + + +def fuse_to_bdv_h5( + fused_image_shape: tuple[int, int, int], + group: h5py.Group, + output_path: Path, + tile_names: list, + translations: list, + xml_path: Path, +): num_tiles = len(tile_names) output_file = h5py.File(output_path, mode="w") @@ -114,7 +222,6 @@ def fuse_image( shape=fused_image_shape, dtype="i2", ) - subdivisions = np.array( [[32, 32, 16], [32, 32, 16], [32, 32, 16], [32, 32, 16]], dtype=np.int16, @@ -123,7 +230,6 @@ def fuse_image( [[1, 1, 1], [2, 2, 1], [4, 4, 1], [8, 8, 1]], dtype=np.int16, ) - output_file.require_dataset( "s00/resolutions", data=resolutions, @@ -137,7 +243,6 @@ def fuse_image( shape=subdivisions.shape, ) ds_list = [ds] - for i in range(1, resolutions.shape[0]): ds_list.append( output_file.require_dataset( @@ -172,9 +277,7 @@ def fuse_image( ] write_bdv_xml(Path("testing.xml"), xml_path, output_path, ds.shape) - output_file.close() - input_file.close() def write_ome_zarr(output_path: Path, image: da, overwrite: bool): @@ -205,7 +308,7 @@ def write_ome_zarr(output_path: Path, image: da, overwrite: bool): [{"type": "scale", "scale": [10.0, 41.6, 41.6]}], ], storage_options=dict( - chunks=(4, image.shape[1], image.shape[2]), + chunks=(2, image.shape[1], image.shape[2]), compressor=compressor, ), ) diff --git a/mesospim_stitcher/workflow.py b/mesospim_stitcher/workflow.py index 13b07e0..8b5ecec 100644 --- a/mesospim_stitcher/workflow.py +++ b/mesospim_stitcher/workflow.py @@ -2,28 +2,22 @@ import numpy as np -from mesospim_stitcher.big_stitcher_bridge import run_big_stitcher -from mesospim_stitcher.file_utils import ( - create_pyramid_bdv_h5, - write_big_stitcher_tile_config, -) from mesospim_stitcher.fuse import fuse_image XML_PATH = ( - "D:/TiledDataset/2.5x_tile/" - "2.5x_tile_igor_rightonly_Mag2.5x_ch488_ch561_ch647_bdv3.xml" + "C:/Users/Igor/Documents/NIU-dev/stitching/One_Channel/" + "2.5x_tile_igor_rightonly_Mag2.5x_ch488_ch561_ch647_bdv.xml" ) -H5_PATH = ( - "D:/TiledDataset/2.5x_tile/" - "2.5x_tile_igor_rightonly_Mag2.5x_ch488_ch561_ch647_bdv.h5" +H5_PATH = "C:/Users/Igor/Documents/NIU-dev/stitching/One_Channel/test.h5" +OUT_PATH = ( + "C:/Users/Igor/Documents/NIU-dev/stitching/One_Channel/test_out.zarr" ) META_PATH = ( - "D:/TiledDataset/2.5x_tile/" + "C:/Users/Igor/Documents/NIU-dev/stitching/One_Channel/" "2.5x_tile_igor_rightonly_Mag2.5x_ch488_ch561_ch647_bdv.h5_meta.txt" ) TILE_CONFIG_PATH = ( - "D:/TiledDataset/" - "2.5x_tile/" + "C:/Users/Igor/Documents/NIU-dev/stitching/One_Channel/" "2.5x_tile_igor_rightonly_Mag2.5x" "_ch488_ch561_ch647_bdv_tile_config.txt" ) @@ -35,24 +29,22 @@ [[32, 32, 16], [32, 32, 16], [32, 32, 16], [32, 32, 16], [32, 32, 16]] ) -create_pyramid_bdv_h5(Path(H5_PATH), DOWNSAMPLE_ARRAY, SUBDIVISION_ARRAY) - -tile_metadata = write_big_stitcher_tile_config(Path(TILE_CONFIG_PATH)) - -result_big_stitcher = run_big_stitcher( - Path("C:/Users/Igor/Documents/Fiji.app/ImageJ-win64.exe"), - Path(XML_PATH), - Path(TILE_CONFIG_PATH), -) - -with open("big_stitcher_output.txt", "w") as f: - f.write(result_big_stitcher.stdout) - f.write(result_big_stitcher.stderr) +# create_pyramid_bdv_h5(Path(H5_PATH), DOWNSAMPLE_ARRAY, SUBDIVISION_ARRAY) +# +# tile_metadata = write_big_stitcher_tile_config(Path(TILE_CONFIG_PATH)) +# +# result_big_stitcher = run_big_stitcher( +# Path("C:/Users/Igor/Documents/Fiji.app/ImageJ-win64.exe"), +# Path(XML_PATH), +# Path(TILE_CONFIG_PATH), +# ) +# +# with open("big_stitcher_output.txt", "w") as f: +# f.write(result_big_stitcher.stdout) +# f.write(result_big_stitcher.stderr) fuse_image( Path(XML_PATH), Path(H5_PATH), - Path( - "D:/TiledDataset/2.5x_tile/2.5x_tile_igor_rightonly_Mag2.5x_ch488_ch561_ch647_bdv_fused.h5" - ), + Path(OUT_PATH), )