From 6cc40b2c8309919bc007a3c80dc4801cf74808a9 Mon Sep 17 00:00:00 2001 From: Valentin Boussot Date: Thu, 27 Aug 2026 18:59:25 +0200 Subject: [PATCH] feat(geometry): bound a displacement by where it reaches, not how far TransformBound carried a symmetric residual, so a map was bounded by affine(p) +/- r. That is the right shape for a wiggle and the wrong one for a field solved between two frames, which carries the offset between them in every voxel: priced as a radius such a field opens a region's window both ways by twice its largest value, where the truth is a window of the same size sitting that far along. It becomes an interval. WorldBox.extended moves each end by its own end of it, DisplacementStage keeps min and max where it kept the larger magnitude -- the same two reductions, measured 79 ms against 80 on a 184 M-value window -- and after() folds through the outer affine by splitting the matrix into its non-negative and non-positive parts, because |A| @ r loses which end a sign flip sends where. Measured on an ExaSPIM field whose z runs [-28.1, -22.2] mm on a volume 20.6 mm thick: as a radius every region pulls the whole volume and the fold refuses, 23.57 GiB held against a 19.01 GiB budget; as an interval a 24-row region pulls 175 source rows of 514. residual_xyz stays as the symmetric envelope, which is what sup |v| already was, so coverage and anything else reading one number per axis is unchanged. --- konfai/data/geometry.py | 105 ++++++++++++++++++++++++++++-------- konfai/utils/ome_zarr.py | 91 ++++++++++++++++++++++++++++--- tests/unit/test_geometry.py | 57 +++++++++++++++++++- tests/unit/test_warp.py | 23 +++++--- 4 files changed, 239 insertions(+), 37 deletions(-) diff --git a/konfai/data/geometry.py b/konfai/data/geometry.py index 02c67e9f..cb48747b 100644 --- a/konfai/data/geometry.py +++ b/konfai/data/geometry.py @@ -120,6 +120,18 @@ def grown(self, radius_xyz: np.ndarray | float) -> WorldBox: radius = np.broadcast_to(np.asarray(radius_xyz, dtype=np.float64), self.low_xyz.shape) return WorldBox(self.low_xyz - radius, self.high_xyz + radius) + def extended(self, low_xyz: np.ndarray, high_xyz: np.ndarray) -> WorldBox: + """This box plus a per-component interval: each end moved by its own end of it. + + The asymmetric form of :meth:`grown`, and the one a signed displacement bound needs. An + interval that does not straddle zero MOVES the box instead of widening it, which is the + difference between a field's reach and twice its largest value. + """ + return WorldBox( + self.low_xyz + np.asarray(low_xyz, dtype=np.float64), + self.high_xyz + np.asarray(high_xyz, dtype=np.float64), + ) + def image_under(self, affine: AffineMap) -> WorldBox: """The axis-aligned hull of this box's image under ``affine``. @@ -361,37 +373,74 @@ def sub_grid(self, region_zyx: tuple[slice, ...]) -> Grid: @dataclass(frozen=True) class TransformBound: - """What a stored transform is guaranteed to do: an exact affine part and a bounded residual. - - ``T(p)`` lies in ``affine(p) ± residual_xyz`` for every ``p``, per world component. For a - linear transform the residual is zero and the statement is exact; for a BSpline it is the - sup-norm of the coefficients (non-negative basis functions summing to one make every - displacement a convex combination of them); for a dense field it is the recorded or declared - per-component bound. The affine part is read structurally off the transform, never probed: - a probe measures a local gradient and extrapolates it, which under-bounds (measured). + """What a stored transform is guaranteed to do: an exact affine part and a bounded interval. + + ``T(p)`` lies in ``affine(p) + [low_xyz, high_xyz]`` for every ``p``, per world component. For a + linear transform the interval is empty and the statement is exact; for a BSpline it is the range + of the coefficients (non-negative basis functions summing to one make every displacement a + convex combination of them, so it lies between their smallest and largest); for a dense field it + is the range of its values. The affine part is read structurally off the transform, never + probed: a probe measures a local gradient and extrapolates it, which under-bounds (measured). + + SIGNED, NOT A RADIUS. A displacement field solved between two frames carries the offset between + them in its values, and an interval that does not straddle zero MOVES a region's window rather + than widening it. Measured on an ExaSPIM field whose z component runs [-28.1, -22.2] mm on a + volume 20.6 mm thick: as a radius it reaches 28.1 mm either way, so every region pulls the whole + volume and the fold refuses (23.57 GiB held against a 19.01 GiB budget); as an interval it + reaches 5.9 mm, and a 24-row region pulls 175 source rows of 514. The same two reductions + produce either (:attr:`DisplacementStage.range_xyz`), so the tighter one is free. """ affine: AffineMap - residual_xyz: np.ndarray + low_xyz: np.ndarray + high_xyz: np.ndarray @staticmethod def exact(affine: AffineMap) -> TransformBound: - return TransformBound(affine, np.zeros(affine.rank)) + return TransformBound(affine, np.zeros(affine.rank), np.zeros(affine.rank)) + + @staticmethod + def interval(low_xyz: np.ndarray, high_xyz: np.ndarray) -> TransformBound: + """A pure displacement whose value lies in ``[low_xyz, high_xyz]`` per component.""" + low = np.asarray(low_xyz, dtype=np.float64) + return TransformBound(AffineMap.identity(int(low.size)), low, np.asarray(high_xyz, dtype=np.float64)) @staticmethod def shift(residual_xyz: np.ndarray) -> TransformBound: - return TransformBound(AffineMap.identity(int(residual_xyz.size)), np.asarray(residual_xyz, dtype=np.float64)) + """A pure displacement bounded in magnitude only, ``± residual_xyz``. + + For a caller that knows a radius and not a range. Anything that can state both ends should + say so with :meth:`interval`: this one is twice as wide wherever the range is one-sided. + """ + radius = np.asarray(residual_xyz, dtype=np.float64) + return TransformBound.interval(-radius, radius) + + @property + def residual_xyz(self) -> np.ndarray: + """The symmetric envelope of the interval, for a caller that wants one number per axis.""" + return np.maximum(np.abs(self.low_xyz), np.abs(self.high_xyz)) def after(self, inner: TransformBound) -> TransformBound: - """The bound of ``self(inner(p))``: interval arithmetic through the outer affine.""" + """The bound of ``self(inner(p))``: interval arithmetic through the outer affine. + + NOT ``|A| @ residual``. That is right for an interval centred on zero and wrong for one that + is not: a negative entry of ``A`` sends the inner interval's low end to the outer's high, + and taking absolute values first loses which end went where -- so a rotation folded onto a + one-sided field would be bounded by a box that does not contain it. Splitting the matrix + into its non-negative and non-positive parts is the same arithmetic written to hold either + way, and it reduces to ``|A| @ r`` when ``low = -high``. + """ + matrix = self.affine.matrix + rise, fall = np.maximum(matrix, 0.0), np.minimum(matrix, 0.0) return TransformBound( inner.affine.then(self.affine), - np.abs(self.affine.matrix) @ inner.residual_xyz + self.residual_xyz, + rise @ inner.low_xyz + fall @ inner.high_xyz + self.low_xyz, + rise @ inner.high_xyz + fall @ inner.low_xyz + self.high_xyz, ) def map_box(self, box: WorldBox) -> WorldBox: """Where the image of ``box`` is guaranteed to lie.""" - return box.image_under(self.affine).grown(self.residual_xyz) + return box.image_under(self.affine).extended(self.low_xyz, self.high_xyz) @dataclass(frozen=True) @@ -467,28 +516,40 @@ def tensor_by_voxel(self, device: torch.device, dtype: torch.dtype) -> torch.Ten def __getstate__(self) -> dict: state = dict(self.__dict__) state["_tensors"] = {} + state.pop("range_xyz", None) state.pop("bound_xyz", None) return state @cached_property - def bound_xyz(self) -> np.ndarray: - """``sup |values|`` per component: one pass over the field, kept for the stage's life. + def range_xyz(self) -> tuple[np.ndarray, np.ndarray]: + """``(min, max)`` per component: one pass over the field, kept for the stage's life. Every pull map asks for it (per patch, per plan block, per pushed slab), so it is kept: a thousand patches recomputing one constant 3-vector cost 73 s on a 3x160x256x256 field. ``cached_property`` writes through ``__dict__``, which a frozen dataclass allows; the entry is dropped from the pickle beside ``_tensors``. - Two reductions and no temporary: ``sup |v| = max(max v, -min v)`` is exact for real values, - where ``np.abs(...).max()`` first writes a values-sized copy and then walks it again. On the - 31 M-voxel field the difference is 72 ms and 250 MiB; on a native ExaSPIM field window - (3 x 141 x 1331 x 1775) the copy alone is 4 GiB, and the fold spent 26.4 s of its 179 s here. + Two reductions and no temporary, which is what the pair costs and what its magnitude cost + before it: ``np.abs(...).max()`` would first write a values-sized copy and then walk it + again. On the 31 M-voxel field that copy is 72 ms and 250 MiB; on a native ExaSPIM field + window (3 x 141 x 1331 x 1775) it is 4 GiB, and the fold spent 26.4 s of its 179 s here. + Keeping both ends instead of the larger magnitude measured 79 ms against 80 on a 184 + M-value window: the range is free, and it is the one a region can be sized from. """ flat = self.values.reshape(self.values.shape[0], -1) - return np.maximum(flat.max(axis=1), -flat.min(axis=1)) + return flat.min(axis=1), flat.max(axis=1) + + @cached_property + def bound_xyz(self) -> np.ndarray: + """``sup |values|`` per component, the range's symmetric envelope.""" + low, high = self.range_xyz + return np.maximum(np.abs(low), np.abs(high)) def bound(self) -> TransformBound: - return TransformBound.shift(self.bound_xyz) + # The interval is clamped to include zero: the stage applies NO displacement outside its + # grid, so a target region past the field's edge still needs its identity-mapped samples. + low, high = self.range_xyz + return TransformBound.interval(np.minimum(low, 0.0), np.maximum(high, 0.0)) #: A decoded stored transform: stages in APPLICATION order (first applied first). SimpleITK's diff --git a/konfai/utils/ome_zarr.py b/konfai/utils/ome_zarr.py index 6d162bc5..94363cf4 100644 --- a/konfai/utils/ome_zarr.py +++ b/konfai/utils/ome_zarr.py @@ -39,6 +39,7 @@ import itertools import operator import shutil +import tempfile import threading from collections import OrderedDict from collections.abc import Sequence @@ -48,6 +49,7 @@ import numpy as np +from konfai.utils import uri from konfai.utils.errors import DatasetManagerError from konfai.utils.runtime import map_over_rank_pool @@ -177,6 +179,22 @@ def _konfai_attributes(store_path: str) -> dict[str, Any]: return {} +def _from_ngff_zarr(store_path: str | Path) -> Any: + """ngff-zarr's multiscales for ``store_path``. A remote root goes in as a key-to-bytes mapping + over its own filesystem (ngff-zarr >= 0.44 reads a remote string as a local path) and as its URL + when ngff-zarr refuses the mapping, which older releases resolve themselves.""" + if not uri.is_uri(store_path): + return ngff_zarr.from_ngff_zarr(str(store_path)) + # Through uri.filesystem, so a missing fsspec backend or configuration is the structured + # DatasetManagerError, never a raw dependency error. + filesystem = uri.filesystem(store_path) + _, target = uri.split_scheme(str(store_path)) + try: + return ngff_zarr.from_ngff_zarr(filesystem.get_mapper(target)) + except ValueError: + return ngff_zarr.from_ngff_zarr(str(store_path)) + + @lru_cache(maxsize=8) def _load_image(store_path: str, level: int) -> Any: """Return the ``NgffImage`` for ``level`` of an OME-Zarr store, memoised per (store, level). @@ -194,7 +212,7 @@ def _load_image(store_path: str, level: int) -> Any: """ _require_ngff_zarr() try: - multiscales = ngff_zarr.from_ngff_zarr(str(store_path)) + multiscales = _from_ngff_zarr(store_path) except (KeyError, IndexError, OSError, TypeError, ValueError) as exc: raise DatasetManagerError( f"Cannot open OME-Zarr store '{store_path}' (level {level}).", @@ -258,7 +276,7 @@ def is_displacement_field(store_path: str | Path) -> bool: if not _NGFF_ZARR_AVAILABLE: return False try: - metadata = ngff_zarr.from_ngff_zarr(str(store_path)).metadata + metadata = _from_ngff_zarr(store_path).metadata except Exception: # "Not a displacement field" is the only answer this owes: it is asked purely to decide HOW to # read an entry, and an absent or unreadable store is not one either. @@ -608,7 +626,7 @@ def place(coords: tuple) -> None: def _level_path(store_path: str, level: int) -> str | None: """The zarr path of one level, from the store's multiscales metadata, memoised beside the image.""" try: - datasets = ngff_zarr.from_ngff_zarr(store_path).metadata.datasets + datasets = _from_ngff_zarr(store_path).metadata.datasets return str(datasets[level if len(datasets) > 1 else 0].path) except Exception: return None @@ -624,7 +642,32 @@ def _read_level_window(store_path: str, level: int, image: Any, index: tuple) -> array = _level_array(store_path, level_path) if tuple(array.shape) == tuple(image.data.shape): return _read_chunked(store_path, level_path, array, index) - return np.asarray(image.data[index]) + return _lazy_window(image.data, index) + + +def _lazy_window(data: Any, index: tuple) -> np.ndarray: + """The plain lazy read. A lazy array that refuses a stepped slice (ngff-zarr >= 0.44 wraps the + level in an adapter that takes unit steps only) is read over the unit-step span and stepped here.""" + try: + return np.asarray(data[index]) + except NotImplementedError: + # Each slice normalized against the shape, then its ascending unit-step span; a negative + # step reads that span backwards from its own end, which lands on the indices the original + # slice named. + bounds = tuple( + slice(*item.indices(size)) if isinstance(item, slice) else item + for item, size in zip(index, data.shape, strict=True) + ) + span = tuple( + slice(item.stop + 1, item.start + 1) + if isinstance(item, slice) and item.step < 0 + else slice(item.start, item.stop) + if isinstance(item, slice) + else item + for item in bounds + ) + steps = tuple(slice(None, None, item.step) for item in bounds if isinstance(item, slice)) + return np.asarray(data[span])[steps] def _store_index( @@ -813,6 +856,11 @@ def write_ome_zarr( exist only from 0.6, so a caller passing both could only ever pass them consistently: an invariant worth removing rather than documenting. """ + if scale_factors and uri.is_uri(store_path): + raise DatasetManagerError( + f"Cannot append pyramid levels to the remote store '{store_path}'.", + "Levels are derived in place through local paths; write the store locally and upload it.", + ) array_data = np.asarray(data) # The one write path: the store described and created empty (ngff-zarr's metadata, the # caller's chunking), filled by zarr itself, its levels grafted beside level 0. Handing @@ -876,6 +924,30 @@ def _type_component_axis(multiscales: Any, axis_type: str) -> None: axis.type = axis_type +def _write_skeleton(store_path: str | Path, multiscales: Any, version: str) -> None: + """ngff-zarr's metadata for the store, written in place. ngff-zarr (>= 0.44) writes local + directories only, so a remote root gets the skeleton written locally and uploaded through the + root's own filesystem: a few bytes of metadata, before the array is created underneath it.""" + if not uri.is_uri(store_path): + ngff_zarr.to_ngff_zarr(str(store_path), multiscales, overwrite=True, version=version) + return + filesystem = uri.filesystem(store_path) + _, target = uri.split_scheme(str(store_path)) + if "/" not in target.strip("/"): + raise DatasetManagerError( + f"Refusing to create the store at the filesystem root '{store_path}'.", + "Name a key under the root (e.g. '.../dataset.ome.zarr'): creating a store replaces what its path holds.", + ) + if filesystem.exists(target): + filesystem.rm(target, recursive=True) + filesystem.makedirs(target, exist_ok=True) + with tempfile.TemporaryDirectory() as scratch: + local = Path(scratch) / "skeleton" + ngff_zarr.to_ngff_zarr(str(local), multiscales, overwrite=True, version=version) + for file in sorted(path for path in local.rglob("*") if path.is_file()): + filesystem.put_file(str(file), uri.join(target, file.relative_to(local).as_posix())) + + def create_ome_zarr_store( store_path: str | Path, shape: Sequence[int], @@ -934,7 +1006,7 @@ def create_ome_zarr_store( _type_component_axis(multiscales, _DISPLACEMENT_AXIS_TYPE) version = _RFC5_VERSION # version is explicit because to_ngff_zarr defaults to 0.5, which zarr-python 2 cannot write. - ngff_zarr.to_ngff_zarr(str(store_path), multiscales, overwrite=True, version=version) + _write_skeleton(store_path, multiscales, version) # The level-0 key comes from the metadata rather than a literal: ngff-zarr builds it from the # image name, so "scale0/image" is its convention to change, not ours to hardcode. @@ -1077,13 +1149,18 @@ def append_ome_zarr_levels( KonfAI attribute sidecar is untouched, being a key beside theirs; a displacement field keeps its typed component axis through the same call that types it at creation. """ + if uri.is_uri(store_path): + raise DatasetManagerError( + f"Cannot append pyramid levels to the remote store '{store_path}'.", + "Levels are derived in place through local paths; write the store locally and upload it.", + ) _require_ngff_zarr() if not scale_factors: return store = Path(store_path) clear_ome_zarr_cache(store) field = is_displacement_field(store) - base = ngff_zarr.from_ngff_zarr(str(store)).images[0] + base = _from_ngff_zarr(store).images[0] stored_chunks = tuple(int(size) for size in base.data.chunksize) if downsample_method in (None, "ITKWASM_BIN_SHRINK"): multiscales = _bin_shrink_multiscales(base, _level_zero_scale_factors(scale_factors), stored_chunks) @@ -1152,7 +1229,7 @@ def get_ome_zarr_info(store_path: str | Path, level: int = 0) -> dict[str, Any]: image = _load_image(str(store_path), level) dims = [str(axis).lower() for axis in image.dims] try: - n_levels = len(ngff_zarr.from_ngff_zarr(str(store_path)).images) + n_levels = len(_from_ngff_zarr(store_path).images) except (OSError, TypeError, ValueError): n_levels = 1 return { diff --git a/tests/unit/test_geometry.py b/tests/unit/test_geometry.py index f721ad55..c73b350d 100644 --- a/tests/unit/test_geometry.py +++ b/tests/unit/test_geometry.py @@ -174,7 +174,7 @@ def test_map_box_contains_the_true_image_of_a_nonlinear_map(self): rng = np.random.RandomState(3) affine = AffineMap(np.eye(3) + rng.randn(3, 3) * 0.1, rng.randn(3) * 10.0) residual = np.array([4.0, 2.0, 1.0]) - bound = TransformBound(affine, residual) + bound = TransformBound(affine, -residual, residual) box = WorldBox(np.array([-20.0, -10.0, 0.0]), np.array([15.0, 25.0, 30.0])) mapped = bound.map_box(box) points = rng.uniform(box.low_xyz, box.high_xyz, size=(500, 3)) @@ -182,6 +182,39 @@ def test_map_box_contains_the_true_image_of_a_nonlinear_map(self): images = affine.apply(points) + rng.uniform(-1.0, 1.0, size=(500, 3)) * residual assert np.all(images >= mapped.low_xyz - 1e-9) and np.all(images <= mapped.high_xyz + 1e-9) + def test_a_one_sided_interval_moves_the_box_instead_of_widening_it(self): + # The whole reason the bound is a range and not a radius. A field that displaces every point + # by the same 30 mm is an offset, not a spread: the region it reads sits 30 mm away and is + # no larger. Priced as a radius it would be 60 mm wider AND still centred where it started, + # which on a volume thinner than 60 mm is every voxel there is. + box = WorldBox(np.zeros(3), np.array([10.0, 10.0, 10.0])) + offset = TransformBound.interval(np.full(3, -30.0), np.full(3, -30.0)) + moved = offset.map_box(box) + np.testing.assert_array_equal(moved.low_xyz, [-30.0, -30.0, -30.0]) + np.testing.assert_array_equal(moved.high_xyz, [-20.0, -20.0, -20.0]) + # The radius spelling of the same field is the one that cannot say that. + radius = TransformBound.shift(np.full(3, 30.0)).map_box(box) + assert np.all(radius.high_xyz - radius.low_xyz > moved.high_xyz - moved.low_xyz) + + def test_after_carries_a_one_sided_interval_through_a_sign_flip(self): + # |A| @ residual is right for an interval centred on zero and wrong for one that is not: a + # negative entry sends the inner low end to the outer high end. Here the outer affine + # mirrors x, so a field that only ever pushes -x must come out only ever pushing +x. The + # symmetric arithmetic would answer [-4, +4] on every axis: containing, but four times the + # width, and it hides which side the map actually reaches. + mirror = TransformBound.exact(AffineMap(np.diag([-1.0, 1.0, 1.0]), np.zeros(3))) + pushes = TransformBound.interval(np.array([-4.0, 0.0, 0.0]), np.array([-2.0, 0.0, 0.0])) + folded = mirror.after(pushes) + np.testing.assert_array_equal(folded.low_xyz, [2.0, 0.0, 0.0]) + np.testing.assert_array_equal(folded.high_xyz, [4.0, 0.0, 0.0]) + # And it still contains the map it bounds, which is the property the sign trick must keep. + rng = np.random.RandomState(11) + points = rng.uniform(-50.0, 50.0, size=(400, 3)) + displaced = points + np.stack([rng.uniform(-4.0, -2.0, 400), np.zeros(400), np.zeros(400)], axis=-1) + images = mirror.affine.apply(displaced) + box = folded.map_box(WorldBox(points.min(axis=0), points.max(axis=0))) + assert np.all(images >= box.low_xyz - 1e-9) and np.all(images <= box.high_xyz + 1e-9) + class TestDisplacementStageBound: def test_the_bound_is_one_pass_per_stage_and_stays_out_of_the_pickle(self): @@ -207,6 +240,26 @@ def test_the_bound_is_one_pass_per_stage_and_stays_out_of_the_pickle(self): skewed_stage = DisplacementStage(_grid(), skewed, 1) np.testing.assert_array_equal(skewed_stage.bound_xyz, np.abs(skewed).reshape(3, -1).max(axis=1)) # The cache is per instance, not per pickle: a rank rebuilds it from the values it receives. - assert "bound_xyz" not in stage.__getstate__() + assert "bound_xyz" not in stage.__getstate__() and "range_xyz" not in stage.__getstate__() again = pickle.loads(pickle.dumps(stage)) np.testing.assert_array_equal(again.bound_xyz, first) + + def test_the_stage_bounds_a_one_sided_field_by_where_it_reaches(self): + from konfai.data.geometry import DisplacementStage + + # A field that only ever pushes one way -- which is what a registration between two frames + # writes, the offset between them baked into every voxel. The stage must say where it + # reaches, not how far: as a magnitude this field is worth 12 either side, and it never + # sends a point anywhere but 8 to 12 units below where it started. + values = np.random.RandomState(7).uniform(-12.0, -8.0, (3, 5, 6, 7)) + stage = DisplacementStage(_grid(), values, 1) + low, high = stage.range_xyz + np.testing.assert_allclose(low, values.reshape(3, -1).min(axis=1)) + np.testing.assert_allclose(high, values.reshape(3, -1).max(axis=1)) + bound = stage.bound() + # The BOUND still includes zero: outside its grid the stage displaces nothing, so a region + # past the field's edge keeps its identity-mapped samples in the source window. + np.testing.assert_array_equal(bound.low_xyz, low) + np.testing.assert_array_equal(bound.high_xyz, np.zeros(3)) + # The envelope is still there for whoever wants one number, and it is the old sup |v|. + np.testing.assert_array_equal(stage.bound_xyz, np.abs(values).reshape(3, -1).max(axis=1)) diff --git a/tests/unit/test_warp.py b/tests/unit/test_warp.py index 3b8afc92..5632bc32 100644 --- a/tests/unit/test_warp.py +++ b/tests/unit/test_warp.py @@ -88,13 +88,18 @@ def _recorded(warp: Resample, attribute: Attribute | None = None, shape: tuple[i return warp -def test_the_source_region_is_the_target_grown_by_the_field_reach(tmp_path: Path) -> None: - """A warp is a regrid onto the case's own grid, and its window is the field's reach in voxels. +def test_the_source_region_is_the_target_moved_by_the_field_range(tmp_path: Path) -> None: + """A warp is a regrid onto the case's own grid, and its window is where the field points. Spacing is (x=2, y=1, z=1), so in array order (z, y, x) 4 um of displacement is 4, 4 and 2 voxels: plus the one voxel the linear taps reach. Declared as REGRID and not HALO because the window is derived from the case's GEOMETRY: see the oblique case below, which a per-axis halo cannot express at all. + + MOVED, not grown. This field displaces every point by the same +4 um, so the region reads a + window of its own size sitting 4 um along, not one twice as wide centred where it started. The + two spellings differ by the whole displacement on each face, which is nothing on a wiggle and + everything on a field that carries the offset between two frames. """ _source, _fields, _volume = _fixture(tmp_path, shift_um=(4.0, 4.0, 4.0)) warp = _recorded(Resample(field=f"{tmp_path / 'dvf'}:h5", field_group="DVF")) @@ -104,14 +109,18 @@ def test_the_source_region_is_the_target_grown_by_the_field_reach(tmp_path: Path window = warp.measured_region_source("CASE_000", target, [10, 12, 14], _attributes()) # The rule, written out: the region's OUTER faces (start - 0.5 .. stop - 0.5) in world units, - # grown by the field's 4 um, back to indices, floor/ceil, one voxel of margin for the taps. + # extended by the field's range CLAMPED TO INCLUDE the identity -- here [0, +4] um, so the + # window keeps its own faces and reaches one-sidedly where the field points -- back to indices, + # floor/ceil, one voxel of margin for the taps. extents, per_voxel = (10, 12, 14), (1.0, 1.0, 2.0) # array order (z, y, x) expected = [] for axis, extent in enumerate(extents): reach = 4.0 / per_voxel[axis] - low, high = 4 - 0.5 - reach, 6 - 0.5 + reach + low, high = 4 - 0.5, 6 - 0.5 + reach expected.append((max(0, int(np.floor(low)) - 1), min(extent, int(np.ceil(high)) + 2))) assert [(part.start, part.stop) for part in window] == expected + # And it is strictly less than what the same field priced as a radius would have read. + assert all(part.stop - part.start < 2 * reach + 5 for part, reach in zip(window, (4.0, 4.0, 2.0), strict=True)) def test_an_oblique_case_grows_its_window_on_every_axis(tmp_path: Path) -> None: @@ -169,8 +178,10 @@ def test_a_field_with_no_bound_still_streams_with_windows_measured_at_run(tmp_pa # The plan prices as if the field were zero: the target's outer faces plus the taps' voxel. assert [(part.start, part.stop) for part in priced] == [(2, 8), (2, 8), (2, 8)] - # The run pays the shift the values actually hold: 4 um at spacing 2 is 2 voxels, on x alone. - assert [(part.start, part.stop) for part in measured] == [(2, 8), (2, 8), (0, 10)] + # The run pays the shift the values actually hold: 4 um at spacing 2 is 2 voxels, on x alone, + # and only where the field points -- the window keeps its identity faces (the field displaces + # nothing outside its grid) and extends one-sidedly along x. + assert [(part.start, part.stop) for part in measured] == [(2, 8), (2, 8), (2, 10)] def test_sizing_and_sampling_share_one_field_read(tmp_path: Path, monkeypatch: pytest.MonkeyPatch) -> None: