Rewritten TiffSource using dask/zarr functionality

Added high level Dask function
Completed Image Generator
Improved ome-zarr metadata for napari support

CHANGELOG.md

@@ -1,3 +1,9 @@
+#### Version 0.6.6
+- Rewritten TiffSource using dask/zarr functionality
+- Added high level Dask function
+- Completed Image Generator
+- Improved ome-zarr metadata for napari support
+
 #### Version 0.6.5
 - Expanded Omero image selection
 - Improved image info format, including pyramid sizes

OmeSliCC/GeneratorSource.py

@@ -1,7 +1,4 @@
-import dask
-import dask.array as da
 import numpy as np
-from tqdm import tqdm
 
 from OmeSliCC.OmeSource import OmeSource
 from OmeSliCC.OmeZarr import OmeZarr
@@ -11,15 +8,22 @@ class GeneratorSource(OmeSource):
     def __init__(self, size, tile_size, dtype=np.uint8, source_pixel_size: list = None, seed=None):
         super().__init__()
 
-        size_xyzct = size
-        while len(size_xyzct) < 5:
-            size_xyzct = list(size_xyzct) + [1]
-        while len(tile_size) < len(size):
-            tile_size = list(tile_size) + [1]
+        size_xyzct = list(size)
+        tile_shape = list(np.flip(tile_size))
+        if len(size_xyzct) < 3:
+            size_xyzct += [1]
+        if len(tile_shape) < 3:
+            tile_shape += [1]
+        size_xyzct += [3]
+        tile_shape += [3]
+        if len(size_xyzct) < 5:
+            size_xyzct += [1]
+        if len(tile_shape) < 5:
+            tile_shape += [1]
         self.size = size
         self.sizes = [size[:2]]
         self.sizes_xyzct = [size_xyzct]
-        self.tile_size = tile_size
+        self.tile_shape = tile_shape
         dtype = np.dtype(dtype)
         self.dtype = dtype
         self.pixel_types = [dtype]
@@ -32,72 +36,62 @@ def __init__(self, size, tile_size, dtype=np.uint8, source_pixel_size: list = No
         else:
             self.max_val = 1
 
-        ranges = np.flip(np.ceil(np.divide(size, tile_size)).astype(int))
-        self.tile_indices = list(np.ndindex(tuple(ranges)))
-
         if seed is not None:
             np.random.seed(seed)
 
         self.color_value_table = [np.sin(np.divide(range(dim), dim, dtype=np.float32) * np.pi)
                                   for dim in np.flip(size)]
 
-        # self.noise = np.random.random(size=shape) - 0.5     # uniform
-        self.noise = np.random.normal(loc=0, scale=0.1, size=np.flip(tile_size))  # gaussian
-
     def _find_metadata(self):
         self._get_ome_metadate()
 
-    def calc_color(self, *args):
+    def calc_color(self, *args, **kwargs):
         channels = []
         channel = None
+        range0 = kwargs['range0']
         for index, value in enumerate(reversed(args)):
-            #channel = np.sin((value + self.range0[index]) / self.size[index] * np.pi)
-            channel = self.color_value_table[index][value + self.range0[index]]
+            #channel = np.sin((value + range0[index]) / self.size[index] * np.pi)
+            channel = self.color_value_table[index][value + range0[index]]
             channels.append(channel)
         while len(channels) < 3:
             channels.append(channel)
         return np.stack(channels, axis=-1)
 
-    def get_tile(self, indices, channels_last=False):
+    def get_tile(self, indices, tile_size=None):
         # tile in (z,),y,x,c
-        self.range0 = np.flip(indices[1:]) * self.tile_size
-        self.range1 = np.min([self.range0 + self.tile_size, self.size], 0)
-        shape = list(reversed(self.range1 - self.range0))
-        tile = np.fromfunction(self.calc_color, shape, dtype=int)
+        if not tile_size:
+            tile_size = np.flip(self.tile_shape)
+        range0 = np.flip(indices)
+        range1 = np.min([range0 + tile_size, self.size], 0)
+        shape = list(reversed(range1 - range0))
+        tile = np.fromfunction(self.calc_color, shape, dtype=int, range0=range0)
         # apply noise to each channel separately
         for channeli in range(3):
-            tile[..., channeli] = np.clip(tile[..., channeli] + self.noise, 0, 1)
+            noise = np.random.random(size=shape) - 0.5
+            tile[..., channeli] = np.clip(tile[..., channeli] + noise, 0, 1)
         if self.dtype.kind != 'f':
             tile *= self.max_val
         tile = tile.astype(self.dtype)
-        if not channels_last:
-            tile = np.moveaxis(tile, -1, 0)
         return tile
 
-    def get_tiles(self, channels_last=False):
-        for indices in tqdm(self.tile_indices):
-            yield self.get_tile(indices, channels_last)
-
-    def _get_output_dask(self):
-        # TODO: fix dimensions / shapes
-
-        #data = da.fromfunction(lambda indices: self.get_tile(indices), self.tile_indices, shape=self.tile_size, dtype=self.dtype)
-
-        delayed_reader = dask.delayed(self.get_tile)
-        tile_shape = [3] + list(np.flip(self.tile_size))
-        dask_tiles = []
-        for indices in self.tile_indices:
-            indices_tzyx = indices
-            while len(indices_tzyx) < 4:
-                indices_tzyx = [0] + list(indices_tzyx)
-            dask_tile = da.from_delayed(delayed_reader(indices_tzyx), shape=tile_shape, dtype=self.dtype)
-            dask_tiles.append(dask_tile)
-        dask_data = da.block(dask_tiles)
-        return dask_data
-
     def _asarray_level(self, level: int, x0: float = 0, y0: float = 0, x1: float = -1, y1: float = -1,
                        c: int = None, z: int = None, t: int = None) -> np.ndarray:
-        return self.get_output_dask()
+        # ignore c
+        indices = [y0, x0]
+        tile_size = [y1 - y0, x1 - x0]
+        if z:
+            indices = [0] + indices
+            tile_size = [1] + tile_size
+        if t:
+            indices = [0] + indices
+            tile_size = [1] + tile_size
+        data = self.get_tile(indices, tile_size)
+        if not z:
+            data = np.expand_dims(data, 0)
+        data = np.moveaxis(data, -1, 0)
+        if not t:
+            data = np.expand_dims(data, 0)
+        return data
 
 
 if __name__ == '__main__':
@@ -115,8 +109,8 @@ def _asarray_level(self, level: int, x0: float = 0, y0: float = 0, x1: float = -
     generator = GeneratorSource(size, tile_size, dtype, pixel_size, seed)
     print('init done')
 
-    #save_tiff('D:/slides/test.ome.tiff', generator.get_tiles(channels_last=True), shape, dtype, tile_size=tile_shape)
+    #save_tiff('D:/slides/test.ome.tiff', generator.asarray(), shape, dtype, tile_size=tile_shape)
     zarr = OmeZarr('D:/slides/test.ome.zarr')
-    zarr.write(generator.get_output_dask(), generator, tile_size=tile_size, npyramid_add=3, pyramid_downsample=2)
+    zarr.write(generator.asarray(), generator, tile_size=tile_size, npyramid_add=3, pyramid_downsample=2)
 
     print('done')

OmeSliCC/OmeSource.py

@@ -1,3 +1,4 @@
+import dask.array
 import logging
 import numpy as np
 
@@ -152,13 +153,6 @@ def _init_sizes(self):
     def get_source_dask(self):
         raise NotImplementedError('Implement method in subclass')
 
-    def get_output_dask(self):
-        data = self._get_output_dask()
-        return redimension_data(data, self.dimension_order, self.get_dimension_order())
-
-    def _get_output_dask(self):
-        raise NotImplementedError('Implement method in subclass')
-
     def get_mag(self) -> float:
         mag = self.source_mag
         # get effective mag at target pixel size
@@ -228,7 +222,7 @@ def get_channel_window(self, channeli):
         min_quantile = 0.001
         max_quantile = 0.999
 
-        if 'window' in self.channels[channeli]:
+        if channeli < len(self.channels) and 'window' in self.channels[channeli]:
             return self.channels[channeli].get('window')
 
         dtype = self.get_pixel_type()
@@ -241,7 +235,7 @@ def get_channel_window(self, channeli):
         nsizes = len(self.sizes)
         if nsizes > 1:
             image = self._asarray_level(nsizes - 1)
-            image = image[:, channeli:channeli+1, ...]
+            image = np.asarray(image[:, channeli:channeli+1, ...])
             min, max = get_image_quantile(image, min_quantile), get_image_quantile(image, max_quantile)
         else:
             # do not query full size image
@@ -343,6 +337,48 @@ def asarray(self, x0: float = 0, y0: float = 0, x1: float = -1, y1: float = -1,
             image = image0
         return image
 
+    def asarray_um(self, x_um: float, y_um: float, w_um: float, h_um: float,
+                   c: int = None, z: int = None, t: int = None):
+        pixel_size = self.get_pixel_size_micrometer()[:2]
+        x0, y0 = np.divide([x_um, y_um], pixel_size)
+        w, h = np.divide([w_um, h_um], pixel_size)
+        x1, y1 = x0 + w, y0 + h
+        return self.asarray(x0, y0, x1, y1, c, z, t)
+
+    def asdask(self, chunk_size: tuple) -> da.Array:
+        chunk_shape = list(np.flip(chunk_size))
+        while len(chunk_shape) < 3:
+            chunk_shape = [1] + chunk_shape
+        chunk_shape = [self.get_nchannels()] + chunk_shape
+        while len(chunk_shape) < 5:
+            chunk_shape = [1] + chunk_shape
+        chunks = np.ceil(np.flip(self.get_size_xyzct()) / chunk_shape).astype(int)
+
+        delayed_reader = dask.delayed(self.asarray)
+        dtype = self.get_pixel_type()
+
+        dask_times = []
+        for ti in range(chunks[0]):
+            dask_planes = []
+            for zi in range(chunks[2]):
+                dask_rows = []
+                for yi in range(chunks[3]):
+                    dask_tiles = []
+                    for xi in range(chunks[4]):
+                        x0, x1 = xi * chunk_shape[4], (xi + 1) * chunk_shape[4]
+                        y0, y1 = yi * chunk_shape[3], (yi + 1) * chunk_shape[3]
+                        z = zi * chunk_shape[2]
+                        t = ti * chunk_shape[0]
+                        dask_tile = da.from_delayed(delayed_reader(x0, y0, x1, y1, z=z, t=t), shape=chunk_shape,
+                                                    dtype=dtype)
+                        dask_tiles.append(dask_tile)
+                    dask_rows.append(da.concatenate(dask_tiles, axis=4))
+                dask_planes.append(da.concatenate(dask_rows, axis=3))
+            dask_times.append(da.concatenate(dask_planes, axis=1))
+        dask_data = da.concatenate(dask_times, axis=0)
+
+        return dask_data
+
     def clone_empty(self) -> np.ndarray:
         return np.zeros(self.get_shape(), dtype=self.get_pixel_type())
 

OmeSliCC/Omero.py

@@ -97,9 +97,9 @@ def get_annotation_image_ids(self) -> dict:
                 for image_id in self._get_dataset_images(dataset.getId()):
                     images.pop(image_id, None)
 
+        # regex
         for image_id, image in images.items():
             name = image.getName()
-            print(name)
             include = True
             if include_regex:
                 include = False

OmeSliCC/OmeroSource.py

@@ -37,18 +37,26 @@ def __init__(self,
         image_object = self.omero.get_image_object(image_id)
         self.image_object = image_object
 
-        zsize = get_default(image_object.getPixelSizeZ(), 1)
+        zsize = get_default(image_object.getSizeZ(), 1)
         nchannels = np.sum([channel.getLogicalChannel().getSamplesPerPixel() for channel in image_object.getChannels()])
         pixel_type = np.dtype(image_object.getPixelsType())
         # currently only support/output yxc - allow default value
         #self.dimension_order = image_object.getPrimaryPixels().getDimensionOrder().getValue().lower()
+
         self.pixels_store = self.omero.create_pixels_store(image_object)
         for resolution in self.pixels_store.getResolutionDescriptions():
             self.sizes.append((resolution.sizeX, resolution.sizeY))
             self.sizes_xyzct.append((resolution.sizeX, resolution.sizeY, zsize, nchannels, 1))
             self.pixel_types.append(pixel_type)
             self.pixel_nbits.append(pixel_type.itemsize * 8)
 
+        if not self.sizes:
+            xsize, ysize = image_object.getSizeX(), image_object.getSizeY()
+            self.sizes.append((xsize, ysize))
+            self.sizes_xyzct.append((xsize, ysize, zsize, nchannels, 1))
+            self.pixel_types.append(pixel_type)
+            self.pixel_nbits.append(pixel_type.itemsize * 8)
+
         # Omero API issue: pixel store level order not guaranteed
         default_level = self.pixels_store.getResolutionLevel()
         nlevels = self.pixels_store.getResolutionLevels()