Add spline_filter and spline_filter1d

[grlee77]

This PR adds the spline prefiltering functions that are correspond to:
scipy.ndimage.spline_filter
scipy.ndimage.spline_filter1d

The periodic modes "wrap" and "grid-wrap" are not currently supported, but all other modes are.

The one empirical nature of this implementation is the equation I came up with for estimating the overlap (depth) needed based on the interpolation order. In practice these are infinite impulse response (IIR) filters, but in practice the error terms decay pretty quickly with distance from an edge and the rate of decay is related to the magnitude of the filter poles. The current settings are enough for tests to pass at the tolerance of 1e-6 used in the tests. With the current default, the depth currently ends up being 11, 14, 19 and 22 for spline order of 2, 3, 4 and 5, respectively. I exposed a depth keyword-only argument in case users want to override this setting, but am not certain if that is desirable.

This PR also enables use of prefilter=True in the existing affine_transform function when mode='constant'. This is needed to get good results with order > 1. Although reflect and mirror are also supported by the spline_filter functions added here, they are not currently supported in affine_transform itself. I have not yet tried to support mode='nearest' or mode='grid-constant' for use within affine_transform because those two cases do not have an analytic solution implemented and involve the need for temporary internal padding in SciPy's own implementation, making it more difficult to adapt here.

[grlee77]

cc @m-albert regarding previous discussion related to this

[grlee77]

This is a simple performance benchmark on a somewhat large 3D volume (512 x 512 x 512)

from time import time

import dask.array as da
import numpy as np
import scipy.ndimage as ndi

from dask_image.ndinterp import spline_filter


# define test image
a = 512
ndim = 3
shape = (a,) * ndim
dtype = np.float64
im = np.random.standard_normal(shape).astype(dtype, copy=False)

# transform into dask array
chunksize = (128,) * ndim
dim = da.from_array(im, chunks=chunksize)

kwargs = dict(order=3, mode='constant', output=dtype)

tstart = time()
y_cpu = ndi.spline_filter(im, **kwargs)
dur = time() - tstart
print(f"Duration (scipy.ndimage) = {dur}")

tstart = time()
y = spline_filter(dim, **kwargs)
y = y.compute()
dur = time() - tstart
print(f"Duration (dask-image) = {dur}")

On my system with a 10-core Intel Skylake X CPU the result is:

for order=3

Duration (scipy.ndimage) = 7.603273630142212
Duration (dask-image) = 1.2285776138305664

or for order=5

Duration (scipy.ndimage) = 11.27594804763794
Duration (dask-image) = 1.9986212253570557

[grlee77]

I see there are quite a few test failures. I think a number of these tests cases will have to be skipped when SciPy < 1.6. I can fix that locally and then update here, later.

Currently all cases are passing in my local environment with SciPy 1.6.3 installed.

[m-albert]

Hi @grlee77, this is great!

The speedup on the CPU is considerable and of course this allows applying spline_filter on the GPU.

Really nice that you added prefilter functionality to affine_transform. Actually, the otherwise resulting blur has just been mentioned by @martinschorb. The implementation and tests for this look really good to me.

Currently all cases are passing in my local environment with SciPy 1.6.3 installed.

Locally all tests passing for me, too.

I see there are quite a few test failures. I think a number of these tests cases will have to be skipped when SciPy < 1.6. I can fix that locally and then update here, later.

Looking at the CI tests it seems there could be an error unrelated to the scipy version. In spline_filter1d, map_overlap seems unable to deal with a too large depth:

Call of map_overlap with:

x = dask.array<array, shape=(48, 48, 48), dtype=float64, chunksize=(33, 33, 33), chunktype=numpy.ndarray>
depth = {0: 0, 1: 0, 2: 22}
boundary = {0: 'none', 1: 'none', 2: 'none'}

resulting in:

E ValueError: The overlapping depth 22 is larger than your
E smallest chunk size 15. Rechunk your array
E with a larger chunk size or a chunk size that
E more evenly divides the shape of your array.

Probably the affected chunk is the last one along the array, which is smaller than the previous ones. Could this potentially be avoided by previous padding? Strange though that this doesn't occur locally...

[grlee77]

Probably the affected chunk is the last one along the array, which is smaller than the previous ones. Could this potentially be avoided by previous padding? Strange though that this doesn't occur locally...

I am guessing this is related to the Dask version. I have the recent 2021.4.1 release locally, but the CI job I looked at is using 2.8.1 (released Nov 2019).

[grlee77]

local testing indicates that the test failures don't occur for dask>=2020.12.0. Specifically, it is this PR that fixed it dask/dask#6708

[grlee77]

The one case that is expected to still fail on current dask is if the depth is larger than any of the overall array dimensions.

I could add a test for that case if desired (e.g. we cannot spline_filter a small image of shape (16, 16) with order=5).

edit: actually we should be able to support that case. We probably just need to make sure to set depth=0 on any axis where there is only a single chunk.

[m-albert]

I am guessing this is related to the Dask version. I have the recent 2021.4.1 release locally, but the CI job I looked at is using 2.8.1 (released Nov 2019).

I am guessing this is related to the Dask version. I have the recent 2021.4.1 release locally, but the CI job I looked at is using 2.8.1 (released Nov 2019).

Of course, that makes sense! Thanks for checking this.

[GenevieveBuckley]

The coveralls check shouldn't block this PR. I've just had a look at it & things there are mostly fine.

[GenevieveBuckley]

Thank you for the comment explaining the source of these magic numbers 😄

[GenevieveBuckley]

This looks good to me. I've made a few suggestions, and will leave this open for a bit for other reviewers to make comments too.