A utility for interpolation and resampling

[navjotk]

Over time the work on resampling input data in time and the spatial interpolation of off-grid points seem to have a lot of overlap and this could be cleanly set out in a set of utility functions/classes. I was looking at some existing implementations that we should try and use wherever possible, with @ggorman a few weeks ago, and I am now putting these up in an issue to enable discussion to make this more concrete.

I understand that scipy.signal.firwin can be used to precompute the coefficients required to use the PrecomputedSparseFunction class in devito. This can also be used to resample incoming data in time. I also believe that scipy.signal.lfilter is (functionally) scipy's counterpart of PrecomputedSparseFunction - although only in a single dimension.

Scipy also provides full resampling implementations in the form of scipy.signal.resample and scipy.signal.resample_poly. IIRC, resample only works for integer multiples because it uses FFT. I'm not sure why we can't use resample_poly, though.

Scipy also provides direct functions for computing Kaiser-function [1] based interpolation like scipy.signal.kaiserord and scipy.signal.kaiser_beta.

I didn't find any ready implementations of Lanczos-filter based interpolation in scipy. However, I did find an implementation in ObsPy.

References:
[1] Hicks, Graham J. "Arbitrary source and receiver positioning in finite-difference schemes using Kaiser windowed sinc functions." Geophysics 67, no. 1 (2002): 156-165.

[FabioLuporini]

@ggorman @navjotk @mloubout is this something we're still interested in?

[mloubout]

I would say so yes, Could be a good short project.

[YanivHollander]

Trying to work my way in this issue. I'm new to Devito framework and code base. It seems to me that a proper place to generalize to more interpolator types, like the ones suggested by @navjotk, is in https://github.com/devitocodes/devito/blob/master/devito/operations/interpolators.py. More classes can be added there that implement other types of interpolation, on top of the currently implemented linear interpolation. Any thoughts?

[FabioLuporini]

yes, I agree @YanivHollander .

[YanivHollander]

Inject operation is unclear to me. (bi)linear interpolation from regular to sparse grid seems to be working fine.

Injection back from sparse grid to regular one is less convincing

Which algorithm is employed for injection under LinearInterpolator?

Python:
Issue595.py.zip

[mloubout]

It's linear, the injection is the adjoint of interpolation, so same bilenear coefficietns except it's a scatter instead of a gather.

[YanivHollander]

I need some guidance on the following:

1. The current implementation of SparseTimeFunction or PrecomputedSparseTimeFunction does not have a resample in time functionality. I could only find bi/tri-linear interpolation in space. There aren’t any examples that use resampling in time using these class ‘interpolate’ function.

2. In at least one place in the examples I found a resampling in time function that uses scipy B spline (example/seismic/source.py, line 182). Do we want to bring this functionality to either SparseTimeFunction or PrecomputedSparseTimeFunction classes?

3. Scipy resampling scipy.signal.resample and scipy.signal.resample_poly can be used for upsampling. However, the interpolation is to a regular grid. Those function cannot, to the best of my understanding, be used to interpolate to an irregular grid, unlike bi-linear interpolation, for example. Any thoughts on this point?

4. Have you considered precomputing coefficients for spline in space?

[mloubout]

1. No we didn't use it for that but would be a great addition.
2. I guess it could be moved in there. Would need to move TimeAxis or an equivalent somewhere too.
3. Yes this is why we moved to interpl as we rarely have nicely sampled time serie to use resample.
4. I am not sure, could you elaborate?

[navjotk]

So I think the rationale behind this issue was to combine the interpolation in space and time and move to a common place and @YanivHollander seems to be going there in the last comment.

[YanivHollander]

Hi,
Thanks for your response. I think I have enough for now to start designing my additions. I will share the design once completed for you to comment, before implementing. Do you prefer I would continue commenting under this issue, or move to Slack?

Out of all the options to interpolate in space, I didn't find any consideration of spline in space. That's what I meant, @mloubout, in point number 4.