From cf2519ef8f90a7c2b9b84a464cbcd520cd0790cc Mon Sep 17 00:00:00 2001 From: "Gabriel A. Devenyi" Date: Sun, 26 Jul 2026 14:38:26 -0400 Subject: [PATCH] Expose the slice step that affine_of_package already computes Fixes #177 affine_of_package builds the centre-to-centre slice step, uses it as the third column of the affine, and discards it. A consumer that needs the number -- a NIfTI writer's third zoom, a DICOM SpacingBetweenSlices -- had to re-derive it, and the two obvious parameter sources are both wrong for a large share of real data. Over the 1599 reconstructions of our corpus that have geometry, VisuCoreSlicePacksSliceDist falling back to VisuCoreFrameThickness disagrees with the affine for 197 and resolution[2] for 82. Both fail structurally, not marginally. VisuCoreFrameThickness is the slab thickness of a 3-D acquisition, so T1_FLASH_3D_iso reports 12.0 mm where the plane step is 0.125 mm. resolution[2] is 0.0 where the third dimension is not spatial (13 parametric maps here) and a cross-package diagonal where the reconstruction holds several packages -- 20*sqrt(2) between two perpendicular packages whose spacings are 2.0 mm (61 here). slice_distance is the length of the affine's slice column, per package, so it inherits the affine's handling of all three cases without special-casing any of them: finite and positive for every one of those 1599, and unavailable for the 27 with no geometry -- the same availability as affine_of_package itself. Spacing is not thickness. They differ for a gapped acquisition, and VisuCoreFrameThickness remains the source for the thickness; a consumer writing BIDS needs both, separately. SlicePackageSplitter, which had the same norm inline, now asks for it. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01NuK1cZi8U54WXAdXMmGpzy --- brukerapi/dataset.py | 25 ++++++++++ brukerapi/splitters.py | 3 +- docs/source/compatibility.rst | 7 +++ test/test_geometry.py | 91 ++++++++++++++++++++++++++++++++++- 4 files changed, 123 insertions(+), 3 deletions(-) diff --git a/brukerapi/dataset.py b/brukerapi/dataset.py index e9db08c..7506461 100644 --- a/brukerapi/dataset.py +++ b/brukerapi/dataset.py @@ -1275,6 +1275,31 @@ def affine(self): ) return self.affine_of_package(0) + @property + def slice_distance(self): + """Centre-to-centre slice spacing of each slice package, in mm. + + This is the length of the affine's slice column, so it reports the step + the geometry is actually built from: the measured distance between two + ``VisuCorePosition`` centres where there are two, the encoded plane step + of a 3-D volume, and the declared ``VisuCoreSlicePacksSliceDist`` only + where nothing was measured. + + Neither parameter a caller would reach for otherwise gives that number. + ``VisuCoreFrameThickness`` is the slab thickness of a 3-D acquisition, + not a plane step, and ``resolution[2]`` is zero where the third + dimension is not spatial and a cross-package diagonal where the + reconstruction holds several packages. + + Spacing is not thickness: a gapped acquisition has slices thinner than + the distance between them, and ``VisuCoreFrameThickness`` remains the + source for the thickness. One value per package, because packages can + be spaced differently (spec 7.10). + + :raise: :UnsupportedDatasetType: if the frames carry no geometry + """ + return [float(np.linalg.norm(self.affine_of_package(package)[:3, 2])) for package in range(len(self.slice_packages_index()))] + def get_slice_packages(self): """Return one in-memory 2dseq dataset per slice package. diff --git a/brukerapi/splitters.py b/brukerapi/splitters.py index 85d66f2..4a25425 100644 --- a/brukerapi/splitters.py +++ b/brukerapi/splitters.py @@ -424,8 +424,7 @@ def _split_VisuCoreSlicePacksSlices(self, visu_pars, frame_count): def _split_VisuCoreSlicePacksSliceDist(self, dataset, visu_pars, sp_index): if "VisuCoreSlicePacksSliceDist" not in visu_pars: - distance = float(np.linalg.norm(dataset.affine_of_package(sp_index)[:3, 2])) - self._synthesise_parameter(visu_pars, "VisuCoreSlicePacksSliceDist", str(distance)) + self._synthesise_parameter(visu_pars, "VisuCoreSlicePacksSliceDist", str(dataset.slice_distance[sp_index])) return VisuCoreSlicePacksSliceDist = visu_pars["VisuCoreSlicePacksSliceDist"] # spec 7.10: the inter-slice distance is a double, and truncating it diff --git a/docs/source/compatibility.rst b/docs/source/compatibility.rst index f3b1f04..169831c 100644 --- a/docs/source/compatibility.rst +++ b/docs/source/compatibility.rst @@ -99,3 +99,10 @@ Data contract and limitations in PV5.1), contiguous frames with a common orientation are treated as one package. Consequently two adjacent inferred packages with the same orientation cannot be distinguished. +* ``Dataset.slice_distance`` is the centre-to-centre slice spacing of each + slice package, in mm -- the length of that same slice column, so it is the + step the geometry is built from. It is not the slice *thickness*: a gapped + acquisition has slices thinner than the distance between them, and + ``VisuCoreFrameThickness`` stays the source for the thickness (for a 3-D + acquisition that parameter is the whole slab, not a plane step). It is + available wherever ``affine_of_package`` is. diff --git a/test/test_geometry.py b/test/test_geometry.py index 97c978b..064627e 100644 --- a/test/test_geometry.py +++ b/test/test_geometry.py @@ -13,7 +13,7 @@ from brukerapi.dataset import LOAD_STAGES, Dataset from brukerapi.exceptions import UnsupportedDatasetType -from test.synthetic import axial_orientation, stacked_positions, write_2dseq +from test.synthetic import axial_orientation, stacked_positions, visu_pars_records, write_2dseq, write_binary, write_jcampdx SAGITTAL_ORIENTATION = np.array([0.0, 1.0, 0.0, 0.0, 0.0, 1.0, -1.0, 0.0, 0.0]) PV360_NIFTI_ROOT = Path("test/test_data/PV360_StdData") @@ -281,3 +281,92 @@ def test_geometry_follows_the_data_when_frames_are_stored_in_reverse(tmp_path): # the data array is the last frame on disk (spec 7.2/7.3). assert np.allclose(dataset.affine[:3, 3], positions[-1]) assert np.allclose(dataset.affine[:3, 2], [0.0, 0.0, -1.5]) + + +def test_slice_distance_is_the_step_the_affine_is_built_from(tmp_path): + """A gapped stack: the spacing is not the slice thickness (spec 7.10). + + `VisuCoreFrameThickness` stays the source for the thickness; nothing + reported the spacing, so a consumer had to re-derive it. + """ + dataset = load( + tmp_path, + frame_groups=(("FG_SLICE", 4),), + positions=stacked_positions((-20.0, -20.0, -7.5), (0.0, 0.0, 5.0), 4), + frame_thickness=2.0, + ) + + assert dataset.slice_distance == [5.0] + assert np.isclose(dataset["VisuCoreFrameThickness"].array[0], 2.0) + + +def test_slice_distance_of_a_3d_volume_is_the_plane_step_not_the_slab(tmp_path): + """`VisuCoreFrameThickness` of a 3-D acquisition is the whole slab. + + Reporting it as a slice step overstates the spacing by the partition count + -- 16 mm instead of 2 mm here. + """ + dataset = load( + tmp_path, + dim=3, + dim_desc=("spatial", "spatial", "spatial"), + size=(4, 4, 8), + extent=(40.0, 40.0, 16.0), + frame_groups=(), + positions=np.array([[-20.0, -20.0, -8.0]]), + orientations=axial_orientation(1), + frame_thickness=16.0, + ) + + assert dataset.slice_distance == [2.0] + + +def test_slice_distance_is_reported_per_package(tmp_path): + """Packages can be spaced differently, and a single number cannot say so. + + `resolution[2]` measures across the whole stack, so for two packages that + do not share an orientation it returns the diagonal between them rather + than either spacing. + """ + positions = np.vstack( + [ + stacked_positions((-20.0, -20.0, -3.0), (0.0, 0.0, 1.5), 5), + stacked_positions((0.0, -20.0, -20.0), (-1.0, 0.0, 0.0), 3), + ] + ) + orientations = np.vstack([axial_orientation(5), np.tile(SAGITTAL_ORIENTATION, (3, 1))]) + dataset = load( + tmp_path, + frame_groups=(("FG_SLICE", 8),), + positions=positions, + orientations=orientations, + slice_packs=(0, [(0, 5), (5, 3)]), + slice_pack_distance=[1.5, 1.0], + ) + + assert dataset.slice_distance == [1.5, 1.0] + + +def test_slice_distance_needs_geometry_and_nothing_else(tmp_path): + """Availability is that of `affine_of_package`: geometry, not spatial frames. + + A parametric map whose third axis is a frame group still has a real slice + spacing, and it is reported; a reconstruction with no VisuCorePosition has + none, and says so. + """ + parametric = load( + tmp_path / "isa", + frame_groups=(("FG_ISA", 2),), + positions=stacked_positions((-20.0, -20.0, 0.0), (0.0, 0.0, 0.0), 2), + frame_thickness=0.8, + slice_packs=(0, [(0, 2)]), + slice_pack_distance=0.8, + ) + assert parametric.slice_distance == [0.8] + + flat = tmp_path / "flat" / "9" / "pdata" / "1" + write_jcampdx(flat / "visu_pars", {key: value for key, value in visu_pars_records().items() if not key.startswith(("VisuCorePosition", "VisuCoreOrientation"))}) + write_binary(flat / "2dseq", np.zeros(4 * 4 * 3), np.dtype("