FIX: Interpolation of B-Spline grid on the target EPI grid

[oesteban]

This PR fixes #289 and will address #218, although no solution (computationally tractable) has been found yet.

[oesteban]

@mgxd here's where your data break now. I'm 99% that this is the problem, the lazy (fast) evaluation of B-Spline weights assumes the alignment of the two coordinate systems (i.e., all direction cosines are collinear).

cc/ @effigies

[oesteban]

This line addresses the problem stated in #289 in just one shot.

That's why at "fit" time, when the estimated fieldmap is applied to generate reportlets (and validate everything), the interpolation works (and yields the same output of topup, saving some details). But when it reaches the correction of the EPI, the B-Spline grid is still a grid, but it's not aligned with the EPI grid.

[oesteban]

Cython improvement:

In [1]: import bspline

In [2]: from sdcflows.transform import _cubic_bspline
221111-11:56:02,394 nipype.utils WARNING:
	 A newer version (1.8.4) of nipy/nipype is available. You are using 1.8.4.dev0

In [3]: from numpy import frompyfunc

In [4]: bs = frompyfunc(bspline.bs_eval, 1, 1)

In [5]: %timeit bs(data)
699 ns ± 0.701 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)

In [6]: %timeit _cubic_bspline(data)
23.5 µs ± 102 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

[oesteban]

I have split this PR into smaller ones (#302, #301, #304, and #303). Since I had not managed to resolve the rotation issue, I'm going to open a fresh branch with the Cython code finally resolving this issue.

(I believe that, in the presence of this rotation, there's no way around calculating the b-spline for all combinations of voxels and control points in three dimensions, and hence, Cython)

[oesteban]

@effigies, same test, with scipy.signal.cubic

In [1]: import bspline

In [2]: from scipy.signal import cubic

In [3]: from numpy import frompyfunc

In [4]: bs = frompyfunc(bspline.bs_eval, 1, 1)

In [5]: import numpy as np

In [6]: data = np.array([1.4, 2.1, 0.7, -1, 3.0])

In [7]: %timeit bs(data)
699 ns ± 0.701 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)

In [6]: %timeit cubic(data)
22.2 µs ± 50 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

So basically the same as our _cubic_bspline implementation with np.piecewise.

[effigies]

Could be worth pushing the cython implementation upstream into scipy?

[oesteban]

that'd be interesting :)

I'm not sure I have the bandwidth:

• it will require more thorough testing than ours,
• they will probably want a more comprehensive implementation (e.g., allowing several bspline degrees, and/or derivatives)

so, happy to consider it if we find it necessary to go around our problems, but I don't know how much workload we will be able to take on.

[oesteban]

BTW:

In [1]: import numpy as np
   ...: from scipy.signal import cubic
   ...: from sdcflows.lib.bspline import bs_eval

In [2]: our_cubic = np.frompyfunc(bs_eval, 1, 1)

In [3]: data = np.array([1.4, 2.1, 0.7, -1, 3.0])

In [4]: np.allclose(cubic(data), our_cubic(data).astype("float32"))
Out[4]: True

(for some weird reason, our_cubic returns dtype object)

[effigies]

I believe cubic() is a special implementation from bspline(), so I don't think you would need to generalize it. And I would hope they have tests.

Anyway, agreed that it could be a deep rabbit hole and we should not block any of our process on it. I wonder if just posting your implementation as an issue on scipy could start the ball rolling and somebody else might eventually pick it up and do the integration work.

[oesteban]

Anyway, agreed that it could be a deep rabbit hole and we should not block any of our process on it. I wonder if just posting your implementation as an issue on scipy could start the ball rolling and somebody else might eventually pick it up and do the integration work.

This is something I can definitely do. I will try to first finish up the sdcflows problem and then think of how to make it more widely available.