Add default to bilinear interpolation

[jungerm2]

The first commit fixes the typo bilinier to bilinear and can be easily merged in. Although, you might want to add a type alias for backward compatibility and eventually deprecate it.

The second enables the 2d bilinear interpolation strategy to return a default value if the query is out of bounds, like scipy's fill_value. This works as is, but there are a few limitations to the current implementation that I would like your feedback on:

• Returning a default value is currently handled by the allow_default boolean flag, meaning that there are some checks to ensure that both extrapolate and allow_default aren't both enabled simultaneously. Maybe a BoundaryConditions enum would be better?
• There's no way to specify what the default value should be, it simply uses Default::default() at the moment, but maybe a user would want the default to be something else. I tried enabling this by adding the Sd generic to the Bilinear struct and adding in a new field to capture what the default should be. I'm fairly new to Rust and have not yet been able to figure out all the correct trait bounds needed to pull this off.
• Finally, this currently causes jaggies as there's no smooth interpolation between the data and default out-of-bounds value. This is tricky to solve given how the crate is implemented. It would require finding the neighboring pixels/data points to the out-of-bounds query and interpolate between the default value and any in-bounds neighbors. Since the actual coordinates are stored as Vecs this indirection layer makes this not as trivial as I originally thought. However, for default coordinates (where the index to coord mapping is just a range) this could be done more easily.

[jonasBoss]

First of all thank you for the contributions.

can you create a separate PR for the typo fix? I would merge that right away, but the extrapolation behavior needs some more work which we can address in this PR.

Maybe a BoundaryConditions enum would be better?

yes! But I think it should be named ExtrapolationBehavior or something similar. So they don't get confused with the boundary conditions I am adding for CubicSpline strategy in #5

enum OutOfBoundsBehavior<T, D: Dimension>{
    /// Return an error
    Error,
   
    /// Extrapolate the function
    Extrapolate,

    /// Return nearest known value 
    Nearest,

    /// Return a default value
    Default(T),
    // or maybe:
    /// Return a default value for each data row
    Default(Array<T, D>)
}

It would require finding the neighboring pixels/data points to the out-of-bounds query

maybe we rename the is_in_X_range fn to X_bound_check and return an enum which also contains the index of the boundary:

enum BoundsCheck{
   InRange,
   BelowBounds(usize),
   AboveBounds(usize),
}

But that also does not solve the jaggies properly, as we need some out of bounds x value at which we assume the default value is reached. so there would be three different cases:

1. in bounds
2. out of bounds but not yet the default value
3. out of bounds and the default value

at which value do we transition from 2. to 3.?

[jungerm2]

Made a separate PR for the typo fix, see here.

For the OutOfBoundsBehavior:

• The error case is problematic as it currently errors if a single query point is out of bounds, which kills the whole operation. having it return Result<T> would be nice.
• I'm not sure how Extrapolate currently works but it seems to be a constant extrapolation, just repeating edge values. Maybe a finer distinction needs to be made here.
• It would be very useful to be able to specify a global default. A good use case here is image warping, where default would be a background color.

Here's how I'm currently solving the jaggies issue. I had to switch away from using this crate for the moment because of this problem:

pub fn warp_array3_into(
    mapping: &Mapping,
    data: &Array3<f32>,
    out: &mut Array3<f32>,
    background: Option<Array1<f32>>
) {
    let [out_c, out_h, out_w] = out.shape() else {
        unreachable!("Data is Array3, should have 3 dims!")
    };
    let [_data_c, data_h, data_w] = data.shape() else {
        unreachable!("Data is Array3, should have 3 dims!")
    };

    // Points is a Nx2 array of xy pairs
    let points = Array::from_shape_fn((out_w * out_h, 2), |(i, j)| {
        if j == 0 {
            i % out_w
        } else {
            i / out_w
        }
    });

    // Warp all points and determine indices of in-bound ones
    let (background, padding) = if background.is_some() {
        (background.unwrap(), 1.0)
    } else {
        (Array1::<f32>::zeros(*out_c), 0.0)
    };
    let warpd = mapping.warp_points(&points);
    let in_range_x = |x: f32| -padding <= x && x <= (*data_w as f32) - 1.0 + padding;
    let in_range_y = |y: f32| -padding <= y && y <= (*data_h as f32) - 1.0 + padding;
    let left_top = warpd.mapv(|v| v.floor());
    let right_bottom = left_top.mapv(|v| v + 1.0);
    let right_bottom_weights = &warpd - &left_top;
    
    // Perform interpolation into buffer
    let get_pix_or_bkg = |x: f32, y: f32| {
        if x < 0f32 || x >= *data_w as f32 || y < 0f32 || y >= *data_h as f32 {
            background.to_owned()
        } else {
            data.slice(s![.., y as i32, x as i32]).to_owned()
        }
    };

    let mut out_view = ArrayViewMut::from_shape(
        (*out_c, out_h*out_w), 
        out.as_slice_mut().unwrap()
    ).unwrap();

    (
        out_view.axis_iter_mut(Axis(1)),
        points.axis_iter(Axis(0)),
        warpd.axis_iter(Axis(0)),
        left_top.axis_iter(Axis(0)),
        right_bottom.axis_iter(Axis(0)),
        right_bottom_weights.axis_iter(Axis(0))
    ).into_par_iter().for_each(|(mut out_slice, _dst_p, src_p, lt, rb, rb_w)| {
        let [src_x, src_y] = src_p.to_vec()[..] else {unreachable!("XY pair has only two values")};

        if !in_range_x(src_x) || !in_range_y(src_y) {
            return
        }

        let [l, t] = lt.to_vec()[..] else {unreachable!("XY pair has only two values")};
        let [r, b] = rb.to_vec()[..] else {unreachable!("XY pair has only two values")};
        let [rw, bw] = rb_w.to_vec()[..] else {unreachable!("XY pair has only two values")};

        // Do the interpolation
        let (tl, tr, bl, br) = (
            get_pix_or_bkg(l, t),
            get_pix_or_bkg(r, t),
            get_pix_or_bkg(l, b),
            get_pix_or_bkg(r, b),
        );
        let top = (1.0 - rw) * tl + rw * tr;
        let bottom = (1.0 - rw) * bl + rw * br;
        let value = (1.0 - bw) * top + bw * bottom;
        out_slice.assign(&value)
    });
}

This is probably not the best way of doing this (again, I'm fairly new to rust) and it's more tailored and less general than this crate, but here the "in-bounds" detection actually depends on the background value. If background is specified then I perform the bilinear interpolation for points that are +/-1 out of bounds too (by making the padding value 1) as these may have neighbors that are in bounds. This solves the jaggies issue here.

The above code works well but is a bit slower than ndarray-interp. I'm not sure where the performance difference comes from or how the above could be improved. Do let me know if you have ideas.

[jungerm2]

One more thing, I noticed a few more typos in #5, specifically the changelog you wrote CubicSpline as QubicSpline. Cheers!

[jonasBoss]

• The error case is problematic as it currently errors if a single query point is out of bounds, which kills the whole operation. having it return Result<T> would be nice.

I am not sure if I understand you correctly. All the interpolation functions return Result<_, InterpolateError>. Of course this invalidates the whole query because we can not create a valid solution. Even if only one error occurs. We certainly do not want the return type to be Array<Result<Sd::Elem, InterpolateError>, D> that would preserve successful query points but wastes memory and is cumbersome to use.
Having a default out of bounds value is very useful in this regard. It allows a user to set a sentinel values e.g. f64::NAN as a default and then check for the errors after the fact.

• I'm not sure how Extrapolate currently works but it seems to be a constant extrapolation, just repeating edge values. Maybe a finer distinction needs to be made here.

It just uses the last two in bounds values and calculates the linear function with those. So it just extends the last in bounds segment.

From what I can tell about your issue I feel like the following is a good approach:
Add the default value for out of bounds queries with an enum as I suggested above.
To address the jagged edges you add padding of one pixel with the default value to the source data before creating the interpolator. By setting the x and y data manually you can now also choose how large the padding will be.

[jungerm2]

Agreed, Array<Result<Sd::Elem, InterpolateError>, D> would be cumbersome, and while a sentinel value could be a good idea it would be problematic in the jaggies/padding case as it would mean interpolating between data and say NAN. Maybe a simpler solution is to return the data alongside a boolean array of in/out-of bounds?

[jonasBoss]

I intend to keep the return type as is. While a boolean array can be useful I do not think it should be the default case. Given the ability to set a default value for out-of-bounds query points the user can construct the boolean array himself.

and also your problem can be solved without to much work:

let original: Array3<f32> = .. ;// the data
let mut padded_dim = original.raw_dim(); // use raw_dim to preserve the type information
padded_dim[0] += 2;
padded_dim[1] += 2;
let mut padded_original = Array::zeros(padded_dim); 

// set whatever background color we like

padded_original.slice_mut(s![1..-1, 1..-1, ..]).assign(original);

// we need to set axis manually as the padding changed the coordinates
// the coordinates now go from -1 to n+1
let mut x = Array::linspace(padded_dim[0], -1.0, padded_dim[0] as f32 - 1.0);
let mut y = Array::linspace(padded_dim[1], -1.0, padded_dim[1] as f32 - 1.0);

// by changing the first and last value of the axis we can change how big the padding is
// that changes how fast we fade to black at the edge of the image
x[0] -= 2.0; // fade to black over 3 pixel
x[padded_dim[0] - 1] += 2.0; 

let strat = Bilinear::new().extrapolate(OutOfBoundsBehavior::Default(0.0)); // or maybe set a unique sentinel like f32::NAN
let interp = Interp2D::builder(padded_original)
    .strategy(strat)
    .x(x)
    .y(y)
    .build()
    .unwrap();

let new_x_coords: Array2<f32> = .. ; // transformed x coordinates
let new_y_coords: Array2<f32> = .. ;

let transformed_image = interp.interp_array(&new_x_coords, &new_y_coords);

This might not be 100% correct, and it requires this PR to implement some way to set the out-of-bounds behavior.

Alternatively we can do it without this PR by setting the padding to two pixels on each side and setting extrapolate(true) this way the first pixel is the fade to black (or whatever background you like) and the second pixel makes sure we stay at that when extrapolating.

[jungerm2]

This is a standard use case to me, but clearly, I'm biased... The above solution requires allocating a whole new array which might be very costly depending on the size. It would be nice to have this be part of the API.

[jonasBoss]

I am happy to review Proposals and changes to better facilitate your use-case. As I mentioned above I support the idea of customizing the out of bounds behavior. But I won't be implementing these changes myself.
Even if I decide that your changes are not a good fit for the crate, you can still use them in a custom interpolation strategy and extend the functionality of Interp2D with traits.

---

Alternatively we can do it without this PR by setting the padding to two pixels on each side and setting extrapolate(true) this way the first pixel is the fade to black (or whatever background you like) and the second pixel makes sure we stay at that when extrapolating.

Does that proposal not solve your problem?

The above solution requires allocating a whole new array which might be very costly depending on the size.

If you are concerned about performance I assume you have thousands of images to process. I am also assuming they all have the same size, or can be batched into batches with the same size.

• directly load them into an buffer that already contains the padding
• you can't process them all at once, so reuse the buffer by constructing Interp2D from a ArrayView of the buffer and an ArrayView of the x and y axis.
• use Interp2D::interp_array_into to also reuse the result buffer (assuming that also always has the same size.)
• still not fast enough? use Interp2D::new_unchecked to avoid runtime checks for data lengths and monotonic axis.