diff --git a/docs/api.md b/docs/api.md index 6dd8a136..104b636c 100644 --- a/docs/api.md +++ b/docs/api.md @@ -28,6 +28,22 @@ The full public surface of the `damacy` package. ::: damacy.CudaExecutor +## Spatial queries + +::: damacy.NgffImage + +::: damacy.NgffAxis + +::: damacy.NgffLevel + +::: damacy.SpatialQuery + +::: damacy.ResolvedSpatialQuery + +::: damacy.Sampler + +::: damacy.NgffLimits + ## Output and limits ::: damacy.BatchSpec @@ -84,6 +100,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..0de80168 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](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..aad6dde8 --- /dev/null +++ b/docs/spatial.md @@ -0,0 +1,302 @@ +# 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[Resolve query] + 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") +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 = image.resolve(query, shape=output.shape) + +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=3 << 30), +) +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. `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. 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 + +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 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. +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. 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, 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 + +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. 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 +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. + +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: + +```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. +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 + +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 = {0}; +struct damacy_ngff_limits limits = { + .max_levels = 16, .max_metadata_bytes = 4 << 20 +}; +const int64_t output_shape[] = {64, 64}; +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, &query, 2, output_shape, &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_clear(&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. +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. +`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..9e0d9ee6 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,16 @@ "PoolStarved", "QueueLimits", "RankMismatch", + "ResolvedSpatialQuery", "Sample", + "Sampler", "ShutdownError", + "SpatialQuery", "Stats", "Status", "StorageError", "TryAgain", + "UnsupportedOperation", "ZarrMetadata", "max_concurrency", "set_log_level", @@ -241,6 +249,7 @@ class Status(IntEnum): OOM = _native.STATUS_OOM BUDGET = _native.STATUS_BUDGET SHUTDOWN = _native.STATUS_SHUTDOWN + UNSUPPORTED = _native.STATUS_UNSUPPORTED # ---- exceptions --------------------------------------------------------- @@ -305,6 +314,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 +345,7 @@ class PoolStarved(DamacyError): _native.STATUS_OOM: OutOfMemory, _native.STATUS_BUDGET: BudgetExceeded, _native.STATUS_SHUTDOWN: ShutdownError, + _native.STATUS_UNSUPPORTED: UnsupportedOperation, } @@ -1609,8 +1623,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 +1707,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 +1933,14 @@ 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, +) 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..1a08f07f 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 ---------------------------------------------- @@ -51,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): @@ -252,3 +262,22 @@ 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, ...], + transform: tuple[tuple[float, ...], ...], + filter: int, + boundary: int, + constant_value: float, + level: int, + /, +) -> dict[str, Any]: ... diff --git a/python/damacy/_spatial.c b/python/damacy/_spatial.c new file mode 100644 index 00000000..bc9b6dc5 --- /dev/null +++ b/python/damacy/_spatial.c @@ -0,0 +1,328 @@ +#define PY_SSIZE_T_CLEAN +#include + +#include "_api.h" +#include "damacy_spatial.h" + +static const char image_name[] = "damacy.NgffImage"; + +static void +image_destroy(PyObject* capsule) +{ + damacy_ngff_image_destroy(PyCapsule_GetPointer(capsule, image_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* +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; + int filter, boundary, level; + double value; + if (!PyArg_ParseTuple(args, + "OOOiidi", + &capsule, + &shape_object, + &transform, + &filter, + &boundary, + &value, + &level)) + return NULL; + struct damacy_ngff_image* image = PyCapsule_GetPointer(capsule, image_name); + if (!image) + return NULL; + int64_t output_shape[DAMACY_MAX_RANK]; + 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; + } + for (Py_ssize_t i = 0; i < rank; ++i) { + output_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, (uint8_t)rank, &query.output_to_reference)) + return NULL; + struct damacy_spatial_resolution resolution; + enum damacy_status status; + 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 = resolution_value(&resolution); + damacy_spatial_resolution_clear(&resolution); + 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 }, + { NULL, NULL, 0, NULL } +}; + +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 new file mode 100644 index 00000000..d2336ac7 --- /dev/null +++ b/python/damacy/_spatial.py @@ -0,0 +1,274 @@ +from __future__ import annotations + +import math +import operator +import os +from collections.abc import Iterable +from dataclasses import dataclass +from typing import Any, Literal, NoReturn + +from . import ( + FileMetadataReader, + Status, + UnsupportedOperation, + _component, + _native, + _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): + 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, ...] + + +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(_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. + It finishes before returning. + ``multiscale_index`` explicitly selects an entry in ``ome.multiscales``. + """ + + axes: tuple[NgffAxis, ...] + levels: tuple[NgffLevel, ...] + data_type: str + + 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"]) + + 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: + """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, + _FILTERS[query.sampler.filter], + _BOUNDARIES[query.sampler.boundary], + query.sampler.constant_value, + _native.LEVEL_AUTO if query.level == "auto" else query.level, + ) + return ResolvedSpatialQuery._from_native(info, query.sampler) + + +@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 geometry returned by NgffImage.resolve(). + + Aligned results can be pushed to a Pipeline. Other results raise + UnsupportedOperation when pushed; their geometry remains inspectable. + """ + + 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 + + def __init__(self) -> None: + raise TypeError("use NgffImage.resolve() to create a ResolvedSpatialQuery") + + @classmethod + def _from_native( + cls, info: dict[str, Any], sampler: Sampler + ) -> ResolvedSpatialQuery: + result = object.__new__(cls) + object.__setattr__(result, "sampler", sampler) + for name, value in info.items(): + object.__setattr__(result, name, value) + return result + + def _to_native(self) -> dict[str, Any]: + 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 new file mode 100644 index 00000000..5f4b6cba --- /dev/null +++ b/python/tests/test_spatial.py @@ -0,0 +1,843 @@ +from __future__ import annotations + +import copy +import ctypes +import dataclasses +import gc +import json +import pickle +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 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, + reader=damacy.FileMetadataReader(concurrency=2), + multiscale_index=0, + **kwargs, + ) + + +def query( + scale: float = 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)) + shape = output.shape + factor = 2**level + resolved = image.resolve(query(factor, (factor * 2, factor * 3)), shape=shape) + assert resolved.level == level + assert not resolved.requires_resampling + assert resolved.source_bounds_index == ((2, 6), (3, 7)) + assert resolved.output_to_source == ((1, 0, 2), (0, 1, 3)) + del image + gc.collect() + with pipeline(spatial_executor, output, reader) as p: + 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( + 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)) + 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 pipeline(spatial_executor, output) as p: + with pytest.raises(damacy.UnsupportedOperation, match="resampling required"): + p.push([rotated]) + 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]) + + +def test_fixed_shape_is_checked_at_push(tmp_path, spatial_executor): + root = tmp_path / "image" + write_image(root) + resolved = load(root).resolve(query(), shape=(3, 3)) + 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) + shape = (4, 4) + resolved = image.resolve(query(2), shape=shape) + assert resolved.level == 1 and not resolved.requires_resampling + 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( + 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) + 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 = 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) + 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 ( + image.resolve( + damacy.SpatialQuery( + output_to_reference=transform, sampler=damacy.Sampler() + ), + shape=shape, + ).level + == 0 + ) + transform[0, 0] = 2 + with pytest.raises(damacy.InvalidArgument): + 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) + image = load(root) + shape = (2, 2) + q = damacy.SpatialQuery( + output_to_reference=((3, 2, 0), (2, 3, 0)), + sampler=damacy.Sampler(filter="nearest"), + ) + assert image.resolve(q, shape=shape).level == 0 + forced = dataclasses.replace(q, level=1) + assert image.resolve(forced, shape=shape).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) + 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 + 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) + image = load(root) + shape = (4, 4) + with pytest.raises(damacy.InvalidArgument): + image.resolve( + damacy.SpatialQuery( + output_to_reference=transform, sampler=damacy.Sampler() + ), + shape=shape, + ) + + +@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) + + +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", + [ + 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[:-1] + ",}", + ], +) +def test_invalid_consumed_metadata_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) + 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(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_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()) + 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)) + image = load(root) + shape = output.shape + transform = np.column_stack((np.diag([1, 1, 1, 2, 4]), [1, 1, 4, 4, 0])) + 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: + 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()) + 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))) + 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 = 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 + + +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): + 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)) + 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) + + +def test_collapsed_3d_transforms_are_rejected(tmp_path): + root = tmp_path / "image" + write_image(root, shape=(16, 16, 16)) + image = load(root) + shape = (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): + 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) + 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 = 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/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..2182b775 --- /dev/null +++ b/src/damacy_spatial.h @@ -0,0 +1,125 @@ +#pragma once + +#include "damacy_pipeline.h" + +#ifdef __cplusplus +extern "C" +{ +#endif + + struct damacy_ngff_image; + + 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_resolution + { + 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_spatial_query* query, + 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_clear( + 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..447088c1 --- /dev/null +++ b/src/ngff/load.c @@ -0,0 +1,127 @@ +#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; + 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); + 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..8fae9706 --- /dev/null +++ b/src/ngff/metadata.c @@ -0,0 +1,510 @@ +#include "ngff/ngff.h" + +#include "zarr/zarr_metadata.h" + +#include +#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) { + const char* text = NULL; + enum damacy_status status = read_string(name, &text); + if (status == DAMACY_OOM) + return JSON_ERR_OOM; + int matches = status == DAMACY_OK && !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_object(struct cslice src, struct json_node* root) +{ + return json_resolve(src, NULL, 0, root, NULL) || root->type != JSON_OBJECT; +} + +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); +} + +// 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) +{ + 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; + 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; + 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_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) || + 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; + enum json_err error; + 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_object(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..e8e128bc --- /dev/null +++ b/src/query/spatial.c @@ -0,0 +1,301 @@ +#include "damacy_spatial.h" + +#include "ngff/ngff.h" + +#include +#include +#include +#include + +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; + for (uint8_t i = 0; i < rank; ++i) + 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]]; + 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; + for (uint8_t i = 0; i < rank; ++i) + for (uint8_t j = 0; j < rank; ++j) + 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 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 = (b - a) / (2 * cross); + double t = copysign(1, tau) / (fabs(tau) + hypot(1, tau)); + double c = 1 / hypot(1, t), s = t * c; + 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); + } + return scale * minimum; +} + +static enum damacy_status +validate_query(const struct damacy_ngff_image* image, + const struct damacy_spatial_query* query, + uint8_t rank, + const int64_t* output_shape) +{ + if (!image || !output_shape || !query) + return DAMACY_INVAL; + if (rank != image->info.rank) + return DAMACY_RANK; + 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_shape[i] <= 0 || + output_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; + } + } + if (!linear_is_invertible(image, &query->output_to_reference)) + 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_resolution* 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_spatial_query* query, + uint8_t rank, + const int64_t* output_shape, + struct damacy_spatial_resolution* out) +{ + if (!out) + return DAMACY_INVAL; + *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; + 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 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]; + 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, &result); + if (status != DAMACY_OK) + return status; + char* uri = malloc(strlen(level->uri) + 1); + if (!uri) + return DAMACY_OOM; + strcpy(uri, level->uri); + result.uri = uri; + *out = result; + return DAMACY_OK; +} + +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 || !resolution->uri || !resolution->rank || + resolution->rank > DAMACY_MAX_RANK) + return DAMACY_INVAL; + if (resolution->requires_resampling) + return DAMACY_UNSUPPORTED; + 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_clear(struct damacy_spatial_resolution* resolution) +{ + if (resolution) { + free((void*)resolution->uri); + *resolution = (struct damacy_spatial_resolution){ 0 }; + } +} 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/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 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..4fc94646 --- /dev/null +++ b/tests/test_spatial.c @@ -0,0 +1,507 @@ +#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 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) +{ + 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 const int64_t output_shape[] = { 4, 4 }; + +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, &query, 2, output_shape, &resolved) == + DAMACY_OK); + EXPECT(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_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, &query, 2, output_shape, &resolved) == + DAMACY_OK); + 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(!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_clear(&resolved); + 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) +{ + 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, &query, 2, output_shape, &resolved) == + DAMACY_OK); + 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) == + DAMACY_UNSUPPORTED); + EXPECT(!sample.uri && !sample.rank); + 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, &query, 2, output_shape, &resolved) == + DAMACY_OK); + 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, &query, 2, output_shape, &resolved) == + DAMACY_OK); + 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, &query, 2, output_shape, &resolved) == + DAMACY_OK); + EXPECT(resolved.level == 0); + damacy_spatial_resolution_clear(&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, &query, 2, output_shape, &resolved) == + DAMACY_OK); + 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, &query, 2, output_shape, &resolved) == + DAMACY_INVAL); + EXPECT(!resolved.uri && !resolved.rank); + query.sampler.boundary = DAMACY_BOUNDARY_CONSTANT; + query.sampler.constant_value = 17; + EXPECT(damacy_spatial_resolve(image, &query, 2, output_shape, &resolved) == + DAMACY_OK); + 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, &query, 2, output_shape, &resolved) == + DAMACY_OK); + 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, &query, 2, output_shape, &resolved) == + DAMACY_OK); + 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; +} + +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, &query, 2, output_shape, &resolved) == + DAMACY_INVAL); + EXPECT(!resolved.uri && !resolved.rank); + } + query = identity_query(); + query.level = -2; + EXPECT(damacy_spatial_resolve(image, &query, 2, output_shape, &resolved) == + DAMACY_INVAL); + query.level = 3; + 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, &query, 2, output_shape, &resolved) == + DAMACY_INVAL); + query.sampler.constant_value = 0; + query.sampler.filter = 0; + EXPECT(damacy_spatial_resolve(image, &query, 2, output_shape, &resolved) == + DAMACY_INVAL); + query = identity_query(); + 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_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) +{ + 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; +} + +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 + }; + const int64_t shape[] = { 2, 2, 2 }; + 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; + EXPECT(damacy_spatial_resolve(&image, &query, 3, shape, &resolved) == + DAMACY_INVAL); + 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, &query, 3, shape, &resolved) == + DAMACY_OK); + EXPECT(!resolved.requires_resampling); + damacy_spatial_resolution_clear(&resolved); + return 0; +} + +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); + RUN(test_collapsed_volume); + RUN(test_resolution_shape_and_clear); + RUN(test_unrelated_keys_are_skipped); + RUN(test_metadata_limits_and_load); + return 0; +}