NF - Particle Filtering Tractography (merge) #1384
Conversation
dipy/tracking/local/localtrack.pyx
Outdated
| streamline[streamline_i + s, j] = particle_paths[0, p, s, j] | ||
| directions[streamline_i + s, j] = particle_dirs[0, p, s, j] | ||
|
|
||
| ###TODO fix this |
There was a problem hiding this comment.
@nilgoyette, do you know why cython doesn't let me do tissue_class[0] = TissueClass(particle_states[0, p, 0])
NameError: name 'TissueClass' is not defined
I do have
from .tissue_classifier cimport TissueClass
There was a problem hiding this comment.
It seems that a cython enum doesn't "exists". It's simply defining the constants in the enum. You can cast it to the right type though:
tissue_class[0] = <TissueClass>particle_states[0, p, 0]
I tested and it compiles and all the tests pass. There's no way to check if the particle_states returns a valid TissueClass but hey, it's cython!
Codecov Report
@@ Coverage Diff @@
## master #1384 +/- ##
==========================================
+ Coverage 87.07% 87.19% +0.11%
==========================================
Files 227 230 +3
Lines 29063 29497 +434
Branches 3125 3171 +46
==========================================
+ Hits 25307 25719 +412
- Misses 3046 3061 +15
- Partials 710 717 +7
Continue to review full report at Codecov.
|
|
Thanks for the quick feedback @nilgoyette. That's was the line I needed :) @StongeEtienne it would be awesome if you could test this PFT implementation with SET. You will notice the speedup! |
|
No @gabknight and It is weird that it appears only on one build. I will restart this build to be sure that it is reproducible before investigating it |
|
Thank again @nilgoyette for noticing this bug. I much happier to find it now :p Here is the whole brain tracking (stanford_hardi) with and without PFT, same seeding, same tracking parameters.
|
Codecov Report
@@ Coverage Diff @@
## master #1384 +/- ##
==========================================
- Coverage 87.07% 86.75% -0.33%
==========================================
Files 227 237 +10
Lines 29063 30008 +945
Branches 3125 3232 +107
==========================================
+ Hits 25307 26033 +726
- Misses 3046 3239 +193
- Partials 710 736 +26
Continue to review full report at Codecov.
|
dipy/tracking/local/localtracking.py
Outdated
| / step_size)) | ||
| self.pft_max_nbr_front_steps = int(np.ceil(pft_front_tracking_dist | ||
| / step_size)) | ||
| self.pft_max_steps = (self.pft_max_nbr_back_steps |
There was a problem hiding this comment.
pft_max_steps is never used outsite of this constructor. Should not be in self
Or use it. You have a back_steps + front_steps in your call to _pft.
There was a problem hiding this comment.
back_steps + front_steps is variable (depending on the current tracking point). I removed it from self. thx.
dipy/tracking/local/localtrack.pyx
Outdated
| double point[3] | ||
| double dir[3] | ||
| double voxdir[3] | ||
| double eps = np.finfo(np.float).eps |
There was a problem hiding this comment.
Just to make sure, do you want float32.eps or float64.eps?
There was a problem hiding this comment.
float64 should be fine as we always use double in trilinear_interpolate4d. It might be safer to put eps slightly higher as we do multiple additions in trilinear_interpolate4d.
dipy/tracking/local/localtrack.pyx
Outdated
|
|
||
| sum_weights = 0 | ||
| for p in range(particle_count): | ||
| sum_weights += particle_weights[p] |
There was a problem hiding this comment.
Just to be sure, are you aware that sometimes your weights all equal NaN?
There was a problem hiding this comment.
I wasn't. This happens when all the weighs of the particle are equal to 0. In this case, all particle paths are equally bad and we return the path of the first particle, with the TissueClass INVALIDPOINT. i.e. PFT failed to find an alternative valid path. Having weights = 0 or nan produces the same output. I now removed the division that made the nans. thx.
dipy/tracking/local/localtrack.pyx
Outdated
| sum_squared += particle_weights[p] * particle_weights[p] | ||
|
|
||
| # Resample the particles if the weights are too uneven. | ||
| # Particles with negligable weights are replaced by duplicates of |
dipy/tracking/local/localtrack.pyx
Outdated
| # Particles with negligable weights are replaced by duplicates of | ||
| # those with high weigths through resampling | ||
| N_effective = 1. / sum_squared | ||
| if N_effective < particle_count / 10.: |
There was a problem hiding this comment.
Using your test (without the fixed random), I couldn't get this to happen. In fact, everytime I tested, all elements in the arrays (paths, dirs, states, weights) are perfectly equal in the p dimension. Looking at the code, I don't see how they could be different. If they are actually always equal, then this algo is much too complicated!
I probably don't understand well enough. Can you please explain how it can happen? Is it only because I'm using test data?
There was a problem hiding this comment.
On my side, the resampling happens and the weights, dirs, states, path are different.
try running nosetests --nocapture -v in dipy/tracking/local/
and adding a print at line localtrack.pyx:463 print("Resampling")
add adding a print at line localtracking.py:239 print(self.particle_particle_weights)
You can also try running the example: python doc/example/particle_filtering_tractography.py to check if pft_streamlines and 'local_prob_streamlines are different.
In _pft(.), we do probabilistic tractography for P particles and N steps. at each tracking step, the P particle weights are updates following the information provided by the TissueClassifier (localtrack.pyx:line 441). After the N steps, we sample 1 particle from the P particles using their weight (higher weight, more likely to be sampled). The path of this particle is the one used in the final streamline.
The thing here is that P particles might be very few particles to properly explore a given configuration. To overcome that, we focus the exploration in regions more promising (i.e. particles with higher weights). To do that, we discard particles with low weight and duplicate particles with high weights. This is done when the weight distribution is unbalanced (e.g. when 1. / sum_squared < particle_count / 10).
dipy/tracking/local/localtrack.pyx
Outdated
| fixed_step(point, voxdir, step_size) | ||
| copypoint(point, &particle_paths[0, p, s + 1, 0]) | ||
| copypoint(dir, &particle_dirs[0, p, s + 1, 0]) | ||
| particle_states[0, p, 0] = tc.check_point_c(point) |
There was a problem hiding this comment.
Any reason why you shared the particle tissue_class with the particle length? IMO, it's harder to understand for no real gain. Shouldn't they be 2 arrays?
There was a problem hiding this comment.
I thought that would be more clear to group them to reduce the number of arrays... Let me try with 2 arrays.
| if (self.pft_max_nbr_front_steps < 0 | ||
| or self.pft_max_nbr_back_steps < 0 | ||
| or pft_max_steps < 1): | ||
| raise ValueError("The number of PFT steps must be greater than 0.") |
There was a problem hiding this comment.
Wouldn't <= 0.0 be better for pft_max_nbr_front_steps and pft_max_nbr_back_steps? It seems to me that pft_back_tracking_dist and pft_front_tracking_dist should be > 0.0. No?
There was a problem hiding this comment.
At least one of the 2 must be >0, otherwise, there is no PFT. Putting pft_max_nbr_front_steps = 0 is conceivable, with a pft_back_tracking_dist > 0, or the other way around.
With CMC (probabilistic stopping criteria), pft_back_tracking_dist can be 0 without issue, that wouldn't work for ACT (deterministic stopping criteria), i.e. with ACT and pft_back_tracking_dist = 0, the output would be the same as without PFT.
|
Any other review or request? If not, I will merge it this Friday. |
|
Retested integration in our normal work, works well. It is a +1 from me. Good job @gabknight and everyone! |
|
Give me time to review this until this Sunday. Sorry busy week. |
| be used in conjunction with PFT. In this example, we used CMC. | ||
| """ | ||
|
|
||
| import matplotlib.pyplot as plt |
| step_size=step_size, | ||
| maxlen=1000, | ||
| return_all=False) | ||
| probabilistic_streamlines = [s for s in probabilistic_streamlines] |
There was a problem hiding this comment.
This is the old standard for generating streamlines which has memory issues what you should use is
from dipy.tracking.streamline import Streamlines
streamlines = Streamlines(probabilistic_streamlines)
See example tracking_quick_start.py for more.
| maxlen=1000, | ||
| return_all=False) | ||
| probabilistic_streamlines = [s for s in probabilistic_streamlines] | ||
| save_trk("probabilistic_streamlines.trk", probabilistic_streamlines, affine, |
| pft_front_tracking_dist=1, | ||
| particle_count=15, | ||
| return_all=False) | ||
| streamlines = [s for s in pft_streamlines] |
|
Apart from a few comments. I think this is ready. |
|
Thank you for the feedback @Garyfallidis, I made the changes. Thank you all for your help and feedback on this PR! |
|
Merging 🌟 ! Many years on the making for this PR. Good luck! And @gabknight lets make a workflow :) |
|
🍺 thanks @gabknight for the hard work and everyone (@nilgoyette ) for the reviews! |
NF - Particle Filtering Tractography (merge)
This PR replaces #1340, and follows the changes of PR #1342.
This PR add a new tissue classifier (CMC) and a new tractography algorithm (PFT). Both methods are described in [Girard et al., NImg, 2014].
CMC - new tissue Classifier. CmcTissueClassifier used probabilities and the tissue partial volume fraction maps to determine when the tracking stops. As it is very similar to ActTissueClassifier, both classes now extend an a new abstract class ConstrainedTissueClassifier.
PFT - new tractography method. ParticleFilteringTracking is inherited from the LocalTracking class. The particle filter is used to find a valid streamline path when the streamline being traced reach the state INVALIDPOINT.
An example is giving in doc/examples/particle_filtering_fiber_tracking.py
References
Girard, G., Whittingstall, K., Deriche, R., & Descoteaux, M. Towards quantitative connectivity analysis: reducing tractography biases. NeuroImage, 98, 266-278, 2014.