From 22f95b006b18b2729670c6b7afd025d3b338d435 Mon Sep 17 00:00:00 2001 From: Marc Aubreville Date: Tue, 27 Oct 2026 20:43:14 +0100 Subject: [PATCH] Register the pyramid levels of slideio formats (CZI, VSI) Levels of slideio images were displaced against each other by up to ~0.65 px, which shows as wobbling while the viewer blends zoom levels: - slideio reports 1 / scale (16.0007 for a 16:1 level) but the levels are exact decimation grids; reads at positions that are not multiples of the downsample are displaced by up to a level pixel. SlideIOSlide now reports the exact integer downsample, snaps positions to it and requests size * downsample. - the scale from a slide level to a Deep Zoom level is derived from the downsample (level 0 size / downsample) instead of the stored level size. - CZI pyramids store levels that are displaced against each other (0.2-0.8 px). util.registration measures the displacement once per image (each level against the level below, three textured spots, parabolic refinement, chained to level 0, a level is corrected only if two spots agree) and get_tile applies it as an offset of the fractional sampling box. Level 0 is never shifted. The offsets are stored in the new table images_pyramidregistration (with the file size/mtime and the algorithm version they were measured on). The table is shared with the Rust port of the server, which creates it with the same DDL; the migration uses CREATE TABLE IF NOT EXISTS so either side can create it first. Fixes #165 Co-Authored-By: Claude Sonnet 5.5 --- exact/exact/images/api_views.py | 8 +- .../migrations/0041_pyramidregistration.py | 47 +++++ exact/exact/images/models.py | 17 +- .../exact/images/test_pyramid_registration.py | 57 ++++++ exact/exact/images/views.py | 6 +- exact/util/registration.py | 162 ++++++++++++++++++ exact/util/slide_server.py | 24 ++- exact/util/slideio.py | 34 +++- exact/util/test_registration.py | 94 ++++++++++ 9 files changed, 433 insertions(+), 16 deletions(-) create mode 100644 exact/exact/images/migrations/0041_pyramidregistration.py create mode 100644 exact/exact/images/test_pyramid_registration.py create mode 100644 exact/util/registration.py create mode 100644 exact/util/test_registration.py diff --git a/exact/exact/images/api_views.py b/exact/exact/images/api_views.py index e878375e..7eb71735 100644 --- a/exact/exact/images/api_views.py +++ b/exact/exact/images/api_views.py @@ -44,7 +44,7 @@ from util.cellvizio import ReadableCellVizioMKTDataset # just until data access is pip installable from PIL import Image as PIL_Image -from util.slide_server import SlideCache, SlideFile, PILBytesIO, image_cache +from util.slide_server import SlideCache, SlideFile, PILBytesIO, image_cache, slide_for_image logger = logging.getLogger('django') cache = caches['default'] @@ -157,7 +157,7 @@ def view_image_tile(self, request, pk, z_dimension, frame, level, tile_path): file_path = os.path.join(settings.IMAGE_PATH, image.path(z_dimension, frame)) try: - slide = image_cache.get(file_path) + slide = slide_for_image(image, file_path) tile = slide.get_tile(level, (col, row), frame=frame) @@ -186,7 +186,7 @@ def view_image(self, request, pk, z_dimension, frame): return HttpResponse(value, content_type='application/xml') file_path = os.path.join(settings.IMAGE_PATH, image.path(z_dimension, frame)) - slide = image_cache.get(file_path) + slide = slide_for_image(image, file_path) value = slide.get_dzi("jpeg") if hasattr(cache, "delete_pattern"): @@ -203,7 +203,7 @@ def update_image_cache(self, request, pk=None): z_dimension = int(request.data.get("z_dimension", 1)) frame = int(request.data.get("frame", 1)) - slide = image_cache.get(image.path()) + slide = slide_for_image(image, image.path()) keys = {} results = {} diff --git a/exact/exact/images/migrations/0041_pyramidregistration.py b/exact/exact/images/migrations/0041_pyramidregistration.py new file mode 100644 index 00000000..be780246 --- /dev/null +++ b/exact/exact/images/migrations/0041_pyramidregistration.py @@ -0,0 +1,47 @@ +import django.db.models.deletion +from django.db import migrations, models + +# The Rust port of the server (exact-rust) creates the same table with the same DDL (crates/exact-db/migrations/ +# 0005_pyramid_registration.sql). Whichever side runs first creates it, the other is a no-op: the statement is +# CREATE TABLE IF NOT EXISTS and the model state is added independently of the database operation. +CREATE_TABLE = """ +CREATE TABLE IF NOT EXISTS images_pyramidregistration ( + id bigint GENERATED BY DEFAULT AS IDENTITY PRIMARY KEY, + image_id bigint NOT NULL UNIQUE REFERENCES images_image(id) ON DELETE CASCADE DEFERRABLE INITIALLY DEFERRED, + offsets jsonb NOT NULL, + algorithm smallint NOT NULL, + file_size bigint NOT NULL, + file_mtime bigint NOT NULL, + updated timestamp with time zone NOT NULL DEFAULT now() +); +""" + + +class Migration(migrations.Migration): + + dependencies = [ + ('images', '0040_backfill_crop_source'), + ] + + operations = [ + migrations.SeparateDatabaseAndState( + database_operations=[ + migrations.RunSQL(CREATE_TABLE, reverse_sql='DROP TABLE IF EXISTS images_pyramidregistration;'), + ], + state_operations=[ + migrations.CreateModel( + name='PyramidRegistration', + fields=[ + ('id', models.BigAutoField(auto_created=True, primary_key=True, serialize=False, verbose_name='ID')), + ('offsets', models.JSONField()), + ('algorithm', models.SmallIntegerField()), + ('file_size', models.BigIntegerField()), + ('file_mtime', models.BigIntegerField()), + ('updated', models.DateTimeField(auto_now=True)), + ('image', models.OneToOneField(on_delete=django.db.models.deletion.CASCADE, + related_name='pyramid_registration', to='images.image')), + ], + ), + ], + ), + ] diff --git a/exact/exact/images/models.py b/exact/exact/images/models.py index 2c59db95..1b931d83 100644 --- a/exact/exact/images/models.py +++ b/exact/exact/images/models.py @@ -886,4 +886,19 @@ def convert_coodinates(self, vector): result_vector[f"x{id+1}"] = self.transformation_matrix["t_00"] * x + self.transformation_matrix["t_01"] * y + self.transformation_matrix["t_02"] result_vector[f"y{id+1}"] = self.transformation_matrix["t_10"] * x + self.transformation_matrix["t_11"] * y + self.transformation_matrix["t_12"] - return result_vector \ No newline at end of file + return result_vector + +class PyramidRegistration(models.Model): + """ + Displacement of the levels of an image pyramid against each other (util.registration), measured once per image. + + ``offsets`` holds one ``[dx, dy]`` per pyramid level in pixels of that level; level 0 is always ``[0, 0]``. + Rows are only used for the file they were measured on (``file_size``, ``file_mtime`` in seconds) and for the + current ``algorithm`` version. The table is shared with the Rust port of the server. + """ + image = models.OneToOneField(Image, on_delete=models.CASCADE, related_name='pyramid_registration') + offsets = models.JSONField() + algorithm = models.SmallIntegerField() + file_size = models.BigIntegerField() + file_mtime = models.BigIntegerField() + updated = models.DateTimeField(auto_now=True) diff --git a/exact/exact/images/test_pyramid_registration.py b/exact/exact/images/test_pyramid_registration.py new file mode 100644 index 00000000..fcad8321 --- /dev/null +++ b/exact/exact/images/test_pyramid_registration.py @@ -0,0 +1,57 @@ +from django.test import TestCase + +from exact.administration.models import Product +from exact.images.models import Image, ImageSet, PyramidRegistration +from exact.users.models import Team, User +from util import registration + + +class _FakeSlide: + """A slide that needs registration; counts how often it was measured.""" + needs_registration = True + + def __init__(self): + self.reads = 0 + self.level_dimensions = [(2000, 2000), (1000, 1000)] + self.level_downsamples = [1.0, 2.0] + + +class PyramidRegistrationStorageTest(TestCase): + + def setUp(self): + user = User.objects.create_user('registration-test', password='x') + team = Team.objects.create(name='registration-test') + image_set = ImageSet.objects.create(name='set', team=team, creator=user, path='registration-test') + self.image = Image.objects.create(name='a.czi', filename='a.czi', image_set=image_set, checksum=b'x') + self.path = __file__ # any existing file: its size and mtime identify the "file version" + self.calls = [] + self._estimate = registration.estimate_level_offsets + registration.estimate_level_offsets = lambda osr: self.calls.append(1) or [(0.0, 0.0), (0.4, -0.2)] + + def tearDown(self): + registration.estimate_level_offsets = self._estimate + + def test_measured_once_then_read_from_the_database(self): + slide = _FakeSlide() + self.assertEqual(registration.level_offsets_for(self.image, slide, self.path), [(0.0, 0.0), (0.4, -0.2)]) + self.assertEqual(registration.level_offsets_for(self.image, slide, self.path), [(0.0, 0.0), (0.4, -0.2)]) + self.assertEqual(len(self.calls), 1) + row = PyramidRegistration.objects.get(image=self.image) + self.assertEqual(row.offsets, [[0.0, 0.0], [0.4, -0.2]]) + self.assertEqual(row.algorithm, registration.ALGORITHM) + + def test_measured_again_when_the_file_or_the_algorithm_changed(self): + slide = _FakeSlide() + registration.level_offsets_for(self.image, slide, self.path) + PyramidRegistration.objects.filter(image=self.image).update(file_mtime=1) # file replaced + registration.level_offsets_for(self.image, slide, self.path) + PyramidRegistration.objects.filter(image=self.image).update(algorithm=registration.ALGORITHM - 1) + registration.level_offsets_for(self.image, slide, self.path) + self.assertEqual(len(self.calls), 3) + + def test_formats_with_registered_levels_are_not_measured_or_stored(self): + slide = _FakeSlide() + slide.needs_registration = False + self.assertEqual(registration.level_offsets_for(self.image, slide, self.path), []) + self.assertEqual(self.calls, []) + self.assertFalse(PyramidRegistration.objects.exists()) diff --git a/exact/exact/images/views.py b/exact/exact/images/views.py index 0a0ddf0f..db5fff08 100644 --- a/exact/exact/images/views.py +++ b/exact/exact/images/views.py @@ -55,7 +55,7 @@ AnnotationType, Verification, LogImageAction from exact.tagger_messages.models import Message, TeamMessage, GlobalMessage -from util.slide_server import SlideCache, SlideFile, PILBytesIO, image_cache +from util.slide_server import SlideCache, SlideFile, PILBytesIO, image_cache, slide_for_image from util.cellvizio import ReadableCellVizioMKTDataset # just until data access is pip installable from plugins.pluginFinder import PluginFinder @@ -604,7 +604,7 @@ def view_image(request, image_id, z_dimension:int=1, frame:int=1): return HttpResponse(value, content_type='application/xml') file_path = os.path.join(settings.IMAGE_PATH, image.path(effective_z, frame)) - slide = image_cache.get(file_path, plane=plane) + slide = slide_for_image(image, file_path, plane=plane) value = slide.get_dzi("jpeg") if hasattr(cache, "delete_pattern"): @@ -818,7 +818,7 @@ def view_image_tile(request, image_id, z_dimension, frame, level, tile_path): file_path = os.path.join(settings.IMAGE_PATH, image.path(effective_z, frame)) try: - slide = image_cache.get(file_path, plane=plane) + slide = slide_for_image(image, file_path, plane=plane) n_frames = slide.nFrames if slide.nFrames is not None else image.frames tile = slide.get_tile(level, (col, row), frame=min(frame, n_frames - 1), plane=plane) diff --git a/exact/util/registration.py b/exact/util/registration.py new file mode 100644 index 00000000..38ba372f --- /dev/null +++ b/exact/util/registration.py @@ -0,0 +1,162 @@ +""" +Registration of the levels of an image pyramid. + +Some scanner formats store pyramid levels that are displaced against each other by a fraction of a pixel (seen: Zeiss +CZI read through slideio, 0.2-0.8 px between neighbouring levels at every position). The viewer blends tiles of two +Deep Zoom levels while it loads or zooms, so such a displacement shows as wobbling. + +The displacement is measured once per image: each level is compared with the level below it on up to three textured +spots (integer shift search plus parabolic refinement), the results are chained down to level 0 and a level is only +corrected if two spots agree. The offsets are kept in ``images.PyramidRegistration`` together with the size and +modification time of the file and applied when tiles are generated (``zDeepZoomGenerator.get_tile``). Level 0, the +full resolution, is the reference and is never shifted. + +The same table and algorithm are used by the Rust port of the server (exact-rust); keep ``ALGORITHM`` in sync. +""" +import logging +import os + +import numpy as np + +logger = logging.getLogger(__name__) + +# bump when the estimation changes: stored results of other versions are measured again +ALGORITHM = 2 + +WINDOW = 384 # side of the window under test, in pixels of the level +MARGIN = 8 +MAX_SHIFT = 3 +MAX_SPREAD = 0.15 # the two closest spots must agree within this (px) +MIN_OFFSET = 0.04 # smaller offsets are measurement noise + + +def _gray(image): + rgb = np.asarray(image.convert('RGB'), dtype=np.float32) + return rgb[..., 0] * 0.299 + rgb[..., 1] * 0.587 + rgb[..., 2] * 0.114 + + +def _reduce(array, factor): + h, w = array.shape[0] // factor * factor, array.shape[1] // factor * factor + return array[:h, :w].reshape(h // factor, factor, w // factor, factor).mean(axis=(1, 3)) + + +def _best_shift(level, reference): + """(dx, dy) with level(x + dx, y + dy) ~ reference(x, y), or None if the minimum is not interior.""" + h, w = reference.shape + core = reference[MARGIN:h - MARGIN, MARGIN:w - MARGIN] + n = 2 * MAX_SHIFT + 1 + grid = np.empty((n, n)) + for dy in range(-MAX_SHIFT, MAX_SHIFT + 1): + for dx in range(-MAX_SHIFT, MAX_SHIFT + 1): + patch = level[MARGIN + dy:h - MARGIN + dy, MARGIN + dx:w - MARGIN + dx] + grid[dy + MAX_SHIFT, dx + MAX_SHIFT] = float(((patch - core) ** 2).sum()) + iy, ix = np.unravel_index(np.argmin(grid), grid.shape) + by, bx = int(iy) - MAX_SHIFT, int(ix) - MAX_SHIFT + if abs(bx) >= MAX_SHIFT or abs(by) >= MAX_SHIFT: + return None + + def refine(left, mid, right): + d = left - 2 * mid + right + return 0.0 if d <= 0 else 0.5 * (left - right) / d + + return (bx + refine(grid[iy, ix - 1], grid[iy, ix], grid[iy, ix + 1]), + by + refine(grid[iy - 1, ix], grid[iy, ix], grid[iy + 1, ix])) + + +def _offsets_at(osr, centre): + """Chained offsets of every level at one spot; NaN where a level could not be measured.""" + dims, downs = osr.level_dimensions, osr.level_downsamples + out = [(0.0, 0.0)] + for k in range(1, len(dims)): + ratio = downs[k] / downs[k - 1] + r = round(ratio) + w, h = int(dims[k][0]), int(dims[k][1]) + side = min(WINDOW, w - 2 * MARGIN, h - 2 * MARGIN) + if abs(ratio - r) > 0.02 or r < 1 or side < 64: + out.append((0.0, 0.0)) + continue + a = int(min(max(centre[0] / downs[k] - side // 2, MARGIN), w - side - MARGIN)) + b = int(min(max(centre[1] / downs[k] - side // 2, MARGIN), h - side - MARGIN)) + location = (int(a * downs[k]), int(b * downs[k])) + try: + level = _gray(osr.read_region(location, k, (side, side))) + below = _gray(osr.read_region(location, k - 1, (side * r, side * r))) + except Exception: # noqa: BLE001 - a level that cannot be read is simply not corrected + out.append((float('nan'), float('nan'))) + continue + reference = _reduce(below, r) + measured = None if reference.std() < 6.0 else _best_shift(level, reference) + if measured is None: + out.append((float('nan'), float('nan'))) + continue + # the level below is itself displaced by out[k - 1] of its own pixels, i.e. that / r pixels of this level + out.append((measured[0] + out[k - 1][0] / r, measured[1] + out[k - 1][1] / r)) + return out + + +def estimate_level_offsets(osr): + """ + Offsets per level, in pixels of that level: sampling level k at (x + dx, y + dy) gives the content of level 0 at the + matching position. Empty if there is a single level or no texture to measure on. + """ + dims, downs = osr.level_dimensions, osr.level_downsamples + n = len(dims) + if n < 2: + return [] + last = n - 1 + thumb = _gray(osr.read_region((0, 0), last, (int(dims[last][0]), int(dims[last][1])))) + lh, lw = thumb.shape + ch, cw = max(lh // 4, 1), max(lw // 4, 1) + cells = sorted(((thumb[cy * ch:(cy + 1) * ch, cx * cw:(cx + 1) * cw].std(), cx, cy) + for cy in range(4) for cx in range(4)), reverse=True) + spots = [_offsets_at(osr, ((cx + 0.5) * cw * downs[last], (cy + 0.5) * ch * downs[last])) + for std, cx, cy in cells if std >= 8.0][:3] + if not spots: + return [] + offsets = [(0.0, 0.0)] # level 0 is the reference + for k in range(1, n): + values = [s[k] for s in spots if np.isfinite(s[k][0]) and np.isfinite(s[k][1])] + if not values: + offsets.append((0.0, 0.0)) + elif len(values) == 1: + offsets.append(tuple(0.0 if abs(v) < MIN_OFFSET else v for v in values[0])) + else: + # the closest pair of spots must agree, otherwise nothing is applied + best = min(((max(abs(a[0] - b[0]), abs(a[1] - b[1])), a, b) + for i, a in enumerate(values) for b in values[i + 1:]), key=lambda t: t[0]) + if best[0] <= MAX_SPREAD: + mean = ((best[1][0] + best[2][0]) / 2, (best[1][1] + best[2][1]) / 2) + offsets.append(tuple(0.0 if abs(v) < MIN_OFFSET else v for v in mean)) + else: + offsets.append((0.0, 0.0)) + return offsets + + +def level_offsets_for(image, osr, path): + """ + The registration offsets of an image: from the database if they were measured on this version of the file with this + version of the algorithm, otherwise measured now and stored. [] for formats whose levels are registered. + """ + if not getattr(osr, 'needs_registration', False) or image.depth != 1: + return [] + from exact.images.models import PyramidRegistration + + stat = os.stat(path) + size, mtime = stat.st_size, int(stat.st_mtime) + stored = PyramidRegistration.objects.filter(image=image).first() + if stored and stored.algorithm == ALGORITHM and stored.file_size == size and stored.file_mtime == mtime: + return [tuple(o) for o in stored.offsets] + offsets = estimate_level_offsets(osr) + if offsets: + offsets[0] = (0.0, 0.0) + logger.info('pyramid level offsets of image %s: %s', image.pk, ['(%.2f, %.2f)' % o for o in offsets]) + PyramidRegistration.objects.update_or_create(image=image, defaults={ + 'offsets': [list(o) for o in offsets], 'algorithm': ALGORITHM, 'file_size': size, 'file_mtime': mtime}) + return offsets + + +def attach_level_offsets(slide, image, path): + """Sets ``slide.level_offsets`` (once per cached slide) and returns the slide.""" + if getattr(slide, 'level_offsets', None) is None: + slide.level_offsets = level_offsets_for(image, slide._osr, path) + return slide diff --git a/exact/util/slide_server.py b/exact/util/slide_server.py index 6aac6bae..5c0cc6f9 100644 --- a/exact/util/slide_server.py +++ b/exact/util/slide_server.py @@ -50,20 +50,28 @@ def get_tile(self, level, address, frame=0, plane=PlaneType.AXIAL): _, z_size = self._get_tile_info(level, address) # validates level and address slide_level = self._slide_from_dz_level[level] z_w, z_h = self._z_dimensions[level] + # A slide level covers ``level 0 size / downsample`` pixels of the image, which is not necessarily its stored size + # (levels made by decimation with rounding down cover slightly less): the scale from the slide level to the + # Deep Zoom level is derived from the downsample, not from the stored size. + downsample = self._l0_l_downsamples[slide_level] l_w, l_h = self._l_dimensions[slide_level] - scale_x, scale_y = z_w / l_w, z_h / l_h # slide level pixels -> Deep Zoom pixels + scale_x = z_w / (self._l_dimensions[0][0] / downsample) # slide level pixels -> Deep Zoom pixels + scale_y = z_h / (self._l_dimensions[0][1] / downsample) + # registration offsets of the pyramid levels (util.registration); level 0 is the reference and never shifted + offsets = getattr(self, 'level_offsets', None) + off_x, off_y = offsets[slide_level] if offsets and 0 < slide_level < len(offsets) else (0.0, 0.0) col, row = address overlap = self._z_overlap # tile rectangle in Deep Zoom pixels, converted to (fractional) slide level pixels - box_x = (self._z_t_downsample * col - (overlap if col else 0)) / scale_x - box_y = (self._z_t_downsample * row - (overlap if row else 0)) / scale_y + box_x = max(0.0, (self._z_t_downsample * col - (overlap if col else 0)) / scale_x + off_x) + box_y = max(0.0, (self._z_t_downsample * row - (overlap if row else 0)) / scale_y + off_y) box_w, box_h = z_size[0] / scale_x, z_size[1] / scale_y # integer region that covers the box x0, y0 = max(0, math.floor(box_x)), max(0, math.floor(box_y)) x1, y1 = min(l_w, math.ceil(box_x + box_w)), min(l_h, math.ceil(box_y + box_h)) r_w, r_h = max(1, x1 - x0), max(1, y1 - y0) - l0_location = (int(self._l0_l_downsamples[slide_level] * x0) + self._l0_offset[0], - int(self._l0_l_downsamples[slide_level] * y0) + self._l0_offset[1]) + l0_location = (int(downsample * x0) + self._l0_offset[0], + int(downsample * y0) + self._l0_offset[1]) args = (l0_location, slide_level, (r_w, r_h)) # Read tile — pass plane only to handlers that declare MPR support @@ -807,6 +815,12 @@ def stats(self): image_cache = SlideCache(cache_size=20) +def slide_for_image(image, path, plane=PlaneType.AXIAL): + """``image_cache.get`` plus the pyramid registration offsets of the image (see util.registration).""" + from util.registration import attach_level_offsets + return attach_level_offsets(image_cache.get(path, plane=plane), image, path) + + class SlideFile(object): def __init__(self, relpath): self.name = os.path.basename(relpath) diff --git a/exact/util/slideio.py b/exact/util/slideio.py index bfd18b3d..f326c04a 100644 --- a/exact/util/slideio.py +++ b/exact/util/slideio.py @@ -65,9 +65,29 @@ def __init__(self,filename): def seriesInstanceUID(self) -> str: return '' + @staticmethod + def _exact_downsample(downsample): + """ + The integer downsample if ``downsample`` is within 1 % of one, otherwise None. + + slideio reports ``1 / scale`` for a level, and the scale is derived from the (rounded) level size, so a 16:1 + level is reported as 16.0007. The levels themselves are exact decimation grids: a level pixel covers exactly + ``round(downsample)`` pixels of level 0 (reading a level at positions that are multiples of it matches level 0 + reduced by exactly that factor to 0.00 px; any other position is displaced by up to a level pixel). + """ + exact = round(downsample) + if exact >= 1 and abs(downsample - exact) <= 0.01 * exact: + return float(exact) + return None + @property def level_downsamples(self): - return [1/self.scene.get_zoom_level_info(k).scale for k in range(self.scene.num_zoom_levels)] + downsamples = [1/self.scene.get_zoom_level_info(k).scale for k in range(self.scene.num_zoom_levels)] + return [self._exact_downsample(d) or d for d in downsamples] + + # Zeiss CZI pyramids store levels that are displaced against each other by a fraction of a pixel; the offsets + # are measured once per image (util.registration) and applied when tiles are generated + needs_registration = True @property def level_dimensions(self): @@ -211,8 +231,16 @@ def get_thumbnail(self, size): def read_region(self, location: tuple, level:int, size:tuple, frame:int=0): - ds = self.scene.get_zoom_level_info(level).scale - img = self.scene.read_block(rect=[*location, int(size[0]/ds), int(size[1]/ds)], size=size, slices=(frame, frame+1)) + downsample = self._exact_downsample(1/self.scene.get_zoom_level_info(level).scale) + if downsample: + # slideio snaps a region to its level grid: reading at a position that is not a multiple of the downsample + # displaces the result by up to a level pixel. Snap the position and request exactly size * downsample. + location = tuple(int(round(v / downsample) * downsample) for v in location) + rect = [*location, int(round(size[0] * downsample)), int(round(size[1] * downsample))] + else: + ds = self.scene.get_zoom_level_info(level).scale + rect = [*location, int(size[0]/ds), int(size[1]/ds)] + img = self.scene.read_block(rect=rect, size=size, slices=(frame, frame+1)) img_4ch = np.zeros([size[1], size[0],4], dtype=np.uint8) img_4ch[:,:,3] = 255 if img.ndim == 2: diff --git a/exact/util/test_registration.py b/exact/util/test_registration.py new file mode 100644 index 00000000..16cf3086 --- /dev/null +++ b/exact/util/test_registration.py @@ -0,0 +1,94 @@ +"""Estimation of the displacement between pyramid levels (util.registration) and its use by the tile generator.""" +import unittest + +import numpy as np +import openslide +from PIL import Image + +from util import registration +from util.slide_server import zDeepZoomGenerator + + +def _texture(width, height, seed=3): + rng = np.random.RandomState(seed) + y, x = np.mgrid[0:height, 0:width].astype(np.float64) + img = np.full((height, width), 128.0) + for _ in range(90): + cx, cy, s, a = rng.uniform(0, width), rng.uniform(0, height), rng.uniform(3, 35), rng.uniform(-110, 110) + img += a * np.exp(-(((x - cx) ** 2 + (y - cy) ** 2) / (2 * s * s))) + return np.clip(img, 0, 255) + + +def _shifted_half(big, sx, sy, size): + """The ``size`` x ``size`` window of ``big`` at (sx, sy) (fractional), halved: exactly 2:1, only displaced.""" + image = Image.fromarray(big.astype(np.float32), mode='F') + return np.asarray(image.resize((size // 2, size // 2), Image.BOX, box=(sx, sy, sx + size, sy + size))) + + +class _Pyramid: + """Minimal slide object with levels of exactly 2^k, as the handlers in slide_server provide.""" + + def __init__(self, levels): + self._levels = [Image.fromarray(np.clip(l, 0, 255).astype(np.uint8)).convert('RGBA') for l in levels] + self.level_count = len(levels) + self.level_dimensions = [im.size for im in self._levels] + self.level_downsamples = [2.0 ** k for k in range(len(levels))] + self.dimensions = self.level_dimensions[0] + self.properties = {} + self.needs_registration = True + + def get_best_level_for_downsample(self, downsample): + return max(k for k, d in enumerate(self.level_downsamples) if d <= downsample + 1e-9) + + def read_region(self, location, level, size, frame=0): + ds = int(self.level_downsamples[level]) + x, y = location[0] // ds, location[1] // ds + tile = Image.new('RGBA', size, (0, 0, 0, 0)) + tile.paste(self._levels[level], (-x, -y)) + return tile + + +N = 1024 + + +class RegistrationEstimateTest(unittest.TestCase): + + def test_finds_a_displaced_level(self): + big = _texture(N + 16, N + 16) + level0 = big[8:8 + N, 8:8 + N] + # level 1 shows the content displaced by (+1.2, -0.8) source px = (+0.6, -0.4) of its own pixels + level1 = _shifted_half(big, 8 + 1.2, 8 - 0.8, N) + offsets = registration.estimate_level_offsets(_Pyramid([level0, level1])) + self.assertEqual(offsets[0], (0.0, 0.0), 'level 0 is the reference and never shifted') + # sampling level 1 at (x + dx) must give the reference: dx = -0.6, dy = +0.4 + self.assertAlmostEqual(offsets[1][0], -0.6, delta=0.15) + self.assertAlmostEqual(offsets[1][1], 0.4, delta=0.15) + + def test_a_registered_pyramid_gets_no_offsets(self): + big = _texture(N + 16, N + 16) + offsets = registration.estimate_level_offsets(_Pyramid([big[8:8 + N, 8:8 + N], _shifted_half(big, 8, 8, N)])) + self.assertEqual(offsets, [(0.0, 0.0), (0.0, 0.0)], 'noise below the threshold is not applied') + + def test_a_single_level_has_nothing_to_register(self): + self.assertEqual(registration.estimate_level_offsets(_Pyramid([_texture(512, 512)])), []) + + +class GeneratorUsesOffsetsTest(unittest.TestCase): + + def _tiles(self, offsets): + big = _texture(N + 16, N + 16) + slide = _Pyramid([big[8:8 + N, 8:8 + N], _shifted_half(big, 8, 8, N)]) + generator = zDeepZoomGenerator(slide) + generator.level_offsets = offsets + top = generator.level_count - 1 + return generator.get_tile(top, (1, 1)), generator.get_tile(top - 1, (1, 1)) + + def test_level_zero_is_not_shifted_and_other_levels_are(self): + plain0, plain1 = self._tiles(None) + shifted0, shifted1 = self._tiles([(7.0, -5.0), (0.6, 0.3)]) # an absurd offset for level 0 must be ignored + self.assertEqual(plain0.tobytes(), shifted0.tobytes()) + self.assertNotEqual(plain1.tobytes(), shifted1.tobytes()) + + +if __name__ == '__main__': + unittest.main()