Curve Fitting - add calculation and chart plotting

Test files added to ensure correctness.

GitHub issue #268.

Cargo.toml

@@ -29,6 +29,7 @@ num = "0.4.3"
 num-traits = "0.2"
 num_cpus = "1.16.0"
 petgraph = { version = "0.6", default-features = false, features = ["stable_graph"] }
+plotters = { version = "0.3.7", default-features = false, features = ["image", "bitmap_encoder", "bitmap_backend", "line_series", "ttf"] }
 rand = { version = "0.8.5", default-features = false, features = ["std", "alloc", "small_rng"] }
 rayon = "1.10.0"
 rustc-hash = "2.0.0"

lib/rust/mmscenegraph/Cargo.toml

@@ -15,14 +15,16 @@ name = "bench"
 harness = false
 
 [dependencies]
+anyhow = { workspace = true }
 approx = { workspace = true }
 fastapprox = { workspace = true }
-rustc-hash = { workspace = true }
 log = { workspace = true }
+nalgebra = { workspace = true }
 num-traits = { workspace = true }
-rand = { workspace = true }
 petgraph = { workspace = true }
-nalgebra = { workspace = true }
+rand = { workspace = true }
+rustc-hash = { workspace = true }
 
 [dev-dependencies]
 criterion = { workspace = true }
+plotters = { workspace = true }

lib/rust/mmscenegraph/src/constant.rs

@@ -24,6 +24,7 @@ pub type AttrIndex = usize;
 pub type HashValue = u64;
 pub type NodeIndex = usize;
 pub type FrameValue = u32;
+pub type FrameTime = f64;
 pub type Real = f64;
 pub type Matrix44 = nalgebra::Matrix4<Real>;
 pub type Matrix33 = nalgebra::Matrix3<Real>;

lib/rust/mmscenegraph/src/math/curve_fit.rs

@@ -0,0 +1,94 @@
+//
+// Copyright (C) 2024 David Cattermole.
+//
+// This file is part of mmSolver.
+//
+// mmSolver is free software: you can redistribute it and/or modify it
+// under the terms of the GNU Lesser General Public License as
+// published by the Free Software Foundation, either version 3 of the
+// License, or (at your option) any later version.
+//
+// mmSolver is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with mmSolver.  If not, see <https://www.gnu.org/licenses/>.
+// ====================================================================
+//
+
+use anyhow::Result;
+use log::debug;
+
+use crate::constant::Real;
+use crate::math::line::curve_fit_linear_regression_type1;
+
+#[derive(Copy, Clone, Debug, PartialEq)]
+pub struct Point2 {
+    x: Real,
+    y: Real,
+}
+
+impl Point2 {
+    pub fn x(self) -> Real {
+        self.x
+    }
+
+    pub fn y(self) -> Real {
+        self.y
+    }
+}
+
+#[derive(Copy, Clone, Debug, PartialEq)]
+pub struct AngleRadian {
+    inner: Real,
+}
+
+impl AngleRadian {
+    pub fn value(self) -> Real {
+        self.inner
+    }
+
+    pub fn as_degrees(self) -> Real {
+        self.value().to_degrees()
+    }
+
+    pub fn as_direction(self) -> (Real, Real) {
+        (self.value().cos(), self.value().sin())
+    }
+}
+
+pub fn linear_regression(
+    x_values: &[Real],
+    y_values: &[Real],
+) -> Result<(Point2, AngleRadian)> {
+    let mut point_x = 0.0;
+    let mut point_y = 0.0;
+    let mut angle = 0.0;
+    curve_fit_linear_regression_type1(
+        &x_values,
+        &y_values,
+        &mut point_x,
+        &mut point_y,
+        &mut angle,
+    );
+    debug!("angle={angle}");
+
+    let dir_x = angle.cos();
+    let dir_y = angle.sin();
+    debug!("point_x={point_x} point_y={point_y}");
+    debug!("dir_x={dir_x} dir_y={dir_y}");
+
+    let point = Point2 {
+        x: point_x,
+        y: point_y,
+    };
+
+    let angle = AngleRadian {
+        inner: dir_y.atan2(dir_x),
+    };
+    debug!("angle={angle:?}");
+
+    Ok((point, angle))
+}

lib/rust/mmscenegraph/src/math/line.rs

@@ -18,6 +18,8 @@
 // ====================================================================
 //
 /// Line Regression.
+use log::debug;
+
 use crate::constant::Real;
 
 const EPSILON: Real = 1.0e-9;
@@ -87,6 +89,77 @@ fn impl_fit_line_to_points_type1(
     }
 }
 
+fn impl_combine_values_logic(
+    ok_a: bool,
+    ok_b: bool,
+    intercept_a: Real,
+    slope_a: Real,
+    slope_b: Real,
+    mean_x: Real,
+    mean_y: Real,
+    out_point_x: &mut Real,
+    out_point_y: &mut Real,
+    out_dir_x: &mut Real,
+    out_dir_y: &mut Real,
+) -> bool {
+    debug!("ok_a={ok_a} ok_b={ok_b}");
+
+    if ok_a && ok_b {
+        // Both linear regressions returned.
+
+        // Transpose the values, because of the swapped values given
+        // to the impl_fit_line_to_points_type1() function.
+        let slope_b = slope_b.recip().min(std::f64::MAX).copysign(slope_a);
+        assert_eq!(slope_a.signum(), slope_b.signum());
+
+        *out_point_x = mean_x;
+        *out_point_y = mean_y;
+
+        let slope = slope_a.signum() * (slope_a * slope_b).sqrt();
+        let angle = slope.atan();
+        *out_dir_x = angle.sin();
+        *out_dir_y = angle.cos();
+
+        debug!("slope_a={slope_a} slope_b={slope_b}");
+        debug!("angle={angle} slope={slope}");
+
+        true
+    } else if ok_a && !ok_b {
+        // The special case that the line is entirely vertical.
+        *out_point_x = intercept_a;
+        *out_point_y = mean_y;
+
+        let angle = slope_a.atan();
+        *out_dir_x = angle.sin();
+        *out_dir_y = angle.cos();
+
+        debug!("slope_a={slope_a} intercept_a={intercept_a}");
+        debug!("angle={angle}");
+
+        true
+    } else if !ok_a && ok_b {
+        // The special case that the line is entirely horizontal.
+
+        // Transpose the values, because of the swapped values given
+        // to the impl_fit_line_to_points_type1() function.
+        let slope_b = slope_b.recip().min(std::f64::MAX);
+
+        *out_point_x = mean_x;
+        *out_point_y = mean_y;
+
+        let angle = (-slope_b).atan();
+        *out_dir_x = angle.sin();
+        *out_dir_y = angle.cos();
+
+        debug!("slope_b={slope_b}");
+        debug!("angle={angle}");
+
+        true
+    } else {
+        false
+    }
+}
+
 pub fn fit_line_to_points_type1(
     x: &[Real],
     y: &[Real],
@@ -140,6 +213,60 @@ pub fn fit_line_to_points_type1(
     ok
 }
 
+pub fn curve_fit_linear_regression_type1(
+    x: &[Real],
+    y: &[Real],
+    out_mean_x: &mut Real,
+    out_mean_y: &mut Real,
+    out_angle: &mut Real,
+) -> bool {
+    let mut sum_x: Real = 0.0;
+    let mut sum_y: Real = 0.0;
+    let mut sum_xy: Real = 0.0;
+    let mut sum_x2: Real = 0.0;
+    let mut sum_y2: Real = 0.0;
+    let mut mean_x: Real = 0.0;
+    let mut mean_y: Real = 0.0;
+
+    let ok = impl_precompute_line_fit_data(
+        &x,
+        &y,
+        &mut sum_x,
+        &mut sum_y,
+        &mut sum_xy,
+        &mut sum_x2,
+        &mut sum_y2,
+        &mut mean_x,
+        &mut mean_y,
+    );
+    if !ok {
+        return false;
+    }
+
+    let mut intercept: Real = 0.0;
+    let mut slope: Real = 0.0;
+    let ok = impl_fit_line_to_points_type1(
+        sum_x,
+        sum_xy,
+        sum_x2,
+        mean_x,
+        mean_y,
+        &mut intercept,
+        &mut slope,
+    );
+    if ok {
+        *out_mean_x = mean_x;
+        *out_mean_y = mean_y;
+
+        *out_angle = slope.atan();
+        let angle_degree = out_angle.to_degrees();
+        println!(
+            "fit_line_to_points_type1_curve angle={out_angle} degrees={angle_degree}"
+        );
+    }
+    ok
+}
+
 /// Fits a Line to X and Y data, using "Linear Regression Type II".
 ///
 /// This function will correctly fit a data set regardless if the line
@@ -206,51 +333,19 @@ pub fn fit_line_to_points_type2(
         &mut slope_b,
     );
 
-    if ok_a && ok_b {
-        // Both linear regressions returned.
-
-        // Transpose the values, because of the swapped values given
-        // to the fit_line_to_points_type1() function.
-        slope_b = slope_b.recip().min(std::f64::MAX).copysign(slope_a);
-        assert_eq!(slope_a.signum(), slope_b.signum());
-
-        *out_point_x = mean_x;
-        *out_point_y = mean_y;
-
-        let slope = slope_a.signum() * (slope_a * slope_b).sqrt();
-        let angle = slope.atan();
-        *out_dir_x = angle.sin();
-        *out_dir_y = angle.cos();
-
-        true
-    } else if ok_a && !ok_b {
-        // The special case that the line is entirely vertical.
-        *out_point_x = intercept_a;
-        *out_point_y = mean_y;
-
-        let angle = slope_a.atan();
-        *out_dir_x = angle.sin();
-        *out_dir_y = angle.cos();
-
-        true
-    } else if !ok_a && ok_b {
-        // The special case that the line is entirely horizontal.
-
-        // Transpose the values, because of the swapped values given
-        // to the fit_line_to_points_type1() function.
-        slope_b = slope_b.recip().min(std::f64::MAX);
-
-        *out_point_x = mean_x;
-        *out_point_y = mean_y;
-
-        let angle = (-slope_b).atan();
-        *out_dir_x = angle.sin();
-        *out_dir_y = angle.cos();
-
-        true
-    } else {
-        false
-    }
+    impl_combine_values_logic(
+        ok_a,
+        ok_b,
+        intercept_a,
+        slope_a,
+        slope_b,
+        mean_x,
+        mean_y,
+        out_point_x,
+        out_point_y,
+        out_dir_x,
+        out_dir_y,
+    )
 }
 
 /// Approximates a perfectly straight line from a set of (ordered) line

lib/rust/mmscenegraph/src/math/mod.rs

@@ -19,6 +19,7 @@
 //
 
 pub mod camera;
+pub mod curve_fit;
 pub mod dag;
 pub mod line;
 pub mod line_intersect;