Merge remote-tracking branch 'upstream/master' into atheb/s0_div_0

.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/documentation/construct_participants_tsv_file.rst

@@ -1,4 +1,4 @@
-Instructions to write the tsv files "participants.tsv" for the script scil_group_comparison.py
+Instructions to write the tsv files "participants.tsv" for the script scil_stats_group_comparison.py
 ===============================================================================================
 
 The TSV file should follow the BIDS `specification <https://bids-specification.readthedocs.io/en/stable/03-modality-agnostic-files.html#participants-file>`_.
@@ -12,7 +12,7 @@ participant_id categorical_var_1 categorical_var_2 ...
 
 (ex: participant_id sex nb_children)
 
-The categorical variable name are the "group_by" variable that can be called by scil_group_comparison.py
+The categorical variable name are the "group_by" variable that can be called by scil_stats_group_comparison.py
 
 Specific row
 ------------

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

@@ -364,3 +364,70 @@ def get_next_n_pos(self, random_generator, shuffled_indices,
                                 random_generator)
 
         return seeds
+
+
+class CustomSeedsDispenser(SeedGenerator):
+    """
+    Adaptation of the scilpy.tracking.seed.SeedGenerator interface for
+    using already generated, custom seeds.
+    """
+    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.
+        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.origin = origin
+        self.space = space
+        self.seeds = custom_seeds
+        self.i = 0
+
+    def init_generator(self, rng_seed, numbers_to_skip):
+        """
+        Does not do anything. Simulates SeedGenerator's implementation for
+        retro-compatibility.
+
+        Parameters
+        ----------
+        rng_seed : int
+            The "seed" for the random generator.
+        numbers_to_skip : int
+            Number of seeds (i.e. voxels) to skip. Useful if you want to
+            continue sampling from the same generator as in a first experiment
+            (with a fixed rng_seed).
+
+        Return
+        ------
+        random_generator : numpy random generator
+            Initialized numpy number generator.
+        indices : ndarray
+            Empty list for interface retro-compatibility
+        """
+        self.i = numbers_to_skip
+
+        return np.random.default_rng(rng_seed), []
+
+    def get_next_pos(self, random_generator: np.random.Generator,
+                     shuffled_indices, which_seed):
+        seed = self.seeds[self.i]
+        self.i += 1
+
+        return seed[0], seed[1], seed[2]
+
+    def get_next_n_pos(self, random_generator, shuffled_indices,
+                       which_seed_start, n):
+        seeds = self.seeds[self.i:self.i+n]
+        self.i += n
+
+        return seeds

scilpy/tracking/utils.py

@@ -126,6 +126,12 @@ def add_seeding_options(p):
                                     help='Number of seeds per voxel.')
     seed_sub_exclusive.add_argument('--nt', type=int,
                                     help='Total number of seeds to use.')
+    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 \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/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))

scilpy/tractanalysis/mrds_along_streamlines.py

@@ -2,7 +2,8 @@
 
 import numpy as np
 
-from scilpy.tractanalysis.grid_intersections import grid_intersections
+from scilpy.tractanalysis.voxel_boundary_intersection import\
+    subdivide_streamlines_at_voxel_faces
 
 
 def mrds_metrics_along_streamlines(sft, mrds_pdds,
@@ -79,9 +80,10 @@ def mrds_metric_sums_along_streamlines(sft, mrds_pdds, metrics,
     weight_map = np.zeros(metrics[0].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.
@@ -91,7 +93,7 @@ def mrds_metric_sums_along_streamlines(sft, mrds_pdds, metrics,
         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)