[help needed] How to do 2D (coordinates in a grid) -> 1D interpolation (elevation of said points)?

[robinmoussu]

Hello, I am trying to interpolate the (1D) elevations of points in a 2D grid. I cannot figure what I am supposed to do. I am pretty sure that this crate can do it, but the example section doesn't have any 2D -> 1D interpolation example that I could adapt.

Note: My maths background with splines is really low, so it may be the source of some non-understanding.

[robinmoussu]

As far as I understand the following code should more or less work

use cgmath::Vector2;

let coordinates: Vec<Vector2<f64>> = vec![
    Vector2::new(0.0, 0.0), Vector2::new(0.0, 0.1), Vector2::new(0.0, 0.2), Vector2::new(0.0, 0.3),
    Vector2::new(1.0, 0.0), Vector2::new(1.0, 0.1), Vector2::new(1.0, 0.2), Vector2::new(1.0, 0.3),
    Vector2::new(2.0, 0.0), Vector2::new(2.0, 0.1), Vector2::new(2.0, 0.2), Vector2::new(2.0, 0.3),
    Vector2::new(3.0, 0.0), Vector2::new(3.0, 0.1), Vector2::new(3.0, 0.2), Vector2::new(3.0, 0.3),
];  

let alti_array: Vec<f64> = vec![
    00.0, 01.0, 02.0, 03.0,
    10.0, 11.0, 12.0, 13.0,
    20.0, 21.0, 22.0, 23.0,
    30.0, 31.0, 32.0, 33.0,
];

let keys: Vec<Key<_, _>>;
for (c, altitude) in Iterator::zip(coordinates.iter(), alti_array.iter()) {
    keys.push(Key::new(c, altitude, Interpolation::Linear));
}
let spline = Spline::from_vec(keys);

assert_eq!(00.0, spline.clamped_sample(Vector2::new(0.0, 0.0)));
assert_eq!(33.0, spline.clamped_sample(Vector2::new(3.0, 0.3)));

assert_eq!(11.5, spline.clamped_sample(Vector2::new(1.0, 0.15)));
assert_eq!(16.5, spline.clamped_sample(Vector2::new(1.5, 0.15)));
assert_eq!(21.5, spline.clamped_sample(Vector2::new(2.0, 0.15)));

assert_eq!(16.0, spline.clamped_sample(Vector2::new(1.5, 0.1)));
assert_eq!(16.5, spline.clamped_sample(Vector2::new(1.5, 0.15)));
assert_eq!(17.0, spline.clamped_sample(Vector2::new(1.5, 0.2)));

I have enable cgmath support in my Cargo.toml:

splines = { version = "3.3.0", features = ["std", "impl-cgmath"] }

But I get the following error:

error[E0277]: can't compare `cgmath::vector::Vector2<f64>` with `cgmath::vector::Vector2<f64>`
   --> src/main.rs:216:35
    |
216 |     let spline = Spline::from_vec(keys);
    |                                   ^^^^ no implementation for `cgmath::vector::Vector2<f64> < cgmath::vector::Vector2<f64>` and `cgmath::vector::Vector2<f64> > cgmath::vector::Vector2<f64>`
    |
    = help: the trait `std::cmp::PartialOrd` is not implemented for `cgmath::vector::Vector2<f64>`
    = note: required because of the requirements on the impl of `std::cmp::PartialOrd` for `&cgmath::vector::Vector2<f64>`
    = note: required by `splines::spline::Spline::<T, V>::from_vec`

It feels very strange that cgmath::Vector2 doesn't works out-of-the box. What I am missing?

[hadronized]

let keys: Vec<Key<_, _>>;

Where do you assign the value for keys?

[robinmoussu]

It's a typo. I assume it should have been (I'm not on my computer atm):

let keys: Vec<Key<_, _>> = Vec::new();

The keys are pushed in it just bellow.

[DefinitelyNotKornflex]

I might have not read the doc thoroughly enough, but I don't see how this library could handle this use-case. The sample method calls search_lower_cp which performs a binary search of highest lower key than the argument (thus finding which segment of the spline is relevant), which is why Key::T must implement PartialOrd. This trait is not implemented, probably because there a several ways of ordering two dimentional vectors, none of which : you might for instance say (a,b) > (c,d) if a>c or a = c and b>c, or compare the norms of the vectors ; either way, I doubt the behavior of the spline would be relevant.

If you have control points on each point of a regular grid, you might be able to implement a bilinear or bicubic interpolation by implementing the Interpolate trait for Spline (I'm not sure that'd be pretty, you'd probably better look for a bilinear or bicubic interpolation implementation elsewhere).

Otherwise, you should look into bivariate interpolation algorithms. The csaps crate might have what you need.

[hadronized]

Yeah, that’s an interesting need. Right now keys must implement a 1D-like interface. It would be interesting to generalize it so that we can build iso-surface in 2D / 3D / 4D.

[robinmoussu]

What puzzled me was that traits are implemented for cgmath::Vector{2,3,4}, but it feels that currently they are unusable.

[hadronized]

What do you mean unusable? They are, I use them. However, the traits that need to be implemented is on value types. The key type right now doesn’t require much to be able to be an interpolator.

Right now, key types are very limited and the traits used are a bit messy. For instance, the Trigo trait only has two implementors: f32 and f64.

But I see the need to get a better interface for keys. Right now, for sampling, we need:

• Interpolator types must implement Additive + One + Trigo + Mul<T, Output = T> + Div<T, Output = T> + PartialOrd. The situation is messy, as all this should be hidden behind a trait (or trait alias when supported).
• Value types must implement Additive + Interpolate<T>.

Maybe a small cleanup is required so that we can support nD sampling. :)

[DefinitelyNotKornflex]

@robinmoussu The (at least one of the) unimplemented trait on cgmath vectors is PartialOrd, which prevents it from being used as an interpolator (key). The traits implemented by the cgmath feature let you use these vectors as values for splines (Additive + Interpolate).

@phaazon Another problem is, with the current interface, if I have a type for which Trigo doesn't make much sense, and want to use a bézier spline for interpolation, I still have to provide a Trigo implementation (I mean, it doesn't hurt to define the methods and use the unreachable macro, but that's kinda awful).

I don't see a better solution to that than define different versions of the sample methods depending on the traits available, which isn't very cool either :/

[hadronized]

Hm. Yeah I see what you mean. That property arises from the possibility to change the type of interpolation on a segment basis.

I’ll try to think about it a bit and cleanup the code. I think it should be doable.

[hadronized]

I’m moving the project to sourcehut. Feel free to reopen there if you want to — issues / mailing list. Thank you.