From 9c1394501f87e3a218f8451e9951485b1fa54dc9 Mon Sep 17 00:00:00 2001 From: arnaudbore Date: Tue, 3 Feb 2026 16:03:01 -0500 Subject: [PATCH 1/6] first test to implement DTI propagator --- src/scilpy/cli/scil_tracking_local_dev.py | 103 ++++++++++---- src/scilpy/tracking/propagator.py | 159 ++++++++++++++++++++++ 2 files changed, 234 insertions(+), 28 deletions(-) diff --git a/src/scilpy/cli/scil_tracking_local_dev.py b/src/scilpy/cli/scil_tracking_local_dev.py index 8904efad2..d25e1839d 100755 --- a/src/scilpy/cli/scil_tracking_local_dev.py +++ b/src/scilpy/cli/scil_tracking_local_dev.py @@ -64,8 +64,11 @@ assert_inputs_exist, assert_outputs_exist, parse_sh_basis_arg, verify_compression_th, load_matrix_in_any_format) +from scilpy.io.tensor import (convert_tensor_to_dipy_format, + supported_tensor_formats, + tensor_format_description) from scilpy.image.volume_space_management import DataVolume -from scilpy.tracking.propagator import ODFPropagator +from scilpy.tracking.propagator import ODFPropagator, TensorPropagator from scilpy.tracking.rap import RAPContinue from scilpy.tracking.seed import SeedGenerator, CustomSeedsDispenser from scilpy.tracking.tracker import Tracker @@ -83,14 +86,37 @@ def _build_arg_parser(): formatter_class=argparse.RawTextHelpFormatter, epilog=version_string) + # Input data options + data_g = p.add_argument_group('Input data options') + data_group = data_g.add_mutually_exclusive_group(required=True) + data_group.add_argument('--in_odf', + help='Path to the ODF SH coefficient file (.nii.gz).\n' + 'Use this for ODF-based tracking.') + data_group.add_argument('--in_tensor', + help='Path to the DTI tensor file (.nii.gz).\n' + 'Use this for tensor-based tracking.') + p.add_argument('in_seed', + help='Seeding mask (.nii.gz).') + p.add_argument('in_mask', + help='Tracking mask (.nii.gz).\n' + 'Tracking will stop outside this mask. The last point ' + 'of each \nstreamline (triggering the stopping ' + 'criteria) IS added to the streamline.') + p.add_argument('out_tractogram', + help='Tractogram output file (must be .trk or .tck).') + # Options common to both scripts - add_mandatory_options_tracking(p) - track_g = add_tracking_options(p) add_seeding_options(p) + track_g = add_tracking_options(p) # Options only for here. track_g.add_argument('--algo', default='prob', choices=['det', 'prob'], help='Algorithm to use. [%(default)s]') + track_g.add_argument('--tensor_format', type=str, default='dipy', + choices=supported_tensor_formats, + help="Format of the input tensor file.\n" + "Only used with --in_tensor. [%(default)s]\n" + + tensor_format_description) add_sphere_arg(track_g, symmetric_only=False) track_g.add_argument('--sub_sphere', type=int, default=0, @@ -178,7 +204,11 @@ def main(): parser.error('Invalid output streamline file format (must be trk or ' + 'tck): {0}'.format(args.out_tractogram)) - inputs = [args.in_odf, args.in_seed, args.in_mask] + inputs = [args.in_seed, args.in_mask] + if args.in_odf: + inputs.append(args.in_odf) + if args.in_tensor: + inputs.append(args.in_tensor) assert_inputs_exist(parser, inputs) assert_outputs_exist(parser, args, args.out_tractogram) @@ -205,7 +235,8 @@ def main(): max_nbr_pts = int(args.max_length / args.step_size) min_nbr_pts = max(int(args.min_length / args.step_size), 1) - assert_same_resolution([args.in_mask, args.in_odf, args.in_seed]) + input_data_file = args.in_odf if args.in_odf else args.in_tensor + assert_same_resolution([args.in_mask, input_data_file, args.in_seed]) # Choosing our space and origin for this tracking # If save_seeds, space and origin must be vox, center. Choosing those @@ -254,29 +285,45 @@ def main(): mask = DataVolume(mask_data, mask_res, args.mask_interp) # ------- INSTANTIATING PROPAGATOR ------- - logging.info("Loading ODF SH data.") - odf_sh_img = nib.load(args.in_odf) - odf_sh_data = odf_sh_img.get_fdata(caching='unchanged', dtype=float) - odf_sh_res = odf_sh_img.header.get_zooms()[:3] - dataset = DataVolume(odf_sh_data, odf_sh_res, args.sh_interp) - - logging.info("Instantiating propagator.") - # Converting step size to vox space - # We only support iso vox for now but allow slightly different vox 1e-3. - assert np.allclose(np.mean(odf_sh_res[:3]), - odf_sh_res, atol=1e-03) - voxel_size = odf_sh_img.header.get_zooms()[0] - vox_step_size = args.step_size / voxel_size - - # Using space and origin in the propagator: vox and center, like - # in dipy. - sh_basis, is_legacy = parse_sh_basis_arg(args) - - propagator = ODFPropagator( - dataset, vox_step_size, args.rk_order, args.algo, sh_basis, - args.sf_threshold, args.sf_threshold_init, theta, args.sphere, - sub_sphere=args.sub_sphere, - space=our_space, origin=our_origin, is_legacy=is_legacy) + if args.in_odf: + logging.info("Loading ODF SH data.") + odf_sh_img = nib.load(args.in_odf) + odf_sh_data = odf_sh_img.get_fdata(caching='unchanged', dtype=float) + odf_sh_res = odf_sh_img.header.get_zooms()[:3] + dataset = DataVolume(odf_sh_data, odf_sh_res, args.sh_interp) + + logging.info("Instantiating ODF propagator.") + voxel_size = odf_sh_img.header.get_zooms()[0] + vox_step_size = args.step_size / voxel_size + sh_basis, is_legacy = parse_sh_basis_arg(args) + + propagator = ODFPropagator( + dataset, vox_step_size, args.rk_order, args.algo, sh_basis, + args.sf_threshold, args.sf_threshold_init, theta, args.sphere, + sub_sphere=args.sub_sphere, + space=our_space, origin=our_origin, is_legacy=is_legacy) + + else: # args.in_tensor + logging.info("Loading DTI tensor data.") + tensor_img = nib.load(args.in_tensor) + tensor_data = tensor_img.get_fdata(caching='unchanged', dtype=float) + + # Convert tensor to dipy format if needed + if args.tensor_format != 'dipy': + logging.info(f"Converting tensor from {args.tensor_format} to dipy format.") + tensor_data = convert_tensor_to_dipy_format(tensor_data, args.tensor_format) + + logging.info(f"Tensor data shape: {tensor_data.shape}") + tensor_res = tensor_img.header.get_zooms()[:3] + dataset = DataVolume(tensor_data, tensor_res, args.mask_interp) + + logging.info("Instantiating tensor propagator.") + voxel_size = tensor_img.header.get_zooms()[0] + vox_step_size = args.step_size / voxel_size + + propagator = TensorPropagator( + dataset, vox_step_size, args.rk_order, args.algo, + theta, space=our_space, origin=our_origin) # ------- INSTANTIATING RAP OBJECT ------- if args.rap_mask: diff --git a/src/scilpy/tracking/propagator.py b/src/scilpy/tracking/propagator.py index 777f46c53..ececc569a 100644 --- a/src/scilpy/tracking/propagator.py +++ b/src/scilpy/tracking/propagator.py @@ -7,6 +7,7 @@ from dipy.data import get_sphere from dipy.io.stateful_tractogram import Space, Origin from dipy.reconst.shm import sh_to_sf_matrix +from dipy.reconst.dti import eig_from_lo_tri from scilpy.reconst.utils import (get_sphere_neighbours, get_sh_order_and_fullness) @@ -691,3 +692,161 @@ def _get_possible_next_dirs(self, pos, v_in): valid_volumes = np.array(valid_volumes) return valid_dirs, valid_volumes + + +class TensorPropagator(AbstractPropagator): + """ + Propagator for DTI tensor tracking. Tracks along the principal + eigenvector of the diffusion tensor. + """ + def __init__(self, datavolume, step_size, rk_order, algo, theta, + space=Space('vox'), origin=Origin('center')): + """ + Parameters + ---------- + datavolume: scilpy.image.volume_space_management.DataVolume + Trackable DataVolume object containing tensor data in lower + triangular format (6 coefficients: Dxx, Dxy, Dyy, Dxz, Dyz, Dzz). + step_size: float + The step size for tracking. + rk_order: int + Order for the Runge Kutta integration. + algo: string + Type of algorithm. Choices are 'det' or 'prob' + theta: float + Maximum angle (radians) between two steps. + space: dipy Space + Space of the streamlines during tracking. Default: VOX. + origin: dipy Origin + Origin of the streamlines during tracking. Default: center. + """ + super().__init__(datavolume, step_size, rk_order, space, origin) + + if self.space == Space.RASMM: + raise NotImplementedError( + "This version of the propagator on tensors is not ready to work " + "in RASMM space.") + + self.algo = algo + self.theta = theta + self.normalize_directions = True + self.line_rng_generator = None + + def reset_data(self, new_data=None): + return super().reset_data(new_data) + + def prepare_forward(self, seeding_pos, random_generator): + """Get initial direction from tensor at seeding position.""" + self.line_rng_generator = random_generator + + # Get tensor at seeding position + tensor_data = self.datavolume.get_value_at_coordinate( + *seeding_pos, space=self.space, origin=self.origin) + + if tensor_data is None: + logging.debug(f"Seed at {seeding_pos}: tensor_data is None") + return PropagationStatus.ERROR + + if np.all(tensor_data == 0): + logging.debug(f"Seed at {seeding_pos}: tensor_data is all zeros") + return PropagationStatus.ERROR + + # Get principal eigenvector + direction = self._get_direction_from_tensor(tensor_data) + + if direction is None: + logging.debug(f"Seed at {seeding_pos}: failed to extract direction from tensor") + return PropagationStatus.ERROR + + return TrackingDirection(direction) + + def prepare_backward(self, line, forward_dir): + """Flip direction for backward tracking.""" + # forward_dir is a TrackingDirection (which is a list) + return TrackingDirection(-np.array(forward_dir)) + + def finalize_streamline(self, last_pos, v_in): + return super().finalize_streamline(last_pos, v_in) + + def propagate(self, line, v_in): + """Propagate using Runge-Kutta integration.""" + return super().propagate(line, v_in) + + def _sample_next_direction(self, pos, v_in): + """Sample next direction from tensor.""" + tensor_data = self.datavolume.get_value_at_coordinate( + *pos, space=self.space, origin=self.origin) + + if tensor_data is None or np.all(tensor_data == 0): + return None + + # Get principal eigenvector + direction = self._get_direction_from_tensor(tensor_data) + + if direction is None: + return None + + # Check angle constraint + cosine = np.dot(v_in, direction) / (np.linalg.norm(v_in) * np.linalg.norm(direction)) + cosine = np.clip(cosine, -1, 1) + + # Flip if needed to maintain direction continuity + if cosine < 0: + direction = -direction + cosine = abs(cosine) + + if np.arccos(cosine) > self.theta: + return None + + # For probabilistic tracking, add some noise + if self.algo == 'prob' and self.line_rng_generator is not None: + # Add gaussian noise to direction + noise = self.line_rng_generator.normal(0, 0.1, 3) + direction = direction + noise + direction = direction / np.linalg.norm(direction) + + return direction + + def _get_direction_from_tensor(self, tensor_data): + """ + Extract principal eigenvector and FA from tensor data. + + Parameters + ---------- + tensor_data : ndarray + Tensor coefficients in lower triangular format (6 values). + + Returns + ------- + direction : ndarray or None + Principal eigenvector (3D direction). + """ + if len(tensor_data) != 6: + logging.warning(f"Expected 6 tensor coefficients, got {len(tensor_data)}") + return None, 0.0 + + logging.debug(f"Tensor data: {tensor_data}") + + # Compute eigenvalues and eigenvectors + # eig_from_lo_tri returns a flat array of 12 values: + # [eval1, eval2, eval3, evec1_x, evec1_y, evec1_z, evec2_x, evec2_y, evec2_z, evec3_x, evec3_y, evec3_z] + try: + result = eig_from_lo_tri(tensor_data) + evals = result[:3] # First 3 values are eigenvalues + evecs = result[3:].reshape(3, 3) # Next 9 values are eigenvectors (3x3 matrix) + logging.debug(f"Eigenvalues: {evals}, Eigenvectors shape: {evecs.shape}") + except Exception as e: + logging.debug(f"Exception in eig_from_lo_tri: {e}") + return None + + # Sort by eigenvalue magnitude (largest first) + order = np.argsort(evals)[::-1] + evals = evals[order] + evecs = evecs[:, order] + + # Principal eigenvector (associated with largest eigenvalue) + principal_evec = evecs[:, 0] + + logging.debug(f"Principal eigenvector: {principal_evec}, Sorted eigenvalues: {evals}") # Calculate FA + + return principal_evec From bcd4317883f7e2d92e649dd801ffa00d116f5730 Mon Sep 17 00:00:00 2001 From: arnaudbore Date: Wed, 4 Feb 2026 10:25:24 -0500 Subject: [PATCH 2/6] fix test --- src/scilpy/cli/tests/test_tracking_local_dev.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/scilpy/cli/tests/test_tracking_local_dev.py b/src/scilpy/cli/tests/test_tracking_local_dev.py index b7d95bdee..9e8fdda54 100644 --- a/src/scilpy/cli/tests/test_tracking_local_dev.py +++ b/src/scilpy/cli/tests/test_tracking_local_dev.py @@ -25,7 +25,7 @@ def test_execution_tracking_fodf(script_runner, monkeypatch): 'fodf.nii.gz') in_mask = os.path.join(SCILPY_HOME, 'tracking', 'seeding_mask.nii.gz') - ret = script_runner.run(['scil_tracking_local_dev', in_fodf, + ret = script_runner.run(['scil_tracking_local_dev', '--in_odf', in_fodf, in_mask, in_mask, 'local_prob.trk', '--nt', '10', '--compress', '0.1', '--sh_basis', 'descoteaux07', '--min_length', '20', '--max_length', '200', @@ -44,7 +44,7 @@ def test_execution_tracking_rap(script_runner, monkeypatch): in_rap_mask = os.path.join(SCILPY_HOME, 'tracking', 'seeding_mask.nii.gz') - ret = script_runner.run(['scil_tracking_local_dev', in_fodf, + ret = script_runner.run(['scil_tracking_local_dev', '--in_odf', in_fodf, in_mask, in_mask, 'local_prob_rap.trk', '--nt', '10', '--compress', '0.1', '--sh_basis', 'descoteaux07', @@ -68,7 +68,7 @@ def test_execution_tracking_fodf_custom_seeds(script_runner, monkeypatch): custom_seeds = [[1., 1., 1.], [2., 2., 2.], [3., 3., 3.]] np.save(in_custom_seeds, custom_seeds) - ret = script_runner.run(['scil_tracking_local_dev', in_fodf, + ret = script_runner.run(['scil_tracking_local_dev', '--in_odf', in_fodf, in_mask, in_mask, 'local_prob2.trk', '--in_custom_seeds', in_custom_seeds, '--compress', '0.1', '--sh_basis', 'descoteaux07', From 6c4061100f168b938db1bb3617d14ed267d920cb Mon Sep 17 00:00:00 2001 From: arnaudbore Date: Fri, 17 Apr 2026 12:16:18 -0400 Subject: [PATCH 3/6] add docstring + reference for log euclidean interpolation --- src/scilpy/image/volume_space_management.py | 45 +++++++++++++++++++++ 1 file changed, 45 insertions(+) diff --git a/src/scilpy/image/volume_space_management.py b/src/scilpy/image/volume_space_management.py index 85e11e67c..b895927f2 100644 --- a/src/scilpy/image/volume_space_management.py +++ b/src/scilpy/image/volume_space_management.py @@ -16,6 +16,25 @@ from dipy.reconst.shm import sf_to_sh +""" +Tensor interpolation notes +-------------------------- +The tensor interpolation path supports a Log-Euclidean mode for 6-coefficient +DTI volumes. This follows the idea of interpolating symmetric positive-definite +matrices in the log-domain and then mapping the result back with the matrix +exponential. + +References +~~~~~~~~~~ +Arsigny, V., Fillard, P., Pennec, X., Ayache, N. (2007). Geometric Means in a +Novel Vector Space Structure on Symmetric Positive-Definite Matrices. SIAM +Journal on Matrix Analysis and Applications, 29(1), 328-347. + +Pennec, X., Fillard, P., Ayache, N. (2006). A Riemannian Framework for Tensor +Computing. International Journal of Computer Vision, 66, 41-66. +""" + + class DataVolume(object): """ Class to access/interpolate data from nibabel object @@ -296,6 +315,32 @@ def _spd_expm(sym_tensor): return eigvecs @ np.diag(np.exp(eigvals)) @ eigvecs.T def _log_euclidean_interpolate4d(self, coord): + """ + Interpolate a 4D tensor volume using Log-Euclidean interpolation. + + The input volume must store diffusion tensors in lower-triangular + 6-coefficient form: [Dxx, Dxy, Dyy, Dxz, Dyz, Dzz]. Each tensor is + converted to a 3x3 SPD matrix, mapped to the log-domain, interpolated + with trilinear weights in that domain, and mapped back with the matrix + exponential. + + Parameters + ---------- + coord: ndarray shape (3,) + Coordinate in voxel space, already converted to center-origin + indexing. + + Returns + ------- + ndarray shape (6,) + Interpolated tensor coefficients in lower-triangular form. + + Notes + ----- + This implementation assumes SPD tensors. Eigenvalues are clamped to a + small positive epsilon before taking the logarithm to avoid numerical + issues on nearly-singular tensors. + """ if self.data.shape[-1] != 6: raise ValueError("log_euclidean interpolation requires 6 tensor " "coefficients per voxel.") From fef49138b2dde5beed503f3879c803748d2af63d Mon Sep 17 00:00:00 2001 From: arnaudbore Date: Fri, 17 Apr 2026 13:04:05 -0400 Subject: [PATCH 4/6] fix tests --- src/scilpy/cli/scil_tracking_local.py | 9 +++------ src/scilpy/cli/scil_tracking_local_dev.py | 6 +++++- src/scilpy/tracking/utils.py | 5 +---- 3 files changed, 9 insertions(+), 11 deletions(-) diff --git a/src/scilpy/cli/scil_tracking_local.py b/src/scilpy/cli/scil_tracking_local.py index bed9e5e7b..0f6c3bac3 100755 --- a/src/scilpy/cli/scil_tracking_local.py +++ b/src/scilpy/cli/scil_tracking_local.py @@ -77,6 +77,7 @@ verify_compression_th, load_matrix_in_any_format) from scilpy.tracking.tracker import GPUTracker from scilpy.tracking.utils import (add_mandatory_options_tracking, + add_tracking_sh_options, add_out_options, add_seeding_options, add_tracking_options, add_tracking_ptt_options, @@ -100,6 +101,7 @@ def _build_arg_parser(): # Options that are the same in this script and scil_tracking_local_dev: add_mandatory_options_tracking(p) track_g = add_tracking_options(p) + add_tracking_sh_options(p) add_seeding_options(p) # Other options, only available in this script: @@ -110,11 +112,7 @@ def _build_arg_parser(): track_g.add_argument('--algo', default='prob', choices=['det', 'prob', 'ptt', 'eudx'], help='Algorithm to use. [%(default)s]') - add_sphere_arg(track_g, symmetric_only=False) - track_g.add_argument('--sub_sphere', - type=int, default=0, - help='Subdivides each face of the sphere into 4^s new' - ' faces. [%(default)s]') + add_tracking_ptt_options(p) gpu_g = p.add_argument_group('GPU options') gpu_g.add_argument('--use_gpu', action='store_true', @@ -131,7 +129,6 @@ def _build_arg_parser(): ' [{}]'.format(DEFAULT_BATCH_SIZE)) out_g = add_out_options(p) - out_g.add_argument('--seed', type=int, help='Random number generator seed.') diff --git a/src/scilpy/cli/scil_tracking_local_dev.py b/src/scilpy/cli/scil_tracking_local_dev.py index 6bd49f444..2c626e28c 100755 --- a/src/scilpy/cli/scil_tracking_local_dev.py +++ b/src/scilpy/cli/scil_tracking_local_dev.py @@ -162,7 +162,11 @@ def _build_arg_parser(): "given the \ninitial seed. The direction is " "randomly drawn from the ODF/Tensor.") - add_tracking_sh_options(p) + sh_options = add_tracking_sh_options(p) + sh_options.add_argument('--sh_interp', default='trilinear', + choices=['nearest', 'trilinear'], + help="Spherical harmonic interpolation: " + "nearest-neighbor \nor trilinear. [%(default)s]") add_tracking_tensor_options(p) rap_g = p.add_argument_group('Region-Adaptive Propagation options') diff --git a/src/scilpy/tracking/utils.py b/src/scilpy/tracking/utils.py index beda2b322..74ee42934 100644 --- a/src/scilpy/tracking/utils.py +++ b/src/scilpy/tracking/utils.py @@ -114,15 +114,12 @@ def add_tracking_sh_options(p): "for the step function.\n" "For more information, refer to the note in the" " script description. [%(default)s]") - sh_options.add_argument('--sh_interp', default='trilinear', - choices=['nearest', 'trilinear'], - help="Spherical harmonic interpolation: " - "nearest-neighbor \nor trilinear. [%(default)s]") sh_options.add_argument('--sfthres', dest='sf_threshold', metavar='sf_th', type=float, default=0.1, help='Spherical function relative threshold. ' '[%(default)s]') add_sh_basis_args(sh_options) + return sh_options def add_tracking_ptt_options(p): track_g = p.add_argument_group('PTT options') From 782094495f3fef49823212df0e03fc836f5aaef8 Mon Sep 17 00:00:00 2001 From: arnaudbore Date: Fri, 17 Apr 2026 13:30:40 -0400 Subject: [PATCH 5/6] fix pep8 --- src/scilpy/cli/scil_tracking_local.py | 9 ++-- src/scilpy/cli/scil_tracking_local_dev.py | 26 +++++----- .../cli/tests/test_tracking_local_dev.py | 2 +- src/scilpy/tracking/tracker.py | 2 +- src/scilpy/tracking/utils.py | 50 ++++++++++--------- 5 files changed, 46 insertions(+), 43 deletions(-) diff --git a/src/scilpy/cli/scil_tracking_local.py b/src/scilpy/cli/scil_tracking_local.py index 0f6c3bac3..c1b07d1bf 100755 --- a/src/scilpy/cli/scil_tracking_local.py +++ b/src/scilpy/cli/scil_tracking_local.py @@ -41,9 +41,10 @@ * Forward tracking: For GPU tracking, the `--forward_only` flag can be used to disable backward tracking. This option isn't available for CPU tracking. - * Random number generator seed (RNG): CPU and GPU use different RNG implementations,< - so the same `--seed` is reproducible within a backend but does not guarantee - identical streamlines across CPU vs GPU tracking. + * Random number generator seed (RNG): CPU and GPU use different RNG + implementations, so the same `--seed` is reproducible within a + backend but does not guarantee identical streamlines + across CPU vs GPU tracking. All the input nifti files must be in isotropic resolution. @@ -71,7 +72,7 @@ from dipy.tracking.stopping_criterion import BinaryStoppingCriterion from dipy.tracking.tracker import eudx_tracking from scilpy.io.image import get_data_as_mask -from scilpy.io.utils import (add_sphere_arg, add_verbose_arg, +from scilpy.io.utils import (add_verbose_arg, assert_headers_compatible, assert_inputs_exist, assert_outputs_exist, parse_sh_basis_arg, verify_compression_th, load_matrix_in_any_format) diff --git a/src/scilpy/cli/scil_tracking_local_dev.py b/src/scilpy/cli/scil_tracking_local_dev.py index 2c626e28c..3fe99dc80 100755 --- a/src/scilpy/cli/scil_tracking_local_dev.py +++ b/src/scilpy/cli/scil_tracking_local_dev.py @@ -166,28 +166,28 @@ def _build_arg_parser(): sh_options.add_argument('--sh_interp', default='trilinear', choices=['nearest', 'trilinear'], help="Spherical harmonic interpolation: " - "nearest-neighbor \nor trilinear. [%(default)s]") + "nearest-neighbor \nor trilinear. [%(default)s]") add_tracking_tensor_options(p) rap_g = p.add_argument_group('Region-Adaptive Propagation options') rap_mode = rap_g.add_mutually_exclusive_group() rap_mode.add_argument('--rap_mask', default=None, help='Region-Adaptive Propagation mask (.nii.gz).\n' - 'Region-Adaptive Propagation tractography will start within ' - 'this mask.') + 'Region-Adaptive Propagation tractography will start within ' + 'this mask.') rap_mode.add_argument('--rap_labels', default=None, help='Region-Adaptive Propagation label volume (.nii.gz) .\n' - 'Voxel values are integer labels (0=background, 1..N=regions) .\n' - 'Used with --rap_method switch to select policies per label.') + 'Voxel values are integer labels (0=background, 1..N=regions) .\n' + 'Used with --rap_method switch to select policies per label.') rap_g.add_argument('--rap_method', default='None', choices=['None', 'continue', 'switch'], help="Region-Adaptive Propagation tractography method.\n" - "'continue': continues tracking with same params,\n" - "'switch': switches tracking params inside RAP mask.\n" - " [%(default)s]") + "'continue': continues tracking with same params,\n" + "'switch': switches tracking params inside RAP mask.\n" + " [%(default)s]") rap_g.add_argument('--rap_save_entry_exit', default=None, help='Save RAP entry/exit coordinates as a binary mask.\n' - 'Provide output filename (.nii.gz).') + 'Provide output filename (.nii.gz).') r_g = p.add_argument_group('Random seeding options') r_g.add_argument('--rng_seed', type=int, default=None, @@ -230,8 +230,8 @@ def main(): with open(args.rap_params, 'r') as f: rap_params = json.load(f) models = [cfg['filename'] for cfg in rap_params.get('methods', {}).values() - if 'filename' in cfg] - + if 'filename' in cfg] + assert_inputs_exist(parser, inputs + models) assert_outputs_exist(parser, args, args.out_tractogram) @@ -332,12 +332,12 @@ def main(): args.sf_threshold, args.sf_threshold_init, theta, args.sphere, sub_sphere=args.sub_sphere, space=our_space, origin=our_origin, is_legacy=is_legacy) - + elif args.in_tensor: # args.in_tensor logging.info("Loading DTI tensor data.") tensor_img = nib.load(args.in_tensor) tensor_data = tensor_img.get_fdata(caching='unchanged', dtype=float) - + # Convert tensor to dipy format if needed if args.tensor_format != 'dipy': logging.info(f"Converting tensor from {args.tensor_format} to dipy format.") diff --git a/src/scilpy/cli/tests/test_tracking_local_dev.py b/src/scilpy/cli/tests/test_tracking_local_dev.py index 553e93f52..f3ffeca45 100644 --- a/src/scilpy/cli/tests/test_tracking_local_dev.py +++ b/src/scilpy/cli/tests/test_tracking_local_dev.py @@ -79,4 +79,4 @@ def test_execution_tracking_fodf_custom_seeds(script_runner, monkeypatch): '--save_seeds', '--rng_seed', '0', '--sub_sphere', '2', '--rk_order', '4']) - assert ret.success \ No newline at end of file + assert ret.success diff --git a/src/scilpy/tracking/tracker.py b/src/scilpy/tracking/tracker.py index 183f40168..a3b4c8918 100644 --- a/src/scilpy/tracking/tracker.py +++ b/src/scilpy/tracking/tracker.py @@ -754,4 +754,4 @@ def _track(self): # output is yielded so that we can use LazyTractogram. # seed and strl with origin center (same as DIPY) - yield strl - 0.5, seed - 0.5 \ No newline at end of file + yield strl - 0.5, seed - 0.5 diff --git a/src/scilpy/tracking/utils.py b/src/scilpy/tracking/utils.py index 74ee42934..83b738eae 100644 --- a/src/scilpy/tracking/utils.py +++ b/src/scilpy/tracking/utils.py @@ -82,45 +82,47 @@ def add_tracking_options(p): def add_tracking_tensor_options(p): tensor_options = p.add_argument_group('Tensor options') tensor_options.add_argument('--tensor_format', type=str, default='dipy', - choices=supported_tensor_formats, - help="Format of the input tensor file.\n" - "Only used with --in_tensor. [%(default)s]\n" + - tensor_format_description) + choices=supported_tensor_formats, + help="Format of the input tensor file.\n" + "Only used with --in_tensor. [%(default)s]\n" + + tensor_format_description) tensor_options.add_argument('--tensor_interp', type=str, - default='log_euclidean', - choices=['nearest', 'trilinear', - 'log_euclidean'], - help="Tensor interpolation method. " - "Only used with --in_tensor. " - "[%(default)s]") + default='log_euclidean', + choices=['nearest', 'trilinear', + 'log_euclidean'], + help="Tensor interpolation method. " + "Only used with --in_tensor. " + "[%(default)s]") tensor_options.add_argument('--std', type=float, default=0.1, help="Standard deviation of the noise added " "to the direction for prob tensor-based tracking.") + def add_tracking_sh_options(p): sh_options = p.add_argument_group('Spherical harmonics options') add_sphere_arg(sh_options, symmetric_only=False) sh_options.add_argument('--sub_sphere', - type=int, default=0, - help='Subdivides each face of the sphere into 4^s new' - ' faces. [%(default)s]') + type=int, default=0, + help='Subdivides each face of the sphere into 4^s new' + ' faces. [%(default)s]') sh_options.add_argument('--sfthres_init', metavar='sf_th', type=float, - default=0.5, dest='sf_threshold_init', - help="Spherical function relative threshold value " - "for the \ninitial direction. [%(default)s]") + default=0.5, dest='sf_threshold_init', + help="Spherical function relative threshold value " + "for the \ninitial direction. [%(default)s]") sh_options.add_argument('--rk_order', metavar="K", type=int, default=1, - choices=[1, 2, 4], - help="The order of the Runge-Kutta integration used " - "for the step function.\n" - "For more information, refer to the note in the" - " script description. [%(default)s]") + choices=[1, 2, 4], + help="The order of the Runge-Kutta integration used " + "for the step function.\n" + "For more information, refer to the note in the" + " script description. [%(default)s]") sh_options.add_argument('--sfthres', dest='sf_threshold', metavar='sf_th', - type=float, default=0.1, - help='Spherical function relative threshold. ' - '[%(default)s]') + type=float, default=0.1, + help='Spherical function relative threshold. ' + '[%(default)s]') add_sh_basis_args(sh_options) return sh_options + def add_tracking_ptt_options(p): track_g = p.add_argument_group('PTT options') track_g.add_argument('--probe_length', dest='probe_length', From 2a2596611354bea003cd5f6c63d7f750c10c37c7 Mon Sep 17 00:00:00 2001 From: arnaudbore Date: Fri, 17 Apr 2026 22:38:53 -0400 Subject: [PATCH 6/6] answer Charles comments --- src/scilpy/image/volume_space_management.py | 39 ++++++----- src/scilpy/tracking/propagator.py | 72 ++++++++++++++++++--- src/scilpy/tracking/rap.py | 14 +++- 3 files changed, 94 insertions(+), 31 deletions(-) diff --git a/src/scilpy/image/volume_space_management.py b/src/scilpy/image/volume_space_management.py index b895927f2..6f815fb4f 100644 --- a/src/scilpy/image/volume_space_management.py +++ b/src/scilpy/image/volume_space_management.py @@ -16,28 +16,25 @@ from dipy.reconst.shm import sf_to_sh -""" -Tensor interpolation notes --------------------------- -The tensor interpolation path supports a Log-Euclidean mode for 6-coefficient -DTI volumes. This follows the idea of interpolating symmetric positive-definite -matrices in the log-domain and then mapping the result back with the matrix -exponential. - -References -~~~~~~~~~~ -Arsigny, V., Fillard, P., Pennec, X., Ayache, N. (2007). Geometric Means in a -Novel Vector Space Structure on Symmetric Positive-Definite Matrices. SIAM -Journal on Matrix Analysis and Applications, 29(1), 328-347. - -Pennec, X., Fillard, P., Ayache, N. (2006). A Riemannian Framework for Tensor -Computing. International Journal of Computer Vision, 66, 41-66. -""" - - class DataVolume(object): """ Class to access/interpolate data from nibabel object + + Tensor interpolation notes + -------------------------- + The tensor interpolation path supports a Log-Euclidean mode for 6-coefficient + DTI volumes. This follows the idea of interpolating symmetric positive-definite + matrices in the log-domain and then mapping the result back with the matrix + exponential. + + References + ~~~~~~~~~~ + Arsigny, V., Fillard, P., Pennec, X., Ayache, N. (2007). Geometric Means in a + Novel Vector Space Structure on Symmetric Positive-Definite Matrices. SIAM + Journal on Matrix Analysis and Applications, 29(1), 328-347. + + Pennec, X., Fillard, P., Ayache, N. (2006). A Riemannian Framework for Tensor + Computing. International Journal of Computer Vision, 66, 41-66. """ def __init__(self, data, voxres, interpolation=None, must_be_3d=False): @@ -49,8 +46,8 @@ def __init__(self, data, voxres, interpolation=None, must_be_3d=False): voxres: np.array(3,) The pixel resolution, ex, using img.header.get_zooms()[:3]. interpolation: str or None - The interpolation choice amongst "trilinear" or "nearest". If - None, functions getting a coordinate in mm instead of voxel + The interpolation choice amongst "trilinear" or "nearest" or "log_euclidean". + If None, functions getting a coordinate in mm instead of voxel coordinates are not available. must_be_3d: bool If True, dataset can't be 4D. diff --git a/src/scilpy/tracking/propagator.py b/src/scilpy/tracking/propagator.py index 337e6ef1c..738afa499 100644 --- a/src/scilpy/tracking/propagator.py +++ b/src/scilpy/tracking/propagator.py @@ -716,6 +716,10 @@ def __init__(self, datavolume, step_size, rk_order, algo, theta, Type of algorithm. Choices are 'det' or 'prob' theta: float Maximum angle (radians) between two steps. + std: float or None + Standard deviation for the Gaussian noise added to the principal + eigenvector in probabilistic tracking. If None, a default value of + 0.1 is used. Ignored for deterministic tracking. space: dipy Space Space of the streamlines during tracking. Default: VOX. origin: dipy Origin @@ -792,8 +796,9 @@ def _sample_next_direction(self, pos, v_in): if tensor_data is None or np.all(tensor_data == 0): return None - # Get principal eigenvector - direction = self._get_direction_from_tensor(tensor_data) + # Get principal eigenvector and local eigenvalues. + direction, evals = self._get_direction_from_tensor( + tensor_data, return_evals=True) if direction is None: return None @@ -812,30 +817,39 @@ def _sample_next_direction(self, pos, v_in): # For probabilistic tracking, add some noise if self.algo == 'prob' and self.line_rng_generator is not None: - # Add gaussian noise to direction - noise = self.line_rng_generator.normal(0, self.std, 3) + # Scale angular perturbation by local anisotropy (FA). + # Rationale: principal direction is less stable in low-FA regions + # (crossing/isotropic tissue) and more stable in high-FA regions, + # so we use std_local = std_global * (1 - FA). + local_std = self._compute_local_std_from_evals(evals) + noise = self.line_rng_generator.normal(0, local_std, 3) direction = direction + noise direction = direction / np.linalg.norm(direction) return direction - def _get_direction_from_tensor(self, tensor_data): + def _get_direction_from_tensor(self, tensor_data, return_evals=False): """ - Extract principal eigenvector and FA from tensor data. + Extract principal eigenvector from tensor data. Parameters ---------- tensor_data : ndarray Tensor coefficients in lower triangular format (6 values). + return_evals: bool + If True, also return sorted eigenvalues. Returns ------- direction : ndarray or None Principal eigenvector (3D direction). + evals : ndarray or None + Sorted eigenvalues (descending), returned only when + return_evals=True. """ if len(tensor_data) != 6: logging.warning(f"Expected 6 tensor coefficients, got {len(tensor_data)}") - return None, 0.0 + return (None, None) if return_evals else None logging.debug(f"Tensor data: {tensor_data}") @@ -849,7 +863,7 @@ def _get_direction_from_tensor(self, tensor_data): logging.debug(f"Eigenvalues: {evals}, Eigenvectors shape: {evecs.shape}") except Exception as e: logging.debug(f"Exception in eig_from_lo_tri: {e}") - return None + return (None, None) if return_evals else None # Sort by eigenvalue magnitude (largest first) order = np.argsort(evals)[::-1] @@ -859,6 +873,46 @@ def _get_direction_from_tensor(self, tensor_data): # Principal eigenvector (associated with largest eigenvalue) principal_evec = evecs[:, 0] - logging.debug(f"Principal eigenvector: {principal_evec}, Sorted eigenvalues: {evals}") # Calculate FA + logging.debug(f"Principal eigenvector: {principal_evec}, Sorted eigenvalues: {evals}") + if return_evals: + return principal_evec, evals return principal_evec + + def _compute_local_std_from_evals(self, evals): + """ + Compute voxel-wise noise std from local eigenvalues. + + Uses FA-derived scaling so that isotropic tensors get larger + perturbations and highly anisotropic tensors get smaller ones: + + std_local = std_global * (1 - FA) + + This is a pragmatic uncertainty proxy (not a full posterior model): + FA captures how strongly diffusion is oriented, so low FA implies less + directional confidence and therefore wider perturbation. + + References + ---------- + Basser, P. J., Mattiello, J., & LeBihan, D. (1994). MR diffusion + tensor spectroscopy and imaging. Biophysical Journal, 66(1), 259-267. + + Pierpaoli, C., & Basser, P. J. (1996). Toward a quantitative + assessment of diffusion anisotropy. Magnetic Resonance in Medicine, + 36(6), 893-906. + """ + if evals is None: + return self.std + + evals = np.asarray(evals, dtype=np.float64) + evals = np.maximum(evals, 0.0) + denom = np.sum(evals * evals) + + if denom <= 0.0: + fa = 0.0 + else: + mean_eval = np.mean(evals) + fa = np.sqrt(1.5 * np.sum((evals - mean_eval) ** 2) / denom) + fa = float(np.clip(fa, 0.0, 1.0)) + + return self.std * (1.0 - fa) diff --git a/src/scilpy/tracking/rap.py b/src/scilpy/tracking/rap.py index 9884556ea..59a8b6f2e 100644 --- a/src/scilpy/tracking/rap.py +++ b/src/scilpy/tracking/rap.py @@ -163,7 +163,19 @@ def rap_multistep_propagate(self, line, prev_direction): self.propagator = new_propagator logging.debug(f"RAP propagator switched to default label {self._propagators.keys()[0]}") - # Switch from ODF to Tensor or vice versa if needed + # Normalize previous direction representation when switching + # propagator families. + # + # ODF propagators rely on TrackingDirection.index to lookup angular + # neighborhoods on the discrete sphere (tracking_neighbours). Tensor + # propagators only need the Cartesian direction and can operate on a + # plain ndarray. + # + # Therefore: + # - Tensor -> ODF: wrap/quantize direction using prepare_backward so + # an index is available on the target ODF sphere. + # - ODF -> Tensor: drop the index and keep only Cartesian components + # to avoid carrying stale sphere metadata across models. prev_direction_has_index = getattr(prev_direction, 'index', None) is not None if hasattr(self.propagator, 'tracking_neighbours') and not prev_direction_has_index: prev_direction = self.propagator.prepare_backward(line, prev_direction)