diff --git a/geo/src/algorithm/mod.rs b/geo/src/algorithm/mod.rs index bba8a95d2..9436beae4 100644 --- a/geo/src/algorithm/mod.rs +++ b/geo/src/algorithm/mod.rs @@ -176,6 +176,11 @@ pub use lines_iter::LinesIter; pub mod map_coords; pub use map_coords::{MapCoords, MapCoordsInPlace}; +/// Offset the edges of a geometry perpendicular to the edge direction, either +/// to the left or to the right depending on the sign of the specified distance. +pub mod offset_curve; +pub use offset_curve::OffsetCurve; + /// Orient a `Polygon`'s exterior and interior rings. pub mod orient; pub use orient::Orient; diff --git a/geo/src/algorithm/offset_curve/line_intersection.rs b/geo/src/algorithm/offset_curve/line_intersection.rs new file mode 100644 index 000000000..eaaebdaba --- /dev/null +++ b/geo/src/algorithm/offset_curve/line_intersection.rs @@ -0,0 +1,370 @@ +/// Definitions used in documentation for this module; +/// +/// - **line**: +/// - The straight path on a plane that +/// - extends infinitely in both directions. +/// - defined by two distinct points `a` and `b` and +/// - has the direction of the vector from `a` to `b` +/// +/// - **segment**: +/// - A finite portion of a `line` which +/// - lies between the points `a` and `b` +/// - has the direction of the vector from `a` to `b` +/// +/// - **ray**: +/// - A segment which extends infinitely in the forward direction +/// +/// - `Line`: the type [crate::Line] which is actually a **segment**. + +use super::vector_extensions::VectorExtensions; +use crate::{ + Coord, + CoordFloat, + CoordNum, + // algorithm::kernels::Kernel, + // algorithm::kernels::RobustKernel, + // Orientation +}; + +/// Used to encode the relationship between a **segment** and an intersection +/// point. See documentation for [LineIntersectionResultWithRelationships] +#[derive(PartialEq, Eq, Debug, Clone, Copy)] +pub(super) enum LineSegmentIntersectionType { + /// The intersection point lies between the start and end of the **segment** + /// + /// Abbreviated to `TIP` in original paper + TrueIntersectionPoint, + /// The intersection point is 'false' or 'virtual': it lies on the same + /// **line** as the **segment**, but not between the start and end points of + /// the **segment**. + /// + /// Abbreviated to `FIP` in original paper + + // Note: Rust does not permit nested enum declaration, so + // FalseIntersectionPointType has to be declared below. + FalseIntersectionPoint(FalseIntersectionPointType), +} + +/// These are the variants of [LineSegmentIntersectionType::FalseIntersectionPoint] +#[derive(PartialEq, Eq, Debug, Clone, Copy)] +pub(super) enum FalseIntersectionPointType { + /// The intersection point is 'false' or 'virtual': it lies on the same + /// **line** as the **segment**, and before the start of the **segment**. + /// + /// Abbreviated to `NFIP` in original paper (Negative) + /// (also referred to as `FFIP` in Figure 6, but i think this is an + /// error?) + BeforeStart, + /// The intersection point is 'false' or 'virtual': it lies on the same + /// **line** as the **segment**, and after the end of the **segment**. + /// + /// Abbreviated to `PFIP` in original paper (Positive) + AfterEnd, +} + + +/// Struct to contain the result for [line_segment_intersection_with_relationships] +#[derive(Clone, Debug)] +pub(super) struct LineIntersectionResultWithRelationships +where + T: CoordNum, +{ + pub ab: LineSegmentIntersectionType, + pub cd: LineSegmentIntersectionType, + pub intersection: Coord, +} + +/// Computes the intersection between two line segments; +/// a to b (`ab`), and c to d (`cd`) +/// +/// > note: looks like there is already `cartesian_intersect` as a private +/// > method in simplifyvw.rs. It uses the orient2d method of [Kernel], +/// > however it only gives a true/false answer and does not return the +/// > intersection point or parameters needed. +/// +/// We already have LineIntersection trait BUT we need a function that also +/// returns the parameters for both lines described below. The LineIntersection +/// trait uses some fancy unrolled code it seems unlikely it could be adapted +/// for this purpose. +/// +/// Returns the intersection point **and** parameters `t_ab` and `t_cd` +/// described below +/// +/// The intersection of segments can be expressed as a parametric equation +/// where `t_ab` and `t_cd` are unknown scalars : +/// +/// ```text +/// a + ab · t_ab = c + cd · t_cd +/// ``` +/// +/// > note: a real intersection can only happen when `0 <= t_ab <= 1` and +/// > `0 <= t_cd <= 1` but this function will find intersections anyway +/// > which may lay outside of the line segments +/// +/// This can be rearranged as follows: +/// +/// ```text +/// ab · t_ab - cd · t_cd = c - a +/// ``` +/// +/// Collecting the scalars `t_ab` and `-t_cd` into the column vector `T`, +/// and by collecting the vectors `ab` and `cd` into matrix `M`: +/// we get the matrix form: +/// +/// ```text +/// [ab_x cd_x][ t_ab] = [ac_x] +/// [ab_y cd_y][-t_cd] [ac_y] +/// ``` +/// +/// or +/// +/// ```text +/// M·T=ac +/// ``` +/// +/// 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. The lines are either +/// a) parallel or +/// b) collinear +/// +/// Pre-multiplying both sides by the inverted 2x2 matrix we get: +/// +/// ```text +/// [ t_ab] = 1/(ab×cd) · [ cd_y -cd_x][ac_x] +/// [-t_cd] [-ab_y ab_x][ac_y] +/// ``` +/// +/// or +/// +/// ```text +/// T = M⁻¹·ac +/// ``` +/// +/// Expands to: +/// +/// ```text +/// [ t_ab] = 1/(ab_x·cd_y - ab_y·cd_x)·[ cd_y·ac_x - cd_x·ac_y] +/// [-t_cd] [-ab_y·ac_x + ab_x·ac_y] +/// ``` +/// +/// Since it is tidier to write cross products, observe that the above is +/// equivalent to: +/// +/// ```text +/// [t_ab] = [ ac×cd / ab×cd ] +/// [t_cd] = [ - ab×ac / ab×cd ] +/// ``` + +fn line_segment_intersection_with_parameters( + a: Coord, + b: Coord, + c: Coord, + d: Coord, +) -> Option<(T, T, Coord)> +where + T: CoordFloat, +{ + let ab = b - a; + let cd = d - c; + let ac = c - a; + + let ab_cross_cd = ab.cross_product_2d(cd); + if T::is_zero(&ab_cross_cd) { + // Segments are exactly parallel or colinear + None + } else { + // Division by zero is prevented, but testing is needed to see what + // happens for near-parallel sections of line. + let t_ab = ac.cross_product_2d(cd) / ab_cross_cd; + let t_cd = -ab.cross_product_2d(ac) / ab_cross_cd; + let intersection = a + ab * t_ab; + + Some((t_ab, t_cd, intersection)) + } + + // TODO: + // The above could be replaced with the following, but at the cost of + // repeating some computation. + + // match RobustKernel::orient2d(*a, *b, *d) { + // Orientation::Collinear => None, + // _ => { + // 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((t_ab, t_cd, intersection)) + // } + // } +} + +/// This is a simple wrapper for [line_segment_intersection_with_parameters]; +/// Returns the intersection point as well as the relationship between the point +/// and each of the input line segments. See [LineSegmentIntersectionType] +pub(super) fn line_segment_intersection_with_relationships( + a: Coord, + b: Coord, + c: Coord, + d: Coord, +) -> Option> +where + T: CoordFloat, +{ + line_segment_intersection_with_parameters(a, b, c, d).map(|(t_ab, t_cd, intersection)| { + use FalseIntersectionPointType::{AfterEnd, BeforeStart}; + use LineSegmentIntersectionType::{FalseIntersectionPoint, TrueIntersectionPoint}; + LineIntersectionResultWithRelationships { + ab: if T::zero() <= t_ab && t_ab <= T::one() { + TrueIntersectionPoint + } else if t_ab < T::zero() { + FalseIntersectionPoint(BeforeStart) + } else { + FalseIntersectionPoint(AfterEnd) + }, + cd: if T::zero() <= t_cd && t_cd <= T::one() { + TrueIntersectionPoint + } else if t_cd < T::zero() { + FalseIntersectionPoint(BeforeStart) + } else { + FalseIntersectionPoint(AfterEnd) + }, + intersection, + } + }) +} + +#[cfg(test)] +mod test { + use super::{ + line_segment_intersection_with_parameters, line_segment_intersection_with_relationships, + FalseIntersectionPointType, LineIntersectionResultWithRelationships, + LineSegmentIntersectionType, + }; + use crate::{Coord}; + use FalseIntersectionPointType::{AfterEnd, BeforeStart}; + use LineSegmentIntersectionType::{FalseIntersectionPoint, TrueIntersectionPoint}; + + #[test] + fn test_line_segment_intersection_with_parameters() { + let a = Coord { x: 0f64, y: 0f64 }; + let b = Coord { x: 2f64, y: 2f64 }; + let c = Coord { x: 0f64, y: 1f64 }; + let d = Coord { x: 1f64, y: 0f64 }; + if let Some((t_ab, t_cd, intersection)) = + line_segment_intersection_with_parameters(a, b, c, d) + { + assert_eq!(t_ab, 0.25f64); + assert_eq!(t_cd, 0.50f64); + assert_eq!( + intersection, + Coord { + x: 0.5f64, + y: 0.5f64 + } + ); + } else { + assert!(false) + } + } + + #[test] + fn test_line_segment_intersection_with_parameters_parallel() { + let a = Coord { x: 3f64, y: 4f64 }; + let b = Coord { x: 6f64, y: 8f64 }; + let c = Coord { x: 9f64, y: 9f64 }; + let d = Coord { x: 12f64, y: 13f64 }; + assert_eq!( + line_segment_intersection_with_parameters(a, b, c, d), + None + ) + } + #[test] + fn test_line_segment_intersection_with_parameters_colinear() { + let a = Coord { x: 1f64, y: 2f64 }; + let b = Coord { x: 2f64, y: 4f64 }; + let c = Coord { x: 3f64, y: 6f64 }; + let d = Coord { x: 5f64, y: 10f64 }; + assert_eq!( + line_segment_intersection_with_parameters(a, b, c, d), + None + ) + } + + #[test] + fn test_line_segment_intersection_with_relationships() { + let a = Coord { x: 1f64, y: 2f64 }; + let b = Coord { x: 2f64, y: 3f64 }; + let c = Coord { x: 0f64, y: 2f64 }; + let d = Coord { x: -2f64, y: 6f64 }; + + let expected_intersection_point = Coord { x: 1f64 / 3f64, y: 4f64 / 3f64 }; + + if let Some(LineIntersectionResultWithRelationships { + ab, + cd, + intersection, + }) = line_segment_intersection_with_relationships(a, b, c, d) + { + assert_eq!(ab, FalseIntersectionPoint(BeforeStart)); + assert_eq!(cd, FalseIntersectionPoint(BeforeStart)); + assert_relative_eq!(intersection, expected_intersection_point); + } else { + assert!(false); + } + + if let Some(LineIntersectionResultWithRelationships { + ab, + cd, + intersection, + }) = line_segment_intersection_with_relationships(b, a, c, d) + { + assert_eq!(ab, FalseIntersectionPoint(AfterEnd)); + assert_eq!(cd, FalseIntersectionPoint(BeforeStart)); + assert_relative_eq!(intersection, expected_intersection_point); + } else { + assert!(false); + } + + if let Some(LineIntersectionResultWithRelationships { + ab, + cd, + intersection, + }) = line_segment_intersection_with_relationships(a, b, d, c) + { + assert_eq!(ab, FalseIntersectionPoint(BeforeStart)); + assert_eq!(cd, FalseIntersectionPoint(AfterEnd)); + assert_relative_eq!(intersection, expected_intersection_point); + } else { + assert!(false); + } + } + + #[test] + fn test_line_segment_intersection_with_relationships_true() { + let a = Coord { x: 0f64, y: 1f64 }; + let b = Coord { x: 2f64, y: 3f64 }; + let c = Coord { x: 0f64, y: 2f64 }; + let d = Coord { x: -2f64, y: 6f64 }; + + let expected_intersection_point = Coord { x: 1f64 / 3f64, y: 4f64 / 3f64 }; + + if let Some(LineIntersectionResultWithRelationships { + ab, + cd, + intersection, + }) = line_segment_intersection_with_relationships(a, b, c, d) + { + assert_eq!(ab, TrueIntersectionPoint); + assert_eq!(cd, FalseIntersectionPoint(BeforeStart)); + assert_relative_eq!(intersection, expected_intersection_point); + } else { + assert!(false); + } + } +} diff --git a/geo/src/algorithm/offset_curve/line_measured.rs b/geo/src/algorithm/offset_curve/line_measured.rs new file mode 100644 index 000000000..c5d95c265 --- /dev/null +++ b/geo/src/algorithm/offset_curve/line_measured.rs @@ -0,0 +1,30 @@ + +use crate::{Line, CoordNum}; + +// Note: Previously I had a struct called "OffsetLineRaw" which turned out to +// be a line and it's length. This new struct [LineMeasured] is basically the +// same, but +// - has wider use cases +// - has a name that communicate's it's content better + +/// A struct containing a [Line] and it's precalculated length. +/// +/// TODO: I always Assume that it is faster to store calculated lengths than +/// simply recalculate them every time they are needed. This is likely the case +/// if the length is recalculated many times for the same Line. But in practice +/// there may be no benefit if it is only used once or twice in the same +/// function, where the compiler might have made the same optimization on our +/// behalf +/// +/// TODO: I would like to mark this type as immutable so that it can only be +/// instantiated using the full struct constructor and never mutated. Apparently +/// that isn't possible in rust without making all members private. This is sad +/// because we give up nice destructuring syntax if we make members private. For +/// the time being, I am leaving members public. I think ultimately this type +/// might go away in refactoring and integrating with other parts of the geo +/// crate. +#[derive(Clone, PartialEq, Eq, Debug)] +pub(super)struct LineMeasured where T:CoordNum { + pub line:Line, + pub length:T +} \ No newline at end of file diff --git a/geo/src/algorithm/offset_curve/mod.rs b/geo/src/algorithm/offset_curve/mod.rs new file mode 100644 index 000000000..b5662e925 --- /dev/null +++ b/geo/src/algorithm/offset_curve/mod.rs @@ -0,0 +1,9 @@ + +mod vector_extensions; +mod line_intersection; +mod offset_segments_iterator; +mod offset_line_raw; +mod line_measured; + +mod offset_curve_trait; +pub use offset_curve_trait::OffsetCurve; diff --git a/geo/src/algorithm/offset_curve/offset_curve_trait.rs b/geo/src/algorithm/offset_curve/offset_curve_trait.rs new file mode 100644 index 000000000..d7e87a188 --- /dev/null +++ b/geo/src/algorithm/offset_curve/offset_curve_trait.rs @@ -0,0 +1,324 @@ +// TODO: Should I be doing `use crate ::{...}` or `use geo_types::{...}` +use crate::{Coord, CoordFloat, Line, LineString, MultiLineString}; + +use super::line_intersection::{ + FalseIntersectionPointType::AfterEnd, + LineIntersectionResultWithRelationships, + LineSegmentIntersectionType::{FalseIntersectionPoint, TrueIntersectionPoint}, +}; + +use super::line_measured::LineMeasured; + +use super::vector_extensions::VectorExtensions; + +use super::offset_line_raw::offset_line_raw; + +use super::offset_segments_iterator::{LineStringOffsetSegmentPairs, OffsetSegmentsIteratorItem}; + + +/// The OffsetCurve trait is implemented for geometries where the edges of the +/// geometry can be offset perpendicular to the direction of the edges by some +/// positive or negative distance. For example, an offset [Line] will become a +/// [Line], and an offset [LineString] will become a [LineString]. +/// Geometry with no length ([geo_types::Point]) cannot be offset as it has no +/// directionality. +/// +/// > NOTE: The [OffsetCurve::offset_curve()] function is different to a `buffer` operation. +/// > A buffer (or inset / outset operation) would normally produce an enclosed +/// > shape; For example a [geo_types::Point] would become a circular +/// > [geo_types::Polygon], a [geo_types::Line] would become a capsule shaped +/// > [geo_types::Polygon]. + +pub trait OffsetCurve +where + T: CoordFloat, + Self: Sized, +{ + /// Offset the edges of the geometry by `distance`. + /// + /// In a coordinate system where positive is up and to the right; + /// when facing the direction of increasing coordinate index: + /// + /// - Positive `distance` will offset the edges of a geometry to the left + /// - Negative `distance` will offset the edges of a geometry to the right + /// + /// If you are using 'screen coordinates' where the y axis is often flipped + /// then the offset direction described above will be reversed. + /// + /// # Examples + /// + /// ``` + /// #use crate::{line_string, Coord}; + /// let input = line_string![ + /// Coord { x: 0f64, y: 0f64 }, + /// Coord { x: 0f64, y: 2f64 }, + /// Coord { x: 2f64, y: 2f64 }, + /// ]; + /// let output_expected = line_string![ + /// Coord { x: 1f64, y: 0f64 }, + /// Coord { x: 1f64, y: 1f64 }, + /// Coord { x: 2f64, y: 1f64 }, + /// ]; + /// let output_actual = input.offset_curve(-1f64).unwrap(); + /// assert_eq!(output_actual, output_expected); + /// ``` + fn offset_curve(&self, distance: T) -> Option; +} + +impl OffsetCurve for Line +where + T: CoordFloat, +{ + fn offset_curve(&self, distance: T) -> Option { + if distance == T::zero() { + // prevent unnecessary work; + Some(self.clone()) + // TODO: for typical use cases the offset would rarely be zero; + // This check may add unnecessary branching when there are a lot of + // Lines. It makes more sense to do this performance check for + // LineStrings... + } else { + let Line { start: a, end: b } = *self; + match offset_line_raw(a, b, distance) { + Some(LineMeasured { line, .. }) => Some(line), + _ => None, + } + } + } +} + +impl OffsetCurve for LineString +where + T: CoordFloat, +{ + fn offset_curve(&self, distance: T) -> Option { + // Loosely follows the algorithm described by + // [Xu-Zheng Liu, Jun-Hai Yong, Guo-Qin Zheng, Jia-Guang Sun. An offset + // algorithm for polyline curves. Computers in Industry, Elsevier, 2007, + // 15p. inria-00518005] + // (https://hal.inria.fr/inria-00518005/document) + + // Handle trivial cases; + // Note: Docs say LineString is valid "if it is either empty or contains + // two or more coordinates" + + // TODO: is `self.into_inner()` rather than `self.0` preferred? The + // contents of the tuple struct are public. + // Issue #816 seems to suggest that `self.0` is to be deprecated + match self.0.len() { + 0 => return Some(self.clone()), + 1 => return None, + 2 => { + return match Line::new(self.0[0], self.0[1]).offset_curve(distance) { + Some(line) => Some(line.into()), + None => None, + } + } + _ => (), + } + + // Prevent unnecessary work: + if T::is_zero(&distance) { + return Some(self.clone()); + } + + // TODO: Parameterize miter limit, and miter limit distance / factor. + let mitre_limit_factor = T::from(2.0).unwrap(); + let mitre_limit_distance = distance.abs() * mitre_limit_factor; + let mitre_limit_distance_squared = mitre_limit_distance * mitre_limit_distance; + + let mut offset_points = Vec::with_capacity(self.0.len()); + + for item in self.iter_offset_segment_pairs(distance) { + println!("{item:?}"); + if let OffsetSegmentsIteratorItem { + ab_offset: Some(LineMeasured { line, .. }), + first: true, + .. + } = item + { + offset_points.push(line.start) + }; + match item { + OffsetSegmentsIteratorItem { + ab_offset: + Some(LineMeasured { + line: Line { start: m, end: n }, + length: ab_len, + }), + bc_offset: + Some(LineMeasured { + line: Line { start: o, end: p }, + length: bc_len, + }), + i: + Some(LineIntersectionResultWithRelationships { + ab, + cd, + intersection, + }), + .. + } => match (ab, cd) { + (TrueIntersectionPoint, TrueIntersectionPoint) => { + // Inside elbow + // No mitre limit needed + offset_points.push(intersection) + } + (FalseIntersectionPoint(AfterEnd), FalseIntersectionPoint(_)) => { + // Outside elbow + // Check for Mitre Limit + let elbow_length_squared = (intersection - n).magnitude_squared(); + if elbow_length_squared > mitre_limit_distance_squared { + // Mitre Limited / Truncated Corner + let mn: Coord = n - m; + let op: Coord = p - o; + offset_points.push(n + mn / ab_len * mitre_limit_distance); + offset_points.push(o - op / bc_len * mitre_limit_distance); + } else { + // Sharp Corner + offset_points.push(intersection) + } + } + _ => { + // Inside pinched elbow + // (ie forearm curled back through bicep 🙃) + //println!("CASE 3 - bridge"); + offset_points.push(n); + offset_points.push(o); + } + }, + OffsetSegmentsIteratorItem { + ab_offset: + Some(LineMeasured { + line: Line { end: n, .. }, + .. + }), + i: None, + .. + } => { + // Collinear + // TODO: this is not an elegant way to handle colinear + // input: in some (all?) cases this produces a redundant + // colinear point in the output. It might be easier to + // eliminate this redundant point in a pre-processing step + // rather than try do it here. + offset_points.push(n) + } + _ => { + // Several ways to end up here... probably one of the + // segments could not be offset + return None; + } + } + if let OffsetSegmentsIteratorItem { + bc_offset: Some(LineMeasured { line, .. }), + last: true, + .. + } = item + { + offset_points.push(line.end) + }; + } + + Some(offset_points.into()) + } +} + +impl OffsetCurve for MultiLineString +where + T: CoordFloat, +{ + fn offset_curve(&self, distance: T) -> Option { + self.iter() + .map(|item| item.offset_curve(distance)) + .collect() + } +} + +#[cfg(test)] +mod test { + + use crate::{ + line_string, + Coord, + Line, + //LineString, + MultiLineString, + OffsetCurve, + }; + + #[test] + fn test_offset_line() { + let input = Line::new(Coord { x: 1f64, y: 1f64 }, Coord { x: 1f64, y: 2f64 }); + let output_actual = input.offset_curve(-1.0); + let output_expected = Some(Line::new( + Coord { x: 2f64, y: 1f64 }, + Coord { x: 2f64, y: 2f64 }, + )); + assert_eq!(output_actual, output_expected); + } + #[test] + fn test_offset_line_negative() { + let input = Line::new(Coord { x: 1f64, y: 1f64 }, Coord { x: 1f64, y: 2f64 }); + let output_actual = input.offset_curve(1.0); + let output_expected = Some(Line::new( + Coord { x: 0f64, y: 1f64 }, + Coord { x: 0f64, y: 2f64 }, + )); + assert_eq!(output_actual, output_expected); + } + + #[test] + fn test_offset_line_string() { + let input = line_string![ + Coord { x: 0f64, y: 0f64 }, + Coord { x: 0f64, y: 2f64 }, + Coord { x: 2f64, y: 2f64 }, + ]; + let output_actual = input.offset_curve(-1f64); + let output_expected = Some(line_string![ + Coord { x: 1f64, y: 0f64 }, + Coord { x: 1f64, y: 1f64 }, + Coord { x: 2f64, y: 1f64 }, + ]); + assert_eq!(output_actual, output_expected); + } + + #[test] + fn test_offset_line_string_invalid() { + let input = line_string![Coord { x: 0f64, y: 0f64 },]; + let output_actual = input.offset_curve(-1f64); + let output_expected = None; + assert_eq!(output_actual, output_expected); + } + + #[test] + fn test_offset_multi_line_string() { + let input = MultiLineString::new(vec![ + line_string![ + Coord { x: 0f64, y: 0f64 }, + Coord { x: 0f64, y: 2f64 }, + Coord { x: 2f64, y: 2f64 }, + ], + line_string![ + Coord { x: 0f64, y: 0f64 }, + Coord { x: 0f64, y: -2f64 }, + Coord { x: -2f64, y: -2f64 }, + ], + ]); + let output_actual = input.offset_curve(-1f64); + let output_expected = Some(MultiLineString::new(vec![ + line_string![ + Coord { x: 1f64, y: 0f64 }, + Coord { x: 1f64, y: 1f64 }, + Coord { x: 2f64, y: 1f64 }, + ], + line_string![ + Coord { x: -1f64, y: 0f64 }, + Coord { x: -1f64, y: -1f64 }, + Coord { x: -2f64, y: -1f64 }, + ], + ])); + assert_eq!(output_actual, output_expected); + } +} diff --git a/geo/src/algorithm/offset_curve/offset_line_raw.rs b/geo/src/algorithm/offset_curve/offset_line_raw.rs new file mode 100644 index 000000000..c693dc710 --- /dev/null +++ b/geo/src/algorithm/offset_curve/offset_line_raw.rs @@ -0,0 +1,84 @@ +use crate::{Coord, CoordFloat, Line}; +use super::{vector_extensions::VectorExtensions, line_measured::LineMeasured}; + + +/// The result of the [offset_line_raw()] function +// #[derive(Clone)] +// pub(super) struct OffsetLineRawResult where T:CoordFloat { +// pub a_offset:Coord, +// pub b_offset:Coord, +// pub ab_len:T, +// } + + +/// Offset a line defined by [Coord]s `a` and `b` by `distance`. +/// +/// In a coordinate system where positive is up and to the right; +/// Positive `distance` will offset the line to the left (when standing +/// at `a` and facing `b`) +/// +/// This could be implemented on [geo_types::Line]... +/// +/// There are 2 reasons +/// +/// 1. I am trying to localize my changes to the offset_curve module for now. +/// 2. I am trying to do is avoid repeated calculation of segment length. +/// This function has a special return type which also yields the length. +/// +/// TODO: In future it may be preferable to create new types called +/// `LineMeasured` and `LineStringMeasured` which store pre-computed length. +/// +/// - Confirm if significant performance benefit to using a bigger structs to +/// avoid recomputing the line segment length? +/// - I think there certainly might be in future parts of the algorithm which +/// need the length repeatedly) +/// +/// +pub(super) fn offset_line_raw( + a: Coord, + b: Coord, + distance: T, +) -> Option> +where + T: CoordFloat, +{ + let ab = b - a; + let length = ab.magnitude(); + if length == T::zero() { + return None; + } + let ab_offset = ab.left() / length * distance; + + Some(LineMeasured { + line: Line{ + start:a + ab_offset, + end: b + ab_offset + }, + length, + }) +} + + +// TODO: test + +#[cfg(test)] +mod test { + + use crate::{ + Coord, + Line, + offset_curve::{offset_line_raw::offset_line_raw, line_measured::LineMeasured}, + }; + + #[test] + fn test_offset_line_raw() { + let a = Coord { x: 0f64, y: 0f64 }; + let b = Coord { x: 0f64, y: 1f64 }; + let output_actual = offset_line_raw(a, b, 1f64); + let output_expected = Some(LineMeasured{ + line:Line { start: Coord { x: 1f64, y: 0f64 }, end: Coord { x: 1f64, y: 1f64 } }, + length:1f64, + }); + assert_eq!(output_actual, output_expected); + } +} \ No newline at end of file diff --git a/geo/src/algorithm/offset_curve/offset_segments_iterator.rs b/geo/src/algorithm/offset_curve/offset_segments_iterator.rs new file mode 100644 index 000000000..704a3cb5d --- /dev/null +++ b/geo/src/algorithm/offset_curve/offset_segments_iterator.rs @@ -0,0 +1,214 @@ +/// +/// My requirements are +/// +/// - Facilitate iterating over `Line`s in a LineString in a pairwise fashion +/// - Offset the `Line` inside the iterator +/// - Avoid repeatedly calculating length for each line +/// - Make iterator lazier (don't keep all offset `Line`s in memory) +/// - Iterator should provide +/// - the offset points +/// - the intersection point ([LineIntersectionResultWithRelationships]) +/// - the pre-calculated length of offset line segments (for miter limit +/// calculation) +/// - support wrapping over to the first segment at the end to simplify +/// closed shapes +/// +use crate::{Coord, CoordFloat, CoordNum, LineString}; + +use super::{ + line_intersection::{ + line_segment_intersection_with_relationships, LineIntersectionResultWithRelationships, + }, + line_measured::LineMeasured, + offset_line_raw::offset_line_raw, +}; + +/// Bring this into scope to imbue [LineString] with +/// [iter_offset_segment_pairs()] +pub(super) trait LineStringOffsetSegmentPairs +where + T: CoordFloat, +{ + /// Loop over the segments of a [LineString] in a pairwise fashion, + /// offsetting and intersecting them as we go. + /// + /// Returns an [OffsetSegmentsIterator] + fn iter_offset_segment_pairs(&self, distance: T) -> OffsetSegmentsIterator; +} + +pub(super) struct OffsetSegmentsIterator<'a, T> +where + T: CoordFloat, +{ + line_string: &'a LineString, + distance: T, + previous_offset_segment: Option>, + index: usize, +} + +impl LineStringOffsetSegmentPairs for LineString +where + T: CoordFloat, +{ + fn iter_offset_segment_pairs(&self, distance: T) -> OffsetSegmentsIterator + where + T: CoordNum, + { + if self.0.len() < 3 { + // LineString is not long enough, therefore return an iterator that + // will return None as first result + OffsetSegmentsIterator { + line_string: self, + distance, + previous_offset_segment: None, + index: 0, + } + } else { + // TODO: Length check above prevents panic; use + // unsafe get_unchecked for performance? + let a = self.0[0]; + let b = self.0[1]; + OffsetSegmentsIterator { + line_string: self, + distance, + previous_offset_segment: offset_line_raw(a, b, distance), + index: 0, + } + } + } +} + + +/// - [LineString] `a---b---c` is offset to form +/// - [LineMeasured] `ab_offset` (`a'---b'`) and +/// - [LineMeasured] `bc_offset` (`b'---c'`) +/// - [LineIntersectionResultWithRelationships] `i` is the intersection point. +/// +/// ```text +/// a +/// a' \ +/// \ \ +/// \ b---------c +/// b' +/// +/// i b'--------c' +/// ``` +#[derive(Clone, Debug)] +pub(super) struct OffsetSegmentsIteratorItem +where + T: CoordNum, +{ + /// This is true for the first result + pub first: bool, + + // this is true for the last result + pub last: bool, + + // TODO: seems a,b,c are unused... + pub a: Coord, + pub b: Coord, + pub c: Coord, + + pub ab_offset: Option>, + pub bc_offset: Option>, + + /// Intersection [Coord] between segments `mn` and `op` + pub i: Option>, +} + +impl<'a, T> Iterator for OffsetSegmentsIterator<'a, T> +where + T: CoordFloat, +{ + /// Option since each step of the iteration may fail. + type Item = OffsetSegmentsIteratorItem; + + /// Return type is confusing; `Option>>` + /// + /// TODO: Revise + /// + /// The outer Option is required by the Iterator trait, and indicates if + /// iteration is finished, (When this iterator is used via `.map()` or + /// similar the user does not see the outer Option.) + /// The inner Option indicates if the result of each iteration is valid. + /// Returning None will halt iteration, returning Some(None) will not, + /// but the user should stop iterating. + /// + fn next(&mut self) -> Option { + if self.index + 3 > self.line_string.0.len() { + // Iteration is complete + return None; + } else { + // TODO: Length check above prevents panic; use + // unsafe get_unchecked for performance? + let a = self.line_string[self.index]; + let b = self.line_string[self.index + 1]; + let c = self.line_string[self.index + 2]; + + self.index += 1; + + // Fetch previous offset segment + let ab_offset = self.previous_offset_segment.clone(); + + // Compute next offset segment + self.previous_offset_segment = offset_line_raw(b, c, self.distance); + + Some(OffsetSegmentsIteratorItem { + first: self.index == 1, + last: self.index + 3 > self.line_string.0.len(), + a, + b, + c, + i: match (&ab_offset, &self.previous_offset_segment) { + (Some(ab_offset), Some(bc_offset)) => { + line_segment_intersection_with_relationships( + ab_offset.line.start, + ab_offset.line.end, + bc_offset.line.start, + bc_offset.line.end, + ) + } + _ => None, + }, + ab_offset, + bc_offset: self.previous_offset_segment.clone(), + }) + } + } +} + +#[cfg(test)] +mod test { + use super::{LineStringOffsetSegmentPairs, OffsetSegmentsIteratorItem}; + use crate::{ + line_string, + offset_curve::{ + line_intersection::LineIntersectionResultWithRelationships, line_measured::LineMeasured, + }, + Coord, + }; + + #[test] + fn test_iterator() { + let input = line_string![ + Coord { x: 1f64, y: 0f64 }, + Coord { x: 1f64, y: 1f64 }, + Coord { x: 2f64, y: 1f64 }, + ]; + + // TODO: this test is a bit useless after recent changes + let result: Option> = input + .iter_offset_segment_pairs(1f64) + .map(|item| match item { + OffsetSegmentsIteratorItem { + ab_offset: Some(LineMeasured { .. }), + bc_offset: Some(LineMeasured { .. }), + i: Some(LineIntersectionResultWithRelationships { .. }), + .. + } => Some(()), + _ => None, + }) + .collect(); + assert!(result.unwrap().len() == 1); + } +} diff --git a/geo/src/algorithm/offset_curve/vector_extensions.rs b/geo/src/algorithm/offset_curve/vector_extensions.rs new file mode 100644 index 000000000..f50b9bbcd --- /dev/null +++ b/geo/src/algorithm/offset_curve/vector_extensions.rs @@ -0,0 +1,284 @@ +use crate::{Coord, CoordFloat}; + +/// Extends the `Coord` struct with some common vector operations; +/// +/// - [VectorExtensions::cross_product_2d], +/// - [VectorExtensions::magnitude], +/// - [VectorExtensions::magnitude_squared], +/// - [VectorExtensions::dot_product], +/// - [VectorExtensions::left], +/// - [VectorExtensions::right], +/// +/// TODO: I implemented these functions here because I am trying to localize +/// my changes to the [crate::algorithm::offset_curve] module at the moment. +/// +/// The [crate::algorithm::Kernel] trait has some functions which overlap with +/// this trait. I have realized I am re-inventing the wheel here. + +pub(super) trait VectorExtensions +where + Self: Copy, +{ + type NumericType; + /// The 2D cross product is the signed magnitude of the 3D "Cross Product" + /// assuming z ordinates are zero. + /// + /// ## Other names for this function: + /// + /// - In exterior algebra, it is called the wedge product. + /// - If the inputs are packed into a 2x2 matrix, this is the determinant. + /// + /// ## Other appearances in this library + /// + /// 1. [geo_types::Point::cross_prod()] is already defined on + /// [geo_types::Point]... but that it seems to be some other + /// operation on 3 points?? + /// + /// 2. Note: The [geo_types::Line] struct also has a + /// [geo_types::Line::determinant()] function which is the same as + /// `cross_product_2d(line.start, line.end)` + /// + /// 3. The [crate::algorithm::Kernel::orient2d()] trait default + /// implementation uses cross product to compute orientation. It returns + /// an enum, not the numeric value which is needed for line segment + /// intersection. + /// + /// + /// ## Properties + /// + /// - The absolute value of the cross product is the area of the + /// parallelogram formed by the operands + /// - The sign of the output is reversed if the operands are reversed + /// - The sign can be used to check if the operands are clockwise / + /// anti-clockwise orientation with respect to the origin; + /// or phrased differently: + /// "is b to the left of the line between the origin and a"? + /// - If the operands are colinear with the origin, the magnitude is zero + /// + /// + /// ## Derivation + /// + /// From basis vectors `i`,`j`,`k` and the axioms on wikipedia + /// [Cross product](https://en.wikipedia.org/wiki/Cross_product#Computing); + /// + /// ```text + /// i×j = k + /// j×k = i + /// k×i = j + /// + /// j×i = -k + /// k×j = -i + /// i×k = -j + /// + /// i×i = j×j = k×k = 0 + /// ``` + /// + /// We can define the 2D cross product as the magnitude of the 3D cross + /// product as follows + /// + /// ```text + /// |a × b| = |(a_x·i + a_y·j + 0·k) × (b_x·i + b_y·j + 0·k)| + /// = |a_x·b_x·(i×i) + a_x·b_y·(i×j) + a_y·b_x·(j×i) + a_y·b_y·(j×j)| + /// = |a_x·b_x·( 0 ) + a_x·b_y·( k ) + a_y·b_x·(-k ) + a_y·b_y·( 0 )| + /// = | (a_x·b_y - a_y·b_x)·k | + /// = a_x·b_y - a_y·b_x + /// ``` + fn cross_product_2d(self, other: Rhs) -> Self::NumericType; + + /// The inner product of the coordinate components + /// + /// ```ignore + /// self.x*other.x + self.y*other.y + /// ``` + fn dot_product(self, other: Rhs) -> Self::NumericType; + + /// The euclidean distance between this coordinate and the origin + /// + /// ``` + /// (self.x*self.x + self.y*self.y).sqrt() + /// ``` + fn magnitude(self) -> Self::NumericType; + + /// The squared distance between this coordinate and the origin. + /// (Avoids the square root calculation when it is not needed) + /// + /// ``` + /// self.x*self.x + self.y*self.y + /// ``` + fn magnitude_squared(self) -> Self::NumericType; + + /// In a coordinate system where positive is up and to the right; + /// Rotate this coordinate around the origin in the xy plane 90 degrees + /// anti-clockwise (Consistent with [crate::algorithm::rotate::Rotate]). + fn left(self) -> Self; + + /// In a coordinate system where positive is up and to the right; + /// Rotate this coordinate around the origin in the xy plane 90 degrees + /// clockwise (Consistent with [crate::algorithm::rotate::Rotate]). + fn right(self) -> Self; +} + +impl VectorExtensions for Coord +where + T: CoordFloat, +{ + type NumericType = T; + + fn cross_product_2d(self, right: Coord) -> Self::NumericType { + self.x * right.y - self.y * right.x + } + + fn dot_product(self, other: Self) -> Self::NumericType { + self.x * other.x + self.y * other.y + } + + fn magnitude(self) -> Self::NumericType { + (self.x * self.x + self.y * self.y).sqrt() + } + + fn magnitude_squared(self) -> Self::NumericType { + self.x * self.x + self.y * self.y + } + + fn left(self) -> Self { + Self { + x: -self.y, + y: self.x, + } + } + + fn right(self) -> Self { + Self { + x: self.y, + y: -self.x, + } + } +} + +#[cfg(test)] +mod test { + // crate dependencies + use crate::Coord; + + // private imports + use super::VectorExtensions; + + #[test] + fn test_cross_product() { + // perpendicular unit length + let a = Coord { x: 1f64, y: 0f64 }; + let b = Coord { x: 0f64, y: 1f64 }; + + // expect the area of parallelogram + assert_eq!(a.cross_product_2d(b), 1f64); + // expect swapping will result in negative + assert_eq!(b.cross_product_2d(a), -1f64); + + // Add skew; expect results should be the same + let a = Coord { x: 1f64, y: 0f64 }; + let b = Coord { x: 1f64, y: 1f64 }; + + // expect the area of parallelogram + assert_eq!(a.cross_product_2d(b), 1f64); + // expect swapping will result in negative + assert_eq!(b.cross_product_2d(a), -1f64); + + // Make Colinear; expect zero + let a = Coord { x: 2f64, y: 2f64 }; + let b = Coord { x: 1f64, y: 1f64 }; + assert_eq!(a.cross_product_2d(b), 0f64); + } + + #[test] + fn test_dot_product() { + // perpendicular unit length + let a = Coord { x: 1f64, y: 0f64 }; + let b = Coord { x: 0f64, y: 1f64 }; + // expect zero for perpendicular + assert_eq!(a.dot_product(b), 0f64); + + // Parallel, same direction + let a = Coord { x: 1f64, y: 0f64 }; + let b = Coord { x: 2f64, y: 0f64 }; + // expect +ive product of magnitudes + assert_eq!(a.dot_product(b), 2f64); + // expect swapping will have same result + assert_eq!(b.dot_product(a), 2f64); + + // Parallel, opposite direction + let a = Coord { x: 3f64, y: 4f64 }; + let b = Coord { x: -3f64, y: -4f64 }; + // expect -ive product of magnitudes + assert_eq!(a.dot_product(b), -25f64); + // expect swapping will have same result + assert_eq!(b.dot_product(a), -25f64); + } + + #[test] + fn test_magnitude() { + let a = Coord { x: 1f64, y: 0f64 }; + assert_eq!(a.magnitude(), 1f64); + + let a = Coord { x: 0f64, y: 0f64 }; + assert_eq!(a.magnitude(), 0f64); + + let a = Coord { x: -3f64, y: 4f64 }; + assert_eq!(a.magnitude(), 5f64); + } + + #[test] + fn test_magnitude_squared() { + let a = Coord { x: 1f64, y: 0f64 }; + assert_eq!(a.magnitude_squared(), 1f64); + + let a = Coord { x: 0f64, y: 0f64 }; + assert_eq!(a.magnitude_squared(), 0f64); + + let a = Coord { x: -3f64, y: 4f64 }; + assert_eq!(a.magnitude_squared(), 25f64); + } + + #[test] + fn test_left_right() { + let a = Coord { x: 1f64, y: 0f64 }; + let a_left = Coord { x: 0f64, y: 1f64 }; + let a_right = Coord { x: 0f64, y: -1f64 }; + + assert_eq!(a.left(), a_left); + assert_eq!(a.right(), a_right); + assert_eq!(a.left(), -a.right()); + } + + #[test] + fn test_left_right_match_rotate() { + use crate::algorithm::rotate::Rotate; + use crate::Point; + // The aim of this test is to confirm that wording in documentation is + // consistent. + + // when the user is in a coordinate system where the y axis is flipped + // (eg screen coordinates in a HTML canvas), then rotation directions + // will be different to those described in the documentation. + + // The documentation for the Rotate trait says: 'Positive angles are + // counter-clockwise, and negative angles are clockwise rotations' + + let counter_clockwise_rotation_degrees = 90.0; + let clockwise_rotation_degrees = -counter_clockwise_rotation_degrees; + + let a: Point = Coord { x: 1.0, y: 0.0 }.into(); + let origin:Point = Coord::::zero().into(); + + // left is anti-clockwise + assert_relative_eq!( + Point::from(a.0.left()), + a.rotate_around_point(counter_clockwise_rotation_degrees, origin), + ); + // right is clockwise + assert_relative_eq!( + Point::from(a.0.right()), + a.rotate_around_point(clockwise_rotation_degrees, origin), + ); + + } +} diff --git a/rfcs/2022-11-11-offset.md b/rfcs/2022-11-11-offset.md new file mode 100644 index 000000000..3d8345926 --- /dev/null +++ b/rfcs/2022-11-11-offset.md @@ -0,0 +1,182 @@ +# Offset Curve + +- [1. Introduction](#1-introduction) +- [2. OffsetCurve Trait](#2-offsetcurve-trait) +- [3. Trait Implementations](#3-trait-implementations) +- [4. To Do / Limitations](#4-to-do--limitations) +- [5. Algorithm](#5-algorithm) + - [5.1. References](#51-references) + - [5.2. Definitions (For the psudo-code in this readme only)](#52-definitions-for-the-psudo-code-in-this-readme-only) + - [5.3. Algorithm Part 0 - Pre-Treatment](#53-algorithm-part-0---pre-treatment) + - [5.4. Algorithm Part 0.1 - Segment Offset](#54-algorithm-part-01---segment-offset) + - [5.5. Algorithm Part 1 - Line Extension](#55-algorithm-part-1---line-extension) + - [5.6. Algorithm Part 4.1 - Dual Clipping - **(TODO: not yet implemented)**](#56-algorithm-part-41---dual-clipping---todo-not-yet-implemented) + - [5.7. Algorithm Part 4.1.2 - Cookie Cutter - **(TODO: not yet implemented)**](#57-algorithm-part-412---cookie-cutter---todo-not-yet-implemented) + - [5.8. Algorithm Part 4.1.3 - Proximity Clipping **(TODO: not yet implemented)**](#58-algorithm-part-413---proximity-clipping-todo-not-yet-implemented) + +## 1. Introduction + +Proposal to add an Offset Curve operation which allows linear geometry to be offset perpendicular to its direction based on the paper outlined below. + +Note that the proposed `offset_curve` is different to a `buffer` operation. + +## 2. OffsetCurve Trait + +Create a Trait called `OffsetCurve` in the `algorithms` module. + +## 3. Trait Implementations + +- [X] `Line` +- [X] `LineString` +- [X] `MultiLineString` +- [ ] `Triangle` +- [ ] `Rectangle` +- [ ] `Polygon` +- [ ] `MultiPolygon` +- [ ] ~~`Point`~~
+ - [ ] ~~`MultiPoint`~~ + - [ ] ~~`Geometry`~~ + - [ ] ~~`GeometryCollection`~~ + +> Note: +> Offset for `Point` Can't be defined without confusion. +> The user may expect the result to be one of the following: +> +> 1. Return the input unchanged, +> 1. a translation (but in what direction?) +> 1. a circle (like a buffer, but as previously explained, +> the `offset` operation is different to the `buffer` operation). +> 1. an empty `LineString` (wierd, yes, but this is actually what GEOS does 🤔) + + +For coordinate types + +- [X] `where T:CoordFloat` +- [ ] `where T:CoordNum` (??? seems tricky) + +## 4. To Do / Limitations + +Some may be removed during development + +- [X] Change name from `Offset` to `OffsetCurve` to match terminology in `GEOS` + and `jts` +- [ ] Make proper use of SimpleKernel / RobustKernel +- [X] Prevent execution when the specified offset distance is zero +- [X] Mitre-limit is implemented so that a LineString which doubles-back on + itself will not produce an elbow at infinity + - [ ] Make Mitre Limit configurable? +- [ ] Option to clip output such that non-adjacent line segments in the output + are not self-intersecting. +- [ ] Handle closed shapes + +## 5. Algorithm + +### 5.1. References + +Loosely follows the algorithm described by +[Xu-Zheng Liu, Jun-Hai Yong, Guo-Qin Zheng, Jia-Guang Sun. An offset algorithm for polyline curves. Computers in Industry, Elsevier, 2007, 15p. inria-00518005](https://hal.inria.fr/inria-00518005/document) + +### 5.2. Definitions (For the psudo-code in this readme only) + +Type definitions +```python +from typing import List, Tuple +from nicks_line_tools.Vector2 import Vector2 + +# define type aliases for convenience: +LineString = List[Vector2] +MultiLineString = List[LineString] +LineSegment = Tuple[Vector2, Vector2] +LineSegmentList = List[Tuple[Vector2, Vector2]] +Parameter = float + +# declare type of variables used: +input_linestring: LineString +offset_ls: LineString +offset_segments: LineSegmentList +raw_offset_ls: LineString +``` + +Function Type Definitions (pseudocode) + +```python +intersect = (tool: LineString, target: LineString) -> (point_of_intersection: Optional[Vector2], distance_along_target: List[Parameter]) +project = (tool: Vector2, target: LineString) -> (nearest_point_on_target_to_tool: Vector2, distance_along_target: Parameter) +interpolate = (distance_along_target: Parameter, target: LineString) -> (point_on_target: Vector2) +``` + +### 5.3. Algorithm Part 0 - Pre-Treatment + +1. Pretreatment steps from the paper are not implemented... these mostly deal with arcs and malformed input geometry + + +### 5.4. Algorithm Part 0.1 - Segment Offset + +1. Create an empty `LineSegmentList` called `offset_segments` +1. For each `LineSegment` of `input_linestring` + 1. Take each segment `(a,b)` of `input_linestring` and compute the vector from the start to the end of the segment
+ `ab = b - a` + 1. rotate this vector by -90 degrees to obtain the 'left normal'
+ `ab_left = Vector2(-ab.y, ab.x)` + 1. normalise `ab_left` by dividing each component by the magnitude of `ab_left`. + 1. multiply the vector by the scalar `d` to obtain the `segment_offset_vector` + 1. add the `segment_offset_vector` to `a` and `b` to get `offset_a` and `offset_b` + 1. append `(offset_a, offset_b)` to `offset_segments` + + +### 5.5. Algorithm Part 1 - Line Extension + +4. Create an empty `LineString` called `raw_offset_ls` +1. Append `offset_segments[0][0]` to `raw_offset_ls` +1. For each pair of consecutive segments `(a,b),(c,d)` in `offset_segments` + 1. If `(a,b)` is co-linear with `(c,d)` then append `b` to raw_offset_ls, and go to the next pair of segments. + 1. Otherwise, find the intersection point `p` of the infinite lines that are collinear with `(a,b)` and `(c,d)`;
+ Pay attention to the location of `p` relative to each of the segments: + 1. if `p` is within the segment it is a *True Intersection Point* or **TIP** + 1. If `p` is outside the segment it is called a *False Intersection Point* or **FIP**.
+ **FIP**s are further classified as follows: + - **Positive FIP** or **PFIP** if `p` is after the end of a segment, or + - **Negative FIP** or **NFIP** if `p` is before the start of a segment. + 1. If `p` is a **TIP** for both `(a,b)` and `(c,d)` append `p` to `raw_offset_ls` + 1. If `p` is a **FIP** for both `(a,b)` and `(c,d)` + 1. If `p` is a **PFIP** for `(a,b)` then append `p`to `raw_offset_ls` + 1. Otherwise, append `b` then `c` to `raw_offset_ls` + 1. Otherwise, append `b` then `c` to `raw_offset_ls` +1. Remove zero length segments in `raw_offset_ls` + +### 5.6. Algorithm Part 4.1 - Dual Clipping - **(TODO: not yet implemented)** + +8. Find `raw_offset_ls_twin` by repeating Algorithms 0.1 and 1 but offset the `input_linestring` in the opposite direction (`-d`) +1. Find `intersection_points` between + 1. `raw_offset_ls` and `raw_offset_ls` + 1. `raw_offset_ls` and `raw_offset_ls_twin` + +1. Find `split_offset_mls` by splitting `raw_offset_ls` at each point in `intersection_points` +1. If `intersection_points` was empty, then add `raw_offset_ls` to `split_offset_mls` and skip to Algorithm 4.2. +1. Delete each `LineString` in `split_offset_mls` if it intersects the `input_linestring`
+ unless the intersection is with the first or last `LineSegment` of `input_linestring` + 1. If we find such an intersection point that *is* on the first or last `LineSegment` of `input_linestring`
+ then add the intersection point to a list called `cut_targets` + +### 5.7. Algorithm Part 4.1.2 - Cookie Cutter - **(TODO: not yet implemented)** + +13. For each point `p` in `cut_targets` + 1. construct a circle of diameter `d` with its center at `p` + 1. delete all parts of any linestring in `split_offset_mls` which falls within this circle +1. Empty the `cut_targets` list + +### 5.8. Algorithm Part 4.1.3 - Proximity Clipping **(TODO: not yet implemented)** + +17. For each linestring `item_ls` in `split_offset_mls` + 1. For each segment `(a,b)` in `item_ls` + 1. For each segment `(u,v)` of `input_linestring` + - Find `a_proj` and `b_proj` by projecting `a` and `b` onto segment `(u,v)` + - Adjust the projected points such that `a_proj` and `b_proj` lie **at** or **between** `u` and `v` + - Find `a_dist` by computing `magnitude(a_proj - a)` + - Find `b_dist` by computing `magnitude(b_proj - b)` + - if either `a_dist` or `b_dist` is less than the absolute value of `d` then + - if `a_dist > b_dist`add `a_proj` to `cut_targets` + - Otherwise, add `b_proj` to `cut_targets` +1. Repeat Algorithm 4.1.2 +1. Remove zero length segments and empty linestrings etc +1. Join remaining linestrings that are touching to form new linestring(s) \ No newline at end of file