rustic data structure

This commit is contained in:
randogoth 2024-03-17 19:18:53 +02:00
parent 486f3246fc
commit 13d53c36e3
3 changed files with 401 additions and 67 deletions

181
Cargo.lock generated
View file

@ -35,6 +35,12 @@ version = "1.1.0"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "d468802bab17cbc0cc575e9b053f41e72aa36bfa6b7f55e3529ffa43161b97fa"
[[package]]
name = "bytemuck"
version = "1.15.0"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "5d6d68c57235a3a081186990eca2867354726650f42f7516ca50c28d6281fd15"
[[package]]
name = "byteorder"
version = "1.5.0"
@ -87,7 +93,7 @@ version = "0.4.3"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "79127ed59a85d7687c409e9978547cffb7dc79675355ed22da6b66fd5f6ead01"
dependencies = [
"itertools",
"itertools 0.11.0",
"num-traits",
]
@ -190,6 +196,21 @@ dependencies = [
"either",
]
[[package]]
name = "itertools"
version = "0.12.1"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "ba291022dbbd398a455acf126c1e341954079855bc60dfdda641363bd6922569"
dependencies = [
"either",
]
[[package]]
name = "lazy_static"
version = "1.4.0"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "e2abad23fbc42b3700f2f279844dc832adb2b2eb069b2df918f455c4e18cc646"
[[package]]
name = "libc"
version = "0.2.153"
@ -208,6 +229,74 @@ version = "0.4.21"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "90ed8c1e510134f979dbc4f070f87d4313098b704861a105fe34231c70a3901c"
[[package]]
name = "matrixmultiply"
version = "0.3.8"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "7574c1cf36da4798ab73da5b215bbf444f50718207754cb522201d78d1cd0ff2"
dependencies = [
"autocfg",
"rawpointer",
]
[[package]]
name = "nalgebra"
version = "0.29.0"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "d506eb7e08d6329505faa8a3a00a5dcc6de9f76e0c77e4b75763ae3c770831ff"
dependencies = [
"approx",
"matrixmultiply",
"nalgebra-macros",
"num-complex",
"num-rational",
"num-traits",
"rand",
"rand_distr",
"simba",
"typenum",
]
[[package]]
name = "nalgebra-macros"
version = "0.1.0"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "01fcc0b8149b4632adc89ac3b7b31a12fb6099a0317a4eb2ebff574ef7de7218"
dependencies = [
"proc-macro2",
"quote",
"syn 1.0.109",
]
[[package]]
name = "num-complex"
version = "0.4.5"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "23c6602fda94a57c990fe0df199a035d83576b496aa29f4e634a8ac6004e68a6"
dependencies = [
"num-traits",
]
[[package]]
name = "num-integer"
version = "0.1.46"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "7969661fd2958a5cb096e56c8e1ad0444ac2bbcd0061bd28660485a44879858f"
dependencies = [
"num-traits",
]
[[package]]
name = "num-rational"
version = "0.4.1"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "0638a1c9d0a3c0914158145bc76cff373a75a627e6ecbfb71cbe6f453a5a19b0"
dependencies = [
"autocfg",
"num-integer",
"num-traits",
]
[[package]]
name = "num-traits"
version = "0.2.18"
@ -224,6 +313,12 @@ version = "1.19.0"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "3fdb12b2476b595f9358c5161aa467c2438859caa136dec86c26fdd2efe17b92"
[[package]]
name = "paste"
version = "1.0.14"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "de3145af08024dea9fa9914f381a17b8fc6034dfb00f3a84013f7ff43f29ed4c"
[[package]]
name = "ppv-lite86"
version = "0.2.17"
@ -278,6 +373,22 @@ dependencies = [
"getrandom",
]
[[package]]
name = "rand_distr"
version = "0.4.3"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "32cb0b9bc82b0a0876c2dd994a7e7a2683d3e7390ca40e6886785ef0c7e3ee31"
dependencies = [
"num-traits",
"rand",
]
[[package]]
name = "rawpointer"
version = "0.2.1"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "60a357793950651c4ed0f3f52338f53b2f809f32d83a07f72909fa13e4c6c1e3"
[[package]]
name = "rayon"
version = "1.9.0"
@ -321,6 +432,15 @@ dependencies = [
"smallvec",
]
[[package]]
name = "safe_arch"
version = "0.7.1"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "f398075ce1e6a179b46f51bd88d0598b92b00d3551f1a2d4ac49e771b56ac354"
dependencies = [
"bytemuck",
]
[[package]]
name = "serde"
version = "1.0.197"
@ -338,7 +458,20 @@ checksum = "7eb0b34b42edc17f6b7cac84a52a1c5f0e1bb2227e997ca9011ea3dd34e8610b"
dependencies = [
"proc-macro2",
"quote",
"syn",
"syn 2.0.53",
]
[[package]]
name = "simba"
version = "0.6.0"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "f0b7840f121a46d63066ee7a99fc81dcabbc6105e437cae43528cea199b5a05f"
dependencies = [
"approx",
"num-complex",
"num-traits",
"paste",
"wide",
]
[[package]]
@ -365,6 +498,30 @@ version = "1.2.0"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "a8f112729512f8e442d81f95a8a7ddf2b7c6b8a1a6f509a95864142b30cab2d3"
[[package]]
name = "statrs"
version = "0.16.0"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "2d08e5e1748192713cc281da8b16924fb46be7b0c2431854eadc785823e5696e"
dependencies = [
"approx",
"lazy_static",
"nalgebra",
"num-traits",
"rand",
]
[[package]]
name = "syn"
version = "1.0.109"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "72b64191b275b66ffe2469e8af2c1cfe3bafa67b529ead792a6d0160888b4237"
dependencies = [
"proc-macro2",
"quote",
"unicode-ident",
]
[[package]]
name = "syn"
version = "2.0.53"
@ -382,10 +539,18 @@ version = "0.1.0"
dependencies = [
"delaunator",
"geo",
"itertools 0.12.1",
"rand",
"rayon",
"statrs",
]
[[package]]
name = "typenum"
version = "1.17.0"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "42ff0bf0c66b8238c6f3b578df37d0b7848e55df8577b3f74f92a69acceeb825"
[[package]]
name = "unicode-ident"
version = "1.0.12"
@ -404,6 +569,16 @@ version = "0.11.0+wasi-snapshot-preview1"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "9c8d87e72b64a3b4db28d11ce29237c246188f4f51057d65a7eab63b7987e423"
[[package]]
name = "wide"
version = "0.7.15"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "89beec544f246e679fc25490e3f8e08003bc4bf612068f325120dad4cea02c1c"
dependencies = [
"bytemuck",
"safe_arch",
]
[[package]]
name = "zerocopy"
version = "0.7.32"
@ -421,5 +596,5 @@ checksum = "9ce1b18ccd8e73a9321186f97e46f9f04b778851177567b1975109d26a08d2a6"
dependencies = [
"proc-macro2",
"quote",
"syn",
"syn 2.0.53",
]

View file

@ -8,5 +8,7 @@ edition = "2021"
[dependencies]
delaunator = "1.0.2"
geo = "0.28.0"
itertools = "0.12.1"
rand = "0.8.5"
rayon = "1.9.0"
statrs = "0.16.0"

View file

@ -4,6 +4,83 @@ use std::cmp::{min, max};
use rand::Rng;
use rayon::prelude::*;
use delaunator::{triangulate, Point as DelaunatorPoint};
use itertools::Itertools;
use std::sync::{Arc, Mutex};
#[derive(Debug, Clone, Copy, Hash, PartialEq, Eq)]
struct Edge(usize, usize);
#[derive(Debug)]
struct TriangleData {
index: usize,
area: f32,
terminal_edge: Edge,
edges_with_lengths: Vec<(Edge, f32)>,
node_connections: HashSet<usize>,
}
#[derive(Debug)]
struct GeometryData {
triangles: Vec<TriangleData>,
edge_to_triangles: HashMap<Edge, Vec<usize>>, // Maps an edge to triangle indices
edge_lengths: HashMap<Edge, f32>, // Edge lengths
vertex_to_triangles: HashMap<usize, Vec<usize>>, // Maps a vertex to connected triangle indices
}
impl GeometryData {
fn new() -> Self {
GeometryData {
triangles: Vec::new(),
edge_to_triangles: HashMap::new(),
edge_lengths: HashMap::new(),
vertex_to_triangles: HashMap::new(),
}
}
// Function to add a triangle to the GeometryData
fn add_triangle(&mut self, index: usize, points: &[Point<f32>], tri_idx: &[usize]) {
let point_a: Point<f32> = points[tri_idx[0]];
let point_b: Point<f32> = points[tri_idx[1]];
let point_c: Point<f32> = 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<usize> = tri_idx.iter().cloned().collect::<HashSet<_>>();
self.triangles.push(TriangleData {
index,
area,
terminal_edge,
edges_with_lengths: edges_sorted,
node_connections,
});
for &(edge, length) in &edges_with_lengths {
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);
}
}
}
pub fn random_points(center: (f32, f32), radius: f32, num_points: usize) -> Vec<Point<f32>> {
let mut rng: rand::prelude::ThreadRng = rand::thread_rng();
@ -37,80 +114,160 @@ pub fn delaunay(points: &Vec<Point<f32>>) -> Vec<usize> {
result.triangles
}
fn preprocess(points: &[Point<f32>], triangles: &[usize]) -> (
HashMap<usize, HashSet<usize>>,
HashMap<usize, f32>,
HashMap<(usize, usize), HashSet<usize>>,
HashMap<(usize, usize), f32>,
HashMap<usize, HashSet<usize>>,
) {
// Process each set of triangle indices in parallel
let triangle_calculation: Vec<_> = triangles.par_chunks(3).map(|tri_idx| {
let point_a: Point<f32> = points[tri_idx[0]];
let point_b: Point<f32> = points[tri_idx[1]];
let point_c: Point<f32> = points[tri_idx[2]];
fn preprocess(points: &[Point<f32>], triangles: &[usize]) -> GeometryData {
let geometry_data: Arc<Mutex<GeometryData>> = Arc::new(Mutex::new(GeometryData::new()));
// Calculate the lengths of each edge and pair them with their vertex indices
let edges_with_lengths: [((usize, usize), f32); 3] = [
((min(tri_idx[0], tri_idx[1]), max(tri_idx[0], tri_idx[1])), point_a.euclidean_distance(&point_b)),
((min(tri_idx[1], tri_idx[2]), max(tri_idx[1], tri_idx[2])), point_b.euclidean_distance(&point_c)),
((min(tri_idx[2], tri_idx[0]), max(tri_idx[2], tri_idx[0])), point_c.euclidean_distance(&point_a)),
];
triangles.par_chunks(3).enumerate().for_each(|(index, tri_idx)| {
let points_clone: Vec<Point<f32>> = points.to_vec(); // Clone points to avoid borrowing issues
let gd: Arc<Mutex<GeometryData>> = geometry_data.clone(); // Clone Arc for use in each thread
// Sort edges by length to ensure the longest edge is first
let mut edges_sorted: Vec<((usize, usize), f32)> = edges_with_lengths.to_vec();
edges_sorted.sort_by(|a: &((usize, usize), f32), b: &((usize, usize), f32)| b.1.partial_cmp(&a.1).unwrap());
let terminal_edges: HashSet<usize> = [edges_sorted[0].0 .0, edges_sorted[0].0 .1].iter().cloned().collect::<HashSet<_>>();
gd.lock().unwrap().add_triangle(index, &points_clone, tri_idx);
});
// Collect 'area_map' with area for each triangle
let poly: Polygon<f32> = 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![]);
let area: f32 = poly.unsigned_area();
// Extract the GeometryData from the Arc<Mutex<>>. This is safe to do here because
// the par_iter has completed, and we know no other threads are accessing it.
Arc::try_unwrap(geometry_data).unwrap().into_inner().unwrap()
}
// Generate the node connections
let mut node_connections: HashMap<usize, HashSet<usize>> = HashMap::new();
for &idx in tri_idx.iter() {
let connected_nodes: HashSet<usize> = tri_idx.iter().filter(|&&x| x != idx).cloned().collect::<HashSet<_>>();
node_connections.insert(idx, connected_nodes);
fn mean_std(dataset: Vec<f32>) -> (f32, f32) {
let mean: f32 = dataset.iter().sum::<f32>() / dataset.len() as f32;
let std: f32 = (dataset.iter().map(|&length| {
let diff = length - mean;
diff * diff}
).sum::<f32>() / dataset.len() as f32).sqrt();
(mean, std)
}
fn delfin(
geometry_data: &GeometryData,
min_voidness: f32,
min_distance: f32,
) -> Vec<HashSet<usize>> {
// Sort all triangles by the longest terminal edge
let triangles_sorted: Vec<usize> = 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)
})
.sorted_by(|a, b| b.1.partial_cmp(&a.1).unwrap()) // Sort in descending order by edge length
.map(|(index, _)| index) // Extract triangle indices
.collect();
// Calculate densities based on reverse area
let densities: Vec<f32> = geometry_data.triangles.iter()
.map(|triangle_data: &TriangleData| 1.0 / triangle_data.area)
.collect();
// Calculate mean and standard deviation of terminal edges lengths
let terminal_edge_lengths: Vec<f32> = geometry_data.triangles.iter()
.map(|triangle_data: &TriangleData| geometry_data.edge_lengths[&triangle_data.terminal_edge])
.collect();
let (mean_terminal_edge, std_terminal_edge) = mean_std(terminal_edge_lengths);
let (mean_density, std_density) = mean_std(densities);
let mut void_polygons: Vec<HashSet<usize>> = Vec::new();
let mut processed_triangles: HashSet<usize> = HashSet::new();
for &triangle_index in &triangles_sorted {
// Skip if this triangle has already been processed
if processed_triangles.contains(&triangle_index) {
continue;
}
// Return all calculated data for this triangle
(tri_idx[0] / 3, terminal_edges, area, edges_sorted, node_connections)
}).collect();
// Initialize shared data structures
let mut terminal_map: HashMap<usize, HashSet<usize>> = HashMap::new(); // terminal edge for each triangle
let mut area_map: HashMap<usize, f32> = HashMap::new(); // area for each triangle
let mut wing_map: HashMap<(usize, usize), HashSet<usize>> = HashMap::new(); // triangle index for each edge
let mut edge_map: HashMap<(usize, usize), f32> = HashMap::new(); // length for each edge
let mut node_map: HashMap<usize, HashSet<usize>> = HashMap::new(); // vertices connected to each vertex
// Merge all triangle results
for (triangle_index, terminal_edges, area, edges_sorted, node_connections) in triangle_calculation {
terminal_map.insert(triangle_index, terminal_edges);
area_map.insert(triangle_index, area);
for &(edge, length) in &edges_sorted {
wing_map.entry(edge).or_insert_with(HashSet::new).insert(triangle_index);
edge_map.insert(edge, length);
// Retrieve the terminal edge for the current triangle
let triangle_data: &TriangleData = &geometry_data.triangles[triangle_index];
let terminal_edge: Edge = triangle_data.terminal_edge;
// Calculate the Z-score for the terminal edge length
let terminal_edge_length: f32 = geometry_data.edge_lengths[&terminal_edge];
let z_score: f32 = (terminal_edge_length - mean_terminal_edge) / std_terminal_edge;
// println!("Terminal Length: {}", terminal_edge_length);
// println!("Terminal Length Z-Score: {}", z_score);
// Continue if the Z-score is below the minimum distance threshold
if z_score < min_distance {
continue;
}
for (idx, connections) in node_connections {
node_map.entry(idx).or_insert_with(HashSet::new).extend(connections);
// Retrieve triangles that share the terminal edge, continue if less than 2 triangles share it
if let Some(connected_triangles) = geometry_data.edge_to_triangles.get(&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<usize> = 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<usize> = triangle_set.clone();
while let Some(current_idx) = triangles_to_expand.iter().next().cloned() {
triangles_to_expand.remove(&current_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);
}
}
}
// 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;
}
}
(terminal_map, area_map, wing_map, edge_map, node_map)
// Filter out void polygon sets
void_polygons.retain(|poly_set: &HashSet<usize>| {
// 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)
.sum();
// Calculate the polygon density
let polygon_density: f32 = if total_area > 0.0 { 1.0 / total_area } else { 0.0 };
// Calculate the density Z-score
let density_z_score: f32 = (polygon_density - mean_density) / std_density;
// println!("Density Z-Score: {}", density_z_score.abs());
// Filter based on the density Z-score and the minimum number of triangles
density_z_score.abs() >= min_voidness && poly_set.len() >= 3
});
return void_polygons;
}
fn main() {
let points: Vec<Point<f32>> = random_points((0.0, 0.0), 300.0, 2000);
let triangles: Vec<usize> = delaunay(&points);
let result: (HashMap<usize, HashSet<usize>>, HashMap<usize, f32>, HashMap<(usize, usize), HashSet<usize>>, HashMap<(usize, usize), f32>, HashMap<usize, HashSet<usize>>) = preprocess(&points, &triangles);
println!("{:?}", result);
let points: Vec<Point<f32>> = random_points((0.0, 0.0), 1000.0, 100000);
let triangles_indices: Vec<usize> = delaunay(&points);
// println!("{:?}", triangles_indices);
// Preprocess to create GeometryData
let geometry_data: GeometryData = preprocess(&points, &triangles_indices);
// Define minimum voidness and minimum distance for delfin function
let min_voidness: f32 = 0.2; // Example threshold for voidness
let min_distance: f32 = 0.0; // Example threshold for minimum distance (Z-score)
// Execute delfin function with the generated GeometryData
let void_polygons: Vec<HashSet<usize>> = delfin(&geometry_data, min_voidness, min_distance);
// To display the result, let's just print the count of void polygons found
println!("Void Polygons Found: {}", void_polygons.len());
}