From b4da925330b11f7107eceb856983a28fb6170ce0 Mon Sep 17 00:00:00 2001 From: David Manthey Date: Mon, 18 Dec 2023 13:36:24 -0500 Subject: [PATCH] Support multi-source affine transforms --- CHANGELOG.md | 8 + large_image/tilesource/utilities.py | 28 ++- .../large_image_source_multi/__init__.py | 187 +++++++++++++++--- test/test_files/multi_affine.yml | 72 +++++++ test/test_source_multi.py | 34 ++++ 5 files changed, 295 insertions(+), 34 deletions(-) create mode 100644 test/test_files/multi_affine.yml diff --git a/CHANGELOG.md b/CHANGELOG.md index 81530e1e9..ced9fcad2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,13 @@ # Change Log +## 1.27.0 + +### Features +- Support affine transforms in the multi-source ([#1415](../../pull/1415)) + +### Improvements +- Remove NaN values from band information ([#1414](../../pull/1414)) + ## 1.26.3 ### Improvements diff --git a/large_image/tilesource/utilities.py b/large_image/tilesource/utilities.py index 20e70bc93..146fe4ccf 100644 --- a/large_image/tilesource/utilities.py +++ b/large_image/tilesource/utilities.py @@ -704,14 +704,28 @@ def getAvailableNamedPalettes(includeColors=True, reduced=False): return sorted(palettes) +def fullAlphaValue(arr): + """ + Given a numpy array, return the value that should be used for a fully + opaque alpha channel. For uint variants, this is the max value. + + :param arr: a numpy array. + :returns: the value for the alpha channel. + """ + if arr.dtype.kind == 'u': + return np.iinfo(arr.dtype).max + return 1 + + def _makeSameChannelDepth(arr1, arr2): """ Given two numpy arrays that are either two or three dimensions, make the third dimension the same for both of them. Specifically, if they are two - dimensions, first convert to trhee dimensions with a single final value. + dimensions, first convert to three dimensions with a single final value. Otherwise, the dimensions are assumed to be channels of L, LA, RGB, RGBA, or . If L is needed to change to RGB, it is repeated (LLL). - Missing A channels are filled with 1. + Missing A channels are filled with 255, 65535, or 1 depending on if the + dtype is uint8, uint16, or something else. :param arr1: one array to compare. :param arr2: a second array to compare. @@ -729,7 +743,7 @@ def _makeSameChannelDepth(arr1, arr2): for key, arr in arrays.items(): other = arrays['arr1' if key == 'arr2' else 'arr2'] if arr.shape[2] < 3 and other.shape[2] >= 3: - newarr = np.ones((arr.shape[0], arr.shape[1], arr.shape[2] + 2)) + newarr = np.ones((arr.shape[0], arr.shape[1], arr.shape[2] + 2), dtype=arr.dtype) newarr[:, :, 0:1] = arr[:, :, 0:1] newarr[:, :, 1:2] = arr[:, :, 0:1] newarr[:, :, 2:3] = arr[:, :, 0:1] @@ -741,9 +755,11 @@ def _makeSameChannelDepth(arr1, arr2): for key, arr in arrays.items(): other = arrays['arr1' if key == 'arr2' else 'arr2'] if arr.shape[2] < other.shape[2]: - newarr = np.ones((arr.shape[0], arr.shape[1], other.shape[2])) - newarr[:, :, :arr.shape[2]] = arr - arrays[key] = newarr + arrays[key] = np.pad( + arr, + ((0, 0), (0, 0), (0, other.shape[2] - arr.shape[2])), + mode='constant', + constant_values=fullAlphaValue(arr)) return arrays['arr1'], arrays['arr2'] diff --git a/sources/multi/large_image_source_multi/__init__.py b/sources/multi/large_image_source_multi/__init__.py index 3317d2a1d..c69c07e2e 100644 --- a/sources/multi/large_image_source_multi/__init__.py +++ b/sources/multi/large_image_source_multi/__init__.py @@ -19,7 +19,7 @@ from large_image.constants import TILE_FORMAT_NUMPY, SourcePriority from large_image.exceptions import TileSourceError, TileSourceFileNotFoundError from large_image.tilesource import FileTileSource -from large_image.tilesource.utilities import _makeSameChannelDepth +from large_image.tilesource.utilities import _makeSameChannelDepth, fullAlphaValue try: __version__ = _importlib_version(__name__) @@ -752,6 +752,8 @@ def _collectFrames(self): ts = self._openSource(source) self.tileWidth = self.tileWidth or ts.tileWidth self.tileHeight = self.tileHeight or ts.tileHeight + if not hasattr(self, '_firstdtype'): + self._firstdtype = ts.dtype if not numChecked: tsMag = ts.getNativeMagnification() for key in self._nativeMagnification: @@ -920,26 +922,157 @@ def _mergeTiles(self, base, tile, x, y): return tile if base is None: base = np.zeros((0, 0, tile.shape[2]), dtype=tile.dtype) + if tile.shape[2] in {2, 4}: + base[:, :, -1] = fullAlphaValue(tile) base, tile = _makeSameChannelDepth(base, tile) if base.shape[0] < tile.shape[0] + y: vfill = np.zeros( (tile.shape[0] + y - base.shape[0], base.shape[1], base.shape[2]), dtype=base.dtype) if base.shape[2] == 2 or base.shape[2] == 4: - vfill[:, :, -1] = 1 + vfill[:, :, -1] = fullAlphaValue(base) base = np.vstack((base, vfill)) if base.shape[1] < tile.shape[1] + x: hfill = np.zeros( (base.shape[0], tile.shape[1] + x - base.shape[1], base.shape[2]), dtype=base.dtype) if base.shape[2] == 2 or base.shape[2] == 4: - hfill[:, :, -1] = 1 + hfill[:, :, -1] = fullAlphaValue(base) base = np.hstack((base, hfill)) if base.flags.writeable is False: base = base.copy() - base[y:y + tile.shape[0], x:x + tile.shape[1], :] = tile + if base.shape[2] == 2 or base.shape[2] == 4: + baseA = base[y:y + tile.shape[0], x:x + tile.shape[1], -1].astype( + float) / fullAlphaValue(base) + tileA = tile[:, :, -1].astype(float) / fullAlphaValue(tile) + outA = tileA + baseA * (1 - tileA) + base[y:y + tile.shape[0], x:x + tile.shape[1], :-1] = ( + tile[:, :, :-1] * tileA[..., np.newaxis] + + base[y:y + tile.shape[0], x:x + tile.shape[1], :-1] * baseA[..., np.newaxis] * + (1 - tileA[..., np.newaxis]) + ) / outA[..., np.newaxis] + base[y:y + tile.shape[0], x:x + tile.shape[1], -1] = outA * fullAlphaValue(base) + else: + base[y:y + tile.shape[0], x:x + tile.shape[1], :] = tile return base + def _getTransformedTile(self, ts, transform, corners, scale, frame, crop=None): + """ + Determine where the target tile's corners are located on the source. + Fetch that so that we have at least sqrt(2) more resolution, then use + scikit-image warp to transform it. scikit-image does a better job than + scipy.ndimage.affine_transform. + + :param ts: the source of the image to transform. + :param transform: a 3x3 affine 2d matrix for transforming the source + image at full resolution to the target tile at full resolution. + :param corners: corners of the destination in full res coordinates. + corner 0 must be the upper left, 2 must be the lower right. + :param scale: scaling factor from full res to the target resolution. + :param frame: frame number of the source image. + :param crop: an optional dictionary to crop the source image in full + resolution, untransformed coordinates. This may contain left, top, + right, and bottom values in pixels. + :returns: a numpy array tile or None, x, y coordinates within the + target tile for the placement of the numpy tile array. + """ + try: + import skimage.transform + except ImportError: + msg = 'scikit-image is required for affine transforms.' + raise TileSourceError(msg) + # From full res source to full res destination + transform = transform.copy() + # Scale dest corners to actual size; adjust transform for the same + corners = np.array(corners) + corners[:, :2] //= scale + transform[:2, :] /= scale + # Offset so our target is the actual destination array we use + transform[0][2] -= corners[0][0] + transform[1][2] -= corners[0][1] + corners[:, :2] -= corners[0, :2] + outw, outh = corners[2][0], corners[2][1] + if not outh or not outw: + return None, 0, 0 + srccorners = np.dot(np.linalg.inv(transform), np.array(corners).T).T.tolist() + minx = min(c[0] for c in srccorners) + maxx = max(c[0] for c in srccorners) + miny = min(c[1] for c in srccorners) + maxy = max(c[1] for c in srccorners) + srcscale = max((maxx - minx) / outw, (maxy - miny) / outh) + # we only need every 1/srcscale pixel. + srcscale = int(2 ** math.log2(max(1, srcscale))) + # Pad to reduce edge effects at tile boundaries + border = int(math.ceil(2 * srcscale)) + region = { + 'left': int(max(0, minx - border) // srcscale) * srcscale, + 'top': int(max(0, miny - border) // srcscale) * srcscale, + 'right': int((min(ts.sizeX, maxx + border) + srcscale - 1) // srcscale) * srcscale, + 'bottom': int((min(ts.sizeY, maxy + border) + srcscale - 1) // srcscale) * srcscale, + } + if crop: + region['left'] = max(region['left'], crop.get('left', region['left'])) + region['top'] = max(region['top'], crop.get('top', region['top'])) + region['right'] = min(region['right'], crop.get('right', region['right'])) + region['bottom'] = min(region['bottom'], crop.get('bottom', region['bottom'])) + output = { + 'maxWidth': (region['right'] - region['left']) // srcscale, + 'maxHeight': (region['bottom'] - region['top']) // srcscale, + } + if output['maxWidth'] <= 0 or output['maxHeight'] <= 0: + return None, 0, 0 + srcImage = ts.getRegion( + format=TILE_FORMAT_NUMPY, region=region, output=output, frame=frame)[0] + # This is the region we actually took in our source coordinates, scaled + # for if we took a low res version + regioncorners = np.array([ + [region['left'] // srcscale, region['top'] // srcscale, 1], + [region['right'] // srcscale, region['top'] // srcscale, 1], + [region['right'] // srcscale, region['bottom'] // srcscale, 1], + [region['left'] // srcscale, region['bottom'] // srcscale, 1]], dtype=float) + # adjust our transform if we took a low res version of the source + transform[:2, :2] *= srcscale + # Find where the source corners land on the destination. + preshiftcorners = (np.dot(transform, regioncorners.T).T).tolist() + regioncorners[:, :2] -= regioncorners[0, :2] + destcorners = (np.dot(transform, regioncorners.T).T).tolist() + offsetx, offsety = None, None + for idx in range(4): + if offsetx is None or destcorners[idx][0] < offsetx: + x = preshiftcorners[idx][0] + offsetx = destcorners[idx][0] - (x - math.floor(x)) + if offsety is None or destcorners[idx][1] < offsety: + y = preshiftcorners[idx][1] + offsety = destcorners[idx][1] - (y - math.floor(y)) + transform[0][2] -= offsetx + transform[1][2] -= offsety + x, y = int(math.floor(x)), int(math.floor(y)) + # Recompute where the source corners will land + destcorners = (np.dot(transform, regioncorners.T).T).tolist() + destsize = (max(math.ceil(c[0]) for c in destcorners), + max(math.ceil(c[1]) for c in destcorners)) + # Add an alpha band if needed + if srcImage.shape[2] in {1, 3}: + _, srcImage = _makeSameChannelDepth(np.zeros((1, 1, srcImage.shape[2] + 1)), srcImage) + # Add enough space to warp the source image in place + srcImage = np.pad(srcImage, ( + (0, max(0, destsize[1] - srcImage.shape[0])), + (0, max(0, destsize[0] - srcImage.shape[1])), + (0, 0)), mode='constant') + destImage = skimage.transform.warp( + srcImage.astype(float), + skimage.transform.AffineTransform(np.linalg.inv( + transform))).astype(srcImage.dtype) + # Trim to target size and location + destImage = destImage[max(0, -y):, max(0, -x):, :] + x += max(0, -x) + y += max(0, -y) + if x >= outw or y >= outh: + return None, None, None + destImage = destImage[:min(destImage.shape[0], outh - y), + :min(destImage.shape[1], outw - x), :] + return destImage, x, y + def _addSourceToTile(self, tile, sourceEntry, corners, scale): """ Add a source to the current tile. @@ -964,17 +1097,20 @@ def _addSourceToTile(self, tile, sourceEntry, corners, scale): corners[2][1] <= bbox['top'] or corners[0][1] >= bbox['bottom']): return tile transform = bbox.get('transform') - srccorners = ( - list(np.dot(bbox['inverse'], np.array(corners).T).T) - if transform is not None else corners) x = y = 0 # If there is no transform or the diagonals are positive and there is - # no sheer, use getRegion with an appropriate size (be wary of edges) - if (transform is None or + # no sheer and integer pixel alignment, use getRegion with an + # appropriate size + scaleX = transform[0][0] if transform is not None else 1 + scaleY = transform[1][1] if transform is not None else 1 + if ((transform is None or ( transform[0][0] > 0 and transform[0][1] == 0 and - transform[1][0] == 0 and transform[1][1] > 0): - scaleX = transform[0][0] if transform is not None else 1 - scaleY = transform[1][1] if transform is not None else 1 + transform[1][0] == 0 and transform[1][1] > 0 and + transform[0][2] % scaleX == 0 and transform[1][2] % scaleY == 0)) and + (scaleX % scale) == 0 and (scaleY % scale) == 0): + srccorners = ( + list(np.dot(bbox['inverse'], np.array(corners).T).T) + if transform is not None else corners) region = { 'left': srccorners[0][0], 'top': srccorners[0][1], 'right': srccorners[2][0], 'bottom': srccorners[2][1], @@ -1005,21 +1141,12 @@ def _addSourceToTile(self, tile, sourceEntry, corners, scale): sourceTile, _ = ts.getRegion( region=region, output=output, frame=sourceEntry.get('frame', 0), format=TILE_FORMAT_NUMPY) - # Otherwise, determine where the target tile's corners are located on - # the source; fetch that so that we have at least sqrt(2) more - # resolution, then use scikit-image warp to transform it. scikit-image - # does a better job than scipy.ndimage.affine_transform. else: - # try: - # import skimge.transform - # except ImportError: - # msg = 'scikit-image is required for affine transforms.' - # raise TileSourceError(msg) - msg = 'Not implemented' - raise TileSourceError(msg) - # Crop - # TODO - tile = self._mergeTiles(tile, sourceTile, x, y) + sourceTile, x, y = self._getTransformedTile( + ts, transform, corners, scale, sourceEntry.get('frame', 0), + source.get('position', {}).get('crop')) + if sourceTile is not None and all(dim > 0 for dim in sourceTile.shape): + tile = self._mergeTiles(tile, sourceTile, x, y) return tile @methodcache() @@ -1060,7 +1187,9 @@ def getTile(self, x, y, z, pilImageAllowed=False, numpyAllowed=False, **kwargs): if fill: colors = self._info.get('backgroundColor') if colors: - tile = np.full((self.tileHeight, self.tileWidth, len(colors)), colors) + tile = np.full((self.tileHeight, self.tileWidth, len(colors)), + colors, + dtype=getattr(self, '_firstdtype', np.uint8)) # Add each source to the tile for sourceEntry in sourceList: tile = self._addSourceToTile(tile, sourceEntry, corners, scale) @@ -1068,7 +1197,9 @@ def getTile(self, x, y, z, pilImageAllowed=False, numpyAllowed=False, **kwargs): # TODO number of channels? colors = self._info.get('backgroundColor', [0]) if colors: - tile = np.full((self.tileHeight, self.tileWidth, len(colors)), colors) + tile = np.full((self.tileHeight, self.tileWidth, len(colors)), + colors, + dtype=getattr(self, '_firstdtype', np.uint8)) # We should always have a tile return self._outputTile(tile, TILE_FORMAT_NUMPY, x, y, z, pilImageAllowed, numpyAllowed, **kwargs) diff --git a/test/test_files/multi_affine.yml b/test/test_files/multi_affine.yml new file mode 100644 index 000000000..b2ff82283 --- /dev/null +++ b/test/test_files/multi_affine.yml @@ -0,0 +1,72 @@ +--- +name: Multi orientation +description: A test multi file +scale: + mm_x: 0.0005 + mm_y: 0.0005 +width: 2048 +height: 1540 +tileWidth: 256 +tileHeight: 256 +backgroundColor: + - 0 + - 0 + - 255 +sources: + - path: ./test_orient1.tif + z: 0 + - path: ./test_orient2.tif + z: 1 + position: + scale: 4 + - path: ./test_orient3.tif + z: 2 + position: + scale: 0.25 + - path: ./test_orient4.tif + z: 3 + position: + x: 137 + y: 32 + scale: 4 + - path: ./test_orient5.tif + z: 4 + position: + x: 524 + y: 768 + s11: 0.866 + s12: 0.5 + s21: -0.5 + s22: 0.866 + - path: ./test_orient6.tif + z: 5 + position: + x: 1024 + y: 768 + scale: 2 + s11: 0.866 + s12: 0.5 + s21: -0.5 + s22: 0.866 + - path: ./test_orient7.tif + z: 6 + position: + x: 1024 + y: 768 + s11: 2.598 + s12: 1.5 + s21: -1 + s22: 1.732 + - path: ./test_orient8.tif + z: 7 + position: + x: 1024 + y: 768 + s11: 0.866 + s12: -0.5 + s21: 0.5 + s22: 0.866 + crop: + left: 4 + top: 6 + right: 170 diff --git a/test/test_source_multi.py b/test/test_source_multi.py index 93da891e1..e47b8242c 100644 --- a/test/test_source_multi.py +++ b/test/test_source_multi.py @@ -1,8 +1,10 @@ +import io import json import os import large_image_source_multi import numpy as np +import PIL.Image import pytest import large_image @@ -262,3 +264,35 @@ def testStyleFrameBase(): assert image == imageB imageC = source.getTile(0, 0, 2, frame=0) assert image == imageC + + +def testMultiAffineTransform(): + testDir = os.path.dirname(os.path.realpath(__file__)) + imagePath = os.path.join(testDir, 'test_files', 'multi_affine.yml') + source = large_image_source_multi.open(imagePath) + frameval = [ + [243, 243, 49065], + [3920, 3920, 47672], + [13, 13, 49149], + [3943, 3943, 47695], + [246, 246, 49063], + [972, 972, 48785], + [1455, 1455, 48597], + [204, 204, 49059], + ] + for frame in range(8): + thumbs = {} + for res in {256, 512, 1024, 2048}: + img = PIL.Image.open(io.BytesIO(source.getThumbnail( + frame=frame, width=res, encoding='PNG')[0])) + img = np.asarray(img.resize((256, 192), PIL.Image.LANCZOS).convert('RGB')) + thumbs[res] = img + r = np.sum(img[:, :, 0] >= 128) + g = np.sum(img[:, :, 1] >= 128) + b = np.sum(img[:, :, 2] >= 128) + assert abs(r - frameval[frame][0]) < 200 + assert abs(g - frameval[frame][1]) < 200 + assert abs(b - frameval[frame][2]) < 200 + assert abs(np.average(thumbs[2048] - thumbs[256])) < 5 + assert abs(np.average(thumbs[1024] - thumbs[256])) < 5 + assert abs(np.average(thumbs[512] - thumbs[256])) < 5