diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 42c99cb2f..d5afe5c89 100644 --- a/.github/workflows/test.yml +++ b/.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 diff --git a/docs/source/fake_files/grid_intersections.py b/docs/source/fake_files/grid_intersections.py deleted file mode 100644 index 272c3e90c..000000000 --- a/docs/source/fake_files/grid_intersections.py +++ /dev/null @@ -1,5 +0,0 @@ -# -*- coding: utf-8 -*- - - -def grid_intersections(): - pass diff --git a/docs/source/fake_files/voxel_boundary_intersection.py b/docs/source/fake_files/voxel_boundary_intersection.py new file mode 100644 index 000000000..113345ea6 --- /dev/null +++ b/docs/source/fake_files/voxel_boundary_intersection.py @@ -0,0 +1,5 @@ +# -*- coding: utf-8 -*- + + +def subdivide_streamlines_at_voxel_faces(): + pass diff --git a/docs/source/modules/scilpy.tractanalysis.rst b/docs/source/modules/scilpy.tractanalysis.rst index 26af90350..3a5ea8905 100644 --- a/docs/source/modules/scilpy.tractanalysis.rst +++ b/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: diff --git a/scilpy/connectivity/connectivity.py b/scilpy/connectivity/connectivity.py index dd6b75901..11ed51a39 100644 --- a/scilpy/connectivity/connectivity.py +++ b/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) diff --git a/scilpy/io/fetcher.py b/scilpy/io/fetcher.py index 6876d6e6e..945a1e2fb 100644 --- a/scilpy/io/fetcher.py +++ b/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 diff --git a/scilpy/io/hdf5.py b/scilpy/io/hdf5.py index e31b9dad8..aa3766ae6 100644 --- a/scilpy/io/hdf5.py +++ b/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: diff --git a/scilpy/io/streamlines.py b/scilpy/io/streamlines.py index cf26c85cc..0fc196f6a 100644 --- a/scilpy/io/streamlines.py +++ b/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 diff --git a/scilpy/tracking/seed.py b/scilpy/tracking/seed.py index 8f869548c..6fa1247a0 100644 --- a/scilpy/tracking/seed.py +++ b/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): """ diff --git a/scilpy/tracking/utils.py b/scilpy/tracking/utils.py index 15ffd2e5a..853ebd4bc 100644 --- a/scilpy/tracking/utils.py +++ b/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): diff --git a/scilpy/tractanalysis/afd_along_streamlines.py b/scilpy/tractanalysis/afd_along_streamlines.py index e5f8d6aa8..ad15ed3e2 100644 --- a/scilpy/tractanalysis/afd_along_streamlines.py +++ b/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) diff --git a/scilpy/tractanalysis/bingham_metric_along_streamlines.py b/scilpy/tractanalysis/bingham_metric_along_streamlines.py index f20fa27e9..91eaf7a5a 100644 --- a/scilpy/tractanalysis/bingham_metric_along_streamlines.py +++ b/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) diff --git a/scilpy/tractanalysis/fibertube_scoring.py b/scilpy/tractanalysis/fibertube_scoring.py index c10b39a64..53bc1b23b 100644 --- a/scilpy/tractanalysis/fibertube_scoring.py +++ b/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) diff --git a/scilpy/tractanalysis/fixel_density.py b/scilpy/tractanalysis/fixel_density.py index 3773ace84..42d5d90dc 100644 --- a/scilpy/tractanalysis/fixel_density.py +++ b/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)) diff --git a/scilpy/tractanalysis/mrds_along_streamlines.py b/scilpy/tractanalysis/mrds_along_streamlines.py index 91bdeca14..15b648726 100644 --- a/scilpy/tractanalysis/mrds_along_streamlines.py +++ b/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) diff --git a/scilpy/tractanalysis/grid_intersections.pyx b/scilpy/tractanalysis/voxel_boundary_intersection.pyx similarity index 95% rename from scilpy/tractanalysis/grid_intersections.pyx rename to scilpy/tractanalysis/voxel_boundary_intersection.pyx index 7edcc18c2..05eb63c08 100644 --- a/scilpy/tractanalysis/grid_intersections.pyx +++ b/scilpy/tractanalysis/voxel_boundary_intersection.pyx @@ -29,7 +29,23 @@ cdef struct Pointers: @cython.boundscheck(False) @cython.wraparound(False) @cython.cdivision(True) -def grid_intersections(streamlines): +def subdivide_streamlines_at_voxel_faces(streamlines): + """ + Cut streamlines segments into smaller segments such that a segment covering + more than one voxel is split into smaller segments that either end or start + at voxel boundaries. + + Parameters + ---------- + streamlines: list of ndarray + Streamlines coordinates in voxel space, corner origin. + + Returns + ------- + split_coordinates: list of ndarray + Updated streamline coordinates with added coordinate points + at voxel boundaries. + """ cdef: cnp.npy_intp nb_streamlines = len(streamlines._lengths) cnp.npy_intp at_point = 0 diff --git a/scilpy/tractograms/dps_and_dpp_management.py b/scilpy/tractograms/dps_and_dpp_management.py index bf40e4350..5390bf063 100644 --- a/scilpy/tractograms/dps_and_dpp_management.py +++ b/scilpy/tractograms/dps_and_dpp_management.py @@ -217,7 +217,7 @@ def project_dpp_to_map(sft, dpp_key, sum_lines=False, endpoints_only=False): for p in points: x, y, z = sft.streamlines[s][p, :].astype(int) # Or floor count[x, y, z] += 1 - the_map[x, y, z] += sft.data_per_point[dpp_key][s][p] + the_map[x, y, z] += np.squeeze(sft.data_per_point[dpp_key][s][p]) if not sum_lines: count = np.maximum(count, 1e-6) # Avoid division by 0 diff --git a/scilpy/tractograms/streamline_operations.py b/scilpy/tractograms/streamline_operations.py index 5249b80a8..6e244af92 100644 --- a/scilpy/tractograms/streamline_operations.py +++ b/scilpy/tractograms/streamline_operations.py @@ -408,11 +408,8 @@ def filter_streamlines_by_length(sft, min_length=0., max_length=np.inf, valid_length_ids = np.logical_and(lengths >= min_length, lengths <= max_length) filtered_sft = sft[valid_length_ids] - - if return_rejected: - rejected_sft = sft[~valid_length_ids] else: - valid_length_ids = [] + valid_length_ids = np.array([], dtype=bool) filtered_sft = sft # Return to original space @@ -420,7 +417,7 @@ def filter_streamlines_by_length(sft, min_length=0., max_length=np.inf, filtered_sft.to_space(orig_space) if return_rejected: - rejected_sft.to_space(orig_space) + rejected_sft = sft[~valid_length_ids] return filtered_sft, valid_length_ids, rejected_sft else: return filtered_sft, valid_length_ids diff --git a/scilpy/tractograms/tests/test_streamline_operations.py b/scilpy/tractograms/tests/test_streamline_operations.py index ccddcc972..d381da976 100644 --- a/scilpy/tractograms/tests/test_streamline_operations.py +++ b/scilpy/tractograms/tests/test_streamline_operations.py @@ -7,6 +7,7 @@ import pytest from dipy.io.streamline import load_tractogram from dipy.tracking.streamlinespeed import length +from dipy.io.stateful_tractogram import StatefulTractogram from scilpy import SCILPY_HOME from scilpy.io.fetcher import fetch_data, get_testing_files_dict @@ -174,6 +175,17 @@ def test_filter_streamlines_by_length(): # Test that streamlines shorter than 100 and longer than 120 were removed. assert np.all(lengths >= min_length) and np.all(lengths <= max_length) + # === 4. Return rejected streamlines with empty sft === + empty_sft = short_sft[[]] # Empty sft from short_sft (chosen arbitrarily) + filtered_sft, _, rejected = \ + filter_streamlines_by_length(empty_sft, min_length=min_length, + max_length=max_length, + return_rejected=True) + assert isinstance(filtered_sft, StatefulTractogram) + assert isinstance(rejected, StatefulTractogram) + assert len(filtered_sft) == 0 + assert len(rejected) == 0 + def test_filter_streamlines_by_total_length_per_dim(): long_sft = load_tractogram(in_long_sft, in_ref) diff --git a/scilpy/viz/backends/fury.py b/scilpy/viz/backends/fury.py index 3940dacfb..ab930cb6d 100644 --- a/scilpy/viz/backends/fury.py +++ b/scilpy/viz/backends/fury.py @@ -148,7 +148,7 @@ def set_viewport(scene, orientation, slice_index, volume_shape, aspect_ratio): Ratio between viewport's width and height. """ - scene.projection('parallel') + scene.projection(proj_type='parallel') camera = initialize_camera( orientation, slice_index, volume_shape, aspect_ratio) scene.set_camera(position=camera[CamParams.VIEW_POS], diff --git a/scripts/scil_NODDI_maps.py b/scripts/scil_NODDI_maps.py index 96be2f227..2d067c5ff 100755 --- a/scripts/scil_NODDI_maps.py +++ b/scripts/scil_NODDI_maps.py @@ -120,12 +120,13 @@ def main(): shells_centroids, indices_shells = identify_shells(bvals, args.tolerance, round_centroids=True) - nb_shells = len(shells_centroids) + non_b0_shells = shells_centroids[shells_centroids > args.tolerance] + nb_shells = len(non_b0_shells) if nb_shells <= 1: raise ValueError("Amico's NODDI works with data with more than one " "shell, but you seem to have single-shell data (we " - "found shells {}). Change tolerance if necessary." - .format(np.sort(shells_centroids))) + "found shell {}). Change tolerance if necessary." + .format(non_b0_shells[0])) logging.info('Will compute NODDI with AMICO on {} shells at found at {}.' .format(len(shells_centroids), np.sort(shells_centroids))) diff --git a/scripts/scil_fibertube_score_tractogram.py b/scripts/scil_fibertube_score_tractogram.py index 18fef9926..b8272c148 100644 --- a/scripts/scil_fibertube_score_tractogram.py +++ b/scripts/scil_fibertube_score_tractogram.py @@ -6,6 +6,12 @@ tracking, computes metrics about the quality of individual fiber reconstruction. +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. +TODO: Add the seed's segment index as dps, to allow different seeding methods +and forward_only=False. + Each streamline is associated with an "Arrival fibertube segment", which is the closest fibertube segment to its before-last coordinate. We then define the following terms: diff --git a/scripts/scil_fibertube_tracking.py b/scripts/scil_fibertube_tracking.py index f14df22ab..2ef4b48a9 100644 --- a/scripts/scil_fibertube_tracking.py +++ b/scripts/scil_fibertube_tracking.py @@ -111,10 +111,6 @@ def _build_arg_parser(): 'tracking sphere or max angle is reached.\n' 'Default: 0, i.e. do not add points following ' 'an invalid direction.') - track_g.add_argument( - '--forward_only', action='store_true', - help='If set, tracks in one direction only (forward) \n' - 'given the \ninitial seed.') track_g.add_argument( '--keep_last_out_point', action='store_true', help='If set, keep the last point (once out of the \n' @@ -248,7 +244,7 @@ def main(): args.max_invalid_nb_points, 0, args.nbr_processes, True, 'r+', rng_seed=args.rng_seed, - track_forward_only=args.forward_only, + track_forward_only=True, skip=args.skip, verbose=args.verbose, append_last_point=args.keep_last_out_point) diff --git a/scripts/scil_freewater_maps.py b/scripts/scil_freewater_maps.py index f24c3a94f..854eb30ab 100755 --- a/scripts/scil_freewater_maps.py +++ b/scripts/scil_freewater_maps.py @@ -2,7 +2,7 @@ # -*- coding: utf-8 -*- """ -Compute Free Water maps [1] using AMICO. +Compute Free Water maps [1] using the AMICO framework [2]. This script supports both single and multi-shell data. Formerly: scil_compute_freewater.py @@ -35,6 +35,9 @@ [1] Pasternak 0, Sochen N, Gur Y, Intrator N, Assaf Y. Free water elimination and mapping from diffusion mri. Magn Reson Med. 62 (3) (2009) 717-730. + [2] Daducci A, et al. Accelerated microstructure imaging + via convex optimization (AMICO) from diffusion MRI data. + Neuroimage 105 (2015) 32-44. """ diff --git a/scripts/scil_tracking_local_dev.py b/scripts/scil_tracking_local_dev.py index 889655558..d250c6f55 100755 --- a/scripts/scil_tracking_local_dev.py +++ b/scripts/scil_tracking_local_dev.py @@ -206,7 +206,7 @@ def main(): seed_res = seed_img.header.get_zooms()[:3] if args.in_custom_seeds: seeds = np.squeeze(load_matrix_in_any_format(args.in_custom_seeds)) - seed_generator = CustomSeedsDispenser(seeds, seed_res, space=our_space, + seed_generator = CustomSeedsDispenser(seeds, space=our_space, origin=our_origin) nbr_seeds = len(seeds) else: diff --git a/scripts/scil_tractogram_dpp_math.py b/scripts/scil_tractogram_dpp_math.py index 9303c1d32..5949fddb6 100755 --- a/scripts/scil_tractogram_dpp_math.py +++ b/scripts/scil_tractogram_dpp_math.py @@ -19,7 +19,8 @@ If endpoints_only and dps mode is set operation will be calculated across the data at the endpoints and stored as a single value (or array in the 4D case) -per streamline. +per streamline. If you wish to perform operations on dps values, please use +scil_tractogram_dps_math.py. Endpoint only operation: correlation: correlation calculated between arrays extracted from streamline diff --git a/scripts/scil_tractogram_dps_math.py b/scripts/scil_tractogram_dps_math.py index 30432b093..02837e07d 100755 --- a/scripts/scil_tractogram_dps_math.py +++ b/scripts/scil_tractogram_dps_math.py @@ -2,19 +2,30 @@ # -*- coding: utf-8 -*- """ -Add, extract or delete information to each streamline from a tractogram +Import, extract or delete dps (data_per_streamline) information to a tractogram file. Can be for example SIFT2 weights, processing information, bundle IDs, tracking seeds, etc. -Input and output tractograms must always be .trk to simplify filetype checks -for each operation. +This script is not the same as the dps mode of scil_tractogram_dpp_math.py, +which performs operations on dpp (data_per_point) and saves the result as dps. +Instead this script performs operations directly on dps values. + +Input and output tractograms must be .trk, unless you are using the 'import' +operation, in which case a .tck input tractogram is accepted. + +Usage examples: + > scil_tractogram_dps_math.py tractogram.trk import "bundle_ids" + --in_dps_file my_bundle_ids.txt + > scil_tractogram_dps_math.py tractogram.trk export "seeds" + --out_dps_file seeds.npy """ -import os +import nibabel as nib import argparse import logging from dipy.io.streamline import save_tractogram, load_tractogram +from scilpy.io.streamlines import load_tractogram_with_reference import numpy as np from scilpy.io.utils import (add_overwrite_arg, @@ -33,32 +44,29 @@ def _build_arg_parser(): formatter_class=argparse.RawTextHelpFormatter) p.add_argument('in_tractogram', - help='Input tractogram (.trk).') + help='Input tractogram (.trk for all operations,' + '.tck accepted for import).') p.add_argument('operation', metavar='OPERATION', - choices=['add', 'delete', 'export'], + choices=['import', 'delete', 'export'], help='The type of operation to be performed on the\n' 'tractogram\'s data_per_streamline at the given\n' 'key. Must be one of the following: [%(choices)s].\n' 'The additional arguments required for each\n' 'operation are specified under each group below.') p.add_argument('dps_key', type=str, - help='Where to find the data to be exported,\n' - 'in the tractogram.') + help='Key name used for the operation.') p.add_argument('--out_tractogram', - help='Output tractogram (.trk). Required for any mutation.') + help='Output tractogram (.trk). Required for "import" and\n' + '"delete" operations.') - add_args = p.add_argument_group('Add operation', - 'Requires the out_tractogram argument.') - add_args.add_argument('--add_dps_file', - help='File containing the data to add to\n' - 'streamlines (.txt, .npy or .mat).') + import_args = p.add_argument_group('Operation "import" mandatory options') + import_args.add_argument('--in_dps_file', + help='File containing the data to import to\n' + 'streamlines (.txt, .npy or .mat).') - _ = p.add_argument_group('Delete operation', - 'Requires the out_tractogram argument.') - - export_args = p.add_argument_group('Export operation') - export_args.add_argument('--export_dps_file', + export_args = p.add_argument_group('Operation "export" mandatory options') + export_args.add_argument('--out_dps_file', help='File in which the extracted data will be\n' 'saved (.txt or .npy).') @@ -73,18 +81,31 @@ def main(): args = parser.parse_args() logging.getLogger().setLevel(logging.getLevelName(args.verbose)) - check_tract_trk(parser, args.in_tractogram) + if args.operation == 'import': + if not nib.streamlines.is_supported(args.in_tractogram): + parser.error('Invalid input streamline file format (must be trk ' + + 'or tck): {0}'.format(args.in_tractogram)) + else: + check_tract_trk(parser, args.in_tractogram) + if args.out_tractogram: check_tract_trk(parser, args.out_tractogram) - # I/O assertions - assert_inputs_exist(parser, args.in_tractogram, args.add_dps_file) - assert_outputs_exist(parser, args, [], optional=[args.export_dps_file, + assert_inputs_exist(parser, args.in_tractogram, args.in_dps_file) + assert_outputs_exist(parser, args, [], optional=[args.out_dps_file, args.out_tractogram]) - sft = load_tractogram(args.in_tractogram, 'same') + sft = load_tractogram_with_reference(parser, args, args.in_tractogram) + + if args.operation == 'import': + if args.in_dps_file is None: + parser.error('The --in_dps_file option is required for ' + + 'the "import" operation.') + + if args.out_tractogram is None: + parser.error('The --out_tractogram option is required for ' + + 'the "import" operation.') - if args.operation == 'add': # Make sure the user is not unwillingly overwritting dps if (args.dps_key in sft.get_data_per_streamline_keys() and not args.overwrite): @@ -92,7 +113,7 @@ def main(): ' overwriting.'.format(args.dps_key)) # Load data and remove extraneous dimensions - data = np.squeeze(load_matrix_in_any_format(args.add_dps_file)) + data = np.squeeze(load_matrix_in_any_format(args.in_dps_file)) # Quick check as the built-in error from sft is not too explicit if len(sft) != data.shape[0]: @@ -105,11 +126,19 @@ def main(): save_tractogram(sft, args.out_tractogram) if args.operation == 'delete': + if args.out_tractogram is None: + parser.error('The --out_tractogram option is required for ' + + 'the "delete" operation.') + del sft.data_per_streamline[args.dps_key] save_tractogram(sft, args.out_tractogram) if args.operation == 'export': + if args.out_dps_file is None: + parser.error('The --out_dps_file option is required for ' + + 'the "export" operation.') + # Extract data and reshape if args.dps_key not in sft.data_per_streamline.keys(): raise ValueError('Data does not have any data_per_streamline' @@ -117,7 +146,7 @@ def main(): .format(args.dps_key)) data = np.squeeze(sft.data_per_streamline[args.dps_key]) - save_matrix_in_any_format(args.export_dps_file, data) + save_matrix_in_any_format(args.out_dps_file, data) if __name__ == '__main__': diff --git a/scripts/scil_tractogram_filter_by_length.py b/scripts/scil_tractogram_filter_by_length.py index 1aed05283..9b3ba5871 100755 --- a/scripts/scil_tractogram_filter_by_length.py +++ b/scripts/scil_tractogram_filter_by_length.py @@ -48,8 +48,9 @@ def _build_arg_parser(): help='Do not write file if there is no streamline.') p.add_argument('--display_counts', action='store_true', help='Print streamline count before and after filtering') - p.add_argument('--save_rejected', action='store_true', - help='Save rejected streamlines to output tractogram.') + p.add_argument('--out_rejected', + help='If specified, save rejected streamlines to this ' + 'file.') add_json_args(p) add_reference_arg(p) add_verbose_arg(p) @@ -75,8 +76,17 @@ def main(): sft = load_tractogram_with_reference(parser, args, args.in_tractogram) # Processing - new_sft, _, outliers_sft = filter_streamlines_by_length( - sft, args.minL, args.maxL, return_rejected=True) + return_rejected = args.out_rejected is not None + + # Returns (new_sft, _, [rejected_sft]) in filtered_result + filtered_result = filter_streamlines_by_length( + sft, args.minL, args.maxL, return_rejected=return_rejected) + + new_sft = filtered_result[0] + + if return_rejected: + rejected_sft = filtered_result[2] + save_tractogram(rejected_sft, args.out_rejected, args.no_empty) if args.display_counts: sc_bf = len(sft.streamlines) diff --git a/scripts/tests/test_NODDI_maps.py b/scripts/tests/test_NODDI_maps.py index 105bb72e6..ed7e66b9f 100644 --- a/scripts/tests/test_NODDI_maps.py +++ b/scripts/tests/test_NODDI_maps.py @@ -8,7 +8,8 @@ from scilpy.io.fetcher import fetch_data, get_testing_files_dict # If they already exist, this only takes 5 seconds (check md5sum) -fetch_data(get_testing_files_dict(), keys=['commit_amico.zip']) +fetch_data(get_testing_files_dict(), keys=['commit_amico.zip', + 'processing.zip']) tmp_dir = tempfile.TemporaryDirectory() @@ -32,5 +33,22 @@ def test_execution_commit_amico(script_runner, monkeypatch): '--out_dir', 'noddi', '--tol', '30', '--para_diff', '0.0017', '--iso_diff', '0.003', '--lambda1', '0.5', '--lambda2', '0.001', - '--processes', '1') + '--processes', '1', '-f') assert ret.success + + +def test_single_shell_fail(script_runner, monkeypatch): + monkeypatch.chdir(os.path.expanduser(tmp_dir.name)) + in_dwi = os.path.join(SCILPY_HOME, 'processing', + 'dwi_crop_1000.nii.gz') + in_bval = os.path.join(SCILPY_HOME, 'processing', + '1000.bval') + in_bvec = os.path.join(SCILPY_HOME, 'processing', + '1000.bvec') + ret = script_runner.run('scil_NODDI_maps.py', in_dwi, + in_bval, in_bvec, + '--out_dir', 'noddi', '--tol', '30', + '--para_diff', '0.0017', '--iso_diff', '0.003', + '--lambda1', '0.5', '--lambda2', '0.001', + '--processes', '1', '-f') + assert not ret.success diff --git a/scripts/tests/test_tractogram_dps_math.py b/scripts/tests/test_tractogram_dps_math.py index 6a4e5e007..03b758de2 100644 --- a/scripts/tests/test_tractogram_dps_math.py +++ b/scripts/tests/test_tractogram_dps_math.py @@ -21,7 +21,7 @@ def test_help_option(script_runner): assert ret.success -def test_execution_dps_math_add(script_runner, monkeypatch): +def test_execution_dps_math_import(script_runner, monkeypatch): monkeypatch.chdir(os.path.expanduser(tmp_dir.name)) in_bundle = os.path.join(SCILPY_HOME, 'filtering', 'bundle_4.trk') @@ -30,14 +30,15 @@ def test_execution_dps_math_add(script_runner, monkeypatch): outname = 'out.trk' np.save(filename, np.arange(len(sft))) ret = script_runner.run('scil_tractogram_dps_math.py', - in_bundle, 'add', 'key', - '--add_dps_file', filename, + in_bundle, 'import', 'key', + '--in_dps_file', filename, '--out_tractogram', outname, '-f') assert ret.success -def test_execution_dps_math_add_with_missing_vals(script_runner, monkeypatch): +def test_execution_dps_math_import_with_missing_vals(script_runner, + monkeypatch): monkeypatch.chdir(os.path.expanduser(tmp_dir.name)) in_bundle = os.path.join(SCILPY_HOME, 'filtering', 'bundle_4.trk') @@ -46,14 +47,15 @@ def test_execution_dps_math_add_with_missing_vals(script_runner, monkeypatch): outname = 'out.trk' np.save(filename, np.arange(len(sft) - 10)) ret = script_runner.run('scil_tractogram_dps_math.py', - in_bundle, 'add', 'key', - '--add_dps_file', filename, + in_bundle, 'import', 'key', + '--in_dps_file', filename, '--out_tractogram', outname, '-f') assert ret.stderr -def test_execution_dps_math_add_with_existing_key(script_runner, monkeypatch): +def test_execution_dps_math_import_with_existing_key(script_runner, + monkeypatch): monkeypatch.chdir(os.path.expanduser(tmp_dir.name)) in_bundle = os.path.join(SCILPY_HOME, 'filtering', 'bundle_4.trk') @@ -63,19 +65,19 @@ def test_execution_dps_math_add_with_existing_key(script_runner, monkeypatch): outname2 = 'out_2.trk' np.save(filename, np.arange(len(sft))) ret = script_runner.run('scil_tractogram_dps_math.py', - in_bundle, 'add', 'key', - '--add_dps_file', filename, + in_bundle, 'import', 'key', + '--in_dps_file', filename, '--out_tractogram', outname, '-f') assert ret.success ret = script_runner.run('scil_tractogram_dps_math.py', - outname, 'add', 'key', - '--add_dps_file', filename, + outname, 'import', 'key', + '--in_dps_file', filename, '--out_tractogram', outname2,) assert not ret.success -def test_execution_dps_math_tck(script_runner, monkeypatch): +def test_execution_dps_math_tck_output(script_runner, monkeypatch): monkeypatch.chdir(os.path.expanduser(tmp_dir.name)) in_bundle = os.path.join(SCILPY_HOME, 'filtering', 'bundle_4.trk') @@ -84,8 +86,8 @@ def test_execution_dps_math_tck(script_runner, monkeypatch): outname = 'out.tck' np.save(filename, np.arange(len(sft))) ret = script_runner.run('scil_tractogram_dps_math.py', - in_bundle, 'add', 'key', - '--add_dps_file', filename, + in_bundle, 'import', 'key', + '--in_dps_file', filename, '--out_tractogram', outname, '-f') assert not ret.success @@ -134,7 +136,7 @@ def test_execution_dps_math_export(script_runner, monkeypatch): filename = 'out.txt' ret = script_runner.run('scil_tractogram_dps_math.py', in_bundle, 'export', 'key', - '--export_dps_file', filename, + '--out_dps_file', filename, '-f') assert ret.success @@ -146,6 +148,6 @@ def test_execution_dps_math_export_no_key(script_runner, monkeypatch): filename = 'out.txt' ret = script_runner.run('scil_tractogram_dps_math.py', in_bundle, 'export', 'key', - '--export_dps_file', filename, + '--out_dps_file', filename, '-f') assert not ret.success diff --git a/scripts/tests/test_tractogram_filter_by_length.py b/scripts/tests/test_tractogram_filter_by_length.py index fdac121ef..d5f4eda7d 100644 --- a/scripts/tests/test_tractogram_filter_by_length.py +++ b/scripts/tests/test_tractogram_filter_by_length.py @@ -6,6 +6,7 @@ from scilpy import SCILPY_HOME from scilpy.io.fetcher import fetch_data, get_testing_files_dict +from scilpy.io.streamlines import load_tractogram # If they already exist, this only takes 5 seconds (check md5sum) fetch_data(get_testing_files_dict(), keys=['filtering.zip']) @@ -20,9 +21,55 @@ def test_help_option(script_runner): def test_execution_filtering(script_runner, monkeypatch): monkeypatch.chdir(os.path.expanduser(tmp_dir.name)) + # Effectively, this doesn't filter anything, since bundle_4.trk has + # all streamlines with lengths of 130mm. This is just to test the + # script execution. in_bundle = os.path.join(SCILPY_HOME, 'filtering', 'bundle_4.trk') ret = script_runner.run('scil_tractogram_filter_by_length.py', in_bundle, 'bundle_4_filtered.trk', '--minL', '125', '--maxL', '130') + + sft = load_tractogram('bundle_4_filtered.trk', 'same') + assert len(sft) == 52 + + assert ret.success + + +def test_rejected_filtering(script_runner, monkeypatch): + monkeypatch.chdir(os.path.expanduser(tmp_dir.name)) + in_bundle = os.path.join(SCILPY_HOME, 'filtering', + 'bundle_all_1mm.trk') + ret = script_runner.run('scil_tractogram_filter_by_length.py', + in_bundle, 'bundle_all_1mm_filtered.trk', + '--minL', '125', '--maxL', '130', + '--out_rejected', 'bundle_all_1mm_rejected.trk') + assert ret.success + assert os.path.exists('bundle_all_1mm_rejected.trk') + assert os.path.exists('bundle_all_1mm_rejected.trk') + + sft = load_tractogram('bundle_all_1mm_filtered.trk', 'same') + rejected_sft = load_tractogram('bundle_all_1mm_rejected.trk', 'same') + + assert len(sft) == 266 + assert len(rejected_sft) == 2824 + + +def test_rejected_filtering_no_rejection(script_runner, monkeypatch): + monkeypatch.chdir(os.path.expanduser(tmp_dir.name)) + in_bundle = os.path.join(SCILPY_HOME, 'filtering', + 'bundle_4.trk') + ret = script_runner.run('scil_tractogram_filter_by_length.py', + in_bundle, 'bundle_4_filtered_no_rejection.trk', + '--minL', '125', '--maxL', '130', + '--out_rejected', 'bundle_4_rejected.trk') assert ret.success + + # File should be created even though there are no rejected streamlines + assert os.path.exists('bundle_4_rejected.trk') + + sft = load_tractogram('bundle_4_filtered_no_rejection.trk', 'same') + rejected_sft = load_tractogram('bundle_4_rejected.trk', 'same') + + assert len(sft) == 52 + assert len(rejected_sft) == 0 diff --git a/scripts/tests/test_viz_volume_screenshot.py b/scripts/tests/test_viz_volume_screenshot.py index 1e4db3509..12cb79aab 100644 --- a/scripts/tests/test_viz_volume_screenshot.py +++ b/scripts/tests/test_viz_volume_screenshot.py @@ -12,7 +12,8 @@ tmp_dir = tempfile.TemporaryDirectory() -def test_screenshot(script_runner): +def test_screenshot(script_runner, monkeypatch): + monkeypatch.chdir(os.path.expanduser(tmp_dir.name)) in_fa = os.path.join(SCILPY_HOME, 'bst', 'fa.nii.gz') ret = script_runner.run("scil_viz_volume_screenshot.py", in_fa, 'fa.png') diff --git a/setup.py b/setup.py index 07a2a7ad1..0a03ecc40 100644 --- a/setup.py +++ b/setup.py @@ -20,11 +20,14 @@ def get_extensions(): ['scilpy/tractograms/uncompress.pyx']) quick_tools = Extension('scilpy.tractanalysis.quick_tools', ['scilpy/tractanalysis/quick_tools.pyx']) - grid_intersections = Extension('scilpy.tractanalysis.grid_intersections', - ['scilpy/tractanalysis/grid_intersections.pyx']) - streamlines_metrics = Extension('scilpy.tractanalysis.streamlines_metrics', - ['scilpy/tractanalysis/streamlines_metrics.pyx']) - return [uncompress, quick_tools, grid_intersections, streamlines_metrics] + voxel_boundary_intersection =\ + Extension('scilpy.tractanalysis.voxel_boundary_intersection', + ['scilpy/tractanalysis/voxel_boundary_intersection.pyx']) + streamlines_metrics =\ + Extension('scilpy.tractanalysis.streamlines_metrics', + ['scilpy/tractanalysis/streamlines_metrics.pyx']) + return [uncompress, quick_tools, + voxel_boundary_intersection, streamlines_metrics] class CustomBuildExtCommand(build_ext):