Test the pure computational functions

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))

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
+    )

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),
+        )