Skip to content
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

WIP - PTT PFT #11

Open
wants to merge 6 commits into
base: master
Choose a base branch
from
+94 −43
Open

WIP - PTT PFT #11

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
2a6166b
Eliminate the use of variable directions
JesusMda Oct 16, 2024
0140cee
Using the generate_streamline function in pft_trakcer
JesusMda Nov 1, 2024
cdd1ce7
An ENDPOINT was added in the generate_streamline of PTT.
JesusMda Nov 8, 2024
61b65f7
The parameters directions and particles_dirs are removed.
JesusMda Nov 9, 2024
0294872
BF - streamstatus
gabknight Nov 22, 2024
d308ad2
Merge pull request #1 from gabknight/BF_streamstatus
JesusMda Nov 26, 2024
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
12 changes: 6 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,16 @@ def __init__(
)

self.min_wm_pve_before_stopping = min_wm_pve_before_stopping
self.directions = np.empty((maxlen + 1, 3), dtype=float)
# 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_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 +434,14 @@ def _tracker(self, seed, first_step, streamline):
first_step,
self._voxel_size,
streamline,
self.directions,
# 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_dirs,
self.particle_weights,
self.particle_steps,
self.particle_stream_statuses,
Expand Down
125 changes: 88 additions & 37 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,14 @@ def pft_tracker(
cnp.float_t[:] first_step,
cnp.float_t[:] voxel_size,
cnp.float_t[:, :] streamline,
cnp.float_t[:, :] directions,
# 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_dirs,
cnp.float_t[:] particle_weights,
cnp.npy_intp[:, :] particle_steps,
cnp.npy_intp[:, :] particle_stream_statuses,
Expand Down Expand Up @@ -170,10 +170,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 +187,114 @@ cdef _pft_tracker(DirectionGetter dg,
double[::1] direction,
double* voxel_size,
cnp.float_t[:, :] streamline,
cnp.float_t[:, :] directions,
# 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_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], &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:
# print(stream_status[0])
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_point(&streamline[i-1,0], point)
copy_point(&streamline[i-2,0], prev_point)
for j in range(3):
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
# print(stream_status[0])
# if stream_status[0] == TRACKPOINT:
# The tracking continues normally
# print("Se sale con status TRACKPOINT", i_sub)
#if dg.get_direction_c(point, direction):
# no valid diffusion direction to follow
#stream_status[0] = INVALIDPOINT
#else:
#for j in range(3):
# step forward
point[j] += direction[j] / voxel_size[j] * step_size
#point[j] += direction[j] / voxel_size[j] * step_size

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

i += i_sub-1
copy_point(&streamline[i-1,0], point)
# stream_status[0] = sc.check_point_c(point)

copy_point(point, &streamline[i, 0])
copy_point(&direction[0], &directions[i, 0])
stream_status[0] = sc.check_point_c(point)
i += 1
# update max_wm_pve
for j in range(i_sub):
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)
Copy link
Owner Author

Choose a reason for hiding this comment

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

this block needs to check all generated points, no only le last point.

      # update max_wm_pve
        for j in range(i_sub):
            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

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])

# update max_wm_pve following pft
for j in range(i):
Expand Down Expand Up @@ -291,7 +328,7 @@ cdef _pft_tracker(DirectionGetter dg,
@cython.cdivision(True)
cdef _pft(cnp.float_t[:, :] streamline,
int streamline_i,
cnp.float_t[:, :] directions,
# cnp.float_t[:, :] directions,
DirectionGetter dg,
AnatomicalStoppingCriterion sc,
double* voxel_size,
Expand All @@ -300,7 +337,7 @@ 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_dirs,
cnp.float_t[:] particle_weights,
cnp.npy_intp[:, :] particle_steps,
cnp.npy_intp[:, :] particle_stream_statuses):
Expand All @@ -311,12 +348,20 @@ 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])
# 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 +371,16 @@ 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
# 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)
# copy_point(&particle_dirs[0, p, s, 0], dir)
# Update the last direction
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 +391,7 @@ 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])
# 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 +421,8 @@ 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])
# 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 +439,8 @@ 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])
# 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 +460,6 @@ 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])
# 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]