Merge branch 'features/spline' into dev

Author: Axect

example_data/bspline_calculus.png

@@ -0,0 +1,69 @@
+use peroxide::fuga::*;
+use anyhow::Result;
+
+fn main() -> Result<()> {
+    // 1. Define the B-spline
+    let degree = 3;
+    let knots = vec![0f64, 1f64, 2f64, 3f64, 4f64];
+    let control_points = vec![
+        vec![0f64, 0f64],
+        vec![1f64, 2f64],
+        vec![2f64, -1f64],
+        vec![3f64, 1f64],
+        vec![4f64, -1f64],
+        vec![5f64, 2f64],
+        vec![6f64, 0f64],
+    ];
+
+    let spline = BSpline::clamped(degree, knots, control_points.clone())?;
+    let deriv_spline = spline.derivative()?;
+
+    // 2. Define the integrand for ∫y dx = ∫ y(t) * (dx/dt) dt
+    let integrand = |t: f64, s: &BSpline, ds: &BSpline| -> f64 {
+        let p = s.eval(t);
+        let dp = ds.eval(t);
+        let y_t = p.1;
+        let dx_dt = dp.0;
+        y_t * dx_dt
+    };
+
+    // 3. Perform numerical integration using self-contained Gauss-Legendre quadrature
+    let t_start = 0f64;
+    let t_end = 4f64;
+    
+    //let area = gauss_legendre_integrate(
+    //    |t| integrand(t, &spline, &deriv_spline),
+    //    t_start,
+    //    t_end,
+    //    64,
+    //);
+    let area = integrate(|t| integrand(t, &spline, &deriv_spline), (t_start, t_end), G7K15R(1e-4, 20));
+
+    // 4. Print the result
+    println!("The area under the B-spline curve (∫y dx) is: {}", area);
+
+    // For verification, let's plot the original curve
+    #[cfg(feature = "plot")]
+    {
+        let t = linspace(t_start, t_end, 200);
+        let (x, y): (Vec<f64>, Vec<f64>) = spline.eval_vec(&t).into_iter().unzip();
+        
+        let mut plt = Plot2D::new();
+        plt.set_title(&format!("Original B-Spline (Area approx. {:.4})", area))
+            .set_xlabel("x")
+            .set_ylabel("y")
+            .insert_pair((x.clone(), y.clone()))
+            .insert_pair((control_points.iter().map(|p| p[0]).collect(), control_points.iter().map(|p| p[1]).collect()))
+            .set_plot_type(vec![(0, PlotType::Line), (1, PlotType::Scatter)])
+            .set_color(vec![(0, "blue"), (1, "red")])
+            .set_legend(vec!["Spline", "Control Points"])
+            .set_style(PlotStyle::Nature)
+            .set_path("example_data/bspline_with_area.png")
+            .set_dpi(600)
+            .savefig()?;
+        
+        println!("Generated plot: example_data/bspline_with_area.png");
+    }
+
+    Ok(())
+}

example_data/bspline_derivative.png

@@ -0,0 +1,92 @@
+use peroxide::fuga::*;
+use anyhow::Result;
+
+fn main() -> Result<()> {
+    // 1. Define the B-spline
+    let degree = 3;
+    let knots = vec![0f64, 1f64, 2f64, 3f64, 4f64];
+    let control_points = vec![
+        vec![0f64, 0f64],
+        vec![1f64, 2f64],
+        vec![2f64, -1f64],
+        vec![3f64, 1f64],
+        vec![4f64, -1f64],
+        vec![5f64, 2f64],
+        vec![6f64, 0f64],
+    ];
+
+    let spline = BSpline::clamped(degree, knots.clone(), control_points.clone())?;
+
+    // 2. Get the derivative and integral
+    let deriv_spline = spline.derivative()?;
+    let integ_spline = spline.integral()?;
+
+    // 3. Evaluate points on all splines
+    let t = linspace(0f64, 4f64, 200);
+    let (x, y): (Vec<f64>, Vec<f64>) = spline.eval_vec(&t).into_iter().unzip();
+    let (dx, dy): (Vec<f64>, Vec<f64>) = deriv_spline.eval_vec(&t).into_iter().unzip();
+    let (ix, iy): (Vec<f64>, Vec<f64>) = integ_spline.eval_vec(&t).into_iter().unzip();
+
+    // 4. Plotting
+    #[cfg(feature = "plot")]
+    {
+        // Plot Original Spline and its control points
+        let control_x: Vec<f64> = control_points.iter().map(|p| p[0]).collect();
+        let control_y: Vec<f64> = control_points.iter().map(|p| p[1]).collect();
+
+        let mut plt = Plot2D::new();
+        plt.set_title("Original B-Spline")
+            .set_xlabel("x")
+            .set_ylabel("y")
+            .insert_pair((x.clone(), y.clone()))
+            .insert_pair((control_x, control_y))
+            .set_plot_type(vec![(0, PlotType::Line), (1, PlotType::Scatter)])
+            .set_color(vec![(0, "darkblue"), (1, "darkred")])
+            .set_legend(vec!["Spline", "Control Points"])
+            .set_style(PlotStyle::Nature)
+            .set_path("example_data/bspline_original.png")
+            .set_dpi(600)
+            .savefig()?;
+
+        // Plot Derivative (Hodograph) and its control points
+        let deriv_control_x: Vec<f64> = deriv_spline.control_points.iter().map(|p| p[0]).collect();
+        let deriv_control_y: Vec<f64> = deriv_spline.control_points.iter().map(|p| p[1]).collect();
+
+        let mut plt_deriv = Plot2D::new();
+        plt_deriv.set_title("B-Spline Derivative (Hodograph)")
+            .set_xlabel("dx/dt")
+            .set_ylabel("dy/dt")
+            .insert_pair((dx, dy))
+            .insert_pair((deriv_control_x, deriv_control_y))
+            .set_plot_type(vec![(0, PlotType::Line), (1, PlotType::Scatter)])
+            .set_color(vec![(0, "darkblue"), (1, "darkred")])
+            .set_legend(vec!["Derivative", "Derivative Control Points"])
+            .set_style(PlotStyle::Nature)
+            .set_path("example_data/bspline_derivative.png")
+            .set_dpi(600)
+            .savefig()?;
+
+        // Plot Integral and its control points
+        let integ_control_x: Vec<f64> = integ_spline.control_points.iter().map(|p| p[0]).collect();
+        let integ_control_y: Vec<f64> = integ_spline.control_points.iter().map(|p| p[1]).collect();
+
+        let mut plt_integ = Plot2D::new();
+        plt_integ.set_title("B-Spline Integral")
+            .set_xlabel("Integral of x")
+            .set_ylabel("Integral of y")
+            .insert_pair((ix, iy))
+            .insert_pair((integ_control_x, integ_control_y))
+            .set_plot_type(vec![(0, PlotType::Line), (1, PlotType::Scatter)])
+            .set_color(vec![(0, "darkblue"), (1, "darkred")])
+            .set_legend(vec!["Integral", "Integral Control Points"])
+            .set_style(PlotStyle::Nature)
+            .set_path("example_data/bspline_integral.png")
+            .set_dpi(600)
+            .savefig()?;
+    }
+
+    println!("Generated plots for original, derivative, and integral splines.");
+    println!("Please check 'example_data/bspline_original.png', 'example_data/bspline_derivative.png', and 'example_data/bspline_integral.png'.");
+
+    Ok(())
+}

example_data/bspline_integral.png

@@ -1206,6 +1206,76 @@ impl BSpline {
 
         B[0][p]
     }
+
+    /// Returns the derivative of the B-spline.
+    /// The derivative of a B-spline of degree p is a B-spline of degree p-1.
+    pub fn derivative(&self) -> Result<Self> {
+        if self.degree == 0 {
+            bail!("Cannot take the derivative of a B-spline with degree 0.");
+        }
+
+        let p = self.degree;
+        let c = self.control_points.len();
+
+        // The new control points (Q_i) are calculated from the old ones (P_i)
+        // Q_i = p * (P_{i+1} - P_i) / (t_{i+p+1} - t_{i+1})
+        let mut new_control_points = Vec::with_capacity(c - 1);
+        for i in 0..c - 1 {
+            let p_i = &self.control_points[i];
+            let p_i1 = &self.control_points[i + 1];
+            let denominator = self.knots[i + p + 1] - self.knots[i + 1];
+            if denominator == 0.0 {
+                new_control_points.push(vec![0.0; p_i.len()]);
+            } else {
+                let factor = p as f64 / denominator;
+                let q_i = p_i1
+                    .iter()
+                    .zip(p_i)
+                    .map(|(a, b)| factor * (a - b))
+                    .collect();
+                new_control_points.push(q_i);
+            }
+        }
+
+        // The new knot vector is the old one with the first and last knots removed.
+        let new_knots = self.knots[1..self.knots.len() - 1].to_vec();
+        let new_degree = self.degree - 1;
+
+        BSpline::open(new_degree, new_knots, new_control_points)
+    }
+
+    /// Returns the integral of the B-spline.
+    /// The integral of a B-spline of degree p is a B-spline of degree p+1.
+    /// The first control point of the integral is set to zero (integration constant).
+    pub fn integral(&self) -> Result<Self> {
+        let p = self.degree;
+        let c = self.control_points.len();
+        let dim = self.control_points[0].len();
+
+        // The new knot vector is the old one with the first and last knots duplicated.
+        let mut new_knots = self.knots.clone();
+        new_knots.insert(0, self.knots[0]);
+        new_knots.push(self.knots[self.knots.len() - 1]);
+
+        // The new control points (R_i) are calculated recursively.
+        // R_j = R_{j-1} + P_{j-1} * (t_{j+p} - t_{j-1}) / (p+1)
+        let mut new_control_points = vec![vec![0.0; dim]; c + 1]; // R_0 is the zero vector
+
+        for i in 1..=c {
+            let factor = (self.knots[i + p] - self.knots[i - 1]) / ((p + 1) as f64);
+            let prev_r = new_control_points[i - 1].clone();
+            let p_im1 = &self.control_points[i - 1];
+            let mut r_i = Vec::with_capacity(dim);
+            for j in 0..dim {
+                r_i.push(prev_r[j] + factor * p_im1[j]);
+            }
+            new_control_points[i] = r_i;
+        }
+
+        let new_degree = self.degree + 1;
+
+        BSpline::open(new_degree, new_knots, new_control_points)
+    }
 }
 
 impl Spline<(f64, f64)> for BSpline {