From 3b448b413f8394b209aebfd82ca34e82033a7b60 Mon Sep 17 00:00:00 2001 From: Nathan Gouwens Date: Fri, 28 Aug 2026 17:08:18 -0700 Subject: [PATCH 1/2] Test the pure computational functions Covers every publicly documented name in coordinates, linestring3d, processing, dataset, metrics, morphology and angle, plus the shared voxel-matching helper in projection. These need no fixtures beyond small arrays, so expectations are derived from the domain rather than recorded from the current implementation: a 3-4-5 triangle for segment lengths, a streamline perpendicular to a plane at 90 degrees and one lying in it at 0, one voxel per layer giving 10-micron layers at 10-micron resolution, scale invariance where doubling both inputs must not move the answer. Golden-master values are deliberately absent. One of the two open contributions is a behaviour change whose results match neither architecture's previous answer, so golden values taken from the current branch would fail a correct contribution. Large-allocation helpers are always called with small shapes; their defaults allocate against the full 1320x800x1140 atlas. Two defects are recorded as xfail(strict) tests asserting correct behaviour: - issue #8: load_swc_as_dataframe uses sep=" " and cannot read tab- or column-aligned SWC. It does not fail -- every numeric column comes back NaN, silently. - find_closest_streamline accepts `resolution` and uses it to scale the returned coordinates, but does not forward it to coordinates_to_voxels, so at any resolution other than (10, 10, 10) it voxelises the input against the wrong grid and looks up the wrong streamline. Two behaviours are documented as characterization tests rather than asserted as correct, since neither is knowably a bug: LineString3D.rotation_to_vector divides by 1 + dot and so is undefined for an antiparallel target (cortical streamlines never hit this), and _matching_voxel_indices silently returns wrong answers for an unsorted lookup with no sorter, because np.searchsorted does not check. Co-Authored-By: Claude Opus 5 --- tests/test_angle.py | 218 +++++++++++++++++++++++++++ tests/test_coordinates.py | 74 +++++++++ tests/test_dataset.py | 77 ++++++++++ tests/test_linestring3d.py | 152 +++++++++++++++++++ tests/test_matching_voxel_indices.py | 153 +++++++++++++++++++ tests/test_metrics.py | 142 +++++++++++++++++ tests/test_morphology.py | 204 +++++++++++++++++++++++++ tests/test_processing.py | 86 +++++++++++ 8 files changed, 1106 insertions(+) create mode 100644 tests/test_angle.py create mode 100644 tests/test_dataset.py create mode 100644 tests/test_linestring3d.py create mode 100644 tests/test_matching_voxel_indices.py create mode 100644 tests/test_metrics.py create mode 100644 tests/test_morphology.py create mode 100644 tests/test_processing.py diff --git a/tests/test_angle.py b/tests/test_angle.py new file mode 100644 index 0000000..b9e19e2 --- /dev/null +++ b/tests/test_angle.py @@ -0,0 +1,218 @@ +"""Affine construction, nearest-streamline lookup, and streamline-plane angle. + +The angle cases are anchored at values that can be verified by inspection: +a streamline perpendicular to the plane is 90 degrees, one lying in the plane +is 0, and one at 45 degrees is 45. +""" + +import numpy as np +import pytest + +from ccf_streamlines.angle import ( + determine_angle_between_streamline_and_plane, + find_closest_streamline, + vector_to_3d_affine_matrix, +) + +RESOLUTION_NOT_FORWARDED = ( + "`find_closest_streamline` accepts `resolution` but does not forward it to " + "`coordinates_to_voxels`, so a non-default resolution voxelises against " + "(10, 10, 10); remove this marker when it is fixed" +) + +#: Maps the unit square onto the xy-plane, so the plane normal is +z. +XY_PLANE = vector_to_3d_affine_matrix([1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0]) + + +# -- vector_to_3d_affine_matrix -------------------------------------------- + + +def test_affine_matrix_layout(): + """The first nine entries are the 3x3 basis, the last three the translation.""" + M = vector_to_3d_affine_matrix(list(range(12))) + + assert M.shape == (3, 4) + assert np.array_equal(M[:, :3], np.arange(9).reshape(3, 3)) + assert np.array_equal(M[:, 3], np.array([9, 10, 11])) + + +def test_affine_matrix_translates_the_origin(): + M = vector_to_3d_affine_matrix([1, 0, 0, 0, 1, 0, 0, 0, 1, 5, 6, 7]) + assert np.array_equal(M @ np.array([0, 0, 0, 1]), np.array([5, 6, 7])) + + +# -- determine_angle_between_streamline_and_plane -------------------------- + + +def test_streamline_perpendicular_to_the_plane_is_ninety_degrees(): + """Running along +z, the xy-plane's normal.""" + streamline = np.array([[0.0, 0.0, 10.0], [0.0, 0.0, 0.0]]) + angle = determine_angle_between_streamline_and_plane(streamline, XY_PLANE) + assert angle == pytest.approx(90.0) + + +@pytest.mark.parametrize( + "direction", [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [1.0, 1.0, 0.0]] +) +def test_streamline_lying_in_the_plane_is_zero_degrees(direction): + streamline = np.array([direction, [0.0, 0.0, 0.0]]) + angle = determine_angle_between_streamline_and_plane(streamline, XY_PLANE) + assert angle == pytest.approx(0.0, abs=1e-9) + + +def test_streamline_at_forty_five_degrees(): + streamline = np.array([[1.0, 0.0, 1.0], [0.0, 0.0, 0.0]]) + angle = determine_angle_between_streamline_and_plane(streamline, XY_PLANE) + assert angle == pytest.approx(45.0) + + +def test_only_the_endpoints_of_the_streamline_matter(): + """The function uses the pia and white-matter ends, not the path between.""" + straight = np.array([[0.0, 0.0, 10.0], [0.0, 0.0, 0.0]]) + wiggly = np.array([[0.0, 0.0, 10.0], [5.0, 5.0, 5.0], [0.0, 0.0, 0.0]]) + assert determine_angle_between_streamline_and_plane( + straight, XY_PLANE + ) == pytest.approx(determine_angle_between_streamline_and_plane(wiggly, XY_PLANE)) + + +def test_a_rotated_plane_rotates_the_angle(): + """Swap the plane to the xz-plane; a streamline along +z now lies in it.""" + xz_plane = vector_to_3d_affine_matrix([1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0]) + streamline = np.array([[0.0, 0.0, 10.0], [0.0, 0.0, 0.0]]) + assert determine_angle_between_streamline_and_plane( + streamline, xz_plane + ) == pytest.approx(0.0, abs=1e-9) + + +def test_plane_translation_does_not_change_the_angle(): + translated = vector_to_3d_affine_matrix([1, 0, 0, 0, 1, 0, 0, 0, 1, 100, 200, 300]) + streamline = np.array([[0.0, 0.0, 10.0], [0.0, 0.0, 0.0]]) + assert determine_angle_between_streamline_and_plane( + streamline, translated + ) == pytest.approx(90.0) + + +# -- find_closest_streamline ----------------------------------------------- + + +def test_a_coordinate_on_a_streamline_returns_that_streamline(mini_ccf): + path_index = 0 + coord = mini_ccf.coord_on_path(path_index, 3) + + result = find_closest_streamline( + coord, + mini_ccf.closest_surface_voxel_file, + mini_ccf.surface_paths_file, + resolution=mini_ccf.resolution, + volume_shape=mini_ccf.volume_shape, + ) + + assert np.array_equal(result, mini_ccf.path_microns(path_index)) + + +def test_the_reference_may_be_a_preloaded_array(mini_ccf): + """The documented alternative to passing a file path.""" + import h5py + + with h5py.File(mini_ccf.closest_surface_voxel_file, "r") as f: + closest = f["closest surface voxel"][:] + + coord = mini_ccf.coord_on_path(0, 3) + from_array = find_closest_streamline( + coord, closest, mini_ccf.surface_paths_file, + resolution=mini_ccf.resolution, volume_shape=mini_ccf.volume_shape, + ) + from_path = find_closest_streamline( + coord, mini_ccf.closest_surface_voxel_file, mini_ccf.surface_paths_file, + resolution=mini_ccf.resolution, volume_shape=mini_ccf.volume_shape, + ) + assert np.array_equal(from_array, from_path) + + +def test_surface_paths_may_be_an_open_h5py_file(mini_ccf): + import h5py + + coord = mini_ccf.coord_on_path(0, 3) + with h5py.File(mini_ccf.surface_paths_file, "r") as f: + from_handle = find_closest_streamline( + coord, mini_ccf.closest_surface_voxel_file, f, + resolution=mini_ccf.resolution, volume_shape=mini_ccf.volume_shape, + ) + assert np.array_equal(from_handle, mini_ccf.path_microns(0)) + + +def test_a_coordinate_outside_cortex_returns_an_empty_array(mini_ccf, caplog): + """The dorso-ventral planes at each end have no streamline voxels.""" + outside = np.array([3.0, 0.0, 1.0]) * np.array(mini_ccf.resolution) + + result = find_closest_streamline( + outside, + mini_ccf.closest_surface_voxel_file, + mini_ccf.surface_paths_file, + resolution=mini_ccf.resolution, + volume_shape=mini_ccf.volume_shape, + ) + + assert result.size == 0 + assert "not within isocortex" in caplog.text + + +def test_a_right_hemisphere_coordinate_comes_back_on_the_right(mini_ccf): + """Reference data exists only on the left, so the lookup reflects, then + reflects the answer back.""" + z_size = mini_ccf.volume_shape[2] + left_voxel = mini_ccf.path_voxels(0)[3] + right_voxel = left_voxel.copy() + right_voxel[2] = z_size - left_voxel[2] + coord = right_voxel * np.array(mini_ccf.resolution) + + result = find_closest_streamline( + coord, + mini_ccf.closest_surface_voxel_file, + mini_ccf.surface_paths_file, + resolution=mini_ccf.resolution, + volume_shape=mini_ccf.volume_shape, + ) + + expected = mini_ccf.path_voxels(0).copy() + expected[:, 2] = z_size - expected[:, 2] + assert np.array_equal(result, expected * np.array(mini_ccf.resolution)) + + +def test_a_coordinate_may_be_given_as_a_flat_triple(mini_ccf): + coord = mini_ccf.coord_on_path(0, 3) + flat = find_closest_streamline( + coord, mini_ccf.closest_surface_voxel_file, mini_ccf.surface_paths_file, + resolution=mini_ccf.resolution, volume_shape=mini_ccf.volume_shape, + ) + nested = find_closest_streamline( + coord.reshape(1, 3), mini_ccf.closest_surface_voxel_file, + mini_ccf.surface_paths_file, + resolution=mini_ccf.resolution, volume_shape=mini_ccf.volume_shape, + ) + assert np.array_equal(flat, nested) + + +@pytest.mark.xfail(strict=True, reason=RESOLUTION_NOT_FORWARDED) +def test_resolution_is_honoured_when_finding_the_streamline(mini_ccf): + """The same physical point, expressed at a coarser voxel size. + + ``resolution`` is used to scale the *returned* coordinates but not to + convert the *input* coordinate to a voxel, so at any resolution other than + (10, 10, 10) the wrong voxel is looked up. Here it lands outside the + lookup entirely and an empty array comes back. + """ + resolution = (20, 20, 20) + voxel = mini_ccf.path_voxels(0)[3] + coord = voxel * np.array(resolution) + + result = find_closest_streamline( + coord, + mini_ccf.closest_surface_voxel_file, + mini_ccf.surface_paths_file, + resolution=resolution, + volume_shape=mini_ccf.volume_shape, + ) + + assert result.size > 0 + assert np.array_equal(result, mini_ccf.path_voxels(0) * np.array(resolution)) diff --git a/tests/test_coordinates.py b/tests/test_coordinates.py index 99b969f..cbaace5 100644 --- a/tests/test_coordinates.py +++ b/tests/test_coordinates.py @@ -48,3 +48,77 @@ def test_coords_to_voxels(): double_resolution ) == expected_voxels ) + + +def test_default_resolution_is_ten_microns(): + test_coords = np.array([[0., 0., 0.], [15., 25., 35.]]) + + assert np.all( + coordinates.coordinates_to_voxels(test_coords) == + coordinates.coordinates_to_voxels(test_coords, (10, 10, 10)) + ) + + +def test_anisotropic_resolution_is_applied_per_axis(): + test_coords = np.array([[100., 100., 100.]]) + expected_voxels = np.array([[10, 5, 1]]) + + assert np.all( + coordinates.coordinates_to_voxels(test_coords, (10, 20, 100)) == + expected_voxels + ) + + +def test_coordinates_are_floored_not_rounded(): + # 19.9 microns is still inside the second 10-micron voxel + test_coords = np.array([[9.9, 10.0, 19.9]]) + expected_voxels = np.array([[0, 1, 1]]) + + assert np.all( + coordinates.coordinates_to_voxels(test_coords, (10, 10, 10)) == + expected_voxels + ) + + +def test_negative_coordinates_floor_away_from_zero(): + # Flooring is toward negative infinity, so -0.1 microns is voxel -1, not 0. + test_coords = np.array([[-0.1, -10., -25.]]) + expected_voxels = np.array([[-1, -1, -3]]) + + assert np.all( + coordinates.coordinates_to_voxels(test_coords, (10, 10, 10)) == + expected_voxels + ) + + +def test_result_is_an_integer_array(): + test_coords = np.array([[0., 0., 0.], [15., 25., 35.]]) + voxels = coordinates.coordinates_to_voxels(test_coords, (10, 10, 10)) + + assert np.issubdtype(voxels.dtype, np.integer) + assert voxels.shape == test_coords.shape + + +def test_non_numeric_dtype(): + test_coords = np.array([ + ["0", "0", "0"], + ["1", "1", "1"], + ]) + resolution = (10, 10, 10) + + with pytest.raises(ValueError, match="numeric dtype"): + coordinates.coordinates_to_voxels( + test_coords, + resolution) + + +def test_two_dimensional_coordinates_work_with_a_two_tuple(): + # The function is not hardcoded to three dimensions; it only requires that + # `resolution` match the second dimension of `coords`. + test_coords = np.array([[0., 0.], [15., 25.]]) + expected_voxels = np.array([[0, 0], [1, 2]]) + + assert np.all( + coordinates.coordinates_to_voxels(test_coords, (10, 10)) == + expected_voxels + ) diff --git a/tests/test_dataset.py b/tests/test_dataset.py new file mode 100644 index 0000000..655ab97 --- /dev/null +++ b/tests/test_dataset.py @@ -0,0 +1,77 @@ +"""The ISH upsampling helper. + +Always called with a small ``target_volume_shape``: the default allocates a +1320 x 800 x 1140 float array, several gigabytes, which no test should do. +""" + +import numpy as np +import pytest + +from ccf_streamlines.dataset import upscale_ish_volume + + +@pytest.fixture +def small_volume(): + """Shape (3, 1, 2) with distinct values, so any axis mixup is visible.""" + return np.arange(6, dtype=float).reshape(3, 1, 2) + + +def test_each_target_voxel_takes_its_downscaled_source_value(small_volume): + """With a 2x ratio, target voxel (i, j, k) comes from source (i//2, j//2, k//2).""" + # rotate_axes swaps 0 and 2, so a (3, 1, 2) input becomes (2, 1, 3). + result = upscale_ish_volume( + small_volume, + orig_voxel_size=20, + target_voxel_size=10, + target_volume_shape=(4, 2, 6), + ) + swapped = np.swapaxes(small_volume, 0, 2) + + assert result.shape == (4, 2, 6) + for i in range(4): + for j in range(2): + for k in range(6): + assert result[i, j, k] == swapped[i // 2, j // 2, k // 2] + + +def test_rotate_axes_swaps_the_first_and_last_axes(small_volume): + """Volumes from the ISH atlas API have anterior-posterior in z and + left-right in x; the CCF has those swapped.""" + rotated = upscale_ish_volume( + small_volume, orig_voxel_size=10, target_voxel_size=10, + target_volume_shape=(2, 1, 3), rotate_axes=True, + ) + unrotated = upscale_ish_volume( + np.swapaxes(small_volume, 0, 2), orig_voxel_size=10, target_voxel_size=10, + target_volume_shape=(2, 1, 3), rotate_axes=False, + ) + assert np.array_equal(rotated, unrotated) + + +def test_without_rotation_a_matching_shape_round_trips(small_volume): + """A 1:1 ratio and a matching target shape is the identity.""" + result = upscale_ish_volume( + small_volume, orig_voxel_size=10, target_voxel_size=10, + target_volume_shape=small_volume.shape, rotate_axes=False, + ) + assert np.array_equal(result, small_volume) + + +def test_upscaling_repeats_each_source_voxel_ratio_times(): + """A single source voxel fills a ratio-cubed block of the target.""" + volume = np.array([[[7.0]]]) + result = upscale_ish_volume( + volume, orig_voxel_size=30, target_voxel_size=10, + target_volume_shape=(3, 3, 3), rotate_axes=False, + ) + assert np.array_equal(result, np.full((3, 3, 3), 7.0)) + + +def test_target_larger_than_the_scaled_source_raises(small_volume): + """Asking for more target voxels than the source can cover is an + out-of-bounds index, not a silent zero-fill.""" + with pytest.raises(IndexError): + upscale_ish_volume( + small_volume, orig_voxel_size=20, target_voxel_size=10, + target_volume_shape=(100, 2, 6), + ) diff --git a/tests/test_linestring3d.py b/tests/test_linestring3d.py new file mode 100644 index 0000000..1794ce2 --- /dev/null +++ b/tests/test_linestring3d.py @@ -0,0 +1,152 @@ +"""Streamline geometry, against analytically known cases. + +Every depth calculation in the package is built on this class, so the cases +here are chosen so the right answer can be worked out by hand: a 3-4-5 +triangle, an axis-aligned path, a point whose perpendicular offset is exact. +""" + +import numpy as np +import pytest + +from ccf_streamlines.linestring3d import LineString3D + + +@pytest.fixture +def straight_path(): + """Ten units long, along +y, so distances along it are just y.""" + return LineString3D(np.array([[0.0, y, 0.0] for y in range(11)])) + + +@pytest.fixture +def bent_path(): + """Two segments: 3 along x, then 4 along y. Total length 7.""" + return LineString3D(np.array([[0.0, 0.0, 0.0], [3.0, 0.0, 0.0], [3.0, 4.0, 0.0]])) + + +def test_segment_lengths(bent_path): + assert np.array_equal(bent_path.segment_lengths(), np.array([3.0, 4.0])) + + +def test_length_is_the_sum_of_segments(bent_path, straight_path): + assert bent_path.length == 7.0 + assert straight_path.length == 10.0 + + +def test_segment_lengths_of_a_diagonal(): + """A 3-4-5 triangle, so the hypotenuse is exactly 5.""" + path = LineString3D(np.array([[0.0, 0.0, 0.0], [3.0, 4.0, 0.0]])) + assert path.segment_lengths() == pytest.approx([5.0]) + assert path.length == pytest.approx(5.0) + + +# -- project --------------------------------------------------------------- + + +@pytest.mark.parametrize("y", [0.0, 1.0, 5.5, 10.0]) +def test_project_a_point_lying_on_the_path(straight_path, y): + assert straight_path.project(np.array([0.0, y, 0.0])) == pytest.approx(y) + + +def test_project_a_point_beside_the_path(straight_path): + """Distance along is unaffected by perpendicular displacement.""" + assert straight_path.project(np.array([3.0, 4.0, 2.0])) == pytest.approx(4.0) + + +def test_project_normalized_is_the_fraction_of_total_length(straight_path): + assert straight_path.project(np.array([0.0, 2.5, 0.0]), normalized=True) == pytest.approx(0.25) + assert straight_path.project(np.array([0.0, 10.0, 0.0]), normalized=True) == pytest.approx(1.0) + + +def test_project_across_a_bend(bent_path): + """Three units along the first segment, then two into the second.""" + assert bent_path.project(np.array([3.0, 2.0, 0.0])) == pytest.approx(5.0) + + +def test_project_before_the_start_clamps_to_zero(straight_path): + """A point behind the pia end does not project onto any segment.""" + assert straight_path.project(np.array([0.0, -5.0, 0.0])) == pytest.approx(0.0) + + +def test_project_past_the_end_clamps_to_the_last_vertex(straight_path): + """Past the white-matter end, the nearest vertex is the last one. + + Note the returned value is the length up to that vertex, which for the + final vertex is the whole path. + """ + assert straight_path.project(np.array([0.0, 50.0, 0.0])) == pytest.approx(10.0) + + +# -- offset_of_point ------------------------------------------------------- + + +def test_offset_of_a_point_on_the_path_is_zero(straight_path): + """The property that makes exact integer assertions possible elsewhere.""" + offset = straight_path.offset_of_point(np.array([0.0, 4.0, 0.0])) + assert np.array_equal(offset, np.zeros(3)) + + +def test_offset_is_the_perpendicular_displacement(straight_path): + offset = straight_path.offset_of_point(np.array([2.0, 4.0, -3.0])) + assert offset == pytest.approx([2.0, 0.0, -3.0]) + + +def test_offset_of_a_point_beyond_the_end_is_measured_from_the_last_vertex(straight_path): + offset = straight_path.offset_of_point(np.array([1.0, 20.0, 0.0])) + assert offset == pytest.approx([1.0, 10.0, 0.0]) + + +def test_offset_across_a_bend(bent_path): + offset = bent_path.offset_of_point(np.array([3.0, 2.0, 5.0])) + assert offset == pytest.approx([0.0, 0.0, 5.0]) + + +# -- rotation_to_vector ---------------------------------------------------- + + +def test_rotation_aligns_the_path_with_the_target_vector(): + """A path along +x, rotated to point along +y.""" + path = LineString3D(np.array([[0.0, 0.0, 0.0], [5.0, 0.0, 0.0]])) + rot = path.rotation_to_vector(np.array([0.0, 1.0, 0.0])) + + rotated_end = rot @ path.coords[-1, :] + assert rotated_end == pytest.approx([0.0, 5.0, 0.0]) + + +def test_rotation_is_orthonormal_and_preserves_length(): + path = LineString3D(np.array([[0.0, 0.0, 0.0], [1.0, 2.0, 3.0]])) + rot = path.rotation_to_vector(np.array([0.0, 1.0, 0.0])) + + assert rot @ rot.T == pytest.approx(np.identity(3)) + assert np.linalg.det(rot) == pytest.approx(1.0) + assert np.linalg.norm(rot @ path.coords[-1, :]) == pytest.approx(np.sqrt(14)) + + +def test_rotation_to_an_already_aligned_vector_is_the_identity(): + """The mini-CCF's streamlines run straight down +y, so this is the case + the coordinate projector actually hits, and it must be exact.""" + path = LineString3D(np.array([[0.0, 0.0, 0.0], [0.0, 7.0, 0.0]])) + rot = path.rotation_to_vector(np.array([0.0, 1.0, 0.0])) + assert np.array_equal(rot, np.identity(3)) + + +def test_rotation_target_need_not_be_unit_length(): + path = LineString3D(np.array([[0.0, 0.0, 0.0], [5.0, 0.0, 0.0]])) + assert np.array_equal( + path.rotation_to_vector(np.array([0.0, 1.0, 0.0])), + path.rotation_to_vector(np.array([0.0, 9.0, 0.0])), + ) + + +def test_rotation_to_an_antiparallel_vector_is_undefined(): + """Characterization, not an assertion that this is correct. + + The Rodrigues construction divides by ``1 + dot``, which is zero when the + path points exactly opposite the target. Cortical streamlines always run + roughly pia-to-white-matter, so the coordinate projector never asks for + this -- but a caller using `LineString3D` directly can, and gets + non-finite values rather than an error. + """ + path = LineString3D(np.array([[0.0, 0.0, 0.0], [0.0, -1.0, 0.0]])) + with np.errstate(divide="ignore", invalid="ignore"): + rot = path.rotation_to_vector(np.array([0.0, 1.0, 0.0])) + assert not np.all(np.isfinite(rot)) diff --git a/tests/test_matching_voxel_indices.py b/tests/test_matching_voxel_indices.py new file mode 100644 index 0000000..28c140b --- /dev/null +++ b/tests/test_matching_voxel_indices.py @@ -0,0 +1,153 @@ +"""The shared voxel-matching helper. + +This is the common core beneath the closest-streamline search, surface-voxel +collapsing, and 2-D coordinate calculation, and both open contributions reduce +to its behaviour. It is pure and takes small arrays, so it can be covered +exhaustively. + +Its contract, stated plainly: given a lookup whose ``lookup_ind`` column is +sorted (or made sorted by a ``sorter``), return the ``ref_ind`` column of the +matching row for each query, and ``missing_value`` for queries with no match. +Ties in the key column resolve to whichever tied row the ordering placed first. +""" + +import numpy as np +import pytest + +from ccf_streamlines.projection import _matching_voxel_indices + + +@pytest.fixture +def lookup(): + """Key column sorted and unique; value column deliberately unsorted.""" + return np.array([ + [10, 700], + [20, 500], + [30, 900], + [40, 100], + ]) + + +def test_every_query_returns_its_matching_row(lookup): + result = _matching_voxel_indices(np.array([10, 30, 40]), lookup) + assert np.array_equal(result, np.array([700, 900, 100])) + + +def test_query_order_is_preserved(lookup): + result = _matching_voxel_indices(np.array([40, 10, 30, 20]), lookup) + assert np.array_equal(result, np.array([100, 700, 900, 500])) + + +def test_a_repeated_query_returns_the_same_answer_each_time(lookup): + result = _matching_voxel_indices(np.array([20, 20, 20]), lookup) + assert np.array_equal(result, np.array([500, 500, 500])) + + +def test_missing_keys_get_the_missing_value(lookup): + """Zero by default, which is why flat index 0 must never be a real voxel.""" + result = _matching_voxel_indices(np.array([10, 15, 40]), lookup) + assert np.array_equal(result, np.array([700, 0, 100])) + + +def test_the_missing_value_is_configurable(lookup): + result = _matching_voxel_indices(np.array([15]), lookup, missing_value=-1) + assert result[0] == -1 + + +def test_all_missing(lookup): + result = _matching_voxel_indices(np.array([1, 2, 3]), lookup, missing_value=-1) + assert np.array_equal(result, np.array([-1, -1, -1])) + + +def test_an_empty_query_returns_an_empty_result(lookup): + result = _matching_voxel_indices(np.array([], dtype=int), lookup) + assert result.shape == (0,) + + +def test_columns_can_be_swapped(lookup): + """The 2-D coordinate path searches column 1 and returns column 0.""" + result = _matching_voxel_indices( + np.array([900, 100]), lookup, lookup_ind=1, ref_ind=0, missing_value=-1, + sorter=np.argsort(lookup[:, 1]), + ) + assert np.array_equal(result, np.array([30, 40])) + + +def test_a_sorter_makes_an_unsorted_key_column_searchable(): + """The view-lookup files are not sorted on their volume-index column.""" + unsorted = np.array([ + [1, 300], + [2, 100], + [3, 400], + [4, 200], + ]) + sorter = np.argsort(unsorted[:, 1]) + + result = _matching_voxel_indices( + np.array([100, 400]), unsorted, lookup_ind=1, ref_ind=0, + missing_value=-1, sorter=sorter, + ) + assert np.array_equal(result, np.array([2, 3])) + + +def test_without_a_sorter_an_unsorted_lookup_returns_wrong_answers_silently(): + """Characterization, not an endorsement. + + ``np.searchsorted`` assumes a sorted array and does not check. Nothing in + the package validates the reference file's ordering, so an unsorted lookup + yields wrong voxels rather than an error. Documented here so the assumption + is visible. + """ + unsorted = np.array([ + [1, 300], + [2, 100], + [3, 400], + [4, 200], + ]) + + without_sorter = _matching_voxel_indices( + np.array([100]), unsorted, lookup_ind=1, ref_ind=0, missing_value=-1 + ) + with_sorter = _matching_voxel_indices( + np.array([100]), unsorted, lookup_ind=1, ref_ind=0, missing_value=-1, + sorter=np.argsort(unsorted[:, 1]), + ) + + assert with_sorter[0] == 2 + assert without_sorter[0] != with_sorter[0] + + +def test_a_tie_resolves_to_the_first_row_in_sorter_order(): + """The behaviour issue #12 is about, in its smallest possible form. + + ``np.searchsorted`` returns the *left* insertion point, so a tied key picks + whichever tied row the ordering placed first. The helper is not wrong; the + caller choosing an unstable ordering is. + """ + tied = np.array([ + [7, 100], + [8, 100], + [9, 200], + ]) + + stable = _matching_voxel_indices( + np.array([100]), tied, lookup_ind=1, ref_ind=0, missing_value=-1, + sorter=np.argsort(tied[:, 1], kind="stable"), + ) + reversed_ties = _matching_voxel_indices( + np.array([100]), tied, lookup_ind=1, ref_ind=0, missing_value=-1, + sorter=np.lexsort((-np.arange(len(tied)), tied[:, 1])), + ) + + assert stable[0] == 7 + assert reversed_ties[0] == 8 + + +def test_the_result_dtype_is_integer(lookup): + result = _matching_voxel_indices(np.array([10]), lookup) + assert np.issubdtype(result.dtype, np.integer) + + +def test_the_result_has_the_shape_of_the_query(lookup): + query = np.array([10, 20, 30, 40, 50]) + assert _matching_voxel_indices(query, lookup).shape == query.shape diff --git a/tests/test_metrics.py b/tests/test_metrics.py new file mode 100644 index 0000000..5d8d3bb --- /dev/null +++ b/tests/test_metrics.py @@ -0,0 +1,142 @@ +"""Layer thickness measurement, against streamlines whose layers we chose. + +The layer memberships are assigned voxel by voxel here, so the returned start, +end and thickness are all derivable by hand: a straight path of unit-spaced +voxels at 10 micron resolution puts layer *i* between 10*i and 10*(i+1). +""" + +import numpy as np +import pytest + +from ccf_streamlines.metrics import measure_streamline_layer_thicknesses + +# Structure set IDs from the mouse ontology, in cortical depth order. These are +# also the values `measure_streamline_layer_thicknesses` looks for, and their +# numeric order matches their depth order -- which the function relies on when +# it sorts annotations to stop layers being intercalated. +LAYER_IDS = [667481440, 667481441, 667481445, 667481446, 667481449, 667481450] +LAYER_NAMES = [ + "Isocortex layer 1", + "Isocortex layer 2/3", + "Isocortex layer 4", + "Isocortex layer 5", + "Isocortex layer 6a", + "Isocortex layer 6b", +] +RESOLUTION = (10, 10, 10) +SHAPE = (4, 10, 3) + + +def _straight_path(x, z, layer_ids, padded_length=8): + """A path down +y at (x, z), one voxel per entry of ``layer_ids``. + + Returns the padded flat-index row and a volume labelled to match. + """ + volume = np.zeros(SHAPE, dtype=np.int64) + ys = np.arange(1, 1 + len(layer_ids)) + for y, layer in zip(ys, layer_ids): + volume[x, y, z] = layer + flat = np.ravel_multi_index((np.full_like(ys, x), ys, np.full_like(ys, z)), SHAPE) + row = np.zeros(padded_length, dtype=np.int64) + row[: len(flat)] = flat + return row, volume + + +def test_one_voxel_per_layer_gives_ten_micron_layers(): + """Six voxels, one per layer, so every layer is exactly one voxel thick.""" + row, volume = _straight_path(1, 1, LAYER_IDS) + paths = row.reshape(1, -1) + + result = measure_streamline_layer_thicknesses(volume, paths, RESOLUTION) + + assert set(result) == set(LAYER_NAMES) + for i, name in enumerate(LAYER_NAMES): + start, end, thickness = result[name][0] + assert start == pytest.approx(10.0 * i) + assert end == pytest.approx(10.0 * (i + 1)) + assert thickness == pytest.approx(10.0) + + +def test_thicknesses_follow_the_number_of_voxels_in_each_layer(): + """Layer 2/3 given three voxels must come back three times as thick.""" + layers = [ + LAYER_IDS[0], + LAYER_IDS[1], LAYER_IDS[1], LAYER_IDS[1], + LAYER_IDS[2], + LAYER_IDS[3], + LAYER_IDS[4], + LAYER_IDS[5], + ] + row, volume = _straight_path(1, 1, layers, padded_length=10) + + result = measure_streamline_layer_thicknesses(volume, row.reshape(1, -1), RESOLUTION) + + assert result["Isocortex layer 1"][0][2] == pytest.approx(10.0) + assert result["Isocortex layer 2/3"][0][2] == pytest.approx(30.0) + assert result["Isocortex layer 2/3"][0][0] == pytest.approx(10.0) + assert result["Isocortex layer 2/3"][0][1] == pytest.approx(40.0) + assert result["Isocortex layer 4"][0][0] == pytest.approx(40.0) + + +def test_an_absent_layer_reports_all_zeros(): + """This is how the projectors detect that a layer is missing.""" + layers = [i for i in LAYER_IDS if i != LAYER_IDS[2]] # no layer 4 + row, volume = _straight_path(1, 1, layers) + + result = measure_streamline_layer_thicknesses(volume, row.reshape(1, -1), RESOLUTION) + + assert np.array_equal(result["Isocortex layer 4"][0], np.zeros(3)) + # ...and the layers around it stay contiguous + assert result["Isocortex layer 2/3"][0][1] == pytest.approx(20.0) + assert result["Isocortex layer 5"][0][0] == pytest.approx(20.0) + + +def test_total_thickness_equals_the_number_of_annotated_voxels(): + row, volume = _straight_path(1, 1, LAYER_IDS) + result = measure_streamline_layer_thicknesses(volume, row.reshape(1, -1), RESOLUTION) + + total = sum(result[name][0][2] for name in LAYER_NAMES) + assert total == pytest.approx(60.0) + + +def test_several_streamlines_are_measured_independently(): + layers_a = LAYER_IDS + layers_b = [LAYER_IDS[0], LAYER_IDS[0], LAYER_IDS[1], LAYER_IDS[2], + LAYER_IDS[3], LAYER_IDS[4], LAYER_IDS[5]] + row_a, volume = _straight_path(1, 1, layers_a, padded_length=10) + row_b, volume_b = _straight_path(2, 2, layers_b, padded_length=10) + volume = volume + volume_b + + paths = np.vstack([row_a, row_b]) + result = measure_streamline_layer_thicknesses(volume, paths, RESOLUTION) + + assert result["Isocortex layer 1"][0][2] == pytest.approx(10.0) + assert result["Isocortex layer 1"][1][2] == pytest.approx(20.0) + + +def test_duplicate_consecutive_voxels_are_collapsed_before_measuring(): + """A path that lingers on a voxel must not double-count its layer.""" + row, volume = _straight_path(1, 1, LAYER_IDS, padded_length=10) + with_dupes = np.zeros(12, dtype=np.int64) + values = row[row > 0] + doubled = np.repeat(values[:1], 2) # visit the first voxel twice + packed = np.concatenate([doubled, values[1:]]) + with_dupes[: len(packed)] = packed + + plain = measure_streamline_layer_thicknesses(volume, row.reshape(1, -1), RESOLUTION) + duped = measure_streamline_layer_thicknesses( + volume, with_dupes.reshape(1, -1), RESOLUTION + ) + + for name in LAYER_NAMES: + assert plain[name][0] == pytest.approx(duped[name][0]) + + +def test_resolution_scales_the_measured_thicknesses(): + row, volume = _straight_path(1, 1, LAYER_IDS) + + at_10 = measure_streamline_layer_thicknesses(volume, row.reshape(1, -1), (10, 10, 10)) + at_20 = measure_streamline_layer_thicknesses(volume, row.reshape(1, -1), (20, 20, 20)) + + for name in LAYER_NAMES: + assert at_20[name][0][2] == pytest.approx(2 * at_10[name][0][2]) diff --git a/tests/test_morphology.py b/tests/test_morphology.py new file mode 100644 index 0000000..5fce01a --- /dev/null +++ b/tests/test_morphology.py @@ -0,0 +1,204 @@ +"""Morphology loading and coordinate-to-volume transformation. + +``transform_swc_to_volume`` and ``transform_coordinates_to_volume`` default to +the full CCF shape, which allocates a 1320 x 800 x 1140 array. Every call here +passes a small shape explicitly. +""" + +import numpy as np +import pandas as pd +import pytest + +from ccf_streamlines.morphology import ( + find_topological_point_coordinates, + load_swc_as_dataframe, + transform_coordinates_to_volume, + transform_swc_to_volume, +) + +ISSUE_8 = ( + "pinned to AllenInstitute/ccf_streamlines#8: load_swc_as_dataframe uses " + "sep=' ' and cannot read the delimiters real reconstruction tools emit; " + "remove this marker when it is fixed" +) + +# id type x y z r parent_id +SPACE_DELIMITED = ( + "# generated by a reconstruction tool\n" + "1 1 0.0 0.0 0.0 1.0 -1\n" + "2 3 10.0 0.0 0.0 1.0 1\n" + "3 3 20.0 0.0 0.0 1.0 2\n" + "4 3 10.0 10.0 0.0 1.0 2\n" +) + + +@pytest.fixture +def swc_file(tmp_path): + path = tmp_path / "neuron.swc" + path.write_text(SPACE_DELIMITED) + return str(path) + + +def test_space_delimited_swc_loads(swc_file): + df = load_swc_as_dataframe(swc_file) + + assert list(df.columns) == ["id", "type", "x", "y", "z", "r", "parent_id"] + assert len(df) == 4 + assert df["id"].tolist() == [1, 2, 3, 4] + assert df["x"].tolist() == [0.0, 10.0, 20.0, 10.0] + assert df["parent_id"].tolist() == [-1, 1, 2, 2] + + +def test_comment_lines_are_skipped(swc_file): + df = load_swc_as_dataframe(swc_file) + assert not df["id"].astype(str).str.startswith("#").any() + + +@pytest.mark.xfail(strict=True, reason=ISSUE_8) +def test_tab_delimited_swc_loads(tmp_path): + """MorphIO v3.3.3 and other tools emit tab-separated SWC. + + With ``sep=" "`` the whole line lands in the first column and every + numeric column is NaN -- silently, so the caller gets a dataframe of the + right shape full of nothing. + """ + path = tmp_path / "tabs.swc" + path.write_text(SPACE_DELIMITED.replace(" ", "\t").replace("#\tgenerated", "# generated")) + + df = load_swc_as_dataframe(str(path)) + + assert df["x"].tolist() == [0.0, 10.0, 20.0, 10.0] + assert not df["x"].isna().any() + + +@pytest.mark.xfail(strict=True, reason=ISSUE_8) +def test_column_aligned_swc_loads(tmp_path): + """Runs of spaces used for alignment are equally common, and equally fatal.""" + path = tmp_path / "aligned.swc" + path.write_text( + "1 1 0.0 0.0 0.0 1.0 -1\n" + "2 3 10.0 0.0 0.0 1.0 1\n" + ) + + df = load_swc_as_dataframe(str(path)) + + assert df["x"].tolist() == [0.0, 10.0] + assert not df["x"].isna().any() + + +# -- transform_coordinates_to_volume --------------------------------------- + + +def test_coordinates_are_counted_into_their_voxels(): + coords = np.array([ + [0.0, 0.0, 0.0], + [5.0, 5.0, 5.0], # same voxel as the first + [10.0, 0.0, 0.0], + ]) + + volume = transform_coordinates_to_volume( + coords, volume_shape=(3, 2, 2), resolution=(10, 10, 10) + ) + + assert volume.shape == (3, 2, 2) + assert volume[0, 0, 0] == 2 + assert volume[1, 0, 0] == 1 + assert volume.sum() == 3 + + +def test_counts_are_invariant_to_scaling_coordinates_and_resolution_together(): + coords = np.array([[0.0, 0.0, 0.0], [10.0, 10.0, 10.0]]) + + at_10 = transform_coordinates_to_volume(coords, (3, 3, 3), (10, 10, 10)) + at_20 = transform_coordinates_to_volume(coords * 2, (3, 3, 3), (20, 20, 20)) + + assert np.array_equal(at_10, at_20) + + +def test_empty_coordinates_give_an_empty_volume(): + volume = transform_coordinates_to_volume( + np.zeros((0, 3)), volume_shape=(2, 2, 2), resolution=(10, 10, 10) + ) + assert volume.sum() == 0 + + +# -- transform_swc_to_volume ----------------------------------------------- + + +def test_swc_nodes_are_counted_into_a_volume(swc_file): + volume = transform_swc_to_volume( + swc_file, volume_shape=(3, 2, 2), resolution=(10, 10, 10) + ) + + assert volume.sum() == 4 + assert volume[0, 0, 0] == 1 # the soma at the origin + assert volume[1, 0, 0] == 1 + assert volume[2, 0, 0] == 1 + assert volume[1, 1, 0] == 1 + + +def test_compartments_filter_which_nodes_are_counted(swc_file): + soma_only = transform_swc_to_volume( + swc_file, volume_shape=(3, 2, 2), resolution=(10, 10, 10), compartments=[1] + ) + dendrite_only = transform_swc_to_volume( + swc_file, volume_shape=(3, 2, 2), resolution=(10, 10, 10), compartments=[3] + ) + + assert soma_only.sum() == 1 + assert dendrite_only.sum() == 3 + assert soma_only[0, 0, 0] == 1 + assert dendrite_only[0, 0, 0] == 0 + + +# -- find_topological_point_coordinates ------------------------------------ + + +def test_branch_and_termination_nodes_are_found(): + """Node 2 has two children, so it branches; 3 and 4 are nobody's parent.""" + df = pd.DataFrame({ + "id": [1, 2, 3, 4], + "type": [1, 3, 3, 3], + "x": [0.0, 10.0, 20.0, 10.0], + "y": [0.0, 0.0, 0.0, 10.0], + "z": [0.0, 0.0, 0.0, 0.0], + "r": [1.0, 1.0, 1.0, 1.0], + "parent_id": [-1, 1, 2, 2], + }) + + coords = find_topological_point_coordinates(df) + + # Branch nodes first, then terminations. + assert coords.shape == (3, 3) + assert np.array_equal(coords[0], np.array([10.0, 0.0, 0.0])) # branch: node 2 + assert {tuple(c) for c in coords[1:]} == { + (20.0, 0.0, 0.0), # termination: node 3 + (10.0, 10.0, 0.0), # termination: node 4 + } + + +def test_an_unbranched_neurite_has_only_a_termination(): + df = pd.DataFrame({ + "id": [1, 2, 3], + "type": [1, 3, 3], + "x": [0.0, 10.0, 20.0], + "y": [0.0, 0.0, 0.0], + "z": [0.0, 0.0, 0.0], + "r": [1.0, 1.0, 1.0], + "parent_id": [-1, 1, 2], + }) + + coords = find_topological_point_coordinates(df) + + assert coords.shape == (1, 3) + assert np.array_equal(coords[0], np.array([20.0, 0.0, 0.0])) + + +def test_topological_points_can_be_transformed_to_a_volume(swc_file): + """The documented pairing of the two functions.""" + df = load_swc_as_dataframe(swc_file) + coords = find_topological_point_coordinates(df) + + volume = transform_coordinates_to_volume(coords, (3, 2, 2), (10, 10, 10)) + + assert volume.sum() == coords.shape[0] diff --git a/tests/test_processing.py b/tests/test_processing.py new file mode 100644 index 0000000..bfae781 --- /dev/null +++ b/tests/test_processing.py @@ -0,0 +1,86 @@ +"""Path cleaning, against paths with known repeats.""" + +import numpy as np +import pytest + +from ccf_streamlines.processing import remove_duplicate_voxels_from_paths + + +def test_consecutive_duplicates_are_removed_and_the_row_is_refilled(): + """Duplicates collapse and the row is left-packed, padded back to width.""" + paths = np.array([ + [5, 5, 7, 7, 9, 0, 0, 0], + [2, 3, 4, 0, 0, 0, 0, 0], + ]) + + result = remove_duplicate_voxels_from_paths(paths) + + assert result.shape == paths.shape + assert np.array_equal(result[0], np.array([5, 7, 9, 0, 0, 0, 0, 0])) + assert np.array_equal(result[1], np.array([2, 3, 4, 0, 0, 0, 0, 0])) + + +def test_a_path_with_no_duplicates_is_unchanged(): + paths = np.array([[11, 12, 13, 14, 0, 0]]) + assert np.array_equal(remove_duplicate_voxels_from_paths(paths), paths) + + +def test_only_consecutive_duplicates_are_removed(): + """A voxel repeated non-consecutively is a genuine revisit, and is kept.""" + paths = np.array([[4, 4, 8, 4, 0, 0]]) + result = remove_duplicate_voxels_from_paths(paths) + assert np.array_equal(result[0], np.array([4, 8, 4, 0, 0, 0])) + + +def test_a_long_run_collapses_to_one_voxel(): + paths = np.array([[6, 6, 6, 6, 6, 0]]) + result = remove_duplicate_voxels_from_paths(paths) + assert np.array_equal(result[0], np.array([6, 0, 0, 0, 0, 0])) + + +def test_rows_are_cleaned_independently(): + paths = np.array([ + [1, 1, 1, 2, 0, 0], + [3, 4, 5, 6, 7, 0], + ]) + result = remove_duplicate_voxels_from_paths(paths) + assert np.array_equal(result[0], np.array([1, 2, 0, 0, 0, 0])) + assert np.array_equal(result[1], np.array([3, 4, 5, 6, 7, 0])) + + +def test_input_must_be_zero_padded(): + """A row that fills its full width loses its last voxel. + + The function detects voxels by looking at where the row *changes*, so the + final voxel is only seen because the padding after it differs. The real + ``paths`` dataset is always zero-padded, so this never bites in practice -- + but it is a precondition, not a free choice, and a caller building paths by + hand needs to know. + """ + unpadded = np.array([[5, 5, 7, 7, 9, 9]]) + result = remove_duplicate_voxels_from_paths(unpadded) + assert 9 not in result[0].tolist() + + padded = np.array([[5, 5, 7, 7, 9, 9, 0]]) + assert 9 in remove_duplicate_voxels_from_paths(padded)[0].tolist() + + +def test_matches_the_inline_deduplication_in_metrics(): + """`metrics.measure_streamline_layer_thicknesses` re-implements this with a + Python loop instead of calling it. The two must agree, or layer depths are + measured against different paths than everything else.""" + paths = np.array([ + [5, 5, 7, 7, 9, 0, 0, 0], + [2, 3, 3, 0, 0, 0, 0, 0], + ]) + + vectorized = remove_duplicate_voxels_from_paths(paths.copy()) + + # The loop from metrics.py, verbatim. + inline = np.zeros_like(paths) + paths_diff = np.diff(paths, axis=1) + for i in range(paths.shape[0]): + unique_inds = np.flatnonzero(paths_diff[i, :]) + inline[i, :len(unique_inds)] = paths[i, :][unique_inds] + + assert np.array_equal(vectorized, inline) From 5bbc79dd88c74cf23ba5ec50ac855ab376279bfc Mon Sep 17 00:00:00 2001 From: Nathan Gouwens Date: Mon, 31 Aug 2026 12:39:50 -0700 Subject: [PATCH 2/2] Cite the tracker issue in the xfail reason The defect now has an issue (#23), so the marker names it. A contributor looking at a reported unexpected pass can find the report, and closing the issue is verifiable by CI rather than by reading the diff. Co-Authored-By: Claude Opus 5 --- tests/test_angle.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/tests/test_angle.py b/tests/test_angle.py index b9e19e2..b4b38b5 100644 --- a/tests/test_angle.py +++ b/tests/test_angle.py @@ -15,9 +15,10 @@ ) RESOLUTION_NOT_FORWARDED = ( - "`find_closest_streamline` accepts `resolution` but does not forward it to " - "`coordinates_to_voxels`, so a non-default resolution voxelises against " - "(10, 10, 10); remove this marker when it is fixed" + "AllenInstitute/ccf_streamlines#23: `find_closest_streamline` accepts " + "`resolution` but does not forward it to `coordinates_to_voxels`, so a " + "non-default resolution voxelises against (10, 10, 10); remove this marker " + "when it is fixed" ) #: Maps the unit square onto the xy-plane, so the plane normal is +z.