Add support functions for spreading and interpolating with scaling in cufinufft

[chaithyagr]

This resolves #306
Also I have a bunch of fixes which did not make through for the earlier refactoring PR.

[chaithyagr]

This support is only added for cufinufft at the moment, adding to finufft soon.

[chaithyagr]

I have rebased all my changes to the latest version. I still will need the cufinufft_spread and interp functions for exposing them. @blackwer can you please comment on this and if the use of cuda_array_interface essentially makes this PR null?
Is it possible to call spread and interpolate functions only from python side?

[blackwer]

I have rebased all my changes to the latest version. I still will need the cufinufft_spread and interp functions for exposing them. @blackwer can you please comment on this and if the use of cuda_array_interface essentially makes this PR null? Is it possible to call spread and interpolate functions only from python side?

I'm not sure how the _cuda_array_interface_ stuff plays into this. We could certainly expose the spread/interp functionality as you've done here, and that functionality shouldn't get too much in the way of it. It's more an issue of if we want to publicly expose these functions, which might be annoying to maintain a stable API for. If we maintain a stable public API, there's nothing tying it to any particular language, since we expose the API via C. We're meeting this week to discuss some things and I'll bring this up. It might be tagged as unstable/unsupported, but I'm not opposed to exposing it.

[chaithyagr]

If it does come up in the meeting do mention that spread and interpolate functions are very important for the below 2 important reasons:

1. it is crucially needed for estimating good density compensators which is repeatedly used by folks who use NUFFT to reduce their lipschitz constant (to better shape the problem by weighting out the Fourier terms) which can significantly help in speeding up many reconstruction algorithms
2. spread and interp functions are very good approximate candidates for gridding and interpolation of data in general.

[turbotage]

I second @chaithyagr's message about the need for spreading/interp interface exposure. Also for me, the interest is in density compensation MRI. Perhaps it also nicely ties into any wishes for NUFFT inverse documentation or examples, #264, maybe a Pipe Menon example or similar. As a start simply exposing the spreading/interp would be greatly appreciated!

Great work!

[paquiteau]

Hello there,

more pragmatically, the exposed API can consist of a new argument/ flags in the (cu)finufft plan struct (e.g., spreadinterp_only) that tells that this plan (be it type 1 or type2) will only perform the spreading/interpolating; No FFT or deapodization. I had try such an approach a while back (mimicking what is being done in gpunufft in fact), but @chaithyagr was faster with his PR :D

[ahbarnett]

I'm confused where this code came from (surely it already exists in cufinufft?)

[chaithyagr]

It did indeed, it was however removed during cleanup as a part of 27c3621

This was perhaps done as it was used to test internal modules, but then it was this very functions which I chose to expose (with a wrapper)

[ahbarnett]

We have no idea what this code does without a mathematical description, and documentation of inputs/outputs.

[ahbarnett]

What is n here? This routine looks like a crude quadrature approximation to the integral of the spread kernel. But it is not accurate. A Gauss-Legendre rule would be needed.

[ahbarnett]

rank would need to be documented. This all suggests that this code should reside outside of (cu)FINUFFT, since it's not part of the NUFFT task.

[ahbarnett]

Dear MRI users,

We are glad you find (cu)FINUFFT useful, and want that to remain the case! But we (or at least I) am confused what additional features are wanted from our library going forward. Hence I moved #306 to a Discussion.
If it is exposing a C/C++ interface to spread and interp, in 1,2,3 dims, that is maybe doable, but is going to remain in the unsupported/undocumented category for a long time... and it adds confusion. The main reason is: our library performs a very well-specified task, namely a NUFFT, which has a clear mathematical definition:
https://finufft.readthedocs.io/en/latest/math.html
(Indeed, without this definition, the library would be unusable, and unstable.) How we choose to do spreading/interpolation is internal to the library, and could change at any time (upsampling factor, window function, nspread width, etc, which we choose carefully according to the user-requested tolerance) while still performing the fixed required task and passing its tests.

Thus what you're requesting (spread/interp and deconvolution) is not part of the interface to a NUFFT library, is unstable, and should remain under the hood. (Eg, one would not expect access to the internals of an FFT library.)
In fact, it's unclear precisely what the mathematical definition is for what you seek (what spread/interp function do you seek, how many grid-points wide, and why would it have the same optimality properties as what we happen to need to perform the NUFFT well?)
I tried to read the Pipe-Menon 1999 MRM paper but found it too domain-specific to understand. "Scale factors" are certainly too MRI-specific for us to even understand how to support (or what they are, what unit tests for them would be, etc).
[Another practical way to look at it: If we "switch off" the FFT and deconvolve steps inside (cu)FINUFFT, then one is left with a upsampfac^d.N sized array - and how should we pass that to the user? This is a whole new interface.]

Digging into the MRI application, here's why I'm confused re the requests in this PR. In standard iterative recon papers (eg Fessler and Noll, ISMRM 2007, or more recent papers of Lustig et al.), the imaging forward model is a 2D type 2 NUFFT, and usually given the matrix symbol A. Thus one wants to minimize ||Ax-y|| in some norm plus some regularization term. In these works, the algorithms all involve applications of A or A^H, ie type 2 and type 1 NUFFTs. Those are what you should use FINUFFT for. (I should check that those are in fact what you are using it for - right?) In contrast, density compensation weights appear to affect convergence rates, and are computed in some distinct offline step (which could also happen to involve NUFFTs, eg sinc^2 weights of Greengard, Inati, et al).

BTW, I should add: if you take the output of the type 1 (ie, the regular array) and send that to a sized-N FFT (which is very cheap), you get an interpolator from scattered data that makes no reference to window functions or density compensation. Maybe that is what you want?

Sorry about the length of this, but we have to all communicate about it successfully to know how to proceed. In conclusion: if you feel strongly about wanting these features that fall outside of the NUFFT interface, if would be great if you could please give us a literature reference that explains the precise mathematical formula desired, in a similar way to the A matrix in the above iterative MRI papers.

[ahbarnett]

I second @chaithyagr's message about the need for spreading/interp interface exposure. Also for me, the interest is in density compensation MRI. Perhaps it also nicely ties into any wishes for NUFFT inverse documentation or examples, #264, maybe a Pipe Menon example or similar. As a start simply exposing the spreading/interp would be greatly appreciated!

Great work!

Please see the 2.2.0 docs which now have a simple 1D inverse type 2 tutorial:
https://finufft.readthedocs.io/en/latest/tutorial/inv1d2.html
2D inversion to follow in next few days I hope. Best, Alex

[ahbarnett]

Apologies to the users who want interfaces from spreadinterp.cpp exposed directly - see discussion for why we're not going to be able to do it. Hwoever, it will be easy for you to do in your forks. PS spreadinterp just got 10-80% faster due to XSIMD vecotrization... Best, Alex