diff --git a/Cargo.lock b/Cargo.lock
index 893b4c3..28abaf7 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -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",
 ]
diff --git a/Cargo.toml b/Cargo.toml
index 69e3a2b..41f371e 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -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"
diff --git a/src/main.rs b/src/main.rs
index 28ef627..f8e996f 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -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());
 }