fix bug in aabb2 edge2 intersection length

del-geo-core/src/edge2.rs

@@ -1,5 +1,6 @@
 //! methods for 2D edge (line segment)
 
+use crate::aabb;
 use num_traits::AsPrimitive;
 
 pub fn length<T>(ps: &[T; 2], pe: &[T; 2]) -> T
@@ -142,36 +143,19 @@ fn test_nearest_point2() {
 
 pub fn intersection_length_against_aabb2(ps: &[f32; 2], pe: &[f32; 2], aabb2: &[f32; 4]) -> f32 {
     // 0 min, 1 max
-    let seg_range_x = [ps[0].min(pe[0]), ps[0].max(pe[0])];
+    let edge_range_x = [ps[0].min(pe[0]), ps[0].max(pe[0])];
     let aabb_range_x = [aabb2[0], aabb2[2]];
-    let seg_range_y = [ps[1].min(pe[1]), ps[1].max(pe[1])];
-    let aabb_range_y = [aabb2[1], aabb2[3]];
-
-    let dx = range_intersection_length(&aabb_range_x, &seg_range_x);
-    let dy = range_intersection_length(&aabb_range_y, &seg_range_y);
-    (dx * dx + dy * dy).sqrt()
-}
-
-pub fn range_intersection_length(r1: &[f32; 2], r2: &[f32; 2]) -> f32 {
-    // separeted
-    if r1[1] <= r2[0] || r1[0] >= r2[1] {
+    let Some(dx) = crate::range::intersection_length(&aabb_range_x, &edge_range_x) else {
         return 0f32;
     };
-
-    // contains
-    if r1[0] >= r2[0] && r1[1] <= r2[1] {
-        return r1[1] - r1[0];
-    } else if r1[0] <= r2[0] && r1[1] >= r2[1] {
-        return r2[1] - r2[0];
+    //
+    let edge_range_y = [ps[1].min(pe[1]), ps[1].max(pe[1])];
+    let aabb_range_y = [aabb2[1], aabb2[3]];
+    let Some(dy) = crate::range::intersection_length(&aabb_range_y, &edge_range_y) else {
+        return 0f32;
     };
-
-    if r1[1] > r2[1] {
-        r2[1] - r1[0]
-    } else if r2[1] > r1[1] {
-        r1[1] - r2[0]
-    } else {
-        (r1[1] - r2[0]).min(r1[1] - r1[0])
-    }
+    //
+    (dx * dx + dy * dy).sqrt()
 }
 
 #[test]
@@ -191,20 +175,20 @@ fn test_intersection_length_against_aabb2() {
     let length = intersection_length_against_aabb2(&ps, &pe, &aabb);
     assert_eq!(length, 0.0);
 
+    // outside AABB
+    let ps = [2.0, 0.0];
+    let pe = [3.0, 0.5];
+    let aabb = [2.0, 2.0, 3.0, 3.0];
+    let length = intersection_length_against_aabb2(&ps, &pe, &aabb);
+    assert_eq!(length, 0.0);
+
     // partially inside AABB
     let ps = [0.0, 0.0];
     let pe = [2.0, 2.0];
     let aabb = [1.0, 1.0, 3.0, 3.0];
     let length = intersection_length_against_aabb2(&ps, &pe, &aabb);
     assert!((length - SQRT_2).abs() < 1e-6);
 
-    // touching border
-    let ps = [1.0, 0.0];
-    let pe = [1.0, 3.0];
-    let aabb = [1.0, 1.0, 2.0, 2.0];
-    let length = intersection_length_against_aabb2(&ps, &pe, &aabb);
-    assert!((length - 1.0).abs() < 1e-6);
-
     // segment with coordinates in reverse order
     let ps = [3.0, 3.0];
     let pe = [1.0, 1.0];
@@ -226,3 +210,124 @@ fn test_intersection_length_against_aabb2() {
     let length = intersection_length_against_aabb2(&ps, &pe, &aabb);
     assert!((length - SQRT_2).abs() < 1e-6);
 }
+
+pub fn area_left_side_of_edge2(aabb2: &[f32; 4], ps: &[f32; 2], pe: &[f32; 2]) -> f32 {
+    let Some((tmin, tmax)) = aabb::intersections_against_line(aabb2, ps, &crate::vec2::sub(pe, ps))
+    else {
+        return 0f32;
+    };
+    // dbg!(tmin, tmax);
+    let qs = crate::vec2::axpy(tmin, &crate::vec2::sub(pe, ps), ps);
+    let qe = crate::vec2::axpy(tmax, &crate::vec2::sub(pe, ps), ps);
+    // dbg!(qe, qs);
+    if (qs[0] - aabb2[0]).abs() < f32::EPSILON {
+        // left in
+        assert!(aabb2[1] <= qs[1] && qs[1] < aabb2[3]);
+        if (qe[1] - aabb2[1]).abs() < f32::EPSILON {
+            // bottom out
+            let l_s_y = qs[1] - aabb2[1];
+            let l_e_x = qe[0] - aabb2[0];
+            let a = (aabb2[3] - aabb2[1]) * (aabb2[2] - aabb2[0]) - l_s_y * l_e_x * 0.5;
+            return a;
+        } else if (qe[0] - aabb2[2]).abs() < f32::EPSILON {
+            // right out
+            assert!(aabb2[1] <= qe[1] && qe[1] < aabb2[3]);
+            let l_s_y = aabb2[3] - qs[1];
+            let l_e_y = aabb2[3] - qe[1];
+            let a = (aabb2[2] - aabb2[0]) * (l_s_y + l_e_y) * 0.5;
+            return a;
+        } else if (qe[1] - aabb2[3]).abs() < f32::EPSILON {
+            // top out
+            let l_s_y = aabb2[3] - qs[1];
+            let l_e_x = qe[0] - aabb2[0];
+            return l_s_y * l_e_x * 0.5;
+        }
+    } else if (qs[1] - aabb2[1]).abs() < f32::EPSILON {
+        // bottom in
+        assert!(aabb2[0] <= qs[0] && qs[0] < aabb2[2]);
+        if (qe[0] - aabb2[0]).abs() < f32::EPSILON {
+            // left out
+            let l_s_x = qs[0] - aabb2[0];
+            let l_e_y = qe[1] - aabb2[1];
+            let a = l_s_x * l_e_y * 0.5;
+            return a;
+        }
+        if (qe[1] - aabb2[3]).abs() < f32::EPSILON {
+            // top out
+            assert!(aabb2[0] <= qe[0] && qe[0] < aabb2[2]);
+            let l_s_x = qs[0] - aabb2[0];
+            let l_e_x = qe[0] - aabb2[0];
+            let a = (aabb2[3] - aabb2[1]) * (l_s_x + l_e_x) * 0.5;
+            return a;
+        }
+        if (qe[0] - aabb2[2]).abs() < f32::EPSILON {
+            // right out
+            assert!(aabb2[1] <= qe[1] && qe[1] < aabb2[3]);
+            let l_s_x = aabb2[2] - qs[0];
+            let l_e_y = qe[1] - aabb2[1];
+            let a = (aabb2[3] - aabb2[1]) * (aabb2[2] - aabb2[0]) - (l_s_x * l_e_y) * 0.5;
+            return a;
+        }
+    } else if (qs[0] - aabb2[2]).abs() < f32::EPSILON {
+        // right in
+        assert!(aabb2[1] <= qs[1] && qs[1] < aabb2[3]);
+        if (qe[0] - aabb2[0]).abs() < f32::EPSILON {
+            // left out
+            assert!(aabb2[1] <= qe[1] && qe[1] < aabb2[3]);
+            let l_s_y = qs[1] - aabb2[1];
+            let l_e_y = qe[1] - aabb2[1];
+            let a = (aabb2[2] - aabb2[0]) * (l_s_y + l_e_y) * 0.5;
+            return a;
+        }
+    } else if (qs[1] - aabb2[3]).abs() < f32::EPSILON && (qe[1] - aabb2[1]).abs() < f32::EPSILON {
+        // up to down
+        assert!(aabb2[0] <= qs[0] && qs[0] < aabb2[2]);
+        assert!(aabb2[0] <= qe[0] && qe[0] < aabb2[2]);
+        let l_s_x = aabb2[2] - qs[0];
+        let l_e_x = aabb2[2] - qe[0];
+        let a = (aabb2[3] - aabb2[1]) * (l_s_x + l_e_x) * 0.5;
+        return a;
+    }
+    0f32
+}
+
+#[test]
+fn test_area_of_aabb2_left_side_of_edge2() {
+    let aabb2 = [0f32, 0.0, 1.0, 1.0];
+    {
+        let ps = [0.1f32, -0.5];
+        let pe = [0.7f32, 1.5];
+        let a = area_left_side_of_edge2(&aabb2, &ps, &pe); // down to up
+        assert!((a - 0.4).abs() < 1.0e-7);
+        let a = area_left_side_of_edge2(&aabb2, &pe, &ps); // up to down
+        assert!((a - 0.6).abs() < 1.0e-7);
+    }
+    {
+        let ps = [-0.5, 0.1];
+        let pe = [1.5, 0.7f32];
+        let a = area_left_side_of_edge2(&aabb2, &ps, &pe); // left to right
+        assert!((a - 0.6).abs() < 1.0e-7);
+        let a = area_left_side_of_edge2(&aabb2, &pe, &ps); // right to left
+        assert!((a - 0.4).abs() < 1.0e-7);
+    }
+    {
+        // left in
+        let ps = [-0.2, 0.5];
+        let a = area_left_side_of_edge2(&aabb2, &ps, &[0.5, -0.2]); // bottom out
+        assert!((a - 0.955).abs() < 1.0e-7); // 0.955=1.0-0.3*0.3*0.5
+        let a = area_left_side_of_edge2(&aabb2, &ps, &[0.5, 1.2]); // top out
+        assert!((a - 0.045).abs() < 1.0e-7); // // 0.045=0.3*0.3*0.5
+        let a = area_left_side_of_edge2(&aabb2, &ps, &[1.2, 0.7]); // right out
+        assert!((a - 0.4).abs() < 1.0e-7);
+    }
+    {
+        // bottom in
+        let ps = [0.5f32, -0.2];
+        let a = area_left_side_of_edge2(&aabb2, &ps, &[0.7, 1.2]); // top out
+        assert!((a - 0.6).abs() < 1.0e-7);
+        let a = area_left_side_of_edge2(&aabb2, &ps, &[-0.2, 0.5]); // left out
+        assert!((a - 0.045).abs() < 1.0e-7);
+        let a = area_left_side_of_edge2(&aabb2, &ps, &[1.2, 0.5]); // right out
+        assert!((a - 0.955).abs() < 1.0e-7);
+    }
+}

del-geo-core/src/range.rs

@@ -1,8 +1,11 @@
 /// Find the distance of two ranges, return None if they are overlapped
+///
 pub fn distance_to_range<Real>(a: (Real, Real), b: (Real, Real)) -> Option<Real>
 where
     Real: num_traits::Float,
 {
+    debug_assert!(a.0 <= a.1);
+    debug_assert!(b.0 <= b.1);
     if a.0 > b.1 {
         return Some(a.0 - b.1);
     }
@@ -11,3 +14,15 @@ where
     }
     None
 }
+
+pub fn intersection_length(r1: &[f32; 2], r2: &[f32; 2]) -> Option<f32> {
+    debug_assert!(r1[0] <= r1[1]);
+    debug_assert!(r2[0] <= r2[1]);
+    // separated
+    if r1[1] <= r2[0] || r1[0] >= r2[1] {
+        return None;
+    };
+    let vmin = r1[0].max(r2[0]);
+    let vmax = r1[1].min(r2[1]);
+    Some(vmax - vmin)
+}

del-geo-core/src/vec2.rs

@@ -112,6 +112,13 @@ pub fn orthogonalize(u: &[f32; 2], v: &[f32; 2]) -> [f32; 2] {
     [v[0] - t * u[0], v[1] - t * u[1]]
 }
 
+pub fn axpy<Real>(alpha: Real, x: &[Real; 2], y: &[Real; 2]) -> [Real; 2]
+where
+    Real: num_traits::Float,
+{
+    [alpha * x[0] + y[0], alpha * x[1] + y[1]]
+}
+
 // -------------------------------
 // below: about the Vec2 class
 pub struct XY<'a, Real> {