rustyqlib/core/interpolation/
bilinear.rs1use super::linear::bracket;
5use crate::core::errors::RustyQLibError;
6
7#[derive(Debug, Clone)]
11pub struct BilinearGrid {
12 xs: Vec<f64>,
13 ys: Vec<f64>,
14 z: Vec<Vec<f64>>,
15}
16
17impl BilinearGrid {
18 pub fn new(xs: &[f64], ys: &[f64], z: &[Vec<f64>]) -> Result<Self, RustyQLibError> {
19 if xs.len() < 2 || ys.len() < 2 {
20 return Err(RustyQLibError::invalid_input("bilinear", "need at least a 2 x 2 grid"));
21 }
22 if xs.windows(2).any(|w| w[1] <= w[0]) || ys.windows(2).any(|w| w[1] <= w[0]) {
23 return Err(RustyQLibError::invalid_input("bilinear", "grid axes must be strictly increasing"));
24 }
25 if z.len() != xs.len() || z.iter().any(|row| row.len() != ys.len()) {
26 return Err(RustyQLibError::invalid_input("bilinear", "z must be an xs.len() x ys.len() matrix"));
27 }
28 Ok(BilinearGrid { xs: xs.to_vec(), ys: ys.to_vec(), z: z.to_vec() })
29 }
30
31 pub fn eval(&self, x: f64, y: f64) -> f64 {
32 let (i, wx) = bracket(&self.xs, x);
33 let (j, wy) = bracket(&self.ys, y);
34 let z00 = self.z[i - 1][j - 1];
35 let z01 = self.z[i - 1][j];
36 let z10 = self.z[i][j - 1];
37 let z11 = self.z[i][j];
38 z00 * (1.0 - wx) * (1.0 - wy)
39 + z01 * (1.0 - wx) * wy
40 + z10 * wx * (1.0 - wy)
41 + z11 * wx * wy
42 }
43}
44
45#[cfg(test)]
46mod tests {
47 use super::*;
48
49 #[test]
50 fn reproduces_a_bilinear_function_exactly() {
51 let f = |x: f64, y: f64| 2.0 + 3.0 * x + 4.0 * y + 5.0 * x * y;
53 let xs = [0.0, 0.7, 1.5, 3.0];
54 let ys = [-1.0, 0.5, 2.0];
55 let z: Vec<Vec<f64>> =
56 xs.iter().map(|&x| ys.iter().map(|&y| f(x, y)).collect()).collect();
57 let grid = BilinearGrid::new(&xs, &ys, &z).unwrap();
58 for i in 0..=30 {
59 for j in 0..=30 {
60 let (x, y) = (i as f64 * 0.1, -1.0 + j as f64 * 0.1);
61 assert!((grid.eval(x, y) - f(x, y)).abs() < 1e-12, "({x}, {y})");
62 }
63 }
64 }
65
66 #[test]
67 fn clamps_outside_the_grid() {
68 let grid = BilinearGrid::new(
69 &[0.0, 1.0],
70 &[0.0, 1.0],
71 &[vec![1.0, 2.0], vec![3.0, 4.0]],
72 )
73 .unwrap();
74 assert_eq!(grid.eval(-5.0, -5.0), 1.0);
75 assert_eq!(grid.eval(9.0, 9.0), 4.0);
76 assert_eq!(grid.eval(0.5, -3.0), 2.0); }
78}