diff --git a/src/main.rs b/src/main.rs index 6a0a957..d6bae7a 100644 --- a/src/main.rs +++ b/src/main.rs @@ -14,10 +14,10 @@ struct Edge(usize, usize); #[derive(Debug)] struct TriangleData { index: usize, - area: f32, - terminal_edge: Edge, - edges_with_lengths: Vec<(Edge, f32)>, - node_connections: HashSet, + area: Option, + terminal_edge: Option, + edges_with_lengths: Option>, + node_connections: Option>, } #[derive(Debug)] @@ -38,48 +38,61 @@ impl GeometryData { } } - // Function to add a triangle to the GeometryData - fn add_triangle(&mut self, index: usize, points: &[Point], tri_idx: &[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 edges_with_lengths = [ - (Edge(min(tri_idx[0], tri_idx[1]), max(tri_idx[0], tri_idx[1])), point_a.euclidean_distance(&point_b)), - (Edge(min(tri_idx[1], tri_idx[2]), max(tri_idx[1], tri_idx[2])), point_b.euclidean_distance(&point_c)), - (Edge(min(tri_idx[2], tri_idx[0]), max(tri_idx[2], tri_idx[0])), point_c.euclidean_distance(&point_a)), - ]; - - let mut edges_sorted = edges_with_lengths.to_vec(); - edges_sorted.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap()); - - let terminal_edge = edges_sorted[0].0; - let area = 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(); - - let node_connections: HashSet = tri_idx.iter().cloned().collect::>(); - + fn add_triangle(&mut self, index: usize, points: &[Point], tri_idx: &[usize], types: usize) { + let point_a = points[tri_idx[0]]; + let point_b = points[tri_idx[1]]; + let point_c = points[tri_idx[2]]; + + let edges_with_lengths: Option> = if types == 0 || types == 2 { + Some([ + (Edge(min(tri_idx[0], tri_idx[1]), max(tri_idx[0], tri_idx[1])), point_a.euclidean_distance(&point_b)), + (Edge(min(tri_idx[1], tri_idx[2]), max(tri_idx[1], tri_idx[2])), point_b.euclidean_distance(&point_c)), + (Edge(min(tri_idx[2], tri_idx[0]), max(tri_idx[2], tri_idx[0])), point_c.euclidean_distance(&point_a)), + ].to_vec().into_iter().sorted_by(|a, b| b.1.partial_cmp(&a.1).unwrap()).collect()) + } else { + None + }; + + let terminal_edge = edges_with_lengths.as_ref().map(|edges| edges[0].0); + let area = 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 + }; + + let node_connections: Option> = if types == 0 || types == 1 { + Some(tri_idx.iter().cloned().collect()) + } else { + None + }; + self.triangles.push(TriangleData { index, area, terminal_edge, - edges_with_lengths: edges_sorted, - node_connections, + edges_with_lengths: edges_with_lengths.clone(), + node_connections: node_connections.clone(), }); - - for &(edge, length) in &edges_with_lengths { - self.edge_lengths.insert(edge, length); - self.edge_to_triangles.entry(edge).or_default().push(index); + + if let Some(edges) = &edges_with_lengths { + for &(edge, length) in edges { + self.edge_lengths.insert(edge, length); + self.edge_to_triangles.entry(edge).or_default().push(index); + } } - - for &vertex in tri_idx { - self.vertex_to_triangles.entry(vertex).or_default().push(index); + + if let Some(nodes) = &node_connections { + for &vertex in nodes { + self.vertex_to_triangles.entry(vertex).or_default().push(index); + } } } + } pub fn random_points(center: (f32, f32), radius: f32, num_points: usize) -> Vec> { @@ -114,14 +127,14 @@ pub fn delaunay(points: &Vec>) -> Vec { result.triangles } -fn preprocess(points: &[Point], triangles: &[usize]) -> GeometryData { +fn preprocess(points: &[Point], triangles: &[usize], types: usize) -> GeometryData { let geometry_data: Arc> = Arc::new(Mutex::new(GeometryData::new())); triangles.par_chunks(6).enumerate().for_each(|(index, tri_idx)| { let points_clone: Vec> = points.to_vec(); // Clone points to avoid borrowing issues let gd: Arc> = geometry_data.clone(); // Clone Arc for use in each thread - gd.lock().unwrap().add_triangle(index, &points_clone, tri_idx); + gd.lock().unwrap().add_triangle(index, &points_clone, tri_idx, types); }); // Extract the GeometryData from the Arc>. This is safe to do here because @@ -147,21 +160,29 @@ fn delfin( // Sort all triangles by the longest terminal edge let triangles_sorted: Vec<(usize, f32)> = geometry_data.triangles.iter() - .map(|triangle_data: &TriangleData| { - let terminal_edge_length: f32 = geometry_data.edge_lengths[&triangle_data.terminal_edge]; - (triangle_data.index, terminal_edge_length) + .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(); - // Calculate densities based on reverse area + // Calculate areas for triangles that have an area calculated let areas: Vec = geometry_data.triangles.par_iter() - .map(|triangle_data: &TriangleData| triangle_data.area) - .collect(); + .filter_map(|triangle_data| triangle_data.area) + .map(|area| area) + .collect(); // Calculate mean and standard deviation of terminal edges lengths let terminal_edge_lengths: Vec = geometry_data.triangles.par_iter() - .map(|triangle_data: &TriangleData| geometry_data.edge_lengths[&triangle_data.terminal_edge]) + .filter_map(|triangle_data| { + triangle_data.terminal_edge.and_then(|edge| geometry_data.edge_lengths.get(&edge)) + }) + .cloned() .collect(); let (mean_terminal_edge, std_terminal_edge) = mean_std(terminal_edge_lengths); @@ -185,58 +206,69 @@ fn delfin( // Retrieve triangles that share the terminal edge, continue if less than 2 triangles share it let triangle_data: &TriangleData = &geometry_data.triangles[triangle_index]; - let terminal_edge: Edge = triangle_data.terminal_edge; - if let Some(connected_triangles) = geometry_data.edge_to_triangles.get(&terminal_edge) { - if connected_triangles.len() < 2 { - continue; - } + 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; + } - // 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() { - triangles_to_expand.remove(¤t_idx); - - // For each triangle, check its edges against the edges of the neighbors - for &neighbor_idx in connected_triangles { - if triangle_set.contains(&neighbor_idx) || processed_triangles.contains(&neighbor_idx) { - continue; - } - - let neighbor_data = &geometry_data.triangles[neighbor_idx]; - // Check if neighbor shares a terminal edge - if neighbor_data.terminal_edge == terminal_edge { - triangle_set.insert(neighbor_idx); - processed_triangles.insert(neighbor_idx); - triangles_to_expand.insert(neighbor_idx); + // 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) = 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; + // 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 + // 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)) - .map(|triangle_data: &TriangleData| triangle_data.area) + .filter_map(|&idx| geometry_data.triangles.get(idx).and_then(|td| td.area)) .sum(); - - // Calculate the area Z-score - let area_z_score: f32 = (total_area - mean_area) / std_area; - - // Filter based on the area Z-score and the minimum number of triangles + + // Calculate the area Z-score if std_area is non-zero to avoid division by zero. + let area_z_score: f32 = if std_area != 0.0 { + (total_area - mean_area) / std_area + } else { + -5.0 + }; + + // Filter based on the area Z-score and the minimum number of triangles. area_z_score >= min_area && poly_set.len() >= 3 }); @@ -248,9 +280,9 @@ fn main() { let triangles_indices: Vec = delaunay(&points); // Preprocess to create GeometryData - let geometry_data: GeometryData = preprocess(&points, &triangles_indices); + let geometry_data: GeometryData = preprocess(&points, &triangles_indices, 0); - // Define minimum voidness and minimum distance for delfin function + // Define minimum area and minimum distance for delfin function let min_area: f32 = 4.0; // Example threshold for voidness let min_distance: f32 = 1.0; // Example threshold for minimum distance (Z-score)