diff --git a/tests/test_angle.py b/tests/test_angle.py new file mode 100644 index 0000000..b4b38b5 --- /dev/null +++ b/tests/test_angle.py @@ -0,0 +1,219 @@ +"""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 = ( + "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. +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)