From 6c14ee18b73e5b1427c2cd6b27088e87e63f294c Mon Sep 17 00:00:00 2001 From: Nathan Clack Date: Tue, 15 Sep 2026 00:54:46 +0000 Subject: [PATCH 01/10] Add NGFF spatial query resolution --- docs/api.md | 20 + docs/pipeline.md | 24 +- docs/spatial.md | 279 +++++++++++++ mkdocs.yml | 1 + python/CMakeLists.txt | 2 + python/damacy/__init__.py | 37 +- python/damacy/_api.c | 19 +- python/damacy/_api.h | 5 + python/damacy/_components.c | 6 + python/damacy/_native.c | 3 +- python/damacy/_native.pyi | 24 ++ python/damacy/_spatial.c | 373 +++++++++++++++++ python/damacy/_spatial.py | 255 ++++++++++++ python/tests/test_spatial.py | 659 +++++++++++++++++++++++++++++++ src/CMakeLists.txt | 6 +- src/damacy.h | 1 + src/damacy_spatial.h | 127 ++++++ src/damacy_status.c | 2 + src/ngff/load.c | 131 ++++++ src/ngff/metadata.c | 514 ++++++++++++++++++++++++ src/ngff/ngff.h | 23 ++ src/query/spatial.c | 273 +++++++++++++ src/store/metadata_store_async.c | 23 +- src/store/metadata_store_async.h | 6 + src/util/json.c | 48 +++ src/util/json.h | 12 + tests/CMakeLists.txt | 2 + tests/test_json.c | 29 ++ tests/test_spatial.c | 335 ++++++++++++++++ 29 files changed, 3214 insertions(+), 25 deletions(-) create mode 100644 docs/spatial.md create mode 100644 python/damacy/_spatial.c create mode 100644 python/damacy/_spatial.py create mode 100644 python/tests/test_spatial.py create mode 100644 src/damacy_spatial.h create mode 100644 src/ngff/load.c create mode 100644 src/ngff/metadata.c create mode 100644 src/ngff/ngff.h create mode 100644 src/query/spatial.c create mode 100644 tests/test_spatial.c diff --git a/docs/api.md b/docs/api.md index 6dd8a136..90b13d81 100644 --- a/docs/api.md +++ b/docs/api.md @@ -28,6 +28,24 @@ The full public surface of the `damacy` package. ::: damacy.CudaExecutor +## Spatial queries + +::: damacy.NgffImage + +::: damacy.NgffAxis + +::: damacy.NgffLevel + +::: damacy.SpatialResolver + +::: damacy.SpatialQuery + +::: damacy.ResolvedSpatialQuery + +::: damacy.Sampler + +::: damacy.NgffLimits + ## Output and limits ::: damacy.BatchSpec @@ -84,6 +102,8 @@ The full public surface of the `damacy` package. ::: damacy.BudgetExceeded +::: damacy.UnsupportedOperation + ::: damacy.ShutdownError ::: damacy.PoolStarved diff --git a/docs/pipeline.md b/docs/pipeline.md index 38c07a99..bce339e1 100644 --- a/docs/pipeline.md +++ b/docs/pipeline.md @@ -53,8 +53,10 @@ with damacy.Pipeline( del array ``` -The paths must name little-endian numeric Zarr v3 **arrays**. An NGFF group is -not resolved to a level automatically. The examples assume decoded chunks no larger than 2 MiB. +The paths in `Sample` and `IndexQuery` must name little-endian numeric Zarr v3 +**arrays**. Use an [NGFF image and spatial resolver](spatial.md) to select a +level from a group before pushing an aligned crop. The examples assume decoded +chunks no larger than 2 MiB. `FileMetadataReader` supplies small asynchronous metadata reads. `ZarrMetadata` supplies Zarr interpretation and cache capacities; `MetadataCache` is a capacity @@ -321,15 +323,11 @@ raise `InvalidArgument`. On Linux CUDA defaults on, builds both executors, and requires the CUDA toolkit, nvCOMP, and a runtime NVIDIA driver. GDS requires a CUDA build. -## Future queries +## Spatial queries -Rectangular and indexed queries copy or gather existing voxels. Spatial -queries will describe a fixed output tensor, a transform -from output coordinates to source space, and sampler settings including -interpolation, antialiasing, and boundary handling. - -For NGFF images, resolving a spatial query will choose an appropriate source -level from the requested sampling scale before enumerating chunks. The plan -will record the chosen array and transform. Resampling also needs interpolation -halos and dependencies across chunks. These are later planner and executor -operations; they are not implemented by the current box query. +[Spatial queries](spatial.md) resolve a fixed output grid, transform, and sampler +against loaded NGFF metadata. The resolver chooses a source level before chunk +planning. Aligned crops can use either executor now. Requests that still need +resampling expose their geometry and bounds, but fail when converted or pushed +for decoding. Future resampling operations will carry these dependencies across +source chunks in the prepared plan. diff --git a/docs/spatial.md b/docs/spatial.md new file mode 100644 index 00000000..03c53bb5 --- /dev/null +++ b/docs/spatial.md @@ -0,0 +1,279 @@ +# NGFF images and spatial queries + +A spatial query describes a fixed output grid mapped into a reference image. +Metadata loading and query resolution are separate from chunk planning and decoding: + +```mermaid +flowchart LR + reader[Metadata reader] --> image[Loaded NGFF image] + image --> resolver[Spatial resolver] + query[Transform and sampler] --> resolver + output[Fixed output shape] --> resolver + resolver --> resolved[Source array, transform, and bounds] + resolved --> planner[Chunk planner] + planner --> executor[CPU or CUDA executor] +``` + +The resolver chooses a source level before the chunk planner reads data. It +also reports whether that source grid can supply the output by copying voxels. +**Execution currently requires an aligned crop.** Rotation, reflection, shear, +fractional offsets, and scaling that remains after level selection require a +future resampling operation. Those queries can be resolved and inspected now; +converting or submitting them raises `UnsupportedOperation` (`DAMACY_UNSUPPORTED`). + +## Assemble the stages + +This example uses a three-dimensional image with axes ordered `z, y, x`. The +reference array must contain the requested crop. See +[pipeline composition](pipeline.md) for executor limits and batch ownership. + +```python +import damacy +import numpy as np + +metadata_reader = damacy.FileMetadataReader(concurrency=8) +image = damacy.NgffImage( + "/data/image.zarr", + reader=metadata_reader, + multiscale_index=0, + limits=damacy.NgffLimits(max_levels=16, max_metadata_bytes=4 << 20), +) +output = damacy.BatchSpec(samples=1, shape=(16, 64, 64), dtype="f32") +resolver = damacy.SpatialResolver(image=image, output=output) +query = damacy.SpatialQuery( + output_to_reference=( + (1, 0, 0, 8), + (0, 1, 0, 16), + (0, 0, 1, 16), + ), + sampler=damacy.Sampler(filter="nearest", boundary="error"), +) +resolved = resolver.resolve(query) + +planner = damacy.ChunkPlanner( + metadata=damacy.ZarrMetadata( + reader=metadata_reader, cache=damacy.MetadataCache() + ), + limits=damacy.PlanLimits(), +) +executor = damacy.CpuExecutor( + reader=damacy.FileReader(workers=4), + limits=damacy.CpuLimits(max_memory_bytes=256 << 20), +) +with damacy.Pipeline( + planner=planner, + executor=executor, + output=output, + queues=damacy.QueueLimits(lookahead_samples=8), +) as pipeline: + pipeline.push([resolved]) + with pipeline.pop() as batch: + array = np.from_dlpack(batch) + assert array.shape == (1, 16, 64, 64) +``` + +`NgffImage` loads the group's `zarr.json` and each level's array metadata once. +It owns an immutable description. `SpatialResolver.resolve()` performs no I/O +and can be called from multiple threads. A result owns its source URI and +geometry, so it remains usable after the image and resolver are released. +Source metadata and data must remain unchanged while the image is in use. + +Loading is synchronous and uses the supplied reader exclusively for the +duration of the call. A reader already serving a pipeline or another load is +rejected. The example reuses its reader after loading finishes. Use a separate +reader if loading more images while a pipeline is running. + +## Coordinates and the output grid + +All public spatial coordinates use **voxel corners**. Voxel index `i` covers +`[i, i + 1)` and has its center at `i + 0.5`. The volume covered by voxel zero +starts at coordinate zero, matching interval and indexed queries. Coordinates +in `output_to_reference` use the highest-resolution level, `image.levels[0]`, +as the reference; they are not physical coordinates. + +For rank `N`, the Python transform is an `N × (N + 1)` matrix. Its last column +is the offset. In C, `damacy_affine.linear` and `.offset` represent the same map: + +```text +reference_corner = linear * output_corner + offset +``` + +Rows and columns follow the metadata's array axis order, including time and +channel dimensions. The `BatchSpec` supplies the rank and fixed output shape. +Output centers are the grid points `(j[0] + 0.5, ..., j[N-1] + 0.5)`. +Only spatial axes may mix, rotate, reflect, or scale. Time and channel axes +require identity rows/columns and integral offsets, and must remain in bounds. +Combining spatial resampling with indexed time/channel selections is a later +extension; use `IndexQuery` for current arbitrary index selections. + +The spatial linear map must be finite and numerically nonsingular. Shapes and +transformed sample coordinates are limited to `2**52 - 1` index units so +half-voxel centers remain representable. Rank, dtype conversion, output size, +and every sampler parameter are validated before a resolution is returned. + +## NGFF interpretation + +The loader supports the [OME-Zarr 0.5](https://ngff.openmicroscopy.org/0.5/) +`attributes.ome` layout on Zarr v3. Select the `multiscales` entry explicitly +with `multiscale_index`. Supported images have two or three spatial axes and +optional time/channel axes, in the declared order. `dimension_names` must +match those axes at every level. Arrays must use supported numeric Zarr codecs +and a common dtype. Custom/unspecified axis types and transforms stored in +external arrays are currently unsupported. + +Dataset transforms must have one positive scale and an optional following +translation. Optional transforms shared by all levels follow the same rule. +Levels must have nondecreasing scales and nonincreasing shapes along each +axis; nonspatial transforms and shapes must be unchanged. Paths must be +relative child paths. Axis units are retained as declared; no physical-unit +conversion or inference is needed for reference-level queries. + +NGFF uses center-origin coordinates, as described in its +[coordinate convention](https://ngff.openmicroscopy.org/rfc/5/#coordinate-convention). +The adapter derives a corner-origin map once while loading metadata. If the +NGFF center-index map at level `l` is `p = S_l * i_l + T_l`, then: + +```text +R_l = S_l / S_0 +origin_l = (T_l - T_0) / S_0 + (1 - R_l) / 2 +reference_corner = R_l * level_corner + origin_l +``` + +These calculations are componentwise. A transform shared by every level +cancels when mapping to reference-level indices. Level zero is exactly the +identity. For a 2× level with a half-reference-voxel center translation, +`origin_l` is zero. A scale-only 2× level instead has `origin_l = -0.5`; +Damacy preserves that declared alignment. The API never asks callers to add +or remove NGFF center offsets themselves. + +`NgffLevel.scale_to_reference` and `.origin_reference_index` expose this +adapted map. Given query matrix `A` and offset `b`, resolution computes: + +```text +output_to_source.linear = inverse(R_l) * A +output_to_source.offset = inverse(R_l) * (b - origin_l) +``` + +## Choosing a level + +`SpatialQuery.level` defaults to `"auto"` (`DAMACY_LEVEL_AUTO` in C). The +resolver selects the coarsest level whose spacing is no coarser than the +requested output spacing in any spatial direction. It tests the smallest +singular value of `inverse(R_l) * A`, restricted to spatial axes, against one. +The comparison allows `64 * DBL_EPSILON` for numerical error. If no level +qualifies, it selects level zero for upsampling. Rotation, anisotropic scale, +and shear are included in this calculation. Source position, boundary mode, +and decoder capabilities do not change the level choice. + +An explicit nonnegative level index overrides that policy. A bad index fails; +there is no fallback from a requested level. For an aligned crop of level `l` +beginning at index `beg`, use the level's scale as the query's diagonal linear +map and `origin_l + R_l * beg` as its offset. + +The current copy path requires the resolved linear map to equal identity, +its offsets to be integral, and its source bounds to be inside the array. +Fractional values are not rounded into a crop. Small floating-point differences +in metadata can therefore make a result require resampling. `requires_resampling` +reports this, and `as_sample()` checks it before returning an interval query. +The pipeline still checks the resulting shape against its own `BatchSpec`. + +## Sampler and bounds + +Each spatial query carries a `Sampler`. Its current point filters are: + +| Filter | Source indices contributing at corner coordinate `s` | +| --- | --- | +| `nearest` | `floor(s)`; a tie on a voxel boundary selects the larger index | +| `linear` | `floor(s - 0.5)` and `ceil(s - 0.5)`, combined across spatial axes | + +Linear interpolation omits neighbors with zero weight. These definitions also +apply to negative coordinates. No additional antialias filter is currently +specified; choosing a pyramid level uses the downsampling already present in +the image, whose construction is outside the sampler's control. + +| Boundary | Interpretation | +| --- | --- | +| `error` | Reject if any contributing index is outside the source | +| `constant` | Extend the source with the finite `constant_value` | +| `clamp` | Replace each out-of-range spatial index with the nearest edge index | + +`constant_value` must be zero unless the boundary is `constant`. Filter and +boundary enums have no valid zero value in C. Unknown settings fail validation. +Boundary rules apply to the chosen source level, not the reference array. + +`source_bounds_index` is a half-open integer box enclosing all source indices +that contribute to the output, before boundary handling. `read_bounds_index` +accounts for boundary handling: constant extension clips to the array, while +clamp includes the touched edge voxels. Empty read intervals are valid for +constant-only outputs. The bounds include interpolation support. They are +conservative boxes around transformed crops, not a list of exact touched chunks. +Boundary sampling currently requires resampling and cannot be decoded yet. + +## Public C API + +Include `damacy_spatial.h` alongside the pipeline API. A two-dimensional +identity crop can be resolved and pushed as follows; `pipeline` must already +have matching output geometry and an executor: + +```c +struct damacy_metadata_reader* reader = NULL; +struct damacy_ngff_image* image = NULL; +struct damacy_spatial_resolution* resolved = NULL; +struct damacy_ngff_limits limits = { + .max_levels = 16, .max_metadata_bytes = 4 << 20 +}; +struct damacy_batch_spec output = { + .dtype = DAMACY_F32, .sample_shape = {64, 64}, + .sample_rank = 2, .samples_per_batch = 1 +}; +struct damacy_spatial_query query = { + .output_to_reference = { + .linear = {{1, 0}, {0, 1}}, .offset = {16, 24} + }, + .sampler = { + .filter = DAMACY_FILTER_NEAREST, .boundary = DAMACY_BOUNDARY_ERROR + }, + .level = DAMACY_LEVEL_AUTO +}; +struct damacy_sample sample; +enum damacy_status status = + damacy_file_metadata_reader_create(8, NULL, &reader); +if (status == DAMACY_OK) + status = damacy_ngff_image_load(reader, "/data/image.zarr", 0, &limits, &image); +if (status == DAMACY_OK) + status = damacy_spatial_resolve(image, &output, &query, &resolved); +if (status == DAMACY_OK) + status = damacy_spatial_resolution_sample(resolved, &sample); +if (status == DAMACY_OK) { + struct damacy_push_result pushed = damacy_push( + pipeline, (struct damacy_sample_slice){.beg = &sample, .end = &sample + 1}); + status = pushed.status; +} +damacy_spatial_resolution_destroy(resolved); +damacy_ngff_image_destroy(image); +damacy_metadata_reader_destroy(reader); +``` + +A production caller must retry the unconsumed suffix when `damacy_push` returns +`DAMACY_AGAIN`, keeping the sample's owner alive until acceptance or abandonment. +The example releases it on return. Accepted pushes copy the URI and selectors. +`damacy_spatial_resolution_sample()` borrows its URI from `resolved` and zeroes +its output on failure. Image and resolution info accessors return borrowed, +immutable views; their owners must remain alive while reading those views. +Resolution copies the information it needs from the image. Destruction accepts +null pointers. C affine entries outside the configured rank are unused. + +`max_metadata_bytes` bounds the sum of input JSON file sizes; files exceeding +the remaining budget are rejected before allocating or reading their contents. +`max_levels` separately bounds the number of level records. Exhaustion returns +`DAMACY_BUDGET`. Metadata loading performs no chunk reads and requires no GPU. + +## Next execution step + +A later planner operation can consume the owned resolution directly, group +source chunks needed by each output tile, and preserve the transform and +sampler in the prepared plan. CPU resampling needs a bounded decoded working +set and all contributing chunks available before writing a tile. CUDA can +implement the same operation after that contract is tested. Neither executor +should interpret NGFF or choose a pyramid level. Additional filters and +antialiasing parameters belong on the query's sampler when implemented. diff --git a/mkdocs.yml b/mkdocs.yml index 9461ea66..ede077c4 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -78,6 +78,7 @@ markdown_extensions: nav: - Home: index.md - Pipeline composition: pipeline.md + - Spatial queries: spatial.md - GPU memory budget: budget.md - Distributed: distributed.md - Async prefetch: prefetch.md diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index 5aec7972..75c2ec52 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -12,6 +12,7 @@ python_add_library( damacy/_log_sink.c damacy/_api.c damacy/_components.c + damacy/_spatial.c ) target_link_libraries(_native PRIVATE damacy log warnings) @@ -26,6 +27,7 @@ set_target_properties( PROPERTIES LIBRARY_OUTPUT_DIRECTORY "${CMAKE_CURRENT_BINARY_DIR}/damacy" ) configure_file(damacy/__init__.py damacy/__init__.py COPYONLY) +configure_file(damacy/_spatial.py damacy/_spatial.py COPYONLY) if(BUILD_TESTING) execute_process( diff --git a/python/damacy/__init__.py b/python/damacy/__init__.py index 83539869..8f448500 100644 --- a/python/damacy/__init__.py +++ b/python/damacy/__init__.py @@ -77,6 +77,10 @@ "MetadataCache", "Metric", "NativeCudaError", + "NgffAxis", + "NgffImage", + "NgffLevel", + "NgffLimits", "NotFound", "NumaStrategy", "OutOfMemory", @@ -85,12 +89,17 @@ "PoolStarved", "QueueLimits", "RankMismatch", + "ResolvedSpatialQuery", "Sample", + "Sampler", "ShutdownError", + "SpatialQuery", + "SpatialResolver", "Stats", "Status", "StorageError", "TryAgain", + "UnsupportedOperation", "ZarrMetadata", "max_concurrency", "set_log_level", @@ -241,6 +250,7 @@ class Status(IntEnum): OOM = _native.STATUS_OOM BUDGET = _native.STATUS_BUDGET SHUTDOWN = _native.STATUS_SHUTDOWN + UNSUPPORTED = _native.STATUS_UNSUPPORTED # ---- exceptions --------------------------------------------------------- @@ -305,6 +315,10 @@ class ShutdownError(DamacyError): """Pipeline destroyed or in failed state.""" +class UnsupportedOperation(DamacyError): + """The requested metadata or execution operation is not supported.""" + + class PoolStarved(DamacyError): """Raised when :meth:`Pipeline.pop` waits longer than ``Config.pop_timeout_s`` for the next batch. @@ -332,6 +346,7 @@ class PoolStarved(DamacyError): _native.STATUS_OOM: OutOfMemory, _native.STATUS_BUDGET: BudgetExceeded, _native.STATUS_SHUTDOWN: ShutdownError, + _native.STATUS_UNSUPPORTED: UnsupportedOperation, } @@ -1609,8 +1624,10 @@ def __init__( # _pending_buf is the head iterator's already-pulled-but-not-yet- # pushed samples; held flat to avoid wrapping `it` in successive # itertools.chain() layers under sustained backpressure. - self._pending: deque[Iterator[Sample | IndexQuery]] = deque() - self._pending_buf: list[Sample | IndexQuery] = [] + self._pending: deque[Iterator[Sample | IndexQuery | ResolvedSpatialQuery]] = ( + deque() + ) + self._pending_buf: list[Sample | IndexQuery | ResolvedSpatialQuery] = [] # damacy_pop has no timed variant; on timeout the worker stays # parked inside it and the next pop() adopts the same thread. self._pop_lock = threading.Lock() @@ -1691,7 +1708,9 @@ def __exit__( # ---- pipeline ---------------------------------------------------- - def push(self, samples: Iterable[Sample | IndexQuery]) -> None: + def push( + self, samples: Iterable[Sample | IndexQuery | ResolvedSpatialQuery] + ) -> None: """Queue samples for processing. Accepts any iterable (list, generator, infinite generator, …); large or unbounded sources are pulled lazily as :meth:`pop` frees space. @@ -1915,3 +1934,15 @@ def stats_reset(self) -> None: # Avoid "no handler for damacy" warnings in apps that don't configure # logging; users opt in by attaching their own handler / level. logging.getLogger(__name__).addHandler(logging.NullHandler()) + + +from ._spatial import ( # noqa: E402 + NgffAxis, + NgffImage, + NgffLevel, + NgffLimits, + ResolvedSpatialQuery, + Sampler, + SpatialQuery, + SpatialResolver, +) diff --git a/python/damacy/_api.c b/python/damacy/_api.c index a9793c25..b2eff758 100644 --- a/python/damacy/_api.c +++ b/python/damacy/_api.c @@ -1322,12 +1322,19 @@ api_register_types(PyObject* m) const char* name; int value; } statuses[] = { - { "STATUS_OK", DAMACY_OK }, { "STATUS_AGAIN", DAMACY_AGAIN }, - { "STATUS_INVAL", DAMACY_INVAL }, { "STATUS_NOTFOUND", DAMACY_NOTFOUND }, - { "STATUS_DTYPE", DAMACY_DTYPE }, { "STATUS_RANK", DAMACY_RANK }, - { "STATUS_IO", DAMACY_IO }, { "STATUS_DECODE", DAMACY_DECODE }, - { "STATUS_CUDA", DAMACY_CUDA }, { "STATUS_OOM", DAMACY_OOM }, - { "STATUS_BUDGET", DAMACY_BUDGET }, { "STATUS_SHUTDOWN", DAMACY_SHUTDOWN }, + { "STATUS_OK", DAMACY_OK }, + { "STATUS_AGAIN", DAMACY_AGAIN }, + { "STATUS_INVAL", DAMACY_INVAL }, + { "STATUS_NOTFOUND", DAMACY_NOTFOUND }, + { "STATUS_DTYPE", DAMACY_DTYPE }, + { "STATUS_RANK", DAMACY_RANK }, + { "STATUS_IO", DAMACY_IO }, + { "STATUS_DECODE", DAMACY_DECODE }, + { "STATUS_CUDA", DAMACY_CUDA }, + { "STATUS_OOM", DAMACY_OOM }, + { "STATUS_BUDGET", DAMACY_BUDGET }, + { "STATUS_SHUTDOWN", DAMACY_SHUTDOWN }, + { "STATUS_UNSUPPORTED", DAMACY_UNSUPPORTED }, }; for (size_t i = 0; i < sizeof statuses / sizeof statuses[0]; ++i) { if (PyModule_AddIntConstant(m, statuses[i].name, statuses[i].value) < 0) diff --git a/python/damacy/_api.h b/python/damacy/_api.h index 1478e108..2526342d 100644 --- a/python/damacy/_api.h +++ b/python/damacy/_api.h @@ -15,3 +15,8 @@ api_pipeline_from_components(struct damacy_planner* planner, const struct damacy_batch_spec* output, const struct damacy_queue_limits* queues, PyObject* dependencies); + +int +spatial_register(PyObject* module); +struct damacy_metadata_reader* +api_metadata_reader(PyObject* capsule); diff --git a/python/damacy/_components.c b/python/damacy/_components.c index 2753a2c2..e4486f0e 100644 --- a/python/damacy/_components.c +++ b/python/damacy/_components.c @@ -321,3 +321,9 @@ components_register(PyObject* module) { return PyModule_AddFunctions(module, methods); } + +struct damacy_metadata_reader* +api_metadata_reader(PyObject* capsule) +{ + return component_value(capsule, METADATA_READER); +} diff --git a/python/damacy/_native.c b/python/damacy/_native.c index 7f62e524..35e4e726 100644 --- a/python/damacy/_native.c +++ b/python/damacy/_native.c @@ -258,7 +258,8 @@ module_exec(PyObject* m) PyErr_SetString(PyExc_RuntimeError, "failed to install damacy log sink"); return -1; } - if (api_register_types(m) != 0 || components_register(m) != 0) + if (api_register_types(m) != 0 || components_register(m) != 0 || + spatial_register(m) != 0) return -1; #ifdef DAMACY_HAS_CUDA if (PyModule_AddIntConstant(m, "CUDA_ENABLED", 1) < 0) diff --git a/python/damacy/_native.pyi b/python/damacy/_native.pyi index 03f432b7..3897ec18 100644 --- a/python/damacy/_native.pyi +++ b/python/damacy/_native.pyi @@ -33,6 +33,7 @@ STATUS_CUDA: Final[int] STATUS_OOM: Final[int] STATUS_BUDGET: Final[int] STATUS_SHUTDOWN: Final[int] +STATUS_UNSUPPORTED: Final[int] # ---- damacy_dtype integers ---------------------------------------------- @@ -252,3 +253,26 @@ def compose_pipeline( prepared_batches: int, /, ) -> Pipeline: ... +def ngff_load( + reader: object, + uri: str, + multiscale_index: int, + max_levels: int, + max_metadata_bytes: int, + /, +) -> object: ... +def ngff_info(image: object, /) -> dict[str, Any]: ... +def spatial_resolve( + image: object, + shape: tuple[int, ...], + samples: int, + dtype: int, + transform: tuple[tuple[float, ...], ...], + filter: int, + boundary: int, + constant_value: float, + level: int, + /, +) -> object: ... +def spatial_info(resolution: object, /) -> dict[str, Any]: ... +def spatial_sample(resolution: object, /) -> dict[str, Any]: ... diff --git a/python/damacy/_spatial.c b/python/damacy/_spatial.c new file mode 100644 index 00000000..0e2f410c --- /dev/null +++ b/python/damacy/_spatial.c @@ -0,0 +1,373 @@ +#define PY_SSIZE_T_CLEAN +#include + +#include "_api.h" +#include "damacy_spatial.h" + +static const char image_name[] = "damacy.NgffImage"; +static const char resolution_name[] = "damacy.SpatialResolution"; + +static void +image_destroy(PyObject* capsule) +{ + damacy_ngff_image_destroy(PyCapsule_GetPointer(capsule, image_name)); +} + +static void +resolution_destroy(PyObject* capsule) +{ + damacy_spatial_resolution_destroy( + PyCapsule_GetPointer(capsule, resolution_name)); +} + +static int +put(PyObject* dict, const char* key, PyObject* value) +{ + if (!value) + return -1; + int result = PyDict_SetItemString(dict, key, value); + Py_DECREF(value); + return result; +} + +static PyObject* +integers(const int64_t* values, uint8_t rank) +{ + PyObject* result = PyTuple_New(rank); + if (!result) + return NULL; + for (uint8_t i = 0; i < rank; ++i) { + PyObject* value = PyLong_FromLongLong(values[i]); + if (!value) { + Py_DECREF(result); + return NULL; + } + PyTuple_SET_ITEM(result, i, value); + } + return result; +} + +static PyObject* +doubles(const double* values, uint8_t rank) +{ + PyObject* result = PyTuple_New(rank); + if (!result) + return NULL; + for (uint8_t i = 0; i < rank; ++i) { + PyObject* value = PyFloat_FromDouble(values[i]); + if (!value) { + Py_DECREF(result); + return NULL; + } + PyTuple_SET_ITEM(result, i, value); + } + return result; +} + +static PyObject* +affine(const struct damacy_affine* transform, uint8_t rank) +{ + PyObject* result = PyTuple_New(rank); + if (!result) + return NULL; + for (uint8_t i = 0; i < rank; ++i) { + double row[DAMACY_MAX_RANK + 1]; + for (uint8_t j = 0; j < rank; ++j) + row[j] = transform->linear[i][j]; + row[rank] = transform->offset[i]; + PyObject* value = doubles(row, rank + 1); + if (!value) { + Py_DECREF(result); + return NULL; + } + PyTuple_SET_ITEM(result, i, value); + } + return result; +} + +static PyObject* +bounds(const struct damacy_aabb* box) +{ + PyObject* result = PyTuple_New(box->rank); + if (!result) + return NULL; + for (uint8_t i = 0; i < box->rank; ++i) { + PyObject* value = Py_BuildValue( + "(LL)", (long long)box->dims[i].beg, (long long)box->dims[i].end); + if (!value) { + Py_DECREF(result); + return NULL; + } + PyTuple_SET_ITEM(result, i, value); + } + return result; +} + +static PyObject* +load_image(PyObject* self, PyObject* args) +{ + (void)self; + PyObject* reader_object; + const char* uri; + unsigned int index, levels; + unsigned long long bytes; + if (!PyArg_ParseTuple( + args, "OsIIK", &reader_object, &uri, &index, &levels, &bytes)) + return NULL; + struct damacy_metadata_reader* reader = api_metadata_reader(reader_object); + if (!reader) + return NULL; + struct damacy_ngff_limits limits = { .max_levels = levels, + .max_metadata_bytes = bytes }; + struct damacy_ngff_image* image = NULL; + enum damacy_status status; + Py_BEGIN_ALLOW_THREADS status = + damacy_ngff_image_load(reader, uri, index, &limits, &image); + Py_END_ALLOW_THREADS if (status != DAMACY_OK) return api_raise_status( + status, "load NGFF metadata"); + PyObject* capsule = PyCapsule_New(image, image_name, image_destroy); + if (!capsule) + damacy_ngff_image_destroy(image); + return capsule; +} + +static PyObject* +image_info(PyObject* self, PyObject* capsule) +{ + (void)self; + struct damacy_ngff_image* image = PyCapsule_GetPointer(capsule, image_name); + if (!image) + return NULL; + const struct damacy_ngff_info* info = damacy_ngff_image_info(image); + PyObject* result = PyDict_New(); + PyObject* axes = PyTuple_New(info->rank); + PyObject* levels = PyTuple_New(info->level_count); + if (!result || !axes || !levels) + goto Fail; + for (uint8_t i = 0; i < info->rank; ++i) { + const struct damacy_ngff_axis* axis = &info->axes[i]; + const char* kind = axis->kind == DAMACY_NGFF_SPACE ? "space" + : axis->kind == DAMACY_NGFF_TIME ? "time" + : "channel"; + PyObject* value = Py_BuildValue("(ssz)", axis->name, kind, axis->unit); + if (!value) + goto Fail; + PyTuple_SET_ITEM(axes, i, value); + } + for (uint32_t i = 0; i < info->level_count; ++i) { + const struct damacy_ngff_level* level = &info->levels[i]; + PyObject* value = PyDict_New(); + if (!value) + goto Fail; + PyTuple_SET_ITEM(levels, i, value); + if (put(value, "uri", PyUnicode_FromString(level->uri)) || + put(value, "shape", integers(level->shape, info->rank)) || + put(value, + "scale_to_reference", + doubles(level->scale_to_reference, info->rank)) || + put(value, + "origin_reference_index", + doubles(level->origin_reference_index, info->rank))) + goto Fail; + } + if (PyDict_SetItemString(result, "axes", axes) || + PyDict_SetItemString(result, "levels", levels) || + put(result, "data_type", PyUnicode_FromString(info->data_type))) + goto Fail; + Py_DECREF(axes); + Py_DECREF(levels); + return result; +Fail: + Py_XDECREF(result); + Py_XDECREF(axes); + Py_XDECREF(levels); + return NULL; +} + +static int +parse_transform(PyObject* object, uint8_t rank, struct damacy_affine* transform) +{ + PyObject* rows = PySequence_Fast(object, "transform must be a sequence"); + if (!rows) + return -1; + if (PySequence_Fast_GET_SIZE(rows) != rank) { + PyErr_SetString(PyExc_ValueError, "transform rank must match output rank"); + Py_DECREF(rows); + return -1; + } + for (uint8_t i = 0; i < rank; ++i) { + PyObject* row = PySequence_Fast(PySequence_Fast_GET_ITEM(rows, i), + "transform rows must be sequences"); + if (!row) { + Py_DECREF(rows); + return -1; + } + if (PySequence_Fast_GET_SIZE(row) != rank + 1) { + PyErr_SetString(PyExc_ValueError, "transform rows need rank + 1 entries"); + Py_DECREF(row); + Py_DECREF(rows); + return -1; + } + for (uint8_t j = 0; j <= rank; ++j) { + double value = PyFloat_AsDouble(PySequence_Fast_GET_ITEM(row, j)); + if (PyErr_Occurred()) { + Py_DECREF(row); + Py_DECREF(rows); + return -1; + } + if (j == rank) + transform->offset[i] = value; + else + transform->linear[i][j] = value; + } + Py_DECREF(row); + } + Py_DECREF(rows); + return 0; +} + +static PyObject* +resolve_query(PyObject* self, PyObject* args) +{ + (void)self; + PyObject *capsule, *shape_object, *transform; + unsigned int samples, dtype; + int filter, boundary, level; + double value; + if (!PyArg_ParseTuple(args, + "OOIIOiidi", + &capsule, + &shape_object, + &samples, + &dtype, + &transform, + &filter, + &boundary, + &value, + &level)) + return NULL; + struct damacy_ngff_image* image = PyCapsule_GetPointer(capsule, image_name); + if (!image) + return NULL; + struct damacy_batch_spec output = { .samples_per_batch = samples, + .dtype = (enum damacy_dtype)dtype }; + PyObject* shape = PySequence_Fast(shape_object, "shape must be a sequence"); + if (!shape) + return NULL; + Py_ssize_t rank = PySequence_Fast_GET_SIZE(shape); + if (rank < 1 || rank > DAMACY_MAX_RANK) { + Py_DECREF(shape); + PyErr_SetString(PyExc_ValueError, "invalid output rank"); + return NULL; + } + output.sample_rank = (uint8_t)rank; + for (Py_ssize_t i = 0; i < rank; ++i) { + output.sample_shape[i] = + PyLong_AsLongLong(PySequence_Fast_GET_ITEM(shape, i)); + if (PyErr_Occurred()) { + Py_DECREF(shape); + return NULL; + } + } + Py_DECREF(shape); + struct damacy_spatial_query query = { + .sampler = { .filter = (enum damacy_filter)filter, + .boundary = (enum damacy_boundary)boundary, + .constant_value = value }, + .level = level + }; + if (parse_transform( + transform, output.sample_rank, &query.output_to_reference)) + return NULL; + struct damacy_spatial_resolution* resolution = NULL; + enum damacy_status status; + Py_BEGIN_ALLOW_THREADS status = + damacy_spatial_resolve(image, &output, &query, &resolution); + Py_END_ALLOW_THREADS if (status != DAMACY_OK) return api_raise_status( + status, "resolve spatial query"); + PyObject* result = + PyCapsule_New(resolution, resolution_name, resolution_destroy); + if (!result) + damacy_spatial_resolution_destroy(resolution); + return result; +} + +static PyObject* +resolution_info(PyObject* self, PyObject* capsule) +{ + (void)self; + struct damacy_spatial_resolution* resolution = + PyCapsule_GetPointer(capsule, resolution_name); + if (!resolution) + return NULL; + const struct damacy_spatial_info* info = + damacy_spatial_resolution_info(resolution); + PyObject* result = PyDict_New(); + if (!result) + return NULL; + if (put(result, "uri", PyUnicode_FromString(info->uri)) || + put(result, "level", PyLong_FromUnsignedLong(info->level)) || + put(result, "shape", integers(info->output_shape, info->rank)) || + put(result, "source_shape", integers(info->source_shape, info->rank)) || + put(result, + "output_to_source", + affine(&info->output_to_source, info->rank)) || + put(result, "source_bounds_index", bounds(&info->source_bounds_index)) || + put(result, "read_bounds_index", bounds(&info->read_bounds_index)) || + put(result, + "requires_resampling", + PyBool_FromLong(info->requires_resampling))) { + Py_DECREF(result); + return NULL; + } + return result; +} + +static PyObject* +resolution_sample(PyObject* self, PyObject* capsule) +{ + (void)self; + struct damacy_spatial_resolution* resolution = + PyCapsule_GetPointer(capsule, resolution_name); + if (!resolution) + return NULL; + struct damacy_sample sample; + enum damacy_status status = + damacy_spatial_resolution_sample(resolution, &sample); + if (status != DAMACY_OK) + return api_raise_status(status, + "submit spatial query: resampling required"); + PyObject* axes = PyList_New(sample.rank); + if (!axes) + return NULL; + for (uint8_t i = 0; i < sample.rank; ++i) { + PyObject* axis = Py_BuildValue("(s(LL))", + "interval", + (long long)sample.axes[i].interval.beg, + (long long)sample.axes[i].interval.end); + if (!axis) { + Py_DECREF(axes); + return NULL; + } + PyList_SET_ITEM(axes, i, axis); + } + PyObject* result = + Py_BuildValue("{s:s,s:O}", "uri", sample.uri, "axes", axes); + Py_DECREF(axes); + return result; +} + +static PyMethodDef methods[] = { + { "ngff_load", load_image, METH_VARARGS, NULL }, + { "ngff_info", image_info, METH_O, NULL }, + { "spatial_resolve", resolve_query, METH_VARARGS, NULL }, + { "spatial_info", resolution_info, METH_O, NULL }, + { "spatial_sample", resolution_sample, METH_O, NULL }, + { NULL, NULL, 0, NULL } +}; + +int +spatial_register(PyObject* module) +{ + return PyModule_AddFunctions(module, methods); +} diff --git a/python/damacy/_spatial.py b/python/damacy/_spatial.py new file mode 100644 index 00000000..3b44ded9 --- /dev/null +++ b/python/damacy/_spatial.py @@ -0,0 +1,255 @@ +from __future__ import annotations + +import math +import operator +import os +from collections.abc import Iterable +from dataclasses import dataclass +from typing import Any, Literal + +from . import BatchSpec, FileMetadataReader, Sample, _component, _native, _positive_int + + +def _index(value: int, name: str) -> int: + if isinstance(value, bool): + raise TypeError(f"{name} must be an integer") + result = operator.index(value) + if not 0 <= result <= (1 << 31) - 1: + raise ValueError(f"{name} must be between 0 and 2**31 - 1") + return result + + +def _finite_number(value: float) -> float: + if isinstance(value, (str, bytes, bool)): + raise TypeError("transform and sampler values must be numbers") + result = float(value) + if not math.isfinite(result): + raise ValueError("transform and sampler values must be finite") + return result + + +@dataclass(frozen=True, slots=True) +class NgffLimits: + """Maximum levels and total JSON bytes read when loading an image.""" + + max_levels: int = 32 + max_metadata_bytes: int = 4 << 20 + + def __post_init__(self) -> None: + _positive_int(self.max_levels, "max_levels", (1 << 31) - 1) + _positive_int(self.max_metadata_bytes, "max_metadata_bytes", (1 << 64) - 1) + + +@dataclass(frozen=True, slots=True) +class NgffAxis: + """An axis in the array's dimension order, with its declared NGFF unit.""" + + name: str + kind: Literal["space", "time", "channel"] + unit: str | None + + +@dataclass(frozen=True, slots=True) +class NgffLevel: + """A source array and its voxel-corner mapping into reference-level indices.""" + + uri: str + shape: tuple[int, ...] + scale_to_reference: tuple[float, ...] + origin_reference_index: tuple[float, ...] + + +@dataclass(frozen=True, slots=True, init=False) +class NgffImage: + """Load an immutable OME-Zarr 0.5 image description through the given reader. + + Loading reads the image group's metadata and each level's array metadata. + It finishes before returning and releases the reader for other uses. + ``multiscale_index`` explicitly selects an entry in ``ome.multiscales``. + """ + + axes: tuple[NgffAxis, ...] + levels: tuple[NgffLevel, ...] + data_type: str + _native: object + + def __init__( + self, + uri: str | os.PathLike[str], + *, + reader: FileMetadataReader, + multiscale_index: int, + limits: NgffLimits | None = None, + ) -> None: + if not isinstance(reader, FileMetadataReader): + raise TypeError("reader must be a FileMetadataReader") + if limits is None: + limits = NgffLimits() + if not isinstance(limits, NgffLimits): + raise TypeError("limits must be NgffLimits") + path = os.fspath(uri) + if not isinstance(path, str) or not path or "\0" in path: + raise ValueError("uri must be a nonempty path without NUL characters") + native = _component( + _native.ngff_load, + reader._native, + path, + _index(multiscale_index, "multiscale_index"), + limits.max_levels, + limits.max_metadata_bytes, + ) + info = _native.ngff_info(native) + object.__setattr__(self, "_native", native) + object.__setattr__( + self, "axes", tuple(NgffAxis(*axis) for axis in info["axes"]) + ) + object.__setattr__( + self, "levels", tuple(NgffLevel(**level) for level in info["levels"]) + ) + object.__setattr__(self, "data_type", info["data_type"]) + + +@dataclass(frozen=True, slots=True) +class Sampler: + """Point interpolation and treatment of source indices outside spatial axes. + + Nearest selects ``floor(source_corner)``. Linear interpolates between + centers at ``index + 0.5``. Constant extends the source with + ``constant_value``; clamp repeats the closest edge value; error rejects + a query needing out-of-bounds samples. No extra antialias filter is applied. + """ + + filter: Literal["nearest", "linear"] = "linear" + boundary: Literal["error", "constant", "clamp"] = "error" + constant_value: float = 0 + + def __post_init__(self) -> None: + if self.filter not in ("nearest", "linear"): + raise ValueError("filter must be 'nearest' or 'linear'") + if self.boundary not in ("error", "constant", "clamp"): + raise ValueError("boundary must be 'error', 'constant', or 'clamp'") + value = _finite_number(self.constant_value) + if self.boundary != "constant" and value != 0: + raise ValueError("constant_value requires boundary='constant'") + object.__setattr__(self, "constant_value", value) + + +@dataclass(frozen=True, slots=True, init=False) +class SpatialQuery: + """Map a fixed output grid into reference-level voxel-corner coordinates. + + ``output_to_reference`` has rank rows and rank + 1 columns. The last + column is the translation; the preceding columns are the linear map. + Rows and columns follow NGFF array axis order, including time/channel + dimensions, which only allow identity plus integer translation. + ``level='auto'`` selects the coarsest level no coarser than the output + spacing in any direction, falling back to level zero for upsampling. + """ + + output_to_reference: tuple[tuple[float, ...], ...] + sampler: Sampler + level: int | Literal["auto"] + + def __init__( + self, + *, + output_to_reference: Iterable[Iterable[float]], + sampler: Sampler, + level: int | Literal["auto"] = "auto", + ) -> None: + if not isinstance(sampler, Sampler): + raise TypeError("sampler must be a Sampler") + matrix = tuple( + tuple(_finite_number(v) for v in row) for row in output_to_reference + ) + rank = len(matrix) + if not 2 <= rank <= _native.MAX_RANK or any( + len(row) != rank + 1 for row in matrix + ): + raise ValueError( + "output_to_reference must have rank rows and rank + 1 columns" + ) + if isinstance(level, str): + if level != "auto": + raise ValueError("level must be 'auto' or a nonnegative integer") + else: + level = _index(level, "level") + object.__setattr__(self, "output_to_reference", matrix) + object.__setattr__(self, "sampler", sampler) + object.__setattr__(self, "level", level) + + +@dataclass(frozen=True, slots=True, init=False) +class ResolvedSpatialQuery: + """Owned source selection and geometry, independent of the loaded image. + + Aligned results can be pushed to a Pipeline or converted with ``as_sample``. + Other results raise UnsupportedOperation when converted or pushed; their + geometry remains available for inspection without invoking a decoder. + """ + + uri: str + level: int + shape: tuple[int, ...] + source_shape: tuple[int, ...] + output_to_source: tuple[tuple[float, ...], ...] + sampler: Sampler + source_bounds_index: tuple[tuple[int, int], ...] + read_bounds_index: tuple[tuple[int, int], ...] + requires_resampling: bool + _native: object + + @classmethod + def _from_native(cls, native: object, sampler: Sampler) -> ResolvedSpatialQuery: + result = object.__new__(cls) + object.__setattr__(result, "_native", native) + object.__setattr__(result, "sampler", sampler) + for name, value in _native.spatial_info(native).items(): + object.__setattr__(result, name, value) + return result + + def _to_native(self) -> dict[str, Any]: + return _component(_native.spatial_sample, self._native) + + def as_sample(self) -> Sample: + """Return an interval sample, or fail if resampling is required.""" + sample = self._to_native() + return Sample(uri=sample["uri"], aabb=tuple(axis[1] for axis in sample["axes"])) + + +@dataclass(frozen=True, slots=True) +class SpatialResolver: + """Resolve queries using injected image metadata and fixed output geometry. + + Resolution does no I/O and can be shared across threads. The same + BatchSpec should be supplied to the downstream Pipeline. + """ + + image: NgffImage + output: BatchSpec + + def __post_init__(self) -> None: + if not isinstance(self.image, NgffImage) or not isinstance( + self.output, BatchSpec + ): + raise TypeError("SpatialResolver requires an NgffImage and a BatchSpec") + if len(self.image.axes) != len(self.output.shape): + raise ValueError("output rank must match the NGFF image") + + def resolve(self, query: SpatialQuery) -> ResolvedSpatialQuery: + """Choose the source level, map coordinates, and bound source reads.""" + if not isinstance(query, SpatialQuery): + raise TypeError("query must be a SpatialQuery") + native = _component( + _native.spatial_resolve, + self.image._native, + self.output.shape, + self.output.samples, + int(self.output.dtype), + query.output_to_reference, + {"nearest": 1, "linear": 2}[query.sampler.filter], + {"error": 1, "constant": 2, "clamp": 3}[query.sampler.boundary], + query.sampler.constant_value, + -1 if query.level == "auto" else query.level, + ) + return ResolvedSpatialQuery._from_native(native, query.sampler) diff --git a/python/tests/test_spatial.py b/python/tests/test_spatial.py new file mode 100644 index 00000000..bc3194c4 --- /dev/null +++ b/python/tests/test_spatial.py @@ -0,0 +1,659 @@ +from __future__ import annotations + +import ctypes +import dataclasses +import gc +import json +from concurrent.futures import ThreadPoolExecutor +from typing import Literal + +import damacy +import numpy as np +import pytest + + +def write_image(path, *, shape=(32, 32), axes=None, scales=None, data=False): + path.mkdir() + rank = len(shape) + axes = axes or [(name, "space") for name in ("z", "y", "x")[-rank:]] + scales = scales or [ + [2**level if kind == "space" else 1 for _, kind in axes] for level in range(3) + ] + datasets = [] + arrays = [] + for level, scale in enumerate(scales): + level_shape = tuple(int(n / s) for n, s in zip(shape, scale, strict=True)) + values = ( + np.arange(np.prod(level_shape), dtype=np.uint16).reshape(level_shape) + + level * 10000 + ) + arrays.append(values) + array_path = path / str(level) + array_path.mkdir() + chunks = [ + min(8, n) if kind == "space" else 1 + for n, (_, kind) in zip(level_shape, axes, strict=True) + ] + metadata = { + "zarr_format": 3, + "node_type": "array", + "shape": level_shape, + "dimension_names": [name for name, _ in axes], + "data_type": "uint16", + "fill_value": 0, + "chunk_grid": {"name": "regular", "configuration": {"chunk_shape": chunks}}, + "chunk_key_encoding": { + "name": "default", + "configuration": {"separator": "/"}, + }, + "codecs": [{"name": "bytes", "configuration": {"endian": "little"}}], + } + (array_path / "zarr.json").write_text(json.dumps(metadata)) + datasets.append( + { + "path": str(level), + "coordinateTransformations": [ + {"type": "scale", "scale": scale}, + { + "type": "translation", + "translation": [(s - 1) / 2 for s in scale], + }, + ], + } + ) + if data: + grid = tuple( + (n + c - 1) // c for n, c in zip(level_shape, chunks, strict=True) + ) + for coordinate in np.ndindex(grid): + slices = tuple( + slice(i * c, min((i + 1) * c, n)) + for i, c, n in zip(coordinate, chunks, level_shape, strict=True) + ) + block = np.zeros(chunks, dtype=np.uint16) + selected = values[slices] + block[tuple(slice(0, n) for n in selected.shape)] = selected + file = array_path.joinpath("c", *map(str, coordinate)) + file.parent.mkdir(parents=True, exist_ok=True) + file.write_bytes(block.tobytes()) + metadata = { + "zarr_format": 3, + "node_type": "group", + "attributes": { + "ome": { + "version": "0.5", + "multiscales": [ + { + "axes": [{"name": name, "type": kind} for name, kind in axes], + "datasets": datasets, + } + ], + } + }, + } + (path / "zarr.json").write_text(json.dumps(metadata)) + return arrays + + +def edit(path, keys, value): + metadata = json.loads(path.read_text()) + target = metadata + for key in keys[:-1]: + target = target[key] + if value is None: + del target[keys[-1]] + else: + target[keys[-1]] = value + path.write_text(json.dumps(metadata)) + + +def load(path, **kwargs): + return damacy.NgffImage( + path, + reader=damacy.FileMetadataReader(concurrency=2), + multiscale_index=0, + **kwargs, + ) + + +def query( + scale=1, offset=(0, 0), *, sampler=None, level: int | Literal["auto"] = "auto" +): + rank = len(offset) + return damacy.SpatialQuery( + output_to_reference=tuple( + (*(scale if i == j else 0 for j in range(rank)), origin) + for i, origin in enumerate(offset) + ), + sampler=sampler or damacy.Sampler(filter="nearest"), + level=level, + ) + + +@pytest.fixture(params=["cpu", "cuda"]) +def spatial_executor(request): + reader = damacy.FileReader(workers=2, max_inflight_reads=2) + if request.param == "cpu": + return damacy.CpuExecutor( + reader=reader, + limits=damacy.CpuLimits( + max_memory_bytes=8 << 20, + decode_workers=2, + max_encoded_chunk_bytes=4096, + max_decoded_chunk_bytes=4096, + ), + ) + request.getfixturevalue("cuda_ctx") + return damacy.CudaExecutor( + reader=reader, + device=0, + limits=damacy.CudaLimits( + max_gpu_memory_bytes=256 << 20, + max_chunk_bytes=4096, + max_chunks_per_wave=2, + max_substreams_per_chunk=8, + ), + ) + + +def pipeline(executor, output, reader=None): + return damacy.Pipeline( + planner=damacy.ChunkPlanner( + metadata=damacy.ZarrMetadata( + reader=reader or damacy.FileMetadataReader(concurrency=2), + cache=damacy.MetadataCache(array_entries=16, shard_index_entries=32), + ), + limits=damacy.PlanLimits( + max_chunks=128, max_chunk_bytes=4096, max_shards_per_sample=4 + ), + ), + executor=executor, + output=output, + queues=damacy.QueueLimits(lookahead_samples=4), + pop_timeout_s=5, + ) + + +def read_batch(batch): + result = np.empty(batch.info.shape, dtype=np.float32) + if batch.info.device_type == damacy.DeviceType.CPU: + ctypes.memmove(result.ctypes.data, batch.info.data, result.nbytes) + else: + driver = ctypes.CDLL("libcuda.so.1") + copy = driver.cuMemcpyDtoH_v2 + copy.argtypes = [ctypes.c_void_p, ctypes.c_uint64, ctypes.c_size_t] + copy.restype = ctypes.c_int + assert copy(result.ctypes.data, batch.info.data, result.nbytes) == 0 + return result + + +@pytest.mark.parametrize("level", [0, 1, 2]) +def test_resolved_crops_decode_at_selected_level(tmp_path, spatial_executor, level): + root = tmp_path / "image" + arrays = write_image(root, data=True) + reader = damacy.FileMetadataReader(concurrency=2) + image = damacy.NgffImage(root, reader=reader, multiscale_index=0) + output = damacy.BatchSpec(2, (4, 4)) + resolver = damacy.SpatialResolver(image=image, output=output) + factor = 2**level + resolved = resolver.resolve(query(factor, (factor * 2, factor * 3))) + assert resolved.level == level + assert not resolved.requires_resampling + assert resolved.as_sample().aabb == ((2, 6), (3, 7)) + assert resolved.output_to_source == ((1, 0, 2), (0, 1, 3)) + del resolver, image + gc.collect() + with pipeline(spatial_executor, output, reader) as p: + p.push([resolved, resolved.as_sample()]) + with p.pop() as batch: + expected = arrays[level][2:6, 3:7] + np.testing.assert_array_equal( + read_batch(batch), np.stack([expected, expected]) + ) + + +def test_resampling_fails_before_decoding_and_pipeline_recovers( + tmp_path, spatial_executor +): + root = tmp_path / "image" + arrays = write_image(root, data=True) + output = damacy.BatchSpec(1, (4, 4)) + resolver = damacy.SpatialResolver(load(root), output) + rotated = resolver.resolve( + damacy.SpatialQuery( + output_to_reference=((1, -1, 12), (1, 1, 4)), + sampler=damacy.Sampler(filter="linear"), + ) + ) + assert rotated.requires_resampling + with pytest.raises(damacy.UnsupportedOperation, match="resampling required"): + rotated.as_sample() + with pipeline(spatial_executor, output) as p: + with pytest.raises(damacy.UnsupportedOperation, match="resampling required"): + p.push([rotated]) + p.push([resolver.resolve(query())]) + with p.pop() as batch: + np.testing.assert_array_equal(read_batch(batch)[0], arrays[0][:4, :4]) + + +def test_fixed_shape_is_checked_at_push(tmp_path, spatial_executor): + root = tmp_path / "image" + write_image(root) + resolved = damacy.SpatialResolver(load(root), damacy.BatchSpec(1, (3, 3))).resolve( + query() + ) + with ( + pipeline(spatial_executor, damacy.BatchSpec(1, (4, 4))) as p, + pytest.raises(damacy.InvalidArgument), + ): + p.push([resolved]) + + +def test_ngff_center_offsets_stay_in_adapter(tmp_path): + root = tmp_path / "image" + write_image(root) + path = root / "zarr.json" + prefix = ("attributes", "ome", "multiscales", 0) + transforms = ( + [ + {"type": "scale", "scale": [0.5, 1]}, + {"type": "translation", "translation": [10, 20]}, + ], + [ + {"type": "scale", "scale": [1, 2]}, + {"type": "translation", "translation": [10.25, 20.5]}, + ], + [ + {"type": "scale", "scale": [2, 4]}, + {"type": "translation", "translation": [10, 20]}, + ], + ) + for i, value in enumerate(transforms): + edit(path, (*prefix, "datasets", i, "coordinateTransformations"), value) + edit( + path, + (*prefix, "coordinateTransformations"), + [ + {"type": "scale", "scale": [3, 7]}, + {"type": "translation", "translation": [900, -100]}, + ], + ) + image = load(root) + assert image.levels[0].origin_reference_index == (0, 0) + assert image.levels[1].origin_reference_index == (0, 0) + assert image.levels[2].origin_reference_index == (-1.5, -1.5) + resolver = damacy.SpatialResolver(image, damacy.BatchSpec(1, (4, 4))) + resolved = resolver.resolve(query(2)) + assert resolved.level == 1 and not resolved.requires_resampling + assert resolver.resolve(query(4)).requires_resampling + + +def test_anisotropic_rotated_level_selection_and_time_channel_axes(tmp_path): + root = tmp_path / "image" + write_image( + root, + shape=(2, 3, 16, 16, 16), + axes=[ + ("t", "time"), + ("c", "channel"), + ("z", "space"), + ("y", "space"), + ("x", "space"), + ], + scales=[[1, 1, 1, 1, 1], [1, 1, 1, 2, 4]], + ) + image = load(root) + resolver = damacy.SpatialResolver(image, damacy.BatchSpec(1, (1, 1, 2, 2, 2))) + transform = np.zeros((5, 6)) + transform[:5, :5] = np.eye(5) + transform[:2, 5] = (1, 2) + transform[2:5, 2:5] = np.array([[0, -1, 0], [2, 0, 0], [0, 0, 4]]) + transform[2:, 5] = (4, 4, 0) + resolved = resolver.resolve( + damacy.SpatialQuery( + output_to_reference=transform, sampler=damacy.Sampler(filter="nearest") + ) + ) + assert resolved.level == 1 + assert resolved.output_to_source[2] == (0, 0, 0, -1, 0, 4) + assert resolved.output_to_source[3] == (0, 0, 1, 0, 0, 2) + assert resolved.source_bounds_index[:2] == ((1, 2), (2, 3)) + transform[2, 3] = -0.5 + assert ( + resolver.resolve( + damacy.SpatialQuery(output_to_reference=transform, sampler=damacy.Sampler()) + ).level + == 0 + ) + transform[0, 0] = 2 + with pytest.raises(damacy.InvalidArgument): + resolver.resolve( + damacy.SpatialQuery(output_to_reference=transform, sampler=damacy.Sampler()) + ) + + +def test_selection_uses_smallest_spacing_not_column_lengths(tmp_path): + root = tmp_path / "image" + write_image(root) + resolver = damacy.SpatialResolver(load(root), damacy.BatchSpec(1, (2, 2))) + q = damacy.SpatialQuery( + output_to_reference=((3, 2, 0), (2, 3, 0)), + sampler=damacy.Sampler(filter="nearest"), + ) + assert resolver.resolve(q).level == 0 + forced = dataclasses.replace(q, level=1) + assert resolver.resolve(forced).level == 1 + + +@pytest.mark.parametrize( + "filter,boundary,offset,source,reads", + [ + ("nearest", "error", -0.25, (0, 4), (0, 4)), + ("linear", "constant", -0.25, (-1, 4), (0, 4)), + ("linear", "clamp", -0.25, (-1, 4), (0, 4)), + ("nearest", "constant", -10, (-10, -6), (0, 0)), + ("nearest", "clamp", -10, (-10, -6), (0, 1)), + ("nearest", "constant", 40, (40, 44), (32, 32)), + ("nearest", "clamp", 40, (40, 44), (31, 32)), + ], +) +def test_source_and_read_bounds(tmp_path, filter, boundary, offset, source, reads): + root = tmp_path / "image" + write_image(root) + resolver = damacy.SpatialResolver(load(root), damacy.BatchSpec(1, (4, 4))) + resolved = resolver.resolve( + query( + offset=(offset, 0), sampler=damacy.Sampler(filter=filter, boundary=boundary) + ) + ) + assert resolved.source_bounds_index[0] == source + assert resolved.read_bounds_index[0] == reads + assert resolved.requires_resampling + + +@pytest.mark.parametrize( + "transform", + [ + ((0, 0, 0), (0, 1, 0)), + ((1, 2, 0), (2, 4, 0)), + ((1, 0, -1), (0, 1, 0)), + ((1e100, 0, 0), (0, 1, 0)), + ((1, 0, 2**52), (0, 1, 0)), + ], +) +def test_invalid_geometry(tmp_path, transform): + root = tmp_path / "image" + write_image(root) + resolver = damacy.SpatialResolver(load(root), damacy.BatchSpec(1, (4, 4))) + with pytest.raises(damacy.InvalidArgument): + resolver.resolve( + damacy.SpatialQuery(output_to_reference=transform, sampler=damacy.Sampler()) + ) + + +@pytest.mark.parametrize( + "keys,value", + [ + (("version",), "0.4"), + (("multiscales", 0, "axes", 0, "name"), "x"), + (("multiscales", 0, "axes", 0, "type"), "custom"), + (("multiscales", 0, "datasets", 1, "path"), "0"), + (("multiscales", 0, "datasets", 1, "path"), "../outside"), + (("multiscales", 0, "datasets", 1, "path"), "/outside"), + (("multiscales", 0, "datasets", 1, "path"), "1//array"), + ( + ("multiscales", 0, "datasets", 1, "coordinateTransformations", 0, "scale"), + [2], + ), + ( + ("multiscales", 0, "datasets", 1, "coordinateTransformations", 0, "scale"), + [0, 2], + ), + ( + ("multiscales", 0, "datasets", 1, "coordinateTransformations", 0, "scale"), + [0.5, 2], + ), + ( + ("multiscales", 0, "datasets", 1, "coordinateTransformations", 0, "scale"), + [-2, 2], + ), + ( + ("multiscales", 0, "datasets", 1, "coordinateTransformations", 0, "type"), + "affine", + ), + (("multiscales", 0, "datasets", 1, "coordinateTransformations"), []), + ], +) +def test_rejects_unsupported_or_ambiguous_metadata(tmp_path, keys, value): + root = tmp_path / "image" + write_image(root) + edit(root / "zarr.json", ("attributes", "ome", *keys), value) + with pytest.raises((damacy.InvalidArgument, damacy.UnsupportedOperation)): + load(root) + + +@pytest.mark.parametrize( + "key,value", + [ + ("dimension_names", ["x", "y"]), + ("dimension_names", None), + ("dimension_names", ["y"]), + ("shape", [16, 16, 16]), + ("shape", [64, 16]), + ("data_type", "float32"), + ], +) +def test_level_arrays_must_agree_with_ngff(tmp_path, key, value): + root = tmp_path / "image" + write_image(root) + edit(root / "1" / "zarr.json", (key,), value) + with pytest.raises(damacy.DamacyError): + load(root) + + +def test_limits_errors_and_reader_reuse(tmp_path): + root = tmp_path / "image" + write_image(root) + reader = damacy.FileMetadataReader(concurrency=2) + metadata_bytes = sum(p.stat().st_size for p in root.rglob("zarr.json")) + for limits in [ + damacy.NgffLimits(max_levels=2), + damacy.NgffLimits(max_metadata_bytes=metadata_bytes - 1), + damacy.NgffLimits(max_metadata_bytes=1), + ]: + with pytest.raises(damacy.BudgetExceeded): + damacy.NgffImage(root, reader=reader, multiscale_index=0, limits=limits) + image = damacy.NgffImage( + root, + reader=reader, + multiscale_index=0, + limits=damacy.NgffLimits(max_metadata_bytes=metadata_bytes), + ) + assert image.levels[0].shape == (32, 32) + with pytest.raises(damacy.NotFound): + load(tmp_path / "missing") + with pytest.raises(damacy.InvalidArgument): + damacy.NgffImage(root, reader=reader, multiscale_index=1) + + +def test_multiscale_selection_and_escaped_strings(tmp_path): + root = tmp_path / "image" + write_image(root, axes=[("μ", "space"), ("x", "space")]) + path = root / "zarr.json" + metadata = json.loads(path.read_text()) + entry = metadata["attributes"]["ome"]["multiscales"][0] + metadata["attributes"]["ome"]["multiscales"].insert(0, {}) + path.write_text( + json.dumps(metadata) + .replace('"version"', '"ver\\u0073ion"') + .replace('"0.5"', '"0.\\u0035"') + ) + image = damacy.NgffImage( + root, reader=damacy.FileMetadataReader(concurrency=2), multiscale_index=1 + ) + assert image.axes[0].name == "μ" + assert entry["axes"][0]["name"] == image.axes[0].name + with pytest.raises(damacy.InvalidArgument): + load(root) + + +@pytest.mark.parametrize( + "mutate", + [ + lambda text: text.replace( + '"version": "0.5"', '"version": "0.5", "version": "0.4"' + ), + lambda text: text.replace( + '"scale": [1, 1]', '"scale": [1, 1], "scale": [2, 2]' + ), + lambda text: text + " garbage", + lambda text: text[:-1] + ",}", + ], +) +def test_malformed_json_is_rejected(tmp_path, mutate): + root = tmp_path / "image" + write_image(root) + path = root / "zarr.json" + path.write_text(mutate(path.read_text())) + with pytest.raises(damacy.DamacyError): + load(root) + + +def test_resolution_is_independent_of_files_and_can_run_in_threads(tmp_path): + root = tmp_path / "image" + write_image(root) + resolver = damacy.SpatialResolver(load(root), damacy.BatchSpec(1, (4, 4))) + for metadata in root.rglob("zarr.json"): + metadata.unlink() + with ThreadPoolExecutor(max_workers=4) as workers: + results = list(workers.map(resolver.resolve, [query(2)] * 16)) + assert all(r.level == 1 and r.as_sample().aabb == ((0, 4), (0, 4)) for r in results) + + +def test_query_values_are_copied_and_validated(): + matrix = [[1, 0, 2], [0, 1, 3]] + q = damacy.SpatialQuery(output_to_reference=matrix, sampler=damacy.Sampler()) + matrix[0][2] = 100 + assert q.output_to_reference[0][2] == 2 + with pytest.raises(dataclasses.FrozenInstanceError): + q.level = 2 # type: ignore[misc] + for invalid in [(), ((1, 0), (0, 1)), ((1, 0, float("nan")), (0, 1, 0))]: + with pytest.raises(ValueError): + damacy.SpatialQuery(output_to_reference=invalid, sampler=damacy.Sampler()) + with pytest.raises(ValueError): + damacy.Sampler(boundary="clamp", constant_value=3) + with pytest.raises(ValueError): + damacy.NgffLimits(max_levels=0) + + +def test_anisotropic_crop_decodes_with_time_and_channel_axes( + tmp_path, spatial_executor +): + root = tmp_path / "image" + arrays = write_image( + root, + shape=(2, 3, 16, 16, 16), + axes=[ + ("t", "time"), + ("c", "channel"), + ("z", "space"), + ("y", "space"), + ("x", "space"), + ], + scales=[[1, 1, 1, 1, 1], [1, 1, 1, 2, 4]], + data=True, + ) + output = damacy.BatchSpec(1, (1, 2, 4, 4, 4)) + resolver = damacy.SpatialResolver(load(root), output) + transform = np.column_stack((np.diag([1, 1, 1, 2, 4]), [1, 1, 4, 4, 0])) + resolved = resolver.resolve( + damacy.SpatialQuery(output_to_reference=transform, sampler=damacy.Sampler()) + ) + assert resolved.level == 1 and not resolved.requires_resampling + with pipeline(spatial_executor, output) as p: + p.push([resolved]) + with p.pop() as batch: + np.testing.assert_array_equal( + read_batch(batch)[0], arrays[1][1:2, 1:3, 4:8, 2:6, :4] + ) + + +def test_automatic_level_selection_matches_singular_values(tmp_path): + root = tmp_path / "image" + scales = np.array([[1, 1, 1], [1, 2, 4], [2, 4, 8]]) + write_image(root, shape=(32, 32, 32), scales=scales.tolist()) + resolver = damacy.SpatialResolver(load(root), damacy.BatchSpec(1, (2, 2, 2))) + rng = np.random.default_rng(532) + for _ in range(64): + left, _ = np.linalg.qr(rng.normal(size=(3, 3))) + right, _ = np.linalg.qr(rng.normal(size=(3, 3))) + matrix = left @ np.diag(2 ** rng.uniform(-2, 4, 3)) @ right + expected = 0 + for level, scale in enumerate(scales): + if np.linalg.svd(matrix / scale[:, None], compute_uv=False)[-1] >= 1: + expected = level + resolved = resolver.resolve( + damacy.SpatialQuery( + output_to_reference=np.column_stack((matrix, [8, 8, 8])), + sampler=damacy.Sampler(boundary="constant"), + ) + ) + assert resolved.level == expected + + +def test_active_metadata_reader_is_not_reused(tmp_path, spatial_executor): + root = tmp_path / "image" + write_image(root) + reader = damacy.FileMetadataReader(concurrency=2) + with ( + pipeline(spatial_executor, damacy.BatchSpec(1, (4, 4)), reader), + pytest.raises(damacy.InvalidArgument), + ): + damacy.NgffImage(root, reader=reader, multiscale_index=0) + assert damacy.NgffImage(root, reader=reader, multiscale_index=0).levels + + +@pytest.mark.parametrize("path", ["zarr.json", "1/zarr.json"]) +def test_empty_metadata_fails_cleanly(tmp_path, path): + root = tmp_path / "image" + write_image(root) + (root / path).write_text("") + with pytest.raises(damacy.InvalidArgument): + load(root) + + +def test_native_query_validation(tmp_path): + from damacy import _native + + root = tmp_path / "image" + write_image(root) + image = load(root) + args = (image._native, (4, 4), 1, 0) + with pytest.raises(ValueError): + _native.spatial_resolve(*args, ((1, 0), (0, 1)), 1, 1, 0, -1) + for filter, boundary, value, level in [ + (0, 1, 0, -1), + (1, 0, 0, -1), + (1, 1, 1, -1), + (1, 2, float("nan"), -1), + (1, 1, 0, -2), + (1, 1, 0, 3), + ]: + with pytest.raises(_native.DamacyError): + _native.spatial_resolve( + *args, ((1, 0, 0), (0, 1, 0)), filter, boundary, value, level + ) + + +@pytest.mark.parametrize( + "key,value", + [("shape", "[16,16]"), ("data_type", '"float32"'), ("node_type", '"group"')], +) +def test_duplicate_array_fields_are_rejected(tmp_path, key, value): + root = tmp_path / "image" + write_image(root) + path = root / "0" / "zarr.json" + text = path.read_text() + path.write_text(text[:-1] + f', "{key}": {value}' + "}") + with pytest.raises(damacy.InvalidArgument): + load(root) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index d2872158..941569d6 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -278,6 +278,10 @@ if(NOT DAMACY_FUZZ) add_src_lib(pipeline_components SOURCES pipeline/components.h pipeline/components.c damacy_status.c LINKS damacy_config store platform log ) + add_src_lib(spatial_query + SOURCES damacy_spatial.h ngff/ngff.h ngff/metadata.c ngff/load.c query/spatial.c + LINKS pipeline_components metadata_store_async zarr json m + ) add_src_lib(metadata_planner SOURCES pipeline/zarr_planner.c LINKS pipeline_components plan_builder prefetcher array_meta shard_index metadata_store_async lookahead @@ -291,7 +295,7 @@ if(NOT DAMACY_FUZZ) damacy_lifecycle.c damacy_push.c damacy_plan.c damacy_pop.c damacy_scheduler.c LINKS pipeline_components metadata_planner cpu_executor scheduler damacy_stats - query_selection + query_selection spatial_query ) if(DAMACY_CUDA) add_cuda_lib(cuda_executor diff --git a/src/damacy.h b/src/damacy.h index 60596987..826c2fb5 100644 --- a/src/damacy.h +++ b/src/damacy.h @@ -55,6 +55,7 @@ extern "C" DAMACY_OOM, // host allocation failed (calloc/malloc returned null) DAMACY_BUDGET, // a configured cap is too small to satisfy the request DAMACY_SHUTDOWN, // pipeline destroyed or in failed state + DAMACY_UNSUPPORTED, }; // Human-readable name for a status code; safe for log/error messages. diff --git a/src/damacy_spatial.h b/src/damacy_spatial.h new file mode 100644 index 00000000..252a6fde --- /dev/null +++ b/src/damacy_spatial.h @@ -0,0 +1,127 @@ +#pragma once + +#include "damacy_pipeline.h" + +#ifdef __cplusplus +extern "C" +{ +#endif + + struct damacy_ngff_image; + struct damacy_spatial_resolution; + + enum damacy_ngff_axis_kind + { + DAMACY_NGFF_SPACE = 1, + DAMACY_NGFF_TIME, + DAMACY_NGFF_CHANNEL, + }; + + struct damacy_ngff_axis + { + const char* name; + const char* unit; + enum damacy_ngff_axis_kind kind; + }; + + struct damacy_ngff_level + { + const char* uri; + int64_t shape[DAMACY_MAX_RANK]; + double scale_to_reference[DAMACY_MAX_RANK]; + double origin_reference_index[DAMACY_MAX_RANK]; + }; + + struct damacy_ngff_info + { + uint8_t rank; + const char* data_type; + struct damacy_ngff_axis axes[DAMACY_MAX_RANK]; + uint32_t level_count; + const struct damacy_ngff_level* levels; + }; + + struct damacy_ngff_limits + { + uint32_t max_levels; + uint64_t max_metadata_bytes; + }; + + struct damacy_affine + { + double linear[DAMACY_MAX_RANK][DAMACY_MAX_RANK]; + double offset[DAMACY_MAX_RANK]; + }; + + enum damacy_filter + { + DAMACY_FILTER_NEAREST = 1, + DAMACY_FILTER_LINEAR, + }; + + enum damacy_boundary + { + DAMACY_BOUNDARY_ERROR = 1, + DAMACY_BOUNDARY_CONSTANT, + DAMACY_BOUNDARY_CLAMP, + }; + + struct damacy_sampler + { + enum damacy_filter filter; + enum damacy_boundary boundary; + double constant_value; + }; + + enum + { + DAMACY_LEVEL_AUTO = -1 + }; + + struct damacy_spatial_query + { + struct damacy_affine output_to_reference; + struct damacy_sampler sampler; + int32_t level; + }; + + struct damacy_spatial_info + { + const char* uri; + uint32_t level; + uint8_t rank; + int64_t output_shape[DAMACY_MAX_RANK]; + int64_t source_shape[DAMACY_MAX_RANK]; + struct damacy_affine output_to_source; + struct damacy_sampler sampler; + struct damacy_aabb source_bounds_index; + struct damacy_aabb read_bounds_index; + int requires_resampling; + }; + + enum damacy_status damacy_ngff_image_load( + struct damacy_metadata_reader* reader, + const char* uri, + uint32_t multiscale_index, + const struct damacy_ngff_limits* limits, + struct damacy_ngff_image** out); + const struct damacy_ngff_info* damacy_ngff_image_info( + const struct damacy_ngff_image* image); + void damacy_ngff_image_destroy(struct damacy_ngff_image* image); + + enum damacy_status damacy_spatial_resolve( + const struct damacy_ngff_image* image, + const struct damacy_batch_spec* output, + const struct damacy_spatial_query* query, + struct damacy_spatial_resolution** out); + const struct damacy_spatial_info* damacy_spatial_resolution_info( + const struct damacy_spatial_resolution* resolution); + enum damacy_status damacy_spatial_resolution_sample( + const struct damacy_spatial_resolution* resolution, + struct damacy_sample* out); + void damacy_spatial_resolution_destroy( + struct damacy_spatial_resolution* resolution); + +#ifdef __cplusplus +} +#endif diff --git a/src/damacy_status.c b/src/damacy_status.c index f100a23b..0a04ba4f 100644 --- a/src/damacy_status.c +++ b/src/damacy_status.c @@ -28,6 +28,8 @@ damacy_status_str(enum damacy_status s) return "configured budget too small"; case DAMACY_SHUTDOWN: return "shutdown"; + case DAMACY_UNSUPPORTED: + return "unsupported operation"; } return "unknown"; } diff --git a/src/ngff/load.c b/src/ngff/load.c new file mode 100644 index 00000000..adac9699 --- /dev/null +++ b/src/ngff/load.c @@ -0,0 +1,131 @@ +#include "ngff/ngff.h" + +#include "pipeline/components.h" +#include "platform/platform.h" +#include "store/metadata_store_async.h" + +#include +#include + +struct metadata_read +{ + struct platform_mutex* mutex; + struct platform_cond* cond; + int done; + enum damacy_status status; + void* data; + size_t bytes; +}; + +static void +read_complete(void* user, enum damacy_status status, void* data, size_t bytes) +{ + struct metadata_read* read = user; + platform_mutex_lock(read->mutex); + read->status = status; + read->data = data; + read->bytes = bytes; + read->done = 1; + platform_cond_broadcast(read->cond); + platform_mutex_unlock(read->mutex); +} + +static enum damacy_status +read_metadata(struct metadata_store_async* store, + struct metadata_read* read, + const char* uri, + size_t* remaining_bytes) +{ + size_t len = strlen(uri); + if (len > SIZE_MAX - sizeof("/zarr.json")) + return DAMACY_BUDGET; + char* path = malloc(len + sizeof("/zarr.json")); + if (!path) + return DAMACY_OOM; + memcpy(path, uri, len); + if (len && path[len - 1] == '/') + --len; + memcpy(path + len, "/zarr.json", sizeof("/zarr.json")); + read->done = 0; + if (metadata_store_async_read_file_bounded( + store, path, *remaining_bytes, read_complete, read)) { + free(path); + return DAMACY_IO; + } + free(path); + platform_mutex_lock(read->mutex); + while (!read->done) + platform_cond_wait(read->cond, read->mutex); + platform_mutex_unlock(read->mutex); + if (read->status == DAMACY_OK) { + if (!read->data || !read->bytes) + return DAMACY_INVAL; + *remaining_bytes -= read->bytes; + } + return read->status; +} + +enum damacy_status +damacy_ngff_image_load(struct damacy_metadata_reader* reader, + const char* uri, + uint32_t multiscale_index, + const struct damacy_ngff_limits* limits, + struct damacy_ngff_image** out) +{ + if (!out) + return DAMACY_INVAL; + *out = NULL; + if (!reader || !uri || !*uri || !limits || !limits->max_levels || + limits->max_levels > INT32_MAX || !limits->max_metadata_bytes || + limits->max_metadata_bytes > SIZE_MAX) + return DAMACY_INVAL; + int expected = 0; + if (!atomic_compare_exchange_strong(&reader->active, &expected, 1)) + return DAMACY_INVAL; + struct metadata_store_async* store = metadata_store_async_create( + (int)reader->concurrency, NULL, &reader->latency); + struct metadata_read read = { .mutex = platform_mutex_new(), + .cond = platform_cond_new() }; + struct damacy_ngff_image* image = NULL; + enum damacy_status status = DAMACY_OOM; + if (!store || !read.mutex || !read.cond) + goto Done; + size_t remaining_bytes = (size_t)limits->max_metadata_bytes; + status = read_metadata(store, &read, uri, &remaining_bytes); + if (status != DAMACY_OK) + goto Done; + status = ngff_parse_group( + (struct cslice){ read.data, (const char*)read.data + read.bytes }, + uri, + multiscale_index, + limits->max_levels, + &image); + free(read.data); + read.data = NULL; + if (status != DAMACY_OK) + goto Done; + for (uint32_t i = 0; i < image->info.level_count; ++i) { + status = + read_metadata(store, &read, image->levels[i].uri, &remaining_bytes); + if (status != DAMACY_OK) + goto Done; + status = ngff_parse_array( + (struct cslice){ read.data, (const char*)read.data + read.bytes }, + image, + i); + free(read.data); + read.data = NULL; + if (status != DAMACY_OK) + goto Done; + } + *out = image; + image = NULL; +Done: + metadata_store_async_destroy(store); + platform_mutex_free(read.mutex); + platform_cond_free(read.cond); + atomic_store(&reader->active, 0); + damacy_ngff_image_destroy(image); + free(read.data); + return status; +} diff --git a/src/ngff/metadata.c b/src/ngff/metadata.c new file mode 100644 index 00000000..04381936 --- /dev/null +++ b/src/ngff/metadata.c @@ -0,0 +1,514 @@ +#include "ngff/ngff.h" + +#include "zarr/zarr_metadata.h" + +#include +#include +#include + +static enum json_err +elements(struct json_node node, struct json_iter* it) +{ + static const struct json_query query = { .kind = QUERY_ITER }; + if (node.type != JSON_ARRAY) + return JSON_ERR_TYPE; + return json_iter_init(node.s, &query, 1, it, NULL); +} + +static int +hex4(const char* p, uint32_t* out) +{ + *out = 0; + for (int i = 0; i < 4; ++i) { + unsigned char c = (unsigned char)p[i]; + uint32_t digit; + if (c >= '0' && c <= '9') + digit = c - '0'; + else if (c >= 'a' && c <= 'f') + digit = c - 'a' + 10; + else if (c >= 'A' && c <= 'F') + digit = c - 'A' + 10; + else + return 1; + *out = (*out << 4) | digit; + } + return 0; +} + +static enum damacy_status +read_string(struct json_node node, const char** out) +{ + if (node.type != JSON_STRING || !cslice_len(node.s)) + return DAMACY_INVAL; + char* text = malloc(cslice_len(node.s) + 1); + if (!text) + return DAMACY_OOM; + char* dst = text; + const char* p = node.s.beg; + while (p < node.s.end) { + unsigned char c = (unsigned char)*p++; + if (c != '\\') { + *dst++ = (char)c; + continue; + } + if (p == node.s.end) + goto Invalid; + c = (unsigned char)*p++; + if (c == 'u') { + uint32_t code; + if (node.s.end - p < 4 || hex4(p, &code)) + goto Invalid; + p += 4; + if (code >= 0xd800 && code <= 0xdbff) { + uint32_t low; + if (node.s.end - p < 6 || p[0] != '\\' || p[1] != 'u' || + hex4(p + 2, &low) || low < 0xdc00 || low > 0xdfff) + goto Invalid; + p += 6; + code = 0x10000 + ((code - 0xd800) << 10) + low - 0xdc00; + } else if (code >= 0xdc00 && code <= 0xdfff) { + goto Invalid; + } + if (!code) + goto Invalid; + if (code < 0x80) { + *dst++ = (char)code; + } else if (code < 0x800) { + *dst++ = (char)(0xc0 | (code >> 6)); + *dst++ = (char)(0x80 | (code & 0x3f)); + } else if (code < 0x10000) { + *dst++ = (char)(0xe0 | (code >> 12)); + *dst++ = (char)(0x80 | ((code >> 6) & 0x3f)); + *dst++ = (char)(0x80 | (code & 0x3f)); + } else { + *dst++ = (char)(0xf0 | (code >> 18)); + *dst++ = (char)(0x80 | ((code >> 12) & 0x3f)); + *dst++ = (char)(0x80 | ((code >> 6) & 0x3f)); + *dst++ = (char)(0x80 | (code & 0x3f)); + } + } else { + const char* escapes = "\"\\/bfnrt"; + const char* values = "\"\\/\b\f\n\r\t"; + const char* found = strchr(escapes, c); + if (!found) + goto Invalid; + *dst++ = values[found - escapes]; + } + } + *dst = 0; + *out = text; + return DAMACY_OK; +Invalid: + free(text); + return DAMACY_INVAL; +} + +static enum json_err +member(struct json_node node, const char* key, struct json_node* out) +{ + struct json_object_iter it; + enum json_err error = json_object_iter_init(node, &it); + if (error != JSON_OK) + return error; + struct json_node name, value; + int found = 0; + while ((error = json_object_iter_next(&it, &name, &value)) == JSON_OK) { + if (!cslice_len(name.s)) + continue; + const char* text = NULL; + enum damacy_status status = read_string(name, &text); + if (status != DAMACY_OK) + return status == DAMACY_OOM ? JSON_ERR_OOM : JSON_ERR_PARSE; + int matches = !strcmp(key, text); + free((void*)text); + if (matches) { + if (found) + return JSON_ERR_PARSE; + *out = value; + found = 1; + } + } + return error == JSON_ERR_NOT_FOUND && found ? JSON_OK : error; +} + +static int +string_equal(struct json_node node, const char* expected) +{ + const char* text = NULL; + if (read_string(node, &text) != DAMACY_OK) + return 0; + int result = !strcmp(text, expected); + free((void*)text); + return result; +} + +static int +root_document(struct cslice src, struct json_node* root) +{ + if (json_resolve(src, NULL, 0, root, NULL) || root->type != JSON_OBJECT) + return 1; + for (const char* p = root->s.end; p < src.end; ++p) + if (*p != ' ' && *p != '\t' && *p != '\n' && *p != '\r') + return 1; + return 0; +} + +static enum damacy_status +read_axes(struct json_node node, struct damacy_ngff_info* info) +{ + struct json_iter it; + if (elements(node, &it)) + return DAMACY_INVAL; + uint8_t spaces = 0; + int time_seen = 0, channel_seen = 0; + struct json_node axis; + enum json_err error; + while ((error = json_iter_next(&it, &axis)) == JSON_OK) { + if (info->rank == 5) + return DAMACY_RANK; + struct damacy_ngff_axis* out = &info->axes[info->rank++]; + struct json_node value; + if (member(axis, "name", &value)) + return DAMACY_INVAL; + enum damacy_status status = read_string(value, &out->name); + if (status != DAMACY_OK) + return status; + for (uint8_t i = 0; i + 1 < info->rank; ++i) + if (!strcmp(out->name, info->axes[i].name)) + return DAMACY_INVAL; + if (member(axis, "type", &value)) + return DAMACY_UNSUPPORTED; + if (string_equal(value, "space")) { + out->kind = DAMACY_NGFF_SPACE; + spaces++; + } else if (string_equal(value, "time")) { + if (info->rank != 1 || time_seen++) + return DAMACY_INVAL; + out->kind = DAMACY_NGFF_TIME; + } else if (string_equal(value, "channel")) { + if (spaces || channel_seen++) + return DAMACY_INVAL; + out->kind = DAMACY_NGFF_CHANNEL; + } else { + return DAMACY_UNSUPPORTED; + } + error = member(axis, "unit", &value); + if (error == JSON_OK) { + status = read_string(value, &out->unit); + if (status != DAMACY_OK) + return status; + } else if (error != JSON_ERR_NOT_FOUND) { + return DAMACY_INVAL; + } + } + if (error != JSON_ERR_NOT_FOUND || spaces < 2 || spaces > 3) + return DAMACY_INVAL; + return DAMACY_OK; +} + +static enum damacy_status +read_vector(struct json_node node, uint8_t rank, double* values) +{ + struct json_iter it; + if (elements(node, &it)) + return DAMACY_INVAL; + struct json_node value; + for (uint8_t d = 0; d < rank; ++d) + if (json_iter_next(&it, &value) || json_as_double(value, &values[d]) || + !isfinite(values[d])) + return DAMACY_INVAL; + return json_iter_next(&it, &value) == JSON_ERR_NOT_FOUND ? DAMACY_OK + : DAMACY_INVAL; +} + +static enum damacy_status +read_transform(struct json_node node, + uint8_t rank, + double* scale, + double* translation) +{ + struct json_iter it; + if (elements(node, &it)) + return DAMACY_INVAL; + int count = 0; + enum json_err error; + struct json_node transform; + while ((error = json_iter_next(&it, &transform)) == JSON_OK) { + struct json_node type, value; + if (member(transform, "type", &type)) + return DAMACY_INVAL; + const char* key = count == 0 ? "scale" : "translation"; + if (count > 1 || !string_equal(type, key)) + return DAMACY_UNSUPPORTED; + if (member(transform, "path", &value) != JSON_ERR_NOT_FOUND) + return DAMACY_UNSUPPORTED; + if (member(transform, key, &value)) + return DAMACY_INVAL; + enum damacy_status status = + read_vector(value, rank, count ? translation : scale); + if (status != DAMACY_OK) + return status; + ++count; + } + if (error != JSON_ERR_NOT_FOUND || !count) + return DAMACY_INVAL; + for (uint8_t d = 0; d < rank; ++d) + if (scale[d] <= 0) + return DAMACY_INVAL; + return DAMACY_OK; +} + +static int +relative_path(const char* path) +{ + if (!*path || *path == '/' || strchr(path, '\\') || strchr(path, ':')) + return 0; + const char* p = path; + do { + const char* end = strchr(p, '/'); + size_t len = end ? (size_t)(end - p) : strlen(p); + if (!len || (len == 1 && p[0] == '.') || + (len == 2 && p[0] == '.' && p[1] == '.')) + return 0; + p = end ? end + 1 : NULL; + } while (p); + return 1; +} + +static enum damacy_status +read_level(struct json_node node, + const char* uri, + uint8_t rank, + struct damacy_ngff_level* level) +{ + struct json_node value; + if (member(node, "path", &value)) + return DAMACY_INVAL; + const char* path = NULL; + enum damacy_status status = read_string(value, &path); + if (status != DAMACY_OK) + return status; + if (!relative_path(path)) { + free((void*)path); + return DAMACY_INVAL; + } + size_t root_len = strlen(uri), path_len = strlen(path); + if (root_len > SIZE_MAX - path_len - 2) { + free((void*)path); + return DAMACY_BUDGET; + } + char* joined = malloc(root_len + path_len + 2); + if (!joined) { + free((void*)path); + return DAMACY_OOM; + } + memcpy(joined, uri, root_len); + if (root_len && uri[root_len - 1] != '/') + joined[root_len++] = '/'; + memcpy(joined + root_len, path, path_len + 1); + free((void*)path); + level->uri = joined; + if (member(node, "coordinateTransformations", &value)) + return DAMACY_INVAL; + return read_transform( + value, rank, level->scale_to_reference, level->origin_reference_index); +} + +static enum damacy_status +normalize_levels(struct damacy_ngff_image* image) +{ + struct damacy_ngff_level reference = image->levels[0]; + for (uint32_t i = image->info.level_count; i-- > 0;) { + struct damacy_ngff_level* level = &image->levels[i]; + for (uint8_t d = 0; d < image->info.rank; ++d) { + double scale = level->scale_to_reference[d]; + double origin = level->origin_reference_index[d]; + double base_scale = reference.scale_to_reference[d]; + double base_origin = reference.origin_reference_index[d]; + if (image->info.axes[d].kind != DAMACY_NGFF_SPACE) { + if (scale != base_scale || origin != base_origin) + return DAMACY_UNSUPPORTED; + level->scale_to_reference[d] = 1; + level->origin_reference_index[d] = 0; + } else { + if (i && scale < image->levels[i - 1].scale_to_reference[d]) + return DAMACY_INVAL; + double ratio = scale / base_scale; + double offset = (origin - base_origin) / base_scale + 0.5 * (1 - ratio); + if (!isfinite(ratio) || !isfinite(offset)) + return DAMACY_INVAL; + level->scale_to_reference[d] = ratio; + level->origin_reference_index[d] = offset; + } + } + } + return DAMACY_OK; +} + +enum damacy_status +ngff_parse_group(struct cslice src, + const char* uri, + uint32_t multiscale_index, + uint32_t max_levels, + struct damacy_ngff_image** out) +{ + *out = NULL; + struct json_node root, node, ome, multiscale; + uint64_t format; + if (root_document(src, &root) || member(root, "zarr_format", &node) || + json_as_uint(node, &format) || format != 3 || + member(root, "node_type", &node) || !string_equal(node, "group") || + member(root, "attributes", &node) || member(node, "ome", &ome) || + member(ome, "version", &node)) + return DAMACY_INVAL; + if (!string_equal(node, "0.5")) + return DAMACY_UNSUPPORTED; + if (member(ome, "multiscales", &node) || node.type != JSON_ARRAY) + return DAMACY_INVAL; + const struct json_query select = { .kind = QUERY_INDEX, + .index = multiscale_index }; + if (json_resolve(node.s, &select, 1, &multiscale, NULL)) + return DAMACY_INVAL; + struct damacy_ngff_image* image = calloc(1, sizeof(*image)); + if (!image) + return DAMACY_OOM; + enum damacy_status status = DAMACY_INVAL; + if (member(multiscale, "axes", &node)) + goto Fail; + status = read_axes(node, &image->info); + if (status != DAMACY_OK) + goto Fail; + double group_scale[DAMACY_MAX_RANK] = { 0 }; + double group_origin[DAMACY_MAX_RANK] = { 0 }; + enum json_err error = member(multiscale, "coordinateTransformations", &node); + if (error == JSON_OK) { + status = read_transform(node, image->info.rank, group_scale, group_origin); + if (status != DAMACY_OK) + goto Fail; + } else if (error != JSON_ERR_NOT_FOUND) { + status = DAMACY_INVAL; + goto Fail; + } + status = DAMACY_INVAL; + if (member(multiscale, "datasets", &node)) + goto Fail; + struct json_iter it; + if (elements(node, &it)) + goto Fail; + struct json_node dataset; + while ((error = json_iter_next(&it, &dataset)) == JSON_OK) { + if (image->info.level_count >= max_levels) { + status = DAMACY_BUDGET; + goto Fail; + } + size_t count = (size_t)image->info.level_count + 1; + if (count > SIZE_MAX / sizeof(*image->levels)) { + status = DAMACY_BUDGET; + goto Fail; + } + void* levels = realloc(image->levels, count * sizeof(*image->levels)); + if (!levels) { + status = DAMACY_OOM; + goto Fail; + } + image->levels = levels; + struct damacy_ngff_level* level = &image->levels[image->info.level_count++]; + memset(level, 0, sizeof(*level)); + status = read_level(dataset, uri, image->info.rank, level); + if (status != DAMACY_OK) + goto Fail; + for (uint32_t i = 0; i + 1 < image->info.level_count; ++i) + if (!strcmp(level->uri, image->levels[i].uri)) { + status = DAMACY_INVAL; + goto Fail; + } + } + if (error != JSON_ERR_NOT_FOUND || !image->info.level_count) { + status = DAMACY_INVAL; + goto Fail; + } + status = normalize_levels(image); + if (status != DAMACY_OK) + goto Fail; + image->info.levels = image->levels; + *out = image; + return DAMACY_OK; +Fail: + damacy_ngff_image_destroy(image); + return status; +} + +enum damacy_status +ngff_parse_array(struct cslice src, + struct damacy_ngff_image* image, + uint32_t index) +{ + struct json_node root, names, value; + if (root_document(src, &root)) + return DAMACY_INVAL; + static const char* keys[] = { "zarr_format", "node_type", + "shape", "data_type", + "chunk_grid", "chunk_key_encoding", + "fill_value", "codecs", + "dimension_names" }; + for (size_t i = 0; i < sizeof(keys) / sizeof(*keys); ++i) + if (member(root, keys[i], &value)) + return DAMACY_INVAL; + struct zarr_metadata array; + if (zarr_metadata_parse(src.beg, cslice_len(src), &array)) + return DAMACY_UNSUPPORTED; + if (array.rank != image->info.rank) + return DAMACY_RANK; + struct json_iter it; + if (member(root, "dimension_names", &names) || elements(names, &it)) + return DAMACY_INVAL; + for (uint8_t d = 0; d < array.rank; ++d) { + const char* name = NULL; + if (json_iter_next(&it, &value)) + return DAMACY_INVAL; + enum damacy_status status = read_string(value, &name); + if (status != DAMACY_OK) + return status; + int match = !strcmp(name, image->info.axes[d].name); + free((void*)name); + if (!match || !array.shape[d] || array.shape[d] > INT64_C(4503599627370495)) + return DAMACY_INVAL; + if (index && image->info.axes[d].kind != DAMACY_NGFF_SPACE && + array.shape[d] != (uint64_t)image->levels[0].shape[d]) + return DAMACY_UNSUPPORTED; + if (index && array.shape[d] > (uint64_t)image->levels[index - 1].shape[d]) + return DAMACY_INVAL; + image->levels[index].shape[d] = (int64_t)array.shape[d]; + } + if (json_iter_next(&it, &value) != JSON_ERR_NOT_FOUND) + return DAMACY_INVAL; + if (!index) { + image->dtype = array.dtype; + if (member(root, "data_type", &value)) + return DAMACY_INVAL; + return read_string(value, &image->info.data_type); + } + return array.dtype == image->dtype ? DAMACY_OK : DAMACY_DTYPE; +} + +const struct damacy_ngff_info* +damacy_ngff_image_info(const struct damacy_ngff_image* image) +{ + return image ? &image->info : NULL; +} + +void +damacy_ngff_image_destroy(struct damacy_ngff_image* image) +{ + if (!image) + return; + for (uint8_t d = 0; d < image->info.rank; ++d) { + free((void*)image->info.axes[d].name); + free((void*)image->info.axes[d].unit); + } + for (uint32_t i = 0; i < image->info.level_count; ++i) + free((void*)image->levels[i].uri); + free((void*)image->info.data_type); + free(image->levels); + free(image); +} diff --git a/src/ngff/ngff.h b/src/ngff/ngff.h new file mode 100644 index 00000000..f8b13634 --- /dev/null +++ b/src/ngff/ngff.h @@ -0,0 +1,23 @@ +#pragma once + +#include "damacy_spatial.h" +#include "dtype/dtype.h" +#include "util/json.h" + +struct damacy_ngff_image +{ + struct damacy_ngff_info info; + struct damacy_ngff_level* levels; + enum dtype dtype; +}; + +enum damacy_status +ngff_parse_group(struct cslice src, + const char* uri, + uint32_t multiscale_index, + uint32_t max_levels, + struct damacy_ngff_image** out); +enum damacy_status +ngff_parse_array(struct cslice src, + struct damacy_ngff_image* image, + uint32_t level); diff --git a/src/query/spatial.c b/src/query/spatial.c new file mode 100644 index 00000000..90768fb3 --- /dev/null +++ b/src/query/spatial.c @@ -0,0 +1,273 @@ +#include "damacy_spatial.h" + +#include "damacy_config.h" +#include "ngff/ngff.h" +#include "pipeline/components.h" + +#include +#include +#include +#include + +struct damacy_spatial_resolution +{ + struct damacy_spatial_info info; + struct damacy_sample sample; +}; + +static double +minimum_spacing(const struct damacy_ngff_image* image, + const struct damacy_affine* transform, + uint32_t level) +{ + uint8_t axes[3], rank = 0; + for (uint8_t d = 0; d < image->info.rank; ++d) + if (image->info.axes[d].kind == DAMACY_NGFF_SPACE) + axes[rank++] = d; + double matrix[3][3] = { { 0 } }, scale = 0; + for (uint8_t i = 0; i < rank; ++i) + for (uint8_t j = 0; j < rank; ++j) { + double value = transform->linear[axes[i]][axes[j]] / + image->levels[level].scale_to_reference[axes[i]]; + if (!isfinite(value)) + return NAN; + matrix[i][j] = value; + scale = fmax(scale, fabs(value)); + } + if (!scale) + return 0; + double gram[3][3] = { { 0 } }; + for (uint8_t i = 0; i < rank; ++i) + for (uint8_t j = 0; j < rank; ++j) + for (uint8_t k = 0; k < rank; ++k) + gram[i][j] += (matrix[i][k] / scale) * (matrix[j][k] / scale); + for (int sweep = 0; sweep < 24; ++sweep) { + for (uint8_t p = 0; p < rank; ++p) { + for (uint8_t q = p + 1; q < rank; ++q) { + double cross = gram[p][q]; + if (fabs(cross) <= DBL_EPSILON * sqrt(gram[p][p] * gram[q][q])) + continue; + double tau = (gram[q][q] - gram[p][p]) / (2 * cross); + double t = copysign(1, tau) / (fabs(tau) + hypot(1, tau)); + double c = 1 / hypot(1, t), s = t * c; + gram[p][p] -= t * cross; + gram[q][q] += t * cross; + gram[p][q] = gram[q][p] = 0; + for (uint8_t k = 0; k < rank; ++k) { + if (k == p || k == q) + continue; + double a = gram[k][p], b = gram[k][q]; + gram[k][p] = gram[p][k] = c * a - s * b; + gram[k][q] = gram[q][k] = s * a + c * b; + } + } + } + } + double minimum = gram[0][0]; + for (uint8_t d = 1; d < rank; ++d) + minimum = fmin(minimum, gram[d][d]); + return scale * sqrt(fmax(0, minimum)); +} + +static enum damacy_status +validate_query(const struct damacy_ngff_image* image, + const struct damacy_batch_spec* output, + const struct damacy_spatial_query* query) +{ + if (!image || !output || !query) + return DAMACY_INVAL; + if (output->sample_rank != image->info.rank) + return DAMACY_RANK; + int64_t shape[DAMACY_MAX_RANK + 1], strides[DAMACY_MAX_RANK + 1]; + uint64_t bytes; + enum damacy_status status = batch_spec_layout(output, shape, strides, &bytes); + if (status != DAMACY_OK) + return status; + if (!cast_path_supported(output->dtype, image->dtype)) + return DAMACY_DTYPE; + if (query->level < DAMACY_LEVEL_AUTO || + (query->level >= 0 && (uint32_t)query->level >= image->info.level_count)) + return DAMACY_INVAL; + const struct damacy_sampler* sampler = &query->sampler; + if ((sampler->filter != DAMACY_FILTER_NEAREST && + sampler->filter != DAMACY_FILTER_LINEAR) || + (sampler->boundary != DAMACY_BOUNDARY_ERROR && + sampler->boundary != DAMACY_BOUNDARY_CONSTANT && + sampler->boundary != DAMACY_BOUNDARY_CLAMP) || + !isfinite(sampler->constant_value) || + (sampler->boundary != DAMACY_BOUNDARY_CONSTANT && + sampler->constant_value)) + return DAMACY_INVAL; + for (uint8_t i = 0; i < image->info.rank; ++i) { + double offset = query->output_to_reference.offset[i]; + int space = image->info.axes[i].kind == DAMACY_NGFF_SPACE; + if (!isfinite(offset) || + output->sample_shape[i] > INT64_C(4503599627370495)) + return DAMACY_INVAL; + if (!space && offset != floor(offset)) + return DAMACY_INVAL; + for (uint8_t j = 0; j < image->info.rank; ++j) { + double value = query->output_to_reference.linear[i][j]; + if (!isfinite(value)) + return DAMACY_INVAL; + if ((!space || image->info.axes[j].kind != DAMACY_NGFF_SPACE) && + value != (double)(i == j)) + return DAMACY_INVAL; + } + } + double spacing = minimum_spacing(image, &query->output_to_reference, 0); + if (!isfinite(spacing) || spacing <= 0) + return DAMACY_INVAL; + return DAMACY_OK; +} + +static int64_t +clamp_index(int64_t value, int64_t end) +{ + return value < 0 ? 0 : value > end ? end : value; +} + +static enum damacy_status +resolve_bounds(const struct damacy_ngff_image* image, + struct damacy_spatial_info* info) +{ + info->source_bounds_index.rank = info->read_bounds_index.rank = info->rank; + for (uint8_t i = 0; i < info->rank; ++i) { + long double beg = info->output_to_source.offset[i], end = beg; + for (uint8_t j = 0; j < info->rank; ++j) { + double value = info->output_to_source.linear[i][j]; + long double first = 0.5L * value; + long double last = ((long double)info->output_shape[j] - 0.5L) * value; + beg += fminl(first, last); + end += fmaxl(first, last); + if (value != (double)(i == j)) + info->requires_resampling = 1; + } + double offset = info->output_to_source.offset[i]; + if (offset != floor(offset)) + info->requires_resampling = 1; + if (!isfinite(beg) || !isfinite(end) || beg < -4503599627370495.0L || + end > 4503599627370495.0L) + return DAMACY_INVAL; + struct damacy_interval span; + if (info->sampler.filter == DAMACY_FILTER_NEAREST) { + span.beg = (int64_t)floorl(beg); + span.end = (int64_t)floorl(end) + 1; + } else { + span.beg = (int64_t)floorl(beg - 0.5L); + span.end = (int64_t)ceill(end - 0.5L) + 1; + } + info->source_bounds_index.dims[i] = span; + int64_t size = info->source_shape[i]; + if (span.beg < 0 || span.end > size) { + if (info->sampler.boundary == DAMACY_BOUNDARY_ERROR || + image->info.axes[i].kind != DAMACY_NGFF_SPACE) + return DAMACY_INVAL; + info->requires_resampling = 1; + } + struct damacy_interval read; + if (info->sampler.boundary == DAMACY_BOUNDARY_CLAMP) { + read.beg = clamp_index(span.beg, size - 1); + read.end = clamp_index(span.end - 1, size - 1) + 1; + } else { + read.beg = clamp_index(span.beg, size); + read.end = clamp_index(span.end, size); + } + info->read_bounds_index.dims[i] = read; + } + return DAMACY_OK; +} + +enum damacy_status +damacy_spatial_resolve(const struct damacy_ngff_image* image, + const struct damacy_batch_spec* output, + const struct damacy_spatial_query* query, + struct damacy_spatial_resolution** out) +{ + if (!out) + return DAMACY_INVAL; + *out = NULL; + enum damacy_status status = validate_query(image, output, query); + if (status != DAMACY_OK) + return status; + uint32_t level_index = query->level < 0 ? 0 : (uint32_t)query->level; + if (query->level == DAMACY_LEVEL_AUTO) + for (uint32_t i = 1; i < image->info.level_count; ++i) + if (minimum_spacing(image, &query->output_to_reference, i) >= + 1 - 64 * DBL_EPSILON) + level_index = i; + const struct damacy_ngff_level* level = &image->levels[level_index]; + struct damacy_spatial_resolution* resolution = calloc(1, sizeof(*resolution)); + if (!resolution) + return DAMACY_OOM; + struct damacy_spatial_info* info = &resolution->info; + info->level = level_index; + info->rank = image->info.rank; + info->sampler = query->sampler; + for (uint8_t i = 0; i < info->rank; ++i) { + info->output_shape[i] = output->sample_shape[i]; + info->source_shape[i] = level->shape[i]; + double scale = level->scale_to_reference[i]; + info->output_to_source.offset[i] = (query->output_to_reference.offset[i] - + level->origin_reference_index[i]) / + scale; + for (uint8_t j = 0; j < info->rank; ++j) + info->output_to_source.linear[i][j] = + query->output_to_reference.linear[i][j] / scale; + } + status = resolve_bounds(image, info); + if (status != DAMACY_OK) { + free(resolution); + return status; + } + char* uri = malloc(strlen(level->uri) + 1); + if (!uri) { + free(resolution); + return DAMACY_OOM; + } + strcpy(uri, level->uri); + info->uri = uri; + if (!info->requires_resampling) { + resolution->sample.uri = uri; + resolution->sample.rank = info->rank; + for (uint8_t i = 0; i < info->rank; ++i) + resolution->sample.axes[i] = + (struct damacy_axis_selection){ .kind = DAMACY_AXIS_INTERVAL, + .interval = + info->source_bounds_index.dims[i] }; + } + *out = resolution; + return DAMACY_OK; +} + +const struct damacy_spatial_info* +damacy_spatial_resolution_info( + const struct damacy_spatial_resolution* resolution) +{ + return resolution ? &resolution->info : NULL; +} + +enum damacy_status +damacy_spatial_resolution_sample( + const struct damacy_spatial_resolution* resolution, + struct damacy_sample* out) +{ + if (!out) + return DAMACY_INVAL; + *out = (struct damacy_sample){ 0 }; + if (!resolution) + return DAMACY_INVAL; + if (resolution->info.requires_resampling) + return DAMACY_UNSUPPORTED; + *out = resolution->sample; + return DAMACY_OK; +} + +void +damacy_spatial_resolution_destroy(struct damacy_spatial_resolution* resolution) +{ + if (resolution) { + free((void*)resolution->info.uri); + free(resolution); + } +} diff --git a/src/store/metadata_store_async.c b/src/store/metadata_store_async.c index 80d34541..fbdf7bec 100644 --- a/src/store/metadata_store_async.c +++ b/src/store/metadata_store_async.c @@ -405,6 +405,17 @@ handle_statx_complete(struct metadata_store_async* s, struct metadata_job* job) return; } + if (job->stx.stx_size > job->requested_len || + job->stx.stx_size > UINT32_MAX) { + job->status = DAMACY_BUDGET; + if (job->open_done) { + if (job->fd >= 0) + finish_job(s, job); + else + complete_job(s, job); + } + return; + } job->len = (size_t)job->stx.stx_size; if (job->open_done && job->fd < 0) { complete_job(s, job); @@ -803,7 +814,17 @@ metadata_store_async_read_file(struct metadata_store_async* s, metadata_store_read_cb cb, void* user) { - return post_read(s, key, 0, 0, REQ_READ_FILE, cb, user); + return post_read(s, key, 0, SIZE_MAX, REQ_READ_FILE, cb, user); +} + +int +metadata_store_async_read_file_bounded(struct metadata_store_async* s, + const char* key, + size_t max_bytes, + metadata_store_read_cb cb, + void* user) +{ + return post_read(s, key, 0, max_bytes, REQ_READ_FILE, cb, user); } int diff --git a/src/store/metadata_store_async.h b/src/store/metadata_store_async.h index 2ce3816a..a7f73baa 100644 --- a/src/store/metadata_store_async.h +++ b/src/store/metadata_store_async.h @@ -92,6 +92,12 @@ extern "C" void metadata_store_async_op_latency_stats_reset( struct metadata_store_async* s); + int metadata_store_async_read_file_bounded(struct metadata_store_async* s, + const char* key, + size_t max_bytes, + metadata_store_read_cb cb, + void* user); + int metadata_store_async_read_file(struct metadata_store_async* s, const char* key, metadata_store_read_cb cb, diff --git a/src/util/json.c b/src/util/json.c index b1481c49..d15aaf11 100644 --- a/src/util/json.c +++ b/src/util/json.c @@ -746,3 +746,51 @@ json_str_eq(struct json_node n, const char* lit) size_t span = cslice_len(n.s); return span == L && memcmp(n.s.beg, lit, L) == 0; } + +enum json_err +json_object_iter_init(struct json_node node, struct json_object_iter* it) +{ + if (!it) + return JSON_ERR_INVALID; + if (node.type != JSON_OBJECT) + return JSON_ERR_TYPE; + if (cslice_len(node.s) < 2 || *node.s.beg != '{' || node.s.end[-1] != '}') + return JSON_ERR_PARSE; + *it = (struct json_object_iter){ .remaining = { node.s.beg + 1, node.s.end }, + .first = 1 }; + return JSON_OK; +} + +enum json_err +json_object_iter_next(struct json_object_iter* it, + struct json_node* key, + struct json_node* value) +{ + if (!it || !key || !value) + return JSON_ERR_INVALID; + struct cslice cursor = it->remaining; + skip_ws(&cursor); + if (cs_at_end(cursor)) + return JSON_ERR_PARSE; + if (*cursor.beg == '}') + return JSON_ERR_NOT_FOUND; + if (!it->first) { + if (*cursor.beg != ',') + return JSON_ERR_PARSE; + cursor.beg++; + skip_ws(&cursor); + } + key->type = JSON_STRING; + if (lex_string_span(&cursor, &key->s, &key->flag)) + return JSON_ERR_PARSE; + skip_ws(&cursor); + if (cs_at_end(cursor) || *cursor.beg++ != ':') + return JSON_ERR_PARSE; + skip_ws(&cursor); + enum json_err error = lex_value(&cursor, value); + if (error != JSON_OK) + return error; + it->remaining = cursor; + it->first = 0; + return JSON_OK; +} diff --git a/src/util/json.h b/src/util/json.h index 80914588..34dada90 100644 --- a/src/util/json.h +++ b/src/util/json.h @@ -66,6 +66,18 @@ extern "C" struct cslice s; // span into src (between quotes for strings) }; + struct json_object_iter + { + struct cslice remaining; + int first; + }; + + enum json_err json_object_iter_init(struct json_node node, + struct json_object_iter* it); + enum json_err json_object_iter_next(struct json_object_iter* it, + struct json_node* key, + struct json_node* value); + // Position info for parse failures. offset is the byte offset into the // src cslice where the lexer gave up. Filled in by json_resolve and // json_iter_init when the caller passes a non-NULL err. diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 34f54287..1704b4bc 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -68,6 +68,8 @@ add_damacy_test(test_scheduler scheduler platform log) # requires `uv` on PATH at test time (provided by flake.nix devShell). if(NOT DAMACY_FUZZ) add_damacy_test(test_damacy_plan damacy) + add_damacy_test(test_spatial damacy test_fixture) + set_tests_properties(test_spatial PROPERTIES LABELS tsan TIMEOUT 120) add_damacy_test(test_cpu_pipeline damacy test_fixture Threads::Threads) add_damacy_test(test_cpu_executor cpu_executor) set_tests_properties(test_cpu_pipeline PROPERTIES TIMEOUT 120) diff --git a/tests/test_json.c b/tests/test_json.c index 151fcd67..247646fd 100644 --- a/tests/test_json.c +++ b/tests/test_json.c @@ -324,6 +324,34 @@ test_err_offset_set(void) return 0; } +static int +test_object_iterator(void) +{ + struct json_node root, key, value; + struct json_object_iter it; + EXPECT(json_resolve( + SRC("{\"name\":\"x\",\"values\":[1,2]}"), NULL, 0, &root, NULL) == + JSON_OK); + EXPECT(json_object_iter_init(root, &it) == JSON_OK); + EXPECT(json_object_iter_next(&it, &key, &value) == JSON_OK); + EXPECT(json_str_eq(key, "name") && json_str_eq(value, "x")); + EXPECT(json_object_iter_next(&it, &key, &value) == JSON_OK); + EXPECT(json_str_eq(key, "values") && value.type == JSON_ARRAY); + EXPECT(json_object_iter_next(&it, &key, &value) == JSON_ERR_NOT_FOUND); + EXPECT(json_object_iter_next(&it, &key, &value) == JSON_ERR_NOT_FOUND); + EXPECT(json_object_iter_init(value, &it) == JSON_ERR_TYPE); + EXPECT(json_object_iter_init(root, NULL) == JSON_ERR_INVALID); + EXPECT(json_object_iter_next(NULL, &key, &value) == JSON_ERR_INVALID); + EXPECT(json_resolve(SRC("{\"a\":1,}"), NULL, 0, &root, NULL) == JSON_OK); + EXPECT(json_object_iter_init(root, &it) == JSON_OK); + EXPECT(json_object_iter_next(&it, &key, &value) == JSON_OK); + EXPECT(json_object_iter_next(&it, &key, &value) == JSON_ERR_PARSE); + EXPECT(json_resolve(SRC("{}"), NULL, 0, &root, NULL) == JSON_OK); + EXPECT(json_object_iter_init(root, &it) == JSON_OK); + EXPECT(json_object_iter_next(&it, &key, &value) == JSON_ERR_NOT_FOUND); + return 0; +} + int main(void) { @@ -344,6 +372,7 @@ main(void) RUN(test_primitive_conversions); RUN(test_null_args); RUN(test_err_offset_set); + RUN(test_object_iterator); log_info("all tests passed"); return 0; } diff --git a/tests/test_spatial.c b/tests/test_spatial.c new file mode 100644 index 00000000..d9781d7d --- /dev/null +++ b/tests/test_spatial.c @@ -0,0 +1,335 @@ +#include "damacy_spatial.h" +#include "expect.h" +#include "fixture.h" +#include "ngff/ngff.h" + +#include +#include +#include +#include +#include +#include + +static const char group[] = + "{\"zarr_format\":3,\"node_type\":\"group\",\"attributes\":{\"ome\":{" + "\"version\":\"0.5\",\"multiscales\":[{\"axes\":[" + "{\"name\":\"y\",\"type\":\"space\",\"unit\":\"micrometer\"}," + "{\"name\":\"x\",\"type\":\"space\",\"unit\":\"micrometer\"}]," + "\"coordinateTransformations\":[{\"type\":\"scale\",\"scale\":[3,5]}," + "{\"type\":\"translation\",\"translation\":[80,-20]}],\"datasets\":[" + "{\"path\":\"0\",\"coordinateTransformations\":[" + "{\"type\":\"scale\",\"scale\":[0.5,1]}," + "{\"type\":\"translation\",\"translation\":[10,20]}]}," + "{\"path\":\"1\",\"coordinateTransformations\":[" + "{\"type\":\"scale\",\"scale\":[1,2]}," + "{\"type\":\"translation\",\"translation\":[10.25,20.5]}]}," + "{\"path\":\"2\",\"coordinateTransformations\":[" + "{\"type\":\"scale\",\"scale\":[2,4]}," + "{\"type\":\"translation\",\"translation\":[10,20]}]}]}]}}}"; + +static void +array_json(char* dst, size_t capacity, int size) +{ + snprintf(dst, + capacity, + "{\"zarr_format\":3,\"node_type\":\"array\",\"shape\":[%d,%d]," + "\"dimension_names\":[\"y\",\"x\"],\"data_type\":\"uint16\"," + "\"chunk_grid\":{\"name\":\"regular\",\"configuration\":{\"chunk_" + "shape\":[8,8]}}," + "\"chunk_key_encoding\":{\"name\":\"default\",\"configuration\":{" + "\"separator\":\"/\"}}," + "\"fill_value\":0,\"codecs\":[{\"name\":\"bytes\",\"configuration\":" + "{\"endian\":\"little\"}}]}", + size, + size); +} + +static struct cslice +text_slice(const char* text) +{ + return (struct cslice){ text, text + strlen(text) }; +} + +static int +image_create(struct damacy_ngff_image** image) +{ + EXPECT(ngff_parse_group(text_slice(group), "volume", 0, 3, image) == + DAMACY_OK); + for (uint32_t i = 0; i < 3; ++i) { + char array[1024]; + array_json(array, sizeof(array), 64 >> i); + EXPECT(ngff_parse_array(text_slice(array), *image, i) == DAMACY_OK); + } + return 0; +} + +static struct damacy_spatial_query +identity_query(void) +{ + return (struct damacy_spatial_query){ + .output_to_reference = { .linear = { { 1, 0 }, { 0, 1 } } }, + .sampler = { .filter = DAMACY_FILTER_NEAREST, + .boundary = DAMACY_BOUNDARY_ERROR }, + .level = DAMACY_LEVEL_AUTO + }; +} + +static struct damacy_batch_spec output = { .dtype = DAMACY_F32, + .sample_shape = { 4, 4 }, + .sample_rank = 2, + .samples_per_batch = 1 }; + +static int +test_corner_conversion_and_copy(void) +{ + struct damacy_ngff_image* image; + EXPECT(!image_create(&image)); + const struct damacy_ngff_info* metadata = damacy_ngff_image_info(image); + EXPECT(metadata->rank == 2 && metadata->level_count == 3); + EXPECT(!strcmp(metadata->data_type, "uint16")); + EXPECT(!strcmp(metadata->axes[0].unit, "micrometer")); + EXPECT(!strcmp(metadata->levels[1].uri, "volume/1")); + for (uint8_t d = 0; d < 2; ++d) { + EXPECT(metadata->levels[0].scale_to_reference[d] == 1); + EXPECT(metadata->levels[0].origin_reference_index[d] == 0); + EXPECT(metadata->levels[1].scale_to_reference[d] == 2); + EXPECT(metadata->levels[1].origin_reference_index[d] == 0); + EXPECT(metadata->levels[2].origin_reference_index[d] == -1.5); + } + struct damacy_spatial_query query = identity_query(); + struct damacy_spatial_resolution* resolved; + EXPECT(damacy_spatial_resolve(image, &output, &query, &resolved) == + DAMACY_OK); + EXPECT(damacy_spatial_resolution_info(resolved)->level == 0); + struct damacy_sample sample; + EXPECT(damacy_spatial_resolution_sample(resolved, &sample) == DAMACY_OK); + EXPECT(sample.axes[0].interval.beg == 0 && sample.axes[0].interval.end == 4); + damacy_spatial_resolution_destroy(resolved); + query.output_to_reference.linear[0][0] = 2; + query.output_to_reference.linear[1][1] = 2; + query.output_to_reference.offset[0] = 4; + query.output_to_reference.offset[1] = 6; + query.sampler.filter = DAMACY_FILTER_LINEAR; + EXPECT(damacy_spatial_resolve(image, &output, &query, &resolved) == + DAMACY_OK); + const struct damacy_spatial_info* info = + damacy_spatial_resolution_info(resolved); + EXPECT(info->level == 1 && !info->requires_resampling); + damacy_ngff_image_destroy(image); + memset(&query, 0, sizeof(query)); + EXPECT(damacy_spatial_resolution_sample(resolved, &sample) == DAMACY_OK); + EXPECT(!strcmp(sample.uri, "volume/1")); + EXPECT(sample.rank == 2 && sample.axes[0].kind == DAMACY_AXIS_INTERVAL); + EXPECT(sample.axes[0].interval.beg == 2 && sample.axes[0].interval.end == 6); + EXPECT(sample.axes[1].interval.beg == 3 && sample.axes[1].interval.end == 7); + damacy_spatial_resolution_destroy(resolved); + return 0; +} + +static int +test_scale_rotation_and_shear(void) +{ + struct damacy_ngff_image* image; + EXPECT(!image_create(&image)); + struct damacy_spatial_query query = identity_query(); + query.output_to_reference.offset[0] = query.output_to_reference.offset[1] = + 12; + double s = sqrt(2.0); + query.output_to_reference.linear[0][0] = s; + query.output_to_reference.linear[0][1] = -s; + query.output_to_reference.linear[1][0] = s; + query.output_to_reference.linear[1][1] = s; + struct damacy_spatial_resolution* resolved; + EXPECT(damacy_spatial_resolve(image, &output, &query, &resolved) == + DAMACY_OK); + const struct damacy_spatial_info* info = + damacy_spatial_resolution_info(resolved); + EXPECT(info->level == 1 && info->requires_resampling); + EXPECT(fabs(info->output_to_source.linear[0][0] - s / 2) < 1e-15); + struct damacy_sample sample; + EXPECT(damacy_spatial_resolution_sample(resolved, &sample) == + DAMACY_UNSUPPORTED); + EXPECT(!sample.uri && !sample.rank); + damacy_spatial_resolution_destroy(resolved); + query.output_to_reference.linear[0][0] = 3; + query.output_to_reference.linear[0][1] = 2; + query.output_to_reference.linear[1][0] = 2; + query.output_to_reference.linear[1][1] = 3; + EXPECT(damacy_spatial_resolve(image, &output, &query, &resolved) == + DAMACY_OK); + EXPECT(damacy_spatial_resolution_info(resolved)->level == 0); + damacy_spatial_resolution_destroy(resolved); + query = identity_query(); + query.output_to_reference.linear[0][0] = 4; + query.output_to_reference.linear[1][1] = 4; + EXPECT(damacy_spatial_resolve(image, &output, &query, &resolved) == + DAMACY_OK); + info = damacy_spatial_resolution_info(resolved); + EXPECT(info->level == 2 && info->requires_resampling); + EXPECT(info->output_to_source.offset[0] == 0.375); + damacy_spatial_resolution_destroy(resolved); + query = identity_query(); + query.output_to_reference.linear[0][0] = 0.5; + query.output_to_reference.linear[1][1] = 0.5; + EXPECT(damacy_spatial_resolve(image, &output, &query, &resolved) == + DAMACY_OK); + EXPECT(damacy_spatial_resolution_info(resolved)->level == 0); + damacy_spatial_resolution_destroy(resolved); + damacy_ngff_image_destroy(image); + return 0; +} + +static int +test_sampler_bounds(void) +{ + struct damacy_ngff_image* image; + EXPECT(!image_create(&image)); + struct damacy_spatial_query query = identity_query(); + query.output_to_reference.offset[0] = -0.25; + struct damacy_spatial_resolution* resolved; + EXPECT(damacy_spatial_resolve(image, &output, &query, &resolved) == + DAMACY_OK); + const struct damacy_spatial_info* info = + damacy_spatial_resolution_info(resolved); + EXPECT(info->source_bounds_index.dims[0].beg == 0); + EXPECT(info->source_bounds_index.dims[0].end == 4); + damacy_spatial_resolution_destroy(resolved); + query.sampler.filter = DAMACY_FILTER_LINEAR; + EXPECT(damacy_spatial_resolve(image, &output, &query, &resolved) == + DAMACY_INVAL); + EXPECT(!resolved); + query.sampler.boundary = DAMACY_BOUNDARY_CONSTANT; + query.sampler.constant_value = 17; + EXPECT(damacy_spatial_resolve(image, &output, &query, &resolved) == + DAMACY_OK); + info = damacy_spatial_resolution_info(resolved); + EXPECT(info->source_bounds_index.dims[0].beg == -1); + EXPECT(info->source_bounds_index.dims[0].end == 4); + EXPECT(info->read_bounds_index.dims[0].beg == 0); + EXPECT(info->read_bounds_index.dims[0].end == 4); + damacy_spatial_resolution_destroy(resolved); + query.output_to_reference.offset[0] = -20; + EXPECT(damacy_spatial_resolve(image, &output, &query, &resolved) == + DAMACY_OK); + info = damacy_spatial_resolution_info(resolved); + EXPECT(info->read_bounds_index.dims[0].beg == 0); + EXPECT(info->read_bounds_index.dims[0].end == 0); + damacy_spatial_resolution_destroy(resolved); + query.sampler.boundary = DAMACY_BOUNDARY_CLAMP; + query.sampler.constant_value = 0; + EXPECT(damacy_spatial_resolve(image, &output, &query, &resolved) == + DAMACY_OK); + info = damacy_spatial_resolution_info(resolved); + EXPECT(info->read_bounds_index.dims[0].beg == 0); + EXPECT(info->read_bounds_index.dims[0].end == 1); + damacy_spatial_resolution_destroy(resolved); + damacy_ngff_image_destroy(image); + return 0; +} + +static int +test_invalid_queries(void) +{ + struct damacy_ngff_image* image; + EXPECT(!image_create(&image)); + struct damacy_spatial_query query = identity_query(); + struct damacy_spatial_resolution* resolved; + double invalid[] = { 0, NAN, INFINITY, 1e100 }; + for (size_t i = 0; i < sizeof(invalid) / sizeof(*invalid); ++i) { + query.output_to_reference.linear[0][0] = invalid[i]; + EXPECT(damacy_spatial_resolve(image, &output, &query, &resolved) == + DAMACY_INVAL); + EXPECT(!resolved); + } + query = identity_query(); + query.level = -2; + EXPECT(damacy_spatial_resolve(image, &output, &query, &resolved) == + DAMACY_INVAL); + query.level = 3; + EXPECT(damacy_spatial_resolve(image, &output, &query, &resolved) == + DAMACY_INVAL); + query.level = 0; + query.sampler.constant_value = 1; + EXPECT(damacy_spatial_resolve(image, &output, &query, &resolved) == + DAMACY_INVAL); + query.sampler.constant_value = 0; + query.sampler.filter = 0; + EXPECT(damacy_spatial_resolve(image, &output, &query, &resolved) == + DAMACY_INVAL); + query = identity_query(); + struct damacy_batch_spec wrong = output; + wrong.sample_rank = 3; + EXPECT(damacy_spatial_resolve(image, &wrong, &query, &resolved) == + DAMACY_RANK); + damacy_ngff_image_destroy(image); + return 0; +} + +static int +test_metadata_limits_and_load(void) +{ + struct damacy_ngff_image* image = NULL; + EXPECT(ngff_parse_group(text_slice(group), "v", 0, 2, &image) == + DAMACY_BUDGET); + EXPECT(!image); + EXPECT(ngff_parse_group(text_slice(group), "v", 1, 3, &image) == + DAMACY_INVAL); + char root[4096]; + const char* scratch = getenv("TMPDIR"); + snprintf( + root, sizeof(root), "%s/damacy-spatial-XXXXXX", scratch ? scratch : "/tmp"); + EXPECT(mkdtemp(root)); + char path[8192], array[1024]; + snprintf(path, sizeof(path), "%s/zarr.json", root); + EXPECT(!fixture_write_file(path, group)); + size_t metadata_bytes = strlen(group); + for (int i = 0; i < 3; ++i) { + snprintf(path, sizeof(path), "%s/%d", root, i); + EXPECT(!mkdir(path, 0700)); + snprintf(path, sizeof(path), "%s/%d/zarr.json", root, i); + array_json(array, sizeof(array), 64 >> i); + EXPECT(!fixture_write_file(path, array)); + metadata_bytes += strlen(array); + } + struct damacy_metadata_reader* reader; + EXPECT(damacy_file_metadata_reader_create(2, NULL, &reader) == DAMACY_OK); + struct damacy_ngff_limits limits = { .max_levels = 3, + .max_metadata_bytes = metadata_bytes }; + EXPECT(damacy_ngff_image_load(reader, root, 0, &limits, &image) == DAMACY_OK); + EXPECT(damacy_ngff_image_info(image)->levels[2].shape[0] == 16); + damacy_ngff_image_destroy(image); + limits.max_metadata_bytes--; + EXPECT(damacy_ngff_image_load(reader, root, 0, &limits, &image) == + DAMACY_BUDGET); + EXPECT(!image); + limits.max_metadata_bytes = 1; + EXPECT(damacy_ngff_image_load(reader, root, 0, &limits, &image) == + DAMACY_BUDGET); + EXPECT(!image); + limits.max_metadata_bytes = metadata_bytes; + EXPECT(damacy_ngff_image_load( + reader, "/damacy-missing-image", 0, &limits, &image) == + DAMACY_NOTFOUND); + damacy_metadata_reader_destroy(reader); + for (int i = 0; i < 3; ++i) { + snprintf(path, sizeof(path), "%s/%d/zarr.json", root, i); + EXPECT(!unlink(path)); + snprintf(path, sizeof(path), "%s/%d", root, i); + EXPECT(!rmdir(path)); + } + snprintf(path, sizeof(path), "%s/zarr.json", root); + EXPECT(!unlink(path)); + EXPECT(!rmdir(root)); + return 0; +} + +int +main(void) +{ + RUN(test_corner_conversion_and_copy); + RUN(test_scale_rotation_and_shear); + RUN(test_sampler_bounds); + RUN(test_invalid_queries); + RUN(test_metadata_limits_and_load); + return 0; +} From 65ec48ae6f737be921dd6be130225908bfe8a0c3 Mon Sep 17 00:00:00 2001 From: Nathan Clack Date: Tue, 15 Sep 2026 01:03:17 +0000 Subject: [PATCH 02/10] Validate spatial matrix rank robustly --- docs/spatial.md | 4 +- python/tests/test_spatial.py | 35 ++++++++++ src/query/spatial.c | 124 +++++++++++++++++++++++++---------- tests/test_spatial.c | 51 ++++++++++++++ 4 files changed, 180 insertions(+), 34 deletions(-) diff --git a/docs/spatial.md b/docs/spatial.md index 03c53bb5..4f06cfee 100644 --- a/docs/spatial.md +++ b/docs/spatial.md @@ -106,7 +106,9 @@ require identity rows/columns and integral offsets, and must remain in bounds. Combining spatial resampling with indexed time/channel selections is a later extension; use `IndexQuery` for current arbitrary index selections. -The spatial linear map must be finite and numerically nonsingular. Shapes and +The spatial linear map must be finite and numerically nonsingular. Validation +normalizes each spatial row by its largest coefficient and eliminates using +the largest remaining pivot. Pivots no larger than `64 * DBL_EPSILON` fail validation. Shapes and transformed sample coordinates are limited to `2**52 - 1` index units so half-voxel centers remain representable. Rank, dtype conversion, output size, and every sampler parameter are validated before a resolution is returned. diff --git a/python/tests/test_spatial.py b/python/tests/test_spatial.py index bc3194c4..c334e07d 100644 --- a/python/tests/test_spatial.py +++ b/python/tests/test_spatial.py @@ -657,3 +657,38 @@ def test_duplicate_array_fields_are_rejected(tmp_path, key, value): path.write_text(text[:-1] + f', "{key}": {value}' + "}") with pytest.raises(damacy.InvalidArgument): load(root) + + +def test_collapsed_3d_transforms_are_rejected(tmp_path): + root = tmp_path / "image" + write_image(root, shape=(16, 16, 16)) + resolver = damacy.SpatialResolver(load(root), damacy.BatchSpec(1, (2, 2, 2))) + rng = np.random.default_rng(832) + for _ in range(128): + matrix = rng.normal(size=(3, 3)) + matrix[2] = matrix[0] + matrix[1] + q = damacy.SpatialQuery( + output_to_reference=np.column_stack((matrix, [4, 4, 4])), + sampler=damacy.Sampler(boundary="constant"), + ) + with pytest.raises(damacy.InvalidArgument): + resolver.resolve(q) + + +@pytest.mark.parametrize("spacing,level", [(0.5, 0), (3.0, 1)]) +def test_level_selection_with_large_scale_difference(tmp_path, spacing, level): + root = tmp_path / "image" + write_image(root) + resolver = damacy.SpatialResolver(load(root), damacy.BatchSpec(1, (2, 2))) + large = 2**28 + matrix = ( + (large + spacing / 2, large - spacing / 2, 0), + (large - spacing / 2, large + spacing / 2, 0), + ) + resolved = resolver.resolve( + damacy.SpatialQuery( + output_to_reference=matrix, + sampler=damacy.Sampler(boundary="constant"), + ) + ) + assert resolved.level == level diff --git a/src/query/spatial.c b/src/query/spatial.c index 90768fb3..a110586c 100644 --- a/src/query/spatial.c +++ b/src/query/spatial.c @@ -15,58 +15,117 @@ struct damacy_spatial_resolution struct damacy_sample sample; }; -static double -minimum_spacing(const struct damacy_ngff_image* image, - const struct damacy_affine* transform, - uint32_t level) +static uint8_t +spatial_matrix(const struct damacy_ngff_image* image, + const struct damacy_affine* transform, + uint32_t level, + double matrix[3][3]) { uint8_t axes[3], rank = 0; for (uint8_t d = 0; d < image->info.rank; ++d) if (image->info.axes[d].kind == DAMACY_NGFF_SPACE) axes[rank++] = d; - double matrix[3][3] = { { 0 } }, scale = 0; for (uint8_t i = 0; i < rank; ++i) - for (uint8_t j = 0; j < rank; ++j) { - double value = transform->linear[axes[i]][axes[j]] / + for (uint8_t j = 0; j < rank; ++j) + matrix[i][j] = transform->linear[axes[i]][axes[j]] / image->levels[level].scale_to_reference[axes[i]]; - if (!isfinite(value)) - return NAN; - matrix[i][j] = value; - scale = fmax(scale, fabs(value)); + return rank; +} + +static int +linear_is_invertible(const struct damacy_ngff_image* image, + const struct damacy_affine* transform) +{ + double matrix[3][3]; + uint8_t rank = spatial_matrix(image, transform, 0, matrix); + for (uint8_t i = 0; i < rank; ++i) { + double scale = 0; + for (uint8_t j = 0; j < rank; ++j) + scale = fmax(scale, fabs(matrix[i][j])); + if (!scale) + return 0; + for (uint8_t j = 0; j < rank; ++j) + matrix[i][j] /= scale; + } + for (uint8_t j = 0; j < rank; ++j) { + uint8_t pivot_row = j, pivot_column = j; + for (uint8_t i = j; i < rank; ++i) + for (uint8_t k = j; k < rank; ++k) + if (fabs(matrix[i][k]) > fabs(matrix[pivot_row][pivot_column])) { + pivot_row = i; + pivot_column = k; + } + if (fabs(matrix[pivot_row][pivot_column]) <= 64 * DBL_EPSILON) + return 0; + for (uint8_t k = j; k < rank; ++k) { + double value = matrix[j][k]; + matrix[j][k] = matrix[pivot_row][k]; + matrix[pivot_row][k] = value; + } + for (uint8_t i = j; i < rank; ++i) { + double value = matrix[i][j]; + matrix[i][j] = matrix[i][pivot_column]; + matrix[i][pivot_column] = value; } + for (uint8_t i = j + 1; i < rank; ++i) { + double factor = matrix[i][j] / matrix[j][j]; + for (uint8_t k = j + 1; k < rank; ++k) + matrix[i][k] -= factor * matrix[j][k]; + } + } + return 1; +} + +static double +minimum_spacing(const struct damacy_ngff_image* image, + const struct damacy_affine* transform, + uint32_t level) +{ + double matrix[3][3]; + uint8_t rank = spatial_matrix(image, transform, level, matrix); + double scale = 0; + for (uint8_t i = 0; i < rank; ++i) + for (uint8_t j = 0; j < rank; ++j) + scale = fmax(scale, fabs(matrix[i][j])); if (!scale) return 0; - double gram[3][3] = { { 0 } }; for (uint8_t i = 0; i < rank; ++i) for (uint8_t j = 0; j < rank; ++j) - for (uint8_t k = 0; k < rank; ++k) - gram[i][j] += (matrix[i][k] / scale) * (matrix[j][k] / scale); - for (int sweep = 0; sweep < 24; ++sweep) { + matrix[i][j] /= scale; + for (int sweep = 0; sweep < 32; ++sweep) { + int changed = 0; for (uint8_t p = 0; p < rank; ++p) { for (uint8_t q = p + 1; q < rank; ++q) { - double cross = gram[p][q]; - if (fabs(cross) <= DBL_EPSILON * sqrt(gram[p][p] * gram[q][q])) + double a = 0, b = 0, cross = 0; + for (uint8_t i = 0; i < rank; ++i) { + a += matrix[i][p] * matrix[i][p]; + b += matrix[i][q] * matrix[i][q]; + cross += matrix[i][p] * matrix[i][q]; + } + if (!cross || fabs(cross) <= DBL_EPSILON * sqrt(a) * sqrt(b)) continue; - double tau = (gram[q][q] - gram[p][p]) / (2 * cross); + double tau = (b - a) / (2 * cross); double t = copysign(1, tau) / (fabs(tau) + hypot(1, tau)); double c = 1 / hypot(1, t), s = t * c; - gram[p][p] -= t * cross; - gram[q][q] += t * cross; - gram[p][q] = gram[q][p] = 0; - for (uint8_t k = 0; k < rank; ++k) { - if (k == p || k == q) - continue; - double a = gram[k][p], b = gram[k][q]; - gram[k][p] = gram[p][k] = c * a - s * b; - gram[k][q] = gram[q][k] = s * a + c * b; + for (uint8_t i = 0; i < rank; ++i) { + double first = matrix[i][p], second = matrix[i][q]; + matrix[i][p] = c * first - s * second; + matrix[i][q] = s * first + c * second; } + changed = 1; } } + if (!changed) + break; + } + double minimum = INFINITY; + for (uint8_t j = 0; j < rank; ++j) { + double length = 0; + for (uint8_t i = 0; i < rank; ++i) + length = hypot(length, matrix[i][j]); + minimum = fmin(minimum, length); } - double minimum = gram[0][0]; - for (uint8_t d = 1; d < rank; ++d) - minimum = fmin(minimum, gram[d][d]); - return scale * sqrt(fmax(0, minimum)); + return scale * minimum; } static enum damacy_status @@ -115,8 +174,7 @@ validate_query(const struct damacy_ngff_image* image, return DAMACY_INVAL; } } - double spacing = minimum_spacing(image, &query->output_to_reference, 0); - if (!isfinite(spacing) || spacing <= 0) + if (!linear_is_invertible(image, &query->output_to_reference)) return DAMACY_INVAL; return DAMACY_OK; } diff --git a/tests/test_spatial.c b/tests/test_spatial.c index d9781d7d..d0e67975 100644 --- a/tests/test_spatial.c +++ b/tests/test_spatial.c @@ -323,6 +323,56 @@ test_metadata_limits_and_load(void) return 0; } +static int +test_collapsed_volume(void) +{ + struct damacy_ngff_level level = { .uri = "volume/0", + .shape = { 8, 8, 8 }, + .scale_to_reference = { 1, 1, 1 } }; + struct damacy_ngff_image image = { + .info = { .rank = 3, + .data_type = "uint16", + .axes = { { .name = "z", .kind = DAMACY_NGFF_SPACE }, + { .name = "y", .kind = DAMACY_NGFF_SPACE }, + { .name = "x", .kind = DAMACY_NGFF_SPACE } }, + .level_count = 1, + .levels = &level }, + .levels = &level, + .dtype = dtype_u16 + }; + struct damacy_batch_spec shape = { .dtype = DAMACY_F32, + .sample_rank = 3, + .sample_shape = { 2, 2, 2 }, + .samples_per_batch = 1 }; + struct damacy_spatial_query query = { + .output_to_reference = { .linear = { { 0.7854589457317591, + 0.4315455905112043, + 0.5705067473246884 }, + { 0.2148316327537341, + 0.7753842829888732, + -1.494658594462818 }, + { 1.0002905784854932, + 1.2069298735000775, + -0.9241518471381295 } }, + .offset = { 4, 4, 4 } }, + .sampler = { .filter = DAMACY_FILTER_NEAREST, + .boundary = DAMACY_BOUNDARY_CONSTANT }, + .level = DAMACY_LEVEL_AUTO + }; + struct damacy_spatial_resolution* resolved = NULL; + EXPECT(damacy_spatial_resolve(&image, &shape, &query, &resolved) == + DAMACY_INVAL); + EXPECT(!resolved); + query.output_to_reference = (struct damacy_affine){ + .linear = { { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } } + }; + EXPECT(damacy_spatial_resolve(&image, &shape, &query, &resolved) == + DAMACY_OK); + EXPECT(!damacy_spatial_resolution_info(resolved)->requires_resampling); + damacy_spatial_resolution_destroy(resolved); + return 0; +} + int main(void) { @@ -330,6 +380,7 @@ main(void) RUN(test_scale_rotation_and_shear); RUN(test_sampler_bounds); RUN(test_invalid_queries); + RUN(test_collapsed_volume); RUN(test_metadata_limits_and_load); return 0; } From 2823100de8180a23e43ff238f36cd0d4a05822b9 Mon Sep 17 00:00:00 2001 From: Nathan Clack Date: Wed, 16 Sep 2026 22:48:59 +0000 Subject: [PATCH 03/10] Simplify spatial query resolution --- docs/api.md | 2 - docs/pipeline.md | 2 +- docs/spatial.md | 63 ++++++++----- python/damacy/__init__.py | 2 - python/damacy/_native.pyi | 6 +- python/damacy/_spatial.c | 126 +++++++------------------ python/damacy/_spatial.py | 103 ++++++++++---------- python/tests/test_spatial.py | 160 +++++++++++++++++++++---------- src/damacy_spatial.h | 12 +-- src/ngff/metadata.c | 25 +---- src/query/spatial.c | 114 +++++++++-------------- tests/test_spatial.c | 176 +++++++++++++++++++---------------- 12 files changed, 381 insertions(+), 410 deletions(-) diff --git a/docs/api.md b/docs/api.md index 90b13d81..104b636c 100644 --- a/docs/api.md +++ b/docs/api.md @@ -36,8 +36,6 @@ The full public surface of the `damacy` package. ::: damacy.NgffLevel -::: damacy.SpatialResolver - ::: damacy.SpatialQuery ::: damacy.ResolvedSpatialQuery diff --git a/docs/pipeline.md b/docs/pipeline.md index bce339e1..0de80168 100644 --- a/docs/pipeline.md +++ b/docs/pipeline.md @@ -54,7 +54,7 @@ with damacy.Pipeline( ``` The paths in `Sample` and `IndexQuery` must name little-endian numeric Zarr v3 -**arrays**. Use an [NGFF image and spatial resolver](spatial.md) to select a +**arrays**. Use an [NGFF image](spatial.md) to select a level from a group before pushing an aligned crop. The examples assume decoded chunks no larger than 2 MiB. diff --git a/docs/spatial.md b/docs/spatial.md index 4f06cfee..0d365932 100644 --- a/docs/spatial.md +++ b/docs/spatial.md @@ -6,7 +6,7 @@ Metadata loading and query resolution are separate from chunk planning and decod ```mermaid flowchart LR reader[Metadata reader] --> image[Loaded NGFF image] - image --> resolver[Spatial resolver] + image --> resolver[Resolve query] query[Transform and sampler] --> resolver output[Fixed output shape] --> resolver resolver --> resolved[Source array, transform, and bounds] @@ -39,7 +39,6 @@ image = damacy.NgffImage( limits=damacy.NgffLimits(max_levels=16, max_metadata_bytes=4 << 20), ) output = damacy.BatchSpec(samples=1, shape=(16, 64, 64), dtype="f32") -resolver = damacy.SpatialResolver(image=image, output=output) query = damacy.SpatialQuery( output_to_reference=( (1, 0, 0, 8), @@ -48,7 +47,7 @@ query = damacy.SpatialQuery( ), sampler=damacy.Sampler(filter="nearest", boundary="error"), ) -resolved = resolver.resolve(query) +resolved = image.resolve(query, shape=output.shape) planner = damacy.ChunkPlanner( metadata=damacy.ZarrMetadata( @@ -73,9 +72,10 @@ with damacy.Pipeline( ``` `NgffImage` loads the group's `zarr.json` and each level's array metadata once. -It owns an immutable description. `SpatialResolver.resolve()` performs no I/O -and can be called from multiple threads. A result owns its source URI and -geometry, so it remains usable after the image and resolver are released. +It owns an immutable description. `NgffImage.resolve(query, shape=...)` performs +no I/O and can be called from multiple threads. A result owns its source URI +and geometry, so it remains usable after the image is released. Python results +contain ordinary immutable values and can be pickled for transfer between processes. Source metadata and data must remain unchanged while the image is in use. Loading is synchronous and uses the supplied reader exclusively for the @@ -99,7 +99,7 @@ reference_corner = linear * output_corner + offset ``` Rows and columns follow the metadata's array axis order, including time and -channel dimensions. The `BatchSpec` supplies the rank and fixed output shape. +channel dimensions. The supplied output shape determines the rank and fixed output grid. Output centers are the grid points `(j[0] + 0.5, ..., j[N-1] + 0.5)`. Only spatial axes may mix, rotate, reflect, or scale. Time and channel axes require identity rows/columns and integral offsets, and must remain in bounds. @@ -110,8 +110,10 @@ The spatial linear map must be finite and numerically nonsingular. Validation normalizes each spatial row by its largest coefficient and eliminates using the largest remaining pivot. Pivots no larger than `64 * DBL_EPSILON` fail validation. Shapes and transformed sample coordinates are limited to `2**52 - 1` index units so -half-voxel centers remain representable. Rank, dtype conversion, output size, -and every sampler parameter are validated before a resolution is returned. +half-voxel centers remain representable. Rank, shape, transforms, and sampler +parameters are checked before a resolution is returned. Batch sizing and dtype +conversion are validated by the downstream pipeline. Resolution allocates only +its small description, so the shape need not fit an output buffer at this stage. ## NGFF interpretation @@ -124,12 +126,19 @@ and a common dtype. Custom/unspecified axis types and transforms stored in external arrays are currently unsupported. Dataset transforms must have one positive scale and an optional following -translation. Optional transforms shared by all levels follow the same rule. +translation. Transforms shared by all levels do not affect reference-level +coordinates and are left uninterpreted. Levels must have nondecreasing scales and nonincreasing shapes along each axis; nonspatial transforms and shapes must be unchanged. Paths must be relative child paths. Axis units are retained as declared; no physical-unit conversion or inference is needed for reference-level queries. +Parsing checks the fields needed for array layout and coordinate interpretation, +including duplicate consumed fields and finite, correctly sized transforms. +Unselected multiscale entries and additional attributes are left uninterpreted. +Unused subtrees and trailing content may remain unexamined; these calls do not +perform whole-document JSON validation. + NGFF uses center-origin coordinates, as described in its [coordinate convention](https://ngff.openmicroscopy.org/rfc/5/#coordinate-convention). The adapter derives a corner-origin map once while loading metadata. If the @@ -176,8 +185,8 @@ The current copy path requires the resolved linear map to equal identity, its offsets to be integral, and its source bounds to be inside the array. Fractional values are not rounded into a crop. Small floating-point differences in metadata can therefore make a result require resampling. `requires_resampling` -reports this, and `as_sample()` checks it before returning an interval query. -The pipeline still checks the resulting shape against its own `BatchSpec`. +reports this, and `Pipeline.push()` rejects such results until a resampler is +available. The pipeline checks the resulting shape against its own `BatchSpec`. ## Sampler and bounds @@ -220,14 +229,11 @@ have matching output geometry and an executor: ```c struct damacy_metadata_reader* reader = NULL; struct damacy_ngff_image* image = NULL; -struct damacy_spatial_resolution* resolved = NULL; +struct damacy_spatial_resolution resolved = {0}; struct damacy_ngff_limits limits = { .max_levels = 16, .max_metadata_bytes = 4 << 20 }; -struct damacy_batch_spec output = { - .dtype = DAMACY_F32, .sample_shape = {64, 64}, - .sample_rank = 2, .samples_per_batch = 1 -}; +const int64_t output_shape[] = {64, 64}; struct damacy_spatial_query query = { .output_to_reference = { .linear = {{1, 0}, {0, 1}}, .offset = {16, 24} @@ -243,15 +249,15 @@ enum damacy_status status = if (status == DAMACY_OK) status = damacy_ngff_image_load(reader, "/data/image.zarr", 0, &limits, &image); if (status == DAMACY_OK) - status = damacy_spatial_resolve(image, &output, &query, &resolved); + status = damacy_spatial_resolve(image, &query, 2, output_shape, &resolved); if (status == DAMACY_OK) - status = damacy_spatial_resolution_sample(resolved, &sample); + status = damacy_spatial_resolution_sample(&resolved, &sample); if (status == DAMACY_OK) { struct damacy_push_result pushed = damacy_push( pipeline, (struct damacy_sample_slice){.beg = &sample, .end = &sample + 1}); status = pushed.status; } -damacy_spatial_resolution_destroy(resolved); +damacy_spatial_resolution_clear(&resolved); damacy_ngff_image_destroy(image); damacy_metadata_reader_destroy(reader); ``` @@ -259,11 +265,18 @@ damacy_metadata_reader_destroy(reader); A production caller must retry the unconsumed suffix when `damacy_push` returns `DAMACY_AGAIN`, keeping the sample's owner alive until acceptance or abandonment. The example releases it on return. Accepted pushes copy the URI and selectors. -`damacy_spatial_resolution_sample()` borrows its URI from `resolved` and zeroes -its output on failure. Image and resolution info accessors return borrowed, -immutable views; their owners must remain alive while reading those views. -Resolution copies the information it needs from the image. Destruction accepts -null pointers. C affine entries outside the configured rank are unused. +The C result is a caller-owned `damacy_spatial_resolution` value with directly +readable fields. It owns its URI independently of the image. Treat the fields as +read-only and call `damacy_spatial_resolution_clear()` before reusing or discarding +the result. A plain struct copy borrows the same URI; clear only the owning +value. Clear frees the URI and zeros the value; repeating it is safe. +Resolution zeros its output on failure. Zero-initialized results can be cleared. + +`damacy_spatial_resolution_sample()` is the compatibility bridge to the existing +C push API. Its sample borrows the result's URI and its output is zeroed on +failure. Keep the result alive until the sample is accepted or abandoned. +The image info accessor returns a borrowed, immutable view. C affine entries +outside the configured rank are unused. `max_metadata_bytes` bounds the sum of input JSON file sizes; files exceeding the remaining budget are rejected before allocating or reading their contents. diff --git a/python/damacy/__init__.py b/python/damacy/__init__.py index 8f448500..9e0d9ee6 100644 --- a/python/damacy/__init__.py +++ b/python/damacy/__init__.py @@ -94,7 +94,6 @@ "Sampler", "ShutdownError", "SpatialQuery", - "SpatialResolver", "Stats", "Status", "StorageError", @@ -1944,5 +1943,4 @@ def stats_reset(self) -> None: ResolvedSpatialQuery, Sampler, SpatialQuery, - SpatialResolver, ) diff --git a/python/damacy/_native.pyi b/python/damacy/_native.pyi index 3897ec18..7eba5228 100644 --- a/python/damacy/_native.pyi +++ b/python/damacy/_native.pyi @@ -265,14 +265,10 @@ def ngff_info(image: object, /) -> dict[str, Any]: ... def spatial_resolve( image: object, shape: tuple[int, ...], - samples: int, - dtype: int, transform: tuple[tuple[float, ...], ...], filter: int, boundary: int, constant_value: float, level: int, /, -) -> object: ... -def spatial_info(resolution: object, /) -> dict[str, Any]: ... -def spatial_sample(resolution: object, /) -> dict[str, Any]: ... +) -> dict[str, Any]: ... diff --git a/python/damacy/_spatial.c b/python/damacy/_spatial.c index 0e2f410c..8c07fb65 100644 --- a/python/damacy/_spatial.c +++ b/python/damacy/_spatial.c @@ -5,7 +5,6 @@ #include "damacy_spatial.h" static const char image_name[] = "damacy.NgffImage"; -static const char resolution_name[] = "damacy.SpatialResolution"; static void image_destroy(PyObject* capsule) @@ -13,13 +12,6 @@ image_destroy(PyObject* capsule) damacy_ngff_image_destroy(PyCapsule_GetPointer(capsule, image_name)); } -static void -resolution_destroy(PyObject* capsule) -{ - damacy_spatial_resolution_destroy( - PyCapsule_GetPointer(capsule, resolution_name)); -} - static int put(PyObject* dict, const char* key, PyObject* value) { @@ -226,20 +218,41 @@ parse_transform(PyObject* object, uint8_t rank, struct damacy_affine* transform) return 0; } +static PyObject* +resolution_value(const struct damacy_spatial_resolution* info) +{ + PyObject* result = PyDict_New(); + if (!result) + return NULL; + if (put(result, "uri", PyUnicode_FromString(info->uri)) || + put(result, "level", PyLong_FromUnsignedLong(info->level)) || + put(result, "shape", integers(info->output_shape, info->rank)) || + put(result, "source_shape", integers(info->source_shape, info->rank)) || + put(result, + "output_to_source", + affine(&info->output_to_source, info->rank)) || + put(result, "source_bounds_index", bounds(&info->source_bounds_index)) || + put(result, "read_bounds_index", bounds(&info->read_bounds_index)) || + put(result, + "requires_resampling", + PyBool_FromLong(info->requires_resampling))) { + Py_DECREF(result); + return NULL; + } + return result; +} + static PyObject* resolve_query(PyObject* self, PyObject* args) { (void)self; PyObject *capsule, *shape_object, *transform; - unsigned int samples, dtype; int filter, boundary, level; double value; if (!PyArg_ParseTuple(args, - "OOIIOiidi", + "OOOiidi", &capsule, &shape_object, - &samples, - &dtype, &transform, &filter, &boundary, @@ -249,8 +262,7 @@ resolve_query(PyObject* self, PyObject* args) struct damacy_ngff_image* image = PyCapsule_GetPointer(capsule, image_name); if (!image) return NULL; - struct damacy_batch_spec output = { .samples_per_batch = samples, - .dtype = (enum damacy_dtype)dtype }; + int64_t output_shape[DAMACY_MAX_RANK]; PyObject* shape = PySequence_Fast(shape_object, "shape must be a sequence"); if (!shape) return NULL; @@ -260,10 +272,8 @@ resolve_query(PyObject* self, PyObject* args) PyErr_SetString(PyExc_ValueError, "invalid output rank"); return NULL; } - output.sample_rank = (uint8_t)rank; for (Py_ssize_t i = 0; i < rank; ++i) { - output.sample_shape[i] = - PyLong_AsLongLong(PySequence_Fast_GET_ITEM(shape, i)); + output_shape[i] = PyLong_AsLongLong(PySequence_Fast_GET_ITEM(shape, i)); if (PyErr_Occurred()) { Py_DECREF(shape); return NULL; @@ -276,84 +286,16 @@ resolve_query(PyObject* self, PyObject* args) .constant_value = value }, .level = level }; - if (parse_transform( - transform, output.sample_rank, &query.output_to_reference)) + if (parse_transform(transform, (uint8_t)rank, &query.output_to_reference)) return NULL; - struct damacy_spatial_resolution* resolution = NULL; + struct damacy_spatial_resolution resolution; enum damacy_status status; - Py_BEGIN_ALLOW_THREADS status = - damacy_spatial_resolve(image, &output, &query, &resolution); + Py_BEGIN_ALLOW_THREADS status = damacy_spatial_resolve( + image, &query, (uint8_t)rank, output_shape, &resolution); Py_END_ALLOW_THREADS if (status != DAMACY_OK) return api_raise_status( status, "resolve spatial query"); - PyObject* result = - PyCapsule_New(resolution, resolution_name, resolution_destroy); - if (!result) - damacy_spatial_resolution_destroy(resolution); - return result; -} - -static PyObject* -resolution_info(PyObject* self, PyObject* capsule) -{ - (void)self; - struct damacy_spatial_resolution* resolution = - PyCapsule_GetPointer(capsule, resolution_name); - if (!resolution) - return NULL; - const struct damacy_spatial_info* info = - damacy_spatial_resolution_info(resolution); - PyObject* result = PyDict_New(); - if (!result) - return NULL; - if (put(result, "uri", PyUnicode_FromString(info->uri)) || - put(result, "level", PyLong_FromUnsignedLong(info->level)) || - put(result, "shape", integers(info->output_shape, info->rank)) || - put(result, "source_shape", integers(info->source_shape, info->rank)) || - put(result, - "output_to_source", - affine(&info->output_to_source, info->rank)) || - put(result, "source_bounds_index", bounds(&info->source_bounds_index)) || - put(result, "read_bounds_index", bounds(&info->read_bounds_index)) || - put(result, - "requires_resampling", - PyBool_FromLong(info->requires_resampling))) { - Py_DECREF(result); - return NULL; - } - return result; -} - -static PyObject* -resolution_sample(PyObject* self, PyObject* capsule) -{ - (void)self; - struct damacy_spatial_resolution* resolution = - PyCapsule_GetPointer(capsule, resolution_name); - if (!resolution) - return NULL; - struct damacy_sample sample; - enum damacy_status status = - damacy_spatial_resolution_sample(resolution, &sample); - if (status != DAMACY_OK) - return api_raise_status(status, - "submit spatial query: resampling required"); - PyObject* axes = PyList_New(sample.rank); - if (!axes) - return NULL; - for (uint8_t i = 0; i < sample.rank; ++i) { - PyObject* axis = Py_BuildValue("(s(LL))", - "interval", - (long long)sample.axes[i].interval.beg, - (long long)sample.axes[i].interval.end); - if (!axis) { - Py_DECREF(axes); - return NULL; - } - PyList_SET_ITEM(axes, i, axis); - } - PyObject* result = - Py_BuildValue("{s:s,s:O}", "uri", sample.uri, "axes", axes); - Py_DECREF(axes); + PyObject* result = resolution_value(&resolution); + damacy_spatial_resolution_clear(&resolution); return result; } @@ -361,8 +303,6 @@ static PyMethodDef methods[] = { { "ngff_load", load_image, METH_VARARGS, NULL }, { "ngff_info", image_info, METH_O, NULL }, { "spatial_resolve", resolve_query, METH_VARARGS, NULL }, - { "spatial_info", resolution_info, METH_O, NULL }, - { "spatial_sample", resolution_sample, METH_O, NULL }, { NULL, NULL, 0, NULL } }; diff --git a/python/damacy/_spatial.py b/python/damacy/_spatial.py index 3b44ded9..8ef25ebf 100644 --- a/python/damacy/_spatial.py +++ b/python/damacy/_spatial.py @@ -7,7 +7,14 @@ from dataclasses import dataclass from typing import Any, Literal -from . import BatchSpec, FileMetadataReader, Sample, _component, _native, _positive_int +from . import ( + FileMetadataReader, + Status, + UnsupportedOperation, + _component, + _native, + _positive_int, +) def _index(value: int, name: str) -> int: @@ -108,6 +115,29 @@ def __init__( ) object.__setattr__(self, "data_type", info["data_type"]) + def resolve( + self, query: SpatialQuery, *, shape: Iterable[int] + ) -> ResolvedSpatialQuery: + """Resolve source geometry without I/O, batching, or output allocation.""" + if not isinstance(query, SpatialQuery): + raise TypeError("query must be a SpatialQuery") + output_shape = tuple( + _positive_int(value, "shape extent", (1 << 52) - 1) for value in shape + ) + if len(output_shape) != len(self.axes): + raise ValueError("output rank must match the NGFF image") + info = _component( + _native.spatial_resolve, + self._native, + output_shape, + query.output_to_reference, + {"nearest": 1, "linear": 2}[query.sampler.filter], + {"error": 1, "constant": 2, "clamp": 3}[query.sampler.boundary], + query.sampler.constant_value, + -1 if query.level == "auto" else query.level, + ) + return ResolvedSpatialQuery._from_native(info, query.sampler) + @dataclass(frozen=True, slots=True) class Sampler: @@ -181,11 +211,10 @@ def __init__( @dataclass(frozen=True, slots=True, init=False) class ResolvedSpatialQuery: - """Owned source selection and geometry, independent of the loaded image. + """Owned source geometry returned by NgffImage.resolve(). - Aligned results can be pushed to a Pipeline or converted with ``as_sample``. - Other results raise UnsupportedOperation when converted or pushed; their - geometry remains available for inspection without invoking a decoder. + Aligned results can be pushed to a Pipeline. Other results raise + UnsupportedOperation when pushed; their geometry remains inspectable. """ uri: str @@ -197,59 +226,27 @@ class ResolvedSpatialQuery: source_bounds_index: tuple[tuple[int, int], ...] read_bounds_index: tuple[tuple[int, int], ...] requires_resampling: bool - _native: object + + def __init__(self) -> None: + raise TypeError("use NgffImage.resolve() to create a ResolvedSpatialQuery") @classmethod - def _from_native(cls, native: object, sampler: Sampler) -> ResolvedSpatialQuery: + def _from_native( + cls, info: dict[str, Any], sampler: Sampler + ) -> ResolvedSpatialQuery: result = object.__new__(cls) - object.__setattr__(result, "_native", native) object.__setattr__(result, "sampler", sampler) - for name, value in _native.spatial_info(native).items(): + for name, value in info.items(): object.__setattr__(result, name, value) return result def _to_native(self) -> dict[str, Any]: - return _component(_native.spatial_sample, self._native) - - def as_sample(self) -> Sample: - """Return an interval sample, or fail if resampling is required.""" - sample = self._to_native() - return Sample(uri=sample["uri"], aabb=tuple(axis[1] for axis in sample["axes"])) - - -@dataclass(frozen=True, slots=True) -class SpatialResolver: - """Resolve queries using injected image metadata and fixed output geometry. - - Resolution does no I/O and can be shared across threads. The same - BatchSpec should be supplied to the downstream Pipeline. - """ - - image: NgffImage - output: BatchSpec - - def __post_init__(self) -> None: - if not isinstance(self.image, NgffImage) or not isinstance( - self.output, BatchSpec - ): - raise TypeError("SpatialResolver requires an NgffImage and a BatchSpec") - if len(self.image.axes) != len(self.output.shape): - raise ValueError("output rank must match the NGFF image") - - def resolve(self, query: SpatialQuery) -> ResolvedSpatialQuery: - """Choose the source level, map coordinates, and bound source reads.""" - if not isinstance(query, SpatialQuery): - raise TypeError("query must be a SpatialQuery") - native = _component( - _native.spatial_resolve, - self.image._native, - self.output.shape, - self.output.samples, - int(self.output.dtype), - query.output_to_reference, - {"nearest": 1, "linear": 2}[query.sampler.filter], - {"error": 1, "constant": 2, "clamp": 3}[query.sampler.boundary], - query.sampler.constant_value, - -1 if query.level == "auto" else query.level, - ) - return ResolvedSpatialQuery._from_native(native, query.sampler) + if self.requires_resampling: + error = UnsupportedOperation("submit spatial query: resampling required") + error.status = Status.UNSUPPORTED + error.what = "submit spatial query: resampling required" + raise error + return { + "uri": self.uri, + "axes": [("interval", span) for span in self.source_bounds_index], + } diff --git a/python/tests/test_spatial.py b/python/tests/test_spatial.py index c334e07d..5fdfd350 100644 --- a/python/tests/test_spatial.py +++ b/python/tests/test_spatial.py @@ -4,6 +4,7 @@ import dataclasses import gc import json +import pickle from concurrent.futures import ThreadPoolExecutor from typing import Literal @@ -194,17 +195,17 @@ def test_resolved_crops_decode_at_selected_level(tmp_path, spatial_executor, lev reader = damacy.FileMetadataReader(concurrency=2) image = damacy.NgffImage(root, reader=reader, multiscale_index=0) output = damacy.BatchSpec(2, (4, 4)) - resolver = damacy.SpatialResolver(image=image, output=output) + shape = output.shape factor = 2**level - resolved = resolver.resolve(query(factor, (factor * 2, factor * 3))) + resolved = image.resolve(query(factor, (factor * 2, factor * 3)), shape=shape) assert resolved.level == level assert not resolved.requires_resampling - assert resolved.as_sample().aabb == ((2, 6), (3, 7)) + assert resolved.source_bounds_index == ((2, 6), (3, 7)) assert resolved.output_to_source == ((1, 0, 2), (0, 1, 3)) - del resolver, image + del image gc.collect() with pipeline(spatial_executor, output, reader) as p: - p.push([resolved, resolved.as_sample()]) + p.push([resolved, pickle.loads(pickle.dumps(resolved))]) with p.pop() as batch: expected = arrays[level][2:6, 3:7] np.testing.assert_array_equal( @@ -218,20 +219,20 @@ def test_resampling_fails_before_decoding_and_pipeline_recovers( root = tmp_path / "image" arrays = write_image(root, data=True) output = damacy.BatchSpec(1, (4, 4)) - resolver = damacy.SpatialResolver(load(root), output) - rotated = resolver.resolve( + image = load(root) + shape = output.shape + rotated = image.resolve( damacy.SpatialQuery( output_to_reference=((1, -1, 12), (1, 1, 4)), sampler=damacy.Sampler(filter="linear"), - ) + ), + shape=shape, ) assert rotated.requires_resampling - with pytest.raises(damacy.UnsupportedOperation, match="resampling required"): - rotated.as_sample() with pipeline(spatial_executor, output) as p: with pytest.raises(damacy.UnsupportedOperation, match="resampling required"): p.push([rotated]) - p.push([resolver.resolve(query())]) + p.push([image.resolve(query(), shape=shape)]) with p.pop() as batch: np.testing.assert_array_equal(read_batch(batch)[0], arrays[0][:4, :4]) @@ -239,9 +240,7 @@ def test_resampling_fails_before_decoding_and_pipeline_recovers( def test_fixed_shape_is_checked_at_push(tmp_path, spatial_executor): root = tmp_path / "image" write_image(root) - resolved = damacy.SpatialResolver(load(root), damacy.BatchSpec(1, (3, 3))).resolve( - query() - ) + resolved = load(root).resolve(query(), shape=(3, 3)) with ( pipeline(spatial_executor, damacy.BatchSpec(1, (4, 4))) as p, pytest.raises(damacy.InvalidArgument), @@ -282,10 +281,10 @@ def test_ngff_center_offsets_stay_in_adapter(tmp_path): assert image.levels[0].origin_reference_index == (0, 0) assert image.levels[1].origin_reference_index == (0, 0) assert image.levels[2].origin_reference_index == (-1.5, -1.5) - resolver = damacy.SpatialResolver(image, damacy.BatchSpec(1, (4, 4))) - resolved = resolver.resolve(query(2)) + shape = (4, 4) + resolved = image.resolve(query(2), shape=shape) assert resolved.level == 1 and not resolved.requires_resampling - assert resolver.resolve(query(4)).requires_resampling + assert image.resolve(query(4), shape=shape).requires_resampling def test_anisotropic_rotated_level_selection_and_time_channel_axes(tmp_path): @@ -303,16 +302,17 @@ def test_anisotropic_rotated_level_selection_and_time_channel_axes(tmp_path): scales=[[1, 1, 1, 1, 1], [1, 1, 1, 2, 4]], ) image = load(root) - resolver = damacy.SpatialResolver(image, damacy.BatchSpec(1, (1, 1, 2, 2, 2))) + shape = (1, 1, 2, 2, 2) transform = np.zeros((5, 6)) transform[:5, :5] = np.eye(5) transform[:2, 5] = (1, 2) transform[2:5, 2:5] = np.array([[0, -1, 0], [2, 0, 0], [0, 0, 4]]) transform[2:, 5] = (4, 4, 0) - resolved = resolver.resolve( + resolved = image.resolve( damacy.SpatialQuery( output_to_reference=transform, sampler=damacy.Sampler(filter="nearest") - ) + ), + shape=shape, ) assert resolved.level == 1 assert resolved.output_to_source[2] == (0, 0, 0, -1, 0, 4) @@ -320,29 +320,36 @@ def test_anisotropic_rotated_level_selection_and_time_channel_axes(tmp_path): assert resolved.source_bounds_index[:2] == ((1, 2), (2, 3)) transform[2, 3] = -0.5 assert ( - resolver.resolve( - damacy.SpatialQuery(output_to_reference=transform, sampler=damacy.Sampler()) + image.resolve( + damacy.SpatialQuery( + output_to_reference=transform, sampler=damacy.Sampler() + ), + shape=shape, ).level == 0 ) transform[0, 0] = 2 with pytest.raises(damacy.InvalidArgument): - resolver.resolve( - damacy.SpatialQuery(output_to_reference=transform, sampler=damacy.Sampler()) + image.resolve( + damacy.SpatialQuery( + output_to_reference=transform, sampler=damacy.Sampler() + ), + shape=shape, ) def test_selection_uses_smallest_spacing_not_column_lengths(tmp_path): root = tmp_path / "image" write_image(root) - resolver = damacy.SpatialResolver(load(root), damacy.BatchSpec(1, (2, 2))) + image = load(root) + shape = (2, 2) q = damacy.SpatialQuery( output_to_reference=((3, 2, 0), (2, 3, 0)), sampler=damacy.Sampler(filter="nearest"), ) - assert resolver.resolve(q).level == 0 + assert image.resolve(q, shape=shape).level == 0 forced = dataclasses.replace(q, level=1) - assert resolver.resolve(forced).level == 1 + assert image.resolve(forced, shape=shape).level == 1 @pytest.mark.parametrize( @@ -360,11 +367,13 @@ def test_selection_uses_smallest_spacing_not_column_lengths(tmp_path): def test_source_and_read_bounds(tmp_path, filter, boundary, offset, source, reads): root = tmp_path / "image" write_image(root) - resolver = damacy.SpatialResolver(load(root), damacy.BatchSpec(1, (4, 4))) - resolved = resolver.resolve( + image = load(root) + shape = (4, 4) + resolved = image.resolve( query( offset=(offset, 0), sampler=damacy.Sampler(filter=filter, boundary=boundary) - ) + ), + shape=shape, ) assert resolved.source_bounds_index[0] == source assert resolved.read_bounds_index[0] == reads @@ -384,10 +393,14 @@ def test_source_and_read_bounds(tmp_path, filter, boundary, offset, source, read def test_invalid_geometry(tmp_path, transform): root = tmp_path / "image" write_image(root) - resolver = damacy.SpatialResolver(load(root), damacy.BatchSpec(1, (4, 4))) + image = load(root) + shape = (4, 4) with pytest.raises(damacy.InvalidArgument): - resolver.resolve( - damacy.SpatialQuery(output_to_reference=transform, sampler=damacy.Sampler()) + image.resolve( + damacy.SpatialQuery( + output_to_reference=transform, sampler=damacy.Sampler() + ), + shape=shape, ) @@ -506,11 +519,10 @@ def test_multiscale_selection_and_escaped_strings(tmp_path): lambda text: text.replace( '"scale": [1, 1]', '"scale": [1, 1], "scale": [2, 2]' ), - lambda text: text + " garbage", lambda text: text[:-1] + ",}", ], ) -def test_malformed_json_is_rejected(tmp_path, mutate): +def test_invalid_consumed_metadata_is_rejected(tmp_path, mutate): root = tmp_path / "image" write_image(root) path = root / "zarr.json" @@ -522,12 +534,17 @@ def test_malformed_json_is_rejected(tmp_path, mutate): def test_resolution_is_independent_of_files_and_can_run_in_threads(tmp_path): root = tmp_path / "image" write_image(root) - resolver = damacy.SpatialResolver(load(root), damacy.BatchSpec(1, (4, 4))) + image = load(root) + shape = (4, 4) for metadata in root.rglob("zarr.json"): metadata.unlink() with ThreadPoolExecutor(max_workers=4) as workers: - results = list(workers.map(resolver.resolve, [query(2)] * 16)) - assert all(r.level == 1 and r.as_sample().aabb == ((0, 4), (0, 4)) for r in results) + results = list( + workers.map(lambda q: image.resolve(q, shape=shape), [query(2)] * 16) + ) + assert all( + r.level == 1 and r.source_bounds_index == ((0, 4), (0, 4)) for r in results + ) def test_query_values_are_copied_and_validated(): @@ -564,10 +581,12 @@ def test_anisotropic_crop_decodes_with_time_and_channel_axes( data=True, ) output = damacy.BatchSpec(1, (1, 2, 4, 4, 4)) - resolver = damacy.SpatialResolver(load(root), output) + image = load(root) + shape = output.shape transform = np.column_stack((np.diag([1, 1, 1, 2, 4]), [1, 1, 4, 4, 0])) - resolved = resolver.resolve( - damacy.SpatialQuery(output_to_reference=transform, sampler=damacy.Sampler()) + resolved = image.resolve( + damacy.SpatialQuery(output_to_reference=transform, sampler=damacy.Sampler()), + shape=shape, ) assert resolved.level == 1 and not resolved.requires_resampling with pipeline(spatial_executor, output) as p: @@ -582,7 +601,8 @@ def test_automatic_level_selection_matches_singular_values(tmp_path): root = tmp_path / "image" scales = np.array([[1, 1, 1], [1, 2, 4], [2, 4, 8]]) write_image(root, shape=(32, 32, 32), scales=scales.tolist()) - resolver = damacy.SpatialResolver(load(root), damacy.BatchSpec(1, (2, 2, 2))) + image = load(root) + shape = (2, 2, 2) rng = np.random.default_rng(532) for _ in range(64): left, _ = np.linalg.qr(rng.normal(size=(3, 3))) @@ -592,11 +612,12 @@ def test_automatic_level_selection_matches_singular_values(tmp_path): for level, scale in enumerate(scales): if np.linalg.svd(matrix / scale[:, None], compute_uv=False)[-1] >= 1: expected = level - resolved = resolver.resolve( + resolved = image.resolve( damacy.SpatialQuery( output_to_reference=np.column_stack((matrix, [8, 8, 8])), sampler=damacy.Sampler(boundary="constant"), - ) + ), + shape=shape, ) assert resolved.level == expected @@ -628,7 +649,7 @@ def test_native_query_validation(tmp_path): root = tmp_path / "image" write_image(root) image = load(root) - args = (image._native, (4, 4), 1, 0) + args = (image._native, (4, 4)) with pytest.raises(ValueError): _native.spatial_resolve(*args, ((1, 0), (0, 1)), 1, 1, 0, -1) for filter, boundary, value, level in [ @@ -662,7 +683,8 @@ def test_duplicate_array_fields_are_rejected(tmp_path, key, value): def test_collapsed_3d_transforms_are_rejected(tmp_path): root = tmp_path / "image" write_image(root, shape=(16, 16, 16)) - resolver = damacy.SpatialResolver(load(root), damacy.BatchSpec(1, (2, 2, 2))) + image = load(root) + shape = (2, 2, 2) rng = np.random.default_rng(832) for _ in range(128): matrix = rng.normal(size=(3, 3)) @@ -672,23 +694,61 @@ def test_collapsed_3d_transforms_are_rejected(tmp_path): sampler=damacy.Sampler(boundary="constant"), ) with pytest.raises(damacy.InvalidArgument): - resolver.resolve(q) + image.resolve(q, shape=shape) @pytest.mark.parametrize("spacing,level", [(0.5, 0), (3.0, 1)]) def test_level_selection_with_large_scale_difference(tmp_path, spacing, level): root = tmp_path / "image" write_image(root) - resolver = damacy.SpatialResolver(load(root), damacy.BatchSpec(1, (2, 2))) + image = load(root) + shape = (2, 2) large = 2**28 matrix = ( (large + spacing / 2, large - spacing / 2, 0), (large - spacing / 2, large + spacing / 2, 0), ) - resolved = resolver.resolve( + resolved = image.resolve( damacy.SpatialQuery( output_to_reference=matrix, sampler=damacy.Sampler(boundary="constant"), - ) + ), + shape=shape, ) assert resolved.level == level + + +@pytest.mark.parametrize("tail", [",", ", garbage", ', {"unused": [garbage,]}']) +def test_unselected_metadata_is_not_validated(tmp_path, tail): + root = tmp_path / "image" + write_image(root) + path = root / "zarr.json" + metadata = json.loads(path.read_text()) + entry = metadata["attributes"]["ome"]["multiscales"][0] + entry["coordinateTransformations"] = {"unused": ["not", "interpreted"]} + metadata["attributes"]["unused"] = {"arbitrary": {"nested": None}} + text = json.dumps(metadata) + position = text.rindex("]") + path.write_text(text[:position] + tail + text[position:] + " ignored suffix") + image = load(root) + resolved = image.resolve(query(2), shape=(4, 4)) + assert resolved.level == 1 + assert not resolved.requires_resampling + + +def test_resolution_only_needs_output_shape(tmp_path): + root = tmp_path / "image" + write_image(root) + image = load(root) + shape = (1 << 40, 1 << 40) + resolved = image.resolve( + query(sampler=damacy.Sampler(boundary="constant")), shape=shape + ) + assert resolved.shape == shape + assert resolved.read_bounds_index == ((0, 32), (0, 32)) + assert resolved.requires_resampling + for invalid in [(), (4,), (0, 4), (-1, 4), (1 << 52, 4), (1.5, 4)]: + with pytest.raises((ValueError, TypeError)): + image.resolve(query(), shape=invalid) + with pytest.raises(TypeError, match=r"NgffImage\.resolve"): + damacy.ResolvedSpatialQuery() diff --git a/src/damacy_spatial.h b/src/damacy_spatial.h index 252a6fde..2182b775 100644 --- a/src/damacy_spatial.h +++ b/src/damacy_spatial.h @@ -8,7 +8,6 @@ extern "C" #endif struct damacy_ngff_image; - struct damacy_spatial_resolution; enum damacy_ngff_axis_kind { @@ -85,7 +84,7 @@ extern "C" int32_t level; }; - struct damacy_spatial_info + struct damacy_spatial_resolution { const char* uri; uint32_t level; @@ -111,15 +110,14 @@ extern "C" enum damacy_status damacy_spatial_resolve( const struct damacy_ngff_image* image, - const struct damacy_batch_spec* output, const struct damacy_spatial_query* query, - struct damacy_spatial_resolution** out); - const struct damacy_spatial_info* damacy_spatial_resolution_info( - const struct damacy_spatial_resolution* resolution); + uint8_t rank, + const int64_t* output_shape, + struct damacy_spatial_resolution* out); enum damacy_status damacy_spatial_resolution_sample( const struct damacy_spatial_resolution* resolution, struct damacy_sample* out); - void damacy_spatial_resolution_destroy( + void damacy_spatial_resolution_clear( struct damacy_spatial_resolution* resolution); #ifdef __cplusplus diff --git a/src/ngff/metadata.c b/src/ngff/metadata.c index 04381936..4a4385a0 100644 --- a/src/ngff/metadata.c +++ b/src/ngff/metadata.c @@ -143,14 +143,9 @@ string_equal(struct json_node node, const char* expected) } static int -root_document(struct cslice src, struct json_node* root) +root_object(struct cslice src, struct json_node* root) { - if (json_resolve(src, NULL, 0, root, NULL) || root->type != JSON_OBJECT) - return 1; - for (const char* p = root->s.end; p < src.end; ++p) - if (*p != ' ' && *p != '\t' && *p != '\n' && *p != '\r') - return 1; - return 0; + return json_resolve(src, NULL, 0, root, NULL) || root->type != JSON_OBJECT; } static enum damacy_status @@ -355,7 +350,7 @@ ngff_parse_group(struct cslice src, *out = NULL; struct json_node root, node, ome, multiscale; uint64_t format; - if (root_document(src, &root) || member(root, "zarr_format", &node) || + if (root_object(src, &root) || member(root, "zarr_format", &node) || json_as_uint(node, &format) || format != 3 || member(root, "node_type", &node) || !string_equal(node, "group") || member(root, "attributes", &node) || member(node, "ome", &ome) || @@ -378,17 +373,7 @@ ngff_parse_group(struct cslice src, status = read_axes(node, &image->info); if (status != DAMACY_OK) goto Fail; - double group_scale[DAMACY_MAX_RANK] = { 0 }; - double group_origin[DAMACY_MAX_RANK] = { 0 }; - enum json_err error = member(multiscale, "coordinateTransformations", &node); - if (error == JSON_OK) { - status = read_transform(node, image->info.rank, group_scale, group_origin); - if (status != DAMACY_OK) - goto Fail; - } else if (error != JSON_ERR_NOT_FOUND) { - status = DAMACY_INVAL; - goto Fail; - } + enum json_err error; status = DAMACY_INVAL; if (member(multiscale, "datasets", &node)) goto Fail; @@ -444,7 +429,7 @@ ngff_parse_array(struct cslice src, uint32_t index) { struct json_node root, names, value; - if (root_document(src, &root)) + if (root_object(src, &root)) return DAMACY_INVAL; static const char* keys[] = { "zarr_format", "node_type", "shape", "data_type", diff --git a/src/query/spatial.c b/src/query/spatial.c index a110586c..e8e128bc 100644 --- a/src/query/spatial.c +++ b/src/query/spatial.c @@ -1,20 +1,12 @@ #include "damacy_spatial.h" -#include "damacy_config.h" #include "ngff/ngff.h" -#include "pipeline/components.h" #include #include #include #include -struct damacy_spatial_resolution -{ - struct damacy_spatial_info info; - struct damacy_sample sample; -}; - static uint8_t spatial_matrix(const struct damacy_ngff_image* image, const struct damacy_affine* transform, @@ -130,20 +122,14 @@ minimum_spacing(const struct damacy_ngff_image* image, static enum damacy_status validate_query(const struct damacy_ngff_image* image, - const struct damacy_batch_spec* output, - const struct damacy_spatial_query* query) + const struct damacy_spatial_query* query, + uint8_t rank, + const int64_t* output_shape) { - if (!image || !output || !query) + if (!image || !output_shape || !query) return DAMACY_INVAL; - if (output->sample_rank != image->info.rank) + if (rank != image->info.rank) return DAMACY_RANK; - int64_t shape[DAMACY_MAX_RANK + 1], strides[DAMACY_MAX_RANK + 1]; - uint64_t bytes; - enum damacy_status status = batch_spec_layout(output, shape, strides, &bytes); - if (status != DAMACY_OK) - return status; - if (!cast_path_supported(output->dtype, image->dtype)) - return DAMACY_DTYPE; if (query->level < DAMACY_LEVEL_AUTO || (query->level >= 0 && (uint32_t)query->level >= image->info.level_count)) return DAMACY_INVAL; @@ -160,8 +146,8 @@ validate_query(const struct damacy_ngff_image* image, for (uint8_t i = 0; i < image->info.rank; ++i) { double offset = query->output_to_reference.offset[i]; int space = image->info.axes[i].kind == DAMACY_NGFF_SPACE; - if (!isfinite(offset) || - output->sample_shape[i] > INT64_C(4503599627370495)) + if (!isfinite(offset) || output_shape[i] <= 0 || + output_shape[i] > INT64_C(4503599627370495)) return DAMACY_INVAL; if (!space && offset != floor(offset)) return DAMACY_INVAL; @@ -187,7 +173,7 @@ clamp_index(int64_t value, int64_t end) static enum damacy_status resolve_bounds(const struct damacy_ngff_image* image, - struct damacy_spatial_info* info) + struct damacy_spatial_resolution* info) { info->source_bounds_index.rank = info->read_bounds_index.rank = info->rank; for (uint8_t i = 0; i < info->rank; ++i) { @@ -238,14 +224,15 @@ resolve_bounds(const struct damacy_ngff_image* image, enum damacy_status damacy_spatial_resolve(const struct damacy_ngff_image* image, - const struct damacy_batch_spec* output, const struct damacy_spatial_query* query, - struct damacy_spatial_resolution** out) + uint8_t rank, + const int64_t* output_shape, + struct damacy_spatial_resolution* out) { if (!out) return DAMACY_INVAL; - *out = NULL; - enum damacy_status status = validate_query(image, output, query); + *out = (struct damacy_spatial_resolution){ 0 }; + enum damacy_status status = validate_query(image, query, rank, output_shape); if (status != DAMACY_OK) return status; uint32_t level_index = query->level < 0 ? 0 : (uint32_t)query->level; @@ -255,56 +242,32 @@ damacy_spatial_resolve(const struct damacy_ngff_image* image, 1 - 64 * DBL_EPSILON) level_index = i; const struct damacy_ngff_level* level = &image->levels[level_index]; - struct damacy_spatial_resolution* resolution = calloc(1, sizeof(*resolution)); - if (!resolution) - return DAMACY_OOM; - struct damacy_spatial_info* info = &resolution->info; - info->level = level_index; - info->rank = image->info.rank; - info->sampler = query->sampler; - for (uint8_t i = 0; i < info->rank; ++i) { - info->output_shape[i] = output->sample_shape[i]; - info->source_shape[i] = level->shape[i]; + struct damacy_spatial_resolution result = { .level = level_index, + .rank = rank, + .sampler = query->sampler }; + for (uint8_t i = 0; i < rank; ++i) { + result.output_shape[i] = output_shape[i]; + result.source_shape[i] = level->shape[i]; double scale = level->scale_to_reference[i]; - info->output_to_source.offset[i] = (query->output_to_reference.offset[i] - - level->origin_reference_index[i]) / - scale; - for (uint8_t j = 0; j < info->rank; ++j) - info->output_to_source.linear[i][j] = + result.output_to_source.offset[i] = (query->output_to_reference.offset[i] - + level->origin_reference_index[i]) / + scale; + for (uint8_t j = 0; j < rank; ++j) + result.output_to_source.linear[i][j] = query->output_to_reference.linear[i][j] / scale; } - status = resolve_bounds(image, info); - if (status != DAMACY_OK) { - free(resolution); + status = resolve_bounds(image, &result); + if (status != DAMACY_OK) return status; - } char* uri = malloc(strlen(level->uri) + 1); - if (!uri) { - free(resolution); + if (!uri) return DAMACY_OOM; - } strcpy(uri, level->uri); - info->uri = uri; - if (!info->requires_resampling) { - resolution->sample.uri = uri; - resolution->sample.rank = info->rank; - for (uint8_t i = 0; i < info->rank; ++i) - resolution->sample.axes[i] = - (struct damacy_axis_selection){ .kind = DAMACY_AXIS_INTERVAL, - .interval = - info->source_bounds_index.dims[i] }; - } - *out = resolution; + result.uri = uri; + *out = result; return DAMACY_OK; } -const struct damacy_spatial_info* -damacy_spatial_resolution_info( - const struct damacy_spatial_resolution* resolution) -{ - return resolution ? &resolution->info : NULL; -} - enum damacy_status damacy_spatial_resolution_sample( const struct damacy_spatial_resolution* resolution, @@ -313,19 +276,26 @@ damacy_spatial_resolution_sample( if (!out) return DAMACY_INVAL; *out = (struct damacy_sample){ 0 }; - if (!resolution) + if (!resolution || !resolution->uri || !resolution->rank || + resolution->rank > DAMACY_MAX_RANK) return DAMACY_INVAL; - if (resolution->info.requires_resampling) + if (resolution->requires_resampling) return DAMACY_UNSUPPORTED; - *out = resolution->sample; + out->uri = resolution->uri; + out->rank = resolution->rank; + for (uint8_t i = 0; i < resolution->rank; ++i) + out->axes[i] = (struct damacy_axis_selection){ + .kind = DAMACY_AXIS_INTERVAL, + .interval = resolution->source_bounds_index.dims[i] + }; return DAMACY_OK; } void -damacy_spatial_resolution_destroy(struct damacy_spatial_resolution* resolution) +damacy_spatial_resolution_clear(struct damacy_spatial_resolution* resolution) { if (resolution) { - free((void*)resolution->info.uri); - free(resolution); + free((void*)resolution->uri); + *resolution = (struct damacy_spatial_resolution){ 0 }; } } diff --git a/tests/test_spatial.c b/tests/test_spatial.c index d0e67975..0ad3fa69 100644 --- a/tests/test_spatial.c +++ b/tests/test_spatial.c @@ -74,10 +74,7 @@ identity_query(void) }; } -static struct damacy_batch_spec output = { .dtype = DAMACY_F32, - .sample_shape = { 4, 4 }, - .sample_rank = 2, - .samples_per_batch = 1 }; +static const int64_t output_shape[] = { 4, 4 }; static int test_corner_conversion_and_copy(void) @@ -97,32 +94,30 @@ test_corner_conversion_and_copy(void) EXPECT(metadata->levels[2].origin_reference_index[d] == -1.5); } struct damacy_spatial_query query = identity_query(); - struct damacy_spatial_resolution* resolved; - EXPECT(damacy_spatial_resolve(image, &output, &query, &resolved) == + struct damacy_spatial_resolution resolved; + EXPECT(damacy_spatial_resolve(image, &query, 2, output_shape, &resolved) == DAMACY_OK); - EXPECT(damacy_spatial_resolution_info(resolved)->level == 0); + EXPECT(resolved.level == 0); struct damacy_sample sample; - EXPECT(damacy_spatial_resolution_sample(resolved, &sample) == DAMACY_OK); + EXPECT(damacy_spatial_resolution_sample(&resolved, &sample) == DAMACY_OK); EXPECT(sample.axes[0].interval.beg == 0 && sample.axes[0].interval.end == 4); - damacy_spatial_resolution_destroy(resolved); + damacy_spatial_resolution_clear(&resolved); query.output_to_reference.linear[0][0] = 2; query.output_to_reference.linear[1][1] = 2; query.output_to_reference.offset[0] = 4; query.output_to_reference.offset[1] = 6; query.sampler.filter = DAMACY_FILTER_LINEAR; - EXPECT(damacy_spatial_resolve(image, &output, &query, &resolved) == + EXPECT(damacy_spatial_resolve(image, &query, 2, output_shape, &resolved) == DAMACY_OK); - const struct damacy_spatial_info* info = - damacy_spatial_resolution_info(resolved); - EXPECT(info->level == 1 && !info->requires_resampling); + EXPECT(resolved.level == 1 && !resolved.requires_resampling); damacy_ngff_image_destroy(image); memset(&query, 0, sizeof(query)); - EXPECT(damacy_spatial_resolution_sample(resolved, &sample) == DAMACY_OK); + EXPECT(damacy_spatial_resolution_sample(&resolved, &sample) == DAMACY_OK); EXPECT(!strcmp(sample.uri, "volume/1")); EXPECT(sample.rank == 2 && sample.axes[0].kind == DAMACY_AXIS_INTERVAL); EXPECT(sample.axes[0].interval.beg == 2 && sample.axes[0].interval.end == 6); EXPECT(sample.axes[1].interval.beg == 3 && sample.axes[1].interval.end == 7); - damacy_spatial_resolution_destroy(resolved); + damacy_spatial_resolution_clear(&resolved); return 0; } @@ -139,42 +134,39 @@ test_scale_rotation_and_shear(void) query.output_to_reference.linear[0][1] = -s; query.output_to_reference.linear[1][0] = s; query.output_to_reference.linear[1][1] = s; - struct damacy_spatial_resolution* resolved; - EXPECT(damacy_spatial_resolve(image, &output, &query, &resolved) == + struct damacy_spatial_resolution resolved; + EXPECT(damacy_spatial_resolve(image, &query, 2, output_shape, &resolved) == DAMACY_OK); - const struct damacy_spatial_info* info = - damacy_spatial_resolution_info(resolved); - EXPECT(info->level == 1 && info->requires_resampling); - EXPECT(fabs(info->output_to_source.linear[0][0] - s / 2) < 1e-15); + EXPECT(resolved.level == 1 && resolved.requires_resampling); + EXPECT(fabs(resolved.output_to_source.linear[0][0] - s / 2) < 1e-15); struct damacy_sample sample; - EXPECT(damacy_spatial_resolution_sample(resolved, &sample) == + EXPECT(damacy_spatial_resolution_sample(&resolved, &sample) == DAMACY_UNSUPPORTED); EXPECT(!sample.uri && !sample.rank); - damacy_spatial_resolution_destroy(resolved); + damacy_spatial_resolution_clear(&resolved); query.output_to_reference.linear[0][0] = 3; query.output_to_reference.linear[0][1] = 2; query.output_to_reference.linear[1][0] = 2; query.output_to_reference.linear[1][1] = 3; - EXPECT(damacy_spatial_resolve(image, &output, &query, &resolved) == + EXPECT(damacy_spatial_resolve(image, &query, 2, output_shape, &resolved) == DAMACY_OK); - EXPECT(damacy_spatial_resolution_info(resolved)->level == 0); - damacy_spatial_resolution_destroy(resolved); + EXPECT(resolved.level == 0); + damacy_spatial_resolution_clear(&resolved); query = identity_query(); query.output_to_reference.linear[0][0] = 4; query.output_to_reference.linear[1][1] = 4; - EXPECT(damacy_spatial_resolve(image, &output, &query, &resolved) == + EXPECT(damacy_spatial_resolve(image, &query, 2, output_shape, &resolved) == DAMACY_OK); - info = damacy_spatial_resolution_info(resolved); - EXPECT(info->level == 2 && info->requires_resampling); - EXPECT(info->output_to_source.offset[0] == 0.375); - damacy_spatial_resolution_destroy(resolved); + EXPECT(resolved.level == 2 && resolved.requires_resampling); + EXPECT(resolved.output_to_source.offset[0] == 0.375); + damacy_spatial_resolution_clear(&resolved); query = identity_query(); query.output_to_reference.linear[0][0] = 0.5; query.output_to_reference.linear[1][1] = 0.5; - EXPECT(damacy_spatial_resolve(image, &output, &query, &resolved) == + EXPECT(damacy_spatial_resolve(image, &query, 2, output_shape, &resolved) == DAMACY_OK); - EXPECT(damacy_spatial_resolution_info(resolved)->level == 0); - damacy_spatial_resolution_destroy(resolved); + EXPECT(resolved.level == 0); + damacy_spatial_resolution_clear(&resolved); damacy_ngff_image_destroy(image); return 0; } @@ -186,43 +178,38 @@ test_sampler_bounds(void) EXPECT(!image_create(&image)); struct damacy_spatial_query query = identity_query(); query.output_to_reference.offset[0] = -0.25; - struct damacy_spatial_resolution* resolved; - EXPECT(damacy_spatial_resolve(image, &output, &query, &resolved) == + struct damacy_spatial_resolution resolved; + EXPECT(damacy_spatial_resolve(image, &query, 2, output_shape, &resolved) == DAMACY_OK); - const struct damacy_spatial_info* info = - damacy_spatial_resolution_info(resolved); - EXPECT(info->source_bounds_index.dims[0].beg == 0); - EXPECT(info->source_bounds_index.dims[0].end == 4); - damacy_spatial_resolution_destroy(resolved); + EXPECT(resolved.source_bounds_index.dims[0].beg == 0); + EXPECT(resolved.source_bounds_index.dims[0].end == 4); + damacy_spatial_resolution_clear(&resolved); query.sampler.filter = DAMACY_FILTER_LINEAR; - EXPECT(damacy_spatial_resolve(image, &output, &query, &resolved) == + EXPECT(damacy_spatial_resolve(image, &query, 2, output_shape, &resolved) == DAMACY_INVAL); - EXPECT(!resolved); + EXPECT(!resolved.uri && !resolved.rank); query.sampler.boundary = DAMACY_BOUNDARY_CONSTANT; query.sampler.constant_value = 17; - EXPECT(damacy_spatial_resolve(image, &output, &query, &resolved) == + EXPECT(damacy_spatial_resolve(image, &query, 2, output_shape, &resolved) == DAMACY_OK); - info = damacy_spatial_resolution_info(resolved); - EXPECT(info->source_bounds_index.dims[0].beg == -1); - EXPECT(info->source_bounds_index.dims[0].end == 4); - EXPECT(info->read_bounds_index.dims[0].beg == 0); - EXPECT(info->read_bounds_index.dims[0].end == 4); - damacy_spatial_resolution_destroy(resolved); + EXPECT(resolved.source_bounds_index.dims[0].beg == -1); + EXPECT(resolved.source_bounds_index.dims[0].end == 4); + EXPECT(resolved.read_bounds_index.dims[0].beg == 0); + EXPECT(resolved.read_bounds_index.dims[0].end == 4); + damacy_spatial_resolution_clear(&resolved); query.output_to_reference.offset[0] = -20; - EXPECT(damacy_spatial_resolve(image, &output, &query, &resolved) == + EXPECT(damacy_spatial_resolve(image, &query, 2, output_shape, &resolved) == DAMACY_OK); - info = damacy_spatial_resolution_info(resolved); - EXPECT(info->read_bounds_index.dims[0].beg == 0); - EXPECT(info->read_bounds_index.dims[0].end == 0); - damacy_spatial_resolution_destroy(resolved); + EXPECT(resolved.read_bounds_index.dims[0].beg == 0); + EXPECT(resolved.read_bounds_index.dims[0].end == 0); + damacy_spatial_resolution_clear(&resolved); query.sampler.boundary = DAMACY_BOUNDARY_CLAMP; query.sampler.constant_value = 0; - EXPECT(damacy_spatial_resolve(image, &output, &query, &resolved) == + EXPECT(damacy_spatial_resolve(image, &query, 2, output_shape, &resolved) == DAMACY_OK); - info = damacy_spatial_resolution_info(resolved); - EXPECT(info->read_bounds_index.dims[0].beg == 0); - EXPECT(info->read_bounds_index.dims[0].end == 1); - damacy_spatial_resolution_destroy(resolved); + EXPECT(resolved.read_bounds_index.dims[0].beg == 0); + EXPECT(resolved.read_bounds_index.dims[0].end == 1); + damacy_spatial_resolution_clear(&resolved); damacy_ngff_image_destroy(image); return 0; } @@ -233,38 +220,69 @@ test_invalid_queries(void) struct damacy_ngff_image* image; EXPECT(!image_create(&image)); struct damacy_spatial_query query = identity_query(); - struct damacy_spatial_resolution* resolved; + struct damacy_spatial_resolution resolved; double invalid[] = { 0, NAN, INFINITY, 1e100 }; for (size_t i = 0; i < sizeof(invalid) / sizeof(*invalid); ++i) { query.output_to_reference.linear[0][0] = invalid[i]; - EXPECT(damacy_spatial_resolve(image, &output, &query, &resolved) == + EXPECT(damacy_spatial_resolve(image, &query, 2, output_shape, &resolved) == DAMACY_INVAL); - EXPECT(!resolved); + EXPECT(!resolved.uri && !resolved.rank); } query = identity_query(); query.level = -2; - EXPECT(damacy_spatial_resolve(image, &output, &query, &resolved) == + EXPECT(damacy_spatial_resolve(image, &query, 2, output_shape, &resolved) == DAMACY_INVAL); query.level = 3; - EXPECT(damacy_spatial_resolve(image, &output, &query, &resolved) == + EXPECT(damacy_spatial_resolve(image, &query, 2, output_shape, &resolved) == DAMACY_INVAL); query.level = 0; query.sampler.constant_value = 1; - EXPECT(damacy_spatial_resolve(image, &output, &query, &resolved) == + EXPECT(damacy_spatial_resolve(image, &query, 2, output_shape, &resolved) == DAMACY_INVAL); query.sampler.constant_value = 0; query.sampler.filter = 0; - EXPECT(damacy_spatial_resolve(image, &output, &query, &resolved) == + EXPECT(damacy_spatial_resolve(image, &query, 2, output_shape, &resolved) == DAMACY_INVAL); query = identity_query(); - struct damacy_batch_spec wrong = output; - wrong.sample_rank = 3; - EXPECT(damacy_spatial_resolve(image, &wrong, &query, &resolved) == + EXPECT(damacy_spatial_resolve(image, &query, 3, output_shape, &resolved) == DAMACY_RANK); damacy_ngff_image_destroy(image); return 0; } +static int +test_resolution_shape_and_clear(void) +{ + struct damacy_ngff_image* image; + EXPECT(!image_create(&image)); + struct damacy_spatial_query query = identity_query(); + query.sampler.boundary = DAMACY_BOUNDARY_CONSTANT; + const int64_t shape[] = { INT64_C(1) << 40, INT64_C(1) << 40 }; + struct damacy_spatial_resolution resolved; + EXPECT(damacy_spatial_resolve(image, &query, 2, shape, &resolved) == + DAMACY_OK); + EXPECT(resolved.output_shape[0] == shape[0]); + EXPECT(resolved.requires_resampling); + damacy_spatial_resolution_clear(&resolved); + damacy_spatial_resolution_clear(&resolved); + EXPECT(!resolved.uri && !resolved.rank); + struct damacy_sample sample; + EXPECT(damacy_spatial_resolution_sample(&resolved, &sample) == DAMACY_INVAL); + EXPECT(!sample.uri && !sample.rank); + const int64_t invalid[] = { 0, -1, INT64_C(1) << 52 }; + for (size_t i = 0; i < sizeof(invalid) / sizeof(*invalid); ++i) { + int64_t bad_shape[] = { invalid[i], 4 }; + EXPECT(damacy_spatial_resolve(image, &query, 2, bad_shape, &resolved) == + DAMACY_INVAL); + EXPECT(!resolved.uri && !resolved.rank); + } + EXPECT(damacy_spatial_resolve(image, &query, 2, NULL, &resolved) == + DAMACY_INVAL); + EXPECT(damacy_spatial_resolve(image, &query, 2, shape, NULL) == DAMACY_INVAL); + damacy_ngff_image_destroy(image); + return 0; +} + static int test_metadata_limits_and_load(void) { @@ -340,10 +358,7 @@ test_collapsed_volume(void) .levels = &level, .dtype = dtype_u16 }; - struct damacy_batch_spec shape = { .dtype = DAMACY_F32, - .sample_rank = 3, - .sample_shape = { 2, 2, 2 }, - .samples_per_batch = 1 }; + const int64_t shape[] = { 2, 2, 2 }; struct damacy_spatial_query query = { .output_to_reference = { .linear = { { 0.7854589457317591, 0.4315455905112043, @@ -359,17 +374,17 @@ test_collapsed_volume(void) .boundary = DAMACY_BOUNDARY_CONSTANT }, .level = DAMACY_LEVEL_AUTO }; - struct damacy_spatial_resolution* resolved = NULL; - EXPECT(damacy_spatial_resolve(&image, &shape, &query, &resolved) == + struct damacy_spatial_resolution resolved; + EXPECT(damacy_spatial_resolve(&image, &query, 3, shape, &resolved) == DAMACY_INVAL); - EXPECT(!resolved); + EXPECT(!resolved.uri && !resolved.rank); query.output_to_reference = (struct damacy_affine){ .linear = { { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } } }; - EXPECT(damacy_spatial_resolve(&image, &shape, &query, &resolved) == + EXPECT(damacy_spatial_resolve(&image, &query, 3, shape, &resolved) == DAMACY_OK); - EXPECT(!damacy_spatial_resolution_info(resolved)->requires_resampling); - damacy_spatial_resolution_destroy(resolved); + EXPECT(!resolved.requires_resampling); + damacy_spatial_resolution_clear(&resolved); return 0; } @@ -381,6 +396,7 @@ main(void) RUN(test_sampler_bounds); RUN(test_invalid_queries); RUN(test_collapsed_volume); + RUN(test_resolution_shape_and_clear); RUN(test_metadata_limits_and_load); return 0; } From b3249041ec5f58f2185e687bd0c7538b6dc78cca Mon Sep 17 00:00:00 2001 From: Nathan Clack Date: Thu, 24 Sep 2026 21:21:53 +0000 Subject: [PATCH 04/10] ngff: round level values at load --- docs/spatial.md | 18 +++++-- python/tests/test_spatial.py | 71 +++++++++++++++++++++++++++- src/ngff/metadata.c | 13 ++++++ tests/test_spatial.c | 91 ++++++++++++++++++++++++++++++++++++ 4 files changed, 188 insertions(+), 5 deletions(-) diff --git a/docs/spatial.md b/docs/spatial.md index 0d365932..d678688f 100644 --- a/docs/spatial.md +++ b/docs/spatial.md @@ -157,6 +157,13 @@ identity. For a 2× level with a half-reference-voxel center translation, Damacy preserves that declared alignment. The API never asks callers to add or remove NGFF center offsets themselves. +Writers usually compute scales and translations in floating point, so these +values can miss by a rounding error. For example, scales `0.1` and `0.3` give +`R_l = 2.9999999999999996`. Loading replaces `R_l` or `origin_l` with the +nearest multiple of 1/256 when it is that close: within `64 * DBL_EPSILON` +times `R_l` for the ratio, or times `(|T_l| + |T_0|) / S_0 + R_l` for the +origin. Whole-number ratios and whole- or half-voxel origins then come out exact. + `NgffLevel.scale_to_reference` and `.origin_reference_index` expose this adapted map. Given query matrix `A` and offset `b`, resolution computes: @@ -183,10 +190,13 @@ map and `origin_l + R_l * beg` as its offset. The current copy path requires the resolved linear map to equal identity, its offsets to be integral, and its source bounds to be inside the array. -Fractional values are not rounded into a crop. Small floating-point differences -in metadata can therefore make a result require resampling. `requires_resampling` -reports this, and `Pipeline.push()` rejects such results until a resampler is -available. The pipeline checks the resulting shape against its own `BatchSpec`. +The resolver checks this exactly and never rounds a query into a crop. Because +loading removes rounding error from the level values, a crop built from them +as above is exact. A query computed another way, such as from physical +coordinates, can still be off by a rounding error and then require resampling. +`requires_resampling` reports this, and `Pipeline.push()` rejects such results +until a resampler is available. The pipeline checks the resulting shape against +its own `BatchSpec`. ## Sampler and bounds diff --git a/python/tests/test_spatial.py b/python/tests/test_spatial.py index 5fdfd350..d3f5a8bc 100644 --- a/python/tests/test_spatial.py +++ b/python/tests/test_spatial.py @@ -108,6 +108,26 @@ def edit(path, keys, value): path.write_text(json.dumps(metadata)) +def set_transforms(root, levels): + for i, (scale, translation) in enumerate(levels): + edit( + root / "zarr.json", + ( + "attributes", + "ome", + "multiscales", + 0, + "datasets", + i, + "coordinateTransformations", + ), + [ + {"type": "scale", "scale": [scale, scale]}, + {"type": "translation", "translation": [translation, translation]}, + ], + ) + + def load(path, **kwargs): return damacy.NgffImage( path, @@ -118,7 +138,11 @@ def load(path, **kwargs): def query( - scale=1, offset=(0, 0), *, sampler=None, level: int | Literal["auto"] = "auto" + scale: float = 1, + offset=(0, 0), + *, + sampler=None, + level: int | Literal["auto"] = "auto", ): rank = len(offset) return damacy.SpatialQuery( @@ -287,6 +311,51 @@ def test_ngff_center_offsets_stay_in_adapter(tmp_path): assert image.resolve(query(4), shape=shape).requires_resampling +def test_float_translations_keep_crops_aligned(tmp_path, spatial_executor): + root = tmp_path / "image" + arrays = write_image(root, data=True) + scales = [0.325 * 2**level for level in range(3)] + set_transforms(root, [(s, 12.3 + (s - scales[0]) / 2) for s in scales]) + image = load(root) + assert image.levels[1].scale_to_reference == (2, 2) + assert image.levels[1].origin_reference_index == (0, 0) + output = damacy.BatchSpec(1, (4, 4)) + resolved = image.resolve(query(2), shape=output.shape) + assert resolved.level == 1 and not resolved.requires_resampling + with pipeline(spatial_executor, output) as p: + p.push([resolved]) + with p.pop() as batch: + np.testing.assert_array_equal(read_batch(batch)[0], arrays[1][:4, :4]) + + +def test_float_scale_ratio_keeps_crops_aligned(tmp_path): + root = tmp_path / "image" + write_image(root, scales=[[1, 1], [3, 3]]) + set_transforms(root, [(0.1, 0), (0.3, 0)]) + image = load(root) + assert image.levels[1].scale_to_reference == (3, 3) + assert image.levels[1].origin_reference_index == (-1, -1) + for beg in range(3): + resolved = image.resolve(query(3, (3 * beg - 1,) * 2), shape=(2, 2)) + assert resolved.level == 1 and not resolved.requires_resampling + assert resolved.source_bounds_index == ((beg, beg + 2),) * 2 + + +def test_documented_level_crop_is_exact(tmp_path): + root = tmp_path / "image" + write_image(root, scales=[[1, 1], [3, 3]]) + set_transforms(root, [(0.1, 12.3), (0.1 * 3, 12.3 + (0.1 * 3 - 0.1) / 2)]) + image = load(root) + scale = image.levels[1].scale_to_reference[0] + origin = image.levels[1].origin_reference_index[0] + for beg in range(1, 4): + resolved = image.resolve( + query(scale, (origin + scale * beg,) * 2, level=1), shape=(2, 2) + ) + assert not resolved.requires_resampling + assert resolved.source_bounds_index == ((beg, beg + 2),) * 2 + + def test_anisotropic_rotated_level_selection_and_time_channel_axes(tmp_path): root = tmp_path / "image" write_image( diff --git a/src/ngff/metadata.c b/src/ngff/metadata.c index 4a4385a0..4a7494b0 100644 --- a/src/ngff/metadata.c +++ b/src/ngff/metadata.c @@ -2,6 +2,7 @@ #include "zarr/zarr_metadata.h" +#include #include #include #include @@ -309,6 +310,14 @@ read_level(struct json_node node, value, rank, level->scale_to_reference, level->origin_reference_index); } +// Removes writers' floating-point error so crops built from levels stay exact. +static double +round_if_close(double value, double error) +{ + double nearest = round(value * 256) / 256; + return fabs(value - nearest) <= error ? nearest : value; +} + static enum damacy_status normalize_levels(struct damacy_ngff_image* image) { @@ -329,7 +338,11 @@ normalize_levels(struct damacy_ngff_image* image) if (i && scale < image->levels[i - 1].scale_to_reference[d]) return DAMACY_INVAL; double ratio = scale / base_scale; + ratio = round_if_close(ratio, 64 * DBL_EPSILON * ratio); double offset = (origin - base_origin) / base_scale + 0.5 * (1 - ratio); + double magnitude = + (fabs(origin) + fabs(base_origin)) / base_scale + ratio; + offset = round_if_close(offset, 64 * DBL_EPSILON * magnitude); if (!isfinite(ratio) || !isfinite(offset)) return DAMACY_INVAL; level->scale_to_reference[d] = ratio; diff --git a/tests/test_spatial.c b/tests/test_spatial.c index 0ad3fa69..ef7b77bf 100644 --- a/tests/test_spatial.c +++ b/tests/test_spatial.c @@ -63,6 +63,45 @@ image_create(struct damacy_ngff_image** image) return 0; } +static int +pyramid_create(const char* scale0, + const char* translation0, + const char* scale1, + const char* translation1, + struct damacy_ngff_image** image) +{ + char text[2048]; + snprintf( + text, + sizeof(text), + "{\"zarr_format\":3,\"node_type\":\"group\",\"attributes\":{\"ome\":{" + "\"version\":\"0.5\",\"multiscales\":[{\"axes\":[" + "{\"name\":\"y\",\"type\":\"space\"}," + "{\"name\":\"x\",\"type\":\"space\"}],\"datasets\":[" + "{\"path\":\"0\",\"coordinateTransformations\":[" + "{\"type\":\"scale\",\"scale\":[%s,%s]}," + "{\"type\":\"translation\",\"translation\":[%s,%s]}]}," + "{\"path\":\"1\",\"coordinateTransformations\":[" + "{\"type\":\"scale\",\"scale\":[%s,%s]}," + "{\"type\":\"translation\",\"translation\":[%s,%s]}]}]}]}}}", + scale0, + scale0, + translation0, + translation0, + scale1, + scale1, + translation1, + translation1); + EXPECT(ngff_parse_group(text_slice(text), "volume", 0, 2, image) == + DAMACY_OK); + for (uint32_t i = 0; i < 2; ++i) { + char array[1024]; + array_json(array, sizeof(array), 64 >> i); + EXPECT(ngff_parse_array(text_slice(array), *image, i) == DAMACY_OK); + } + return 0; +} + static struct damacy_spatial_query identity_query(void) { @@ -121,6 +160,57 @@ test_corner_conversion_and_copy(void) return 0; } +static int +level_one_crop(const struct damacy_ngff_image* image, + double scale, + double offset, + int64_t beg) +{ + struct damacy_spatial_query query = identity_query(); + for (uint8_t d = 0; d < 2; ++d) { + query.output_to_reference.linear[d][d] = scale; + query.output_to_reference.offset[d] = offset; + } + struct damacy_spatial_resolution resolved; + EXPECT(damacy_spatial_resolve(image, &query, 2, output_shape, &resolved) == + DAMACY_OK); + int aligned = resolved.level == 1 && !resolved.requires_resampling && + resolved.source_bounds_index.dims[0].beg == beg && + resolved.source_bounds_index.dims[1].beg == beg; + damacy_spatial_resolution_clear(&resolved); + EXPECT(aligned); + return 0; +} + +static int +test_level_rounding_error_is_removed(void) +{ + struct damacy_ngff_image* image; + EXPECT(!pyramid_create("0.325", "12.3", "0.65", "12.4625", &image)); + const struct damacy_ngff_level* level = + &damacy_ngff_image_info(image)->levels[1]; + EXPECT(level->scale_to_reference[0] == 2); + EXPECT(level->origin_reference_index[0] == 0); + EXPECT(!level_one_crop(image, 2, 0, 0)); + damacy_ngff_image_destroy(image); + EXPECT(!pyramid_create("0.1", "0", "0.3", "0", &image)); + level = &damacy_ngff_image_info(image)->levels[1]; + EXPECT(level->scale_to_reference[0] == 3); + EXPECT(level->origin_reference_index[0] == -1); + for (int64_t beg = 0; beg < 3; ++beg) + EXPECT(!level_one_crop(image, 3, 3.0 * (double)beg - 1, beg)); + damacy_ngff_image_destroy(image); + EXPECT(!pyramid_create("0.1", "12.3", "0.30000000000000004", "12.4", &image)); + level = &damacy_ngff_image_info(image)->levels[1]; + for (int64_t beg = 1; beg < 4; ++beg) { + double scale = level->scale_to_reference[0]; + double offset = level->origin_reference_index[0] + scale * (double)beg; + EXPECT(!level_one_crop(image, scale, offset, beg)); + } + damacy_ngff_image_destroy(image); + return 0; +} + static int test_scale_rotation_and_shear(void) { @@ -392,6 +482,7 @@ int main(void) { RUN(test_corner_conversion_and_copy); + RUN(test_level_rounding_error_is_removed); RUN(test_scale_rotation_and_shear); RUN(test_sampler_bounds); RUN(test_invalid_queries); From b3003958f837a900237285a14347df6d150b596a Mon Sep 17 00:00:00 2001 From: Nathan Clack Date: Thu, 24 Sep 2026 21:22:25 +0000 Subject: [PATCH 05/10] ngff: stop claiming the metadata reader --- docs/spatial.md | 6 ++---- python/damacy/_spatial.py | 2 +- python/tests/test_spatial.py | 10 +++------- src/ngff/load.c | 4 ---- 4 files changed, 6 insertions(+), 16 deletions(-) diff --git a/docs/spatial.md b/docs/spatial.md index d678688f..dc0cab0e 100644 --- a/docs/spatial.md +++ b/docs/spatial.md @@ -78,10 +78,8 @@ and geometry, so it remains usable after the image is released. Python results contain ordinary immutable values and can be pickled for transfer between processes. Source metadata and data must remain unchanged while the image is in use. -Loading is synchronous and uses the supplied reader exclusively for the -duration of the call. A reader already serving a pipeline or another load is -rejected. The example reuses its reader after loading finishes. Use a separate -reader if loading more images while a pipeline is running. +Loading is synchronous. It uses only the reader's settings, so one reader can +serve a running pipeline and other loads at the same time. ## Coordinates and the output grid diff --git a/python/damacy/_spatial.py b/python/damacy/_spatial.py index 8ef25ebf..81d4b6e3 100644 --- a/python/damacy/_spatial.py +++ b/python/damacy/_spatial.py @@ -71,7 +71,7 @@ class NgffImage: """Load an immutable OME-Zarr 0.5 image description through the given reader. Loading reads the image group's metadata and each level's array metadata. - It finishes before returning and releases the reader for other uses. + It finishes before returning. ``multiscale_index`` explicitly selects an entry in ``ome.multiscales``. """ diff --git a/python/tests/test_spatial.py b/python/tests/test_spatial.py index d3f5a8bc..4717e87e 100644 --- a/python/tests/test_spatial.py +++ b/python/tests/test_spatial.py @@ -691,16 +691,12 @@ def test_automatic_level_selection_matches_singular_values(tmp_path): assert resolved.level == expected -def test_active_metadata_reader_is_not_reused(tmp_path, spatial_executor): +def test_metadata_reader_is_shared_with_running_pipeline(tmp_path, spatial_executor): root = tmp_path / "image" write_image(root) reader = damacy.FileMetadataReader(concurrency=2) - with ( - pipeline(spatial_executor, damacy.BatchSpec(1, (4, 4)), reader), - pytest.raises(damacy.InvalidArgument), - ): - damacy.NgffImage(root, reader=reader, multiscale_index=0) - assert damacy.NgffImage(root, reader=reader, multiscale_index=0).levels + with pipeline(spatial_executor, damacy.BatchSpec(1, (4, 4)), reader): + assert damacy.NgffImage(root, reader=reader, multiscale_index=0).levels @pytest.mark.parametrize("path", ["zarr.json", "1/zarr.json"]) diff --git a/src/ngff/load.c b/src/ngff/load.c index adac9699..447088c1 100644 --- a/src/ngff/load.c +++ b/src/ngff/load.c @@ -79,9 +79,6 @@ damacy_ngff_image_load(struct damacy_metadata_reader* reader, limits->max_levels > INT32_MAX || !limits->max_metadata_bytes || limits->max_metadata_bytes > SIZE_MAX) return DAMACY_INVAL; - int expected = 0; - if (!atomic_compare_exchange_strong(&reader->active, &expected, 1)) - return DAMACY_INVAL; struct metadata_store_async* store = metadata_store_async_create( (int)reader->concurrency, NULL, &reader->latency); struct metadata_read read = { .mutex = platform_mutex_new(), @@ -124,7 +121,6 @@ damacy_ngff_image_load(struct damacy_metadata_reader* reader, metadata_store_async_destroy(store); platform_mutex_free(read.mutex); platform_cond_free(read.cond); - atomic_store(&reader->active, 0); damacy_ngff_image_destroy(image); free(read.data); return status; From dc509deca1fa0aa5da891ebd584eeddf4118bab6 Mon Sep 17 00:00:00 2001 From: Nathan Clack Date: Thu, 24 Sep 2026 21:22:31 +0000 Subject: [PATCH 06/10] docs: raise spatial example CPU budget --- docs/spatial.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/spatial.md b/docs/spatial.md index dc0cab0e..aad6dde8 100644 --- a/docs/spatial.md +++ b/docs/spatial.md @@ -57,7 +57,7 @@ planner = damacy.ChunkPlanner( ) executor = damacy.CpuExecutor( reader=damacy.FileReader(workers=4), - limits=damacy.CpuLimits(max_memory_bytes=256 << 20), + limits=damacy.CpuLimits(max_memory_bytes=3 << 30), ) with damacy.Pipeline( planner=planner, From f2861cc798019434b46696e4a2778cdbf2795551 Mon Sep 17 00:00:00 2001 From: Nathan Clack Date: Thu, 24 Sep 2026 21:23:14 +0000 Subject: [PATCH 07/10] ngff: skip keys that do not decode --- python/tests/test_spatial.py | 11 +++++++++++ src/ngff/metadata.c | 8 +++----- tests/test_spatial.c | 14 ++++++++++++++ 3 files changed, 28 insertions(+), 5 deletions(-) diff --git a/python/tests/test_spatial.py b/python/tests/test_spatial.py index 4717e87e..3e4772b1 100644 --- a/python/tests/test_spatial.py +++ b/python/tests/test_spatial.py @@ -579,6 +579,17 @@ def test_multiscale_selection_and_escaped_strings(tmp_path): load(root) +def test_unrelated_keys_are_skipped(tmp_path): + root = tmp_path / "image" + write_image(root) + path = root / "zarr.json" + metadata = json.loads(path.read_text()) + metadata["attributes"]["ome"]["\0"] = 1 + metadata["attributes"]["ome"]["multiscales"][0]["axes"][0]["\ud800"] = 2 + path.write_text(json.dumps(metadata)) + assert load(root).levels[0].shape == (32, 32) + + @pytest.mark.parametrize( "mutate", [ diff --git a/src/ngff/metadata.c b/src/ngff/metadata.c index 4a7494b0..8fae9706 100644 --- a/src/ngff/metadata.c +++ b/src/ngff/metadata.c @@ -114,13 +114,11 @@ member(struct json_node node, const char* key, struct json_node* out) struct json_node name, value; int found = 0; while ((error = json_object_iter_next(&it, &name, &value)) == JSON_OK) { - if (!cslice_len(name.s)) - continue; const char* text = NULL; enum damacy_status status = read_string(name, &text); - if (status != DAMACY_OK) - return status == DAMACY_OOM ? JSON_ERR_OOM : JSON_ERR_PARSE; - int matches = !strcmp(key, text); + if (status == DAMACY_OOM) + return JSON_ERR_OOM; + int matches = status == DAMACY_OK && !strcmp(key, text); free((void*)text); if (matches) { if (found) diff --git a/tests/test_spatial.c b/tests/test_spatial.c index ef7b77bf..4fc94646 100644 --- a/tests/test_spatial.c +++ b/tests/test_spatial.c @@ -373,6 +373,19 @@ test_resolution_shape_and_clear(void) return 0; } +static int +test_unrelated_keys_are_skipped(void) +{ + char text[4096]; + snprintf( + text, sizeof(text), "{\"\\u0000\":1,\"\\ud800\":2,\"\":3,%s", group + 1); + struct damacy_ngff_image* image = NULL; + EXPECT(ngff_parse_group(text_slice(text), "volume", 0, 3, &image) == + DAMACY_OK); + damacy_ngff_image_destroy(image); + return 0; +} + static int test_metadata_limits_and_load(void) { @@ -488,6 +501,7 @@ main(void) RUN(test_invalid_queries); RUN(test_collapsed_volume); RUN(test_resolution_shape_and_clear); + RUN(test_unrelated_keys_are_skipped); RUN(test_metadata_limits_and_load); return 0; } From 8cdb7ecee4883a97c7558838cda1761f6f1b766a Mon Sep 17 00:00:00 2001 From: Nathan Clack Date: Thu, 24 Sep 2026 21:24:20 +0000 Subject: [PATCH 08/10] python: hide native image handle --- python/damacy/_spatial.py | 21 ++++++++++++++++++--- python/tests/test_spatial.py | 13 +++++++++++++ 2 files changed, 31 insertions(+), 3 deletions(-) diff --git a/python/damacy/_spatial.py b/python/damacy/_spatial.py index 81d4b6e3..12f47bb3 100644 --- a/python/damacy/_spatial.py +++ b/python/damacy/_spatial.py @@ -5,7 +5,7 @@ import os from collections.abc import Iterable from dataclasses import dataclass -from typing import Any, Literal +from typing import Any, Literal, NoReturn from . import ( FileMetadataReader, @@ -66,8 +66,15 @@ class NgffLevel: origin_reference_index: tuple[float, ...] +class _NativeImage: + # Keeps the native handle out of the dataclass fields, so repr, ==, hash, + # and asdict() see only the metadata. + __slots__ = ("_native",) + _native: object + + @dataclass(frozen=True, slots=True, init=False) -class NgffImage: +class NgffImage(_NativeImage): """Load an immutable OME-Zarr 0.5 image description through the given reader. Loading reads the image group's metadata and each level's array metadata. @@ -78,7 +85,6 @@ class NgffImage: axes: tuple[NgffAxis, ...] levels: tuple[NgffLevel, ...] data_type: str - _native: object def __init__( self, @@ -115,6 +121,15 @@ def __init__( ) object.__setattr__(self, "data_type", info["data_type"]) + def __copy__(self) -> NgffImage: + return self + + def __deepcopy__(self, memo: dict[int, Any]) -> NgffImage: + return self + + def __reduce__(self) -> NoReturn: + raise TypeError("NgffImage cannot be pickled; load it in each process") + def resolve( self, query: SpatialQuery, *, shape: Iterable[int] ) -> ResolvedSpatialQuery: diff --git a/python/tests/test_spatial.py b/python/tests/test_spatial.py index 3e4772b1..5f4b6cba 100644 --- a/python/tests/test_spatial.py +++ b/python/tests/test_spatial.py @@ -1,5 +1,6 @@ from __future__ import annotations +import copy import ctypes import dataclasses import gc @@ -627,6 +628,18 @@ def test_resolution_is_independent_of_files_and_can_run_in_threads(tmp_path): ) +def test_image_compares_and_copies_by_metadata(tmp_path): + root = tmp_path / "image" + write_image(root) + image = load(root) + assert "_native" not in repr(image) + assert image == load(root) and hash(image) == hash(load(root)) + assert json.loads(json.dumps(dataclasses.asdict(image)))["data_type"] == "uint16" + assert copy.deepcopy(image).resolve(query(), shape=(4, 4)).level == 0 + with pytest.raises(TypeError, match="pickled"): + pickle.dumps(image) + + def test_query_values_are_copied_and_validated(): matrix = [[1, 0, 2], [0, 1, 3]] q = damacy.SpatialQuery(output_to_reference=matrix, sampler=damacy.Sampler()) From b4a286a9a65816595b958c9c7b940f84870ea3c6 Mon Sep 17 00:00:00 2001 From: Nathan Clack Date: Thu, 24 Sep 2026 21:25:35 +0000 Subject: [PATCH 09/10] python: take spatial enums from native --- python/damacy/_native.pyi | 9 +++++++++ python/damacy/_spatial.c | 15 +++++++++++++++ python/damacy/_spatial.py | 13 ++++++++++--- 3 files changed, 34 insertions(+), 3 deletions(-) diff --git a/python/damacy/_native.pyi b/python/damacy/_native.pyi index 7eba5228..1a08f07f 100644 --- a/python/damacy/_native.pyi +++ b/python/damacy/_native.pyi @@ -52,6 +52,15 @@ GDS_AUTO: Final[int] GDS_ON: Final[int] GDS_OFF: Final[int] +# ---- spatial sampler and level integers (mirror damacy_spatial.h) -------- + +FILTER_NEAREST: Final[int] +FILTER_LINEAR: Final[int] +BOUNDARY_ERROR: Final[int] +BOUNDARY_CONSTANT: Final[int] +BOUNDARY_CLAMP: Final[int] +LEVEL_AUTO: Final[int] + # ---- exceptions --------------------------------------------------------- class DamacyError(RuntimeError): diff --git a/python/damacy/_spatial.c b/python/damacy/_spatial.c index 8c07fb65..bc9b6dc5 100644 --- a/python/damacy/_spatial.c +++ b/python/damacy/_spatial.c @@ -309,5 +309,20 @@ static PyMethodDef methods[] = { int spatial_register(PyObject* module) { + static const struct + { + const char* name; + int value; + } constants[] = { + { "FILTER_NEAREST", DAMACY_FILTER_NEAREST }, + { "FILTER_LINEAR", DAMACY_FILTER_LINEAR }, + { "BOUNDARY_ERROR", DAMACY_BOUNDARY_ERROR }, + { "BOUNDARY_CONSTANT", DAMACY_BOUNDARY_CONSTANT }, + { "BOUNDARY_CLAMP", DAMACY_BOUNDARY_CLAMP }, + { "LEVEL_AUTO", DAMACY_LEVEL_AUTO }, + }; + for (size_t i = 0; i < sizeof(constants) / sizeof(*constants); ++i) + if (PyModule_AddIntConstant(module, constants[i].name, constants[i].value)) + return -1; return PyModule_AddFunctions(module, methods); } diff --git a/python/damacy/_spatial.py b/python/damacy/_spatial.py index 12f47bb3..d2336ac7 100644 --- a/python/damacy/_spatial.py +++ b/python/damacy/_spatial.py @@ -16,6 +16,13 @@ _positive_int, ) +_FILTERS = {"nearest": _native.FILTER_NEAREST, "linear": _native.FILTER_LINEAR} +_BOUNDARIES = { + "error": _native.BOUNDARY_ERROR, + "constant": _native.BOUNDARY_CONSTANT, + "clamp": _native.BOUNDARY_CLAMP, +} + def _index(value: int, name: str) -> int: if isinstance(value, bool): @@ -146,10 +153,10 @@ def resolve( self._native, output_shape, query.output_to_reference, - {"nearest": 1, "linear": 2}[query.sampler.filter], - {"error": 1, "constant": 2, "clamp": 3}[query.sampler.boundary], + _FILTERS[query.sampler.filter], + _BOUNDARIES[query.sampler.boundary], query.sampler.constant_value, - -1 if query.level == "auto" else query.level, + _native.LEVEL_AUTO if query.level == "auto" else query.level, ) return ResolvedSpatialQuery._from_native(info, query.sampler) From 504a0ff0c27630f475dabcab471523b5476daa8f Mon Sep 17 00:00:00 2001 From: Nathan Clack Date: Thu, 24 Sep 2026 22:14:31 +0000 Subject: [PATCH 10/10] store: bounded metadata read on macOS --- src/store/metadata_store_async.posix.c | 19 +++++++++++++++++-- 1 file changed, 17 insertions(+), 2 deletions(-) diff --git a/src/store/metadata_store_async.posix.c b/src/store/metadata_store_async.posix.c index 8839869c..a1540571 100644 --- a/src/store/metadata_store_async.posix.c +++ b/src/store/metadata_store_async.posix.c @@ -93,10 +93,15 @@ process_job(struct metadata_store_async* s, struct metadata_job* job) started = metadata_monotonic_ns(); int rc = fstat(fd, &st); metadata_record_op_latency(&s->common, OP_STATX, started); - if (rc || st.st_size < 0 || (uint64_t)st.st_size > SIZE_MAX) { + if (rc || st.st_size < 0) { status = DAMACY_IO; goto Close; } + if ((uint64_t)st.st_size > job->requested_len || + (uint64_t)st.st_size > UINT32_MAX) { + status = DAMACY_BUDGET; + goto Close; + } len = (size_t)st.st_size; } if (!len) @@ -292,7 +297,17 @@ metadata_store_async_read_file(struct metadata_store_async* s, metadata_store_read_cb cb, void* user) { - return post_read(s, key, 0, 0, REQ_READ_FILE, cb, user); + return post_read(s, key, 0, SIZE_MAX, REQ_READ_FILE, cb, user); +} + +int +metadata_store_async_read_file_bounded(struct metadata_store_async* s, + const char* key, + size_t max_bytes, + metadata_store_read_cb cb, + void* user) +{ + return post_read(s, key, 0, max_bytes, REQ_READ_FILE, cb, user); } int