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

Forecast #1322

Merged
merged 14 commits into from Jan 15, 2018

Conversation

Projects
6 participants
@maurozucchelli
Contributor

maurozucchelli commented Aug 23, 2017

Dear DIPY developers,
This pull request is made to add the FORECAST model to DIPY:
http://onlinelibrary.wiley.com/doi/10.1002/mrm.20667/full

FORECAST is a Spherical Deconvolution method that allows estimating voxel specific response function in the form of an axially symmetric tensor. From this tensor, FORECAST has been recently used for estimating crossing invariant diffusivities and fractional anisotropy in:
http://onlinelibrary.wiley.com/doi/10.1002/mrm.25734/full

I'm going to present the comparison between FORECAST and similar techniques at the 2017 MICCAI Computational Diffusion MRI workshop and this implementation is based on the one that I use for that work.

Best wishes,

Mauro

@codecov-io

This comment has been minimized.

codecov-io commented Aug 23, 2017

Codecov Report

Merging #1322 into master will increase coverage by 0.07%.
The diff coverage is 93.37%.

Impacted file tree graph

@@            Coverage Diff             @@
##           master    #1322      +/-   ##
==========================================
+ Coverage   87.07%   87.14%   +0.07%     
==========================================
  Files         227      229       +2     
  Lines       29060    29377     +317     
  Branches     3123     3147      +24     
==========================================
+ Hits        25303    25602     +299     
- Misses       3047     3060      +13     
- Partials      710      715       +5
Impacted Files Coverage Δ
dipy/reconst/forecast.py 90.15% <90.15%> (ø)
dipy/reconst/tests/test_forecast.py 98.38% <98.38%> (ø)
dipy/reconst/csdeconv.py 89.32% <0%> (+0.64%) ⬆️
dipy/core/graph.py 75% <0%> (+1.19%) ⬆️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update ae69540...7641437. Read the comment docs.

@coveralls

This comment has been minimized.

coveralls commented Aug 23, 2017

Coverage Status

Coverage increased (+3.9%) to 89.217% when pulling f0c8138 on maurozucchelli:forecast into c268bd2 on nipy:master.

@arokem

This comment has been minimized.

Member

arokem commented Aug 23, 2017

Hi Mauro - this looks fantastic. Thanks for the contribution!

Before I go into a detailed review, I do have one initial comment: the current Fit method doesn't seem to have a predict method, used to assess the goodness of fit to the signal. Would it be possible to implement the predict method for this model?

@skoudoro

Hi @maurozucchelli , Thank you for this !

I think you should add your example in 2 places : -

  • doc/examples/valid_examples.txt
  • doc/examples_index.rst

After a quick test, everything seems to run/generate without any problems. I just added a couple of PEP8 comment before an algorithm review.

The implementation of FORECAST may require CVXPY (http://www.cvxpy.org/).
"""
def __init__(self,

This comment has been minimized.

@skoudoro

skoudoro Aug 23, 2017

Member

I think you need to initialize the parents class so it miss the call to init of super class. Something like super(ForecastModel, self).__init__(your_arguments).

This comment has been minimized.

@maurozucchelli

maurozucchelli Aug 24, 2017

Contributor

Hi, I followed SHORE and MAPMRI style. What is the difference between this and the simple init that I used?

This comment has been minimized.

@skoudoro

skoudoro Aug 24, 2017

Member

Ok, I see, there is the same mistake on SHORE and MAPMRI :-). ReconstModel should initialize some variable (like gtab). You can see some examples on FWDTI, LIFE or SFM.

You should keep your __init__, this one initialize your class ForecastModel. You just have to initialize ReconstModel since you inherit from this class.

def __init__(self, gtab, sh_order=8, lambda_lb=1e-3,
             optimizer='wls', sphere=None, lambda_csd=1.0)

        # your docstrings
        super(ForecastModel, self).__init__(gtab)  # or ReconstModel.__init__(gtab) 
from dipy.reconst.cache import Cache
from dipy.reconst.multi_voxel import multi_voxel_fit
from dipy.reconst.shm import real_sph_harm
from dipy.core.gradients import gradient_table

This comment has been minimized.

@skoudoro

skoudoro Aug 23, 2017

Member

unused import. I think you can remove it

from dipy.reconst.shm import real_sph_harm
from dipy.core.gradients import gradient_table
from scipy.special import erf
from math import factorial

This comment has been minimized.

@skoudoro

skoudoro Aug 23, 2017

Member

unused import. I think you can remove it

self.bvals = np.round(gtab.bvals/100) * 100
self.bvecs = gtab.bvecs
self.gtab = gtab
if sh_order >= 0 and not(bool(sh_order % 2)) and sh_order<=12:

This comment has been minimized.

@skoudoro

skoudoro Aug 23, 2017

Member

PEP8: missing space around operator <=

data_b0 = data[self.b0s_mask].mean()
data_single_b0 = np.r_[data_b0, data[~self.b0s_mask]] / data_b0
#data_single_b0 = np.clip(data_single_b0, 0, 1.0)

This comment has been minimized.

@skoudoro

skoudoro Aug 23, 2017

Member

PEP8: missing space after block comment

from dipy.reconst.forecast import ForecastModel
from dipy.reconst.shm import sh_to_sf
from dipy.viz import fvtk

This comment has been minimized.

@skoudoro

skoudoro Aug 23, 2017

Member

unused import

from dipy.viz import fvtk
# from dipy.data import fetch_sherbrooke_3shell, read_sherbrooke_3shell, get_sphere
from dipy.data import fetch_cenir_multib, read_cenir_multib, get_sphere
from dipy.core.gradients import gradient_table

This comment has been minimized.

@skoudoro

skoudoro Aug 23, 2017

Member

unused import

mask_small = data_small[..., 0] > 1000
"""
Instantiate the FORECAST Model.

This comment has been minimized.

@skoudoro

skoudoro Aug 23, 2017

Member

missing 1 blank line

Let us consider only a single slice for the FORECAST fitting
"""
data_small = data[18:87,51:52,10:70]

This comment has been minimized.

@skoudoro

skoudoro Aug 23, 2017

Member

PEP8: missing space after ,

Display a part of the fODFs
"""
from dipy.viz import fvtk

This comment has been minimized.

@skoudoro

skoudoro Aug 23, 2017

Member

maybe you should remove this import and keep the one in the top of the file

>>> gtab = get_3shell_gtab()
>>> from dipy.sims.voxel import MultiTensor
>>> mevals = np.array(([0.0017, 0.0003, 0.0003],
>>> [0.0017, 0.0003, 0.0003]))

This comment has been minimized.

@skoudoro

skoudoro Aug 24, 2017

Member

doctest failed, I think it should be ... [0.0017, 0.0003, 0.0003])) here.

>>> angl = [(0, 0), (60, 0)]
>>> data, sticks = MultiTensor(
>>> gtab, mevals, S0=100.0, angles=angl,
>>> fractions=[50, 50], snr=None)

This comment has been minimized.

@skoudoro

skoudoro Aug 24, 2017

Member

same as above

@coveralls

This comment has been minimized.

coveralls commented Aug 25, 2017

Coverage Status

Coverage increased (+0.09%) to 85.427% when pulling 66aa4a4 on maurozucchelli:forecast into c268bd2 on nipy:master.

1 similar comment
@coveralls

This comment has been minimized.

coveralls commented Aug 25, 2017

Coverage Status

Coverage increased (+0.09%) to 85.427% when pulling 66aa4a4 on maurozucchelli:forecast into c268bd2 on nipy:master.

@coveralls

This comment has been minimized.

coveralls commented Aug 25, 2017

Coverage Status

Coverage increased (+0.09%) to 85.427% when pulling 2e6afee on maurozucchelli:forecast into c268bd2 on nipy:master.

1 similar comment
@coveralls

This comment has been minimized.

coveralls commented Aug 25, 2017

Coverage Status

Coverage increased (+0.09%) to 85.427% when pulling 2e6afee on maurozucchelli:forecast into c268bd2 on nipy:master.

@maurozucchelli

This comment has been minimized.

Contributor

maurozucchelli commented Aug 25, 2017

I got an error in travis for "test_csdeconv.py" but I did not modify that file, nor any other files except the FORECAST ones. What should I do with this error?

@arokem

This comment has been minimized.

Member

arokem commented Aug 25, 2017

That error is on the PRE=1 machine. That means that it's detecting errors that will start appearing in upcoming releases of one of our dependencies (in this case, the upcoming numpy 1.14). We should resolve that in a separate PR. I posted an issue here: #1323

@maurozucchelli

This comment has been minimized.

Contributor

maurozucchelli commented Sep 20, 2017

I have updated the branch, merging nipy/dipy master in order to solve the CSD error in Travis. Now everything should be ok

@skoudoro

This comment has been minimized.

Member

skoudoro commented Sep 20, 2017

Thanks for this update @maurozucchelli. New feature for the release in October !

We just need to complete the review and it might be good to have other suggestions. Can you have a look @arokem @Garyfallidis or anyone else ?

@arokem

This comment has been minimized.

Member

arokem commented Sep 20, 2017

Running the example, I get a UserWarning: CSD convergence not reached. Is there some way to avoid this warning -- at least in the example we use? Maybe use a different selection from the data?

import matplotlib.pyplot as plt
from dipy.reconst.forecast import ForecastModel
from dipy.viz import fvtk

This comment has been minimized.

@skoudoro

skoudoro Sep 20, 2017

Member

fvtk will be deprecated in the near future. Can you replace it by the new API ? You can look at the Tracking Quick Start example

This comment has been minimized.

@maurozucchelli

maurozucchelli Sep 25, 2017

Contributor

I'm doing some tests. If I try to do something like:

ren = window.Renderer()
odf_actor = actor.odf_slicer(odf[16:36, :, 30:45], sphere = sphere, colormap='jet')
odf_actor.RotateX(-90)
ren.add(odf_actor)
window.record(ren, out_path='fODFs.png', size=(600, 600), magnification=1)
window.show(ren)

I have a problem: I can see only a line of overlapping ODFs and not a 2D grid of ODFs like while using fvtk. All the visualization parameters are the same for both functions (odf_slicer, and fvtk.sphere_func).

Is it correct what I'm doing? What is the problem here?

fodfs

@maurozucchelli

This comment has been minimized.

Contributor

maurozucchelli commented Sep 25, 2017

@arokem Hi! I replaced my CSD implementation with DIPY csdeconv but we still get the warning. I think that some voxels near the brain border are corrupted. The easiest solution is to mask the data, using median_otsu or some other techniques. What do you think?

@arokem

This comment has been minimized.

Member

arokem commented Sep 25, 2017

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Sep 26, 2017

Hi @maurozucchelli and thanks for the PR.
The odf_slicer works in this way (no need for initial cropping)

odf_actor = actor.odf_slicer(odf, sphere = sphere, colormap='jet')
odf_actor.display(y=20) # shows entire slice at Y=20
or odf_actor.display_extent( ...) to only show a patch of the slice

If your odf is too large for your memory. Just memory map it using np.memmap before giving it as input to odf_slicer.
In summary, the odf_slicer works like a regular slicer. Slices entire volumes. Let me know how it goes.

@skoudoro

This comment has been minimized.

Member

skoudoro commented Sep 27, 2017

Thanks @Garyfallidis ! it works well !

I just wonder why there is some red ball on this result as you can see on the top left in the image below ?

fodfs

odf_actor = actor.odf_slicer(odf, sphere = sphere, colormap='jet', scale=0.25)
odf_actor.display(y=0)
ren = window.Renderer()
ren.add(odf_actor)
window.record(ren, out_path='fODFs.png', size=(600, 600), magnification=1)
window.show(ren)
@maurozucchelli

This comment has been minimized.

Contributor

maurozucchelli commented Oct 6, 2017

Hello! I fixed the example as suggested by @Garyfallidis and @skoudoro and modified a bit the main code to improve stability.

@maurozucchelli

This comment has been minimized.

Contributor

maurozucchelli commented Oct 6, 2017

I got an error in some non-forecast related tests with travis using PRE=1, e.g. vec2vec_rotmat, which I didn't had before. Is there a way to fix these errors?

@skoudoro

This comment has been minimized.

Member

skoudoro commented Oct 6, 2017

we talked about these errors on issue #1346. I plan to fix it today and then, you can rebase.

Thanks for your fix ! I will review it ASAP

@maurozucchelli maurozucchelli force-pushed the maurozucchelli:forecast branch from a1ec448 to ee5d7bf Oct 16, 2017

@maurozucchelli

This comment has been minimized.

Contributor

maurozucchelli commented Oct 17, 2017

Hi I rebased! Now everything seems to work fine

@skoudoro skoudoro added this to Needs a review in Dipy 0.14.0 Nov 8, 2017

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Dec 1, 2017

Hi @maurozucchelli apologies for delaying to review this PR. See my comments. Also you will probably have to rebase again. Will see in your next commit in this PR. We have fixed some issues for upcoming numpy and there will be merging issues if you do not rebase. Thanks for the great PR. Waiting for the minor fixes!! :)

@maurozucchelli maurozucchelli force-pushed the maurozucchelli:forecast branch from ee5d7bf to 083b667 Dec 5, 2017

sh_order is the spherical harmonics order used for the fODF.
optimizer is the algorithm used for the FORECAST basis fitting, in this case

This comment has been minimized.

@Garyfallidis

Garyfallidis Dec 8, 2017

Member

CSD is not an optimizer. CSD is a model. Do you mean something else here @maurozucchelli ?
For example, do you mean the same optimization as it happened in CSD?

This comment has been minimized.

@maurozucchelli

maurozucchelli Dec 11, 2017

Contributor

Exaclty, the same iterative optimization algorithm employed in the CSD model can be used here. We can change the name to "deconvolution_algorithm" or something on this line.

"""
Instantiate the FORECAST Model.
sh_order is the spherical harmonics order used for the fODF.

This comment has been minimized.

@Garyfallidis

Garyfallidis Dec 8, 2017

Member

You may want to add quotes around sh_order.

class ForecastModel(OdfModel, Cache):
r"""Fiber ORientation Estimated using Continuous Axially Symmetric Tensors
(FORECAST) [1,2,3]_. FORECAST is a Spherical Deconvolution reconstruction model

This comment has been minimized.

@Garyfallidis

Garyfallidis Dec 8, 2017

Member

What is interesting in this model is that I don't see any response function. Why is that? Deconvolution without a response? Please explain.

This comment has been minimized.

@maurozucchelli

maurozucchelli Dec 11, 2017

Contributor

The response is estimated in each voxel independently using the mean of the signal on each shell to characterize the diffusivity parallel and perpendicular. This technique works fairly well but is quite sensitive to noise!

r"""Fiber ORientation Estimated using Continuous Axially Symmetric Tensors
(FORECAST) [1,2,3]_. FORECAST is a Spherical Deconvolution reconstruction model
for multi-shell diffusion data which enables the calculation of a voxel
adaptive response function using the Spherical Mean Tecnique (SMT) [2,3]_.

This comment has been minimized.

@Garyfallidis

Garyfallidis Dec 8, 2017

Member

Do we have a comparison anywhere against CSA and and CSD? Do we see important qualitative/quantitative or methodological differences?

This comment has been minimized.

@maurozucchelli

maurozucchelli Dec 11, 2017

Contributor

It can be seen exactly like a CSD with a voxel adaptive response function.

gtab,
sh_order=8,
lambda_lb=1e-3,
optimizer='WLS',

This comment has been minimized.

@Garyfallidis

Garyfallidis Dec 8, 2017

Member

Is the default optimizer the best option? Or the easiest?

This comment has been minimized.

@maurozucchelli

maurozucchelli Dec 11, 2017

Contributor

It is the fastest. I can make the CSD the default because it is more robust!

return v
def Psi_l(l,b):

This comment has been minimized.

@Garyfallidis

Garyfallidis Dec 8, 2017

Member

lower case in function names

def forecast_matrix(sh_order, d_par, d_perp, bvals):
r"""Compute the FORECAST radial matrix
"""
n_c = int((sh_order + 1)*(sh_order + 2)/2)

This comment has been minimized.

@Garyfallidis

Garyfallidis Dec 8, 2017

Member

okay... spaces around operators.

r, theta, phi = cart2sphere(vecs[:, 0], vecs[:, 1], vecs[:, 2])
theta[np.isnan(theta)] = 0
n_c = int((sh_order + 1)*(sh_order + 2)/2)

This comment has been minimized.

@Garyfallidis
def lb_forecast(sh_order):
r"""Returns the Laplace-Beltrami regularization matrix for FORECAST
"""

This comment has been minimized.

@Garyfallidis

Garyfallidis Dec 8, 2017

Member

nice! :)

Display a part of the fODFs
"""
odf_actor = actor.odf_slicer(odf[16:36, :, 30:45], sphere=sphere,

This comment has been minimized.

@Garyfallidis

Garyfallidis Dec 8, 2017

Member

If you want to show the entire slice and the ODF and you don't do that for memory issues you can just create a memmap (using np.memmap) of the odf and give that to the odf_slicer. I believe I have an example of that in the odf_slicer's tests.

This comment has been minimized.

@Garyfallidis

Garyfallidis Dec 8, 2017

Member

This is good work so let's make this as beautiful as possible.

Load an odf reconstruction sphere
"""
sphere = get_sphere('symmetric724')

This comment has been minimized.

@Garyfallidis

Garyfallidis Jan 15, 2018

Member

The sphere called 'repulsion724' is a better sphere. Better evenly distributed. Made by Dr. @mpaquette !

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Jan 15, 2018

OKay! I am merging this. In general., for next time if you want to visualize ODFs see the following lines for ideas.
https://github.com/nipy/dipy/blob/master/dipy/viz/tests/test_actors.py#L406
You can use a mask to mask out as large regions as you want. Or you can use display_extent to show a specific slice patch or both.

@Garyfallidis Garyfallidis merged commit 18fb8c2 into nipy:master Jan 15, 2018

3 checks passed

codecov/patch 93.37% of diff hit (target 87.07%)
Details
codecov/project 87.14% (+0.07%) compared to ae69540
Details
continuous-integration/travis-ci/pr The Travis CI build passed
Details

Dipy 0.14.0 automation moved this from Needs a review to Done Jan 15, 2018

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Jan 15, 2018

Congrats @maurozucchelli!

@maurozucchelli

This comment has been minimized.

Contributor

maurozucchelli commented Jan 16, 2018

Thank you @Garyfallidis! And also @arokem, @skoudoro and @gabknight for the feedbacks! :)

ShreyasFadnavis pushed a commit to ShreyasFadnavis/dipy that referenced this pull request Sep 20, 2018

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