diff --git a/geo/src/algorithm/offset/cross_product.rs b/geo/src/algorithm/offset/cross_product.rs index 6b9f2ef9b..36d3f2d6b 100644 --- a/geo/src/algorithm/offset/cross_product.rs +++ b/geo/src/algorithm/offset/cross_product.rs @@ -1,7 +1,7 @@ use crate::CoordFloat; use geo_types::Coord; -/// 2D "Cross Product" +/// The signed magnitude of the 3D "Cross Product" assuming z ordinates are zero /// /// > Note: `cross_prod` is already defined on `Point`... but that it seems to be /// > some other operation on 3 points @@ -40,7 +40,7 @@ use geo_types::Coord; /// = | (a_x·b_y - a_y·b_x)·k | /// = a_x·b_y - a_y·b_x /// ``` -pub(super) fn cross_product(left: Coord, right: Coord) -> T +pub(super) fn cross_product_2d(left: Coord, right: Coord) -> T where T: CoordFloat, { @@ -53,7 +53,7 @@ mod test { use crate::Coord; // private imports - use super::cross_product; + use super::cross_product_2d; #[test] fn test_cross_product() { @@ -65,9 +65,9 @@ mod test { let ac = c - a; // expect the area of the parallelogram - assert_eq!(cross_product(ac, ab), 1f64); + assert_eq!(cross_product_2d(ac, ab), 1f64); // expect swapping will result in negative - assert_eq!(cross_product(ab, ac), -1f64); + assert_eq!(cross_product_2d(ab, ac), -1f64); // Add skew let a = Coord { x: 0f64, y: 0f64 }; @@ -78,8 +78,8 @@ mod test { let ac = c - a; // expect the area of the parallelogram - assert_eq!(cross_product(ac, ab), 1f64); + assert_eq!(cross_product_2d(ac, ab), 1f64); // expect swapping will result in negative - assert_eq!(cross_product(ab, ac), -1f64); + assert_eq!(cross_product_2d(ab, ac), -1f64); } } diff --git a/geo/src/algorithm/offset/line_intersection.rs b/geo/src/algorithm/offset/line_intersection.rs index 0d2faf653..ead1ec0fa 100644 --- a/geo/src/algorithm/offset/line_intersection.rs +++ b/geo/src/algorithm/offset/line_intersection.rs @@ -1,4 +1,4 @@ -use super::cross_product; +use super::cross_product::cross_product_2d; use crate::{CoordFloat, CoordNum}; use geo_types::Coord; @@ -60,15 +60,15 @@ where /// M·T=ac /// ``` /// -/// The determinant of the matrix `M` is the reciprocal of the cross product -/// of `ab` and `cd`. +/// Inverting the matrix `M` involves taking the reciprocal of the determinant +/// (the determinant is same as the of the [cross_product()] of `ab` and `cd`) /// /// ```text /// 1/(ab×cd) /// ``` /// /// Therefore if `ab×cd = 0` the determinant is undefined and the matrix cannot -/// be inverted this means the lines are either +/// be inverted. The lines are either /// a) parallel or /// b) collinear /// @@ -105,7 +105,7 @@ pub(super) fn line_intersection_with_parameter( b: &Coord, c: &Coord, d: &Coord, -) -> LineIntersectionWithParameterResult +) -> Option> where T: CoordFloat, { @@ -113,24 +113,33 @@ where let cd = *d - *c; let ac = *c - *a; - let ab_cross_cd = cross_product(ab, cd); - - if ab_cross_cd == num_traits::zero() { - // TODO: We can't tolerate this situation as it will cause a divide by - // zero in the next step. Even values close to zero are a problem, - // but I don't know how to deal with that problem jut yet - - // TODO: this is prevented anyway by the only use of this function. - todo!("") - } - - let t_ab = cross_product(ac, cd) / ab_cross_cd; - let t_cd = -cross_product(ab, ac) / ab_cross_cd; - let intersection = *a + ab * t_ab; - LineIntersectionWithParameterResult { - t_ab, - t_cd, - intersection, + // TODO: I'm still confused about how to use Kernel / RobustKernel; + // the following did not work. I need to read more code + // from the rest of this repo to understand. + // if Kernel::orient2d(*a, *b, *d) == Orientation::Collinear { + // note that it is sufficient to check that only one of + // c or d are colinear with ab because of how they are + // related by the original line string. + // TODO: The following line + // - Does not use the Kernel + // - uses an arbitrary threshold value which needs more thought + let ab_cross_cd = cross_product_2d(ab, cd); + if ::from(ab_cross_cd) + .unwrap() + .abs() + < num_traits::cast(0.0000001f64).unwrap() + { + // Segments are parallel or colinear + None + } else { + let t_ab = cross_product_2d(ac, cd) / ab_cross_cd; + let t_cd = -cross_product_2d(ab, ac) / ab_cross_cd; + let intersection = *a + ab * t_ab; + Some(LineIntersectionWithParameterResult { + t_ab, + t_cd, + intersection, + }) } } @@ -144,19 +153,49 @@ mod test { let b = Coord { x: 2f64, y: 2f64 }; let c = Coord { x: 0f64, y: 1f64 }; let d = Coord { x: 1f64, y: 0f64 }; - let LineIntersectionWithParameterResult { + if let Some(LineIntersectionWithParameterResult { t_ab, t_cd, intersection, - } = line_intersection_with_parameter(&a, &b, &c, &d); - assert_eq!(t_ab, 0.25f64); - assert_eq!(t_cd, 0.5f64); - assert_eq!( + }) = line_intersection_with_parameter(&a, &b, &c, &d) + { + assert_eq!(t_ab, 0.25f64); + assert_eq!(t_cd, 0.5f64); + assert_eq!( + intersection, + Coord { + x: 0.5f64, + y: 0.5f64 + } + ); + } else { + assert!(false) + } + } + + #[test] + fn test_intersection_colinear() { + let a = Coord { x: 3f64, y: 4f64 }; + let b = Coord { x: 6f64, y: 8f64 }; + let c = Coord { x: 7f64, y: 7f64 }; + let d = Coord { x: 10f64, y: 9f64 }; + if let Some(LineIntersectionWithParameterResult { + t_ab, + t_cd, intersection, - Coord { - x: 0.5f64, - y: 0.5f64 - } - ); + }) = line_intersection_with_parameter(&a, &b, &c, &d) + { + assert_eq!(t_ab, 0.25f64); + assert_eq!(t_cd, 0.5f64); + assert_eq!( + intersection, + Coord { + x: 0.5f64, + y: 0.5f64 + } + ); + } else { + assert!(false) + } } } diff --git a/geo/src/algorithm/offset/mod.rs b/geo/src/algorithm/offset/mod.rs index c9b7e5f26..1dd3a44d1 100644 --- a/geo/src/algorithm/offset/mod.rs +++ b/geo/src/algorithm/offset/mod.rs @@ -1,7 +1,7 @@ + mod cross_product; +mod slice_itertools; mod line_intersection; -mod offset_trait; -use cross_product::cross_product; -use line_intersection::{line_intersection_with_parameter, LineIntersectionWithParameterResult}; +mod offset_trait; pub use offset_trait::Offset; diff --git a/geo/src/algorithm/offset/offset_trait.rs b/geo/src/algorithm/offset/offset_trait.rs index 2af2c9faa..1603596d7 100644 --- a/geo/src/algorithm/offset/offset_trait.rs +++ b/geo/src/algorithm/offset/offset_trait.rs @@ -1,4 +1,6 @@ -use super::{cross_product, line_intersection_with_parameter, LineIntersectionWithParameterResult}; +use super::line_intersection::{line_intersection_with_parameter, LineIntersectionWithParameterResult}; +use super::slice_itertools::pairwise; + use crate::{ CoordFloat, // Kernel, @@ -6,6 +8,7 @@ use crate::{ Line, LineString, MultiLineString, + Polygon }; use geo_types::Coord; @@ -77,17 +80,6 @@ where } } -/// Iterate over a slice in overlapping pairs -/// -/// ```ignore -/// let items = vec![1, 2, 3, 4, 5]; -/// let actual_result: Vec<(i32, i32)> = pairwise(&items[..]).map(|(a, b)| (*a, *b)).collect(); -/// let expected_result = vec![(1, 2), (2, 3), (3, 4), (4, 5)]; -/// assert_eq!(actual_result, expected_result); -/// ``` -fn pairwise(iterable: &[T]) -> std::iter::Zip, std::slice::Iter> { - iterable.iter().zip(iterable[1..].iter()) -} impl Offset for LineString where @@ -112,52 +104,32 @@ where result.push(first_point); result.extend(pairwise(&offset_segments[..]).flat_map( |(Line { start: a, end: b }, Line { start: c, end: d })| { - let ab = *b - *a; - let cd = *d - *c; - let ab_cross_cd = cross_product(ab, cd); - // TODO: I'm still confused about how to use Kernel / RobustKernel; - // the following did not work. I need to read more code - // from the rest of this repo to understand. - // if Kernel::orient2d(*a, *b, *d) == Orientation::Collinear { - // note that it is sufficient to check that only one of - // c or d are colinear with ab because of how they are - // related by the original line string. - // TODO: The following line - // - Does not use the Kernel - // - uses an arbitrary threshold value which needs more thought - if ::from(ab_cross_cd) - .unwrap() - .abs() - < num_traits::cast(0.0000001f64).unwrap() - { - vec![*b] - } else { - // TODO: if we can inline this function we only need to - // calculate `ab_cross_cd` once - let LineIntersectionWithParameterResult { + match line_intersection_with_parameter(a, b, c, d) { + None => vec![*b], // colinear + Some(LineIntersectionWithParameterResult { t_ab, t_cd, intersection, - } = line_intersection_with_parameter(a, b, c, d); - - let zero = num_traits::zero::(); - let one = num_traits::one::(); - - let tip_ab = zero <= t_ab && t_ab <= one; - let fip_ab = !tip_ab; - let pfip_ab = fip_ab && t_ab > zero; - - let tip_cd = zero <= t_cd && t_cd <= one; - let fip_cd = !tip_cd; - - if tip_ab && tip_cd { - // TODO: test for mitre limit - vec![intersection] - } else if fip_ab && fip_cd && pfip_ab { - // TODO: test for mitre limit - vec![intersection] - } else { - vec![*b, *c] + }) => { + let zero = num_traits::zero::(); + let one = num_traits::one::(); + + let tip_ab = zero <= t_ab && t_ab <= one; + let fip_ab = !tip_ab; + let pfip_ab = fip_ab && t_ab > zero; + + let tip_cd = zero <= t_cd && t_cd <= one; + let fip_cd = !tip_cd; + + if tip_ab && tip_cd { + // TODO: test for mitre limit + vec![intersection] + } else if fip_ab && fip_cd && pfip_ab { + // TODO: test for mitre limit + vec![intersection] + } else { + vec![*b, *c] + } } } }, @@ -178,23 +150,25 @@ where } } + +// impl Offset for Polygon +// where +// T: CoordFloat, +// { +// fn offset(&self, distance: T) -> Self { +// // TODO: not finished yet... need to do interiors +// // self.interiors() +// // TODO: is the winding order configurable? +// self.exterior(); +// todo!("Not finished") +// } +// } + #[cfg(test)] mod test { - // crate dependencies use crate::{line_string, Coord, Line, MultiLineString, Offset}; - // private imports - use super::pairwise; - - #[test] - fn test_pairwise() { - let items = vec![1, 2, 3, 4, 5]; - let actual_result: Vec<(i32, i32)> = pairwise(&items[..]).map(|(a, b)| (*a, *b)).collect(); - let expected_result = vec![(1, 2), (2, 3), (3, 4), (4, 5)]; - assert_eq!(actual_result, expected_result); - } - #[test] fn test_offset_line() { let input = Line::new(Coord { x: 1f64, y: 1f64 }, Coord { x: 1f64, y: 2f64 }); diff --git a/geo/src/algorithm/offset/slice_itertools.rs b/geo/src/algorithm/offset/slice_itertools.rs new file mode 100644 index 000000000..79c77d8b8 --- /dev/null +++ b/geo/src/algorithm/offset/slice_itertools.rs @@ -0,0 +1,45 @@ +/// Iterate over a slice in overlapping pairs +/// +/// ```ignore +/// let items = vec![1, 2, 3, 4, 5]; +/// let actual_result: Vec<(i32, i32)> = pairwise(&items[..]).map(|(a, b)| (*a, *b)).collect(); +/// let expected_result = vec![(1, 2), (2, 3), (3, 4), (4, 5)]; +/// ``` +pub(super) fn pairwise(slice: &[T]) -> std::iter::Zip, std::slice::Iter> { + slice.iter().zip(slice[1..].iter()) +} + +/// Iterate over a slice and repeat the first item at the end +/// +/// ```ignore +/// let items = vec![1, 2, 3, 4, 5]; +/// let actual_result: Vec = wrap_one(&items[..]).cloned().collect(); +/// let expected_result = vec![1, 2, 3, 4, 5, 1]; +/// ``` +pub(super) fn wrap_one( + slice: &[T], +) -> std::iter::Chain, std::slice::Iter> { + slice.iter().chain(slice[..1].iter()) + //.chain::<&T>(std::iter::once(slice[0])) +} + +#[cfg(test)] +mod test { + use super::{pairwise, wrap_one}; + + #[test] + fn test_pairwise() { + let items = vec![1, 2, 3, 4, 5]; + let actual_result: Vec<(i32, i32)> = pairwise(&items[..]).map(|(a, b)| (*a, *b)).collect(); + let expected_result = vec![(1, 2), (2, 3), (3, 4), (4, 5)]; + assert_eq!(actual_result, expected_result); + } + + #[test] + fn test_wrap() { + let items = vec![1, 2, 3, 4, 5]; + let actual_result: Vec = wrap_one(&items[..]).cloned().collect(); + let expected_result = vec![1, 2, 3, 4, 5, 1]; + assert_eq!(actual_result, expected_result); + } +}