use crate::types::{Coord3D, EdgeSources, Line3D};
use crate::utils::parallel::{par_flat_map, par_sort_unstable, par_zip_for_each};
use crate::utils::{compare_angular, z_order_index};
use geo_types::{Coord, LineString};
#[cfg(feature = "parallel")]
use rayon::prelude::*;
use std::cell::RefCell;
use std::cmp::Ordering;
use std::collections::HashMap;
thread_local! {
static NEXT_POINTERS: RefCell<Vec<usize>> = const { RefCell::new(Vec::new()) };
}
pub type NodeId = usize;
pub type EdgeId = usize;
pub type DirEdgeId = usize;
pub(crate) struct ExtractedRing {
pub coords: Vec<Coord3D>,
pub line_ids: Vec<u32>,
pub source_line_ids: Vec<u32>,
pub edge_keys: Vec<(NodeId, NodeId)>,
pub node_ids: Vec<NodeId>,
}
#[derive(Clone, Debug)]
pub(crate) struct Edge {
pub(crate) line: Line3D,
pub(crate) sources: EdgeSources,
pub(crate) dir_edges: [DirEdgeId; 2],
pub(crate) is_marked: bool,
pub(crate) deleted: bool,
}
#[derive(Clone, Debug)]
pub(crate) struct DirectedEdge {
pub(crate) src: NodeId,
pub(crate) dst: NodeId,
pub(crate) edge_idx: EdgeId,
pub(crate) sym_idx: DirEdgeId,
pub(crate) is_visited: bool,
pub(crate) is_marked: bool,
}
#[derive(Clone)]
pub struct PlanarGraph {
pub(crate) nodes_x: Vec<f64>,
pub(crate) nodes_y: Vec<f64>,
pub(crate) nodes_z: Vec<f64>,
pub(crate) nodes_outgoing: Vec<Vec<DirEdgeId>>,
pub(crate) nodes_degree: Vec<usize>,
pub(crate) nodes_marked: Vec<bool>,
pub(crate) edges: Vec<Edge>,
pub(crate) directed_edges: Vec<DirectedEdge>,
pub(crate) node_map: HashMap<NodeKey, NodeId>,
}
#[derive(PartialEq, Eq, Hash, Clone, Copy)]
pub(crate) struct NodeKey(i64, i64);
impl From<Coord<f64>> for NodeKey {
fn from(c: Coord<f64>) -> Self {
NodeKey(c.x.to_bits() as i64, c.y.to_bits() as i64)
}
}
struct NodeEntry {
z_idx: u64, c: Coord3D,
}
impl PartialEq for NodeEntry {
fn eq(&self, other: &Self) -> bool {
self.z_idx == other.z_idx && self.c.x == other.c.x && self.c.y == other.c.y
}
}
impl Eq for NodeEntry {}
impl PartialOrd for NodeEntry {
fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
Some(self.cmp(other))
}
}
impl Ord for NodeEntry {
fn cmp(&self, other: &Self) -> Ordering {
self.z_idx.cmp(&other.z_idx).then_with(|| {
self.c
.x
.total_cmp(&other.c.x)
.then(self.c.y.total_cmp(&other.c.y))
})
}
}
fn create_edge_components(
i: usize,
u: NodeId,
v: NodeId,
line: Line3D,
sources: EdgeSources,
edges_start_len: usize,
dir_edges_start_len: usize,
) -> (
NodeId,
NodeId,
DirEdgeId,
DirEdgeId,
DirectedEdge,
DirectedEdge,
Edge,
) {
let edge_idx = edges_start_len + i;
let de_u_v_idx = dir_edges_start_len + 2 * i;
let de_v_u_idx = dir_edges_start_len + 2 * i + 1;
let de_u_v = DirectedEdge {
src: u,
dst: v,
edge_idx,
sym_idx: de_v_u_idx,
is_visited: false,
is_marked: false,
};
let de_v_u = DirectedEdge {
src: v,
dst: u,
edge_idx,
sym_idx: de_u_v_idx,
is_visited: false,
is_marked: false,
};
let edge = Edge {
line,
sources,
dir_edges: [de_u_v_idx, de_v_u_idx],
is_marked: false,
deleted: false,
};
(u, v, de_u_v_idx, de_v_u_idx, de_u_v, de_v_u, edge)
}
impl Default for PlanarGraph {
fn default() -> Self {
Self::new()
}
}
impl PlanarGraph {
pub fn new() -> Self {
Self {
nodes_x: Vec::new(),
nodes_y: Vec::new(),
nodes_z: Vec::new(),
nodes_outgoing: Vec::new(),
nodes_degree: Vec::new(),
nodes_marked: Vec::new(),
edges: Vec::new(),
directed_edges: Vec::new(),
node_map: HashMap::new(),
}
}
pub(crate) fn clear(&mut self) {
self.nodes_x.clear();
self.nodes_y.clear();
self.nodes_z.clear();
self.nodes_outgoing.clear();
self.nodes_degree.clear();
self.nodes_marked.clear();
self.edges.clear();
self.directed_edges.clear();
self.node_map.clear();
}
pub fn add_node(&mut self, coord: Coord3D) -> NodeId {
let key = NodeKey::from(coord.to_coord_2d());
if let Some(&id) = self.node_map.get(&key) {
return id;
}
let id = self.nodes_x.len();
self.nodes_x.push(coord.x);
self.nodes_y.push(coord.y);
self.nodes_z.push(coord.z);
self.nodes_outgoing.push(Vec::new());
self.nodes_degree.push(0);
self.nodes_marked.push(false);
self.node_map.insert(key, id);
id
}
pub fn bulk_load(&mut self, lines: Vec<Line3D>) {
if lines.is_empty() {
return;
}
let to_entries = |line: &Line3D| {
[
NodeEntry {
z_idx: z_order_index(line.start.to_coord_2d()),
c: line.start,
},
NodeEntry {
z_idx: z_order_index(line.end.to_coord_2d()),
c: line.end,
},
]
};
let mut entries: Vec<NodeEntry> = par_flat_map(&lines, to_entries);
par_sort_unstable(&mut entries);
entries.dedup_by(|a, b| a.c.x == b.c.x && a.c.y == b.c.y);
let start_node_idx = self.nodes_x.len();
self.nodes_x.reserve(entries.len());
self.nodes_y.reserve(entries.len());
self.nodes_z.reserve(entries.len());
self.nodes_outgoing.reserve(entries.len());
self.nodes_degree.reserve(entries.len());
self.nodes_marked.reserve(entries.len());
for entry in &entries {
self.nodes_x.push(entry.c.x);
self.nodes_y.push(entry.c.y);
self.nodes_z.push(entry.c.z);
self.nodes_outgoing.push(Vec::new());
self.nodes_degree.push(0);
self.nodes_marked.push(false);
}
let get_node_id = |pt: Coord3D| -> Option<NodeId> {
let z_pt = z_order_index(pt.to_coord_2d());
let idx_res = entries.binary_search_by(|probe| {
probe
.z_idx
.cmp(&z_pt)
.then_with(|| probe.c.x.total_cmp(&pt.x).then(probe.c.y.total_cmp(&pt.y)))
});
match idx_res {
Ok(i) => Some(start_node_idx + i),
Err(_) => None,
}
};
let mut valid_edges = Vec::with_capacity(lines.len());
let mut degrees = vec![0usize; self.nodes_x.len()];
for line in lines {
let p0 = line.start;
let p1 = line.end;
if p0.x == p1.x && p0.y == p1.y {
continue;
}
let u_opt = get_node_id(p0);
let v_opt = get_node_id(p1);
if let (Some(u), Some(v)) = (u_opt, v_opt) {
valid_edges.push((u, v, line));
}
}
valid_edges.sort_unstable_by(|(u1, v1, l1), (u2, v2, l2)| {
let key1 = ((*u1).min(*v1), (*u1).max(*v1));
let key2 = ((*u2).min(*v2), (*u2).max(*v2));
key1.cmp(&key2)
.then(l1.line_id.cmp(&l2.line_id))
.then(u1.cmp(u2))
.then(v1.cmp(v2))
});
let mut dissolved: Vec<(NodeId, NodeId, Line3D, EdgeSources)> = Vec::new();
for (u, v, line) in valid_edges {
let key = (u.min(v), u.max(v));
if let Some((last_u, last_v, _, sources)) = dissolved.last_mut() {
if ((*last_u).min(*last_v), (*last_u).max(*last_v)) == key {
sources.merge_line_id(line.line_id);
continue;
}
}
dissolved.push((u, v, line, EdgeSources::from_line_id(line.line_id)));
}
for (u, v, _, _) in &dissolved {
degrees[*u] += 1;
degrees[*v] += 1;
}
par_zip_for_each(
&mut self.nodes_outgoing,
°rees,
|adj: &mut Vec<usize>, deg: &usize| {
adj.reserve(*deg);
},
);
self.edges.reserve(dissolved.len());
self.directed_edges.reserve(dissolved.len() * 2);
let edges_start_len = self.edges.len();
let directed_edges_start_len = self.directed_edges.len();
for (i, (u, v, line, sources)) in dissolved.into_iter().enumerate() {
let (u, v, de_u_v_idx, de_v_u_idx, de_u_v, de_v_u, edge) = create_edge_components(
i,
u,
v,
line,
sources,
edges_start_len,
directed_edges_start_len,
);
self.directed_edges.push(de_u_v);
self.directed_edges.push(de_v_u);
self.edges.push(edge);
self.nodes_outgoing[u].push(de_u_v_idx);
self.nodes_degree[u] += 1;
self.nodes_outgoing[v].push(de_v_u_idx);
self.nodes_degree[v] += 1;
}
}
pub fn add_line(&mut self, line: Line3D) -> EdgeId {
let p0 = line.start;
let p1 = line.end;
let u = self.add_node(p0);
let v = self.add_node(p1);
if let Some(edge_idx) = self.nodes_outgoing[u]
.iter()
.find_map(|&directed_edge_idx| {
let directed_edge = &self.directed_edges[directed_edge_idx];
(directed_edge.dst == v && !self.edges[directed_edge.edge_idx].deleted)
.then_some(directed_edge.edge_idx)
})
{
let edge = &mut self.edges[edge_idx];
edge.sources.merge_line_id(line.line_id);
edge.line.line_id = edge.sources.line_ids[0];
return edge_idx;
}
let edge_idx = self.edges.len();
let de_u_v_idx = self.directed_edges.len();
let de_v_u_idx = self.directed_edges.len() + 1;
let de_u_v = DirectedEdge {
src: u,
dst: v,
edge_idx,
sym_idx: de_v_u_idx,
is_visited: false,
is_marked: false,
};
let de_v_u = DirectedEdge {
src: v,
dst: u,
edge_idx,
sym_idx: de_u_v_idx,
is_visited: false,
is_marked: false,
};
self.directed_edges.push(de_u_v);
self.directed_edges.push(de_v_u);
self.edges.push(Edge {
line,
sources: EdgeSources::from_line_id(line.line_id),
dir_edges: [de_u_v_idx, de_v_u_idx],
is_marked: false,
deleted: false,
});
self.nodes_outgoing[u].push(de_u_v_idx);
self.nodes_degree[u] += 1;
self.nodes_outgoing[v].push(de_v_u_idx);
self.nodes_degree[v] += 1;
edge_idx
}
pub fn remove_line_by_id(&mut self, line_id: u32) -> bool {
for edge in &mut self.edges {
if !edge.deleted && edge.sources.remove_line_id(line_id) {
if edge.sources.line_ids.is_empty() {
edge.deleted = true;
} else if edge.line.line_id == line_id {
edge.line.line_id = edge.sources.line_ids[0];
}
return true;
}
}
false
}
pub fn reset_traversal_state(&mut self) {
for d in &mut self.nodes_degree {
*d = 0;
}
for (i, outgoing) in self.nodes_outgoing.iter().enumerate() {
for &de_idx in outgoing {
let de = &self.directed_edges[de_idx];
if !self.edges[de.edge_idx].deleted {
self.nodes_degree[i] += 1;
}
}
}
for m in &mut self.nodes_marked {
*m = false;
}
for de in &mut self.directed_edges {
de.is_visited = false;
de.is_marked = false;
}
for edge in &mut self.edges {
edge.is_marked = false;
}
}
pub fn add_line_string(&mut self, line: LineString<f64>) {
if line.0.is_empty() {
return;
}
let coords = &line.0;
for w in coords.windows(2) {
let p0 = w[0];
let p1 = w[1];
if p0.x == p1.x && p0.y == p1.y {
continue;
}
self.add_line(Line3D::new(p0.into(), p1.into(), 0));
}
}
pub fn sort_edges(&mut self) {
let nodes_x = &self.nodes_x;
let nodes_y = &self.nodes_y;
let directed_edges = &self.directed_edges;
#[cfg(feature = "parallel")]
self.nodes_outgoing
.par_iter_mut()
.enumerate()
.for_each(|(src_idx, adj)| {
adj.retain(|&idx| !self.edges[self.directed_edges[idx].edge_idx].deleted);
let center = Coord {
x: nodes_x[src_idx],
y: nodes_y[src_idx],
};
adj.sort_by(|&a_idx, &b_idx| {
let a_de = &directed_edges[a_idx];
let b_de = &directed_edges[b_idx];
let dst_a_idx = a_de.dst;
let dst_b_idx = b_de.dst;
let target_a = Coord {
x: nodes_x[dst_a_idx],
y: nodes_y[dst_a_idx],
};
let target_b = Coord {
x: nodes_x[dst_b_idx],
y: nodes_y[dst_b_idx],
};
compare_angular(center, target_a, target_b)
});
});
#[cfg(not(feature = "parallel"))]
self.nodes_outgoing
.iter_mut()
.enumerate()
.for_each(|(src_idx, adj)| {
adj.retain(|&idx| !self.edges[self.directed_edges[idx].edge_idx].deleted);
let center = Coord {
x: nodes_x[src_idx],
y: nodes_y[src_idx],
};
adj.sort_by(|&a_idx, &b_idx| {
let a_de = &directed_edges[a_idx];
let b_de = &directed_edges[b_idx];
let dst_a_idx = a_de.dst;
let dst_b_idx = b_de.dst;
let target_a = Coord {
x: nodes_x[dst_a_idx],
y: nodes_y[dst_a_idx],
};
let target_b = Coord {
x: nodes_x[dst_b_idx],
y: nodes_y[dst_b_idx],
};
compare_angular(center, target_a, target_b)
});
});
}
pub fn prune_dangles(&mut self) -> Vec<Vec<Coord3D>> {
let mut dangles = Vec::new();
let mut to_process: Vec<NodeId> = self
.nodes_degree
.iter()
.enumerate()
.filter(|(i, &d)| d == 1 && !self.nodes_marked[*i])
.map(|(i, _)| i)
.collect();
while let Some(node_idx) = to_process.pop() {
if self.nodes_degree[node_idx] != 1 {
continue;
}
self.nodes_marked[node_idx] = true;
self.nodes_degree[node_idx] = 0;
let mut edge_found = false;
let mut neighbor_idx = 0;
let mut found_de_idx = None;
for &de_idx in &self.nodes_outgoing[node_idx] {
let de = &self.directed_edges[de_idx];
if !de.is_marked && !self.edges[de.edge_idx].deleted {
found_de_idx = Some(de_idx);
break;
}
}
if let Some(de_idx) = found_de_idx {
self.directed_edges[de_idx].is_marked = true;
let sym_idx = self.directed_edges[de_idx].sym_idx;
self.directed_edges[sym_idx].is_marked = true;
let edge_idx = self.directed_edges[de_idx].edge_idx;
let line = self.edges[edge_idx].line;
dangles.push(vec![line.start, line.end]);
neighbor_idx = self.directed_edges[de_idx].dst;
edge_found = true;
}
if edge_found && self.nodes_degree[neighbor_idx] > 0 {
self.nodes_degree[neighbor_idx] -= 1;
if self.nodes_degree[neighbor_idx] == 1 && !self.nodes_marked[neighbor_idx] {
to_process.push(neighbor_idx);
}
}
}
dangles
}
pub fn delete_cut_edges(&mut self) -> Vec<Vec<Coord3D>> {
NEXT_POINTERS.with(|cell| {
let mut next_pointers = cell.borrow_mut();
next_pointers.clear();
next_pointers.resize(self.directed_edges.len(), usize::MAX);
self.compute_next_cw_edges(&mut next_pointers);
let mut labels = vec![-1_i64; self.directed_edges.len()];
self.find_and_label_maximal_rings(&next_pointers, &mut labels);
let mut cuts = Vec::new();
for edge in &self.edges {
let [forward, reverse] = edge.dir_edges;
if edge.deleted
|| self.directed_edges[forward].is_marked
|| self.directed_edges[reverse].is_marked
|| labels[forward] != labels[reverse]
{
continue;
}
self.directed_edges[forward].is_marked = true;
self.directed_edges[reverse].is_marked = true;
cuts.push(vec![edge.line.start, edge.line.end]);
}
cuts
})
}
pub fn get_edge_rings(&mut self) -> Vec<(Vec<Coord3D>, Vec<u32>)> {
self.get_edge_rings_with_graph_ids(false, false)
.into_iter()
.map(|ring| (ring.coords, ring.line_ids))
.collect()
}
pub(crate) fn get_edge_rings_with_graph_ids(
&mut self,
include_graph_ids: bool,
include_source_ids: bool,
) -> Vec<ExtractedRing> {
NEXT_POINTERS.with(|cell| {
let mut next_pointers = cell.borrow_mut();
next_pointers.clear();
next_pointers.resize(self.directed_edges.len(), usize::MAX);
let mut labels = vec![-1_i64; self.directed_edges.len()];
self.compute_next_cw_edges(&mut next_pointers);
let maximal_ring_starts =
self.find_and_label_maximal_rings(&next_pointers, &mut labels);
self.convert_maximal_to_minimal_rings(
&maximal_ring_starts,
&mut next_pointers,
&labels,
);
self.extract_valid_rings(&next_pointers, include_graph_ids, include_source_ids)
})
}
fn compute_next_cw_edges(&self, next_pointers: &mut [usize]) {
let mut valid_edges = Vec::new();
for outgoing in &self.nodes_outgoing {
valid_edges.clear();
valid_edges.extend(outgoing.iter().copied().filter(|&idx| {
let de = &self.directed_edges[idx];
!de.is_marked && !self.edges[de.edge_idx].deleted
}));
if valid_edges.is_empty() {
continue;
}
let mut next = *valid_edges.last().unwrap();
for &curr in &valid_edges {
next_pointers[curr] = next;
next = curr;
}
}
}
fn find_and_label_maximal_rings(
&self,
next_pointers: &[usize],
labels: &mut [i64],
) -> Vec<DirEdgeId> {
let mut maximal_ring_starts = Vec::new();
let mut curr_label = 1_i64;
for start_de_idx in 0..self.directed_edges.len() {
if self.directed_edges[start_de_idx].is_marked || labels[start_de_idx] >= 0 {
continue;
}
maximal_ring_starts.push(start_de_idx);
let mut curr = start_de_idx;
loop {
if labels[curr] >= 0 {
break;
}
labels[curr] = curr_label;
let next = next_pointers[self.directed_edges[curr].sym_idx];
if next == usize::MAX || next == start_de_idx {
break;
}
curr = next;
}
curr_label += 1;
}
maximal_ring_starts
}
fn convert_maximal_to_minimal_rings(
&self,
maximal_ring_starts: &[DirEdgeId],
next_pointers: &mut [usize],
labels: &[i64],
) {
let mut intersection_nodes = Vec::<NodeId>::new();
let mut seen_intersection_nodes = vec![false; self.nodes_x.len()];
for &start_de_idx in maximal_ring_starts {
let ring_label = labels[start_de_idx];
if ring_label < 0 {
continue;
}
intersection_nodes.clear();
let mut curr = start_de_idx;
loop {
let node = self.directed_edges[curr].src;
let mut degree_for_label = 0;
for &out_de in &self.nodes_outgoing[node] {
if labels[out_de] == ring_label {
degree_for_label += 1;
}
}
if degree_for_label > 1 && !seen_intersection_nodes[node] {
seen_intersection_nodes[node] = true;
intersection_nodes.push(node);
}
let next = next_pointers[self.directed_edges[curr].sym_idx];
if next == usize::MAX || next == start_de_idx {
break;
}
curr = next;
}
for &node in &intersection_nodes {
let outgoing = &self.nodes_outgoing[node];
let mut first_out: Option<DirEdgeId> = None;
let mut prev_in: Option<DirEdgeId> = None;
for &de_idx in outgoing.iter().rev() {
let sym_idx = self.directed_edges[de_idx].sym_idx;
let out_de = (labels[de_idx] == ring_label).then_some(de_idx);
let in_de = (labels[sym_idx] == ring_label).then_some(sym_idx);
if out_de.is_none() && in_de.is_none() {
continue;
}
if let Some(in_de_idx) = in_de {
prev_in = Some(in_de_idx);
}
if let Some(out_de_idx) = out_de {
if let Some(prev_in_idx) = prev_in.take() {
next_pointers[self.directed_edges[prev_in_idx].sym_idx] = out_de_idx;
}
if first_out.is_none() {
first_out = Some(out_de_idx);
}
}
}
if let (Some(prev_in_idx), Some(first_out_idx)) = (prev_in, first_out) {
next_pointers[self.directed_edges[prev_in_idx].sym_idx] = first_out_idx;
}
seen_intersection_nodes[node] = false;
}
}
}
fn extract_valid_rings(
&mut self,
next_pointers: &[usize],
include_graph_ids: bool,
include_source_ids: bool,
) -> Vec<ExtractedRing> {
for de in &mut self.directed_edges {
de.is_visited = false;
}
let mut ring_edges = Vec::new();
let mut rings = Vec::new();
for start_de_idx in 0..self.directed_edges.len() {
if self.directed_edges[start_de_idx].is_visited
|| self.directed_edges[start_de_idx].is_marked
{
continue;
}
if self.edges[self.directed_edges[start_de_idx].edge_idx].deleted {
continue;
}
ring_edges.clear();
let mut curr_de_idx = start_de_idx;
let mut is_valid_ring = true;
loop {
let curr_de = &mut self.directed_edges[curr_de_idx];
curr_de.is_visited = true;
ring_edges.push(curr_de_idx);
let next_de_idx = next_pointers[self.directed_edges[curr_de_idx].sym_idx];
if next_de_idx == usize::MAX {
is_valid_ring = false;
break;
}
curr_de_idx = next_de_idx;
if curr_de_idx == start_de_idx {
break;
}
if self.directed_edges[curr_de_idx].is_visited {
is_valid_ring = false;
break;
}
}
if is_valid_ring && !ring_edges.is_empty() {
let mut coords = Vec::with_capacity(ring_edges.len() + 1);
let mut ids = Vec::with_capacity(ring_edges.len());
let mut source_ids = Vec::new();
let mut edge_keys = if include_graph_ids {
Vec::with_capacity(ring_edges.len())
} else {
Vec::new()
};
let mut node_ids = if include_graph_ids {
Vec::with_capacity(ring_edges.len())
} else {
Vec::new()
};
let start_node_idx = self.directed_edges[ring_edges[0]].src;
coords.push(Coord3D {
x: self.nodes_x[start_node_idx],
y: self.nodes_y[start_node_idx],
z: self.nodes_z[start_node_idx],
});
for &de_idx in &ring_edges {
let de = &self.directed_edges[de_idx];
let edge_idx = de.edge_idx;
ids.push(self.edges[edge_idx].line.line_id);
if include_source_ids {
source_ids.extend_from_slice(&self.edges[edge_idx].sources.line_ids);
}
if include_graph_ids {
edge_keys.push(if de.src < de.dst {
(de.src, de.dst)
} else {
(de.dst, de.src)
});
node_ids.push(de.src);
}
let dst_idx = de.dst;
coords.push(Coord3D {
x: self.nodes_x[dst_idx],
y: self.nodes_y[dst_idx],
z: self.nodes_z[dst_idx],
});
}
source_ids.sort_unstable();
source_ids.dedup();
rings.push(ExtractedRing {
coords,
line_ids: ids,
source_line_ids: source_ids,
edge_keys,
node_ids,
});
}
}
rings
}
}