SciPy N-D interpolation functions

[jakirkham]

From @jakirkham on May 16, 2017 18:35

The list of functions below are what that this library intends to support with Dask Arrays. These all come from SciPy. In particular they come from scipy.ndimage.interpolation. The intent is to emulate their behavior with Dask Arrays as wrapping them is not likely to work. Will also include tests to verify these retain the behavior that the SciPy functions ordinarily have.

• affine_transform
• geometric_transform
• map_coordinates
• rotate
• shift
• spline_filter
• spline_filter1d
• zoom

Copied from original issue: dask-image/dask-ndinterp#3

[wtbarnes]

@jakirkham I'm very interested in implementations of affine_transform and shift. I'd be happy to put some effort into implementing both of these. Have there been any efforts already to implement either of these functions?

[GenevieveBuckley]

@jakirkham I'm very interested in implementations of affine_transform and shift. I'd be happy to put some effort into implementing both of these. Have there been any efforts already to implement either of these functions?

That's great! No, there aren't any current efforts implementing these so the playing field is wide open for you. I assume you're here at the scipy sprints? Come by the nd-image sprint room tomorrow and we can talk then :)

[wtbarnes]

Great! Unfortunately I am not at SciPy this year. I just watched the dask-image talk on YouTube though. I've recently been looking for a way to do affine transforms with Dask so this is great!

Do you have a suggestion of an existing function that might be good template for getting started? I'm assuming that it is mostly a matter of wrapping these ndimage functions with map_blocks appropriately?

[jni]

@wtbarnes it depends what you mean by "mostly". An affine transform means that your output pixels can come from anywhere in the input array, so it's not a straight mapping from input block to output block. Indeed, in the general case the output blocks will never be aligned with the input blocks.

[m-albert]

Hey @jakirkham @GenevieveBuckley @wtbarnes @jni,

I wrote a function that implements ndimage.affine_transform for dask arrays and takes care of what @jni mentions, that in general input and output blocks are not aligned.

Here's a gist containing the code with example usage and timing.

Would be cool to include this functionality in dask-image!

[jni]

Wow that is awesome, @m-albert! I am not a maintainer here but I'm pretty sure a PR would be very welcome!

[GenevieveBuckley]

It's late here, so I've scrolled very quickly through your gist @m-albert. A pull request would be very welcome!

[grlee77]

Thanks @m-albert, this is interesting! I had been hoping to look at this at some point for GPU applications where memory limits can sometimes be an issue, so it is great to see this here.

In regards to the amount to pad the blocks, it is a little bit tricky, but I will explain what I think is the solution below. The first consideration are the samples that fall within the footprint of the spline kernel. I see you currently extend by 2 samples at each boundary to attempt to account for that which should be enough for all but order 5 where the filter is width 6.

Specifically, the spline interpolation kernel has width order + 1 with starting and (inclusive) ending samples used to interpolate a coordinate, x, given by the following function adapted from the C source

def coord_range(x, order=5):
    """Inclusive coordinate range for spline filter"""
    if order % 1:
        start = math.floor(x) - order // 2
    else:
        start = math.floor(x + 0.5) - order // 2
    stop = start + order
    return (start, stop)

So even for order 5, interpolating a value at location 0.2 for example would involve weights applied to samples from coordinate locations [-2, -1, 0, 1, 2, 3].

However, there is another, possibly more important consideration in the case of order > 1. For those cases, there is also a spline prefiltering step that involves application of an infinite impulse response (IIR) filter. In short this means there will be boundary errors in the prefiltering that decay with increasing distance from the block boundaries and may be noticeable for 10 or more samples. Rather than use a larger padding, the optimal workaround for that should be to first perform prefiltering independently via ndimage.spline_filter and then apply affine_transform with the argument prefilter=False.

spline_filter should be pretty trivial to split up for dask and avoid artifacts because it is separable across axes (i.e. spline_filter just applies spline_filter1d along each axis in turn). For dask-image, when filtering along axis n we should set it to have only 1 chunk on that axis to avoid boundary artifacts, but multiple chunks can be used on any of the other axes. Once you have the prefiltered result, proceed as proposed here, but using prefilter=False.

[m-albert]

@jni @GenevieveBuckley nice, thanks!

@grlee77 Thanks a lot for raising these points and the very useful and detailed explanations!

Regarding the footprint of the spline kernels: Thanks to the code snippet you provided, proper handling of the footprint should be no problem to include in a dask wrapped version.

Regarding the prefiltering, my current understanding would be the following: There should definitely be a straightforward way of implementing the full ndimage.affine_transform functionality by including an iterative 1d prefiltering step with varying chunking as you propose. This would be useful as the entire function can be executed in a way that each step is performed in a chunkwise manner. This would probably be important for e.g. processing large arrays on the GPU. On the other hand, a main functionality people might be looking for in a dask implementation of ndimage.affine_transform consists in processing arrays that do not fit into RAM. However, if I'm not mistaken, the variable chunkings in the prefiltering step would load the entire input array into memory.

Of course, to prevent this, prefilter=False could easily be set. Another option would be to set this (or order=1 instead of scipy's order=3) as the default value for the dask-image implementation and clearly state it in the docstring (not sure what the policy of dask-image would be for changing a wrapped function's API).

In any case, I'll try to open a PR later today, which will provide a better space for this discussion.