RF - Separate direction finder from model fit.

[stefanv]

Here's a concept proposal for pulling out the direction finder from the fit. Suggestions welcome. All tests pass, except for the peaks / indices in GQI / GSI which I had to disable temporarily.

[MrBago]

@stefanv, I can't tell which of the changes here are substantive and which are white-space/file-renames. I'll take a look at this later, but it's going to take me some time to sift though it.

[Garyfallidis]

@MrBago just look at peaks.py and the tests.

[MrBago]

Also a quick note about the direction finder. The purpose of the direction finder, in my eyes at least, is to unify the api for odf-based models and the TensorModel or other tensor based models. By having a unified interface, the user of a fit object doesn't need to worry about whether that fit object exposes it's directions though fit.evecs or directionfinder(fit.odf). Maybe this unification of the api is more trouble then it's worth.

[Garyfallidis]

I think the hole trouble starts even before the direction finder. The trouble starts with the enforcement of the sphere as input to the odf for all odf like models. This is what made designing difficult at the first place. I want you to see the entire picture of the reconstruction problem.

[stefanv]

It seems reasonable, for the sake of consistency, to use the same method of direction finding across all models (including tensors, even when you know exactly what direction should be).

I think we can address @Garyfallidis 's concern by changing the signature of ODF to odf(sphere=default_sphere) (and probably the same for GFA, etc.) When the ODF is requested, we cache it, so that when GFA etc. has to be computed it's available immediately.

[Garyfallidis]

I think this is a good idea. Let me summarize what we discussed on the phone. a) The ReconstructionModel can have as input parameter a sphere if it is needed. b) The odf() will be able to use that sphere if no other sphere is requested. c) The discrete direction finder will also return peak_values and peak_indices as these are required for most tracking algorithms and EuDX.

So, I suppose if we do this then the direction finder separation is not anymore required. Correct @stefanv ?

[MrBago]

Look guys I'm at the point where I'm willing to accept almost anything so we can more past this, but I'd like to point out that you've designed a fit object which is highly optimized for EuDx at some cost to my use-case. I disagree with the statement "peak_values and peak_indices ... are required for most tracking algorithms", but if this is the only use-case you are interested in supporting, we can safely move forward with this design.

[Garyfallidis]

I don't think there will be any cost for your use-case. It should work smoothly with your tracking algorithms you will only use directions, not maxima or indices. Actually, you are correct I haven't got any statistics about how many algorithms discretize their directions. However, DSI Studio, Camino, EuDX and all tractography algorithms developed by Maxime and his collaborators do that. They even made an algorithm that does real time tracking which the surgeon uses in the operation room (incredible?) which again discretizes directions. So, it seems it is more common than what we initially thought. I was surprised too. It makes sense though as you seriously reduce memory size.

Concerning the new DiscreteDirectionFinder with the additional attributes of directions, maxima, indices, we can cover all ODFmodels. For the Tensor case we might need to derive a special DiscreteDirectionFinder i.e TensorDDF which for example returns maxima as the evals etc. We could even have an option to use either the ODFDDF or the TensorDDF as they both make sense.

It is true that if we implement this options it will be very easy to use EuDX with any model which at the moment is quite complicated. By doing that we will at least have some tracking algorithm doing something more or less useful but this adds anyway value to dipy.

Lets move forward then.

[Garyfallidis]

@stefanv and @MrBago, is it okay if I close this PR and PR #97? We can discuss more about this when we will have a better overview of fiber tracking needs. @gabknight has some cool tracking to share. He works together with @mdesco. We could sit after the release and see what extra needs we have to include both Gab's and Bago's methods or even better share code between the two approaches. For the time being I can use peaks_from_models for the examples. I hope this is okay.

[stefanv]

If you don't like these proposals on technical merits, then we should discuss closing, otherwise I'd prefer to find agreement on a way forward. I am concerned about releasing dipy as-is, because once we have put a certain API in place in a release and people start using it, it will be hard to back out. With these two PRs, I think we get closer to the desired outcome, so I'd gladly hear your feedback.

[stefanv]

My gut feel is this: keep all classes minimal. Ideally, the model would store only the model parameters and fitting, the fit would store only fit information and be able to reproduce the ODF and signal. External utility functions operate on these classes.

We have to compromise sometimes for performance, but we keep those compromises to a minimum. E.g., we cache some computations made during fitting on the parent model, so that we don't have to repeat that computation later. This is a hack and ultimately undesirable, but worth it.

For this release, it is feasible to polish the API for single voxel access. We can generalize from single to multi-voxel via utility functions, and also provide utilities for easily looping over a voxel-data/mask combination. What I don't think is achievable is to optimize for both a clean, general, multi-voxel API and speed and memory usage at the same time.

As a controversial illustration of this principle, I think that multi_voxel_model_fit.gfa() should not be provided. Rather, one can compute it with the explicit for (i, j, k), v in voxel(data, mask): gfa[i, j, k] = gfa(model.fit(v.data), sphere) (or similar utility function to broadcast multi-voxel operations).

[Garyfallidis]

Hi @stefanv , sorry for not having been precise in my previous statement. I agree in general with what you suggest above. I thought it would be clear from the example I added in PR #103. I am taking back my initial idea to support caching of gfas etc inside of odfs for dsi-like models. I will have a PR with clean ups asap. Lets move on with what you and Bago have suggested i.e. leave processing of gfas and other properties for a later stage than doing that from inside multivoxel object methods. This is what I meant by saying that we can use peaks_from_model (in odf.py) because this is exactly a function like the gfa(model.fit...) that you propose above. And it can be used to create input variables that they can later be used with EuDX or other algorithms.

My point is: we keep the model/fit objects very simple and not worry about specific non-general properties (like gfa) at the model/fit level. Which I think is the same idea that you propose.

Now the question is do we still need directions to be a method of the single_voxel_modelfit? I think @MrBago suggested that it is not necessary any more and we can use peak_directions which is another utility function similar with peaks_from_model. What do you think @stefanv ?

With all these ideas accepted the only exception is the TensorModel where evecs(directions), and fa are methods of the fit. We could even separate those from the fit and provide utility functions that calculate those variables. That would make the API more general.

As you remember my fear was that we could drive our users crazy as we can write tenfit.fa but not csafit.gfa. Maybe, by taking these properties/attributes out of the tensor we make our API more general. What do you all think?

[MrBago]

I propose to drop the direction_finder on the odfFit, but keep the directions keyword on models that can support it (ie tensor like models).

I'd like to keep fa on the TensorFit for convince, but I've started to move the relevant code into independent functions so it'll be easy to drop if that what you guys think is best.

[Garyfallidis]

I agree with what you propose @MrBago

[stefanv]

Sorry for the slow response due to traveling.

I don't think models should provide directions--they can really only be provided for very few models, and even in those cases I wouldn't use them, because they give that model an unfair advantage in comparisons. If you really want the directions of a tensor fit, why not simply use the directions of the evecs?

[Garyfallidis]

Good point. But what if you use a multitensor or sticks&ball? You will use mevecs rather than evecs or just use directions?

[MrBago]

Do you have a stick&ball or multi-tensor model at the moment?

I think we'll need to write out the logic for going between evecs, mevecs, or sticks to directions (though in most cases the logic might be only a few lines) and because that logic is model specific, putting it in a method or propriety seems like it would be most appropriate. Thought that doesn't have to happen right now, we can punt it for later.

[Garyfallidis]

We have them only for simulations. I introduced them here as a thought experiment to help the discussion.

[MrBago]

I believe this pull request was superseded by #77.