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.

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

Traco Redesign #341

Merged
merged 68 commits into from Dec 8, 2014

Conversation

Projects
None yet
8 participants
@MrBago
Contributor

MrBago commented Mar 25, 2014

I've put together a pull request to help clarify some of the ideas I proposed earlier on this thread.

In this pull request there is an example called newApi.py this example still needs some work, but I hope it will help show how this new tracking framework is meant to be used. In the example I've done some probabilistic tracking with the CSA and CSD models and some deterministic tracking with the CSA model. Be warned that the CSD fitting takes about 45 min.

To describe how this works, in each of the three tracking examples I create a direction getter and a tissue classifier which I pass with seeds and affine and a step_size to LocalTracking.

LocalTracking is an example of what we've been calling a “streamline manager” it tracks each seed and yields a streamline if the streamline doesn't hit any invalid points during tracking.

For the deterministic tracking I've modified the PeaksAndMetrics class to make it a bonefied DirectionGetter. It calls existing cython code that was written for EuDx.

For something to be a DirectionGetter it needs to subclass from dipy.tracking.local.DirectionGetter and implement a cython cdef method called get_direction. The cython requirement is so the inner tracking loop can be compiled. A DirectionGetter needs to also have an initial_direction method. This method is python for now (but we can change this).

For the probabilist tracking the direction getter is ProbabilisticOdfWightedDirectionGetter (I'm open to a more succinct name). I've written the it in python, but used a cython helper type to make it work. Let me know if anyone wants to hear me rant about how this awesomeness works.

I've implemented one example tissue classifier, ThresholdTissueClassifier, which like direction getter has a cython method called check_point. This method returns the type of tissue to the streamline manager/framework which then has to decide what to do next.

This PR needs more tests, a better example and there are some corner cases to deal with (literally the interpolation breaks at the corners of the image). But I promised you guys some code so here it is

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Sep 14, 2014

@MrBago the newApi.py is the only example I see and it does not work. Can you fix it and rebase also this PR? First I got this error

ImportError: cannot import name ProbabilisticOdfWightedDirectionGetter
and then I got this one

AttributeError: 'SphHarmFit' object has no attribute 'get_pmf'

Thx! GGT!

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Sep 14, 2014

Or maybe this should close in favour of @gabknights PR #414 ? That worked for me.

@MrBago

This comment has been minimized.

Contributor

MrBago commented Sep 14, 2014

The two need to be merged, I'll fix the issue tonight when I get home.

Bago
On Sep 13, 2014 7:28 PM, "Eleftherios Garyfallidis" <
notifications@github.com> wrote:

Or maybe this should close in favour of @gabknights PR? That worked for me.


Reply to this email directly or view it on GitHub
#341 (comment).

@MrBago

This comment has been minimized.

Contributor

MrBago commented Sep 14, 2014

I've fixed the example, should work now. Keep in mind that the spherical
deconvolution step takes a long time, most of the run time is in the
model.fit step. I'm open to any ideas that make this faster :).

Bago

On Sat, Sep 13, 2014 at 7:31 PM, Bago mrbago@gmail.com wrote:

The two need to be merged, I'll fix the issue tonight when I get home.

Bago
On Sep 13, 2014 7:28 PM, "Eleftherios Garyfallidis" <
notifications@github.com> wrote:

Or maybe this should close in favour of @gabknights PR? That worked for
me.


Reply to this email directly or view it on GitHub
#341 (comment).

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Sep 19, 2014

@MrBago and @gabknight I got memory error on a machine with 16GBytes just by running the newApi.py example.
Rerunning now to see if that will happen again after recompiling the cython modules.

In the meantime. Please rebase, correct typos, add documentation, add more tests and check coverage. Oh and merge with #414 if you haven't done already.

@gabknight please prepare the tractometer tests for validation so that it is easy to rerun them as the API goes into place.

I will try to think of a way around for the fit part to be faster.

Also I see that you created new interpolation code. It was impossible to use map_coordinates from scipy.ndimage?

Guys, keep the momentum, work hard. Let's have an incredible tracking module in Dipy. This is something we really miss at the moment!

Keep it up! You are the noble tracking team! Such a challenge!

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Sep 19, 2014

Sorry, but got this memory error again. Here is the stdout:

MemoryError Traceback (most recent call last)
/usr/lib/python2.7/dist-packages/IPython/utils/py3compat.pyc in execfile(fname, *where)
202 else:
203 filename = fname
--> 204 builtin.execfile(filename, *where)

/home/eleftherios/Devel/dipy/doc/examples/newApi.py in ()
127 lambda_=np.sqrt(1./2))
128 start = time.time()
--> 129 prob_tracking_example(csdmodel, data, mask, N, hdr, "SphereDeconv.trk")
130 print time.time() - start
131

/home/eleftherios/Devel/dipy/doc/examples/newApi.py in prob_tracking_example(model, data, mask, N, hdr, filename)
61
62 # Create objects to be passed to tracker
---> 63 pdg = ProbabilisticDirectionGetter.fromShmFit(fit, default_sphere, 45.)
64 gfa = fit.gfa
65 gfa = np.where(np.isnan(gfa), 0., gfa)

/home/eleftherios/Devel/dipy/dipy/tracking/local/direction_getter_py.pyc in fromShmFit(klass, shmFit, sphere, max_angle)
94
95 """
---> 96 pmf_gen = ShmFitPmfGen(shmFit, sphere)
97 return klass(pmf_gen, sphere, max_angle)
98

/home/eleftherios/Devel/dipy/dipy/tracking/local/direction_getter_py.pyc in init(self, shmfit, sphere)
35 self.fit = shmfit
36 self.sphere = sphere
---> 37 self._B = shmfit.sampling_matrix(sphere)
38 self._coeff = shmfit.shm_coeff
39

/home/eleftherios/Devel/dipy/dipy/reconst/multi_voxel.pyc in call(self, _args, *_kwargs)
71 if item is not None:
72 result[ijk] = item(_args, *_kwargs)
---> 73 return _squash(result)

/home/eleftherios/Devel/dipy/dipy/reconst/quick_squash.so in dipy.reconst.quick_squash.quick_squash (dipy/reconst/quick_squash.c:2859)()

I have cython version 0.21

@MrBago MrBago force-pushed the MrBago:tracking_WIP branch from 4f89efd to 2eb058c Sep 19, 2014

@MrBago

This comment has been minimized.

Contributor

MrBago commented Sep 19, 2014

Sorry @Garyfallidis it should be working now, try it again. Also the make the fit faster, I'd like to see the peaks from model parallel machinery adapted so it could work for multi voxel fit too, that could speed things up.

def makeNd(array, N):
"""Makes an array that's less than then Nd - Nd

This comment has been minimized.

@stefanv

stefanv Sep 22, 2014

Contributor

typo

"""Makes an array that's less than then Nd - Nd
We need this because numpy 1.6 does not return a "c contiguous array"
when you call ``array(a, order='c', ndmin=N)``

This comment has been minimized.

@stefanv

stefanv Sep 22, 2014

Contributor

Was this a bug? Because that's a requirement of nd.array!

cdef class PAMDirectionGetter(DirectionGetter):

This comment has been minimized.

@stefanv

stefanv Sep 22, 2014

Contributor

Documentation missing

self.ang_thr = 60
self.total_weight = .5
def _initialize(self):

This comment has been minimized.

@stefanv

stefanv Sep 22, 2014

Contributor

Why not use __cinit__?

This comment has been minimized.

@Garyfallidis

Garyfallidis Sep 25, 2014

Member

_initialize is called from outside so it shouldn't be an "internal" function. I suggest to remove the underscore and rethink of the name or put these in cinit as Stefan suggests.

This comment has been minimized.

@MrBago

MrBago Sep 25, 2014

Contributor

_initialize cannot be called until after values/indices and such are computed. Currently those values are computed by "peaks_from_model" and attached to an instance of PeaksAndMetrics. It's common and acceptable for subclasses of a type to use the "internal" methods of their parent classes. Also because peaks_from_model is a de facto constructor for this type, i think it's also reasonalbe to call this method in that function.

That being said I'm open to redesigning peaks_from_model, I think it could use some cleanup, but I don't think we should do it as part of this PR.

This comment has been minimized.

@stefanv

stefanv Sep 25, 2014

Contributor

The reason this caught my attention is because it is being used outside the class (e.g. in peaks). If we can avoid such external calls, then it should be fine as-is.

This comment has been minimized.

@MrBago

MrBago Sep 25, 2014

Contributor

Ya, I know what you mean, but it really needs to be called by the
constructor, and at some point we decided that peaks was going to be the
constructor for this class.

On Wed, Sep 24, 2014 at 6:31 PM, Stefan van der Walt <
notifications@github.com> wrote:

In dipy/reconst/_peakdg.pyx:

+cdef class PAMDirectionGetter(DirectionGetter):
+

  • cdef:
  •    public double qa_thr, ang_thr, total_weight
    
  •    public double[:, :, :, ::1] _qa, _ind
    
  •    public double[:, ::1] _odf_vertices
    
  •    int initialized
    
  • def cinit(self):
  •    self.qa_thr = 0.0239
    
  •    self.ang_thr = 60
    
  •    self.total_weight = .5
    
  • def _initialize(self):

The reason this caught my attention is because it is being used outside
the class (e.g. in peaks). If we can avoid such external calls, then it
should be fine as-is.


Reply to this email directly or view it on GitHub
https://github.com/nipy/dipy/pull/341/files#r18011018.

This comment has been minimized.

@Garyfallidis

Garyfallidis Sep 25, 2014

Member

All I want from you here is to remove the underscore if you keep the method as you use it now. You can also rename it to set_pointers or something like that too to be more different than init and cinit. But if you are calling this method from an external function or other class then this shouldn't be a private function. Doesn't make sense. Please correct.

This comment has been minimized.

@MrBago

MrBago Sep 25, 2014

Contributor

I had an idea about this, it seems a little hackie but it might be the best way to go. What if get_direction or initial_direction call _initialize the first time that one of those methods get used?

This comment has been minimized.

@Garyfallidis

Garyfallidis Sep 29, 2014

Member

Okay that would do the job. Go for it but profile it too please!

self.initialized = 1
def initial_direction(self, double[::1] point):

This comment has been minimized.

@stefanv

stefanv Sep 22, 2014

Contributor

Docstring explaining methodology

@cython.boundscheck(False)
@cython.wraparound(False)
cpdef int get_direction(self, double[::1] point, double[::1] direction):

This comment has been minimized.

@stefanv

stefanv Sep 22, 2014

Contributor

Document return values

class QballBaseModel(SphHarmModel, Cache):
"""To be subclasses by Qball type models."""

This comment has been minimized.

@stefanv

stefanv Sep 22, 2014

Contributor

subclassed

if __name__ == "__main__":
import nose

This comment has been minimized.

@stefanv

stefanv Sep 22, 2014

Contributor

Most of the code seems to use:

from numpy.testing import run_module_suite
run_module_suite()

Perhaps worth using the same everywhere?

odf[r] = 1
return odf
def getDirection(dg, point, dir):

This comment has been minimized.

@stefanv

stefanv Sep 22, 2014

Contributor

get_direction

def _asarray(cython_memview):
# TODO: figure out the best way to get an array from a memory view.
# `np.array(view)` works, but is quite slow. Views are also "array_like",

This comment has been minimized.

@stefanv

stefanv Sep 22, 2014

Contributor

Is np.asarray(view) also slow?

This comment has been minimized.

@MrBago

MrBago Sep 25, 2014

Contributor

yep, asarray is the same as array.

"""
@classmethod
def fromPmf(klass, pmf, sphere, max_angle):

This comment has been minimized.

@stefanv

stefanv Sep 22, 2014

Contributor

camelCase?

@stefanv

This comment has been minimized.

Contributor

stefanv commented Sep 22, 2014

This is a neat PR! I wish I was able to review it better, but with such large changes it feels almost impossible to do well--can we try to split up these things into smaller pieces in general? @matthew-brett, do you know of any progress on the Cython code-coverage front?

@MrBago

This comment has been minimized.

Contributor

MrBago commented Sep 22, 2014

@stefanv I really have tried to keep this RP as focused as I can. The bulk of the PR is the core tracking components, the cython (extension) types needed to make it work and a reference implementation of specific tracking types.

I agree with you that PRs of over 1000 lies are hard to review, and we seem to be getting them somewhat routinely in dipy. That's probably part of the reason that dipy code is not as well reviewed as I think it should be. In this specific case I'm not really sure what to do about it, it's hard to (re)design the whole tracking framework in bits and pieces.

@stefanv

This comment has been minimized.

Contributor

stefanv commented Sep 22, 2014

@MrBago I understand completely. I think the best we can do for big refactorings like this is to have a good idea of test coverage and to document as well as we are able to. Either way, I gave it a read-over and in general this looks like a good step in the right direction.

MrBago added some commits Mar 19, 2014

NF - started working on new framework for local tracking
    Added basic outline of tracking framework in dipy.tracking.localtrack
        including DirectionGetter, TissueClassifier some code to hold it all
        together.
    Added example direction getter for PeaksAndMetrics in dipy.reconst._peakdg
        so PeaksAndMetrics can inherit DirectionGetter attributes.
NF - added LocalTracking class
    This class pulls together everything needed for tracking including:
        direction getter, tissue classifier, seeds, affine, and more
from dipy.direction import ProbabilisticDirectionGetter
prob_dg = ProbabilisticDirectionGetter.from_shcoeff(csd_fit.shm_coeff,

This comment has been minimized.

@Garyfallidis

Garyfallidis Nov 18, 2014

Member

Hey @MrBago! Shouldn't from_shcoeff have a sh basis parameter too?

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Nov 25, 2014

Hey all,

This PR looks ready to be merged on my side. As many have already reviewed the code I think it is safe to say that if nobody has any more comments to address I will merge this by tomorrow night (8 pm EST). 24 hours from now.

This is an important PR as it allows us to create many modular tracking algorithms and hopefully have the best tracking algorithm in the world! So, start using it!!! And play with the different parameters!!

@MrBago MrBago force-pushed the MrBago:tracking_WIP branch from e5f771a to cc87e5d Nov 25, 2014

@mdesco

This comment has been minimized.

Contributor

mdesco commented Nov 25, 2014

I have just spent time reading the PR. Great job and super nice API. This will really allow us to great tracking.

My main concern is using GFA computed from CSD in the examples added. I will point out a few areas where I would not use GFA. GFA from CSD fODF does not work well. The standard deviation of a sharp fODF is too high and the resulting GFA, not useful. Hence, I would not use CSD gfa nor define a classifier based on a thresholded GFA. I recommend using the standard FA or thresholded peaks amplitude as done in MRtrix or in Girard etal NeuroImage 2014. This will avoid confusion in the future and people mis-using GFA with CSD.

@MrBago do you think you can update the examples with that? Also, you might want to refine your CC mask slightly not to grab the fornix, the small structure underneath the CC.

Otherwise, examples are clear and very educational.

Great job!

from dipy.tracking.local import ThresholdTissueClassifier
classifier = ThresholdTissueClassifier(csa_peaks.gfa, .25)

This comment has been minimized.

@mdesco

mdesco Nov 25, 2014

Contributor

GFA on a CSA ODF makes sense and is ok.

probabilistic model of the corpus callosum.
"""
classifier = ThresholdTissueClassifier(csd_fit.gfa, .25)

This comment has been minimized.

@mdesco

mdesco Nov 25, 2014

Contributor

Here. I would not use the GFA from a CSD fODF. You could use a peak_amplitude threshold of 0.1 for instance as in MRtrix or Girard et al NeuroImage 2014. Or simple FA threshold.

I don't think we want Dipy users to start computing GFA from fODFs.

csd_model = ConstrainedSphericalDeconvModel(gtab, None, sh_order=6)
csd_fit = csd_model.fit(data, mask=white_matter)
classifier = ThresholdTissueClassifier(csd_fit.gfa, .25)

This comment has been minimized.

@mdesco

mdesco Nov 25, 2014

Contributor

Same comment as above.

@MrBago

This comment has been minimized.

Contributor

MrBago commented Nov 26, 2014

I agree, here we generally use the directions from CSA or CSD for tracking
but use FA for the tissue classifier. I didn't really want to go into that
for the example, but we also probably shouldn't make examples that are
going to develop bad habits in our users. I'll look over the examples and
try and figure out the best way to fix it.

Bago

On Tue, Nov 25, 2014 at 12:47 PM, Maxime Descoteaux <
notifications@github.com> wrote:

I have just spent time reading the PR. Great job and super nice API. This
will really allow us to great tracking.

My main concern is using GFA computed from CSD in the examples added. I
will point out a few areas where I would not use GFA. GFA from CSD fODF
does not work well. The standard deviation of a sharp fODF is too high and
the resulting GFA, not useful. Hence, I would not use CSD gfa nor define a
classifier based on a thresholded GFA. I recommend using the standard FA or
thresholded peaks amplitude as done in MRtrix or in Girard etal NeuroImage
2014. This will avoid confusion in the future and people mis-using GFA with
CSD.

@MrBago https://github.com/MrBago do you think you can update the
examples with that? Also, you might want to refine your CC mask slightly
not to grab the fornix, the small structure underneath the CC.

Otherwise, examples are clear and very educational.

Great job!


Reply to this email directly or view it on GitHub
#341 (comment).

@arokem

This comment has been minimized.

Member

arokem commented Nov 26, 2014

Take a look at the eudx_tracking: it gets the FA that was already computed
by reconst_dti.

On Tue, Nov 25, 2014 at 4:08 PM, MrBago notifications@github.com wrote:

I agree, here we generally use the directions from CSA or CSD for tracking
but use FA for the tissue classifier. I didn't really want to go into that
for the example, but we also probably shouldn't make examples that are
going to develop bad habits in our users. I'll look over the examples and
try and figure out the best way to fix it.

Bago

On Tue, Nov 25, 2014 at 12:47 PM, Maxime Descoteaux <
notifications@github.com> wrote:

I have just spent time reading the PR. Great job and super nice API.
This
will really allow us to great tracking.

My main concern is using GFA computed from CSD in the examples added. I
will point out a few areas where I would not use GFA. GFA from CSD fODF
does not work well. The standard deviation of a sharp fODF is too high
and
the resulting GFA, not useful. Hence, I would not use CSD gfa nor define
a
classifier based on a thresholded GFA. I recommend using the standard FA
or
thresholded peaks amplitude as done in MRtrix or in Girard etal
NeuroImage
2014. This will avoid confusion in the future and people mis-using GFA
with
CSD.

@MrBago https://github.com/MrBago do you think you can update the
examples with that? Also, you might want to refine your CC mask slightly
not to grab the fornix, the small structure underneath the CC.

Otherwise, examples are clear and very educational.

Great job!


Reply to this email directly or view it on GitHub
#341 (comment).


Reply to this email directly or view it on GitHub
#341 (comment).

@MrBago

This comment has been minimized.

Contributor

MrBago commented Nov 26, 2014

thanks, I'll take a look.
Bago

On Tue, Nov 25, 2014 at 4:10 PM, Ariel Rokem notifications@github.com
wrote:

Take a look at the eudx_tracking: it gets the FA that was already computed
by reconst_dti.

On Tue, Nov 25, 2014 at 4:08 PM, MrBago notifications@github.com wrote:

I agree, here we generally use the directions from CSA or CSD for
tracking
but use FA for the tissue classifier. I didn't really want to go into
that
for the example, but we also probably shouldn't make examples that are
going to develop bad habits in our users. I'll look over the examples
and
try and figure out the best way to fix it.

Bago

On Tue, Nov 25, 2014 at 12:47 PM, Maxime Descoteaux <
notifications@github.com> wrote:

I have just spent time reading the PR. Great job and super nice API.
This
will really allow us to great tracking.

My main concern is using GFA computed from CSD in the examples added.
I
will point out a few areas where I would not use GFA. GFA from CSD
fODF
does not work well. The standard deviation of a sharp fODF is too high
and
the resulting GFA, not useful. Hence, I would not use CSD gfa nor
define
a
classifier based on a thresholded GFA. I recommend using the standard
FA
or
thresholded peaks amplitude as done in MRtrix or in Girard etal
NeuroImage
2014. This will avoid confusion in the future and people mis-using GFA
with
CSD.

@MrBago https://github.com/MrBago do you think you can update the
examples with that? Also, you might want to refine your CC mask
slightly
not to grab the fornix, the small structure underneath the CC.

Otherwise, examples are clear and very educational.

Great job!


Reply to this email directly or view it on GitHub
#341 (comment).


Reply to this email directly or view it on GitHub
#341 (comment).


Reply to this email directly or view it on GitHub
#341 (comment).

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Nov 28, 2014

Hi @MrBago, I would like to move on with the release and this is a crucial PR. Can you fix the CSD/FA issue and then I will gladly go ahead and merge your PR. Keep it up and crack on!!! After that is merged Gab will submit a couple of some small PRs and I believe that will be the end of the tracking PRs for this release. We should also have a short meeting to discuss future actions for the tracking. Maybe next week will be great for that. Let us know when you can be available. -> +Inf !!!!

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Dec 2, 2014

Ping! Ping! There is no update for this PR for the last week although it is at its very last stages. Some feedback please. Thx in advance.

@arokem

This comment has been minimized.

Member

arokem commented Dec 5, 2014

Hey @Garyfallidis! Is this done?

to use the GFA of the CSD FODs to build a tissue classifier, however the GFA
values of these FODs don't classify gray matter and white matter well. We will
therefore use the GFA from the CSA model which we fit for the first section of
this example. Alternatively, one could fit a ``TensoModel`` to the data and use

This comment has been minimized.

@samuelstjean

samuelstjean Dec 5, 2014

Contributor

typo, TensoRModel

The distribution at each point is different and depends on the observed
diffusion data at that point. The distribution of tracking directions at each
point can be represented as a probability mass function (PMF) if the possible
tracking directions are restricted to discrete number of well distributed

This comment has been minimized.

@samuelstjean

samuelstjean Dec 5, 2014

Contributor

discrete numberS

tracking directions are restricted to discrete number of well distributed
points on a sphere.
This example picks up where "introduction to basic tracking" leaves off.

This comment has been minimized.

@samuelstjean

samuelstjean Dec 5, 2014

Contributor

Introduction should be capital, and I'd say left off instead of leaves off, since the tutorial is supposedly finished when one reads this.

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Dec 8, 2014

@arokem and @MrBago I have made a small PR on Bago's side with the typo corrections. MrBago#11

If you merge that one. I will go ahead and merge this one.

@arokem

This comment has been minimized.

Member

arokem commented Dec 8, 2014

I can't merge into Bago's fork, so we'll have to wait for him to merge that.

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Dec 8, 2014

Yes, of course.

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Dec 8, 2014

Okay @arokem just e-mailed me with an alternative suggestion. To merge this PR now and I can make another PR with the small typo changes in the main project. I think it is a good idea. So, let's move forward.

Garyfallidis added a commit that referenced this pull request Dec 8, 2014

@Garyfallidis Garyfallidis merged commit 5d542a6 into nipy:master Dec 8, 2014

1 check passed

continuous-integration/travis-ci The Travis CI build passed
Details
@arokem

This comment has been minimized.

Member

arokem commented Dec 8, 2014

Huzza!

On Mon, Dec 8, 2014 at 7:48 AM, Eleftherios Garyfallidis <
notifications@github.com> wrote:

Merged #341 #341.


Reply to this email directly or view it on GitHub
#341 (comment).

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment