Tractography

Overview

HINEC reconstructs white matter pathways by following the primary eigenvector of the diffusion tensor from seed points throughout the brain. Starting from locations in white matter, the algorithm takes small steps in the direction of maximum water diffusion, building streamlines that approximate fiber tracts. Each streamline grows bidirectionally from its seed point until a termination criterion is met.

Tracking Algorithms

Standard FACT (Euler Integration)

The simplest deterministic tracking method. At each position, the algorithm steps forward along the principal eigenvector by a fixed step size. Fast but can accumulate integration error over long tracks.

HINEC High-Order (RK4)

Fourth-order Runge-Kutta integration evaluates the eigenvector field at four intermediate points per step, producing smoother and more accurate streamlines than Euler integration.

Adaptive RKF45

Runge-Kutta-Fehlberg with adaptive step size control. The algorithm automatically adjusts the step size based on local error estimates — taking smaller steps in regions of high curvature and larger steps in straight sections. This provides the highest accuracy with optimal computational efficiency.

Seed Generation

HINEC uses a hierarchical masking strategy to determine where to place seed points:

1. Primary: White matter mask from the preprocessing pipeline
2. Fallback 1: FA-based white matter mask (FA > 0.2, morphologically eroded)
3. Fallback 2: Eroded brain mask as last resort

% Seed density: number of seeds per voxel
% density = 1: one seed at voxel center
% density = 5: five randomly distributed seeds per voxel (default)
for i = 1:size(base_seeds, 1)
    for j = 1:density
        offset = (rand(1, 3) - 0.5) * 0.8;
        seed_points(idx, :) = base_seeds(i, :) + offset;
        idx = idx + 1;
    end
end

Bidirectional Tracking

From each seed point, HINEC tracks in both the forward (+1) and backward (-1) directions along the eigenvector. The two half-tracks are concatenated to form a complete streamline. This doubles coverage compared to unidirectional tracking.

for direction = [-1, 1]
    track = track_fiber_standard_optimized(nim, seed, direction, options, cos_angle_thresh);
    if size(track, 1) > 1
        track_length_mm = sum(vecnorm(diff(track), 2, 2));
        if track_length_mm >= options.min_length
            tracks{track_count} = track;
        end
    end
end

Direction Interpolation

For sub-voxel accuracy, HINEC uses trilinear interpolation of pre-computed eigenvector components. The critical optimization is separating the 5D eigenvector array into three 3D component arrays before tracking begins, providing a ~40% speedup.

% Pre-compute eigenvector components (KEY OPTIMIZATION)
nim.v1_x = squeeze(nim.evec(:,:,:,1,1));  % X component
nim.v1_y = squeeze(nim.evec(:,:,:,2,1));  % Y component
nim.v1_z = squeeze(nim.evec(:,:,:,3,1));  % Z component
 
% During tracking: trilinear interpolation of each component
x_comp = interp3(nim.v1_x, pos(2), pos(1), pos(3), 'linear', 0);
y_comp = interp3(nim.v1_y, pos(2), pos(1), pos(3), 'linear', 0);
z_comp = interp3(nim.v1_z, pos(2), pos(1), pos(3), 'linear', 0);

Termination Criteria

Tracking terminates when any of the following conditions is met:

Parameters

Output Structure

Results are saved to tractography_results/ with timestamped filenames:

tractography_results/
└── tracks_standard_2025-01-15_14_30_22.mat
    ├── tracks          % Cell array of Nx3 fiber coordinates
    ├── options         % Parameters used
    ├── elapsed_time    % Computation time (seconds)
    └── algorithm       % Algorithm name string