EconForge / interpolation.py

BSD 2-Clause "Simplified" License
123 stars 35 forks source link

Add boundary conditions #69

Open albop opened 4 years ago

albop commented 4 years ago

Currently natural splines is the default for cubic interpolation.

thomasaarholt commented 3 years ago

I was about to create an issue asking for support for masking any pixels outside the boundary as np.nan. Is there any nice way to do that with interpolation.py or shall I write my own script for this?

albop commented 3 years ago

Currently, there is no support for that. Closest you can do is eval_linear( ..., extrap_mode='constant'). It will put zeros outside of the grid. This is an easy feature, though. And it should be easy to add it. The question is how. I see two easy options.

thomasaarholt commented 3 years ago

It was easier than expected - I ended up doing it like this:

mask = np.any((points >= image_shape) | (points < 0.), axis=1)
interpolated_points[mask] = np.nan

So the full function for warping an image became:

def interpolate_image_to_new_position(img: "(Y, X)", points: "(N, 2) or (Y, X, 2)", fill_value=np.nan):
    """Warp an image to new positions given by a list of coordinates that has the same length 
    as the image has pixels

    Parameters
    ----------
    img: Image of shape (Y, X)
    points: Array of points of shape (N, 2), where the last two indices are in traditional (X,Y) convention
    fill_value: Value of points interpolated outside the grid of the image.
        False or float or np.nan
        False follows interpolate.py behaviour.
    """
    # Grid probably becomes a linspace'd array:
    grid = ((0, img.shape[0]-1, img.shape[0]), (0, img.shape[1]-1, img.shape[1]))
    points = points.reshape((-1, 2))
    points = points[:,::-1] # Swap coordinates to (Y, X) convention
    interpolated_values = eval_linear(grid, img, points)
    if fill_value is not False and fill_value is not None:
        mask = np.any((points >= img.shape) | (points < 0.), axis=1)
        interpolated_values[mask] = fill_value
    return interpolated_values.reshape(img.shape)
albop commented 3 years ago

Yes, that would work. It's a workaround with a bit of memory allocation that wouldn't happen if the functionality is included in interpolation.py. So let's keep this issue open.

thomasaarholt commented 3 years ago

I think I suffered some sort of brain-fart during turning my code into the previous comment. While the idea was right, the code I pasted didn't run. I've updated it.