Merge branch 'master' of github.com:VincentBeaud/scilpy into dev

.github/workflows/test.yml

@@ -7,6 +7,7 @@ on:
   pull_request:
     branches:
       - master
+  merge_group:
 
 concurrency:
   group: ${{ github.workflow }}-${{ github.event.pull_request.number || github.ref }}
@@ -70,7 +71,7 @@ jobs:
             .test_reports/
 
   coverage:
-    runs-on: scilus-runners
+    runs-on: ubuntu-latest
     if: github.repository == 'scilus/scilpy'
     needs: test
 

docs/source/fake_files/voxel_boundary_intersection.py

@@ -0,0 +1,5 @@
+# -*- coding: utf-8 -*-
+
+
+def subdivide_streamlines_at_voxel_faces():
+    pass

docs/source/modules/scilpy.tractanalysis.rst

@@ -34,10 +34,10 @@ scilpy.tractanalysis.features module
    :show-inheritance:
 
 
-scilpy.tractanalysis.grid\_intersections module
+scilpy.tractanalysis.voxel\_boundary\_intersection module
 ------------------------------------------------------
 
-.. automodule:: scilpy.tractanalysis.grid_intersections
+.. automodule:: scilpy.tractanalysis.voxel_boundary_intersection
    :members:
    :undoc-members:
    :show-inheritance:

scilpy/connectivity/connectivity.py

@@ -10,14 +10,18 @@
 import h5py
 import nibabel as nib
 import numpy as np
+from scipy.ndimage import map_coordinates
 
 from scilpy.image.labels import get_data_as_labels
 from scilpy.io.hdf5 import reconstruct_streamlines_from_hdf5
 from scilpy.tractanalysis.reproducibility_measures import \
     compute_bundle_adjacency_voxel
 from scilpy.tractanalysis.streamlines_metrics import compute_tract_counts_map
+from scilpy.tractograms.streamline_operations import \
+    resample_streamlines_num_points
 from scilpy.utils.metrics_tools import compute_lesion_stats
 
+
 d = threading.local()
 
 
@@ -31,8 +35,8 @@ def compute_triu_connectivity_from_labels(tractogram, data_labels,
     ----------
     tractogram: StatefulTractogram, or list[np.ndarray]
         Streamlines. A StatefulTractogram input is recommanded.
-        When using directly with a list of streamlines, streamlinee must be in
-        vox space, corner origin.
+        When using directly with a list of streamlines, streamlines must be in
+        vox space, center origin.
     data_labels: np.ndarray
         The loaded nifti image.
     keep_background: Bool
@@ -55,12 +59,12 @@ def compute_triu_connectivity_from_labels(tractogram, data_labels,
         For each streamline, the label at ending point.
     """
     if isinstance(tractogram, StatefulTractogram):
-        #  Vox space, corner origin
-        # = we can get the nearest neighbor easily.
-        # Coord 0 = voxel 0. Coord 0.9 = voxel 0. Coord 1 = voxel 1.
-        tractogram.to_vox()
-        tractogram.to_corner()
-        streamlines = tractogram.streamlines
+        # vox space, center origin: compatible with map_coordinates
+        sfs_2_pts = resample_streamlines_num_points(tractogram, 2)
+        sfs_2_pts.to_vox()
+        sfs_2_pts.to_center()
+        streamlines = sfs_2_pts.streamlines
+
     else:
         streamlines = tractogram
 
@@ -71,23 +75,16 @@ def compute_triu_connectivity_from_labels(tractogram, data_labels,
                   .format(nb_labels))
 
     matrix = np.zeros((nb_labels, nb_labels), dtype=int)
-    start_labels = []
-    end_labels = []
-
-    for s in streamlines:
-        start = ordered_labels.index(
-            data_labels[tuple(np.floor(s[0, :]).astype(int))])
-        end = ordered_labels.index(
-            data_labels[tuple(np.floor(s[-1, :]).astype(int))])
 
-        start_labels.append(start)
-        end_labels.append(end)
+    labels = map_coordinates(data_labels, streamlines._data.T, order=0)
+    start_labels = labels[0::2]
+    end_labels = labels[1::2]
 
-        matrix[start, end] += 1
-        if start != end:
-            matrix[end, start] += 1
+    # sort each pair of labels for start to be smaller than end
+    start_labels, end_labels = zip(*[sorted(pair) for pair in
+                                     zip(start_labels, end_labels)])
 
-    matrix = np.triu(matrix)
+    np.add.at(matrix, (start_labels, end_labels), 1)
     assert matrix.sum() == len(streamlines)
 
     # Rejecting background
@@ -249,7 +246,7 @@ def compute_connectivity_matrices_from_hdf5(
 
     if compute_volume:
         measures_to_return['volume_mm3'] = np.count_nonzero(density) * \
-                                           np.prod(voxel_sizes)
+            np.prod(voxel_sizes)
 
     if compute_streamline_count:
         measures_to_return['streamline_count'] = len(streamlines)

scilpy/io/fetcher.py

@@ -5,9 +5,8 @@
 import logging
 import os
 import pathlib
-import zipfile
-
 import requests
+import zipfile
 
 from scilpy import SCILPY_HOME
 

scilpy/io/hdf5.py

@@ -105,7 +105,7 @@ def reconstruct_sft_from_hdf5(hdf5_handle, group_keys, space=Space.VOX,
             for sub_key in hdf5_handle[group_key].keys():
                 if sub_key not in ['data', 'offsets', 'lengths']:
                     data = hdf5_handle[group_key][sub_key]
-                    if data.shape == hdf5_handle[group_key]['offsets']:
+                    if data.shape == hdf5_handle[group_key]['offsets'].shape:
                         # Discovered dps
                         if load_dps:
                             if i == 0 or not merge_groups:

scilpy/io/streamlines.py

@@ -75,8 +75,8 @@ def load_tractogram_with_reference(parser, args, filepath, arg_name=None):
     filepath: str
         Path of the tractogram file.
     arg_name: str, optional
-        Name of the reference argument. By default the args.ref is used. If
-        arg_name is given, then args.arg_name_ref will be used instead.
+        Name of the reference argument. By default the args.reference is used.
+        If arg_name is given, then args.arg_name_ref will be used instead.
     """
     if is_argument_set(args, 'bbox_check'):
         bbox_check = args.bbox_check

scilpy/tracking/seed.py

@@ -371,21 +371,27 @@ class CustomSeedsDispenser(SeedGenerator):
     Adaptation of the scilpy.tracking.seed.SeedGenerator interface for
     using already generated, custom seeds.
     """
-    def __init__(self, custom_seeds, voxres, space=Space('vox'),
+    def __init__(self, custom_seeds, space=Space('vox'),
                  origin=Origin('center')):
         """
+        Custom seeds need to be in the same space and origin as the ODFs used
+        for tracking.
+
         Parameters
         ----------
         custom_seeds: list
             Custom seeding coordinates.
-        voxres: np.ndarray(3,)
-            The pixel resolution, ex, using img.header.get_zooms()[:3].
+        space: Space (optional)
+            The Dipy space in which the seeds were saved.
+            Default: Space.Vox or 'vox'
+        origin: Origin (optional)
+            The Dipy origin in which the seeds were saved.
+            Default: Origin.NIFTI or 'center'
         """
-        self.voxres = voxres
         self.origin = origin
         self.space = space
         self.seeds = custom_seeds
-        self.count = 0
+        self.i = 0
 
     def init_generator(self, rng_seed, numbers_to_skip):
         """

scilpy/tracking/utils.py

@@ -129,8 +129,9 @@ def add_seeding_options(p):
     seed_sub_exclusive.add_argument(
         '--in_custom_seeds', type=str,
         help='Path to a file containing a list of custom seeding \n'
-             'coordinates (.txt, .mat or .npy). They should be in'
-             'voxel space.')
+             'coordinates (.txt, .mat or .npy). They should be in \n'
+             'voxel space. In the case of a text file, each line should \n'
+             'contain a single seed, written in the format: [x, y, z].')
 
 
 def add_out_options(p):

scilpy/tractanalysis/afd_along_streamlines.py

@@ -6,7 +6,8 @@
 from scipy.special import lpn
 
 from scilpy.reconst.utils import find_order_from_nb_coeff
-from scilpy.tractanalysis.grid_intersections import grid_intersections
+from scilpy.tractanalysis.voxel_boundary_intersection import\
+    subdivide_streamlines_at_voxel_faces
 
 
 def afd_map_along_streamlines(sft, fodf, fodf_basis, length_weighting,
@@ -94,9 +95,10 @@ def afd_and_rd_sums_along_streamlines(sft, fodf, fodf_basis,
     weight_map = np.zeros(shape=fodf_data.shape[:-1])
 
     p_matrix = np.eye(fodf_data.shape[3]) * legendre0_at_n
-    all_crossed_indices = grid_intersections(sft.streamlines)
-    for crossed_indices in all_crossed_indices:
-        segments = crossed_indices[1:] - crossed_indices[:-1]
+    all_split_streamlines =\
+        subdivide_streamlines_at_voxel_faces(sft.streamlines)
+    for split_streamlines in all_split_streamlines:
+        segments = split_streamlines[1:] - split_streamlines[:-1]
         seg_lengths = np.linalg.norm(segments, axis=1)
 
         # Remove points where the segment is zero.
@@ -112,7 +114,7 @@ def afd_and_rd_sums_along_streamlines(sft, fodf, fodf_basis,
         closest_vertex_indices = sorted_angles[:, 0]
 
         # Those starting points are used for the segment vox_idx computations
-        strl_start = crossed_indices[non_zero_lengths]
+        strl_start = split_streamlines[non_zero_lengths]
         vox_indices = (strl_start + (0.5 * segments)).astype(int)
 
         normalization_weights = np.ones_like(seg_lengths)

scilpy/tractanalysis/bingham_metric_along_streamlines.py

@@ -2,7 +2,8 @@
 
 import numpy as np
 from scilpy.reconst.bingham import bingham_to_peak_direction
-from scilpy.tractanalysis.grid_intersections import grid_intersections
+from scilpy.tractanalysis.voxel_boundary_intersection import\
+    subdivide_streamlines_at_voxel_faces
 
 
 def bingham_metric_map_along_streamlines(sft, bingham_coeffs,
@@ -73,9 +74,10 @@ def bingham_metric_sum_along_streamlines(sft, bingham_coeffs, metric,
     weight_map = np.zeros(metric.shape[:-1])
     min_cos_theta = np.cos(np.radians(max_theta))
 
-    all_crossed_indices = grid_intersections(sft.streamlines)
-    for crossed_indices in all_crossed_indices:
-        segments = crossed_indices[1:] - crossed_indices[:-1]
+    all_split_streamlines =\
+        subdivide_streamlines_at_voxel_faces(sft.streamlines)
+    for split_streamlines in all_split_streamlines:
+        segments = split_streamlines[1:] - split_streamlines[:-1]
         seg_lengths = np.linalg.norm(segments, axis=1)
 
         # Remove points where the segment is zero.
@@ -85,7 +87,7 @@ def bingham_metric_sum_along_streamlines(sft, bingham_coeffs, metric,
         seg_lengths = seg_lengths[non_zero_lengths]
 
         # Those starting points are used for the segment vox_idx computations
-        seg_start = crossed_indices[non_zero_lengths]
+        seg_start = split_streamlines[non_zero_lengths]
         vox_indices = (seg_start + (0.5 * segments)).astype(int)
 
         normalization_weights = np.ones_like(seg_lengths)

scilpy/tractanalysis/fibertube_scoring.py

@@ -465,6 +465,10 @@ def endpoint_connectivity(blur_radius, centerlines, centerlines_length,
     segment", which is the closest fibertube segment to its before-last
     coordinate.
 
+    IMPORTANT: Streamlines given as input to be scored should be forward-only,
+    which means they are saved so that [0] is the seeding position and [-1] is
+    the end.
+
     VC: "Valid Connection": A streamline whose arrival fibertube segment is
     the final segment of the fibertube in which is was originally seeded.
 
@@ -492,7 +496,7 @@ def endpoint_connectivity(blur_radius, centerlines, centerlines_length,
         Fixed array containing the number of coordinates of each streamline
     seeds_fiber: list
         Array of the same length as there are streamlines. For every
-        streamline, contains the index of the fiber in which it has been
+        streamline, contains the index of the fibertube in which it has been
         seeded.
 
     Return
@@ -522,43 +526,44 @@ def endpoint_connectivity(blur_radius, centerlines, centerlines_length,
     ic = set()
     nc = set()
 
-    # streamline[1] is the last point with a valid direction
+    # streamline[-2] is the last point with a valid direction
     all_neighbors = tree.query_radius(
-        streamlines[:, 1], blur_radius + max_seg_length / 2 + max_diameter)
+        streamlines[:, -2], blur_radius + max_seg_length / 2 + max_diameter)
 
-    for si, streamline in enumerate(streamlines):
-        seed_fi = seeds_fiber[si]
-        neighbors = all_neighbors[si]
+    for streamline_index, streamline in enumerate(streamlines):
+        seed_fi = seeds_fiber[streamline_index]
+        neighbors = all_neighbors[streamline_index]
 
-        min_dist = np.inf
-        min_seg = 0
+        closest_dist = np.inf
+        closest_seg = 0
 
         # Finding closest segment
         # There will always be a neighbor to override np.inf
-        for neighbor_segi in neighbors:
-            fi = indices[neighbor_segi][0]
-            pi = indices[neighbor_segi][1]
+        for segment_index in neighbors:
+            fibertube_index = indices[segment_index][0]
+            point_index = indices[segment_index][1]
 
-            dist, _, _ = dist_point_segment(centerlines[fi][pi],
-                                            centerlines[fi][pi+1],
-                                            streamline[1])
+            dist, _, _ = dist_point_segment(
+                centerlines[fibertube_index][point_index],
+                centerlines[fibertube_index][point_index+1],
+                streamline[-2])
 
-            if dist < min_dist:
-                min_dist = dist
-                min_seg = neighbor_segi
+            if dist < closest_dist:
+                closest_dist = dist
+                closest_seg = segment_index
 
-        fi = indices[min_seg][0]
-        pi = indices[min_seg][1]
+        fibertube_index = indices[closest_seg][0]
+        point_index = indices[closest_seg][1]
 
         # If the closest segment is the last of its centerlines
-        if pi == centerlines_length[fi]-1:
-            if fi == seed_fi:
-                vc.add(si)
+        if point_index == centerlines_length[fibertube_index]-1:
+            if fibertube_index == seed_fi:
+                vc.add(streamline_index)
             else:
-                ic.add(si)
-        elif pi == 0:
-            ic.add(si)
+                ic.add(streamline_index)
+        elif point_index == 0:
+            ic.add(streamline_index)
         else:
-            nc.add(si)
+            nc.add(streamline_index)
 
     return list(vc), list(ic), list(nc)

scilpy/tractanalysis/fixel_density.py

@@ -3,7 +3,8 @@
 import numpy as np
 
 from dipy.io.streamline import load_tractogram
-from scilpy.tractanalysis.grid_intersections import grid_intersections
+from scilpy.tractanalysis.voxel_boundary_intersection import\
+    subdivide_streamlines_at_voxel_faces
 
 
 def _fixel_density_parallel(args):
@@ -21,9 +22,10 @@ def _fixel_density_single_bundle(bundle, peaks, max_theta, dps_key):
 
     min_cos_theta = np.cos(np.radians(max_theta))
 
-    all_crossed_indices = grid_intersections(sft.streamlines)
-    for i, crossed_indices in enumerate(all_crossed_indices):
-        segments = crossed_indices[1:] - crossed_indices[:-1]
+    all_split_streamlines =\
+        subdivide_streamlines_at_voxel_faces(sft.streamlines)
+    for i, split_streamlines in enumerate(all_split_streamlines):
+        segments = split_streamlines[1:] - split_streamlines[:-1]
         seg_lengths = np.linalg.norm(segments, axis=1)
 
         # Remove points where the segment is zero.
@@ -33,7 +35,7 @@ def _fixel_density_single_bundle(bundle, peaks, max_theta, dps_key):
         seg_lengths = seg_lengths[non_zero_lengths]
 
         # Those starting points are used for the segment vox_idx computations
-        seg_start = crossed_indices[non_zero_lengths]
+        seg_start = split_streamlines[non_zero_lengths]
         vox_indices = (seg_start + (0.5 * segments)).astype(int)
 
         normalized_seg = np.reshape(segments / seg_lengths[..., None], (-1, 3))