added ccd2 in del-geo-nalgebra

del-geo-core/src/aabb3.rs

@@ -24,36 +24,6 @@ where
 // ----------------------------------
 // Below: to method
 
-/// transform aabb to unit square (0,1)^3 while preserving aspect ratio
-/// return 4x4 homogeneous transformation matrix in **column major** order
-pub fn to_mat4_col_major_transf_into_unit_preserve_asp(aabb_world: &[f32; 6]) -> [f32; 16] {
-    let cntr = [
-        (aabb_world[0] + aabb_world[3]) * 0.5,
-        (aabb_world[1] + aabb_world[4]) * 0.5,
-        (aabb_world[2] + aabb_world[5]) * 0.5,
-    ];
-    let size = max_edge_size(aabb_world);
-    let a = 1f32 / size;
-    let b = -a * cntr[0] + 0.5;
-    let c = -a * cntr[1] + 0.5;
-    let d = -a * cntr[2] + 0.5;
-    [a, 0., 0., 0., 0., a, 0., 0., 0., 0., a, 0., b, c, d, 1.]
-}
-
-/// transform aabb to unit square (0,1)^3
-/// return 4x4 homogeneous transformation matrix in **column major** order
-pub fn to_mat4_col_major_transf_into_unit(aabb_world: &[f32; 6]) -> [f32; 16] {
-    let cntr = crate::aabb3::center(aabb_world);
-    let size = crate::aabb3::size(aabb_world);
-    let ax = 1f32 / size[0];
-    let ay = 1f32 / size[1];
-    let az = 1f32 / size[2];
-    let b = -ax * cntr[0] + 0.5;
-    let c = -ay * cntr[1] + 0.5;
-    let d = -az * cntr[2] + 0.5;
-    [ax, 0., 0., 0., 0., ay, 0., 0., 0., 0., az, 0., b, c, d, 1.]
-}
-
 // ----------------------------------
 
 pub fn scale<T>(aabb: &[T; 6], s: T) -> [T; 6]

del-geo-core/src/lib.rs

@@ -1,8 +1,3 @@
-#![allow(clippy::identity_op)]
-#![allow(clippy::erasing_op)]
-#![allow(clippy::needless_borrow)]
-#![allow(dead_code)]
-#![allow(unused_variables)]
 /*
 vec < mat < aabb < obb <
 line < ray < edge <
@@ -25,7 +20,6 @@ pub mod obb3;
 pub mod vec2;
 pub mod vec3;
 //
-pub mod ccd;
 pub mod edge;
 pub mod edge2;
 pub mod edge3;
@@ -35,6 +29,7 @@ pub mod mat3_row_major;
 pub mod mat3_sym;
 pub mod matn;
 pub mod obb2;
+pub mod polynomial_root;
 pub mod quaternion;
 pub mod spherical_harmonics;
 pub mod tet;

del-geo-core/src/mat3_row_major.rs

@@ -114,13 +114,13 @@ pub fn svd(m: &[f64; 9], nitr: usize) {
     let mut u = [0f64; 9];
     if determinant(&v) < 0. {
         // making right hand coordinate
-        v[0 * 3 + 2] *= -1.;
-        v[1 * 3 + 2] *= -1.;
+        v[2] *= -1.; // v[0,2] = v[0*3+2]
+        v[5] *= -1.; // v[1,2] = v[1*3+2]
         v[2 * 3 + 2] *= -1.;
         g[2] *= -1.;
-        u[0 * 3 + 2] *= -1.;
-        u[1 * 3 + 2] *= -1.;
-        u[2 * 3 + 2] *= -1.;
+        u[2] *= -1.; // u[0,2] = u[0*3+2]
+        u[5] *= -1.; // u[1,2] = u[1*3+2]
+        u[8] *= -1.; // u[2,2] = u[2*3+2]
     }
 
     let sql0 = crate::vec3::squared_norm(&u0);

del-geo-core/src/mat4_col_major.rs

@@ -1,5 +1,6 @@
 //! methods for 4x4 matrix
 
+use crate::aabb3::max_edge_size;
 use num_traits::AsPrimitive;
 use std::ops::AddAssign;
 
@@ -85,12 +86,13 @@ where
     ]
 }
 
-/// transformation converting normalized device coordinate (NDC) `[-1,+1]^2` to pixel coordinate
+/// transformation converting normalized device coordinate (NDC) `[-1,+1]^3` to pixel coordinate
 /// depth (-1, +1) is transformed to (0, +1)
 /// for example:
 ///     [-1,-1,-1] becomes (0, H, 0)
 ///     [+1,+1,+1] becomes (W, 0, 1)
-/// * Return
+///
+/// * Arguments
 ///     * `image_shape` - (width, height)
 pub fn from_transform_ndc2pix(img_shape: (usize, usize)) -> [f32; 16] {
     [
@@ -122,7 +124,6 @@ pub fn from_aabb3_fit_into_ndc_preserving_xyasp(aabb: &[f32; 6], asp: f32) -> [f
         let size = crate::aabb3::size(aabb);
         let lenx = size[0] / asp;
         let leny = size[1];
-        dbg!(size, lenx, leny);
         if lenx > leny {
             (lenx / 2f32, size[2] / 2f32)
         } else {
@@ -149,6 +150,36 @@ pub fn from_aabb3_fit_into_ndc_preserving_xyasp(aabb: &[f32; 6], asp: f32) -> [f
     ]
 }
 
+/// transform aabb to unit square (0,1)^3 while preserving aspect ratio
+/// return 4x4 homogeneous transformation matrix in **column major** order
+pub fn from_aabb3_fit_into_unit_preserve_asp(aabb_world: &[f32; 6]) -> [f32; 16] {
+    let cntr = [
+        (aabb_world[0] + aabb_world[3]) * 0.5,
+        (aabb_world[1] + aabb_world[4]) * 0.5,
+        (aabb_world[2] + aabb_world[5]) * 0.5,
+    ];
+    let size = max_edge_size(aabb_world);
+    let a = 1f32 / size;
+    let b = -a * cntr[0] + 0.5;
+    let c = -a * cntr[1] + 0.5;
+    let d = -a * cntr[2] + 0.5;
+    [a, 0., 0., 0., 0., a, 0., 0., 0., 0., a, 0., b, c, d, 1.]
+}
+
+/// transform aabb to unit square (0,1)^3
+/// return 4x4 homogeneous transformation matrix in **column major** order
+pub fn from_aabb3_fit_into_unit(aabb_world: &[f32; 6]) -> [f32; 16] {
+    let cntr = crate::aabb3::center(aabb_world);
+    let size = crate::aabb3::size(aabb_world);
+    let ax = 1f32 / size[0];
+    let ay = 1f32 / size[1];
+    let az = 1f32 / size[2];
+    let b = -ax * cntr[0] + 0.5;
+    let c = -ay * cntr[1] + 0.5;
+    let d = -az * cntr[2] + 0.5;
+    [ax, 0., 0., 0., 0., ay, 0., 0., 0., 0., az, 0., b, c, d, 1.]
+}
+
 // above: from method (making 4x4 matrix)
 // ----------------------------------------
 
@@ -175,9 +206,9 @@ where
     let b = [t[12], t[13], t[14]];
     let d = t[15];
     let c = [t[3], t[7], t[11]];
-    let e = Real::one() / (crate::vec3::dot(&c, &p) + d);
+    let e = Real::one() / (crate::vec3::dot(&c, p) + d);
     let ee = e * e;
-    let f = crate::vec3::add(&crate::mat3_col_major::mult_vec(&a, &p), &b);
+    let f = crate::vec3::add(&crate::mat3_col_major::mult_vec(&a, p), &b);
     [
         a[0] * e - f[0] * c[0] * ee,
         a[1] * e - f[1] * c[0] * ee,

del-geo-core/src/obb3.rs

@@ -146,7 +146,7 @@ pub fn nearest_to_point3<Real>(obb: &[Real; 12], p: &[Real; 3]) -> [Real; 3]
 where
     Real: num_traits::Float,
 {
-    if is_include_point(&obb, p, Real::zero()) {
+    if is_include_point(obb, p, Real::zero()) {
         return *p;
     }
     let (axes, hlen) = unit_axes_and_half_edge_lengths(obb);
@@ -187,8 +187,8 @@ where
     Real: num_traits::Float + std::fmt::Debug,
 {
     let axes = {
-        let (axes_i, _) = unit_axes_and_half_edge_lengths(&obb_i);
-        let (axes_j, _) = unit_axes_and_half_edge_lengths(&obb_j);
+        let (axes_i, _) = unit_axes_and_half_edge_lengths(obb_i);
+        let (axes_j, _) = unit_axes_and_half_edge_lengths(obb_j);
         [
             axes_i[0],
             axes_i[1],

del-geo-core/src/ccd.rs → del-geo-core/src/polynomial_root.rs

@@ -4,7 +4,7 @@ use num_traits::AsPrimitive;
 
 /// find root of quadratic function
 /// f(x) = c0 + c1*x + c2*x^2
-fn quadratic_root<T>(c0: T, c1: T, c2: T) -> Option<(T, T)>
+pub fn quadratic_root<T>(c0: T, c1: T, c2: T) -> Option<[T; 2]>
 where
     T: num_traits::Float + 'static + Copy + std::fmt::Debug,
     i64: AsPrimitive<T>,
@@ -30,7 +30,7 @@ where
         dbg!(c0, c1, c2, det, tmp, x1, x2);
     }
     assert!(x1 <= x2);
-    Some((x1, x2))
+    Some([x1, x2])
 }
 
 #[test]
@@ -46,7 +46,7 @@ fn test_quadratic_root() {
         let c0: f64 = x0 * x1 * c2;
         let res = quadratic_root(c0, c1, c2);
         assert!(res.is_some());
-        let (y0, y1) = res.unwrap();
+        let [y0, y1] = res.unwrap();
         assert!((x0 - y0).abs() < 1.0e-8);
         assert!((x1 - y1).abs() < 1.0e-8);
         //
@@ -89,7 +89,7 @@ where
             }
         } else {
             // quadratic function
-            if let Some((e0, e1)) = quadratic_root(c0, c1, c2) {
+            if let Some([e0, e1]) = quadratic_root(c0, c1, c2) {
                 let (e0, e1) = if e0 < e1 { (e0, e1) } else { (e1, e0) };
                 if e0 >= T::zero() && e0 <= t {
                     result.push(e0);
@@ -130,7 +130,7 @@ where
         Some(r)
     };
     // f'(x) = c1 + 2*c2*x + 3*c3*x^2
-    if let Some((e0, e1)) = quadratic_root(c1, two * c2, three * c3) {
+    if let Some([e0, e1]) = quadratic_root(c1, two * c2, three * c3) {
         assert!(e0 <= e1);
         if T::zero() <= e0 {
             // overlap [0,t] and [-\infty,e0]

del-geo-core/src/spherical_harmonics.rs

@@ -6,7 +6,7 @@ use std::f64::consts::PI;
 
 /// Calculate the normalization of the vector.
 #[inline]
-fn normalize(x: &mut f64, y: &mut f64, z: &mut f64) -> f64 {
+pub fn normalize(x: &mut f64, y: &mut f64, z: &mut f64) -> f64 {
     let r = f64::sqrt(*x * *x + *y * *y + *z * *z);
     let invr = 1.0 / r;
     *x *= invr;
@@ -18,7 +18,7 @@ fn normalize(x: &mut f64, y: &mut f64, z: &mut f64) -> f64 {
 /// Calculate the coefficients of the spherical harmonics for l <= 9 and store them in an array buffer.
 /// Try to access the coefficient Y_l^m by the index: base + l + m, where base = l^2.
 #[inline]
-fn sph_coeff_buffer(n: i8, x: f64, y: f64, z: f64) -> [f64; 100] {
+pub fn sph_coeff_buffer(n: i8, x: f64, y: f64, z: f64) -> [f64; 100] {
     let inv_pi = 1.0 / PI;
     let ep = Complex::new(x, y);
     let mut res = [0.0; 100];
@@ -314,7 +314,7 @@ fn factorial(n: u128) -> u128 {
 
 /// Calculate the coefficient of the term in the Legendre Polynomial.
 /// The term is P_l^m(x) = get_legendre_poly_term_coeff(l, m) * x^m.
-fn get_legendre_poly_term_coeff(func_order: u32, term_order: u32) -> f64 {
+pub fn get_legendre_poly_term_coeff(func_order: u32, term_order: u32) -> f64 {
     if (func_order - term_order) % 2 != 0 {
         0.0
     } else {
@@ -363,7 +363,7 @@ fn calculate_assoc_legendre_poly(l: u64, m: i64, x: f64) -> f64 {
 /// However, u128 is not enough for the factorial calculation, so the maximum l is 11.
 /// Now that it is related to multiply and divide in big factorials, it has large opmitization potential.
 /// But optimization based on combination is hard in both mathematically and programmatically.
-fn get_spherical_harmonics_coeff(l: i64, m: i64, x: f64, y: f64, z: f64) -> f64 {
+pub fn get_spherical_harmonics_coeff(l: i64, m: i64, x: f64, y: f64, z: f64) -> f64 {
     let m_abs = m.abs();
     assert!(l >= 0 && l >= m_abs);
     let r = f64::sqrt(x * x + y * y);
@@ -389,7 +389,7 @@ fn get_spherical_harmonics_coeff(l: i64, m: i64, x: f64, y: f64, z: f64) -> f64
 
 #[test]
 fn test_get_spherical_harmonics_coeff() {
-    let tmp = 1.0 / f64::sqrt(3.0);
+    // let tmp = 1.0 / f64::sqrt(3.0);
     let sph_coeff = sph_coeff_buffer(
         9,
         1.0 / f64::sqrt(3.0),

del-geo-core/src/view_projection.rs

@@ -29,8 +29,8 @@ impl Perspective {
 
     pub fn camera_translation(&mut self, asp: f32, cursor_dx: f32, cursor_dy: f32) {
         let mp = self.mat4_col_major(asp);
-        let sx = (mp[3 + 4 * 3] - mp[0 + 4 * 3]) / mp[0 + 4 * 0];
-        let sy = (mp[3 + 4 * 3] - mp[1 + 4 * 3]) / mp[1 + 4 * 1];
+        let sx = (mp[3 + 4 * 3] - mp[12]) / mp[0]; // (mp[3 + 4 * 3] - mp[0 + 4 * 3]) / mp[0 + 4 * 0];
+        let sy = (mp[3 + 4 * 3] - mp[13]) / mp[5]; // (mp[3 + 4 * 3] - mp[1 + 4 * 3]) / mp[1 + 4 * 1];
         self.cam_pos[0] -= sx * cursor_dx;
         self.cam_pos[1] -= sy * cursor_dy;
     }

del-geo-nalgebra/src/aabb2.rs

@@ -2,14 +2,31 @@
 
 use num_traits::AsPrimitive;
 
+pub fn from_vtx2vec<T>(vtx2vec: &[nalgebra::Vector2<T>]) -> [T; 4]
+    where
+        T: nalgebra::RealField + Copy,
+{
+    let mut aabb = [vtx2vec[0][0], vtx2vec[0][1], vtx2vec[0][0], vtx2vec[0][1]];
+    for xy in vtx2vec.iter().skip(1) {
+        aabb[0] = aabb[0].min(xy[0]);
+        aabb[1] = aabb[1].min(xy[1]);
+        aabb[2] = aabb[2].max(xy[0]);
+        aabb[3] = aabb[3].max(xy[1]);
+    }
+    aabb
+}
+
+// Above: from method
+// ------------------------------
+
 /// signed distance from axis-aligned bounding box
 /// * `pos_in` - where the signed distance is evaluated
 /// * `x_min` - bounding box's x-coordinate minimum
 /// * `x_max` - bounding box's x-coordinate maximum
 /// * `y_min` - bounding box's y-coordinate minimum
 /// * `y_max` - bounding box's y-coordinate maximum
 /// * signed distance (inside is negative)
-pub fn signed_distance_aabb<Real>(
+pub fn signed_distance<Real>(
     pos_in: nalgebra::Vector2<Real>,
     min0: nalgebra::Vector2<Real>,
     max0: nalgebra::Vector2<Real>,
@@ -26,16 +43,3 @@ where
     x_dist.max(y_dist)
 }
 
-pub fn from_vtx2vec<T>(vtx2vec: &[nalgebra::Vector2<T>]) -> [T; 4]
-where
-    T: nalgebra::RealField + Copy,
-{
-    let mut aabb = [vtx2vec[0][0], vtx2vec[0][1], vtx2vec[0][0], vtx2vec[0][1]];
-    for xy in vtx2vec.iter().skip(1) {
-        aabb[0] = aabb[0].min(xy[0]);
-        aabb[1] = aabb[1].min(xy[1]);
-        aabb[2] = aabb[2].max(xy[0]);
-        aabb[3] = aabb[3].max(xy[1]);
-    }
-    aabb
-}