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.
On Tue, Feb 27, 2018 at 4:59 AM, Jaime Fernández del Río < jaime.frio@gmail.com> wrote:
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.
The problem with the factorization is that it assumes infinite data points in all dimensions, i.e, no explicit boundary conditions. With finite data there are edge effects when the data is deconvolved to get the spline coefficients. The way that ndimage deals with that is to extend the data using reflection and start far enough away that the edge effects have died away by the time the "real" data has been reached. How far away that is, is heuristic. I think it should not be too difficult to extend the data in the other ways, but note that since uniform splines are being used, the relevant coefficients lie outside the boundary and need to be picked up from the correct spots in the interior using the symmetries. IIRC, the b-splines are always centered at zero so that they are symmetrical. For odd order splines that will be at the center of a pixel (pixel points), for even order splines at an edge, pixel centers at half integer points. The data can always be considered to be at pixel centers, but the splines are displaced. I don't remember if ndimage treats that correctly.
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!
Chuck
On Thu, Mar 1, 2018 at 9:15 AM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Tue, Feb 27, 2018 at 4:59 AM, Jaime Fernández del Río < jaime.frio@gmail.com> wrote:
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.
The problem with the factorization is that it assumes infinite data points in all dimensions, i.e, no explicit boundary conditions. With finite data there are edge effects when the data is deconvolved to get the spline coefficients. The way that ndimage deals with that is to extend the data using reflection and start far enough away that the edge effects have died away by the time the "real" data has been reached. How far away that is, is heuristic. I think it should not be too difficult to extend the data in the other ways, but note that since uniform splines are being used, the relevant coefficients lie outside the boundary and need to be picked up from the correct spots in the interior using the symmetries.
IIRC, the b-splines are always centered at zero so that they are symmetrical. For odd order splines that will be at the center of a pixel (pixel points), for even order splines at an edge, pixel centers at half integer points. The data can always be considered to be at pixel centers, but the splines are displaced. I don't remember if ndimage treats that correctly.
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!
To explicate a bit more as to why the b-splines are centered, consider quadratic and cubic splines when the data points are considered to occur at the pixel centers. *Quadratic (even order)* 1. If the data points are taken to be halfway between knot points, we need to deconvolve the sequence `array([1, 6, 1])/8`. The Fourier transform of that has no zeros, indeed, is rather smooth, and two extra coefficients, one at each end, are required to interpolate out to the pixel edges. It is all nicely symmetrical. 2. If the data points are taken to correspond to the knot points, we need to deconvolve the sequence `array([1, 1])/2`. The Fourier transform of that sequence has a zero at the Nyquist, not good, and one extra coefficient is needed at some arbitrary end in order to get interpolation out to the pixel edges. *Cubic (odd order)* 1. Taking the data points to correspond to the knot points, we need to deconvolve the sequence `array([1, 4, 1])/6`, whose Fourier transform is nicely behaved. However *four* extra extra coefficients, two at each end, are required to interpolate out to the pixel edges. The "extra" coefficients are needed in order to cover the outer half of the edge pixels, otherwise we are extrapolating. I haven't looked, but I wonder if ndimage gets that right? In both cases, the corner pixels may be a bit of a problem. I haven't thought through that bit. Chuck
On Tue, Feb 27, 2018 at 4:59 AM, Jaime Fernández del Río < jaime.frio@gmail.com> wrote:
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
For "a, b, c", I would write that as "..., a, b, c, b, a, b, c, ..." extended to infinity in both directions because we are applying IIR filters. Fortunately, the IIR filters fall off rapidly, so one need not take too many extended points to start them up. Counter intuitively, I think the interpolation "mode" only applies to the spline coefficients, not to the original image, resulting in unexpected behavior at the edges. If we want to have "mode" apply to the image rather than the filtered results, we can do that, but the result of the filtering would best be a named tuple containing the mode, the spline order, and, in the case of a constant extension, extra edge coefficients. The coefficients for the "wrap" and "constant" modes can then be obtained by appropriately extending the data for starting up the filters. For the "wrap" and "reflect" modes, the coefficents have the same symmetries, so that can be used to get the needed coefficients outside the boundaries. The spline_filter using IIR is just a clever way to solve A*coef = data. For the cubic case and three data points, A depends on the mode as follows: |4 2 0| |1 4 1| x 1/6 (reflect about center of edge pixels -- we do this) |0 2 4| |5 1 0| |1 4 1| x 1/6 (reflect about edge -- we don't do this) |0 1 5| |4 1 1| |1 4 1| x 1/6 (wrap -- we don't do this) |1 1 4| Explicitly inverting such matrices also provides a good test to check that things are done right if we go this way. As you have pointed out, the proper coefficients will vary depending on the mode. The spline coefficients have the same symmetry as the image pixels. Unfortunately, `map_coordinates` reflects about the edges, which doesn't match any of the coefficients that we compute.
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".
Yes. The constant mode is actually the tricky one. I'd be tempted to make it orthogonal to "reflect", "wrap", and "nearest", that way we don't need to compute extra coefficients. Note that if the image isn't prefiltered, we are running a smoothing filter over it.
1. Assuming it is it has the main advantage of doing things in a more general and understandable way once you have enough background knowledge. 2. 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. 3. 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. 1. 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.
Are "mirror" and "reflect" actually different? The function documentation only mentions "reflect".
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!
The docstrings need fixing, they are almost useless. Chuck
participants (2)
-
Charles R Harris -
Jaime Fernández del Río