bsplines 0.0.1-alpha.8

N-dimensional B-spline curves and their derivatives built on top of nalgebra
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71

use nalgebra::SVD;

use crate::{
    curve::{
        knots::Knots,
        parameters::Parameters,
        points::{DataPoints, Points},
    },
    types::MatD,
};

pub fn interpolate(knots: &Knots, points: &DataPoints, params: &Parameters) -> MatD {
    let p = knots.degree();
    let m = points.segments();
    let n = m;

    let Ubar = params.vector();

    let mut Nmat = MatD::zeros(points.count(), points.count());
    for i in 0..=n {
        for g in 0..=m {
            Nmat[(g, i)] = knots.evaluate(0, i, p, Ubar[g]);
        }
    }

    let svd = SVD::new(Nmat, true, true);
    let mat = points.matrix().transpose();
    svd.solve(&mat, f64::EPSILON.sqrt()).unwrap().transpose()
}

#[cfg(test)]
mod tests {
    use approx::assert_relative_eq;
    use nalgebra::dmatrix;

    use crate::curve::{knots, knots::Method::Averaging, parameters, parameters::Method::ChordLength};

    use super::*;

    #[test]
    fn linear() {
        let points = DataPoints::new(dmatrix![
            1., 2., 3., 4.;
            1., 2., 3., 4.;
        ]);

        let params = parameters::generate(&points, ChordLength);
        let knots = knots::generate(1, points.segments(), &params, Averaging).unwrap();

        assert_eq!(interpolate(&knots, &points, &params), points.matrix);
    }

    #[test]
    fn quadratic() {
        let points = DataPoints::new(dmatrix![
            1., 2., 3., 4.;
            1., 2., 3., 4.;
        ]);
        let params = parameters::generate(&points, ChordLength);
        let knots = knots::generate(2, points.segments(), &params, Averaging).unwrap();

        assert_relative_eq!(
            interpolate(&knots, &points, &params),
            dmatrix![
                1., 1.75, 3.25, 4.;
                1., 1.75, 3.25, 4.;
            ],
            epsilon = f64::EPSILON.sqrt()
        );
    }
}