gistools/util/interpolation/
lanczos.rs

1use super::{Interpolatable, get_distance};
2use crate::util::GetInterpolateValue;
3use core::f64::consts::PI;
4use libm::{fabs, sin};
5use s2json::{GetM, GetXY, GetZ};
6
7/// # Lanczos Interpolation
8///
9/// ## Description
10/// Perform interpolation using the Lanczos filter. This method uses a kernel-based approach
11/// to weigh contributions from nearby points, providing a balance between smoothing and sharpness.
12///
13/// ## Usage
14pub fn lanczos_interpolation<
15    M: Clone,
16    P: GetXY + GetZ,
17    R: GetM<M> + GetXY + GetZ,
18    V: Interpolatable,
19>(
20    point: &P,
21    ref_data: &[R],
22    get_value: GetInterpolateValue<R, V>,
23) -> V {
24    if ref_data.is_empty() {
25        return V::default();
26    }
27
28    let mut numerator = V::default();
29    let mut denom = V::default();
30
31    for ref_point in ref_data {
32        let weight = lanczos_kernel(get_distance(point, ref_point), 2.);
33        let mut value = get_value(ref_point);
34        value *= weight;
35        numerator += value;
36        denom += weight;
37    }
38
39    // Avoid division by zero
40    if denom == 0. {
41        return V::default();
42    }
43    numerator /= denom;
44
45    numerator
46}
47
48/// Lanczos kernel function - returns the weight based on the distance from the target point
49/// <https://en.wikipedia.org/wiki/Lanczos_resampling>
50pub fn lanczos_kernel(x: f64, a: f64) -> f64 {
51    if x == 0. {
52        return 1.; // sinc(0) = 1
53    }
54    if fabs(x) >= a {
55        return 0.; // Outside the kernel radius
56    }
57    let pi_x = PI * x;
58    (sin(pi_x) / pi_x) * (sin(pi_x / a) / (pi_x / a))
59}