Simplified rgb source rendering, improved performance of precise resi…

…ze operation using skimage function, improved test functions

CHANGELOG.md

@@ -1,3 +1,7 @@
+#### Version 0.6.10
+- Simplified rgb source rendering
+- Improved performance of precise resize operation using skimage
+
 #### Version 0.6.9
 - Combine action now working: combine separate channel images into single multichannel image(s)
 - Multi-dimensional image generator now working (GeneratorSource)

OmeSliCC/OmeSource.py

@@ -224,12 +224,14 @@ def get_shape(self, dimension_order: str = None, xyzct: tuple = None) -> tuple:
     def get_thumbnail(self, target_size: tuple, precise: bool = False) -> np.ndarray:
         size, index = get_best_size(self.sizes, target_size)
         scale = np.divide(target_size, self.sizes[index])
-        image, _ = self.get_yxc_image(self._asarray_level(index), t=0, z=0)
+        new_dimension_order = 'yxc'
+        image = redimension_data(self._asarray_level(index), self.get_dimension_order(), new_dimension_order, t=0, z=0)
         if precise:
             thumbnail = precise_resize(image, scale)
         else:
-            thumbnail = self.render(image_resize(image, target_size))
-        return thumbnail
+            thumbnail = image_resize(image, target_size)
+        thumbnail_rgb = self.render(thumbnail, new_dimension_order)
+        return thumbnail_rgb
 
     def get_channel_window(self, channeli):
         min_quantile = 0.001
@@ -255,32 +257,13 @@ def get_channel_window(self, channeli):
             min, max = start, end
         return {'start': start, 'end': end, 'min': min, 'max': max}
 
-    def get_yxc_image(self, image, t=None, z=None, c=None):
-        new_dimension_order = self.get_dimension_order()
-        if z is not None:
-            image = image[:, :, z, ...]
-            new_dimension_order = new_dimension_order.replace('z', '')
-
-        new_dimension_order = new_dimension_order.replace('c', '')
-        if c is not None:
-            image = image[:, c, ...]
-        else:
-            image = np.moveaxis(image, 1, -1)
-            new_dimension_order += 'c'
-
-        if t is not None:
-            image = image[t, ...]
-            new_dimension_order = new_dimension_order.replace('t', '')
-        return image, new_dimension_order
-
-    def render(self, image: np.ndarray, t: int = 0, z: int = 0, channels: list = []) -> np.ndarray:
-        if image.ndim == 5:
-            image, _ = self.get_yxc_image(image, t=t, z=z)
+    def render(self, image: np.ndarray, source_dimension_order: str, t: int = 0, z: int = 0, channels: list = []) -> np.ndarray:
+        image = redimension_data(image, source_dimension_order, 'yxc', t=t, z=z)
         new_image = np.zeros(list(image.shape[:2]) + [3], dtype=np.float32)
         tot_alpha = 0
         n = len(self.channels)
 
-        is_rgb = (self.get_nchannels() == 3 and n == 1)
+        is_rgb = (self.get_nchannels() == 3 and (n <= 1 or n == 3))
         do_normalisation = (image.dtype.itemsize == 2)
 
         if not is_rgb:

OmeSliCC/image_util.py

@@ -8,6 +8,7 @@
 import PIL.Image
 from PIL.ExifTags import TAGS
 import tifffile
+from skimage.transform import downscale_local_mean
 from tifffile import TiffFile
 
 from OmeSliCC.util import *
@@ -83,7 +84,7 @@ def convert_image_sign_type(image0: np.ndarray, dtype: np.dtype) -> np.ndarray:
     return image
 
 
-def redimension_data(data, old_order, new_order, **kwargs):
+def redimension_data(data, old_order, new_order, **indices):
     # able to provide optional dimension values e.g. t=0, z=0
     if new_order == old_order:
         return data
@@ -94,7 +95,7 @@ def redimension_data(data, old_order, new_order, **kwargs):
     for o in old_order:
         if o not in new_order:
             index = order.index(o)
-            dim_value = kwargs.get(o, 0)
+            dim_value = indices.get(o, 0)
             new_data = np.take(new_data, indices=dim_value, axis=index)
             order = order.replace(o, '')
     # add
@@ -246,24 +247,11 @@ def image_resize(image: np.ndarray, target_size0: tuple, dimension_order: str =
     return new_image
 
 
-def precise_resize(image: np.ndarray, scale: np.ndarray, use_max: bool = False) -> np.ndarray:
-    h, w = np.ceil(image.shape[:2] * scale).astype(int)
-    shape = list(image.shape).copy()
-    shape[:2] = h, w
-    new_image = np.zeros(shape, dtype=np.float32)
-    step_size = 1 / scale
-    for y in range(h):
-        for x in range(w):
-            y0, y1 = np.round([y * step_size[1], (y + 1) * step_size[1]]).astype(int)
-            x0, x1 = np.round([x * step_size[0], (x + 1) * step_size[0]]).astype(int)
-            image1 = image[y0:y1, x0:x1]
-            if image1.size > 0:
-                if use_max:
-                    value = np.max(image1, axis=(0, 1))
-                else:
-                    value = np.mean(image1, axis=(0, 1))
-                new_image[y, x] = value
-    return new_image.astype(image.dtype)
+def precise_resize(image: np.ndarray, factors) -> np.ndarray:
+    if image.ndim > len(factors):
+        factors = list(factors) + [1]
+    new_image = downscale_local_mean(np.asarray(image), tuple(factors)).astype(image.dtype)
+    return new_image
 
 
 def create_compression_filter(compression: list) -> tuple:
@@ -414,6 +402,15 @@ def compare_image_dist(image0: np.ndarray, image1: np.ndarray) -> tuple:
     return dif, rgb_max, rgb_mean, psnr
 
 
+def compare_image_dist_simple(image0: np.ndarray, image1: np.ndarray) -> dict:
+    dif = cv.absdiff(image0, image1)
+    psnr = cv.PSNR(image0, image1)
+    rgb_dif = np.linalg.norm(dif, axis=2)
+    dif_max = np.max(rgb_dif)
+    dif_mean = np.mean(rgb_dif)
+    return {'dif_max': dif_max, 'dif_mean': dif_mean, 'psnr': psnr}
+
+
 def calc_fraction_used(image: np.ndarray, threshold: float = 0.1) -> float:
     low = int(round(threshold * 255))
     high = int(round((1 - threshold) * 255))

pyproject.toml

@@ -1,6 +1,6 @@
 [project]
 name = "OmeSliCC"
-version = "0.6.9"
+version = "0.6.10"
 description = "Ome(ro) Slide Image Conversion and Compression pipeline"
 readme = "README.md"
 license = {file = "LICENSE.md"}

tests/conversion_test.py

@@ -17,14 +17,14 @@ def load_as_zarr(path, x0_um, x1_um, y0_um, y1_um):
     source = TiffSource(path)
     data = source.asarray_um(x0=x0_um, x1=x1_um, y0=y0_um, y1=y1_um)
     image = np.asarray(data)
-    image = source.render(image)
+    image = source.render(image, source.get_dimension_order())
     show_image(image)
 
 
 def conversion_test():
-    path1 = 'test1.ome.tiff'
-    path2 = 'test1.ome.zarr'
-    path3 = 'test2.ome.tiff'
+    path1 = 'D:/slides/test1.ome.tiff'
+    path2 = 'D:/slides/test1.ome.zarr'
+    path3 = 'D:/slides/test2.ome.tiff'
     size = [1024, 1024]
     tile_size = [256, 256]
     output_params = {'tile_size': tile_size, 'npyramid_add': 3, 'pyramid_downsample': 2}
@@ -46,7 +46,7 @@ def conversion_test():
 
     data = source3.asarray(pixel_size=[10])
     assert data.ndim == 5
-    image, _ = source3.get_yxc_image(data, t=0, z=0)
+    image = redimension_data(data, source3.get_dimension_order(), 'yxc', t=0, z=0)
     assert image.ndim == 3 and image.shape[2] == nchannels
     x0, y0, x1, y1 = size[0] // 2, size[1] // 2, size[0], size[1]
     w, h = np.divide(size, 2).astype(int)
@@ -58,11 +58,12 @@ def conversion_test():
 
 def check_large_tiff_arrays(input):
     source = TiffSource(input, target_pixel_size=[(10, 'um')])
-    image = source.render(source.asarray())
+    dimension_order = source.get_dimension_order()
+    image = source.render(source.asarray(), dimension_order)
     show_image(image)
-    tile = source.render(source.asarray(x0=1100, x1=2100, y0=1200, y1=1400))
+    tile = source.render(source.asarray(x0=1100, x1=2100, y0=1200, y1=1400), dimension_order)
     show_image(tile)
-    tile = source.render(source.asarray(x0=1100, x1=2100, y0=1200, y1=1400, pixel_size=[10]))
+    tile = source.render(source.asarray(x0=1100, x1=2100, y0=1200, y1=1400, pixel_size=[10]), dimension_order)
     show_image(tile)
 
 

tests/resize_test.py

@@ -0,0 +1,55 @@
+import numpy as np
+import os
+
+from OmeSliCC.TiffSource import TiffSource
+from OmeSliCC.image_util import show_image, compare_image_dist_simple, image_resize, precise_resize
+
+
+def precise_resize0(image: np.ndarray, scale: np.ndarray, use_max: bool = False) -> np.ndarray:
+    h, w = np.ceil(image.shape[:2] * scale).astype(int)
+    shape = list(image.shape).copy()
+    shape[:2] = h, w
+    new_image = np.zeros(shape, dtype=np.float32)
+    step_size = 1 / scale
+    for y in range(h):
+        for x in range(w):
+            y0, y1 = np.round([y * step_size[1], (y + 1) * step_size[1]]).astype(int)
+            x0, x1 = np.round([x * step_size[0], (x + 1) * step_size[0]]).astype(int)
+            image1 = image[y0:y1, x0:x1]
+            if image1.size > 0:
+                if use_max:
+                    value = np.max(image1, axis=(0, 1))
+                else:
+                    value = np.mean(image1, axis=(0, 1))
+                new_image[y, x] = value
+    return new_image.astype(image.dtype)
+
+
+def precise_downscale0(image: np.ndarray, patch_size: tuple) -> np.ndarray:
+    scale = 1 / np.array(patch_size)
+    image2 = precise_resize0(image, scale)
+    return image2
+
+
+def precise_downscale(image: np.ndarray, patch_size: tuple) -> np.ndarray:
+    image2 = precise_resize(image, patch_size)
+    return image2
+
+
+if __name__ == '__main__':
+    path = 'E:/Personal/Crick/slides/test_images/H&E K130_PR003.ome.tiff'
+    patch_size = (256, 256)
+
+    source = TiffSource(path)
+    image = source.render(source.asarray(), source.get_dimension_order())
+    image1 = precise_downscale0(image, patch_size)
+    image2 = precise_downscale(image, patch_size)
+    new_size = tuple(np.flip(image2.shape[:2]))
+    image3 = image_resize(image, new_size)
+    new_size = np.divide(source.get_size(), patch_size)
+    image3b = image_resize(image, new_size)
+    print(compare_image_dist_simple(image1[:-1, :-1], image2[:-1, :-1]))
+    show_image(image1)
+    show_image(image2)
+    show_image(image3)
+    show_image(image3b)

tests/test.py

@@ -14,23 +14,19 @@
 
 
 def render_at_pixel_size(filename: str, source_pixel_size: list = None,
-                         target_pixel_size: list = None, pos: tuple = (0, 0, -1, -1)) -> np.ndarray:
+                         target_pixel_size: list = None, **indices) -> np.ndarray:
     if filename.endswith('.zarr'):
         source = OmeZarrSource(filename, source_pixel_size)
     else:
         source = TiffSource(filename, source_pixel_size)
-    image0 = source.asarray(pos[0], pos[1], pos[2], pos[3], pixel_size=target_pixel_size)
-    image = source.render(image0)
+    image0 = source.asarray(pixel_size=target_pixel_size, **indices)
+    image = source.render(image0, source.get_dimension_order())
     return image
 
 
-def test_load(filename: str, pixel_size: tuple = None, position: tuple = None, size: tuple = None) -> np.ndarray:
+def test_load(filename: str, pixel_size: tuple = None, **indices) -> np.ndarray:
     source = TiffSource(filename, pixel_size)
-    if position is None:
-        position = (0, 0)
-    if size is None:
-        size = source.get_size()
-    image = source.asarray(position[0], position[1], position[0] + size[0], position[1] + size[1])
+    image = source.asarray(**indices)
     return image
 
 
@@ -53,14 +49,16 @@ def test_extract_metadata(path: str):
     print(print_dict(metadata))
 
 
-def compare_image_tiles(filename1: str, filename2: str, bits_per_channel: int = 8, tile_size: tuple = (512, 512)) -> tuple:
+def compare_images(filename1: str, filename2: str, tile_size: tuple = None) -> dict:
     source1 = TiffSource(filename1)
     source2 = TiffSource(filename2)
-
     w, h = source1.get_size()
+    if tile_size is None:
+        tile_size = w, h
     tw, th = tile_size
     nx = int(np.ceil(w / tw))
     ny = int(np.ceil(h / th))
+    bits_per_channel = source1.pixel_nbits[0]
 
     difs_max = []
     difs_mean = []
@@ -72,8 +70,8 @@ def compare_image_tiles(filename1: str, filename2: str, bits_per_channel: int =
             ry = y * th
             rx2 = min(rx + tw, w)
             ry2 = min(ry + th, h)
-            patch1 = source1.asarray(rx, ry, rx2, ry2)
-            patch2 = source2.asarray(rx, ry, rx2, ry2)
+            patch1 = source1.asarray(x0=rx, x1=rx2, y0=ry, y1=ry2)
+            patch2 = source2.asarray(x0=rx, x1=rx2, y0=ry, y1=ry2)
             dif, dif_max, dif_mean, _ = compare_image_dist(patch1, patch2)
             mse = np.mean(dif.astype(float) ** 2)
             difs_max.append(dif_max)
@@ -85,7 +83,7 @@ def compare_image_tiles(filename1: str, filename2: str, bits_per_channel: int =
     maxval = 2 ** bits_per_channel - 1
     psnr = 20 * np.log10(maxval / np.sqrt(np.mean(mses)))    # recalculate PSNR based on mean MSE
     print(f'rgb dist max: {dif_max:.1f} mean: {dif_mean:.1f} PSNR: {psnr:.1f}')
-    return dif_max, dif_mean, psnr
+    return {'dif_max': dif_max, 'dif_mean': dif_mean, 'psnr': psnr}
 
 
 def test_read_source(image_filename: str, n: int = 1000):
@@ -112,7 +110,7 @@ def test_read_source(image_filename: str, n: int = 1000):
         yi = random.randrange(ny)
         x = xi * patch_size[0]
         y = yi * patch_size[1]
-        image = source.asarray(x, y, x + patch_size[0], y + patch_size[1])
+        image = source.asarray(x0=x, x1=x+patch_size[0], y0=y, y1=y+patch_size[1])
         image.shape
         #show_image(image)
     elapsed = timer() - start
@@ -168,19 +166,17 @@ def calc_images_fraction(pattern: str):
 
 
 if __name__ == '__main__':
-    image_dir = 'resources/images/'
-    patch_size = (256, 256)
     os.chdir('../')
-
     #path = 'E:/Personal/Crick/slides/test_images/19629.svs'
+    path = 'E:/Personal/Crick/slides/test_images/H&E K130_PR003.ome.tiff'
     #path = 'D:/slides/Pharos_test_images/01-08-23_test2__33.tiff'
     #path = 'D:/slides/Pharos_test_images/Testing7.tiff'
     #path = 'D:/slides/Pharos_test_images/image_navcam.tiff'
     #path = 'C:/temp/sbemimage_test/workspace/OV000.ome.tif'
     #path = 'C:/temp/sbemimage_test/workspace/old OV000.tif'
     #path = 'D:/slides/EM04613/EM04613_04_20x_WF_Reflection-02-Stitching-01.ome.tif'
     #path = 'D:/slides/EM04613/EM04613_04_20x_WF_Reflection-02-Stitching-01.ome.zarr'
-    path = 'D:/slides/EM04573_01t/stitched_norm.ome.zarr'
+    patch_size = (256, 256)
 
     # perform test
     #print(tiff_info(path))