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

MAPMRI #640

Merged
merged 34 commits into from Oct 12, 2015

Conversation

Projects
None yet
8 participants
@maurozucchelli
Contributor

maurozucchelli commented May 1, 2015

Hi everybody, this pull request add the minimal MAPMRI bases for signal reconstruction and ODF estimation. The fit works using cvxopt as an optional package as in shore.py, but since in this method the regularization of the basis is just empirical, adding positivity constraints should be strongly recommended.

After this pull request I will add all the other indexes proposed in Ozarslan et al. Neuroimage 2013 like:
RTOP
RTAP
RTPP
NG
PA
plus a proper regularization.

Thank you, best wishes

Mauro

@arokem

This comment has been minimized.

Member

arokem commented May 2, 2015

Cool! A couple of comments from a general overview, before I go into details:

  • There's a lot of unnecessary white-space in this code. Some functions have an empty line every other line. It would be easier to read the code with less of that.
  • The MapmriModel and MapmriFit need to inherit from reconst.base.ReconstModel and reconst.base.ReconstFit
  • The MapmriFit class needs to have a predict method, which takes a GradientTable as input and produces a prediction as output (for an example, see here: https://github.com/nipy/dipy/blob/master/dipy/reconst/dti.py#L1105). I I understand correctly what is going on here, it should be easy to implement (it's just an inversion of the fitting procedure...), and we need that in order to check whether the model fits the data.
the continuous functions presented in [2]_ but extending it in three
dimension.
The main difference with the SHORE proposed in [3]_ is that MAPMRI 3D
extension is provided using a set of three basis function for the radial

This comment has been minimized.

@Garyfallidis

Garyfallidis May 4, 2015

Member

functions

with respect to the MAPMRI model and compute the real and analytical
ODF.
from dipy.data import get_data,get_sphere

This comment has been minimized.

@Garyfallidis

Garyfallidis May 4, 2015

Member

When you add code in the example section of the docstring you have to start with

>>>

and then you can test if the example is running correctly using

nosetests --with-doctest
@multi_voxel_fit
def fit(self, data):
ind=self.gtab.bvals<=self.bmax_threshold

This comment has been minimized.

@Garyfallidis
ind=self.gtab.bvals<=self.bmax_threshold
gtab2 = gradient_table(self.gtab.bvals[ind], self.gtab.bvecs[ind,:])
tenmodel=dti.TensorModel(gtab2)
tenfit=tenmodel.fit(data[...,ind])

This comment has been minimized.

@Garyfallidis
ind_evals = np.argsort(evals)[::-1]
evals = evals[ind_evals]
R = R[ind_evals,:]

This comment has been minimized.

@Garyfallidis

Garyfallidis May 4, 2015

Member

pep8

Okay you have many pep8 issues. Use spyder or another editor that can automatically tell you were those pep8 are and please correct everywhere that is needed.

evals = np.clip(evals,self.eigenvalue_threshold,evals.max())
if self.anisotropic_scaling:
mu = np.sqrt(evals*2*self.tau)

This comment has been minimized.

@Garyfallidis
mu = np.sqrt(evals*2*self.tau)
else:
mumean=np.sqrt(evals.mean()*2*self.tau)

This comment has been minimized.

@Garyfallidis
else:
mumean=np.sqrt(evals.mean()*2*self.tau)
mu=np.array([mumean,mumean,mumean])

This comment has been minimized.

@Garyfallidis

Garyfallidis May 4, 2015

Member

Okay you need to add spaces between variables (see pep8). Will not make another comment about pep 8 in this PR. Here are some links to refresh your memory

https://github.com/numpy/numpy/blob/master/doc/HOWTO_DOCUMENT.rst.txt
https://www.python.org/dev/peps/pep-0008/

from dipy.reconst.mapmri import MapmriModel
from dipy.viz import fvtk
from dipy.data import fetch_isbi2013_2shell, read_isbi2013_2shell, get_sphere

This comment has been minimized.

@Garyfallidis

Garyfallidis May 4, 2015

Member

I would suggest to use a real dataset with this tutorial. It will attract more attention.

This comment has been minimized.

@Garyfallidis

Garyfallidis May 4, 2015

Member

Try this for example fetch_sherbrooke_3shell()

@arokem

This comment has been minimized.

Member

arokem commented May 11, 2015

Hey @maurozucchelli - have you had a chance to take a look at the comments on this?

@maurozucchelli

This comment has been minimized.

Contributor

maurozucchelli commented May 12, 2015

Yes, sorry for the delay! I almost finish, I just had some trouble with the
example in the code (the one with the >>>), which gave me some error when I
was trying to fetch the gradient table: apparently you cannot fetch in
these examples because it returns a warning when it has finished which is
considered an error by nosetest.

On Mon, May 11, 2015 at 11:00 PM, Ariel Rokem notifications@github.com
wrote:

Hey @maurozucchelli https://github.com/maurozucchelli - have you had a
chance to take a look at the comments on this?


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

@arokem

This comment has been minimized.

Member

arokem commented May 12, 2015

Hey Mauro - what's the warning you were seeing? Maybe we should suppress it.

On Tue, May 12, 2015 at 12:26 AM, maurozucchelli notifications@github.com
wrote:

Yes, sorry for the delay! I almost finish, I just had some trouble with the
example in the code (the one with the >>>), which gave me some error when I
was trying to fetch the gradient table: apparently you cannot fetch in
these examples because it returns a warning when it has finished which is
considered an error by nosetest.

On Mon, May 11, 2015 at 11:00 PM, Ariel Rokem notifications@github.com
wrote:

Hey @maurozucchelli https://github.com/maurozucchelli - have you had a
chance to take a look at the comments on this?


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


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

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented May 12, 2015

No, I don't think we can. He is talking about the message saying that you have already fetched the data. Fetching in a docstring example I think is not a good idea.

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented May 12, 2015

Maybe just don't have an example in the docstring if it creates problems.

@MrBago

This comment has been minimized.

Contributor

MrBago commented May 13, 2015

I think we agreed some time ago that tests should be runnable without network connectivity and without requiring data be pre-fetched. Fetching data in a doc-test violates that rule. Is there a way to leave the example but to tell nose that it should bot be run as a doc-test? Or maybe just replace it with a pointer to the full example, which is much better anyway.

M = mapmri_phi_matrix(self.radial_order, mu, q.T)
ind_mat = mapmri_index_matrix(self.radial_order)

This comment has been minimized.

@MrBago

MrBago May 13, 2015

Contributor

There are a few things here that could be moved outside of the fit function. For example ind_mat and tenmodel don't need to be re-built at each voxel. Instead you could make them in the init and re-use them.

It might make your fit method easier to follow because currently it's kind of long. Having only the steps that are specific to individual voxels might make the method easier to understand.

v = np.dot(v_,self.R)
I_s = mapmri_odf_matrix(self.radial_order,self.mu, s, v)
odf = np.dot(I_s, self._mapmri_coef)
return np.clip(odf,0,odf.max())

This comment has been minimized.

@MrBago

MrBago May 13, 2015

Contributor

clip(odf, 0, None) will be more efficient than using odf.max().

@maurozucchelli

This comment has been minimized.

Contributor

maurozucchelli commented May 15, 2015

Hi, in the example the ODF in the Corpus Callosum are oriented in a strange way! I run the SHORE on it and I get the same result... Do you think that is just a visualization problem?

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented May 15, 2015

@maurozucchelli probably yes. It 's a visualization thing. Or the bvecs need flipping.

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented May 15, 2015

There is a failure in python 2.7 according to Travis.

r"""Mean Apparent Propagator MRI (MAPMRI) [1]_ of the diffusion signal.
The main idea is to model the diffusion signal as a linear combination of
the continuous function presented in [2]_ but extending it in three

This comment has been minimized.

@Garyfallidis

Garyfallidis May 15, 2015

Member

functions

The main idea is to model the diffusion signal as a linear combination of
the continuous function presented in [2]_ but extending it in three
dimension.

This comment has been minimized.

@Garyfallidis

Garyfallidis May 15, 2015

Member

dimensions

the continuous function presented in [2]_ but extending it in three
dimension.
The main difference with the SHORE proposed in [3]_ is that MAPMRI 3D
extension is provided using a set of three basis function for the radial

This comment has been minimized.

@Garyfallidis

Garyfallidis May 15, 2015

Member

functions

uses one basis function to model the radial part and real Spherical
Harmonics to model the angular part.
From the MAPMRI coefficients is possible to use the analytical formulae
to estimate the ODF.

This comment has been minimized.

@Garyfallidis

Garyfallidis May 15, 2015

Member

Correct typos as above.

>>> data, gtab = dsi_voxels()
>>> sphere = get_sphere('symmetric724')
>>> from dipy.sims.voxel import SticksAndBall
>>> data, golden_directions = SticksAndBall(gtab, d=0.0015, S0=1, angles=[(0, 0), (90, 0)], fractions=[50, 50], snr=None)

This comment has been minimized.

@Garyfallidis

Garyfallidis May 15, 2015

Member

I think it should be possible to brake this line so that it doesn't take so much space.

This comment has been minimized.

@maurozucchelli

maurozucchelli May 18, 2015

Contributor

I tried but I couldn't get the indentation right, it always gave me some error.

Parameters
----------
s: unsigned int

This comment has been minimized.

@Garyfallidis

Garyfallidis May 15, 2015

Member

Add a space in front of parameters

s : unsigned int

This comment has been minimized.

@Garyfallidis

Garyfallidis May 15, 2015

Member

Correct this everywhere please.

Returns
-------
B: array, shape (N,)

This comment has been minimized.

@Garyfallidis

Garyfallidis May 15, 2015

Member

Add a space in front of parameters

B : array, shape (N, )

Please correct everywhere in the PR.

--------
In this example, where the data, gradient table and sphere tessellation
used for reconstruction are provided, we model the diffusion signal
with respect to the MAPMRI model and compute the real and analytical

This comment has been minimized.

@Garyfallidis

Garyfallidis May 15, 2015

Member

You say here that you compute the real ODF but you only compute the analytical ODF. Except if you mean something else.

M = mapmri_phi_matrix(self.radial_order, mu, q.T)
# This is a simple empirical regularization, to be replaced
I = np.diag(self.ind_mat.sum(1)**2)

This comment has been minimized.

@Garyfallidis

@maurozucchelli maurozucchelli force-pushed the maurozucchelli:mapmri branch from 4bfb0c2 to e015933 Sep 14, 2015

maurozucchelli
@maurozucchelli

This comment has been minimized.

Contributor

maurozucchelli commented Sep 14, 2015

Ok done.... Let's see what Travis will say about it!

On Mon, Sep 14, 2015 at 4:38 PM, Eleftherios Garyfallidis <
notifications@github.com> wrote:

Hi @maurozucchelli https://github.com/maurozucchelli the interp_rbf
problem should be solved now. Can you rebase again?


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

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Sep 14, 2015

Okay things look much better now. You have only this bug to figure out. Which happens in only one of the Travis bots.
ValueError: Buffer dtype mismatch, expected 'double' but got Python object
Happy hacking!

@maurozucchelli

This comment has been minimized.

Contributor

maurozucchelli commented Sep 15, 2015

Finally Travis accepted it! Apparently clipping the ODF was not very well accepted in a certain version, I removed the clipping and now it is going...

@arokem

This comment has been minimized.

Member

arokem commented Oct 6, 2015

Sorry for the delay here! I see that my comments from back in May are still unaddressed. Any chance to address these comments?

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Oct 7, 2015

@maurozucchelli have you seen @arokem's comments? Let's push forward please.

@maurozucchelli

This comment has been minimized.

Contributor

maurozucchelli commented Oct 7, 2015

Yes but sadly I'm very busy this week, I will probably fix the code this
week-end or Monday.

I thought that I had already addressed everything, what did I miss?

On Wed, Oct 7, 2015 at 2:47 AM, Eleftherios Garyfallidis <
notifications@github.com> wrote:

@maurozucchelli https://github.com/maurozucchelli have you seen @arokem
https://github.com/arokem's comments? Let's push forward please.


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

@arokem

This comment has been minimized.

Member

arokem commented Oct 7, 2015

See comments here: #640 (comment)

@arokem

This comment has been minimized.

Member

arokem commented Oct 7, 2015

Oops - sorry - I see that you have indeed addressed these comments - my bad.

@arokem

This comment has been minimized.

Member

arokem commented Oct 7, 2015

Is this ready to go, otherwise?

@arokem

This comment has been minimized.

Member

arokem commented Oct 7, 2015

OK - the main issue I still see here is test coverage. There are a few other comments here and there otherwise. We also need to think how we deal with the license issue - this code would be GPL licensed. Maybe separate it off into another repo, with its own license, instead of including it as part of the dipy repo? I am happy to help set that up.

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Oct 7, 2015

We have already an explanation for cvxopt in our README file. So, there is no reason to do this now. Let's move on please. But yeah we need to revisit the license issue at a later stage in a more complete way.

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Oct 7, 2015

@maurozucchelli after correcting the few minor things make the next PR with the maps. Mapmri is not so useful without the other maps.

@arokem

This comment has been minimized.

Member

arokem commented Oct 7, 2015

At the very least, something like this needs to be introduced here as well:
https://github.com/nipy/dipy/blob/master/dipy/reconst/shore.py#L234

On Wed, Oct 7, 2015 at 7:58 AM, Eleftherios Garyfallidis <
notifications@github.com> wrote:

@maurozucchelli https://github.com/maurozucchelli after correcting the
few minor things make the next PR with the maps. Mapmri is not so useful
without the other maps.


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

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Oct 7, 2015

Sure.

maurozucchelli
@maurozucchelli

This comment has been minimized.

Contributor

maurozucchelli commented Oct 12, 2015

Ok I added the warning for the cvxopt usage. Is it ok now?

On Wed, Oct 7, 2015 at 5:00 PM, Eleftherios Garyfallidis <
notifications@github.com> wrote:

Sure.


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

maurozucchelli
@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Oct 12, 2015

Good job! @maurozucchelli go ahead with the next PR with the maps please to make this work complete!
We take it that the pinv question will be answered there as you said here.

Garyfallidis added a commit that referenced this pull request Oct 12, 2015

@Garyfallidis Garyfallidis merged commit 11c9e88 into nipy:master Oct 12, 2015

1 check passed

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

This comment has been minimized.

Member

Garyfallidis commented Oct 12, 2015

Thx!

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