Extra NaNs imputed at grid points when interpolating

[Cody-G]

Is this expected behavior?

using Interpolations
img = rand(4,8)
img[:,1:2] = NaN
img[:,7:8] = NaN
@show img
itp = interpolate(img, BSpline(Linear()), OnGrid())
img2 = itp[indices(img)...]
#Notice that column 6 gets set to NaN
@show img2

ref https://github.com/JuliaImages/ImageTransformations.jl/issues/45

[sglyon]

I'd have to take a closer look, but my first guess is that the equations we are using don't specialize the "on the grid point evaluation" case and just implements the Linear B-spline rules.

These rules always use function values at two points for evaluation. This means that the function value at 7 is used whenever interpolating on the interval [6, 8). In this case the NaN at 7 is "poisonous" and infects the interpolated value in that entire interval...

[Cody-G]

I was thinking it may be something like that, but then why doesn't column 3 get set to NaNs as well?

[sglyon]

My guess (again I'd double check this) is that when evaluating at 3 it uses nodes 3 and 4.

If my hypothesis is correct, then the equation for on grid evaluation at grid point i must be (in one dimension for simplicity)

1 * img[i] + 0 * img[i+1]

so for i = 3 you use img[3] and img[4] and neither of those are NaN so you are good

[Cody-G]

I see, that also makes sense. Do you think it makes sense to special-case the grid points? I discovered this behavior when using interpolations with coordinate transformations via ImageTransformations.jl, but I think it would also be unintuitive in other domains. If you're correct I don't see a way around it other than special-casing the grid points. Maybe this is the relevant line for linear BSplines?

https://github.com/JuliaMath/Interpolations.jl/blob/master/src/b-splines/linear.jl#L98

[tomasaschan]

@sglyon's analysis is correct here - the NaNs in the sixth column are there because the seventh is involved in the calculation of their values.

I think special-casing grid points might work for this specific case, but it would make the code much more complex (and probably quite a lot slower, because we'd need a branch in a hot spot) for relatively small gains. Instead, we can utilize features already in the library to work around this. Let me show you how:

This is the data we're working with, as per your example.

using Interpolations
img = rand(4,8)
img[:,1:2] = NaN
img[:,7:8] = NaN
img

4×8 Array{Float64,2}:
 NaN  NaN  0.909054  0.500433   0.803419  0.865767  NaN  NaN
 NaN  NaN  0.444361  0.201766   0.672682  0.136034  NaN  NaN
 NaN  NaN  0.615366  0.741219   0.664258  0.290056  NaN  NaN
 NaN  NaN  0.709181  0.0974666  0.724229  0.735037  NaN  NaN

If we interpolate it naively, we get a bunch of NaN's in the sixth column, because columns 6 and 7 are used to calculate the value there.

itp1 = interpolate(img, BSpline(Linear()), OnGrid())
itp1[indices(img)...]

4×8 Array{Float64,2}:
 NaN  NaN  0.909054  0.500433   0.803419  NaN  NaN  NaN
 NaN  NaN  0.444361  0.201766   0.672682  NaN  NaN  NaN
 NaN  NaN  0.615366  0.741219   0.664258  NaN  NaN  NaN
 NaN  NaN  0.709181  0.0974666  0.724229  NaN  NaN  NaN

The simplest fix is probably to interpolate only non-NaN values, and then use a coordinate transformation that lets you use the original coordinate system. It can be done manually, but by utilizing the features of Interpolations you can actually get the interpolation object to do the work for you.

img_without_nans = img[:,3:6]
itp = interpolate(img_without_nans, BSpline(Linear()), OnGrid())
itp[:,:]

4×4 Array{Float64,2}:
 0.909054  0.500433   0.803419  0.865767
 0.444361  0.201766   0.672682  0.136034
 0.615366  0.741219   0.664258  0.290056
 0.709181  0.0974666  0.724229  0.735037

Next, use a scaling object to re-scale into the original coordinate system, with the second index in 1:8 instead of -1:6:

sitp = scale(itp, 1:4, 3:1:6)

Finally, extrapolate to re-add the NaNs at the edges:

eitp = extrapolate(sitp, NaN)

For some reason, evaluating scaled interpolations doesn't work with sitp[indices(img)...] - this is a bug, and should probably be filed as an issue if it isn't already. (Not 100% sure if the problem is with Interpolations.jl or base, but more likely here...)

But we can still inspect the output, albeit in a more roundabout way:

out = Array(Float64, 4,8)
for x in 1:4, y in 1:8
    out[x,y] = eitp[x,y]
end
out

4×8 Array{Float64,2}:
 NaN  NaN  0.909054  0.500433   0.803419  0.865767  NaN  NaN
 NaN  NaN  0.444361  0.201766   0.672682  0.136034  NaN  NaN
 NaN  NaN  0.615366  0.741219   0.664258  0.290056  NaN  NaN
 NaN  NaN  0.709181  0.0974666  0.724229  0.735037  NaN  NaN

[Cody-G]

Thanks for the detailed response! I see how the scaling trick can help when NaNs are along the edges of the array, but what about when they're in the middle? For example, interpolating this array increases the number of NaNs at grid points from 1 to 4 (1+ 1 per axis + 1 on the diagonal), all surrounded by finite points.

using Interpolations
img = rand(5,5)
img[3,3] = NaN
itp = interpolate(img, BSpline(Linear()), OnGrid())
img2 = itp[indices(img)...]
@show img
@show img2

I share your concern about performance. I wonder if we can eliminate the branching penalty with an ifelse in the line that I pointed out (at least for linear interpolation), but before I look into it I wanted to see if people think it's a good idea, performance aside. If it's not then maybe at least we could document somewhere the asymmetry of the current interpolation scheme and its consequences for NaN elements?

[sglyon]

Sorry for the fly by comment — I’m on my phone..

Could we special case on grid evaluations by defining a getindex method on integer subtypes?

[tomasaschan]

@sglyon I don't think that would help, since the type will be Float64 even if the value is an integer...

Also, we might be able to get it to work for linear interpolation, but the code you're pointing to is today completely generic in order. Special-casing for linear interpolation would make an already quite complex part of the library considerably more complex, with IMO very little gain.

More to the point, I don't think it should be the responsibility of this package to clean the input data - it's much better if you can somehow eliminate the NaNs before interpolating.

@Cody-G, what is your actual use case here? If the grid points are the important points, why use an interpolation at all?

[tomasaschan]

Sorry - I realize now that I was probably adressing @Cody-G already from the second paragraph :)

[sglyon]

@tlycken wouldn't it help to implement getindex(::Interpolations.BSplineInterpolation, is::Integer...) because then we know that we are evaluating on a grid point and can just do a lookup into the underlying array and bypass the B-spline formula entirely? That's fine that the return type of getindex is Float64 -- but I may have misunderstood what you were talking about

[Cody-G]

Thanks again for your replies. The specific domain where this came up for me is in image transformations. The example discussed in https://github.com/JuliaImages/ImageTransformations.jl/issues/45 is a simple one involving a simple translation of an image with no NaNs in interior pixels. That one is definitely solvable with the scale method that you suggest.

But in more complex cases involving rotations and nonlinear deformations (where it seems scale won't help), we usually use NaN values to denote areas of the image where we have no information. If we impute a finite value instead then it can bias some calculations inappropriately, so one can argue that NaNs are actually the cleaner option in some cases. But the dilation of NaN regions that occurs with each interpolation is not ideal. This can be solved by special-casing whole-pixel translations and deformations in the downstream code (as in ImageTransformations where I first created the issue). It's likely that many image transformations will involve sampling some pixels on grid points and others between them, so it may be a pain to implement and I was just hopeful that it could be done here with less difficulty.

@Evizero you may have something to add to this conversation

[tomasaschan]

@sglyon I think it would be really confusing if itp[1:3] and itp[1.:1.:3.] would yield NaN for the data point at 2 in one case but not the others. And there's no way for dispatch to know choose a different method for 1.0 vs 1.1, is there?

@Cody-G let me get back to you when I'm at a computer.

[sglyon]

Ahh I see. Point taken. Thanks for explaining.

I also agree that Interpolations shouldn’t be responsible for data cleaning.

[tomasaschan]

@Cody-G The problem with what you're requesting is that it is fundamentally something different than what Interpolations is meant to do.

Take the following very contrived example, with NaNs in the middle of the data set:

1  2  3   4
5  6  NaN 8
9  10 11  12
13 14 15  16

If you create an interpolation object from this data, and then evaluate again on the grid points, you'll get this result:

1.0   NaN    NaN    NaN  
5.0   NaN    NaN    NaN  
9.0   10.0   11.0   12.0
13.0  14.0   15.0   16.0

However, that's quite misleading. If you would interpolate also on half-points, you'd see that it's not only in the "upward" direction from where 7 would have been you get NaN, but also in the "downward" direction. In fact, for any row-index r and col-index c such that 1 < r < 3 and 2 < r < 4, you'll get NaN (as well as on some of those grid points where the < would be <=, but I'm leaving them out for now). Those indices form a square that's infinitesimally smaller than 3x3 around the original NaN value. The fact that the output above shows the numbers 10 through 12 instead of NaN could just as well be seen as a bug than the other way around.

[Cody-G]

If you would interpolate also on half-points, you'd see that it's not only in the "upward" direction from where 7 would have been you get NaN, but also in the "downward" direction.

I understand that, and I think it's perfectly reasonable that an interpolation at half-point between a NaN and a finite number should yield a NaN. What's less intuitive is that sampling at a finite grid point can also bring a NaN into the computation, and asymmetrically so.

The fact that the output above shows the numbers 10 through 12 instead of NaN could just as well be seen as a bug than the other way around.

I also agree with that statement. This is what motivated me to suggest the more conservative option of documenting the asymmetry if we don't want to change the code.

[tomasaschan]

The thing is, Interpolations.jl was never designed to handle NaN systematically - the asymmetry is just a byproduct of an implementation detail. I think the "most correct" thing to do here would actually be to throw an error if any element in the input is NaN, but since that's not a very useful behavior I'd rather keep the current one.

But as the asymmetry is just that - an implementation detail - it'd be contraproductive to document it (because documented features are features that users will, sooner or later, come to depend on). What if we change something about the implementation that flips the direction of the asymmetry, or even manages to remove it entirely?

If anything, we could document that data with NaN values will propagate NaNs in all calculations, so the output should be expected to contain some NaNs as well - garbage in, garbage out, as they say. However, since this is the case of NaNs everywhere, I'm not sure it makes sense to document it here (should e.g. ImageTransformations.jl document it too? Should everyone that accepts floating number inputs?). But feel free to file a PR for that if you wish :)