diff --git a/docs/source/glossary.rst b/docs/source/glossary.rst index 386935e..b966c48 100644 --- a/docs/source/glossary.rst +++ b/docs/source/glossary.rst @@ -38,8 +38,9 @@ than by a thread. **Flat index 0 is reserved.** It doubles as the right-padding value in ``paths``, so ``paths > 0`` is the validity mask everywhere, and - ``project_volume`` overwrites ``volume.flat[0]`` with a sentinel before - a max or min reduction. + ``project_volume`` substitutes a sentinel for the padded entries -- in + the gathered copy, not in the caller's volume -- before a max or min + reduction. surface voxel The first voxel of a :term:`streamline` - the one at the pia. The diff --git a/src/ccf_streamlines/projection.py b/src/ccf_streamlines/projection.py index 00cd965..37873bb 100644 --- a/src/ccf_streamlines/projection.py +++ b/src/ccf_streamlines/projection.py @@ -184,14 +184,17 @@ def _project_volume_to_view(self, volume, kind="max"): if kind == "max": # The path specification assumes the first point in the volume is not a # valid data point and so should be ignored. Since we are doing a - # maximum projection, we set that to the minimum possible value - # so that it won't be selected - volume.flat[0] = min_val - projected_volume.flat[self.view_lookup[:, 0]] = volume.flat[self.paths].max(axis=1) + # maximum projection, we set the padded entries to the minimum possible + # value so that they won't be selected. The substitution is made in the + # gathered copy, so the caller's volume is not modified. + values = volume.flat[self.paths] + values[self.paths == 0] = min_val + projected_volume.flat[self.view_lookup[:, 0]] = values.max(axis=1) elif kind == "min": # Same thing as above, just set to maximum instead of minimum - volume.flat[0] = max_val - projected_volume.flat[self.view_lookup[:, 0]] = volume.flat[self.paths].min(axis=1) + values = volume.flat[self.paths] + values[self.paths == 0] = max_val + projected_volume.flat[self.view_lookup[:, 0]] = values.min(axis=1) elif kind == "mean" or kind == "average": projected_volume.flat[self.view_lookup[:, 0]] = np.nanmean( np.where(self.paths > 0, volume.flat[self.paths], np.nan), @@ -1255,14 +1258,17 @@ def project_volume(self, volume, kind="max"): if kind == "max": # The path specification assumes the first point in the volume is not a # valid data point and so should be ignored. Since we are doing a - # maximum projection, we set that to the minimum possible value - # so that it won't be selected - volume.flat[0] = min_val - values = volume.flat[self.paths].max(axis=1) + # maximum projection, we set the padded entries to the minimum possible + # value so that they won't be selected. The substitution is made in the + # gathered copy, so the caller's volume is not modified. + path_values = volume.flat[self.paths] + path_values[self.paths == 0] = min_val + values = path_values.max(axis=1) elif kind == "min": # Same thing as above, just set to maximum instead of minimum - volume.flat[0] = max_val - values = volume.flat[self.paths].min(axis=1) + path_values = volume.flat[self.paths] + path_values[self.paths == 0] = max_val + values = path_values.min(axis=1) elif kind == "mean" or kind == "average": values = np.nanmean( np.where(self.paths > 0, volume.flat[self.paths], np.nan), diff --git a/tests/mini_ccf.py b/tests/mini_ccf.py index 1204fd2..f0b3f4d 100644 --- a/tests/mini_ccf.py +++ b/tests/mini_ccf.py @@ -205,9 +205,9 @@ class MiniCCF: def volume(self, dtype=np.float64): """A fresh zero volume of the right shape on *every* call. - ``project_volume`` writes a sentinel into ``volume.flat[0]`` and never - restores it, so tests that share one volume leak state into each other. - Always call this rather than caching the result. + ``project_volume`` no longer writes into its input (issue #20), but a + volume is still mutable caller state, so call this rather than caching + the result and letting tests share one. """ return np.zeros(self.volume_shape, dtype=dtype) diff --git a/tests/test_project_volume_mutation.py b/tests/test_project_volume_mutation.py new file mode 100644 index 0000000..ba518c2 --- /dev/null +++ b/tests/test_project_volume_mutation.py @@ -0,0 +1,227 @@ +"""Regression tests for issue #20 -- ``project_volume`` corrupting its input. + +``paths`` is zero-padded on the right and flat index 0 is reserved as "no +voxel", so a max/min reduction has to keep the padded entries from being +selected. It used to do that by writing the dtype's extreme value into +``volume.flat[0]`` and never restoring it, which modified the caller's array: +projecting twice gave different answers, and with ``hemisphere="both"`` the +second pass wrote through the ``np.flip`` *view*, clobbering ``volume[0, 0, -1]`` +as well. + +The substitution now happens in the gathered copy, so these tests assert both +halves of the contract: the caller's volume is untouched, *and* the reserved +voxel is still kept out of the reduction. +""" + +import numpy as np +import pytest + +from ccf_streamlines.projection import ( + Isocortex2dProjector, + IsocortexEntireProjector, +) + + +@pytest.fixture +def projector_factory(mini_ccf): + """Build an ``Isocortex2dProjector`` over the mini-CCF for a hemisphere.""" + + def build(hemisphere="left"): + return Isocortex2dProjector( + mini_ccf.view_lookup_file, + mini_ccf.surface_paths_file, + hemisphere=hemisphere, + ) + + return build + + +@pytest.fixture +def entire_projector(mini_ccf): + return IsocortexEntireProjector(mini_ccf.surface_paths_file) + + +def ramp_volume(mini, dtype=np.float64): + """A volume whose every voxel holds a distinct value, smallest at index 0. + + Starting the ramp at 1 leaves ``volume.flat[0]`` as the unique global + minimum, which is what makes the ``kind="min"`` value tests below able to + see a leak of the reserved voxel. + """ + volume = mini.volume(dtype=dtype) + volume.flat[:] = np.arange(1, volume.size + 1) + return volume + + +def path_extreme(mini, volume, path_index, kind): + """Reduce one streamline's voxels directly, without touching the library.""" + voxels = mini.path_voxels(path_index) + values = volume[voxels[:, 0], voxels[:, 1], voxels[:, 2]] + return values.max() if kind == "max" else values.min() + + +# --------------------------------------------------------------------------- +# The caller's volume must survive a projection +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize("hemisphere", ["left", "right", "both"]) +@pytest.mark.parametrize("kind", ["max", "min"]) +def test_projection_does_not_modify_the_callers_volume( + mini_ccf, projector_factory, hemisphere, kind +): + """No hemisphere branch may write into the input. + + ``"right"`` and ``"both"`` are the interesting ones: they project + ``np.flip(volume, axis=2)``, a view onto the same buffer, so a write at + flat index 0 of the flipped view lands at ``volume[0, 0, -1]``. + """ + projector = projector_factory(hemisphere) + volume = ramp_volume(mini_ccf) + before = volume.copy() + + projector.project_volume(volume, kind=kind) + + # Named separately from the whole-array check so a regression says which + # of the two cells the old code wrote to. + assert volume.flat[0] == before.flat[0] + assert volume[0, 0, -1] == before[0, 0, -1] + assert np.array_equal(volume, before) + + +@pytest.mark.parametrize("dtype", [np.float64, np.float32, np.uint32, np.int16]) +@pytest.mark.parametrize("kind", ["max", "min"]) +def test_projection_leaves_every_dtype_alone(mini_ccf, projector_factory, dtype, kind): + """The sentinel came from ``np.iinfo``/``np.finfo``; cover both branches.""" + projector = projector_factory("both") + volume = ramp_volume(mini_ccf, dtype=dtype) + before = volume.copy() + + projector.project_volume(volume, kind=kind) + + assert np.array_equal(volume, before) + + +@pytest.mark.parametrize("hemisphere", ["left", "right", "both"]) +@pytest.mark.parametrize("kind", ["max", "min"]) +def test_projecting_the_same_volume_twice_gives_the_same_answer( + mini_ccf, projector_factory, hemisphere, kind +): + """Characterization: this passed before the fix too, and says why. + + Repeatability is the invariant a caller actually notices, so it is worth + pinning -- but it is not what detects this defect. The old code wrote the + same sentinel on every call, and neither cell it clobbered, ``(0, 0, 0)`` + or ``(0, 0, z_size - 1)``, lies on a mini-CCF streamline: the fixture puts + every path in the interior. The corruption was therefore invisible to a + repeated projection and only visible in the volume itself, which is what + the tests above assert. + """ + projector = projector_factory(hemisphere) + volume = ramp_volume(mini_ccf) + + first = projector.project_volume(volume, kind=kind) + second = projector.project_volume(volume, kind=kind) + + assert np.array_equal(first, second) + + +@pytest.mark.parametrize("kind", ["max", "min"]) +def test_a_read_only_volume_can_be_projected(mini_ccf, projector_factory, kind): + """A memory-mapped or otherwise read-only reference volume must work. + + Not merely a nicety: writing the sentinel raised ``ValueError: assignment + destination is read-only`` for anything opened read-only. + """ + projector = projector_factory("both") + volume = ramp_volume(mini_ccf) + volume.setflags(write=False) + + projected = projector.project_volume(volume, kind=kind) + + assert projected.shape[0] > 0 + + +@pytest.mark.parametrize("kind", ["max", "min"]) +def test_entire_projection_does_not_modify_the_callers_volume( + mini_ccf, entire_projector, kind +): + volume = ramp_volume(mini_ccf) + before = volume.copy() + + entire_projector.project_volume(volume, kind=kind) + + assert np.array_equal(volume, before) + + +@pytest.mark.parametrize("kind", ["max", "min"]) +def test_entire_projection_is_repeatable(mini_ccf, entire_projector, kind): + """Characterization, for the reason given above.""" + volume = ramp_volume(mini_ccf) + + first = entire_projector.project_volume(volume, kind=kind) + second = entire_projector.project_volume(volume, kind=kind) + + assert np.array_equal(first, second) + + +# --------------------------------------------------------------------------- +# ...and the reserved voxel must still be excluded from the reduction +# --------------------------------------------------------------------------- + + +def test_min_projection_ignores_the_reserved_first_voxel(mini_ccf, projector_factory): + """The padding must not drag every streamline down to ``volume.flat[0]``. + + Every mini-CCF streamline is padded, and the ramp makes index 0 the unique + global minimum, so a projection that reduced over the padding would return + that value for every pixel. + """ + projector = projector_factory("left") + volume = ramp_volume(mini_ccf) + + projected = projector.project_volume(volume, kind="min") + + covered = [] + for path_index in mini_ccf.in_view_path_indices: + row, col = mini_ccf.view_pixel_for_path(int(path_index)) + assert projected[row, col] == path_extreme( + mini_ccf, volume, int(path_index), "min" + ) + covered.append(projected[row, col]) + + # Pixels no streamline reaches keep the zero the view was built with, so + # the claim is about the covered ones. + assert min(covered) > volume.flat[0] + + +def test_max_projection_ignores_the_reserved_first_voxel(mini_ccf, projector_factory): + """The mirror image: index 0 is made the unique global maximum.""" + projector = projector_factory("left") + volume = ramp_volume(mini_ccf) + volume.flat[0] = volume.size + 1 + + projected = projector.project_volume(volume, kind="max") + + for path_index in mini_ccf.in_view_path_indices: + row, col = mini_ccf.view_pixel_for_path(int(path_index)) + assert projected[row, col] == path_extreme( + mini_ccf, volume, int(path_index), "max" + ) + assert projected.max() < volume.flat[0] + + +@pytest.mark.parametrize("kind", ["max", "min"]) +def test_entire_projection_ignores_the_reserved_first_voxel( + mini_ccf, entire_projector, kind +): + """Same contract for the every-streamline projector, whose result is 1D.""" + volume = ramp_volume(mini_ccf) + if kind == "max": + volume.flat[0] = volume.size + 1 + + values = entire_projector.project_volume(volume, kind=kind) + + assert len(values) == mini_ccf.paths.shape[0] + for path_index in range(mini_ccf.paths.shape[0]): + assert values[path_index] == path_extreme(mini_ccf, volume, path_index, kind)