From 0e4eb257ce36831c30f1cb001a4e390ff9ea8a78 Mon Sep 17 00:00:00 2001 From: Joost de Folter Date: Thu, 15 Feb 2024 23:50:39 +0100 Subject: [PATCH] Simplified rgb source rendering, improved performance of precise resize operation using skimage function, improved test functions --- CHANGELOG.md | 4 +++ OmeSliCC/OmeSource.py | 33 ++++++------------------ OmeSliCC/image_util.py | 37 +++++++++++++-------------- pyproject.toml | 2 +- tests/conversion_test.py | 17 +++++++------ tests/resize_test.py | 55 ++++++++++++++++++++++++++++++++++++++++ tests/test.py | 34 +++++++++++-------------- 7 files changed, 109 insertions(+), 73 deletions(-) create mode 100644 tests/resize_test.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 5a65bb9..3a3ff08 100644 --- a/CHANGELOG.md +++ b/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) diff --git a/OmeSliCC/OmeSource.py b/OmeSliCC/OmeSource.py index 9d5aec8..606b727 100644 --- a/OmeSliCC/OmeSource.py +++ b/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: diff --git a/OmeSliCC/image_util.py b/OmeSliCC/image_util.py index 1bdf35d..5d5e8b8 100644 --- a/OmeSliCC/image_util.py +++ b/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)) diff --git a/pyproject.toml b/pyproject.toml index fa461c3..16f6431 100644 --- a/pyproject.toml +++ b/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"} diff --git a/tests/conversion_test.py b/tests/conversion_test.py index f6d3f72..ea3caa4 100644 --- a/tests/conversion_test.py +++ b/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) diff --git a/tests/resize_test.py b/tests/resize_test.py new file mode 100644 index 0000000..0961b5f --- /dev/null +++ b/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) diff --git a/tests/test.py b/tests/test.py index eea15eb..bb9c296 100644 --- a/tests/test.py +++ b/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,11 +166,9 @@ 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' @@ -180,7 +176,7 @@ def calc_images_fraction(pattern: str): #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))