Merge pull request #1235 from EmmaRenauld/unit_bundles

First unit tests in tractoanalysis module

Author: arnaudbore

src/scilpy/reconst/sh.py

@@ -285,8 +285,12 @@ def peaks_from_sh(shm_coeff, sphere, mask=None, relative_peak_threshold=0.5,
 
     Returns
     -------
-    tuple of np.ndarray
-        peak_dirs, peak_values, peak_indices
+    peak_dirs: np.ndarray
+        Shape: [x, y, z, npeaks * 3]
+    peak_values: np.ndarray
+        Shape: [x, y, z, npeaks]
+    peak_indices: np.ndarray
+        Shape: [x, y, z, npeaks]
     """
     sh_order = order_from_ncoef(shm_coeff.shape[-1],
                                 full_basis=full_basis)

src/scilpy/tractanalysis/afd_along_streamlines.py

@@ -13,9 +13,13 @@
 def afd_map_along_streamlines(sft, fodf, fodf_basis, length_weighting,
                               is_legacy=True):
     """
-    Compute the mean Apparent Fiber Density (AFD) and mean Radial fODF
+    Compute the mean Apparent Fiber Density (AFD) [1] and mean Radial fODF
     (radfODF) maps along a bundle.
 
+    [1] Raffelt et. al (2012). Apparent fibre density: a novel measure for the
+    analysis of diffusion-weighted magnetic resonance images.
+    Neuroimage, 59(4), 3976-3994.
+
     Parameters
     ----------
     sft : StatefulTractogram
@@ -27,6 +31,8 @@ def afd_map_along_streamlines(sft, fodf, fodf_basis, length_weighting,
         Has to be descoteaux07 or tournier07
     length_weighting : bool
         If set, will weigh the AFD values according to segment lengths
+    is_legacy : bool
+        If true, uses legacy basis.
 
     Returns
     -------
@@ -37,9 +43,9 @@ def afd_map_along_streamlines(sft, fodf, fodf_basis, length_weighting,
     """
 
     afd_sum, rd_sum, weights = \
-        afd_and_rd_sums_along_streamlines(sft, fodf, fodf_basis,
-                                          length_weighting,
-                                          is_legacy=is_legacy)
+        _afd_and_rd_sums_along_streamlines(sft, fodf, fodf_basis,
+                                           length_weighting,
+                                           is_legacy=is_legacy)
 
     non_zeros = np.nonzero(afd_sum)
     weights_nz = weights[non_zeros]
@@ -49,8 +55,8 @@ def afd_map_along_streamlines(sft, fodf, fodf_basis, length_weighting,
     return afd_sum, rd_sum
 
 
-def afd_and_rd_sums_along_streamlines(sft, fodf, fodf_basis,
-                                      length_weighting, is_legacy=True):
+def _afd_and_rd_sums_along_streamlines(sft, fodf, fodf_basis,
+                                       length_weighting, is_legacy=True):
     """
     Compute the mean Apparent Fiber Density (AFD) and
     mean Radial fODF (radfODF) maps along a bundle.
@@ -67,7 +73,7 @@ def afd_and_rd_sums_along_streamlines(sft, fodf, fodf_basis,
     length_weighting : bool
         If set, will weigh the AFD values according to segment lengths.
     is_legacy : bool, optional
-        Whether or not the SH basis is in its legacy form.
+        Whether the SH basis is in its legacy form.
 
     Returns
     -------
@@ -89,42 +95,46 @@ def afd_and_rd_sums_along_streamlines(sft, fodf, fodf_basis,
     _, n = sph_harm_ind_list(order)
     legendre0_at_n = legendre_p_all(order, 0)[0][n]
     sphere_norm = np.linalg.norm(sphere.vertices)
+    p_matrix = np.eye(fodf_data.shape[3]) * legendre0_at_n
 
+    # Initializing
     afd_sum_map = np.zeros(shape=fodf_data.shape[:-1])
     rd_sum_map = np.zeros(shape=fodf_data.shape[:-1])
     weight_map = np.zeros(shape=fodf_data.shape[:-1])
 
-    p_matrix = np.eye(fodf_data.shape[3]) * legendre0_at_n
     all_split_streamlines =\
         subdivide_streamlines_at_voxel_faces(sft.streamlines)
     for split_streamlines in all_split_streamlines:
+        # Get the direction of each segment
         segments = split_streamlines[1:] - split_streamlines[:-1]
         seg_lengths = np.linalg.norm(segments, axis=1)
 
-        # Remove points where the segment is zero.
+        # Remove segments of length zero.
         # This removes numpy warnings of division by zero.
         non_zero_lengths = np.nonzero(seg_lengths)[0]
         segments = segments[non_zero_lengths]
         seg_lengths = seg_lengths[non_zero_lengths]
 
+        # Find closest point on sphere
         test = np.dot(segments, sphere.vertices.T)
         test2 = (test.T / (seg_lengths * sphere_norm)).T
         angles = np.arccos(test2)
         sorted_angles = np.argsort(angles, axis=1)
         closest_vertex_indices = sorted_angles[:, 0]
 
-        # Those starting points are used for the segment vox_idx computations
-        strl_start = split_streamlines[non_zero_lengths]
-        vox_indices = (strl_start + (0.5 * segments)).astype(int)
+        # Get the middle voxel of each segment
+        # (They are already cut per voxel, so the middle voxel is probably the
+        # same as the start voxel.)
+        seg_starts = split_streamlines[non_zero_lengths]
+        vox_indices = (seg_starts + (0.5 * segments)).astype(int)
 
         normalization_weights = np.ones_like(seg_lengths)
         if length_weighting:
             normalization_weights = seg_lengths / \
                 np.linalg.norm(fodf.header.get_zooms()[:3])
 
-        for vox_idx, closest_vertex_index, norm_weight in zip(vox_indices,
-                                                              closest_vertex_indices,
-                                                              normalization_weights):
+        for vox_idx, closest_vertex_index, norm_weight in zip(
+                vox_indices, closest_vertex_indices, normalization_weights):
             vox_idx = tuple(vox_idx)
             b_at_idx = b_matrix.T[closest_vertex_index]
             fodf_at_index = fodf_data[vox_idx]

src/scilpy/tractanalysis/tests/test_afd_along_streamlines.py

@@ -0,0 +1,89 @@
+# -*- coding: utf-8 -*-
+import os
+
+import nibabel as nib
+import numpy as np
+from dipy.io.stateful_tractogram import StatefulTractogram, Space, Origin
+from scilpy.tractanalysis.streamlines_metrics import compute_tract_counts_map
+
+from scilpy.tractanalysis.afd_along_streamlines import \
+    afd_map_along_streamlines
+from scilpy import SCILPY_HOME
+from scilpy.io.fetcher import fetch_data, get_testing_files_dict
+
+
+fetch_data(get_testing_files_dict(), keys=['processing.zip'])
+
+
+def test_afd_map_along_streamlines():
+    fodf_img = nib.load(os.path.join(SCILPY_HOME, 'processing', 'fodf.nii.gz'))
+    fodf = fodf_img.get_fdata(dtype=np.float32)
+    fodf_mask = np.sum(fodf, axis=-1) > 0
+
+    peaks = nib.load(os.path.join(SCILPY_HOME, 'processing', 'peaks.nii.gz'))
+    peaks = peaks.get_fdata(dtype=np.float32)
+    peaks = peaks[:, :, :, 0:3]  # Keeping only the first peak
+
+    # Creating a few very short streamlines (2 points) aligned on peaks
+    # in the middle of the image (shape 57x67x56)
+    streamlines = []
+    for i in range(23, 26):
+        start_pt = np.asarray([i + 0.1, i + 0.1, i + 0.1])
+        peak = peaks[tuple(start_pt.astype(int))]
+        peak = np.asarray(peak) / np.linalg.norm(peak)
+        second_pt = np.asarray(start_pt) + peak
+        streamlines.append(np.vstack([start_pt, second_pt], dtype=np.float32))
+    sft = StatefulTractogram(streamlines=streamlines,
+                             reference=fodf_img,
+                             space=Space.VOX, origin=Origin('corner'))
+
+    # Processing
+    afd_map, rd_map = afd_map_along_streamlines(
+        sft, fodf_img, 'descoteaux07',
+        length_weighting=False, is_legacy=False)
+
+    # Should have the same size as fodf
+    assert np.array_equal(fodf_img.shape[0:3], afd_map.shape)
+
+    # Should only have values where we have streamlines + inside fodf mask.
+    # Note. If the streamline starts exactly on the border of the voxel,
+    # we get weird results, because the code for afd uses segments, which means
+    # that the coordinate [i, i, i] may belong to the voxel [i, i, i] or to its
+    # neighbor, depending on the direction of the segment, but the density_map
+    # always considers it to be in voxel [i, i, i]. But ok, streamlines exactly
+    # on a border should not happen too often in real life.
+    density_map = compute_tract_counts_map(sft.streamlines, sft.dimensions)
+    density_mask = density_map > 0
+    mask = np.logical_and(density_mask, fodf_mask)
+    assert np.array_equal(mask, afd_map != 0)
+
+    # Now processing again, but with more of the same streamlines
+    # AFD will be the same (divided by weight map = by number of streamlines)
+    streamlines2 = [s for s in streamlines]
+    for translation in [0.1, 0.11, 0.12]:
+        streamlines2 += [s + translation for s in streamlines]
+    sft.streamlines = streamlines2
+    afd_map2, rd_map2 = afd_map_along_streamlines(
+        sft, fodf_img, 'descoteaux07',
+        length_weighting=False, is_legacy=False)
+    assert np.count_nonzero(afd_map) == np.count_nonzero(afd_map2)
+    assert np.count_nonzero(rd_map) == np.count_nonzero(rd_map2)
+    assert np.allclose(afd_map, afd_map2)
+    assert np.allclose(rd_map, rd_map2)
+
+    # Now processing again, but with rotated peaks. Values should be in the
+    # same voxels, but different
+    streamlines2 = [s for s in streamlines]
+    for s in streamlines:
+        # Same starting point. Moving a little the end point, but we want to
+        # stay in the same voxel
+        modified_s = [s[0, :], s[1, :] + 0.1]
+        streamlines2.append(modified_s)
+    sft.streamlines = streamlines2
+    afd_map2, rd_map2 = afd_map_along_streamlines(
+        sft, fodf_img, 'descoteaux07',
+        length_weighting=False, is_legacy=False)
+    assert np.count_nonzero(afd_map) == np.count_nonzero(afd_map2)
+    assert np.count_nonzero(rd_map) == np.count_nonzero(rd_map2)
+    assert not np.allclose(afd_map, afd_map2)
+    assert not np.allclose(rd_map, rd_map2)

src/scilpy/tractanalysis/tests/test_voxel_boundary_intersection.py

@@ -0,0 +1,40 @@
+# -*- coding: utf-8 -*-
+import numpy as np
+import nibabel as nib
+from dipy.io.stateful_tractogram import StatefulTractogram, Origin, Space
+
+from scilpy.tractanalysis.voxel_boundary_intersection import \
+    subdivide_streamlines_at_voxel_faces
+
+
+def test_subdivide():
+    # Creating a few very short streamlines (2 points).
+
+    # streamline1: stays inside its voxel
+    streamline1 = np.asarray([[23.0, 23.0, 23.0],
+                              [23.5, 23.5, 23.5]], dtype=np.float32)
+
+    # streamline2: changes to next voxel on second dim only. Covers 2 voxels
+    streamline2 = np.asarray([[23.0, 23.1, 23.0],
+                              [23.5, 22.2, 23.5]], dtype=np.float32)
+
+    # streamline3: covers 11 voxels
+    streamline3 = np.asarray([[0.0, 0.0, 0.0],
+                              [10.99, 0.0, 0.0]], dtype=np.float32)
+
+    fake_ref = nib.Nifti1Image(np.zeros((3, 3, 3)), affine=np.eye(4))
+    sft = StatefulTractogram([streamline1, streamline2,
+                              streamline3], reference=fake_ref,
+                             origin=Origin('corner'), space=Space('vox'))
+    split = subdivide_streamlines_at_voxel_faces(sft.streamlines)
+
+    # Streamline 1 should stay the same
+    assert np.array_equal(split[0], streamline1)
+
+    # Streamline 2 should be split into 2 (one more point)
+    assert split[1].shape[0] == 3
+    assert np.array_equal(split[1][0, :], streamline2[0, :])
+    assert np.array_equal(split[1][-1, :], streamline2[-1, :])
+
+    # streamline 3 should be split into 11 (12 points)
+    assert split[2].shape[0] == 12