ActiveAx analytic formulas

[davidrs06]

• Added analytical formula for CylinderGPD, Zeppelin and Ball models (removes dependency on Camino and accelerates the computation of the LUT by roughly 30s in the ActiveAx tutorial).
• Changed Zeppelin and Ball model to dipy, for consistency with StickZeppelinBall and FreeWater models
• N.B.: The output for a Zeppelin given by dipy doesn't match Camino's output, so recovered diameters don't match previous versions of AMICO

[davidrs06]

@daducci : do you remember if there's a reason for the scaling of the b-value for the high resolution scheme to be b_scale=1E6 for the CylinderZeppelinBall model and b_scale=1 for the StickZeppelinBall model ?

[daducci]

Yes, 1e6 was needed to adapt the units for camino, whereas if only tensors were useds (StockZeppelinBall) this was not necessary because all we used models from dipy.

[daducci]

Why the ouput of the zeppelins do not match? They are simple tensors, that's quite strange.

And what abut the cylinders? Did you validate that the signal is the same?

[daducci]

@davidrs06 , it would be great if you could speed up a bit more using cython. This way we might be able to store in files only one response function (in SH space), all the rotation matrices and then perform the rotations on the fly when fitting. This could solve an issue that @barakovic is facing these days in Cardiff. I remember you did it some time ago but we dropped it because it was about 2x slower; but if we use cython we could afford this operation on the fly!

[davidrs06]

Ok for units, for consistency with Camino. But then it doesn't really change anything as the scheme was saved in Stejskaltanner format (without the b-value). The thing is, for clarity, can the scaling be set to b_scale=1.0 now that we don't use Camino anymore, this way we don't have to divide the b-value by 1E6 when creating the gtab for dipy (and for consistency with the StickZeppelinBall and FreeWater models) ?

[davidrs06]

I double checked the output from Camino and CylinderGPD, they match when using np.allclose(). Not sure why Camino and dipy don't match for the Zeppelin, but the difference is small, around the 6th or 7th decimal for each direction, which is enough for np.allclose() to return False

[davidrs06]

Good idea to use cython for the rotation and resampling, I will have a look at it and see how the timing is affected.

[davidrs06]

Hi ! Could we create a separate PR for the rotation on the fly ? Would make more sense to keep the two features separated. Also because I'm messing up a lot with my code right now (the rotation on the fly means you don't have to compute the entire LUT, you have to pass the auxiliary matrices to the fit function, etc... and my code is getting quite ugly) and it would be better to have a separate PR where we can discuss the best way to adapt the behavior of AMICO depending on the status of the flag.

[daducci]

Yes, that is the best way to go!

[davidrs06]

Good, thanks ! Then, regarding the difference between Camino and Dipy for the Zeppelin and Ball, I think the difference might come from the computation of the b-value. Camino computes the b-value from the scheme file, and Dipy uses the b-value given by AMICO right ? If Camino and AMICO don't use the same gyromagnetic ratio to compute the b-value for example, you might get a difference in signal.

[daducci]

That would be strange, indeed it's worth investigating it. However if the bvalues were so different, then also the signals would look very different, no?

[davidrs06]

Yes, the generated signal is different. Camino uses GAMMA=2.6751525e8 while the Scheme class uses GAMMA=267.513e6. The b-value given to dipy's DTI method then generates a zeppelin with slightly different values (starting after the 6th or 7th decimal), but it's enough to give a different fitting in the tutorial for ActiveAx.

Should we change the value for GAMMA in the Scheme class to 2.6751525e8 ? This yields a Zeppelin with identical values as the one generated by Camino, but it would then change a little bit the results for all fittings that have been done with b-values computed with the current Gamma.