[SciPy-Dev] Spline interpolation in ndimage

[Jaime Fernández del Río]

Hi all,

We have been discussing in issue #8465
<https://github.com/scipy/scipy/issues/8465> about the extension modes for
the interpolation module of ndimage. As some other issues linked in that
one show, this is an old problem, that has been recurringly discussed, and
there mostly is agreement on what the correct behavior should be. But as
all things ndimage, implementing any change is hard, in this case not only
because the code is complicated, but also because of the mathematical
subtleties of the interpolation method used.

In trying to come up with a way forward for this I think I can handle the
code complexity, but am having trouble being sure that I come up with a
sound math approach. I think I have more or less figured it out, but I
don't have a good academic background on signal processing, so was hoping
that someone who does (I'm thinking of you, Chuck!) can read my longish
description of things below and validate or correct it.

My main source of knowledge has been:

M. Unser, "Splines: a perfect fit for signal and image processing," in IEEE
Signal Processing Magazine, vol. 16, no. 6, pp. 22-38, Nov 1999.

Recommendations for bigger, better readings on the subject are also welcome!

-----

In what follows I'll assume the input image is 2-D and NxN for simplicity,
but I think everything generalizes easily

If I'm understanding it right, all of our interpolation functions use
B-splines for interpolation. So instead of having NxN discrete values on a
grid, we have NxN B-splines centered at the grid points, and we keep the
NxN coefficients multiplying each spline to reproduce the original image at
the grid points exactly. If the spline order is < 2, the spline
coefficients are the same as the image values. But if order >= 2 (and we
use 3 by default), these have to be calculated in a non-trivial fashion.
This can be done efficiently by applying a separable filter, which is
implemented in ndimage.spline_filter.

Because B-splines have compact support, when using splines of order n we
only need to consider the B-splines on an (n + 1)x(n + 1) grid neighborhood
of the point being interpolated. This is more or less straightforward,
until you move close to the edges and all of a sudden some of the points in
your grid neighborhood fall outside the original image grid. We have our
extend mode, which controls how this points outside should be mapped to
points inside. But here is where things start getting tricky...

When the spline coefficients are computed (i.e. when ndimage.spline_filter
is called), assumptions have to be made about boundary conditions. But
ndimage.spline_filter does not take a mode parameter to control this! So
what I think ndimage does is compute the spline coefficients assuming
"mirror symmetric" boundary conditions, i.e.:

a b c d c b | a b c d | c b a b c d

So if our interpolated point is within the image boundaries, but some of
the grid points in its (n + 1)x(n + 1) neighborhood fall outside the
boundary, the code uses a mirror symmetric extension mode to fill in those
values. This approach has a funny smell to it, but even if it's formally
incorrect I think it would only be marginally so, as the values are
probably going to be very close to the correct ones.

The problem comes when the point being interpolated is outside the
boundaries of the image. We cannot use mirror-symmetric spline coefficients
to extend if e.g. we have been asked to extend using wrap mode. So what
ndimage does is first map the value outside the boundaries to a value
within the boundaries, using the given extension mode, then interpolate it
as before, using mirror-symmetric coefficients if needed because its (n +
1)x(n + 1) neighborhood extends outside. Again, this smells funny, but it
is either correct or very close to correct.

This is mostly all good and well, except for the "wrap" and "reflect"
extension modes: in these cases the area within one pixel of the image
boundaries is different from anything inside the image, so we cannot use
that approach. So what ndimage does is basically make shit up and use
something similar, but not quite right. "reflect" is mostly correct, except
for within that pixel of the boundary, but "wrap" is a surprising and
confusing mess.

So how do we fix this? I see two ways forward:

   1. What seems the more correct approach would be to compute the spline
   coefficients taking into account the extension mode to be used, then use
   the same extension mode to fill in the neighborhood values when
   interpolating for a point outside the boundaries.
   1. First question is whether this is doable? I need to work out the
      math, but for "wrap" it looks like it should be, not 100% sure if also is
      for "reflect".
      2. Assuming it is it has the main advantage of doing things in a more
      general and understandable way once you have enough background knowledge.
      3. It does go a little bit against our API design: you can control
      whether the input is spline-filtered automatically with a parameter, the
      idea being that you may want to do the filtering yourself if you
are going
      to apply several different transformations to the same image. If the mode
      of the filtering has to be synced with the mode of the transformation,
      letting the user do it themselves is a recipe for disaster, because it's
      going to lead to very hard to track bugs.
      4. As elsewhere in ndimage, the current implementation does a lot of
      caching, which works because it always interpolates for a point
within the
      image boundaries. If we started interpolating for points outside the
      boundaries without first mapping to within there may be a performance hit
      which has to be evaluated.
   2. The other approach is,  for "wrap" and "reflect" modes, pad the input
   image with an extra pixel in each direction, then compute our current
   "mirror symmetric" spline coefficients, and leave things as they are right
   now, aside from some changes to the mapping of values to take the extra
   pixels into account.
      1. This looks like a nightmare of special cases everywhere and
      potential off-by-one errors while putting it together, but it
would just go
      along with the ugliness of the existing code.
      2. It's unclear what we would do if we are given an input with
      prefilter=False, so this also breaks the current API, probably
even more so.

Any thoughts or recommendations are very welcome!

Jaime

-- 
(\__/)
( O.o)
( > <) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes
de dominación mundial.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20180227/661d9fbe/attachment-0001.html>