[ENH][WIP] Add DTI tracking

[arnaudbore]

Quick description

DTI Tensor Tracking in scil_tracking_local_dev

Overview

The script now supports tracking from diffusion tensor imaging (DTI) data in addition to ODF-based tracking. This allows streamline tractography using the principal eigenvector of the diffusion tensor.

Usage

Tensor-based tracking

scil_tracking_local_dev.py \
    --in_tensor tensor.nii.gz \
    input_seed.nii.gz \
    input_mask.nii.gz \
    output_tractogram.trk \
    --fa_threshold 0.1 \
    --step_size 0.5 \
    --theta 45 \
    [other tracking options...]

ODF-based tracking (original functionality)

scil_tracking_local_dev.py \
    --in_odf odf.nii.gz \
    input_seed.nii.gz \
    input_mask.nii.gz \
    output_tractogram.trk \
    --sf_threshold 0.1 \
    --step_size 0.5 \
    --theta 45 \
    [other tracking options...]

Key Parameters for Tensor Tracking

--in_tensor

Path to the DTI tensor file (.nii.gz). The tensor must be in lower triangular format with 6 coefficients in the 4th dimension:

• Dxx, Dxy, Dyy, Dxz, Dyz, Dzz

This is the standard format output by most DTI fitting tools including scil_dti_metrics.py.

--algo

Algorithm for direction selection:

• det: Deterministic tracking (follows principal eigenvector exactly)
• prob: Probabilistic tracking (adds gaussian noise to eigenvector)
• Default: prob

How It Works

1. Direction Extraction: At each tracking step, the tensor is interpolated at the current position and decomposed into eigenvalues and eigenvectors.

2. Principal Direction: The eigenvector corresponding to the largest eigenvalue is used as the tracking direction.

3. Angle Constraint: The angle between consecutive steps must not exceed theta (in degrees).

4. Probabilistic Mode: When --algo prob is used, Gaussian noise is added to the direction to introduce variability.

Tensor Format

The input tensor file must have shape (X, Y, Z, 6) with the 6 tensor components in lower triangular order:

[Dxx, Dxy, Dyy, Dxz, Dyz, Dzz]

This corresponds to the full symmetric 3x3 tensor:

| Dxx  Dxy  Dxz |
| Dxy  Dyy  Dyz |
| Dxz  Dyz  Dzz |

Example Workflow

1. Fit DTI model:

scil_dti_metrics.py \
    dwi.nii.gz \
    bvals bvecs \
    --mask brain_mask.nii.gz \
    --tensor tensor.nii.gz \
    --fa fa.nii.gz

2. Track using tensor:

scil_tracking_local_dev.py \
    --in_tensor tensor.nii.gz \
    wm_seed.nii.gz \
    wm_mask.nii.gz \
    tractogram.trk \
    --fa_threshold 0.15 \
    --algo det \
    --step_size 0.5 \
    --theta 45 \
    --nt 1000000

Comparison with ODF Tracking

Feature	Tensor Tracking	ODF Tracking
Input	DTI tensor (6 coefficients)	SH coefficients
Direction	Principal eigenvector	Peak on sphere
Stopping criterion	... (I could put FA)	SF threshold
Best for	Single fiber populations	Crossing fibers
Speed	Faster	Slower

Notes

• Tensor tracking works best in regions with single fiber populations
• For crossing fibers, ODF-based tracking is recommended
• All other tracking parameters (seeding, step size, angle, etc.) work the same for both modes
• RAP (Region-Adaptive Propagation) is compatible with tensor tracking

...

Type of change

Check the relevant options.

• Bug fix (non-breaking change which fixes an issue)
• New feature (non-breaking change which adds functionality)
• Breaking change (fix or feature that would cause existing functionality to not work as expected)
• This change requires a documentation update

Provide data, screenshots, command line to test (if relevant)

...

Checklist

• My code follows the style guidelines of this project (run autopep8)
• I added relevant citations to scripts, modules and functions docstrings and descriptions
• I have performed a self-review of my code
• I have commented my code, particularly in hard-to-understand areas
• I have made corresponding changes to the documentation
• My changes generate no new warnings
• I moved all functions from the script file (except the argparser and main) to scilpy modules
• I have added tests that prove my fix is effective or that my feature works
• New and existing unit tests pass locally with my changes

[codecov]

Codecov Report

❌ Patch coverage is 31.25000% with 77 lines in your changes missing coverage. Please review.
✅ Project coverage is 72.71%. Comparing base (b604b99) to head (34d4fd4).

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1304      +/-   ##
==========================================
- Coverage   72.91%   72.71%   -0.21%     
==========================================
  Files         295      295              
  Lines       25453    25549      +96     
  Branches     3567     3581      +14     
==========================================
+ Hits        18560    18577      +17     
- Misses       5389     5463      +74     
- Partials     1504     1509       +5     

Flag	Coverage Δ
smoketests	69.55% <31.25%> (-0.20%)	⬇️
unittests	14.22% <0.00%> (-0.06%)	⬇️

Flags with carried forward coverage won't be shown. Click here to find out more.

Components	Coverage Δ
Scripts	75.34% <59.52%> (-0.10%)	⬇️
Library	69.57% <14.28%> (-0.34%)	⬇️

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[frheault]

How does the tensor interpolation work? Linear interpolation of the 3x3 matrix?
This does not preserve positive-definiteness and *can) cause swelling artifacts (this would change FA/MD, but I do not know the effect on direction), maybe @mdesco would know.

It is most likely not worth it, but I believe to interpolate a tensor you need Log-Euclidean or Riemannian interpolation.

[gabknight]

Deterministic DTI tracking seem to be working fine. Probabilistic tracking is given the same output as det, not sure why. @arnaudbore do you get the expected tractograms on your side?

scil_tracking_local_dev --in_tensor dti.nii.gz DiSCo1_ROIs-binary.nii.gz DiSCo1_mask.nii.gz dti_tracking_prob.tck --step 0.5 --theta 45 --tensor_format dipy --algo prob

Typo in the PR description:

• --fa_threshold 0.1 \ ---> option not available

• --step_size 0.5 \ --> --step 0.5

[gabknight]

Can you make the normal std a member of the TensorPropagator class, with these default values? Potentially adding it too to scil_tracking_local_dev.py if these is use cases.

[gabknight]

May be less intuitive, but instead of 'det' or 'prob' could be just set the noise STD, in the doc you could suggest 0.1 std provide good prob results. if STD is 0 (default), you skip this step.