Skip to content

NF - Pft accepts any directiongetter #12

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
gabknight wants to merge 1 commit into
base: master
Choose a base branch
from
Open

NF - Pft accepts any directiongetter #12

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 0 additions & 6 deletions dipy/tracking/local_tracking.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -394,16 +394,12 @@ def __init__(
)

self.min_wm_pve_before_stopping = min_wm_pve_before_stopping
self.directions = np.empty((maxlen + 1, 3), dtype=float)
self.pft_max_trial = pft_max_trial
self.particle_count = particle_count
self.particle_paths = np.empty(
(2, self.particle_count, pft_max_steps + 1, 3), dtype=float
)
self.particle_weights = np.empty(self.particle_count, dtype=float)
self.particle_dirs = np.empty(
(2, self.particle_count, pft_max_steps + 1, 3), dtype=float
)
self.particle_steps = np.empty((2, self.particle_count), dtype=np.intp)
self.particle_stream_statuses = np.empty(
(2, self.particle_count), dtype=np.intp
Expand Down Expand Up @@ -434,14 +430,12 @@ def _tracker(self, seed, first_step, streamline):
first_step,
self._voxel_size,
streamline,
self.directions,
self.step_size,
self.pft_max_nbr_back_steps,
self.pft_max_nbr_front_steps,
self.pft_max_trial,
self.particle_count,
self.particle_paths,
self.particle_dirs,
self.particle_weights,
self.particle_steps,
self.particle_stream_statuses,
Expand Down
98 changes: 52 additions & 46 deletions dipy/tracking/localtrack.pyx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@ from dipy.tracking.direction_getter cimport DirectionGetter
from dipy.tracking.stopping_criterion cimport(
StreamlineStatus, StoppingCriterion, AnatomicalStoppingCriterion,
TRACKPOINT, ENDPOINT, OUTSIDEIMAGE, INVALIDPOINT, PYERROR)
from dipy.utils.fast_numpy cimport cumsum, where_to_insert, copy_point
from dipy.utils.fast_numpy cimport cumsum, where_to_insert, copy_point, normalize


def local_tracker(
Expand Down Expand Up @@ -83,14 +83,12 @@ def pft_tracker(
cnp.float_t[:] first_step,
cnp.float_t[:] voxel_size,
cnp.float_t[:, :] streamline,
cnp.float_t[:, :] directions,
double step_size,
int pft_max_nbr_back_steps,
int pft_max_nbr_front_steps,
int pft_max_trials,
int particle_count,
cnp.float_t[:, :, :, :] particle_paths,
cnp.float_t[:, :, :, :] particle_dirs,
cnp.float_t[:] particle_weights,
cnp.npy_intp[:, :] particle_steps,
cnp.npy_intp[:, :] particle_stream_statuses,
Expand All @@ -116,10 +114,6 @@ def pft_tracker(
streamline : array, float, 2d, (N, 3)
Output of tracking will be put into this array. The length of this
array, ``N``, will set the maximum allowable length of the streamline.
directions : array, float, 2d, (N, 3)
Output of tracking directions will be put into this array. The length
of this array, ``N``, will set the maximum allowable length of the
streamline.
step_size : float
Size of tracking steps in mm if ``fixed_step``.
pft_max_nbr_back_steps : int
Expand All @@ -134,8 +128,6 @@ def pft_tracker(
Number of particles to use in the particle filter.
particle_paths : array, float, 4d, (2, particle_count, pft_max_steps, 3)
Temporary array for paths followed by all particles.
particle_dirs : array, float, 4d, (2, particle_count, pft_max_steps, 3)
Temporary array for directions followed by particles.
particle_weights : array, float, 1d (particle_count)
Temporary array for the weights of particles.
particle_steps : array, float, (2, particle_count)
Expand Down Expand Up @@ -170,10 +162,10 @@ def pft_tracker(
copy_point(&seed_pos[0], input_seed_pos)

i = _pft_tracker(dg, sc, input_seed_pos, input_direction, input_voxel_size,
streamline, directions, step_size, &stream_status,
streamline, step_size, &stream_status,
pft_max_nbr_back_steps, pft_max_nbr_front_steps,
pft_max_trials, particle_count, particle_paths,
particle_dirs, particle_weights, particle_steps,
particle_weights, particle_steps,
particle_stream_statuses, min_wm_pve_before_stopping)
return i, stream_status

Expand All @@ -187,77 +179,89 @@ cdef _pft_tracker(DirectionGetter dg,
double[::1] direction,
double* voxel_size,
cnp.float_t[:, :] streamline,
cnp.float_t[:, :] directions,
double step_size,
StreamlineStatus * stream_status,
int pft_max_nbr_back_steps,
int pft_max_nbr_front_steps,
int pft_max_trials,
int particle_count,
cnp.float_t[:, :, :, :] particle_paths,
cnp.float_t[:, :, :, :] particle_dirs,
cnp.float_t[:] particle_weights,
cnp.npy_intp[:, :] particle_steps,
cnp.npy_intp[:, :] particle_stream_statuses,
double min_wm_pve_before_stopping):
cdef:
cnp.npy_intp i, j
cnp.npy_intp i, j, i_sub
int pft_trial, back_steps, front_steps
int strl_array_len
double max_wm_pve, current_wm_pve
double point[3]
void (*step)(double* , double*, double) noexcept nogil

double prev_point[3]
double direction_calc[3]
double input_voxel_size[3]
cnp.float_t[:, :] sub_streamline
StreamlineStatus stream_status_1

copy_point(seed, point)
copy_point(seed, &streamline[0,0])
copy_point(&direction[0], &directions[0, 0])

copy_point(&direction[0], direction_calc)
copy_point(&voxel_size[0], input_voxel_size)
stream_status[0] = TRACKPOINT
pft_trial = 0
max_wm_pve = 0
i = 1
i_sub = 2
strl_array_len = streamline.shape[0]
while i < strl_array_len:
if dg.get_direction_c(point, direction):
# no valid diffusion direction to follow
stream_status[0] = INVALIDPOINT
else:
# If in the previous generate_streamline function at least one new point isn’t obtained, the loop breaks.
if i_sub<2:
break
# Update the point and prev_point variables and get the direction of the last step
stream_status_1 = stream_status[0]
if i>1:
Copy link
Owner Author

@gabknight gabknight Jul 24, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you re-organize the code to start the loop with the call to dg.generate_streamline. This if i>1 and the i_sub = 2 seems unnecessary if you re-order the logic.

copy_point(&streamline[i-1,0], point)
copy_point(&streamline[i-2,0], prev_point)
for j in range(3):
# step forward
point[j] += direction[j] / voxel_size[j] * step_size

copy_point(point, &streamline[i, 0])
copy_point(&direction[0], &directions[i, 0])
stream_status[0] = sc.check_point_c(point)
i += 1
direction_calc[j] = point[j] - prev_point[j]
normalize(&direction_calc[0])
sub_streamline = streamline[i-1:, :]
# Update streamline with generate_streamline
i_sub, stream_status_1 = dg.generate_streamline(point, direction_calc, input_voxel_size, step_size, sc, sub_streamline, stream_status_1, True)
stream_status[0] = stream_status_1
i += i_sub-1
copy_point(&streamline[i-1,0], point)
# update max_wm_pve
for j in range(i_sub):
Copy link
Owner Author

@gabknight gabknight Jul 24, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you add a check here, and only do the for loop if

max_wm_pve < min_wm_pve_before_stopping.

i.e. once we reached the wm, we don't need to update this variable.

pve_j = (1.0 - sc.get_include_c(&sub_streamline[j, 0]) - sc.get_exclude_c(&sub_streamline[j, 0]))
if pve_j > max_wm_pve:
max_wm_pve = pve_j

current_wm_pve = 1.0 - sc.get_include(point) - sc.get_exclude(point)
if current_wm_pve > max_wm_pve:
max_wm_pve = current_wm_pve

if stream_status[0] == TRACKPOINT:
# The tracking continues normally
continue
elif (stream_status[0] == ENDPOINT
if (stream_status[0] == ENDPOINT
and max_wm_pve < min_wm_pve_before_stopping
and current_wm_pve > 0):
# The tracking stopped before reaching the wm
continue
elif stream_status[0] == INVALIDPOINT:
elif stream_status[0] == INVALIDPOINT or stream_status[0] == TRACKPOINT:
if pft_trial < pft_max_trials and i > 1 and i < strl_array_len:
back_steps = min(i - 1, pft_max_nbr_back_steps)
front_steps = min(strl_array_len - i - back_steps - 1,
pft_max_nbr_front_steps)
front_steps = max(0, front_steps)
i = _pft(streamline, i - back_steps, directions, dg, sc,
i = _pft(streamline, i - back_steps, dg, sc,
voxel_size, step_size, stream_status,
back_steps + front_steps, particle_count,
particle_paths, particle_dirs, particle_weights,
particle_paths, particle_weights,
particle_steps, particle_stream_statuses)
pft_trial += 1
# update the current point with the PFT results
copy_point(&streamline[i-1, 0], point)
copy_point(&directions[i-1, 0], &direction[0])
# copy_point(&directions[i-1, 0], &direction[0])
Copy link
Owner Author

@gabknight gabknight Jul 24, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The variable direction is not used, no? If this is correct, could you remove the direction variable completely, and rename direction_calc to direction


# update max_wm_pve following pft
for j in range(i):
Expand Down Expand Up @@ -291,7 +295,6 @@ cdef _pft_tracker(DirectionGetter dg,
@cython.cdivision(True)
cdef _pft(cnp.float_t[:, :] streamline,
int streamline_i,
cnp.float_t[:, :] directions,
DirectionGetter dg,
AnatomicalStoppingCriterion sc,
double* voxel_size,
Expand All @@ -300,7 +303,6 @@ cdef _pft(cnp.float_t[:, :] streamline,
int pft_nbr_steps,
int particle_count,
cnp.float_t[:, :, :, :] particle_paths,
cnp.float_t[:, :, :, :] particle_dirs,
cnp.float_t[:] particle_weights,
cnp.npy_intp[:, :] particle_steps,
cnp.npy_intp[:, :] particle_stream_statuses):
Expand All @@ -311,12 +313,19 @@ cdef _pft(cnp.float_t[:, :] streamline,
double eps = 1e-16
cnp.npy_intp s, p, j

double prev_point[3]

if pft_nbr_steps <= 0:
return streamline_i

# Set the point and prev_point variables and get the direction of the last step
copy_point(&streamline[streamline_i, 0], point)
copy_point(&streamline[streamline_i-1, 0], prev_point)
for j in range(3):
dir[j] = (point[j] - prev_point[j])

for p in range(particle_count):
copy_point(&streamline[streamline_i, 0], &particle_paths[0, p, 0, 0])
copy_point(&directions[streamline_i, 0], &particle_dirs[0, p, 0, 0])
particle_weights[p] = 1. / particle_count
particle_stream_statuses[0, p] = TRACKPOINT
particle_steps[0, p] = 0
Expand All @@ -326,10 +335,13 @@ cdef _pft(cnp.float_t[:, :] streamline,
if particle_stream_statuses[0, p] != TRACKPOINT:
for j in range(3):
particle_paths[0, p, s, j] = 0
particle_dirs[0, p, s, j] = 0
continue # move to the next particle
copy_point(&particle_paths[0, p, s, 0], point)
copy_point(&particle_dirs[0, p, s, 0], dir)
if s>0:
copy_point(&particle_paths[0, p, s , 0], point)
copy_point(&particle_paths[0, p, s-1 , 0], prev_point)
for j in range(3):
dir[j] = (point[j]- prev_point[j])

if dg.get_direction_c(point, dir):
particle_stream_statuses[0, p] = INVALIDPOINT
Expand All @@ -340,7 +352,6 @@ cdef _pft(cnp.float_t[:, :] streamline,
point[j] += dir[j] / voxel_size[j] * step_size

copy_point(point, &particle_paths[0, p, s + 1, 0])
copy_point(dir, &particle_dirs[0, p, s + 1, 0])
particle_stream_statuses[0, p] = sc.check_point_c(point)
particle_steps[0, p] = s + 1
particle_weights[p] *= 1 - sc.get_exclude_c(point)
Expand Down Expand Up @@ -370,8 +381,6 @@ cdef _pft(cnp.float_t[:, :] streamline,
for ss in range(pft_nbr_steps):
copy_point(&particle_paths[0, pp, ss, 0],
&particle_paths[1, pp, ss, 0])
copy_point(&particle_dirs[0, pp, ss, 0],
&particle_dirs[1, pp, ss, 0])
particle_stream_statuses[1, pp] = \
particle_stream_statuses[0, pp]
particle_steps[1, pp] = particle_steps[0, pp]
Expand All @@ -388,8 +397,6 @@ cdef _pft(cnp.float_t[:, :] streamline,
for ss in range(pft_nbr_steps):
copy_point(&particle_paths[1, p_source, ss, 0],
&particle_paths[0, pp, ss, 0])
copy_point(&particle_dirs[1, p_source, ss, 0],
&particle_dirs[0, pp, ss, 0])
particle_stream_statuses[0, pp] = \
particle_stream_statuses[1, p_source]
particle_steps[0, pp] = particle_steps[1, p_source]
Expand All @@ -409,6 +416,5 @@ cdef _pft(cnp.float_t[:, :] streamline,
for s in range(1, particle_steps[0, p]):
copy_point(&particle_paths[0, p, s, 0],
&streamline[streamline_i + s, 0])
copy_point(&particle_dirs[0, p, s, 0], &directions[streamline_i + s, 0])
stream_status[0] = <StreamlineStatus> particle_stream_statuses[0, p]
return streamline_i + particle_steps[0, p]