use delaunator::{triangulate, Point as DelaunatorPoint}; use geo::{Point, Polygon, LineString, Area}; use itertools::Itertools; use rand::Rng; use rayon::prelude::*; use std::cmp::{min, max}; use std::collections::{HashMap, HashSet}; use std::sync::{Arc, Mutex}; #[derive(Debug, Clone, Copy, Hash, PartialEq, Eq)] pub struct Edge(usize, usize); #[derive(Debug)] pub struct TriangleData { pub index: usize, pub area: Option, pub terminal_edge: Option, pub vertices: Vec } #[derive(Debug)] pub struct GeometryData { pub triangles: Vec, pub edge_to_triangles: HashMap>, // Maps an edge to triangle indices pub edge_lengths: HashMap, // Edge lengths pub vertex_connections: HashMap>, // Direct connections between vertices, for DTSCAN } impl GeometryData { fn new() -> Self { GeometryData { triangles: Vec::new(), edge_to_triangles: HashMap::new(), edge_lengths: HashMap::new(), vertex_connections: HashMap::new(), // Adjusted for DTSCAN } } fn add_triangle(&mut self, index: usize, points: &[Point], tri_idx: &[usize], types: usize) { let point_a: Point = points[tri_idx[0]]; let point_b: Point = points[tri_idx[1]]; let point_c: Point = points[tri_idx[2]]; let mut vertices = vec![tri_idx[0], tri_idx[1], tri_idx[2]]; vertices.sort_unstable(); // Temporarily store edges_with_lengths for sorting and determining the terminal_edge. let mut edges_with_lengths_temp = [ (Edge(min(tri_idx[0], tri_idx[1]), max(tri_idx[0], tri_idx[1])), distance(point_a.x(), point_a.y(), point_b.x(), point_b.y())), (Edge(min(tri_idx[1], tri_idx[2]), max(tri_idx[1], tri_idx[2])), distance(point_b.x(), point_b.y(), point_c.x(), point_c.y())), (Edge(min(tri_idx[2], tri_idx[0]), max(tri_idx[2], tri_idx[0])), distance(point_c.x(), point_c.y(), point_a.x(), point_a.y())), ].to_vec(); // Sort edges by length to ensure the longest edge is identified. edges_with_lengths_temp.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap()); let terminal_edge: Option = edges_with_lengths_temp.first().map(|(edge, _)| *edge); let area: Option = if types == 0 || types == 2 { Some(Polygon::new(LineString::from(vec![ (point_a.x(), point_a.y()), (point_b.x(), point_b.y()), (point_c.x(), point_c.y()), (point_a.x(), point_a.y()), ]), vec![]).unsigned_area()) } else { None }; if types == 0 || types == 1 { for &(edge, length) in &edges_with_lengths_temp { self.vertex_connections.entry(edge.0).or_insert_with(HashSet::new).insert(edge.1); self.vertex_connections.entry(edge.1).or_insert_with(HashSet::new).insert(edge.0); self.edge_lengths.insert(edge, length); self.edge_to_triangles.entry(edge).or_default().push(index); } } else { // For types == 2, only update edge_lengths and edge_to_triangles. for &(edge, length) in &edges_with_lengths_temp { self.edge_lengths.insert(edge, length); self.edge_to_triangles.entry(edge).or_default().push(index); } } if types == 0 || types == 2 { self.triangles.push(TriangleData { index, area, terminal_edge, vertices }); } } } fn distance(x1: f32, y1: f32, x2: f32, y2: f32) -> f32 { ((x2 - x1).powi(2) + (y2 - y1).powi(2)).sqrt() } pub struct Xenobalanus { geometry_data: GeometryData, points: Vec>, triangles: Vec, } impl Xenobalanus { pub fn new() -> Self { Xenobalanus { geometry_data: GeometryData::new(), points: Vec::new(), triangles: Vec::new(), } } pub fn points(&self) -> Vec> { self.points.iter() .map(|point| vec![point.x(), point.y()]) .collect() } pub fn points_flat(&self) -> Vec { self.points.iter() .flat_map(|point| vec![point.x(), point.y()]) .collect() } pub fn triangles(&self) -> Vec { self.triangles.clone() } pub fn triangle_vertices(&self) -> Vec> { self.triangles.chunks(3).map(|chunk| { chunk.iter().map(|&index| index).collect() }).collect() } pub fn triangles_coordinates(&self) -> Vec> { self.triangles.chunks(3).map(|chunk| { chunk.iter().flat_map(|&index| { let point = &self.points[index]; vec![point.x(), point.y()] }).collect() }).collect() } // Additional methods moved into GeometryProcessor, operating on self.geometry_data pub fn random_points(&mut self, center: (f32, f32), side_length: f32, num_points: u32) { // generate random points in a square let min_x = center.0 - side_length / 2.0; let max_x = center.0 + side_length / 2.0; let min_y = center.1 - side_length / 2.0; let max_y = center.1 + side_length / 2.0; let mut rng: rand::prelude::ThreadRng = rand::thread_rng(); for _ in 0..num_points { let x = min_x + rng.gen_range(0.0..=1.0) as f32 * ( max_x - min_x); let y: f32 = min_y + rng.gen_range(0.0..=1.0) as f32 * ( max_y - min_y); self.points.push(Point::new(x, y)); } } pub fn delaunay(&mut self) { // Convert geo::Point to delaunator::Point for triangulation let delaunator_points: Vec = self.points.iter() .map(|point: &Point| DelaunatorPoint { x: point.x() as f64, y: point.y() as f64 }) .collect(); // Perform Delaunay triangulation let result: delaunator::Triangulation = triangulate(&delaunator_points); self.triangles = result.triangles } pub fn preprocess(&mut self, types: usize) { let geometry_data = Arc::new(Mutex::new(GeometryData::new())); self.triangles.par_chunks(3).enumerate().for_each(|(index, tri_idx)| { let gd = geometry_data.clone(); // Clone Arc for use in each thread gd.lock().unwrap().add_triangle(index, &self.points, tri_idx, types); }); self.geometry_data = Arc::try_unwrap(geometry_data).unwrap().into_inner().unwrap() } pub fn delfin( &self, min_area: f32, min_distance: f32, ) -> Vec> { // Sort all triangles by the longest terminal edge let triangles_sorted: Vec<(usize, f32)> = self.geometry_data.triangles.iter() .filter_map(|triangle_data| { // Only consider triangles with a terminal edge triangle_data.terminal_edge.map(|terminal_edge| { // Retrieve the length of the terminal edge if it exists self.geometry_data.edge_lengths.get(&terminal_edge) .map(|&length| (triangle_data.index, length)) }).flatten() }) .sorted_by(|a, b| b.1.partial_cmp(&a.1).unwrap()) // Sort in descending order by edge length .collect(); let mut void_polygons: Vec> = Vec::new(); let mut processed_triangles: HashSet = HashSet::new(); for &(triangle_index, terminal_edge_length) in &triangles_sorted { // Skip if this triangle has already been processed if processed_triangles.contains(&triangle_index) { continue; } // Continue if the terminal edge length is below the minimum distance threshold if terminal_edge_length < min_distance { continue; } // Retrieve triangles that share the terminal edge, continue if less than 2 triangles share it let triangle_data: &TriangleData = &self.geometry_data.triangles[triangle_index]; if let Some(terminal_edge) = triangle_data.terminal_edge { if let Some(connected_triangles) = self.geometry_data.edge_to_triangles.get(&terminal_edge) { // Proceed only if there are 2 or more triangles sharing the terminal edge if connected_triangles.len() < 2 { continue; } // Initialize the set with the current triangle and triangles directly connected via their terminal edge let mut triangle_set: HashSet = connected_triangles.iter().cloned().collect(); triangle_set.insert(triangle_index); processed_triangles.extend(&triangle_set); // Dynamically expand the set based on the terminal edge sharing criterion let mut triangles_to_expand: HashSet = triangle_set.clone(); while let Some(current_idx) = triangles_to_expand.iter().next().cloned() { // Remove the current triangle index from the set to avoid reprocessing triangles_to_expand.remove(¤t_idx); // Iterate over each triangle that shares a terminal edge for &neighbor_idx in connected_triangles { // Skip if this triangle has already been considered or processed if triangle_set.contains(&neighbor_idx) || processed_triangles.contains(&neighbor_idx) { continue; } // Safely access the neighbor triangle's data using its index if let Some(neighbor_data) = self.geometry_data.triangles.get(neighbor_idx) { // Check if the neighbor shares the same terminal edge // Directly compare the terminal edges as they are both Option if neighbor_data.terminal_edge == Some(terminal_edge) { // If they share the same terminal edge, include the neighbor in the current void polygon set triangle_set.insert(neighbor_idx); processed_triangles.insert(neighbor_idx); triangles_to_expand.insert(neighbor_idx); } } } } // Add the expanded set to void polygons void_polygons.push(triangle_set); } else { // If no connected triangles are found for the terminal edge, simply skip to the next triangle continue; } } } // Filter out void polygon sets void_polygons.retain(|poly_set: &HashSet| { // Calculate the total area of the polygon set by summing the areas of the triangles it contains. let total_area: f32 = poly_set.iter() .filter_map(|&idx| self.geometry_data.triangles.get(idx).and_then(|td| td.area)) .sum(); // Filter based on the area and the minimum number of triangles. total_area >= min_area && poly_set.len() >= 3 }); return void_polygons; } pub fn dtscan( &self, min_pts: usize, max_closeness: f32, ) -> Vec> { let mut clusters: Vec> = Vec::new(); let mut visited: HashSet = HashSet::new(); for (&vertex_idx, neighbors) in &self.geometry_data.vertex_connections { if visited.contains(&vertex_idx) { continue; } // Check if vertex is a core vertex based on the number of connections and edge lengths if neighbors.len() >= min_pts && neighbors.iter().all(|&n| { if let Some(&length) = self.geometry_data.edge_lengths.get(&Edge(min(vertex_idx, n), max(vertex_idx, n))) { length <= max_closeness } else { false } }) { let mut cluster: Vec = Vec::new(); let mut to_expand: Vec = vec![vertex_idx]; while let Some(current_vertex) = to_expand.pop() { if !visited.insert(current_vertex) { continue; } cluster.push(current_vertex); // Add neighbors that are within max_closeness to to_expand self.geometry_data.vertex_connections.get(¤t_vertex).map(|neighbors: &HashSet| { for &neighbor in neighbors { if let Some(&length) = self.geometry_data.edge_lengths.get(&Edge(min(current_vertex, neighbor), max(current_vertex, neighbor))) { if length <= max_closeness && !visited.contains(&neighbor) { to_expand.push(neighbor); } } } }); } if !cluster.is_empty() { clusters.push(cluster); // Add the constructed cluster to the list of clusters } } } clusters } }