Map coordinate fix

[jjhelmus]

This pull request address Trac Tickets #796, #1378 and #1520 which all stem from the map_coordinate function in ni_interpolation.c

As currently written, this pull causes two tests in test_ndimage to fail (boundary modes and boundary modes 2) which is the reason for the long explaination. In these two tests the 'wrap' mode need further clarification. Take a simplified version of the boundary mode test:

def shift(x):  return (x[0] + 0.5,)
input = numpy.array([2,4,6,8.])
output = ndimage.geometric_transform(input, shift, mode='wrap', output_shape=(5,), order=1)

This generates an index mapping to take from the input for the output of [0.5, 1.5, 2.5, 3.5, 4.5]. With "wrap" mode determining what elements these correspond to is difficult.

The first element is well defined, halfway between the first and second elements in input, 2 and 4, so 3.
The second and third elements are also well defined, 5 and 7.

The fourth element has an index of 3.5, but the input only has indices to 3. Currently in scipy this is "wrapped" an index of 0.5, with a value of 3, but this seems incorrect to me. The fourth element in the input has an index of 3, if the data is wrapped (appended before and after the original data as needed) then the first element would have an index of 4. In this sense an index of 3.5 is halfway between the fourth (last) element and the first element and might be represented as -0.5. The problem come when we try to find a value for this element, we need to extrapolate a half point past either the last element or the first element.

The fifth element in the output has an input index of 4.5, which is currently "wrapped" to an index of 1.5 and takes the value of 5. In my mind, if an index of 4 wraps to the first element of the input, then the 4.5 should be between the first and second element of the input, so 3.

So the question becomes should output be:
1). As scipy has it currently so that 3.5 -> 0.5 and 4.5 -> 1.5. output = [3, 5, 7, 3, 5]
2). As in this pull so that 3.5 -> -0.5, 4.5 -> 0.5. output = [3, 5, 7, ?, 3] and then also what should ? be.

I'm in support of 2, but would be happy if someone can convince me that 1 is correct.

As an addition exame, please see what happens to the third element vs. all other elements of output when shift is defined as:

def shift(x): return (x[0] + 0.9,)
def shift(x): return (x[0] + 1.0,)
def shift(x): return (x[0] + 1.1,)

[rgommers]

Thanks a lot for working on these ndimage issues!

[rgommers]

Can you comment on why you chose the 0.4999 / 0.5001 numbers? Looks like that should be 0.5 in both cases?

[stefanv]

This is a fundamental issue in ndimage that we never explicitly sorted out, and it comes down to two different ways of seeing data.

The first is to see pixels as blocks or bins, with the index value indicating either the start or center (a design decision):

| 0 | 1 | 2 | 3 | ...

The second is to see the pixel values as data points:

0 1 2 3 ...
. . . . .

The latter case is what is currently implemented, whereby "wrap" mode moves past 3 straight into 0. In the earlier model, however, Moving 0.5 past 3 could mean either being at 3.5 or at 4 (i.e. -1) depending on how you choose the center position.

We need to make a call on this, and then have consistent behavior all round.

[rgommers]

Guess the history and motivation for (2) needs some more explanation. I have always thought of (1) when using ndimage and considered wrap mode broken. Both images and gridded data require (1), not (2).

[stefanv]

The question, I guess, is just whether wrap mode should implement:

0 1 2 3/0 1 2 3

or

0 1 2 3 0 1 2 3

My guess would be the latter?

[rgommers]

I'd guess the former. The reason being that the location has meaning in n-D (for images, gridded data). So for

0 1 2 3
4 5 6 7

the 4 is a neighbor of 0 and 5, not of 3 and 5.

If you want the latter behavior, simply reshape to 1-D?

[stefanv]

I'm not sure I follow. Actually, I'm now pretty sure that behavior 1 cannot be implemented, because it would mean that 3 and 0 are the same data-point (but they may have distinct values).

In your example, the padding according to suggestion 2 above would look like this:

6 7 | 4 5 6 7 | 4 5
-------------------
2 3 | 0 1 2 3 | 0 1
6 7 | 4 5 6 7 | 4 5
-------------------
2 3 | 0 1 2 3 | 0 1

with all data spaced one unit apart.

[jjhelmus]

I would vote for the later. This seem to be what the the authors of Ticket #1520 and #796 were expecting. The downside of this is that we would also have to agree on how to extrapolate values in the space between the last index (n) and the first index of the wrap array (n+1).

Maybe we should bring this up on the scipy-dev list to get a few more voices in on the discussion? I can write up a summary of the two options and present them to the list.

As for the 0.4999 / 0.5001, they are present because there are some edge cases where you get indices that are out of range only due to floating point math rounding. Ticket #1378 is a manifestation of this problem, math.sin(180.0 * pi / 180.0) is ~1.2 e -16 and results in an index of 2 + some_tiny_value which is marked out of range.

[stefanv]

The interpolation is not a problem, since you have two data-points with a single unit spacing between them, just like everywhere else.

We should also carefully consider mirror mode, which may currently be implemented as:

0 1 2 3 x 3 2 1

instead of

0 1 2 3 3 2 1 0

(i.e. two unit spacings between the N and N+1 values).

[stefanv]

The out of bounds issue is really tricky. Typically, you want the last data-point to still have the extent of a pixel, just like the rest of the values, but that's not really feasible. The problem could be worked around by making sure (like zoom does) that you never go outside the (0, N-1) interval -- although we should build in towards-zero rounding around a bracket of a certain size (1e-10 or so) there to handle this problem.

[jjhelmus]

I've been trying to understand the various boundary filling modes and created this document to help keep track of the varous modes. Can we agree that this is how we should be doing 'wrap' mode? It is consistent with the short documentation in the Scipy documentation, http://docs.scipy.org/doc/scipy/reference/tutorial/ndimage.html

Boundary filling modes

There are five modes for filling in data outside of borders of inputted data, which are described below.
The value/index examples are for a length 4 array with values 0, 1, 2, 4.

Constant mode

Points with an index, including fraction indices, less than 0 or greater than
len - 1 are replaced with cval.

value: ...| cval | cval | 0 | 1 | 2 | 3 | cval | cval | ...
index: ...|  -2  |  -1  | 0 | 1 | 2 | 3 |   4  |  4   | ...

Nearest mode

Points with an index, including fraction indices, less than 0 or greater than
len - 1 are replaced with the nearest valid point (the first or last point).

value: ... | 0 | 0 | 0 | 1 | 2 | 3 | 3 | 3 | ...
index: ... |-2 |-1 | 0 | 1 | 2 | 3 | 4 | 5 | ... 

Reflect mode

Points with an index less than 0 or greater than len - 1 are reflected with a
repeat of the first/last point. Fraction indicies are interpolated.

value: ... |  3 |  2 |  1 |  0 | 0 | 1 | 2 | 3 | 3 | 2 | 1 | 0 | ...
index: ... | -4 | -3 | -2 | -1 | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | ...

Mirror mode

Points with an index less than 0 or greater than len - 1 are reflected
with no repeat of the first/last point. Fraction indicies are interpolated.

value: ... |  2 |  3 |  2 |  1 | 0 | 1 | 2 | 3 | 2 | 1 | 0 | 1 | ...
index: ... | -4 | -3 | -2 | -1 | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | ...

Wrap Mode

Points with an index less than 0 or greater than len-1 are wrapped by
prepending/appending a copy of the array. Fraction indicies are interpolated.

value: ... |  0 |  1 |  2 |  3 | 0 | 1 | 2 | 3 | 0 | 1 | 2 | 3 | ...
index: ... | -4 | -3 | -2 | -1 | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | ...

Currently the first and last point are taken to be the same data point:

value: ... |  2 |  3 |  1 |  2 | 0 / 3 | 1 | 2 | 3 / 0 | 1 | 2 | 0 | 1 | ...
index: ... | -4 | -3 | -2 | -1 |   0   | 1 | 2 |   3   | 4 | 5 | 6 | 7 | ...

With this pull all of the above modes work as described with integer shifts. Fractional indices fail at the boundaries for the reflect and wrap modes because the repeated values are not used for interpolation. I'm looking into this but if anyone can help I would appreaciate it.

[stefanv]

The description of the modes looks correct; but the 0 / 3 combination of points cannot be (I assume you were just documenting the current state--but if that is the current state it needs to be fixed).

[nfaggian]

Have you considered the behaviour of Matlabs image processing toolbox? Might make sense to be, at least, consistent in that regard. For what it's worth, I agree with Stefans points about the representation of pixels, it's a really common thing to ignore coordinate systems but not always the best thing to do. I guess treating things as an evenly spaced regular grid is an implicit thing here but it would be great to see that be described in code.

Also, consider the possibility of more complex boundary extrapolation algorithms, I think it would be nice to allow the user to define their own whacky way of extrapolating. In weather fields we do things like a local mean extrapolation.

[rgommers]

@Stefan: after reading Jonathan's description I now see that in your option 1, the 0/3 meant 0 and/or 3 for the same point. Since that indeed doesn't make much sense, I had interpreted your slash as a line break. So never mind my earlier comment.

[stefanv]

@rgommers I realised afterwards that my description was vague; sorry about that. So, then it looks like we all agree and can move forward with a solution. @jjhelmus feel free to write this dialogue up, and to use it as motivations for future PRs.

[charris]

Considering pixels as bins with the indexes at the center of the pixels is pretty standard, IDL and Matplotlib do this, for instance. It has some advantages, nice behavior under reflections with the consequent symmetries, centroids behave nicely, and rounding can be used to find which pixel a point is contained in.

[stefanv]

It is, but it also doesn't make much sense in N-dimensions where N > 3. Since both models are fairly common, do you agree that it makes sense in this case to just see the pixels as data points in space?

[charris]

Sure, they are data points in space, but where are they located and what do they represent? For pixels and voxels they are the value of an integral over a region and locating them at the center is in some sense the best estimate of the value at the center.

But apart from that, the choice is fairly arbitrary. The extension modes are related to the descrete cosine transforms and the fft. Reflect mode is DCT-II, mirror mode is DCT-I, the wrap mode is the FFT. For general data points in space that is probably a good way to look at things.

[charris]

And I'll add that the reflect mode, nee pixel centers, is also related to the trapazoidal rule. Generally there are numerical advantages that derive from the symmetry.

[jjhelmus]

I now have a better idea how to fix this issue. map_coordinate is being used to find the location of the shifted point, but not the points surrounding that point which are needed during the spline interpolation. I have to think a bit longer on what the best method for doing this should be.

Until then please hold off on the including this pull request.

[jjhelmus]

Closing this pull request for the moment. I don't have time to work on a proper fix for these issue right now.

[stefanv]

I did some more digging, and will try to summarize the two modes of thinking about wrapping on the boundary:

Do we perceive the first and last datapoints to be on the image boundary, i.e. that the last and first value—when wrapping—lie essentially on top of (very close to) one another. Or do we imagine a regular grid of measurement values that is extended further into a larger regular grid, with values depending on the mode of wrapping. ndimage implements the prior understanding, correctly.

This is in line with viewing the data as a periodic signal, but in contrast with a view of the data as pixels. It basically comes down to where you define the boundary of your image to be: does the boundary run on top of the outside data points of the image, or does it run outside of those data points by a distance of 0.5.

Both concepts are potentially valid, depending on what is desired, although I would argue that the currently implemented version—that the boundary lies on the outside data points, is valid and should remain the default.

We can perhaps consider adding a new wrapping mode, grid-wrap or some such. This notion also clarifies the difference between mirror and reflect, which becomes mirror and grid-mirror.