From acd89cc98838f0cd491bfa0fa6d5ebc9a5f3234c Mon Sep 17 00:00:00 2001 From: Peter Grotz <199680633+peter-grotz@users.noreply.github.com> Date: Tue, 22 Sep 2026 14:39:00 -0700 Subject: [PATCH 1/2] feat(registration): recover registration geometry without the reference volumes The transform needs the geometry of the two volumes the CCF registration worked from. Those are published under ccf_alignment/registration_metadata/ at 1.4-1.8 GB each, but only for some samples: of 60 processed exaSPIM assets with a fused zarr, 40 have them. The rest cannot be transformed at all today. Nothing in the transform reads their voxels -- check_orientation needs input_img.shape and apply_transforms_to_points needs the resampled image's spacing, origin and direction. The only pixel consumer is the QC overlay, which this pipeline does not produce. Two paths, because the samples split cleanly and disjointly. Older samples publish the volumes but carry a v1 processing record with empty parameters; newer ones carry a full v2 record but no volumes. Where the volumes exist their headers are read from a 200 KB byte range rather than the whole file -- exact on all 40. Where the record exists the geometry is derived from the zarr level it names and the sample_scale it recorded. Where both exist, reconcile() checks one against the other and fails rather than mis-registering. Three details are load-bearing, each verified across the 40 and each silently wrong if missed: - the source is fused_ccf_ch.zarr, not fused.zarr; rules from the signal pyramid match 11-27 of 40, the CCF channel matches 40/40 - resampled shapes need exact rationals: 0.0135/0.01 is 1.3499999999999999 in binary float, so a 1270-voxel axis rounds to 1714 where the registration produced 1715 - the direction matrix is ANTs/ITK LPS, which negates the first two rows of the NIfTI RAS srow The written spacing is sample_scale reversed, which does not correspond to the array's physical axes -- the coarse value lands on axis 2 while the coarse axis is z at index 1. Reproducing it is required: the transforms were fit in that space. Also recovers the template-to-CCF version the registration used, so a sample registered against v1.4 is not silently composed with v1.5. --- src/exaspim_swc_processing/registration.py | 402 ++++++++++++++++++ tests/resources/ccf_alignment_processing.json | 265 ++++++++++++ tests/test_registration.py | 254 +++++++++++ 3 files changed, 921 insertions(+) create mode 100644 src/exaspim_swc_processing/registration.py create mode 100644 tests/resources/ccf_alignment_processing.json create mode 100644 tests/test_registration.py diff --git a/src/exaspim_swc_processing/registration.py b/src/exaspim_swc_processing/registration.py new file mode 100644 index 0000000..95027f4 --- /dev/null +++ b/src/exaspim_swc_processing/registration.py @@ -0,0 +1,402 @@ +"""Recover the CCF registration's own parameters and reconstruct its reference geometry. + +Transforming reconstructions into CCF space needs the geometry of the two volumes the +registration worked from. Those are published under +``/ccf_alignment/registration_metadata/`` as NIfTI files of 1.4-1.8 GB each -- +but only sometimes: of 60 processed exaSPIM assets with a fused zarr, 40 have them. The +rest cannot be transformed at all by code that insists on reading the files. + +Nothing in the transform needs the voxels. ``check_orientation`` mirrors coordinates about +``input_img.shape``, and ``apply_transforms_to_points`` converts indices to physical +coordinates using the resampled image's spacing, origin and direction. Both are geometry. +The only consumer of pixel data is the QC overlay rendering, which this pipeline does not +produce. + +So the geometry is rebuilt from metadata the registration already publishes: +``/ccf_alignment/processing.json`` records the zarr it read, the pyramid level, +and the ``sample_scale`` it assigned. Shapes come from that zarr's ``.zarray``. + +Three details are load-bearing, each verified against all 40 samples that publish the real +volumes. Getting any of them wrong yields plausible coordinates rather than an error: + +* The source is ``fused_ccf_ch.zarr``, not ``fused.zarr``. Rules derived from the signal + pyramid match only 11-27 of 40; the CCF channel matches 40/40. +* The resampled shape must be computed with exact rationals. ``0.0135 / 0.01`` evaluates + to ``1.3499999999999999`` in binary floating point, so a shape of 1270 rounds down to + 1714 where the registration produced 1715. +* The direction matrix is expressed in ANTs/ITK LPS, which negates the first two rows of + the NIfTI RAS ``srow``. +""" + +from __future__ import annotations + +import logging +import math +import re +import struct +import zlib +from dataclasses import dataclass +from fractions import Fraction +from pathlib import PurePosixPath + +logger = logging.getLogger(__name__) + +CCF_ALIGNMENT_DIR = "ccf_alignment" +"""Directory under a processed dataset holding the registration's outputs.""" + +PROCESSING_RECORD = f"{CCF_ALIGNMENT_DIR}/processing.json" +"""The registration's own record, relative to the dataset root.""" + +RESAMPLED_SPACING_MM = (0.01, 0.01, 0.01) +"""Isotropic spacing the registration resamples to, in millimetres.""" + +ANTS_DIRECTION = ( + (0.0, 0.0, -1.0), + (1.0, 0.0, 0.0), + (0.0, -1.0, 0.0), +) +"""Direction matrix as ANTs reports it, in LPS. Not the NIfTI RAS ``srow``.""" + +ORIGIN_MM = (0.0, 0.0, 0.0) +"""Origin of both reference volumes.""" + +_TEN_UM = 10 +"""``resolution_um`` of the pass whose geometry the SWC transform consumes.""" + + +class RegistrationRecordError(ValueError): + """Raised when a registration record is missing or does not describe a usable pass.""" + + +@dataclass(frozen=True) +class VolumeGeometry: + """Geometry of one reference volume, without its voxels. + + Attributes + ---------- + shape : tuple[int, int, int] + Array shape in the order the registration wrote it, ``(x, z, y)``. + spacing_mm : tuple[float, float, float] + Voxel size in millimetres, matching ``shape``. + """ + + shape: tuple[int, int, int] + spacing_mm: tuple[float, float, float] + + +@dataclass(frozen=True) +class RegistrationPass: + """One registration pass, as the registration recorded it. + + Attributes + ---------- + input_uri : str + The zarr the registration read, e.g. ``s3://.../fusion/fused_ccf_ch.zarr/``. + level : int + Pyramid level within that zarr. + resolution_um : int + Target resolution of the pass, 10 or 25. + sample_scale_mm : tuple[float, float, float] + Voxel size the registration assigned, in ``(z, y, x)`` order as recorded. + sample_to_template : tuple[str, ...] + Basenames of the sample-to-template transforms, resolved against + ``/ccf_alignment/``. The recorded paths point at Nextflow scratch. + template_to_ccf : tuple[str, ...] + Recorded paths of the template-to-CCF transforms. These are Code Ocean data + assets with no S3 copy, so only the version directory is usable. + """ + + input_uri: str + level: int + resolution_um: int + sample_scale_mm: tuple[float, float, float] + sample_to_template: tuple[str, ...] + template_to_ccf: tuple[str, ...] + + @property + def loaded_spacing_mm(self) -> tuple[float, float, float]: + """Spacing the registration wrote into the loaded volume's header. + + ``sample_scale`` is recorded ``(z, y, x)``, coarse axis first, and the written + header is simply that reversed. Note this does **not** correspond to the array's + physical axes: the volume is stored ``(x, z, y)``, so the coarse value lands on + axis 2 while the coarse physical axis is z at index 1. Verified against the real + headers of 40 samples, which all record ``(0.010125, 0.010125, 0.0135)``. + + Reproducing this exactly is required, not optional -- the transforms were fit in + this space, so a "corrected" spacing would invalidate them. + + Returns + ------- + tuple[float, float, float] + Spacing in millimetres, matching the written array axes. + """ + return tuple(reversed(self.sample_scale_mm)) + + @property + def template_to_ccf_version(self) -> str | None: + """Name of the template-to-CCF asset the registration used. + + Returns + ------- + str | None + The directory name, e.g. ``"reg_exaspim_template_to_ccf_25um_v1.4"``, or + ``None`` if no transform was recorded. + """ + for path in self.template_to_ccf: + parts = PurePosixPath(path).parts + for part in parts: + if part.startswith("reg_exaspim_template_to_ccf"): + return part + return None + + +def _as_float_triple(values: object, field: str) -> tuple[float, float, float]: + """Coerce a recorded three-element sequence to floats. + + Parameters + ---------- + values : object + The recorded value. + field : str + Field name, for the error message. + + Returns + ------- + tuple[float, float, float] + The three values. + + Raises + ------ + RegistrationRecordError + If the value is not a sequence of three numbers. + """ + if not isinstance(values, (list, tuple)) or len(values) != 3: + raise RegistrationRecordError(f"{field} is not a three-element sequence: {values!r}") + return (float(values[0]), float(values[1]), float(values[2])) + + +def parse_registration_record( + record: dict, + resolution_um: int = _TEN_UM, +) -> RegistrationPass: + """Extract one registration pass from ``ccf_alignment/processing.json``. + + The registration runs two passes, coarse then fine. The SWC transform consumes the + fine one. + + Parameters + ---------- + record : dict + The decoded processing record. + resolution_um : int, optional + Which pass to take, by default 10. + + Returns + ------- + RegistrationPass + The recorded parameters of that pass. + + Raises + ------ + RegistrationRecordError + If no pass at that resolution is present, or it omits a required field. + """ + for process in record.get("data_processes") or []: + parameters = ((process.get("code") or {}).get("parameters")) or {} + if parameters.get("resolution_um") != resolution_um: + continue + registration = parameters.get("sample_to_template_registration") or {} + missing = [k for k in ("input_uri", "level") if parameters.get(k) is None] + if missing: + raise RegistrationRecordError( + f"{resolution_um} um pass is missing {', '.join(missing)}" + ) + return RegistrationPass( + input_uri=str(parameters["input_uri"]), + level=int(parameters["level"]), + resolution_um=resolution_um, + sample_scale_mm=_as_float_triple(registration.get("sample_scale"), "sample_scale"), + sample_to_template=tuple( + PurePosixPath(p).name for p in parameters.get("sample_to_template_transforms") or () + ), + template_to_ccf=tuple(parameters.get("template_to_ccf_transforms") or ()), + ) + available = sorted( + { + ((p.get("code") or {}).get("parameters") or {}).get("resolution_um") + for p in record.get("data_processes") or [] + } + - {None} + ) + raise RegistrationRecordError( + f"No {resolution_um} um registration pass in the record; found {available or 'none'}" + ) + + +def loaded_geometry(pass_: RegistrationPass, zarr_shape: tuple[int, int, int]) -> VolumeGeometry: + """Geometry of the volume the registration loaded. + + Parameters + ---------- + pass_ : RegistrationPass + The recorded pass. + zarr_shape : tuple[int, int, int] + Shape of ``/``, in the zarr's ``(z, y, x)`` order. + + Returns + ------- + VolumeGeometry + Shape reordered to the written ``(x, z, y)``, with the recorded spacing. + """ + z, y, x = zarr_shape + return VolumeGeometry(shape=(x, z, y), spacing_mm=pass_.loaded_spacing_mm) + + +def resampled_geometry(loaded: VolumeGeometry) -> VolumeGeometry: + """Geometry of the isotropically resampled volume. + + Shapes are computed with exact rationals and rounded half up. Binary floating point + puts ``0.0135 / 0.01`` just below 1.35, which rounds a 1270-voxel axis down to 1714 + where the registration produced 1715. + + Parameters + ---------- + loaded : VolumeGeometry + Geometry of the loaded volume. + + Returns + ------- + VolumeGeometry + Shape scaled to isotropic spacing, with spacing :data:`RESAMPLED_SPACING_MM`. + """ + factors = [ + Fraction(spacing).limit_denominator(10**6) / Fraction(target).limit_denominator(10**6) + for spacing, target in zip(loaded.spacing_mm, RESAMPLED_SPACING_MM) + ] + shape = tuple( + math.floor(Fraction(extent) * factor + Fraction(1, 2)) + for extent, factor in zip(loaded.shape, factors) + ) + return VolumeGeometry(shape=shape, spacing_mm=RESAMPLED_SPACING_MM) + + +def zarr_level_key(pass_: RegistrationPass) -> str: + """Key of the ``.zarray`` describing the level the registration read. + + Parameters + ---------- + pass_ : RegistrationPass + The recorded pass. + + Returns + ------- + str + Key relative to the bucket, e.g. + ``exaSPIM_..._processed_.../fusion/fused_ccf_ch.zarr/2/.zarray``. + """ + key = re.sub(r"^s3://[^/]+/", "", pass_.input_uri).rstrip("/") + return f"{key}/{pass_.level}/.zarray" + + +NIFTI_HEADER_SIZE = 348 +"""Bytes of a NIfTI-1 header, which carries everything needed here.""" + +_HEADER_FETCH_BYTES = 200_000 +"""Compressed bytes to request. Ample for the header of a gzipped NIfTI.""" + + +def parse_nifti_geometry(compressed_header: bytes) -> VolumeGeometry: + """Read shape and spacing from the leading bytes of a gzipped NIfTI. + + Only the header is needed, so callers can fetch a byte range rather than a volume of + well over a gigabyte. + + Parameters + ---------- + compressed_header : bytes + Leading bytes of the ``.nii.gz``. Need not be the whole file. + + Returns + ------- + VolumeGeometry + The volume's shape and spacing. + + Raises + ------ + RegistrationRecordError + If the bytes do not decompress to a readable NIfTI-1 header. + """ + try: + buffer = zlib.decompressobj(16 + zlib.MAX_WBITS).decompress(compressed_header, 4096) + except zlib.error as error: + raise RegistrationRecordError(f"Could not decompress the NIfTI header: {error}") from error + if len(buffer) < NIFTI_HEADER_SIZE: + raise RegistrationRecordError( + f"Decompressed only {len(buffer)} bytes, need {NIFTI_HEADER_SIZE}" + ) + endian = "<" if struct.unpack("" + dims = struct.unpack(f"{endian}8h", buffer[40:56]) + pixdim = struct.unpack(f"{endian}8f", buffer[76:108]) + if dims[0] < 3: + raise RegistrationRecordError(f"NIfTI header declares {dims[0]} dimensions, need 3") + return VolumeGeometry( + shape=(int(dims[1]), int(dims[2]), int(dims[3])), + spacing_mm=(float(pixdim[1]), float(pixdim[2]), float(pixdim[3])), + ) + + +def registration_volume_key(dataset: str, dataset_id: str, volume: str) -> str: + """Key of a published reference volume, if the registration wrote one. + + Parameters + ---------- + dataset : str + Processed dataset name. + dataset_id : str + Numeric subject id. + volume : str + Either ``"loaded"`` or ``"resampled"``. + + Returns + ------- + str + Key relative to the bucket. + """ + return ( + f"{dataset}/{CCF_ALIGNMENT_DIR}/registration_metadata/" + f"{dataset_id}_10um_{volume}_zarr_img.nii.gz" + ) + + +def reconcile(derived: VolumeGeometry, published: VolumeGeometry, volume: str) -> None: + """Check a derived geometry against the one the registration published. + + Where both are available this turns the derivation from an assumption into a checked + fact. A disagreement means the registration changed and the derivation is stale -- + which would otherwise yield plausible coordinates rather than an error. + + Parameters + ---------- + derived : VolumeGeometry + Geometry reconstructed from metadata. + published : VolumeGeometry + Geometry read from the published volume's header. + volume : str + Which volume, for the error message. + + Raises + ------ + RegistrationRecordError + If shape differs, or spacing differs by more than float32 precision. + """ + if derived.shape != published.shape: + raise RegistrationRecordError( + f"Derived {volume} shape {derived.shape} does not match the published " + f"{published.shape}; the registration's geometry has changed" + ) + if any(abs(a - b) > 1e-6 for a, b in zip(derived.spacing_mm, published.spacing_mm)): + raise RegistrationRecordError( + f"Derived {volume} spacing {derived.spacing_mm} does not match the published " + f"{published.spacing_mm}; the registration's geometry has changed" + ) diff --git a/tests/resources/ccf_alignment_processing.json b/tests/resources/ccf_alignment_processing.json new file mode 100644 index 0000000..6ffd6aa --- /dev/null +++ b/tests/resources/ccf_alignment_processing.json @@ -0,0 +1,265 @@ +{ + "object_type": "Processing", + "describedBy": "https://raw.githubusercontent.com/AllenNeuralDynamics/aind-data-schema/main/src/aind_data_schema/core/processing.py", + "schema_version": "2.2.8", + "data_processes": [ + { + "object_type": "Data process", + "process_type": "Image atlas alignment", + "name": "Image atlas alignment - 25 um", + "stage": "Processing", + "code": { + "object_type": "Code", + "url": "https://codeocean.allenneuraldynamics.org/capsule/6898460/tree", + "name": "aind-exaspim-ccf-registration", + "version": "0.0.1", + "commit_hash": null, + "container": null, + "run_script": "/code/run", + "language": "Python", + "language_version": null, + "input_data": null, + "parameters": { + "input_uri": "s3://aind-open-data/exaSPIM_826509_2026-05-12_18-03-50_processed_2026-05-22_14-04-22/fusion/fused_ccf_ch.zarr/", + "level": 3, + "resolution_um": 25, + "dataset_id": "826509", + "use_fused_mask": true, + "mask_applied": true, + "sample_to_template_registration": { + "sample_scale": [ + 0.027, + 0.02025, + 0.02025 + ], + "exaspim_template_path": "../data/exaSPIM_template_25um/exaspim_template_7sujects_nomask_25um_round6.nii.gz", + "ccf_path": "../data/allen_mouse_ccf/average_template/average_template_25.nii.gz", + "affine_reg_iterations": [ + 200, + 100, + 25, + 3 + ], + "syn_reg_iterations": [ + 200, + 100, + 25, + 3 + ], + "log_level": "ERROR" + }, + "template_to_ccf_transforms": [ + "/data/reg_exaspim_template_to_ccf_25um_v1.4/1Warp.nii.gz", + "/data/reg_exaspim_template_to_ccf_25um_v1.4/0GenericAffine.mat" + ] + }, + "core_dependency": null + }, + "experimenters": [ + "Peter Grotz" + ], + "pipeline_name": "exaspim-data-processing", + "start_date_time": "2026-08-04T19:10:37.639249Z", + "end_date_time": "2026-08-04T19:23:26.784422Z", + "output_path": "ccf_alignment", + "output_parameters": null, + "notes": "Import + template-based registration (25 um): sample -> template -> CCF", + "resources": null + }, + { + "object_type": "Data process", + "process_type": "Image atlas alignment", + "name": "Image atlas alignment - 10 um", + "stage": "Processing", + "code": { + "object_type": "Code", + "url": "https://codeocean.allenneuraldynamics.org/capsule/6898460/tree", + "name": "aind-exaspim-ccf-registration", + "version": "0.0.1", + "commit_hash": null, + "container": null, + "run_script": "/code/run", + "language": "Python", + "language_version": null, + "input_data": null, + "parameters": { + "input_uri": "s3://aind-open-data/exaSPIM_826509_2026-05-12_18-03-50_processed_2026-05-22_14-04-22/fusion/fused_ccf_ch.zarr/", + "level": 2, + "resolution_um": 10, + "dataset_id": "826509_10um", + "sample_to_template_registration": { + "sample_scale": [ + 0.0135, + 0.010125, + 0.010125 + ], + "exaspim_template_path": "../data/exaspim_template_7subjects_nomask_10um_round6_template_only/fixed_median.nii.gz", + "ccf_path": "../data/allen_mouse_ccf/average_template/average_template_10.nii.gz", + "affine_reg_iterations": [ + 300, + 200, + 50, + 0 + ], + "syn_reg_iterations": [ + 400, + 200, + 40, + 0 + ], + "log_level": "ERROR" + }, + "sample_to_template_transforms": [ + "/tmp/nxf.leSOsXxjX4/capsule/results/ccf_alignment/826509_to_exaSPIM_SyN_1Warp.nii.gz", + "/tmp/nxf.leSOsXxjX4/capsule/results/ccf_alignment/826509_to_exaSPIM_SyN_0GenericAffine.mat" + ], + "template_to_ccf_transforms": [ + "/data/reg_exaspim_template_to_ccf_25um_v1.4/1Warp.nii.gz", + "/data/reg_exaspim_template_to_ccf_25um_v1.4/0GenericAffine.mat" + ] + }, + "core_dependency": null + }, + "experimenters": [ + "Peter Grotz" + ], + "pipeline_name": "exaspim-data-processing", + "start_date_time": "2026-08-04T19:23:46.693066Z", + "end_date_time": "2026-08-04T19:51:21.310082Z", + "output_path": "ccf_alignment", + "output_parameters": null, + "notes": "Import + 10um sample->CCF registration using transforms: ['/tmp/nxf.leSOsXxjX4/capsule/results/ccf_alignment/826509_to_exaSPIM_SyN_1Warp.nii.gz', '/tmp/nxf.leSOsXxjX4/capsule/results/ccf_alignment/826509_to_exaSPIM_SyN_0GenericAffine.mat'] and ['/data/reg_exaspim_template_to_ccf_25um_v1.4/1Warp.nii.gz', '/data/reg_exaspim_template_to_ccf_25um_v1.4/0GenericAffine.mat']", + "resources": null + }, + { + "object_type": "Data process", + "process_type": "Image atlas alignment", + "name": "CCF annotation to sample space", + "stage": "Processing", + "code": { + "object_type": "Code", + "url": "https://codeocean.allenneuraldynamics.org/capsule/6898460/tree", + "name": "aind-exaspim-ccf-registration", + "version": "0.0.1", + "commit_hash": null, + "container": null, + "run_script": "/code/run_register_ccf_annotation_10um.sh", + "language": "Python", + "language_version": null, + "input_data": null, + "parameters": { + "subject_id": "826509", + "ccf_annotation_path": "../data/allen_mouse_ccf/annotation/ccf_2017/annotation_10.nii.gz", + "ccf_template_path": "../data/allen_mouse_ccf/average_template/average_template_10.nii.gz", + "exaspim_template_path": "../data/exaspim_template_7subjects_nomask_10um_round6_template_only/fixed_median.nii.gz", + "level": 2 + }, + "core_dependency": null + }, + "experimenters": [ + "Peter Grotz" + ], + "pipeline_name": "exaspim-data-processing", + "start_date_time": "2026-08-04T20:23:03.151214Z", + "end_date_time": "2026-08-04T20:36:13.842546Z", + "output_path": "ccf_alignment/ccf_anno_to_sample", + "output_parameters": { + "standard_outputs": { + "annotation_image": "ccf_alignment/ccf_anno_to_sample/ccf_anno_in_sample_space.nii.gz", + "annotation_zarr": "ccf_alignment/ccf_anno_to_sample/ccf_anno_in_sample_space.zarr" + } + }, + "notes": "Transforms CCF annotation into sample space and writes zarr/precomputed label volumes.", + "resources": null + }, + { + "object_type": "Data process", + "process_type": "Image atlas alignment", + "name": "CCF meshes to sample space", + "stage": "Processing", + "code": { + "object_type": "Code", + "url": "https://codeocean.allenneuraldynamics.org/capsule/6898460/tree", + "name": "aind-exaspim-ccf-registration", + "version": "0.0.1", + "commit_hash": null, + "container": null, + "run_script": "/code/run_invert_ccf_objs.sh", + "language": "Python", + "language_version": null, + "input_data": null, + "parameters": { + "n_meshes": 840, + "mesh_units_um": 1.0, + "ccf_to_template_fwd": [ + "/data/reg_exaspim_template_to_ccf_25um_v1.4/1Warp.nii.gz", + "/data/reg_exaspim_template_to_ccf_25um_v1.4/0GenericAffine.mat" + ], + "template_to_sample_fwd": [ + "/results/ccf_alignment/826509_to_exaSPIM_SyN_1Warp.nii.gz", + "/results/ccf_alignment/826509_to_exaSPIM_SyN_0GenericAffine.mat" + ], + "direction": "reg-FWD(C->T) then SyN-FWD(T->S), two-stage points", + "inv_swaps": [ + [ + 1, + 0 + ], + [ + 2, + 1 + ], + [ + 0, + 2 + ] + ], + "inv_flips": [ + 0, + 1 + ] + }, + "core_dependency": null + }, + "experimenters": [ + "Peter Grotz" + ], + "pipeline_name": "exaspim-data-processing", + "start_date_time": "2026-08-04T21:26:48.231401Z", + "end_date_time": "2026-08-04T21:28:47.269108Z", + "output_path": "ccf_alignment/ccf_mesh_to_sample", + "output_parameters": null, + "notes": "Reverse-transforms CCF region surface meshes into sample space (point transform of vertices through the CCF->template->sample chain).", + "resources": null + } + ], + "pipelines": [ + { + "object_type": "Code", + "url": "https://codeocean.allenneuraldynamics.org/capsule/9578158/tree", + "name": "exaspim-data-processing", + "version": "0.0.0", + "commit_hash": null, + "container": null, + "run_script": null, + "language": null, + "language_version": null, + "input_data": null, + "parameters": null, + "core_dependency": null + } + ], + "notes": "Processes were reordered by start_date_time", + "dependency_graph": { + "CCF annotation to sample space": [ + "Image atlas alignment - 25 um" + ], + "CCF meshes to sample space": [ + "Image atlas alignment - 25 um" + ], + "Image atlas alignment - 10 um": [ + "Image atlas alignment - 25 um" + ], + "Image atlas alignment - 25 um": [] + } +} \ No newline at end of file diff --git a/tests/test_registration.py b/tests/test_registration.py new file mode 100644 index 0000000..99f95fc --- /dev/null +++ b/tests/test_registration.py @@ -0,0 +1,254 @@ +"""Tests for :mod:`exaspim_swc_processing.registration`.""" + +import gzip +import json +import struct +from pathlib import Path + +import pytest + +from exaspim_swc_processing.registration import ( + ANTS_DIRECTION, + NIFTI_HEADER_SIZE, + RegistrationRecordError, + VolumeGeometry, + loaded_geometry, + parse_nifti_geometry, + parse_registration_record, + reconcile, + registration_volume_key, + resampled_geometry, + zarr_level_key, +) + +FIXTURE = Path(__file__).parent / "resources" / "ccf_alignment_processing.json" +DATASET = "exaSPIM_826509_2026-05-12_18-03-50_processed_2026-05-22_14-04-22" + + +def _record() -> dict: + """Load the real registration record for 826509. + + Returns + ------- + dict + The decoded ``ccf_alignment/processing.json``. + """ + return json.loads(FIXTURE.read_text(encoding="utf-8")) + + +def test_the_ten_micron_pass_is_selected() -> None: + """The registration runs coarse then fine; the SWC transform uses the fine pass.""" + result = parse_registration_record(_record()) + assert result.resolution_um == 10 + assert result.level == 2 + assert result.input_uri.endswith("fusion/fused_ccf_ch.zarr/") + + +def test_the_coarse_pass_can_be_selected_too() -> None: + """The 25 um pass is recorded at a different level.""" + result = parse_registration_record(_record(), resolution_um=25) + assert result.level == 3 + assert result.sample_scale_mm == (0.027, 0.02025, 0.02025) + + +def test_written_spacing_is_sample_scale_reversed() -> None: + """The header is ``sample_scale`` reversed, matching all 40 real volumes. + + This puts the coarse value on axis 2 although the coarse physical axis is z at + index 1 -- an upstream quirk that must be reproduced, since the transforms were fit + in this space. + """ + result = parse_registration_record(_record()) + assert result.sample_scale_mm == (0.0135, 0.010125, 0.010125) + assert result.loaded_spacing_mm == (0.010125, 0.010125, 0.0135) + + +def test_transform_basenames_are_extracted() -> None: + """Recorded paths point at Nextflow scratch, so only the basenames are usable.""" + result = parse_registration_record(_record()) + assert result.sample_to_template == ( + "826509_to_exaSPIM_SyN_1Warp.nii.gz", + "826509_to_exaSPIM_SyN_0GenericAffine.mat", + ) + + +def test_the_template_version_is_recovered() -> None: + """826509 registered against v1.4 while the pipeline mounts v1.5.""" + result = parse_registration_record(_record()) + assert result.template_to_ccf_version == "reg_exaspim_template_to_ccf_25um_v1.4" + + +def test_a_missing_template_transform_yields_no_version() -> None: + """A record without the transform does not invent a version.""" + record = _record() + for process in record["data_processes"]: + params = (process.get("code") or {}).get("parameters") or {} + params.pop("template_to_ccf_transforms", None) + assert parse_registration_record(record).template_to_ccf_version is None + + +def test_an_absent_pass_is_reported_with_what_was_found() -> None: + """A record with no matching pass names the resolutions it does have.""" + with pytest.raises(RegistrationRecordError, match="found"): + parse_registration_record(_record(), resolution_um=99) + + +def test_a_pass_missing_required_fields_is_rejected() -> None: + """A truncated record fails loudly rather than yielding partial geometry.""" + record = _record() + for process in record["data_processes"]: + params = (process.get("code") or {}).get("parameters") or {} + if params.get("resolution_um") == 10: + params.pop("level") + with pytest.raises(RegistrationRecordError, match="level"): + parse_registration_record(record) + + +def test_a_malformed_sample_scale_is_rejected() -> None: + """A scale that is not three numbers is an error, not a silent default.""" + record = _record() + for process in record["data_processes"]: + params = (process.get("code") or {}).get("parameters") or {} + if params.get("resolution_um") == 10: + params["sample_to_template_registration"]["sample_scale"] = [0.01, 0.01] + with pytest.raises(RegistrationRecordError, match="sample_scale"): + parse_registration_record(record) + + +def test_zarr_level_key_points_at_the_level_the_registration_read() -> None: + """The key is bucket-relative, so it can be fetched directly.""" + result = parse_registration_record(_record()) + assert zarr_level_key(result) == f"{DATASET}/fusion/fused_ccf_ch.zarr/2/.zarray" + + +def test_loaded_shape_is_the_zarr_shape_reordered() -> None: + """Verified on 40 samples: the written volume is the zarr level as (x, z, y).""" + result = parse_registration_record(_record()) + geometry = loaded_geometry(result, (646, 1267, 2061)) + assert geometry.shape == (2061, 646, 1267) + assert geometry.spacing_mm == (0.010125, 0.010125, 0.0135) + + +def test_resampled_shape_uses_exact_rationals() -> None: + """794492: (2061, 646, 1267) -> (2087, 654, 1710), verified against the real volume.""" + result = parse_registration_record(_record()) + loaded = loaded_geometry(result, (646, 1267, 2061)) + assert resampled_geometry(loaded).shape == (2087, 654, 1710) + + +def test_the_half_boundary_rounds_up_not_to_even() -> None: + """1270 x 27/20 is exactly 1714.5. + + Python's ``round`` is banker's rounding and gives 1714; binary float division makes + ``0.0135 / 0.01`` slightly less than 1.35, which also gives 1714. The registration + produced 1715 for 791116 and 822177. + """ + result = parse_registration_record(_record()) + loaded = loaded_geometry(result, (657, 1270, 2091)) + assert loaded.shape == (2091, 657, 1270) + assert resampled_geometry(loaded).shape == (2117, 665, 1715) + + +def test_resampled_spacing_is_isotropic() -> None: + """The resampled volume is 10 um on every axis.""" + result = parse_registration_record(_record()) + loaded = loaded_geometry(result, (646, 1267, 2061)) + assert resampled_geometry(loaded).spacing_mm == (0.01, 0.01, 0.01) + + +def test_direction_is_lps_not_the_nifti_srow() -> None: + """ANTs negates the first two rows on read; using the srow mirrors x and y.""" + assert ANTS_DIRECTION == ((0.0, 0.0, -1.0), (1.0, 0.0, 0.0), (0.0, -1.0, 0.0)) + + +def _nifti_header( + shape: tuple[int, int, int], + spacing: tuple[float, float, float], + endian: str = "<", +) -> bytes: + """Build a gzipped NIfTI-1 header carrying a given shape and spacing. + + Parameters + ---------- + shape : tuple[int, int, int] + Array dimensions. + spacing : tuple[float, float, float] + Voxel size. + endian : str, optional + Byte order, by default little. + + Returns + ------- + bytes + The gzipped header. + """ + buffer = bytearray(NIFTI_HEADER_SIZE) + struct.pack_into(f"{endian}i", buffer, 0, NIFTI_HEADER_SIZE) + struct.pack_into(f"{endian}8h", buffer, 40, 3, *shape, 1, 1, 1, 1) + struct.pack_into(f"{endian}8f", buffer, 76, 1.0, *spacing, 0, 0, 0, 0) + return gzip.compress(bytes(buffer)) + + +def test_a_published_header_yields_its_geometry() -> None: + """The header alone gives shape and spacing, so no volume need be fetched.""" + raw = _nifti_header((2061, 646, 1267), (0.010125, 0.010125, 0.0135)) + geometry = parse_nifti_geometry(raw) + assert geometry.shape == (2061, 646, 1267) + assert all(abs(a - b) < 1e-6 for a, b in zip(geometry.spacing_mm, (0.010125, 0.010125, 0.0135))) + + +def test_a_big_endian_header_is_read() -> None: + """Byte order is taken from the header's own size field.""" + raw = _nifti_header((10, 20, 30), (1.0, 2.0, 3.0), endian=">") + assert parse_nifti_geometry(raw).shape == (10, 20, 30) + + +def test_undecompressable_bytes_are_rejected() -> None: + """A non-gzip response is an error, not a silent empty geometry.""" + with pytest.raises(RegistrationRecordError, match="decompress"): + parse_nifti_geometry(b"not gzip at all") + + +def test_a_truncated_header_is_rejected() -> None: + """Too few bytes must not be read as a partial header.""" + with pytest.raises(RegistrationRecordError, match="Decompressed only"): + parse_nifti_geometry(gzip.compress(b"\x5c\x01\x00\x00" + b"\x00" * 40)) + + +def test_a_header_with_too_few_dimensions_is_rejected() -> None: + """A 2D volume cannot supply a 3D geometry.""" + buffer = bytearray(NIFTI_HEADER_SIZE) + struct.pack_into(" None: + """Older samples publish the volumes; newer ones do not.""" + assert registration_volume_key(DATASET, "826509", "loaded") == ( + f"{DATASET}/ccf_alignment/registration_metadata/826509_10um_loaded_zarr_img.nii.gz" + ) + + +def test_reconcile_accepts_a_matching_geometry() -> None: + """Where both exist, the derivation is checked rather than assumed.""" + a = VolumeGeometry((2061, 646, 1267), (0.010125, 0.010125, 0.0135)) + b = VolumeGeometry((2061, 646, 1267), (0.010125000029802322, 0.010125, 0.0135)) + reconcile(a, b, "loaded") + + +def test_reconcile_rejects_a_changed_shape() -> None: + """A shape difference means the registration changed; fail rather than mis-register.""" + a = VolumeGeometry((2061, 646, 1267), (0.010125, 0.010125, 0.0135)) + b = VolumeGeometry((2061, 646, 1268), (0.010125, 0.010125, 0.0135)) + with pytest.raises(RegistrationRecordError, match="shape"): + reconcile(a, b, "loaded") + + +def test_reconcile_rejects_a_changed_spacing() -> None: + """Spacing must agree to better than float32 precision.""" + a = VolumeGeometry((10, 10, 10), (0.01, 0.01, 0.01)) + b = VolumeGeometry((10, 10, 10), (0.01, 0.01, 0.02)) + with pytest.raises(RegistrationRecordError, match="spacing"): + reconcile(a, b, "resampled") From 8e613740e170be36b8a06fe2d3f2a2a21ee8ed05 Mon Sep 17 00:00:00 2001 From: Peter Grotz <199680633+peter-grotz@users.noreply.github.com> Date: Tue, 22 Sep 2026 14:44:14 -0700 Subject: [PATCH 2/2] feat(layout): drop the CCF MouseLight JSON from cell assets The resample stage no longer writes them, so the role has no source. --- src/exaspim_swc_processing/layout.py | 7 ------- tests/test_layout.py | 18 ++++++++---------- tests/test_packaging.py | 1 - 3 files changed, 8 insertions(+), 18 deletions(-) diff --git a/src/exaspim_swc_processing/layout.py b/src/exaspim_swc_processing/layout.py index c89b9d4..6b9760c 100644 --- a/src/exaspim_swc_processing/layout.py +++ b/src/exaspim_swc_processing/layout.py @@ -114,13 +114,6 @@ def _spec( ("final/ccf_space_reconstructions/swcs", "ccf_space_reconstructions/swcs"), "ccf_space_reconstructions", ), - _spec( - "ccf_annotation", - ("final/ccf_space_reconstructions/jsons", "ccf_space_reconstructions/jsons"), - "ccf_space_reconstructions", - suffix=".json", - required=False, - ), ) """Artifacts collected into each cell directory, in the order they are placed.""" diff --git a/tests/test_layout.py b/tests/test_layout.py index 0916045..e9d6320 100644 --- a/tests/test_layout.py +++ b/tests/test_layout.py @@ -23,7 +23,6 @@ "refinement/final-voxel": ".swc", "refinement/final-voxel-resampled": ".swc", "final/ccf_space_reconstructions/swcs": ".swc", - "final/ccf_space_reconstructions/jsons": ".json", } @@ -91,21 +90,21 @@ def test_discover_cells_warns_when_a_source_directory_is_absent( stage_root: Path, caplog: pytest.LogCaptureFixture ) -> None: """A role with no source directory is reported and the remaining roles still index.""" - for path in (stage_root / "final/ccf_space_reconstructions/jsons").iterdir(): + for path in (stage_root / "refinement/final-voxel-resampled").iterdir(): path.unlink() - (stage_root / "final/ccf_space_reconstructions/jsons").rmdir() + (stage_root / "refinement/final-voxel-resampled").rmdir() with caplog.at_level(logging.WARNING): cells = discover_cells(stage_root) - assert "ccf_annotation" in caplog.text + assert "specimen_resampled" in caplog.text for cell in cells.values(): - assert "ccf_annotation" not in cell.artifacts + assert "specimen_resampled" not in cell.artifacts def test_missing_roles_reports_only_required_artifacts(stage_root: Path) -> None: """An absent optional artifact does not make a cell incomplete; a required one does.""" for path in (stage_root / "final/ccf_space_reconstructions/swcs").iterdir(): path.unlink() - for path in (stage_root / "final/ccf_space_reconstructions/jsons").iterdir(): + for path in (stage_root / "refinement/final-voxel-resampled").iterdir(): path.unlink() cells = discover_cells(stage_root) for cell in cells.values(): @@ -176,10 +175,10 @@ def test_place_artifact_replaces_an_existing_destination(tmp_path: Path) -> None def test_write_cell_directory_skips_absent_roles(tmp_path: Path, stage_root: Path) -> None: """Only discovered artifacts are written.""" cell = discover_cells(stage_root)[parse_stem(STEMS[0])] - del cell.artifacts["ccf_annotation"] + del cell.artifacts["specimen_resampled"] written = write_cell_directory(cell, tmp_path / "cell") - assert "ccf_annotation" not in written - assert set(written) == {spec.role for spec in ARTIFACT_SPECS} - {"ccf_annotation"} + assert "specimen_resampled" not in written + assert set(written) == {spec.role for spec in ARTIFACT_SPECS} - {"specimen_resampled"} def test_build_cell_layout_writes_the_expected_tree(stage_root: Path, tmp_path: Path) -> None: @@ -193,7 +192,6 @@ def test_build_cell_layout_writes_the_expected_tree(stage_root: Path, tmp_path: for relative in ( f"specimen_space_reconstructions/refined/{stem}.swc", f"ccf_space_reconstructions/{stem}.swc", - f"ccf_space_reconstructions/{stem}.json", ): assert (cell_dir / relative).is_file(), relative diff --git a/tests/test_packaging.py b/tests/test_packaging.py index e012c75..2277b93 100644 --- a/tests/test_packaging.py +++ b/tests/test_packaging.py @@ -25,7 +25,6 @@ "refinement/final-voxel": ".swc", "refinement/final-voxel-resampled": ".swc", "final/ccf_space_reconstructions/swcs": ".swc", - "final/ccf_space_reconstructions/jsons": ".json", } CREATION_TIME = datetime(2026, 8, 19, 22, 10, 2, tzinfo=timezone.utc) PIPELINE = Code(