From e771d9802336476d32079e46ff9a31cf2a053fab Mon Sep 17 00:00:00 2001 From: frheault Date: Thu, 16 May 2024 14:56:37 -0400 Subject: [PATCH 01/12] Added docstring --- scilpy/tractograms/tractogram_operations.py | 364 +++++++++++++++++++- scilpy/utils/spatial.py | 19 +- scripts/scil_tractogram_alter.py | 139 ++++++++ scripts/tests/test_tractogram_alter.py | 6 + 4 files changed, 524 insertions(+), 4 deletions(-) create mode 100644 scripts/scil_tractogram_alter.py create mode 100644 scripts/tests/test_tractogram_alter.py diff --git a/scilpy/tractograms/tractogram_operations.py b/scilpy/tractograms/tractogram_operations.py index 81c4a6567..794a032b0 100644 --- a/scilpy/tractograms/tractogram_operations.py +++ b/scilpy/tractograms/tractogram_operations.py @@ -12,7 +12,8 @@ import logging import random -from dipy.io.stateful_tractogram import StatefulTractogram, Space +from dipy.io.stateful_tractogram import set_sft_logger_level, \ + StatefulTractogram, Space from dipy.io.utils import get_reference_info, is_header_compatible from dipy.segment.clustering import qbx_and_merge from dipy.tracking.streamline import transform_streamlines @@ -24,9 +25,14 @@ from scipy.ndimage import map_coordinates from scipy.spatial import cKDTree +from scilpy.tractanalysis.bundle_operations import uniformize_bundle_sft +from scilpy.tractanalysis.streamlines_metrics import compute_tract_counts_map from scilpy.tractograms.streamline_operations import smooth_line_gaussian, \ resample_streamlines_step_size, parallel_transport_streamline, \ - cut_invalid_streamlines + cut_invalid_streamlines, compress_sft, cut_invalid_streamlines, \ + remove_overlapping_points_streamlines, remove_single_point_streamlines +from scilpy.tractograms.streamline_and_mask_operations import cut_outside_of_mask_streamlines +from scilpy.utils.spatial import generate_rotation_matrix MIN_NB_POINTS = 10 KEY_INDEX = np.concatenate((range(5), range(-1, -6, -1))) @@ -693,7 +699,8 @@ def upsample_tractogram(sft, nb, point_wise_std=None, tube_radius=None, polynomial = Polynomial(poly_coeffs[::-1]) noise_factor = polynomial(x) - vec = s - new_s + tmp_s = s - new_s + vec = tmp_s if not np.any(tmp_s == 0) else np.ones(tmp_s.shape) vec /= np.linalg.norm(vec, axis=0) new_s += vec * np.expand_dims(noise_factor, axis=1) @@ -879,6 +886,357 @@ def split_sft_randomly_per_cluster(orig_sft, chunk_sizes, seed, thresholds): return final_sfts +def subsample_streamlines_alter(sft, min_dice=0.90, epsilon=0.01, + baseline_sft=None): + """ + Function to subsample streamlines based on a dice similarity metric. + The function will keep removing streamlines until the dice similarity + between the original and the subsampled tractogram is close to min_dice. + + Parameters + ---------- + sft: StatefulTractogram + The tractogram to subsample. + min_dice: float + The minimum dice similarity to reach before stopping the subsampling. + epsilon: float + Stopping criteria for convergence. The maximum difference between the + dice similarity and min_dice. + baseline_sft: StatefulTractogram + The tractogram to use as a reference for the dice similarity. If None, + the original tractogram will be used. + + Returns + ------- + new_sft: StatefulTractogram + The tractogram with a subset of streamlines in the same space as the + input tractogram. + """ + # Import in function to avoid circular import error + from scilpy.tractanalysis.reproducibility_measures import compute_dice_voxel + set_sft_logger_level(logging.ERROR) + space = sft.space + origin = sft.origin + + sft.to_vox() + sft.to_corner() + if baseline_sft is None: + original_density_map = compute_tract_counts_map(sft.streamlines, + sft.dimensions) + else: + baseline_sft.to_vox() + baseline_sft.to_corner() + original_density_map = compute_tract_counts_map(baseline_sft.streamlines, + sft.dimensions) + dice = 1.0 + init_pick_min = 0 + init_pick_max = len(sft) + previous_to_pick = None + while dice > min_dice or np.abs(dice - min_dice) > epsilon: + to_pick = init_pick_min + (init_pick_max - init_pick_min) // 2 + if to_pick == previous_to_pick: + logging.warning('No more streamlines to pick, not converging.') + break + previous_to_pick = to_pick + + indices = np.random.choice(len(sft), to_pick, replace=False) + streamlines = sft.streamlines[indices] + curr_density_map = compute_tract_counts_map(streamlines, + sft.dimensions) + dice, _ = compute_dice_voxel(original_density_map, curr_density_map) + logging.debug(f'Subsampled {to_pick} streamlines, dice: {dice}') + + if dice < min_dice: + init_pick_min = to_pick + else: + init_pick_max = to_pick + + new_sft = StatefulTractogram.from_sft(streamlines, sft) + new_sft.to_space(space) + new_sft.to_origin(origin) + return new_sft + + +def cut_streamlines_alter(sft, min_dice=0.90, epsilon=0.01, from_start=True): + """ + Cut streamlines based on a dice similarity metric. + The function will keep removing points from the streamlines until the dice + similarity between the original and the cut tractogram is close to min_dice. + + Parameters + ---------- + sft: StatefulTractogram + The tractogram to cut. + min_dice: float + The minimum dice similarity to reach before stopping the cutting. + epsilon: float + Stopping criteria for convergence. The maximum difference between the + dice similarity and min_dice. + from_start: bool + If True, cut from the start of the streamlines. If False, cut from the + end of the streamlines. + + Returns + ------- + new_sft: StatefulTractogram + The tractogram with cut streamlines in the same space as the input + tractogram. + """ + # Import in function to avoid circular import error + from scilpy.tractanalysis.reproducibility_measures import compute_dice_voxel + set_sft_logger_level(logging.ERROR) + space = sft.space + origin = sft.origin + + # Uniformize endpoints to cut consistently from one end only + uniformize_bundle_sft(sft, swap=not from_start) + sft = resample_streamlines_step_size(sft, 0.5) + sft.to_vox() + sft.to_corner() + original_density_map = compute_tract_counts_map(sft.streamlines, + sft.dimensions) + + # Initialize the dice value and the cut percentage for dichotomic search + dice = 1.0 + init_cut_min = 0 + init_cut_max = 1.0 + previous_to_pick = None + while dice > min_dice or np.abs(dice - min_dice) > epsilon: + to_pick = init_cut_min + (init_cut_max - init_cut_min) / 2 + if to_pick == previous_to_pick: + logging.warning('No more points to pick, not converging.') + break + previous_to_pick = to_pick + + streamlines = [] + for streamline in sft.streamlines: + pos_to_pick = int(len(streamline) * to_pick) + streamline = streamline[:pos_to_pick] + streamlines.append(streamline) + curr_density_map = compute_tract_counts_map(streamlines, + sft.dimensions) + dice, _ = compute_dice_voxel(original_density_map, curr_density_map) + logging.debug(f'Cut {to_pick * 100}% of the streamlines, dice: {dice}') + + if dice < min_dice: + init_cut_min = to_pick + else: + init_cut_max = to_pick + + new_sft = StatefulTractogram.from_sft(streamlines, sft) + new_sft.to_space(space) + new_sft.to_origin(origin) + return compress_sft(new_sft) + + +def replace_streamlines_alter(sft, min_dice=0.90, epsilon=0.01): + """ + Replace streamlines based on a dice similarity metric. + The function will upsamples the streamlines (with parallel transport), + then downsample them until the dice similarity is close to min_dice. + This effectively replaces the streamlines with new ones. + + Parameters + ---------- + sft: StatefulTractogram + The tractogram to replace streamlines from. + min_dice: float + The minimum dice similarity to reach before stopping the replacement. + epsilon: float + Stopping criteria for convergence. The maximum difference between the + dice similarity and min_dice. + + Returns + ------- + new_sft: StatefulTractogram + The tractogram with replaced streamlines in the same space as the input + tractogram. + """ + # Import in function to avoid circular import error + from scilpy.tractanalysis.reproducibility_measures import compute_dice_voxel + set_sft_logger_level(logging.ERROR) + + logging.debug('Upsampling the streamlines by a factor 2x to then ' + 'downsample.') + upsampled_sft = upsample_tractogram(sft, len(sft) * 2, point_wise_std=0.5, + tube_radius=1.0, gaussian=None, + error_rate=0.1, seed=1234) + return subsample_streamlines_alter(upsampled_sft, min_dice, epsilon, + baseline_sft=sft) + + +def trim_streamlines_alter(sft, min_dice=0.90, epsilon=0.01): + # Import in function to avoid circular import error + from scilpy.tractanalysis.reproducibility_measures import compute_dice_voxel + set_sft_logger_level(logging.ERROR) + space = sft.space + origin = sft.origin + + sft.to_vox() + sft.to_corner() + original_density_map = compute_tract_counts_map(sft.streamlines, + sft.dimensions).astype(np.uint64) + thr_density = [1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024] + thr_pos = 0 + voxels_to_remove = np.where( + (original_density_map <= thr_density[thr_pos]) & + (original_density_map > 0)) + + # Initialize the dice value and the number of voxels to pick + dice = 1.0 + previous_dice = 0.0 + init_trim_min = 0 + init_trim_max = np.count_nonzero(voxels_to_remove[0]) + previous_to_pick = None + + while dice > min_dice or np.abs(dice - previous_dice) > epsilon: + to_pick = init_trim_min + (init_trim_max - init_trim_min) // 2 + if to_pick == previous_to_pick or \ + np.abs(dice - previous_dice) < epsilon: + # If too few voxels are picked, increase the threshold + # and reinitialize the picking + + if np.abs(dice - min_dice) > epsilon and \ + thr_pos < len(thr_density) - 1: + thr_pos += 1 + logging.debug(f'Increasing threshold density to ' + f'{thr_density[thr_pos]}.') + + voxels_to_remove = np.where( + (original_density_map <= thr_density[thr_pos]) & + (original_density_map > 0)) + init_trim_min = 0 + init_trim_max = np.count_nonzero(voxels_to_remove[0]) + dice = 1.0 + previous_dice = 0.0 + previous_to_pick = None + continue + else: + break + previous_to_pick = to_pick + + voxel_to_remove = np.where( + (original_density_map <= thr_density[thr_pos]) & + (original_density_map > 0)) + indices = np.random.choice(np.count_nonzero(voxel_to_remove[0]), + to_pick, replace=False) + voxel_to_remove = tuple(np.array(voxel_to_remove).T[indices].T) + mask = original_density_map.copy() + mask[voxel_to_remove] = 0 + + # set logger level to ERROR to avoid logging from cut_outside_of_mask + log_level = logging.getLogger().getEffectiveLevel() + logging.getLogger().setLevel(logging.ERROR) + new_sft = cut_outside_of_mask_streamlines(sft, mask, min_len=10) + # reset logger level + logging.getLogger().setLevel(log_level) + + curr_density_map = compute_tract_counts_map(new_sft.streamlines, + sft.dimensions) + previous_dice = dice + dice, _ = compute_dice_voxel(original_density_map, curr_density_map) + logging.debug(f'Trimmed {to_pick} voxels at density ' + f'{thr_density[thr_pos]}, dice: {dice}') + + if dice < min_dice: + init_trim_max = to_pick + else: + init_trim_min = to_pick + + new_sft.to_space(space) + new_sft.to_origin(origin) + return new_sft + + +def transform_streamlines_alter(sft, min_dice=0.90, epsilon=0.01): + """ + The function will apply random rotations to the streamlines until the dice + similarity between the original and the transformed tractogram is close to + min_dice. + + Start with a large XYZ rotation, then reduce the rotation step by half one + axis at a time until the dice similarity is close to min_dice. + + Parameters + ---------- + sft: StatefulTractogram + The tractogram to transform. + min_dice: float + The minimum dice similarity to reach before stopping the transformation. + epsilon: float + Stopping criteria for convergence. The maximum difference between the + dice similarity and min_dice. + + Returns + ------- + new_sft: StatefulTractogram + The tractogram with transformed streamlines in the same space as the + input tractogram. + """ + # Import in function to avoid circular import error + from scilpy.tractanalysis.reproducibility_measures import compute_dice_voxel + set_sft_logger_level(logging.ERROR) + space = sft.space + origin = sft.origin + + sft.to_vox() + sft.to_corner() + original_density_map = compute_tract_counts_map(sft.streamlines, + sft.dimensions) + + # Initialize the dice value and angles to pick + dice = 1.0 + angle_min = [0.0, 0.0, 0.0] + angle_max = [0.1, 0.1, 0.1] + previous_dice = None + last_pick = np.array([0.0, 0.0, 0.0]) + rand_val = np.random.rand(3) * angle_max[0] + axis_choices = np.random.choice(3, 3, replace=False) + axis = 0 + while dice > min_dice or np.abs(dice - min_dice) > epsilon: + init_angle_min = angle_min[axis] + init_angle_max = angle_max[axis] + to_pick = init_angle_min + (init_angle_max - init_angle_min) / 2 + + # Generate a 4x4 matrix from random euler angles + rand_val = np.array(angle_max) + rand_val[axis] = to_pick + + angles = rand_val * 2 * np.pi + rot_mat = generate_rotation_matrix(angles) + streamlines = transform_streamlines(sft.streamlines, rot_mat) + + # Remove invalid streamlines to avoid numerical issues + curr_sft = StatefulTractogram.from_sft(streamlines, sft) + curr_sft, _ = cut_invalid_streamlines(curr_sft) + curr_sft = remove_single_point_streamlines(curr_sft) + curr_sft = remove_overlapping_points_streamlines(curr_sft) + + curr_density_map = compute_tract_counts_map(curr_sft.streamlines, + sft.dimensions) + dice, _ = compute_dice_voxel(original_density_map, curr_density_map) + logging.debug(f'Transformed {to_pick*360} degree on axis {axis}, ' + f'dice: {dice}') + last_pick[axis] = to_pick + + if dice < min_dice: + angle_max[axis] = to_pick + else: + angle_min[axis] = to_pick + + if (previous_dice is not None) \ + and np.abs(dice - previous_dice) < epsilon / 2: + logging.debug('Not converging, switching axis.\n') + axis_choices = np.roll(axis_choices, 1) + axis = axis_choices[0] + previous_dice = dice + + logging.debug(f'\nFinal angles: {last_pick * 360} at dice: {dice}') + curr_sft.to_space(space) + curr_sft.to_origin(origin) + return curr_sft + + OPERATIONS = { 'difference_robust': difference_robust, 'intersection_robust': intersection_robust, diff --git a/scilpy/utils/spatial.py b/scilpy/utils/spatial.py index 647316632..505ae56e8 100644 --- a/scilpy/utils/spatial.py +++ b/scilpy/utils/spatial.py @@ -73,7 +73,7 @@ def get_coordinate_name(axis_index, affine=np.eye(4)): ------- coordinate_name : str Name of the coordinate suffixed with sign (see RAS_AXES_COORDINATES). - """ # noqa + """ # noqa _ix, _sgn = _any2ras_index(axis_index, affine) return RAS_AXES_COORDINATES[_ix] + _sgn @@ -209,6 +209,7 @@ class WorldBoundingBox(object): voxel_size: np.ndarray Voxel size of the bounding box. """ + def __init__(self, minimums, maximums, voxel_size): self.minimums = minimums self.maximums = maximums @@ -239,3 +240,19 @@ def world_to_voxel(coord, affine): vox_coord = np.dot(iaffine, normalized_coord) vox_coord = np.round(vox_coord).astype(int) return vox_coord[0:3] + + +def generate_rotation_matrix(angles, translation=None): + rotation_matrix = np.eye(4) + x_rot = np.array([[1, 0, 0], + [0, np.cos(angles[0]), -np.sin(angles[0])], + [0, np.sin(angles[0]), np.cos(angles[0])]]) + y_rot = np.array([[np.cos(angles[1]), 0, np.sin(angles[1])], + [0, 1, 0], + [-np.sin(angles[1]), 0, np.cos(angles[1])]]) + z_rot = np.array([[np.cos(angles[2]), -np.sin(angles[2]), 0], + [np.sin(angles[2]), np.cos(angles[2]), 0], + [0, 0, 1]]) + rotation_matrix[:3, :3] = np.dot(np.dot(x_rot, y_rot), z_rot) + rotation_matrix[:3, 3] = translation if translation is not None else 0 + return rotation_matrix diff --git a/scripts/scil_tractogram_alter.py b/scripts/scil_tractogram_alter.py new file mode 100644 index 000000000..1de8f3720 --- /dev/null +++ b/scripts/scil_tractogram_alter.py @@ -0,0 +1,139 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + + +""" +This script is used to alter a bundle to reach specific minimum dice +coefficient from the original bundle. The script will subsample, trim, cut, +replace or transform the streamlines until the minimum dice is reached. +(cannot be combined in one run, use the script multiple times if needed) + +All operation use a dichotomic search to find the parameter that gets as close +as possible to the specified minimum dice coefficient (with an epsilon for +convergence). + +The subsample operation will remove streamlines until the minimum dice is +reached. This affect the whole bundle. + +The trim operation will used the lowest density voxels (starting at 1) to +remove points from the streamlines until the minimum dice is reached. +This typically affect the edge of the bundle. + +The cut operation will remove points from the start (or end) streamlines until +the minimum dice is reached. This affect one end of the bundle. + +The replace operation will upsample the tractogram (generate new streamlines +with noise) and then subsample the tractogram. This effectively replace +streamlines with similar ones until the minimum dice is reached. +This affect the whole bundle. + +The transform operation will apply random transformations to the streamlines +until the minimum dice is reached. This affect the whole bundle. +""" + +import argparse +import logging +import os + +from dipy.io.streamline import load_tractogram, save_tractogram +import numpy as np + +from scilpy.io.streamlines import load_tractogram_with_reference +from scilpy.io.utils import (add_overwrite_arg, + add_reference_arg, + add_verbose_arg, + assert_inputs_exist, + assert_outputs_exist) +from scilpy.tractograms.streamline_operations import remove_overlapping_points_streamlines, \ + remove_single_point_streamlines, cut_invalid_streamlines +from scilpy.tractograms.tractogram_operations import transform_streamlines_alter, \ + trim_streamlines_alter, cut_streamlines_alter, subsample_streamlines_alter, \ + replace_streamlines_alter + + + +def buildArgsParser(): + p = argparse.ArgumentParser(description=__doc__, + formatter_class = argparse.RawTextHelpFormatter) + p.add_argument('in_bundle', + help='Input bundle filename (.trk, .tck).') + p.add_argument('out_bundle', + help='Output bundle filename (.trk, .tck).') + + p.add_argument('--min_dice', type=float, default=0.90, + help='Minimum dice to reach.') + p.add_argument('--epsilon', type=float, default=0.01, + help='Epsilon for the convergence.') + + g = p.add_argument_group(title='Alteration options') + g1 = g.add_mutually_exclusive_group(required=True,) + g1.add_argument('--subsample', action='store_true', + help='Pick a subset of streamlines.') + g1.add_argument('--trim', action='store_true', + help='Trim streamlines using low density voxel.') + g1.add_argument('--cut', action='store_true', + help='Cut streamlines from endpoints.') + g1.add_argument('--replace', action='store_true', + help='Replace streamlines with similar ones.') + g1.add_argument('--transform', action='store_true', + help='Transform streamlines using a random linear ' + 'transformation.') + g.add_argument('--from_end', action='store_true', + help='Cut from the other end of streamlines.\n' + 'Only available with --cut.') + + p.add_argument('--seed', '-s', type=int, default=None, + help='Seed for RNG. Default based on --min_dice.') + add_overwrite_arg(p) + add_verbose_arg(p) + add_reference_arg(p) + + return p + + +def main(): + parser = buildArgsParser() + args = parser.parse_args() + + assert_inputs_exist(parser, args.in_bundle, args.reference) + assert_outputs_exist(parser, args, [args.out_bundle]) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) + + if args.from_end and args.cut is None: + parser.error('The --from_end option is only available with --cut.') + + if args.seed is None: + np.random.seed(int(args.min_dice * 1000)) + else: + np.random.seed(args.seed) + + sft = load_tractogram_with_reference(parser, args, args.in_bundle) + + # Alter the streamlines with only one operation at a time + if args.subsample: + altered_sft = subsample_streamlines_alter(sft, args.min_dice, + epsilon=args.epsilon) + elif args.trim: + altered_sft = trim_streamlines_alter(sft, args.min_dice, + epsilon=args.epsilon) + elif args.cut: + altered_sft = cut_streamlines_alter(sft, args.min_dice, + epsilon=args.epsilon, + from_start=not args.from_end) + elif args.replace: + altered_sft = replace_streamlines_alter(sft, args.min_dice, + epsilon=args.epsilon) + elif args.transform: + altered_sft = transform_streamlines_alter(sft, args.min_dice, + epsilon=args.epsilon) + + # Some operations could have generated invalid streamlines + altered_sft, _ = cut_invalid_streamlines(altered_sft) + altered_sft = remove_single_point_streamlines(altered_sft) + altered_sft = remove_overlapping_points_streamlines(altered_sft) + + save_tractogram(altered_sft, args.out_bundle) + + +if __name__ == "__main__": + main() diff --git a/scripts/tests/test_tractogram_alter.py b/scripts/tests/test_tractogram_alter.py new file mode 100644 index 000000000..5c057d800 --- /dev/null +++ b/scripts/tests/test_tractogram_alter.py @@ -0,0 +1,6 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +def test_help_option(script_runner): + ret = script_runner.run('scil_tractogram_alter.py', '--help') + assert ret.success From 799b9c1b6cf8a379b99a247fc527a60ad958117e Mon Sep 17 00:00:00 2001 From: frheault Date: Thu, 16 May 2024 15:38:36 -0400 Subject: [PATCH 02/12] Fix test after /0 --- scilpy/tractograms/tests/test_tractogram_operations.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/scilpy/tractograms/tests/test_tractogram_operations.py b/scilpy/tractograms/tests/test_tractogram_operations.py index 10868f3ac..928ff4f43 100644 --- a/scilpy/tractograms/tests/test_tractogram_operations.py +++ b/scilpy/tractograms/tests/test_tractogram_operations.py @@ -166,15 +166,15 @@ def test_combining_sft(): def test_upsample_tractogram(): new_sft = upsample_tractogram(sft, 1000, 0.5, 5, False, 0.1, 0) - first_chunk = [[112.64021, 35.409477, 59.42175], - [109.09777, 35.287857, 61.845505], - [110.41855, 37.077374, 56.930523]] + first_chunk = [[112.676895, 35.446144, 59.454517], + [109.087395, 35.27597, 61.876167], + [110.4185, 37.0766, 56.93045]] last_chunk = [[110.40285, 51.036686, 62.419273], [109.698586, 48.330017, 64.50656], [113.04737, 45.89119, 64.778534]] assert len(new_sft) == 1000 - assert len(new_sft.streamlines._data) == 8404 + assert len(new_sft.streamlines._data) == 8402 assert np.allclose(first_chunk, new_sft.streamlines._data[0:30:10]) assert np.allclose(last_chunk, new_sft.streamlines._data[-1:-31:-10]) From 90113f1c6ccf8e96adf70db6737ff34bb16dce73 Mon Sep 17 00:00:00 2001 From: frheault Date: Wed, 22 May 2024 12:08:14 -0400 Subject: [PATCH 03/12] Add shuffle options --- scilpy/tractograms/tractogram_operations.py | 57 ++++++++++++++++++++- scripts/scil_tractogram_alter.py | 11 +++- 2 files changed, 65 insertions(+), 3 deletions(-) diff --git a/scilpy/tractograms/tractogram_operations.py b/scilpy/tractograms/tractogram_operations.py index 794a032b0..8d4f2591a 100644 --- a/scilpy/tractograms/tractogram_operations.py +++ b/scilpy/tractograms/tractogram_operations.py @@ -39,6 +39,20 @@ def shuffle_streamlines(sft, rng_seed=None): + """ + Shuffle the streamlines of a tractogram. + + Parameters + ---------- + sft: StatefulTractogram + The tractogram to shuffle (will slice the streamline, DPS and DPP). + rng_seed: int + + Returns + ------- + shuffled_sft: StatefulTractogram + The shuffled tractogram. + """ indices = np.arange(len(sft.streamlines)) random.shuffle(indices, random=rng_seed) @@ -53,6 +67,47 @@ def shuffle_streamlines(sft, rng_seed=None): return shuffled_sft +def shuffle_streamlines_orientation(sft, rng_seed=None): + """ + Shuffle the orientation of the streamlines. The streamlines are not + reversed, only the orientation of the points within the streamlines is + shuffled. + + Parameters + ---------- + sft: StatefulTractogram + The tractogram that will have its streamlines' orientation shuffled. + rng_seed: int + Random seed. + + Returns + ------- + shuffled_sft: StatefulTractogram + The shuffled tractogram. + """ + if sft.data_per_point is not None and len(sft.data_per_point) > 0: + logging.warning('Shuffling streamlines orientation. DDP will be ' + 'lost.') + + rng = np.random.RandomState(rng_seed) + streamlines = sft.streamlines + shuffled_streamlines = [] + for s in streamlines: + if len(s) < 2: + shuffled_streamlines.append(s) + else: + # flip a coin + if rng.randint(0, 2): + shuffled_streamlines.append(s[::-1]) + else: + shuffled_streamlines.append(s) + + shuffled_sft = StatefulTractogram.from_sft( + shuffled_streamlines, sft, + data_per_streamline=sft.data_per_streamline) + return shuffled_sft + + def get_axis_flip_vector(flip_axes): """ Create a flip vector from a list of axes. @@ -1166,7 +1221,7 @@ def transform_streamlines_alter(sft, min_dice=0.90, epsilon=0.01): epsilon: float Stopping criteria for convergence. The maximum difference between the dice similarity and min_dice. - + Returns ------- new_sft: StatefulTractogram diff --git a/scripts/scil_tractogram_alter.py b/scripts/scil_tractogram_alter.py index 1de8f3720..dde2d5a2b 100644 --- a/scripts/scil_tractogram_alter.py +++ b/scripts/scil_tractogram_alter.py @@ -35,7 +35,7 @@ import logging import os -from dipy.io.streamline import load_tractogram, save_tractogram +from dipy.io.streamline import save_tractogram import numpy as np from scilpy.io.streamlines import load_tractogram_with_reference @@ -48,7 +48,7 @@ remove_single_point_streamlines, cut_invalid_streamlines from scilpy.tractograms.tractogram_operations import transform_streamlines_alter, \ trim_streamlines_alter, cut_streamlines_alter, subsample_streamlines_alter, \ - replace_streamlines_alter + replace_streamlines_alter, shuffle_streamlines, shuffle_streamlines_orientation @@ -84,6 +84,9 @@ def buildArgsParser(): p.add_argument('--seed', '-s', type=int, default=None, help='Seed for RNG. Default based on --min_dice.') + p.add_argument('--shuffle', action='store_true', + help='Shuffle the streamlines and orientation after' + 'alteration.') add_overwrite_arg(p) add_verbose_arg(p) add_reference_arg(p) @@ -132,6 +135,10 @@ def main(): altered_sft = remove_single_point_streamlines(altered_sft) altered_sft = remove_overlapping_points_streamlines(altered_sft) + if args.shuffle: + altered_sft = shuffle_streamlines(altered_sft) + altered_sft = shuffle_streamlines_orientation(altered_sft) + save_tractogram(altered_sft, args.out_bundle) From 3e962dee917aec68da9ce6039130aa68d3cc0f20 Mon Sep 17 00:00:00 2001 From: frheault Date: Wed, 22 May 2024 13:37:46 -0400 Subject: [PATCH 04/12] Change permission --- scripts/scil_tractogram_alter.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) mode change 100644 => 100755 scripts/scil_tractogram_alter.py diff --git a/scripts/scil_tractogram_alter.py b/scripts/scil_tractogram_alter.py old mode 100644 new mode 100755 From 829c90b2bce356c16b1d158a7cbf537e792c2d58 Mon Sep 17 00:00:00 2001 From: frheault Date: Wed, 22 May 2024 16:26:50 -0400 Subject: [PATCH 05/12] Fix pep8 --- scilpy/tractograms/tractogram_operations.py | 9 ++++----- scripts/scil_tractogram_alter.py | 19 +++++++++---------- 2 files changed, 13 insertions(+), 15 deletions(-) diff --git a/scilpy/tractograms/tractogram_operations.py b/scilpy/tractograms/tractogram_operations.py index 8d4f2591a..32abd6f09 100644 --- a/scilpy/tractograms/tractogram_operations.py +++ b/scilpy/tractograms/tractogram_operations.py @@ -31,7 +31,8 @@ resample_streamlines_step_size, parallel_transport_streamline, \ cut_invalid_streamlines, compress_sft, cut_invalid_streamlines, \ remove_overlapping_points_streamlines, remove_single_point_streamlines -from scilpy.tractograms.streamline_and_mask_operations import cut_outside_of_mask_streamlines +from scilpy.tractograms.streamline_and_mask_operations import \ + cut_outside_of_mask_streamlines from scilpy.utils.spatial import generate_rotation_matrix MIN_NB_POINTS = 10 @@ -1107,8 +1108,6 @@ def replace_streamlines_alter(sft, min_dice=0.90, epsilon=0.01): The tractogram with replaced streamlines in the same space as the input tractogram. """ - # Import in function to avoid circular import error - from scilpy.tractanalysis.reproducibility_measures import compute_dice_voxel set_sft_logger_level(logging.ERROR) logging.debug('Upsampling the streamlines by a factor 2x to then ' @@ -1129,8 +1128,8 @@ def trim_streamlines_alter(sft, min_dice=0.90, epsilon=0.01): sft.to_vox() sft.to_corner() - original_density_map = compute_tract_counts_map(sft.streamlines, - sft.dimensions).astype(np.uint64) + original_density_map = compute_tract_counts_map( + sft.streamlines, sft.dimensions).astype(np.uint64) thr_density = [1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024] thr_pos = 0 voxels_to_remove = np.where( diff --git a/scripts/scil_tractogram_alter.py b/scripts/scil_tractogram_alter.py index dde2d5a2b..971c83192 100755 --- a/scripts/scil_tractogram_alter.py +++ b/scripts/scil_tractogram_alter.py @@ -33,7 +33,6 @@ import argparse import logging -import os from dipy.io.streamline import save_tractogram import numpy as np @@ -44,17 +43,17 @@ add_verbose_arg, assert_inputs_exist, assert_outputs_exist) -from scilpy.tractograms.streamline_operations import remove_overlapping_points_streamlines, \ +from scilpy.tractograms.streamline_operations import \ + remove_overlapping_points_streamlines, \ remove_single_point_streamlines, cut_invalid_streamlines from scilpy.tractograms.tractogram_operations import transform_streamlines_alter, \ trim_streamlines_alter, cut_streamlines_alter, subsample_streamlines_alter, \ - replace_streamlines_alter, shuffle_streamlines, shuffle_streamlines_orientation - + replace_streamlines_alter, shuffle_streamlines, shuffle_streamlines_orientation def buildArgsParser(): p = argparse.ArgumentParser(description=__doc__, - formatter_class = argparse.RawTextHelpFormatter) + formatter_class=argparse.RawTextHelpFormatter) p.add_argument('in_bundle', help='Input bundle filename (.trk, .tck).') p.add_argument('out_bundle', @@ -83,10 +82,10 @@ def buildArgsParser(): 'Only available with --cut.') p.add_argument('--seed', '-s', type=int, default=None, - help='Seed for RNG. Default based on --min_dice.') + help='Seed for RNG. Default based on --min_dice.') p.add_argument('--shuffle', action='store_true', - help='Shuffle the streamlines and orientation after' - 'alteration.') + help='Shuffle the streamlines and orientation after' + 'alteration.') add_overwrite_arg(p) add_verbose_arg(p) add_reference_arg(p) @@ -125,11 +124,11 @@ def main(): from_start=not args.from_end) elif args.replace: altered_sft = replace_streamlines_alter(sft, args.min_dice, - epsilon=args.epsilon) + epsilon=args.epsilon) elif args.transform: altered_sft = transform_streamlines_alter(sft, args.min_dice, epsilon=args.epsilon) - + # Some operations could have generated invalid streamlines altered_sft, _ = cut_invalid_streamlines(altered_sft) altered_sft = remove_single_point_streamlines(altered_sft) From 1c25833942d66fdabe10ec8fec6f2fcf8b17e03b Mon Sep 17 00:00:00 2001 From: frheault Date: Tue, 18 Jun 2024 08:25:15 -0400 Subject: [PATCH 06/12] Save matrix --- scilpy/tractograms/tractogram_operations.py | 2 +- scripts/scil_tractogram_alter.py | 15 +++++++++++++-- 2 files changed, 14 insertions(+), 3 deletions(-) diff --git a/scilpy/tractograms/tractogram_operations.py b/scilpy/tractograms/tractogram_operations.py index 32abd6f09..295b58461 100644 --- a/scilpy/tractograms/tractogram_operations.py +++ b/scilpy/tractograms/tractogram_operations.py @@ -1288,7 +1288,7 @@ def transform_streamlines_alter(sft, min_dice=0.90, epsilon=0.01): logging.debug(f'\nFinal angles: {last_pick * 360} at dice: {dice}') curr_sft.to_space(space) curr_sft.to_origin(origin) - return curr_sft + return curr_sft, rot_mat OPERATIONS = { diff --git a/scripts/scil_tractogram_alter.py b/scripts/scil_tractogram_alter.py index 971c83192..8c33e7c61 100755 --- a/scripts/scil_tractogram_alter.py +++ b/scripts/scil_tractogram_alter.py @@ -80,12 +80,16 @@ def buildArgsParser(): g.add_argument('--from_end', action='store_true', help='Cut from the other end of streamlines.\n' 'Only available with --cut.') - + g.add_argument('--save_transform', metavar='FILE', + help='Save the transformation matrix to a file.\n' + 'Only available with --transform.') + p.add_argument('--seed', '-s', type=int, default=None, help='Seed for RNG. Default based on --min_dice.') p.add_argument('--shuffle', action='store_true', help='Shuffle the streamlines and orientation after' 'alteration.') + add_overwrite_arg(p) add_verbose_arg(p) add_reference_arg(p) @@ -104,6 +108,10 @@ def main(): if args.from_end and args.cut is None: parser.error('The --from_end option is only available with --cut.') + if args.save_transform and args.transform is None: + parser.error('The --save_transform option is only available with ' + '--transform.') + if args.seed is None: np.random.seed(int(args.min_dice * 1000)) else: @@ -126,7 +134,7 @@ def main(): altered_sft = replace_streamlines_alter(sft, args.min_dice, epsilon=args.epsilon) elif args.transform: - altered_sft = transform_streamlines_alter(sft, args.min_dice, + altered_sft, matrix = transform_streamlines_alter(sft, args.min_dice, epsilon=args.epsilon) # Some operations could have generated invalid streamlines @@ -138,6 +146,9 @@ def main(): altered_sft = shuffle_streamlines(altered_sft) altered_sft = shuffle_streamlines_orientation(altered_sft) + if args.save_transform: + np.savetxt(args.save_transform, matrix) + save_tractogram(altered_sft, args.out_bundle) From a7e4e982e768cacaa10629666c384caab79011dd Mon Sep 17 00:00:00 2001 From: frheault Date: Tue, 16 Jul 2024 11:53:23 -0400 Subject: [PATCH 07/12] Added new tests, Emmanuelle comments --- scilpy/tractograms/tractogram_operations.py | 50 +++++-- scripts/scil_tractogram_alter.py | 156 -------------------- scripts/tests/test_tractogram_alter.py | 6 - 3 files changed, 38 insertions(+), 174 deletions(-) delete mode 100755 scripts/scil_tractogram_alter.py delete mode 100644 scripts/tests/test_tractogram_alter.py diff --git a/scilpy/tractograms/tractogram_operations.py b/scilpy/tractograms/tractogram_operations.py index ee1bcc33a..5e9bdb69d 100644 --- a/scilpy/tractograms/tractogram_operations.py +++ b/scilpy/tractograms/tractogram_operations.py @@ -70,9 +70,9 @@ def shuffle_streamlines(sft, rng_seed=None): def shuffle_streamlines_orientation(sft, rng_seed=None): """ - Shuffle the orientation of the streamlines. The streamlines are not - reversed, only the orientation of the points within the streamlines is - shuffled. + Shuffle the orientation of the streamlines. Iterate over streamlines + and randomly decide (50/50) if the streamline's head and tail should be + swapped. Parameters ---------- @@ -770,7 +770,11 @@ def upsample_tractogram(sft, nb, point_wise_std=None, tube_radius=None, noise_factor = polynomial(x) vec = s - new_s - vec /= np.linalg.norm(vec, axis=0) + norm = np.linalg.norm(vec, axis=0) + if np.any(norm == 0): + vec = np.ones_like(vec) + norm = np.linalg.norm(vec, axis=0) + vec /= norm new_s += vec * np.expand_dims(noise_factor, axis=1) @@ -1031,7 +1035,7 @@ def subsample_streamlines_alter(sft, min_dice=0.90, epsilon=0.01, return new_sft -def cut_streamlines_alter(sft, min_dice=0.90, epsilon=0.01, from_start=True): +def cut_streamlines_alter(sft, min_dice=0.90, epsilon=0.01, from_end=False): """ Cut streamlines based on a dice similarity metric. The function will keep removing points from the streamlines until the dice @@ -1046,9 +1050,9 @@ def cut_streamlines_alter(sft, min_dice=0.90, epsilon=0.01, from_start=True): epsilon: float Stopping criteria for convergence. The maximum difference between the dice similarity and min_dice. - from_start: bool - If True, cut from the start of the streamlines. If False, cut from the - end of the streamlines. + from_end: bool + If True, the streamlines will be cut from the end. Else, the streamlines + will be cut from the start. Returns ------- @@ -1063,7 +1067,7 @@ def cut_streamlines_alter(sft, min_dice=0.90, epsilon=0.01, from_start=True): origin = sft.origin # Uniformize endpoints to cut consistently from one end only - uniformize_bundle_sft(sft, swap=not from_start) + uniformize_bundle_sft(sft, swap=from_end) sft = resample_streamlines_step_size(sft, 0.5) sft.to_vox() sft.to_corner() @@ -1106,7 +1110,7 @@ def cut_streamlines_alter(sft, min_dice=0.90, epsilon=0.01, from_start=True): def replace_streamlines_alter(sft, min_dice=0.90, epsilon=0.01): """ Replace streamlines based on a dice similarity metric. - The function will upsamples the streamlines (with parallel transport), + The function will upsample the streamlines (with parallel transport), then downsample them until the dice similarity is close to min_dice. This effectively replaces the streamlines with new ones. @@ -1138,6 +1142,28 @@ def replace_streamlines_alter(sft, min_dice=0.90, epsilon=0.01): def trim_streamlines_alter(sft, min_dice=0.90, epsilon=0.01): + """ + Trim streamlines based on a dice similarity metric. + The function will remove low density voxels to trim streamlines until the + similarity between the original and the trimmed tractogram is close to + min_dice. + + Parameters + ---------- + sft: StatefulTractogram + The tractogram to trim. + min_dice: float + The minimum dice similarity to reach before stopping the trimming. + epsilon: float + Stopping criteria for convergence. The maximum difference between the + dice similarity and min_dice. + + Returns + ------- + new_sft: StatefulTractogram + The tractogram with trimmed streamlines in the same space as the input + tract + """ # Import in function to avoid circular import error from scilpy.tractanalysis.reproducibility_measures import compute_dice_voxel set_sft_logger_level(logging.ERROR) @@ -1287,8 +1313,8 @@ def transform_streamlines_alter(sft, min_dice=0.90, epsilon=0.01): curr_density_map = compute_tract_counts_map(curr_sft.streamlines, sft.dimensions) dice, _ = compute_dice_voxel(original_density_map, curr_density_map) - logging.debug(f'Transformed {to_pick*360} degree on axis {axis}, ' - f'dice: {dice}') + logging.debug(f'Transformed {np.round(to_pick * 360, 6)} degree on axis ' + f'{axis}, dice: {dice}') last_pick[axis] = to_pick if dice < min_dice: diff --git a/scripts/scil_tractogram_alter.py b/scripts/scil_tractogram_alter.py deleted file mode 100755 index 63d4b08ff..000000000 --- a/scripts/scil_tractogram_alter.py +++ /dev/null @@ -1,156 +0,0 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- - - -""" -This script is used to alter a bundle to reach specific minimum dice -coefficient from the original bundle. The script will subsample, trim, cut, -replace or transform the streamlines until the minimum dice is reached. -(cannot be combined in one run, use the script multiple times if needed) - -All operation use a dichotomic search to find the parameter that gets as close -as possible to the specified minimum dice coefficient (with an epsilon for -convergence). - -The subsample operation will remove streamlines until the minimum dice is -reached. This affect the whole bundle. - -The trim operation will used the lowest density voxels (starting at 1) to -remove points from the streamlines until the minimum dice is reached. -This typically affect the edge of the bundle. - -The cut operation will remove points from the start (or end) streamlines until -the minimum dice is reached. This affect one end of the bundle. - -The replace operation will upsample the tractogram (generate new streamlines -with noise) and then subsample the tractogram. This effectively replace -streamlines with similar ones until the minimum dice is reached. -This affect the whole bundle. - -The transform operation will apply random transformations to the streamlines -until the minimum dice is reached. This affect the whole bundle. -""" - -import argparse -import logging - -from dipy.io.streamline import save_tractogram -import numpy as np - -from scilpy.io.streamlines import load_tractogram_with_reference -from scilpy.io.utils import (add_overwrite_arg, - add_reference_arg, - add_verbose_arg, - assert_inputs_exist, - assert_outputs_exist) -from scilpy.tractograms.streamline_operations import \ - remove_overlapping_points_streamlines, \ - remove_single_point_streamlines, cut_invalid_streamlines -from scilpy.tractograms.tractogram_operations import transform_streamlines_alter, \ - trim_streamlines_alter, cut_streamlines_alter, subsample_streamlines_alter, \ - replace_streamlines_alter, shuffle_streamlines, shuffle_streamlines_orientation - - -def buildArgsParser(): - p = argparse.ArgumentParser(description=__doc__, - formatter_class=argparse.RawTextHelpFormatter) - p.add_argument('in_bundle', - help='Input bundle filename (.trk, .tck).') - p.add_argument('out_bundle', - help='Output bundle filename (.trk, .tck).') - - p.add_argument('--min_dice', type=float, default=0.90, - help='Minimum dice to reach.') - p.add_argument('--epsilon', type=float, default=0.01, - help='Epsilon for the convergence.') - - g = p.add_argument_group(title='Alteration options') - g1 = g.add_mutually_exclusive_group(required=True,) - g1.add_argument('--subsample', action='store_true', - help='Pick a subset of streamlines.') - g1.add_argument('--trim', action='store_true', - help='Trim streamlines using low density voxel.') - g1.add_argument('--cut', action='store_true', - help='Cut streamlines from endpoints.') - g1.add_argument('--replace', action='store_true', - help='Replace streamlines with similar ones.') - g1.add_argument('--transform', action='store_true', - help='Transform streamlines using a random linear ' - 'transformation.') - g.add_argument('--from_end', action='store_true', - help='Cut from the other end of streamlines.\n' - 'Only available with --cut.') - g.add_argument('--save_transform', metavar='FILE', - help='Save the transformation matrix to a file.\n' - 'Only available with --transform.') - - p.add_argument('--seed', '-s', type=int, default=None, - help='Seed for RNG. Default based on --min_dice.') - p.add_argument('--shuffle', action='store_true', - help='Shuffle the streamlines and orientation after' - 'alteration.') - - add_overwrite_arg(p) - add_verbose_arg(p) - add_reference_arg(p) - - return p - - -def main(): - parser = buildArgsParser() - args = parser.parse_args() - - assert_inputs_exist(parser, args.in_bundle, args.reference) - assert_outputs_exist(parser, args, [args.out_bundle]) - logging.getLogger().setLevel(logging.getLevelName(args.verbose)) - - if args.from_end and args.cut is None: - parser.error('The --from_end option is only available with --cut.') - - if args.save_transform and args.transform is None: - parser.error('The --save_transform option is only available with ' - '--transform.') - - if args.seed is None: - np.random.seed(int(args.min_dice * 1000)) - else: - np.random.seed(args.seed) - - sft = load_tractogram_with_reference(parser, args, args.in_bundle) - - # Alter the streamlines with only one operation at a time - if args.subsample: - altered_sft = subsample_streamlines_alter(sft, args.min_dice, - epsilon=args.epsilon) - elif args.trim: - altered_sft = trim_streamlines_alter(sft, args.min_dice, - epsilon=args.epsilon) - elif args.cut: - altered_sft = cut_streamlines_alter(sft, args.min_dice, - epsilon=args.epsilon, - from_start=not args.from_end) - elif args.replace: - altered_sft = replace_streamlines_alter(sft, args.min_dice, - epsilon=args.epsilon) - elif args.transform: - altered_sft, matrix = transform_streamlines_alter(sft, args.min_dice, - epsilon=args.epsilon) - - # Some operations could have generated invalid streamlines - altered_sft, _ = cut_invalid_streamlines(altered_sft) - altered_sft = remove_single_point_streamlines(altered_sft) - altered_sft = remove_overlapping_points_streamlines(altered_sft) - - if args.shuffle: - altered_sft = shuffle_streamlines(altered_sft) - altered_sft = shuffle_streamlines_orientation(altered_sft) - - if args.save_transform: - np.savetxt(args.save_transform, matrix) - - save_tractogram(altered_sft, args.out_bundle) - - -if __name__ == "__main__": - main() diff --git a/scripts/tests/test_tractogram_alter.py b/scripts/tests/test_tractogram_alter.py deleted file mode 100644 index 5c057d800..000000000 --- a/scripts/tests/test_tractogram_alter.py +++ /dev/null @@ -1,6 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- - -def test_help_option(script_runner): - ret = script_runner.run('scil_tractogram_alter.py', '--help') - assert ret.success From 54cc0cebb8e32b7fc797ed4cc191de34477d25ad Mon Sep 17 00:00:00 2001 From: frheault Date: Tue, 16 Jul 2024 11:53:38 -0400 Subject: [PATCH 08/12] Added new tests, Emmanuelle comments --- scripts/scil_bundle_alter_to_target_dice.py | 159 ++++++++++++++++++ .../tests/test_bundle_alter_to_target_dice.py | 74 ++++++++ 2 files changed, 233 insertions(+) create mode 100755 scripts/scil_bundle_alter_to_target_dice.py create mode 100644 scripts/tests/test_bundle_alter_to_target_dice.py diff --git a/scripts/scil_bundle_alter_to_target_dice.py b/scripts/scil_bundle_alter_to_target_dice.py new file mode 100755 index 000000000..30e91e1a0 --- /dev/null +++ b/scripts/scil_bundle_alter_to_target_dice.py @@ -0,0 +1,159 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + + +""" +This script is used to alter a bundle while keeping it similar enough to the +original version (keeping the dice coefficient above a given threshold). +The script will subsample, trim, cut, replace or transform the streamlines +until the minimum dice is reached. + +Various operations choices cannot be combined in one run, use the script +multiple times if needed. + +All operations use a dichotomic search to find the parameter that gets as close +as possible to the specified minimum dice coefficient (with an epsilon for +convergence). + +- The subsample operation will remove streamlines until the minimum dice is + reached. This affect the whole bundle. + +- The trim operation will used the lowest density voxels (starting at 1) to + remove points from the streamlines until the minimum dice is reached. + This typically affect the edge of the bundle. + +- The cut operation will remove points from the start (or end) streamlines until + the minimum dice is reached. This affect one end of the bundle. + +- The replace operation will upsample the tractogram (generate new streamlines + with noise) and then subsample the tractogram. This effectively replace + streamlines with similar ones until the minimum dice is reached. + This affect the whole bundle. + +- The transform operation will apply random rotations to the streamlines + until the minimum dice is reached. This affect the whole bundle. +""" + +import argparse +import logging + +from dipy.io.streamline import save_tractogram +import numpy as np + +from scilpy.io.streamlines import load_tractogram_with_reference +from scilpy.io.utils import (add_overwrite_arg, + add_reference_arg, + add_verbose_arg, + assert_inputs_exist, + assert_outputs_exist) +from scilpy.tractograms.streamline_operations import \ + remove_overlapping_points_streamlines, \ + remove_single_point_streamlines, cut_invalid_streamlines +from scilpy.tractograms.tractogram_operations import transform_streamlines_alter, \ + trim_streamlines_alter, cut_streamlines_alter, subsample_streamlines_alter, \ + replace_streamlines_alter, shuffle_streamlines, shuffle_streamlines_orientation + + +def buildArgsParser(): + p = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawTextHelpFormatter) + p.add_argument('in_bundle', + help='Input bundle filename (.trk, .tck).') + p.add_argument('out_bundle', + help='Output bundle filename (.trk, .tck).') + + p.add_argument('--min_dice', type=float, default=0.90, + help='Minimum dice to reach [%(default)s].') + p.add_argument('--epsilon', type=float, default=0.01, + help='Epsilon for the convergence [%(default)s].') + + g = p.add_argument_group(title='Alteration options') + g1 = g.add_mutually_exclusive_group(required=True,) + g1.add_argument('--subsample', action='store_true', + help='Pick a subset of streamlines.') + g1.add_argument('--trim', action='store_true', + help='Trim streamlines using low density voxel.') + g1.add_argument('--cut', action='store_true', + help='Cut streamlines from endpoints.') + g1.add_argument('--replace', action='store_true', + help='Replace streamlines with similar ones.') + g1.add_argument('--transform', action='store_true', + help='Transform streamlines using a random linear ' + 'transformation.') + g.add_argument('--from_end', action='store_true', + help='Cut from the other end of streamlines.\n' + 'Only available with --cut.') + g.add_argument('--save_transform', metavar='FILE', + help='Save the transformation matrix to a file.\n' + 'Only available with --transform.') + + p.add_argument('--seed', '-s', type=int, default=None, + help='Seed for RNG. Default based on --min_dice.') + p.add_argument('--shuffle', action='store_true', + help='Shuffle the streamlines and orientation after' + 'alteration.') + + add_overwrite_arg(p) + add_verbose_arg(p) + add_reference_arg(p) + + return p + + +def main(): + parser = buildArgsParser() + args = parser.parse_args() + + assert_inputs_exist(parser, args.in_bundle, args.reference) + assert_outputs_exist(parser, args, [args.out_bundle]) + logging.getLogger().setLevel(logging.getLevelName(args.verbose)) + + if args.from_end and not args.cut: + parser.error('The --from_end option is only available with --cut.') + + if args.save_transform and not args.transform: + parser.error('The --save_transform option is only available with ' + '--transform.') + + if args.seed is None: + np.random.seed(int(args.min_dice * 1000)) + else: + np.random.seed(args.seed) + + sft = load_tractogram_with_reference(parser, args, args.in_bundle) + + # Alter the streamlines with only one operation at a time + if args.subsample: + altered_sft = subsample_streamlines_alter(sft, args.min_dice, + epsilon=args.epsilon) + elif args.trim: + altered_sft = trim_streamlines_alter(sft, args.min_dice, + epsilon=args.epsilon) + elif args.cut: + altered_sft = cut_streamlines_alter(sft, args.min_dice, + epsilon=args.epsilon, + from_end=args.from_end) + elif args.replace: + altered_sft = replace_streamlines_alter(sft, args.min_dice, + epsilon=args.epsilon) + elif args.transform: + altered_sft, matrix = transform_streamlines_alter(sft, args.min_dice, + epsilon=args.epsilon) + + # Some operations could have generated invalid streamlines + altered_sft, _ = cut_invalid_streamlines(altered_sft) + altered_sft = remove_single_point_streamlines(altered_sft) + altered_sft = remove_overlapping_points_streamlines(altered_sft) + + if args.shuffle: + altered_sft = shuffle_streamlines(altered_sft) + altered_sft = shuffle_streamlines_orientation(altered_sft) + + if args.save_transform: + np.savetxt(args.save_transform, matrix) + + save_tractogram(altered_sft, args.out_bundle) + + +if __name__ == "__main__": + main() diff --git a/scripts/tests/test_bundle_alter_to_target_dice.py b/scripts/tests/test_bundle_alter_to_target_dice.py new file mode 100644 index 000000000..c2e476787 --- /dev/null +++ b/scripts/tests/test_bundle_alter_to_target_dice.py @@ -0,0 +1,74 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- + +import os +import tempfile + +from scilpy import SCILPY_HOME +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=['tractometry.zip']) +tmp_dir = tempfile.TemporaryDirectory() + + +def test_help_option(script_runner): + ret = script_runner.run('scil_bundle_alter_to_target_dice.py', '--help') + assert ret.success + + +def test_execution_subsample(script_runner, monkeypatch): + monkeypatch.chdir(os.path.expanduser(tmp_dir.name)) + in_bundle = os.path.join(SCILPY_HOME, 'tractometry', + 'IFGWM_uni.trk') + ret = script_runner.run('scil_bundle_alter_to_target_dice.py', + in_bundle, 'out_tractogram_subsample.trk', + '--min_dice', '0.75', '--epsilon', '0.01', + '--subsample', '--shuffle', '-v') + assert ret.success + + +def test_execution_trim(script_runner, monkeypatch): + monkeypatch.chdir(os.path.expanduser(tmp_dir.name)) + in_bundle = os.path.join(SCILPY_HOME, 'tractometry', + 'IFGWM_uni.trk') + ret = script_runner.run('scil_bundle_alter_to_target_dice.py', + in_bundle, 'out_tractogram_trim.trk', + '--min_dice', '0.75', '--epsilon', '0.01', + '--trim', '-v') + assert ret.success + + +def test_execution_cut(script_runner, monkeypatch): + monkeypatch.chdir(os.path.expanduser(tmp_dir.name)) + in_bundle = os.path.join(SCILPY_HOME, 'tractometry', + 'IFGWM_uni.trk') + ret = script_runner.run('scil_bundle_alter_to_target_dice.py', + in_bundle, 'out_tractogram_cut.trk', + '--min_dice', '0.75', '--epsilon', '0.01', + '--cut', '-v', 'DEBUG') + assert ret.success + + +def test_execution_replace(script_runner, monkeypatch): + monkeypatch.chdir(os.path.expanduser(tmp_dir.name)) + in_bundle = os.path.join(SCILPY_HOME, 'tractometry', + 'IFGWM_uni.trk') + ret = script_runner.run('scil_bundle_alter_to_target_dice.py', + in_bundle, 'out_tractogram_replace.trk', + '--min_dice', '0.75', '--epsilon', '0.01', + '--replace', '-v') + assert ret.success + + +def test_execution_transform(script_runner, monkeypatch): + monkeypatch.chdir(os.path.expanduser(tmp_dir.name)) + in_bundle = os.path.join(SCILPY_HOME, 'tractometry', + 'IFGWM_uni.trk') + ret = script_runner.run('scil_bundle_alter_to_target_dice.py', + in_bundle, 'out_tractogram_transform.trk', + '--min_dice', '0.75', '--epsilon', '0.01', + '--transform', '--save_transform', + 'transform.txt', '-v') + assert ret.success + assert os.path.isfile('transform.txt') From 5735927f710e699f1037e2f9643ced5c773f3dfe Mon Sep 17 00:00:00 2001 From: frheault Date: Thu, 25 Jul 2024 15:53:43 -0400 Subject: [PATCH 09/12] Fix script with new code and Emma comments --- scilpy/tractograms/tractogram_operations.py | 6 +++--- scripts/scil_bundle_alter_to_target_dice.py | 16 ++++++++-------- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/scilpy/tractograms/tractogram_operations.py b/scilpy/tractograms/tractogram_operations.py index 5e9bdb69d..c3d00503a 100644 --- a/scilpy/tractograms/tractogram_operations.py +++ b/scilpy/tractograms/tractogram_operations.py @@ -29,10 +29,10 @@ from scilpy.tractanalysis.streamlines_metrics import compute_tract_counts_map from scilpy.tractograms.streamline_operations import smooth_line_gaussian, \ resample_streamlines_step_size, parallel_transport_streamline, \ - cut_invalid_streamlines, compress_sft, cut_invalid_streamlines, \ + compress_sft, cut_invalid_streamlines, \ remove_overlapping_points_streamlines, remove_single_point_streamlines from scilpy.tractograms.streamline_and_mask_operations import \ - cut_outside_of_mask_streamlines + cut_streamlines_with_mask from scilpy.utils.spatial import generate_rotation_matrix MIN_NB_POINTS = 10 @@ -1225,7 +1225,7 @@ def trim_streamlines_alter(sft, min_dice=0.90, epsilon=0.01): # set logger level to ERROR to avoid logging from cut_outside_of_mask log_level = logging.getLogger().getEffectiveLevel() logging.getLogger().setLevel(logging.ERROR) - new_sft = cut_outside_of_mask_streamlines(sft, mask, min_len=10) + new_sft = cut_streamlines_with_mask(sft, mask, min_len=10) # reset logger level logging.getLogger().setLevel(log_level) diff --git a/scripts/scil_bundle_alter_to_target_dice.py b/scripts/scil_bundle_alter_to_target_dice.py index 30e91e1a0..2c5197443 100755 --- a/scripts/scil_bundle_alter_to_target_dice.py +++ b/scripts/scil_bundle_alter_to_target_dice.py @@ -16,22 +16,22 @@ convergence). - The subsample operation will remove streamlines until the minimum dice is - reached. This affect the whole bundle. + reached. This affects the whole bundle. -- The trim operation will used the lowest density voxels (starting at 1) to +- The trim operation will use the lowest density voxels (starting at 1) to remove points from the streamlines until the minimum dice is reached. This typically affect the edge of the bundle. -- The cut operation will remove points from the start (or end) streamlines until - the minimum dice is reached. This affect one end of the bundle. +- The cut operation will remove points from the start (or end) streamlines + until the minimum dice is reached. This affects one end of the bundle. - The replace operation will upsample the tractogram (generate new streamlines - with noise) and then subsample the tractogram. This effectively replace + with noise) and then subsample the tractogram. This effectively replaces streamlines with similar ones until the minimum dice is reached. - This affect the whole bundle. + This affects the whole bundle. - The transform operation will apply random rotations to the streamlines - until the minimum dice is reached. This affect the whole bundle. + until the minimum dice is reached. This affects the whole bundle. """ import argparse @@ -81,7 +81,7 @@ def buildArgsParser(): help='Transform streamlines using a random linear ' 'transformation.') g.add_argument('--from_end', action='store_true', - help='Cut from the other end of streamlines.\n' + help='Cut streamlines from the tail rather than the head.\n' 'Only available with --cut.') g.add_argument('--save_transform', metavar='FILE', help='Save the transformation matrix to a file.\n' From b57a851e48d50e3d05511a7dc8441e02c660ed87 Mon Sep 17 00:00:00 2001 From: Francois Rheault Date: Mon, 5 Aug 2024 12:46:41 -0400 Subject: [PATCH 10/12] Update scripts/scil_bundle_alter_to_target_dice.py Co-authored-by: Arnaud Bore --- scripts/scil_bundle_alter_to_target_dice.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/scil_bundle_alter_to_target_dice.py b/scripts/scil_bundle_alter_to_target_dice.py index 2c5197443..2ab8ed4ea 100755 --- a/scripts/scil_bundle_alter_to_target_dice.py +++ b/scripts/scil_bundle_alter_to_target_dice.py @@ -54,7 +54,7 @@ replace_streamlines_alter, shuffle_streamlines, shuffle_streamlines_orientation -def buildArgsParser(): +def _build_arg_parser(): p = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawTextHelpFormatter) p.add_argument('in_bundle', From 43f12c1de2b1141666967ca84a0935d0a1cbd326 Mon Sep 17 00:00:00 2001 From: Francois Rheault Date: Mon, 5 Aug 2024 12:46:45 -0400 Subject: [PATCH 11/12] Update scripts/scil_bundle_alter_to_target_dice.py Co-authored-by: Arnaud Bore --- scripts/scil_bundle_alter_to_target_dice.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/scil_bundle_alter_to_target_dice.py b/scripts/scil_bundle_alter_to_target_dice.py index 2ab8ed4ea..c8523d294 100755 --- a/scripts/scil_bundle_alter_to_target_dice.py +++ b/scripts/scil_bundle_alter_to_target_dice.py @@ -101,7 +101,7 @@ def _build_arg_parser(): def main(): - parser = buildArgsParser() + parser = _build_arg_parser() args = parser.parse_args() assert_inputs_exist(parser, args.in_bundle, args.reference) From 20d7a64e69ef3456d49b250c823e2ebfdc3122c4 Mon Sep 17 00:00:00 2001 From: frheault Date: Mon, 5 Aug 2024 13:52:09 -0400 Subject: [PATCH 12/12] Forgotten Em comments --- scilpy/tractograms/tractogram_operations.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/scilpy/tractograms/tractogram_operations.py b/scilpy/tractograms/tractogram_operations.py index c3d00503a..60b67c394 100644 --- a/scilpy/tractograms/tractogram_operations.py +++ b/scilpy/tractograms/tractogram_operations.py @@ -87,13 +87,12 @@ def shuffle_streamlines_orientation(sft, rng_seed=None): The shuffled tractogram. """ if sft.data_per_point is not None and len(sft.data_per_point) > 0: - logging.warning('Shuffling streamlines orientation. DDP will be ' + logging.warning('Shuffling streamlines orientation. DPP will be ' 'lost.') rng = np.random.RandomState(rng_seed) - streamlines = sft.streamlines shuffled_streamlines = [] - for s in streamlines: + for s in sft.streamlines: if len(s) < 2: shuffled_streamlines.append(s) else: