|
| 1 | +# %% Import packages |
| 2 | + |
| 3 | +import itertools |
| 4 | +from pathlib import Path |
| 5 | +from bfio import BioReader |
| 6 | +import numpy as np |
| 7 | +from cloudvolume.lib import touch |
| 8 | +import json |
| 9 | +from cryoet_data_portal_neuroglancer.precompute.segmentation_mask import ( |
| 10 | + create_segmentation_chunk, |
| 11 | + write_metadata, |
| 12 | +) |
| 13 | +from cryoet_data_portal_neuroglancer.precompute.mesh import ( |
| 14 | + generate_multiresolution_mesh_from_segmentation, |
| 15 | +) |
| 16 | + |
| 17 | +# %% Define the path to the files |
| 18 | + |
| 19 | +HERE = Path(__file__).parent |
| 20 | + |
| 21 | +FILEPATH = Path( |
| 22 | + "..." |
| 23 | +) |
| 24 | +OUTPUT_PATH = Path( |
| 25 | + "..." |
| 26 | +) |
| 27 | +OUTPUT_PATH.mkdir(exist_ok=True, parents=True) |
| 28 | +OVERWRITE = False |
| 29 | + |
| 30 | +# %% Load the data |
| 31 | +br = BioReader(str(FILEPATH), backend="bioformats") |
| 32 | + |
| 33 | +# %% Inspect the data |
| 34 | +with open(OUTPUT_PATH / "metadata.txt", "w") as f: |
| 35 | + json.dump(br.metadata.model_dump_json(), f) |
| 36 | + f.write(br.metadata.model_dump_json()) |
| 37 | +print(br.shape) |
| 38 | + |
| 39 | +# %% Extract from the data - units are in nm |
| 40 | +size_x = br.metadata.images[0].pixels.physical_size_x |
| 41 | +if size_x is None: |
| 42 | + size_x = 1 |
| 43 | +size_y = br.metadata.images[0].pixels.physical_size_y |
| 44 | +if size_y is None: |
| 45 | + size_y = 1 |
| 46 | +size_z = br.metadata.images[0].pixels.physical_size_z |
| 47 | +if size_z is None: |
| 48 | + size_z = 1 |
| 49 | +x_unit = br.metadata.images[0].pixels.physical_size_x_unit |
| 50 | +y_unit = br.metadata.images[0].pixels.physical_size_y_unit |
| 51 | +z_unit = br.metadata.images[0].pixels.physical_size_z_unit |
| 52 | + |
| 53 | +# if the the units are um, convert to nm |
| 54 | +if str(x_unit) == "UnitsLength.MICROMETER": |
| 55 | + size_x *= 1000 |
| 56 | +if str(y_unit) == "UnitsLength.MICROMETER": |
| 57 | + size_y *= 1000 |
| 58 | +if str(z_unit) == "UnitsLength.MICROMETER": |
| 59 | + size_z *= 1000 |
| 60 | + |
| 61 | +num_channels = br.shape[-1] |
| 62 | +data_type = "uint16" |
| 63 | +chunk_size = [256, 256, 128] |
| 64 | +volume_size = [br.shape[1], br.shape[0], br.shape[2]] # XYZ |
| 65 | + |
| 66 | +# %% Setup somewhere to hold progress |
| 67 | +progress_dir = OUTPUT_PATH / "progress" |
| 68 | +progress_dir.mkdir(exist_ok=True) |
| 69 | + |
| 70 | +# %% Functions for moving data |
| 71 | +shape = np.array([br.shape[1], br.shape[0], br.shape[2]]) |
| 72 | +chunk_shape = np.array([1024, 1024, 512]) # this is for reading data |
| 73 | +num_chunks_per_dim = np.ceil(shape / chunk_shape).astype(int) |
| 74 | + |
| 75 | + |
| 76 | +def chunked_reader(x_i, y_i, z_i): |
| 77 | + x_start, x_end = x_i * chunk_shape[0], min((x_i + 1) * chunk_shape[0], shape[0]) |
| 78 | + y_start, y_end = y_i * chunk_shape[1], min((y_i + 1) * chunk_shape[1], shape[1]) |
| 79 | + z_start, z_end = z_i * chunk_shape[2], min((z_i + 1) * chunk_shape[2], shape[2]) |
| 80 | + |
| 81 | + # Read the chunk from the BioReader |
| 82 | + chunk = br.read(X=(x_start, x_end), Y=(y_start, y_end), Z=(z_start, z_end)) |
| 83 | + |
| 84 | + # Return the chunk |
| 85 | + return chunk.swapaxes(0, 1) |
| 86 | + |
| 87 | + |
| 88 | +def process(args): |
| 89 | + x_i, y_i, z_i = args |
| 90 | + start = [x_i * chunk_shape[0], y_i * chunk_shape[1], z_i * chunk_shape[2]] |
| 91 | + end = [ |
| 92 | + min((x_i + 1) * chunk_shape[0], shape[0]), |
| 93 | + min((y_i + 1) * chunk_shape[1], shape[1]), |
| 94 | + min((z_i + 1) * chunk_shape[2], shape[2]), |
| 95 | + ] |
| 96 | + f_name = ( |
| 97 | + progress_dir |
| 98 | + / f"{start[0]}-{end[0]}_{start[1]}-{end[1]}_{start[2]}-{end[2]}.done" |
| 99 | + ) |
| 100 | + if f_name.exists() and not OVERWRITE: |
| 101 | + return |
| 102 | + print("Working on", f_name) |
| 103 | + rawdata = chunked_reader(x_i, y_i, z_i) |
| 104 | + # TEMP |
| 105 | + print(np.unique(rawdata)) |
| 106 | + print(rawdata.shape) |
| 107 | + # Create the segmentation mask |
| 108 | + dimensions = [start, end] |
| 109 | + seg_chunk = create_segmentation_chunk(rawdata, dimensions, convert_non_zero_to=None) |
| 110 | + seg_chunk.write_to_directory(OUTPUT_PATH / "data") |
| 111 | + touch(f_name) |
| 112 | + |
| 113 | + |
| 114 | +# %% Try with a single chunk to see if it works |
| 115 | +x_i, y_i, z_i = 0, 0, 0 |
| 116 | +process((x_i, y_i, z_i)) |
| 117 | + |
| 118 | +# %% Loop over all the chunks |
| 119 | +coords = itertools.product( |
| 120 | + range(num_chunks_per_dim[0]), |
| 121 | + range(num_chunks_per_dim[1]), |
| 122 | + range(num_chunks_per_dim[2]), |
| 123 | +) |
| 124 | +# Do it in reverse order because the last chunks are most likely to error |
| 125 | +reversed_coords = list(coords) |
| 126 | +reversed_coords.reverse() |
| 127 | + |
| 128 | +# %% Create the metadata |
| 129 | +metadata = { |
| 130 | + "@type": "neuroglancer_multiscale_volume", |
| 131 | + "data_type": "uint32", |
| 132 | + "num_channels": 1, |
| 133 | + "scales": [ |
| 134 | + { |
| 135 | + "chunk_sizes": [chunk_size], |
| 136 | + "encoding": "compressed_segmentation", |
| 137 | + "copmressed_segmentation_block_size": [8, 8, 8], |
| 138 | + "resolution": [size_x, size_y, size_z], |
| 139 | + "key": "data", |
| 140 | + "size": volume_size, |
| 141 | + } |
| 142 | + ], |
| 143 | + "mesh": "mesh", |
| 144 | + "type": "segmentation", |
| 145 | +} |
| 146 | +write_metadata( |
| 147 | + metadata, |
| 148 | + OUTPUT_PATH, |
| 149 | + overwrite=OVERWRITE, |
| 150 | +) |
| 151 | + |
| 152 | +# %% Now create the mesh |
| 153 | +mesh_shape = np.array(volume_size) |
| 154 | +generate_multiresolution_mesh_from_segmentation(OUTPUT_PATH, "mesh", 3, mesh_shape) |
| 155 | + |
| 156 | + |
| 157 | +# %% Move the data across with a single worker |
| 158 | +for coord in reversed_coords: |
| 159 | + process(coord) |
0 commit comments