#[cfg(test)]
mod tests;
use crate::{LinestringError, linestring_2d::Line2};
use rustc_hash::FxHashSet;
use smallvec::{SmallVec, smallvec};
use std::{cmp::Ordering, collections::BTreeMap, convert::identity, fmt, fmt::Debug};
use vector_traits::{approx::*, num_traits::Float, prelude::*};
#[derive(Clone, Copy)]
struct SiteEventKey<T: GenericVector2> {
pub pos: T,
pub index: Option<u32>,
}
impl<T: GenericVector2> Debug for SiteEventKey<T> {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
f.debug_tuple("")
.field(&self.pos.x())
.field(&self.pos.y())
.field(&self.index)
.finish()
}
}
impl<T: GenericVector2> Ord for SiteEventKey<T> {
#[inline(always)]
fn cmp(&self, other: &Self) -> Ordering {
self.partial_cmp(other).unwrap()
}
}
#[allow(clippy::non_canonical_partial_ord_impl)]
impl<T: GenericVector2> PartialOrd for SiteEventKey<T>
where
T::Scalar: UlpsEq,
{
#[inline(always)]
fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
if ulps_eq!(&self.pos.y(), &other.pos.y()) {
if ulps_eq!(&self.pos.x(), &other.pos.x()) {
Some(Ordering::Equal)
} else {
self.pos.x().partial_cmp(&other.pos.x())
}
} else {
self.pos.y().partial_cmp(&other.pos.y())
}
}
}
impl<T: GenericVector2> Eq for SiteEventKey<T> {}
impl<T: GenericVector2> PartialEq for SiteEventKey<T> {
#[inline(always)]
fn eq(&self, other: &Self) -> bool {
self.pos.is_ulps_eq(
other.pos,
T::Scalar::default_epsilon(),
T::Scalar::default_max_ulps(),
)
}
}
struct MinMax<T: GenericVector2> {
best_left: Option<T::Scalar>,
slope: MinMaxSlope<T>,
best_right: Option<T::Scalar>,
}
impl<T: GenericVector2> MinMax<T>
where
T::Scalar: UlpsEq,
{
fn new() -> Self {
Self {
best_left: None,
best_right: None,
slope: MinMaxSlope::new(),
}
}
fn update(
&mut self,
pivot_x: T::Scalar,
candidate_x: T::Scalar,
candidate_slope: T::Scalar,
candidate_index: u32,
) {
if candidate_x < pivot_x {
if let Some(current_min) = self.best_left {
if ulps_eq!(¤t_min, &candidate_x) {
self.slope
.update_left(false, candidate_slope, candidate_index);
} else if current_min < candidate_x {
self.slope
.update_left(true, candidate_slope, candidate_index);
self.best_left = Some(candidate_x);
} else {
}
} else {
self.best_left = Some(candidate_x);
self.slope
.update_left(false, candidate_slope, candidate_index);
}
} else if candidate_x > pivot_x {
if let Some(current_max) = self.best_right {
if ulps_eq!(¤t_max, &candidate_x) {
self.slope
.update_right(false, candidate_slope, candidate_index);
} else if current_max > candidate_x {
self.slope
.update_right(true, candidate_slope, candidate_index);
self.best_right = Some(candidate_x);
} else {
}
} else {
self.best_right = Some(candidate_x);
self.slope
.update_right(false, candidate_slope, candidate_index);
}
}
}
fn clear(&mut self) {
self.best_left = None;
self.best_right = None;
self.slope.clear();
}
}
struct MinMaxSlope<T: GenericVector2> {
best_left: Option<T::Scalar>, candidates_left: Vec<u32>,
best_right: Option<T::Scalar>, candidates_right: Vec<u32>,
}
impl<T: GenericVector2> MinMaxSlope<T>
where
T::Scalar: UlpsEq,
{
fn new() -> Self {
Self {
best_left: None,
candidates_left: Vec::<_>::new(),
best_right: None,
candidates_right: Vec::<_>::new(),
}
}
fn update_both(&mut self, vertices: &[T], candidate_index: u32, lines: &[(u32, u32)]) {
let line = lines[candidate_index as usize];
let line_start = vertices[line.0 as usize];
let line_end = vertices[line.1 as usize];
let candidate_slope = if ulps_eq!(&line_end.y(), &line_start.y()) {
T::Scalar::INFINITY
} else {
(line_end.x() - line_start.x()) / (line_end.y() - line_start.y())
};
self.update_left(false, candidate_slope, candidate_index);
self.update_right(false, candidate_slope, candidate_index);
}
fn update_left(&mut self, clear: bool, candidate_slope: T::Scalar, candidate_index: u32) {
if clear {
self.candidates_left.clear();
self.best_left = None;
}
if let Some(current_slope) = self.best_left {
if ulps_eq!(¤t_slope, &candidate_slope) {
self.candidates_left.push(candidate_index);
} else if candidate_slope < current_slope {
self.candidates_left.clear();
self.candidates_left.push(candidate_index);
self.best_left = Some(candidate_slope);
} else {
}
} else {
self.best_left = Some(candidate_slope);
self.candidates_left.push(candidate_index);
}
}
fn update_right(&mut self, clear: bool, candidate_slope: T::Scalar, candidate_index: u32) {
if clear {
self.candidates_right.clear();
self.best_right = None;
}
if let Some(current_slope) = self.best_right {
if ulps_eq!(¤t_slope, &candidate_slope) {
self.candidates_right.push(candidate_index);
} else if candidate_slope > current_slope {
self.candidates_right.clear();
self.candidates_right.push(candidate_index);
self.best_right = Some(candidate_slope);
} else {
}
} else {
self.best_right = Some(candidate_slope);
self.candidates_right.push(candidate_index);
}
}
fn clear(&mut self) {
self.best_left = None;
self.candidates_left.clear();
self.best_right = None;
self.candidates_right.clear();
}
}
struct SiteEvent {
drop: Option<Vec<u32>>,
add: Option<Vec<u32>>,
intersection: Option<SmallVec<[u32; 1]>>,
}
impl SiteEvent {
pub(crate) fn with_intersection(i: &[u32]) -> Self {
Self {
drop: None,
add: None,
intersection: Some(i.into()),
}
}
pub(crate) fn with_drop(l: &[u32]) -> Self {
Self {
drop: Some(l.to_vec()),
add: None,
intersection: None,
}
}
pub(crate) fn with_add(l: &[u32]) -> Self {
Self {
drop: None,
add: Some(l.to_vec()),
intersection: None,
}
}
#[allow(dead_code)]
pub fn get_intersections(&self) -> &Option<SmallVec<[u32; 1]>> {
&self.intersection
}
}
fn sweepline_intersection<T: GenericVector2>(
vertices: &[T],
sweepline: T,
other: (u32, u32),
) -> Option<(T::Scalar, T::Scalar)> {
let other_start = vertices[other.0 as usize];
let other_end = vertices[other.1 as usize];
let y1 = other_start.y();
let y2 = other_end.y();
let x1 = other_start.x();
let x2 = other_end.x();
if ulps_eq!(&y1, &y2) {
return (sweepline.x() < x2).then_some((x2, T::Scalar::ZERO));
}
if ulps_eq!(&x1, &x2) {
return Some((x1, T::Scalar::INFINITY));
}
let slope = (y2 - y1) / (x2 - x1);
let d = y1 - slope * x1;
Some(((sweepline.y() - d) / slope, slope))
}
pub struct IntersectionTester<T: GenericVector2> {
sweepline_pos: T,
sweepline_pos_index: Option<u32>,
stop_at_first_intersection: bool,
pub ignore_end_point_intersections: bool,
site_events: Option<BTreeMap<SiteEventKey<T>, SiteEvent>>,
lines: Vec<(u32, u32)>,
vertices: Vec<T>,
}
impl<T: GenericVector2> IntersectionTester<T> {
pub fn new(vertices: Vec<T>) -> Self {
Self {
vertices,
sweepline_pos: T::new_2d(-T::Scalar::max_value(), -T::Scalar::max_value()),
sweepline_pos_index: None,
stop_at_first_intersection: false,
ignore_end_point_intersections: false,
site_events: Some(BTreeMap::new()),
lines: Vec::<(u32, u32)>::new(),
}
}
#[inline]
pub fn get_sweepline_pos(&self) -> (T, Option<u32>) {
(self.sweepline_pos, self.sweepline_pos_index)
}
#[inline]
pub fn get_lines(&self) -> &[(u32, u32)] {
&self.lines
}
pub fn with_stop_at_first_intersection(mut self, value: bool) -> Result<Self, LinestringError> {
self.stop_at_first_intersection = value;
Ok(self)
}
pub fn with_ignore_end_point_intersections(
mut self,
value: bool,
) -> Result<Self, LinestringError> {
self.ignore_end_point_intersections = value;
Ok(self)
}
pub fn with_edges<'a, IT>(mut self, input_iter: IT) -> Result<Self, LinestringError>
where
IT: IntoIterator<Item = &'a (u32, u32)>,
{
let mut site_events = self.site_events.take().ok_or_else(|| {
LinestringError::InternalError("with_lines() could not take 'site_events'".to_string())
})?;
for (index, aline) in input_iter.into_iter().enumerate() {
let index = index as u32;
let mut start_index = aline.0;
let mut start = self.vertices[start_index as usize];
let mut end_index = aline.1;
let mut end = self.vertices[end_index as usize];
if !(start.x().is_finite()
&& start.y().is_finite()
&& end.x().is_finite()
&& end.y().is_finite())
{
return Err(LinestringError::InvalidData(
"Some of the points are infinite".to_string(),
));
}
if !(SiteEventKey {
pos: start,
index: Some(start_index),
})
.lt(&(SiteEventKey {
pos: end,
index: Some(end_index),
})) {
std::mem::swap(&mut start_index, &mut end_index);
std::mem::swap(&mut start, &mut end);
};
self.lines.push(*aline);
let key_start = SiteEventKey {
pos: start,
index: Some(start_index),
};
let key_end = SiteEventKey {
pos: end,
index: Some(end_index),
};
if let Some(event) = site_events.get_mut(&key_start) {
let mut lower = event.add.take().map_or(Vec::<u32>::default(), identity);
lower.push(index);
event.add = Some(lower);
} else {
let event = SiteEvent::with_add(&[index]);
let _ = site_events.insert(key_start, event);
}
if let Some(event) = site_events.get_mut(&key_end) {
let mut upper = event.drop.take().map_or(Vec::<u32>::default(), identity);
upper.push(index);
event.drop = Some(upper);
} else {
let event = SiteEvent::with_drop(&[index]);
let _ = site_events.insert(key_end, event);
}
}
self.site_events = Some(site_events);
Ok(self)
}
#[inline]
fn add_intersection_event(
&self,
site_events: &mut BTreeMap<SiteEventKey<T>, SiteEvent>,
position: &SiteEventKey<T>,
intersecting_lines: &[u32],
) {
if let Some(event) = site_events.get_mut(position) {
let mut intersections_added = 0;
for new_intersection in intersecting_lines.iter() {
let i_line = self.lines[*new_intersection as usize];
if position.pos.is_ulps_eq(
self.vertices[i_line.0 as usize],
T::Scalar::default_epsilon(),
T::Scalar::default_max_ulps(),
) || position.pos.is_ulps_eq(
self.vertices[i_line.1 as usize],
T::Scalar::default_epsilon(),
T::Scalar::default_max_ulps(),
) {
continue;
}
if let Some(prev_intersections) = &mut event.intersection {
prev_intersections.push(*new_intersection);
} else {
let new_vec = smallvec![*new_intersection];
event.intersection = Some(new_vec);
}
intersections_added += 1;
}
if intersections_added > 0 {
let mut intersections = event.intersection.take().unwrap();
intersections.sort_unstable();
intersections.dedup();
event.intersection = Some(intersections);
}
} else {
let event = SiteEvent::with_intersection(intersecting_lines);
let _ = site_events.insert(*position, event);
}
}
#[allow(clippy::type_complexity)]
pub fn compute(
mut self,
) -> Result<(Vec<T>, impl ExactSizeIterator<Item = (u32, Vec<u32>)>), LinestringError> {
let mut site_events = self.site_events.take().ok_or_else(|| {
LinestringError::InternalError(
"compute() could not take ownership of 'site_events'".to_string(),
)
})?;
let mut active_lines =
FxHashSet::<u32>::with_capacity_and_hasher(site_events.len(), Default::default());
let mut result: BTreeMap<SiteEventKey<T>, Vec<u32>> = BTreeMap::default();
let mut neighbour_priority = MinMax::new();
let mut connected_priority = MinMaxSlope::new();
loop {
if let Some((key, event)) = {
site_events.pop_first()
} {
self.handle_event(
&key,
&event,
&mut active_lines,
&mut neighbour_priority,
&mut connected_priority,
&mut site_events,
&mut result,
)?;
if !result.is_empty() && self.stop_at_first_intersection {
break;
}
} else {
self.sweepline_pos = T::new_2d(T::Scalar::max_value(), T::Scalar::max_value());
self.sweepline_pos_index = None;
break;
}
}
Ok((
self.vertices,
result.into_iter().map(|x| (x.0.index.unwrap(), x.1)),
))
}
#[inline(always)]
#[allow(clippy::too_many_arguments)]
fn handle_event(
&mut self,
key: &SiteEventKey<T>,
event: &SiteEvent,
active_lines: &mut FxHashSet<u32>,
neighbour_priority: &mut MinMax<T>,
connected_priority: &mut MinMaxSlope<T>,
site_events: &mut BTreeMap<SiteEventKey<T>, SiteEvent>,
result: &mut BTreeMap<SiteEventKey<T>, Vec<u32>>,
) -> Result<(), LinestringError> {
self.sweepline_pos = key.pos;
self.sweepline_pos_index = key.index;
let removed_active_lines = event.drop.iter().flatten().count();
let added_active_lines = event.add.iter().flatten().count();
let intersections_found = event.intersection.iter().flatten().count();
if self.ignore_end_point_intersections {
if intersections_found > 0 {
if removed_active_lines > 0 {
self.report_intersections_to_result(
result,
self.sweepline_pos_index.unwrap(),
event.drop.iter().flatten(),
);
}
if added_active_lines > 0 {
self.report_intersections_to_result(
result,
self.sweepline_pos_index.unwrap(),
event.add.iter().flatten(),
);
}
self.report_intersections_to_result(
result,
self.sweepline_pos_index.unwrap_or_else(|| {
println!("self.sweepline_pos:{:?}", self.sweepline_pos);
panic!();
}),
event.intersection.iter().flatten(),
);
}
} else if removed_active_lines + added_active_lines + intersections_found > 1 {
if removed_active_lines > 0 {
self.report_intersections_to_result(
result,
self.sweepline_pos_index.unwrap(),
event.drop.iter().flatten(),
);
}
if added_active_lines > 0 {
self.report_intersections_to_result(
result,
self.sweepline_pos_index.unwrap(),
event.add.iter().flatten(),
);
}
if intersections_found > 0 {
self.report_intersections_to_result(
result,
self.sweepline_pos_index.unwrap(),
event.intersection.iter().flatten(),
);
}
}
for line_index in event.drop.iter().flatten() {
let _ = active_lines.remove(line_index);
}
neighbour_priority.clear();
'active_lines: for line_index in active_lines.iter() {
for i in event.intersection.iter().flatten() {
if i == line_index {
continue 'active_lines;
}
}
if let Some((intersection_x, intersection_slope)) = sweepline_intersection(
&self.vertices,
self.sweepline_pos,
*self.lines.get(*line_index as usize).unwrap(),
) {
neighbour_priority.update(
self.sweepline_pos.x(),
intersection_x,
intersection_slope,
*line_index,
);
}
}
for l in event.add.iter().flatten() {
let _ = active_lines.insert(*l);
}
if intersections_found + added_active_lines == 0 {
if !neighbour_priority.slope.candidates_left.is_empty()
&& !neighbour_priority.slope.candidates_right.is_empty()
{
self.find_new_events(
&neighbour_priority.slope.candidates_left,
&neighbour_priority.slope.candidates_right,
site_events,
)?;
}
} else {
connected_priority.clear();
for l in event.add.iter().flatten() {
connected_priority.update_both(&self.vertices, *l, &self.lines);
}
for l in event.intersection.iter().flatten() {
connected_priority.update_both(&self.vertices, *l, &self.lines);
}
if !neighbour_priority.slope.candidates_left.is_empty() {
self.find_new_events(
&neighbour_priority.slope.candidates_left,
&connected_priority.candidates_left,
site_events,
)?;
}
if !neighbour_priority.slope.candidates_right.is_empty() {
self.find_new_events(
&neighbour_priority.slope.candidates_right,
&connected_priority.candidates_right,
site_events,
)?;
}
}
Ok(())
}
fn find_new_events(
&mut self,
left: &[u32],
right: &[u32],
site_events: &mut BTreeMap<SiteEventKey<T>, SiteEvent>,
) -> Result<(), LinestringError> {
for left_i in left.iter() {
for right_i in right.iter() {
let left_l = self.lines[*left_i as usize];
let right_l = self.lines[*right_i as usize];
let left_l_end = self.vertices[left_l.1 as usize];
let right_l_end = self.vertices[right_l.1 as usize];
if left_l_end.is_ulps_eq(
right_l_end,
T::Scalar::default_epsilon(),
T::Scalar::default_max_ulps(),
) {
continue;
}
let left_l_as_line = Line2::new(
self.vertices[left_l.0 as usize],
self.vertices[left_l.1 as usize],
);
let right_l_as_line = Line2::new(
self.vertices[right_l.0 as usize],
self.vertices[right_l.1 as usize],
);
if let Some(intersection_p) = left_l_as_line.intersection_point(right_l_as_line) {
let intersection_p = intersection_p.single();
if !intersection_p.x().is_finite() || !intersection_p.y().is_finite() {
return Err(LinestringError::InvalidData(format!(
"Input data has intersection at invalid position x:{:?}, y:{:?}",
intersection_p.x(),
intersection_p.y()
)));
}
if intersection_p.y() >= self.sweepline_pos.y()
&& !(intersection_p.y() == self.sweepline_pos.y()
&& intersection_p.x() < self.sweepline_pos.x())
&& !intersection_p.is_ulps_eq(
self.sweepline_pos,
T::Scalar::default_epsilon(),
T::Scalar::default_max_ulps(),
)
{
let intersection_p_index = self.vertices.len() as u32;
self.vertices.push(intersection_p);
self.add_intersection_event(
site_events,
&SiteEventKey {
pos: intersection_p,
index: Some(intersection_p_index),
},
&[*left_i, *right_i],
)
}
}
}
}
Ok(())
}
fn report_intersections_to_result<'b, I: Iterator<Item = &'b u32>>(
&mut self,
result: &mut BTreeMap<SiteEventKey<T>, Vec<u32>>,
pos_index: u32,
intersecting_lines: I,
) {
let pos = self.vertices[pos_index as usize];
let key = SiteEventKey {
pos,
index: Some(pos_index),
};
let value = if let Some(value) = result.get_mut(&key) {
value
} else {
let _ = result.insert(key, Vec::default());
result.get_mut(&key).unwrap()
};
for line in intersecting_lines {
value.push(*line);
}
value.sort_unstable();
}
}