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: