Forecast #1322

maurozucchelli · 2017-08-23T14:02:21Z

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 · 2017-08-23T15:59:55Z

Codecov Report

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

@@            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 · 2017-08-23T16:00:08Z

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

arokem · 2017-08-23T17:36:52Z

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.

skoudoro · 2017-08-23T15:32:12Z

dipy/reconst/forecast.py

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


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).

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

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)

skoudoro · 2017-08-23T15:35:12Z

dipy/reconst/forecast.py

+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


unused import. I think you can remove it

skoudoro · 2017-08-23T15:36:03Z

dipy/reconst/forecast.py

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


unused import. I think you can remove it

skoudoro · 2017-08-23T15:44:49Z

dipy/reconst/forecast.py

+        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:


PEP8: missing space around operator <=

skoudoro · 2017-08-23T15:51:02Z

dipy/reconst/forecast.py

+
+        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)


PEP8: missing space after block comment

skoudoro · 2017-08-23T17:35:06Z

doc/examples/reconst_forecast.py

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


unused import

skoudoro · 2017-08-23T17:35:11Z

doc/examples/reconst_forecast.py

+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


unused import

skoudoro · 2017-08-23T17:36:47Z

doc/examples/reconst_forecast.py

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


missing 1 blank line

skoudoro · 2017-08-23T17:38:18Z

doc/examples/reconst_forecast.py

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


PEP8: missing space after ,

skoudoro · 2017-08-23T17:40:20Z

doc/examples/reconst_forecast.py

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


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

skoudoro · 2017-08-24T20:13:48Z

dipy/reconst/forecast.py

+        >>> 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]))


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

skoudoro · 2017-08-24T20:14:05Z

dipy/reconst/forecast.py

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


same as above

coveralls · 2017-08-25T09:09:47Z

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

coveralls · 2017-08-25T09:09:47Z

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

coveralls · 2017-08-25T11:05:33Z

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

coveralls · 2017-08-25T11:05:33Z

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

maurozucchelli · 2017-08-25T12:04:19Z

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 · 2017-08-25T14:31:55Z

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 · 2017-09-20T08:33:23Z

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

skoudoro · 2017-09-20T13:18:12Z

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 · 2017-09-20T17:54:22Z

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?

skoudoro · 2017-09-20T18:05:08Z

doc/examples/reconst_forecast.py

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


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

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?

maurozucchelli · 2017-09-25T12:40:26Z

@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 · 2017-09-25T15:41:56Z

Might there be a 90 degree rotation somewhere relative to the other figure, and we're looking at the slice from the side?

…

On Mon, Sep 25, 2017 at 6:10 AM, maurozucchelli ***@***.***> wrote: ***@***.**** commented on this pull request. ------------------------------ In doc/examples/reconst_forecast.py <#1322 (comment)>: > @@ -0,0 +1,145 @@ +""" +============================================================== +Crossing invariant fiber response function with FORECAST model +============================================================== + +We show how to obtain a voxel specific response function in the form of +axially symmetric tensor and the fODF using the FORECAST model from +[Anderson2005]_ , [Kaden2016]_ and [Zucchelli2017]_. + +First import the necessary modules: +""" + +import matplotlib.pyplot as plt +from dipy.reconst.forecast import ForecastModel +from dipy.viz import fvtk 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? [image: fodfs] <https://user-images.githubusercontent.com/5042250/30809880-913683e4-a203-11e7-917a-000dacdc7980.png> — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#1322 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAHPNvvdNSgw_n-lZxZIGEoZh1-LlGlMks5sl6YsgaJpZM4PAD-t> .

Garyfallidis · 2017-09-26T21:37:22Z

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 · 2017-09-27T15:39:40Z

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 ?

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 · 2017-10-06T12:37:34Z

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

maurozucchelli · 2017-10-06T14:04:00Z

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 · 2017-10-06T14:09:10Z

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 · 2017-10-17T09:03:24Z

Hi I rebased! Now everything seems to work fine

Garyfallidis · 2017-12-01T22:12:12Z

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!! :)

Garyfallidis · 2017-12-01T19:12:32Z

doc/examples/reconst_forecast.py

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


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?

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.

Garyfallidis · 2017-12-01T19:13:11Z

doc/examples/reconst_forecast.py

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


You may want to add quotes around sh_order.

Garyfallidis · 2017-12-01T21:13:04Z

dipy/reconst/forecast.py

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


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

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!

Garyfallidis · 2017-12-01T21:23:19Z

dipy/reconst/forecast.py

+    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]_. 


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

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

Garyfallidis · 2017-12-01T21:45:44Z

dipy/reconst/forecast.py

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


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

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

Garyfallidis · 2017-12-01T22:00:43Z

dipy/reconst/forecast.py

+    return v
+
+
+def Psi_l(l,b):


lower case in function names

Garyfallidis · 2017-12-01T22:01:26Z

dipy/reconst/forecast.py

+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)


okay... spaces around operators.

Garyfallidis · 2017-12-01T22:02:38Z

dipy/reconst/forecast.py

+    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)


Garyfallidis · 2017-12-01T22:03:04Z

dipy/reconst/forecast.py

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


Garyfallidis · 2017-12-01T22:07:51Z

doc/examples/reconst_forecast.py

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


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 is good work so let's make this as beautiful as possible.

Garyfallidis · 2018-01-15T17:46:45Z

doc/examples/reconst_forecast.py

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


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

Garyfallidis · 2018-01-15T18:02:53Z

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 · 2018-01-15T18:03:21Z

Congrats @maurozucchelli!

maurozucchelli · 2018-01-16T10:46:50Z

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

Forecast

skoudoro reviewed Aug 23, 2017

View reviewed changes

skoudoro reviewed Aug 24, 2017

View reviewed changes

arokem mentioned this pull request Aug 25, 2017

Warnings raised in csdeconv for upcoming numpy 1.14 #1323

Closed

skoudoro reviewed Sep 20, 2017

View reviewed changes

maurozucchelli force-pushed the forecast branch from a1ec448 to ee5d7bf Compare October 16, 2017 11:34

maurozucchelli added 2 commits December 4, 2017 14:53

forecast pull request

fe7f28e

final version of test and example

70b1a01

maurozucchelli added 11 commits December 4, 2017 14:53

pep8

dbf2365

change example data

d1968b7

add predict and test_predict

c903979

fixed pep8 comments

7ce5ce2

add forecast example in the index

fbbb534

fix parent class and example indentation

f03c1a7

relax mse test because travis was failing in one version

2a26d2a

substitute my csd with dipy csdeconv

1c5b1c2

add general formula for Psi functions

a9231b3

change odf example

04b37ec

Merge remote-tracking branch 'nipy-dipy/master' into forecast

083b667

maurozucchelli force-pushed the forecast branch from ee5d7bf to 083b667 Compare December 5, 2017 15:03

Garyfallidis requested changes Dec 8, 2017

View reviewed changes

CSD is now the default dec_alg plus pep8 fixes

7641437

Garyfallidis approved these changes Jan 15, 2018

View reviewed changes

Garyfallidis reviewed Jan 15, 2018

View reviewed changes

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

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

Merge pull request dipy#1322 from maurozucchelli/forecast

6ea7e2d

Forecast


		sh_order is the spherical harmonics order used for the fODF.

		optimizer is the algorithm used for the FORECAST basis fitting, in this case

Forecast #1322

Forecast #1322

Conversation

maurozucchelli commented Aug 23, 2017

codecov-io commented Aug 23, 2017 • edited Loading

Codecov Report

coveralls commented Aug 23, 2017

arokem commented Aug 23, 2017

skoudoro left a comment • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

skoudoro Aug 23, 2017 • edited Loading

Choose a reason for hiding this comment

skoudoro Aug 24, 2017 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

coveralls commented Aug 25, 2017

coveralls commented Aug 25, 2017

coveralls commented Aug 25, 2017

coveralls commented Aug 25, 2017

maurozucchelli commented Aug 25, 2017

arokem commented Aug 25, 2017

maurozucchelli commented Sep 20, 2017

skoudoro commented Sep 20, 2017

arokem commented Sep 20, 2017

skoudoro Sep 20, 2017 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

maurozucchelli commented Sep 25, 2017

arokem commented Sep 25, 2017 via email

Garyfallidis commented Sep 26, 2017

skoudoro commented Sep 27, 2017

maurozucchelli commented Oct 6, 2017

maurozucchelli commented Oct 6, 2017

skoudoro commented Oct 6, 2017

maurozucchelli commented Oct 17, 2017

Garyfallidis commented Dec 1, 2017 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Garyfallidis commented Jan 15, 2018

Garyfallidis commented Jan 15, 2018

maurozucchelli commented Jan 16, 2018

codecov-io commented Aug 23, 2017 •

edited

Loading

skoudoro left a comment •

edited

Loading

skoudoro Aug 23, 2017 •

edited

Loading

skoudoro Aug 24, 2017 •

edited

Loading

skoudoro Sep 20, 2017 •

edited

Loading

Garyfallidis commented Dec 1, 2017 •

edited

Loading