From 1b4a538dff7c6af34254327f745bd1a2abd4277f Mon Sep 17 00:00:00 2001 From: randogoth Date: Wed, 27 Mar 2024 11:08:13 +0200 Subject: [PATCH] turned into class --- readme.md | 16 ++- src/lib.rs | 375 +++++++++++++++++++++++++++++------------------------ 2 files changed, 218 insertions(+), 173 deletions(-) diff --git a/readme.md b/readme.md index 3b88dac..dae142b 100644 --- a/readme.md +++ b/readme.md @@ -14,7 +14,7 @@ DTSCAN was developed by Kim, Jongwon, and Jeongho Cho. "[Delaunay triangulation- ## Functions -The combined Xenobalanus implementation is comprised of several key functions: +The Xenobalanus class is comprised of several key methods: - `random_points`: Generates uniformly distributed random points for testing. - `delaunay`: A wrapper of the [Delaunator crate](https://docs.rs/delaunator/latest/delaunator/). Performs Delaunay Triangulation on a given set of points to find their triangular connections. @@ -29,31 +29,33 @@ Below is an example code snippet that demonstrates the workflow. This example ge ```rust use geo::Point; use std::collections::{HashSet}; -use xenobalanus::{delaunay, random_points, preprocess, dtscan, delfin, GeometryData}; +use xenobalanus; fn main() { // Define test area and random points let dots: u32 = 10000; let side_length: f32 = 10000.0; - let points: Vec> = random_points((0.0, 0.0), side_length, dots); + let mut xeno = Xenobalanus::new(); + xeno.random_points((0.0, 0.0), side_length, dots); println!("Generated {:#?} random dots", dots); // Run Delaunay triangulation - let triangles_indices: Vec = delaunay(&points); + xeno.delaunay(); println!("Generated Delaunay triangulation"); // Pre-process triangles - let geometry_data: GeometryData = preprocess(&points, &triangles_indices, 0); + xeno.preprocess(0); // Execute delfin function with the generated GeometryData let min_area: f32 = 1000.0; // threshold for voidness let min_distance: f32 = 200.0; // threshold for minimum distance - let void_polygons: Vec> = delfin(&geometry_data, min_area, min_distance); + let void_polygons: Vec> = xeno.delfin(min_area, min_distance); println!("Found {:#?} Voids", void_polygons.len()); // Execute DTSCAN with the prepared data let min_pts: usize = 5; // threshold for minimum number of points let max_closeness: f32 = 100.5; // threshold for maximum closeness - let clusters: Vec> = dtscan(&geometry_data, min_pts, max_closeness); + let clusters: Vec> = xeno.dtscan(min_pts, max_closeness); println!("Found {:#?} Attractors", clusters.len()); } +``` diff --git a/src/lib.rs b/src/lib.rs index 1beb681..3aeadbf 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -35,7 +35,6 @@ impl GeometryData { 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]]; @@ -89,200 +88,244 @@ impl GeometryData { vertices }); } - } + } + } fn distance(x1: f32, y1: f32, x2: f32, y2: f32) -> f32 { ((x2 - x1).powi(2) + (y2 - y1).powi(2)).sqrt() } -pub fn random_points(center: (f32, f32), side_length: f32, num_points: u32) -> Vec> { - // 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 points: Vec> = Vec::with_capacity(num_points as usize); - 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); - points.push(Point::new(x, y)); - } - - points +pub struct Xenobalanus { + geometry_data: GeometryData, + points: Vec>, + triangles: Vec, } -pub fn delaunay(points: &Vec>) -> Vec { - // Convert geo::Point to delaunator::Point for triangulation - let delaunator_points: Vec = points.iter() +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 + } - // Return the indices of points in the triangles - result.triangles -} + pub fn preprocess(&mut self, types: usize) { + let geometry_data = Arc::new(Mutex::new(GeometryData::new())); -pub fn preprocess(points: &[Point], triangles: &[usize], types: usize) -> GeometryData { - 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 - 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); + }); - gd.lock().unwrap().add_triangle(index, points, tri_idx, types); - }); + self.geometry_data = Arc::try_unwrap(geometry_data).unwrap().into_inner().unwrap() + } - Arc::try_unwrap(geometry_data).unwrap().into_inner().unwrap() -} - -pub fn delfin( - geometry_data: &GeometryData, - min_area: f32, - min_distance: f32, -) -> Vec> { - - // Sort all triangles by the longest terminal edge - let triangles_sorted: Vec<(usize, f32)> = 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 - 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; - } + pub fn delfin( + &self, + min_area: f32, + min_distance: f32, + ) -> Vec> { - // 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 = &geometry_data.triangles[triangle_index]; - if let Some(terminal_edge) = triangle_data.terminal_edge { - if let Some(connected_triangles) = 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; - } + // 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(); - // 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); + 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; + } - // 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) = 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 if the terminal edge length is below the minimum distance threshold + if terminal_edge_length < min_distance { 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| geometry_data.triangles.get(idx).and_then(|td| td.area)) - .sum(); + + // 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; + } - // 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( - geometry_data: &GeometryData, - min_pts: usize, - max_closeness: f32, -) -> Vec> { - let mut clusters: Vec> = Vec::new(); - let mut visited: HashSet = HashSet::new(); - - for (&vertex_idx, neighbors) in &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) = 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 - geometry_data.vertex_connections.get(¤t_vertex).map(|neighbors: &HashSet| { - for &neighbor in neighbors { - if let Some(&length) = 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); + // 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); + } } } } - }); - } - - if !cluster.is_empty() { - clusters.push(cluster); // Add the constructed cluster to the list of clusters + + // 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; + } - clusters + 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 + } } \ No newline at end of file