diff --git a/Cargo.lock b/Cargo.lock
index 7f37377..be479e5 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -20,6 +20,12 @@ version = "0.2.16"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "0942ffc6dcaadf03badf6e6a2d0228460359d5e34b57ccdc720b7382dfbd5ec5"
 
+[[package]]
+name = "anyhow"
+version = "1.0.81"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "0952808a6c2afd1aa8947271f3a60f1a6763c7b912d210184c5149b5cf147247"
+
 [[package]]
 name = "approx"
 version = "0.5.1"
@@ -474,6 +480,16 @@ dependencies = [
  "wide",
 ]
 
+[[package]]
+name = "simple_delaunay_lib"
+version = "0.2.0"
+source = "git+https://github.com/ibujah/simple_delaunay_lib#2812ea85f936723bb9ae3a1cd131589738b41821"
+dependencies = [
+ "anyhow",
+ "log",
+ "robust 1.1.0",
+]
+
 [[package]]
 name = "smallvec"
 version = "1.13.2"
@@ -569,7 +585,7 @@ dependencies = [
 
 [[package]]
 name = "xenobalanus"
-version = "0.1.3"
+version = "0.1.4"
 dependencies = [
  "delaunator",
  "geo",
@@ -577,6 +593,7 @@ dependencies = [
  "itertools 0.12.1",
  "rand",
  "rayon",
+ "simple_delaunay_lib",
  "statrs",
 ]
 
diff --git a/Cargo.toml b/Cargo.toml
index 861b613..19cfd20 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -1,11 +1,15 @@
 [package]
 name = "xenobalanus"
-version = "0.1.3"
+version = "0.1.4"
 edition = "2021"
 
 [lib]
 name = "xenobalanus"
 
+[features]
+default = ["three-d"]
+three-d = ["dep:simple_delaunay_lib"]
+
 [dependencies]
 delaunator = "1.0.2"
 geo = "0.28.0"
@@ -14,3 +18,4 @@ itertools = "0.12.1"
 rand = "0.8.5"
 rayon = "1.9.0"
 statrs = "0.16.0"
+simple_delaunay_lib = { version = "0.2.0", git = "https://github.com/ibujah/simple_delaunay_lib", optional = true }
\ No newline at end of file
diff --git a/src/lib.rs b/src/lib.rs
index 89854b0..5873ad8 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -6,6 +6,10 @@ use std::cmp::{min, max};
 use std::collections::{HashMap, HashSet};
 use std::sync::{Arc, Mutex};
 
+#[cfg(feature = "three-d")]
+use simple_delaunay_lib::delaunay_3d::delaunay_struct_3d::DelaunayStructure3D;
+mod threedee;
+
 #[derive(Debug, Clone, Copy)]
 pub struct Point {
     pub x: f32,
@@ -40,6 +44,12 @@ impl Point {
         bearing
     }
 
+    pub fn bearing_rad(&self, point: Point) -> f32 {
+        let delta_x = point.x - self.x;
+        let delta_y = point.y - self.y;
+        delta_y.atan2(delta_x)
+    }
+
 }
 
 impl From<Point> for Coord<f32> {
@@ -166,6 +176,9 @@ pub struct Xenobalanus {
     geometry_data: GeometryData,
     points: Vec<Point>,
     triangulation: Vec<usize>,
+    #[cfg(feature = "three-d")]
+    nodes: Vec<crate::threedee::Point3D>,
+    tetrahedra: DelaunayStructure3D,
 }
 
 impl Xenobalanus {
@@ -174,6 +187,9 @@ impl Xenobalanus {
             geometry_data: GeometryData::new(),
             points: Vec::new(),
             triangulation: Vec::new(),
+            #[cfg(feature = "three-d")]
+            nodes: Vec::new(),
+            tetrahedra: DelaunayStructure3D::new(),
         }
     }
 
diff --git a/src/threedee.rs b/src/threedee.rs
new file mode 100644
index 0000000..14151a2
--- /dev/null
+++ b/src/threedee.rs
@@ -0,0 +1,255 @@
+use rand::Rng;
+use simple_delaunay_lib::delaunay_3d::delaunay_struct_3d::ExtendedTetrahedron;
+use simple_delaunay_lib::delaunay_3d::simplicial_struct_3d::Node;
+use std::cmp::{min, max};
+use std::collections::HashSet;
+
+use super::{Edge, Xenobalanus};
+
+#[derive(Debug, Clone, Copy)]
+pub struct Point3D {
+    pub x: f32,
+    pub y: f32,
+    pub z: f32,
+}
+
+impl Point3D {
+
+    pub fn new(coords: (f32, f32, f32)) -> Point3D {
+        Point3D {x: coords.0, y: coords.1, z: coords.2 }
+    }
+
+    pub fn arr(&self) -> [f64; 3] {
+        [self.x as f64, self.y as f64, self.z as f64]
+    }
+ 
+    pub fn distance(&self,point: &Point3D) -> f32 {
+        ( 
+            (point.x - &self.x).powi(2) + 
+            (point.y - &self.y).powi(2) + 
+            (point.z - &self.z).powi(2) 
+        ).sqrt()
+    }
+}
+
+impl Xenobalanus {
+
+    pub fn point_3d(&self, index: usize) -> Point3D {
+        self.nodes[index]
+    }
+
+    pub fn points_3d(&self) -> Vec<(f32, f32, f32)> {
+        self.nodes.iter()
+            .map(|point| (point.x, point.y, point.z))
+            .collect()
+    }
+
+    pub fn tetrahedra(&self) -> Vec<[[f32; 3]; 2]> {
+        let total = self.tetrahedra.get_simplicial().get_nb_tetrahedra();
+        let mut edges: Vec<[[f32; 3]; 2]> = vec![];
+        for idx in 0..total {
+            match self.tetrahedra.get_extended_tetrahedron(idx) {
+                Ok(tet) => {
+                    match tet {
+                        ExtendedTetrahedron::Tetrahedron(te) => {
+                            edges.push([
+                                [te[0][0] as f32, te[0][1] as f32, te[0][2] as f32],
+                                [te[1][0] as f32, te[1][1] as f32, te[1][2] as f32]
+                                ]);
+                            edges.push([
+                                [te[1][0] as f32, te[1][1] as f32, te[1][2] as f32],
+                                [te[2][0] as f32, te[2][1] as f32, te[2][2] as f32]
+                                ]);
+                            edges.push([
+                                [te[2][0] as f32, te[2][1] as f32, te[2][2] as f32],
+                                [te[0][0] as f32, te[0][1] as f32, te[0][2] as f32]
+                                ]);
+                            edges.push([
+                                [te[3][0] as f32, te[3][1] as f32, te[3][2] as f32],
+                                [te[0][0] as f32, te[0][1] as f32, te[0][2] as f32]
+                                ]);
+                            edges.push([
+                                [te[3][0] as f32, te[3][1] as f32, te[3][2] as f32],
+                                [te[1][0] as f32, te[1][1] as f32, te[1][2] as f32]
+                                ]);
+                            edges.push([
+                                [te[3][0] as f32, te[3][1] as f32, te[3][2] as f32],
+                                [te[2][0] as f32, te[2][1] as f32, te[2][2] as f32]
+                                ]);
+                        },
+                        ExtendedTetrahedron::Triangle(te) => {
+                            edges.push([
+                                [te[0][0] as f32, te[0][1] as f32, te[0][2] as f32],
+                                [te[1][0] as f32, te[1][1] as f32, te[1][2] as f32]
+                                ]);
+                            edges.push([
+                                [te[1][0] as f32, te[1][1] as f32, te[1][2] as f32],
+                                [te[2][0] as f32, te[2][1] as f32, te[2][2] as f32]
+                                ]);
+                            edges.push([
+                                [te[2][0] as f32, te[2][1] as f32, te[2][2] as f32],
+                                [te[0][0] as f32, te[0][1] as f32, te[0][2] as f32]
+                                ]);
+                        }
+                    }
+                },
+                Err(_) => {
+                    println!("Error getting tetrahedron at index {}", idx);
+                    continue;
+                }
+            }
+        }
+        edges
+    }
+
+    pub fn random_points_3d(&mut self, center: (f32, f32, f32), side_length: f32, num_points: u32) {
+        
+        // generate random points in a cube
+        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 min_z = center.2 - side_length / 2.0;
+        let max_z = center.2 + 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);
+            let z: f32 = min_y + rng.gen_range(0.0..=1.0) as f32 * ( max_z - min_z);
+            self.nodes.push(Point3D {x, y, z});
+        }
+    }
+
+    pub fn delaunay_3d(&mut self) {
+
+        let vertices: Vec<[f64; 3]> = self.nodes.iter().map(|node| {
+            node.arr()
+        }).collect();
+        self.tetrahedra.insert_vertices(&vertices, true).unwrap_or_default();
+    }
+
+    pub fn add_triangle(&mut self, vertex1: usize, vertex2: usize, vertex3: usize) {
+
+        let p1 = self.tetrahedra.get_vertices()[vertex1];
+        let p2 = self.tetrahedra.get_vertices()[vertex2];
+        let p3 = self.tetrahedra.get_vertices()[vertex3];
+
+        let n1 = Point3D::new((p1[0] as f32, p1[1] as f32, p1[2] as f32));
+        let n2 = Point3D::new((p2[0] as f32, p2[1] as f32, p2[2] as f32));
+        let n3 = Point3D::new((p3[0] as f32, p3[1] as f32, p3[2] as f32));
+
+        let l12 = n1.distance(&n2);
+        let l23 = n2.distance(&n3);
+        let l31 = n3.distance(&n1);
+
+        let edge12 = Edge(min(vertex1, vertex2), max(vertex1, vertex2));
+        let edge23 = Edge(min(vertex2, vertex3), max(vertex2, vertex3));
+        let edge31 = Edge(min(vertex3, vertex1), max(vertex3, vertex1));
+
+        self.geometry_data.edge_lengths.insert(edge12, l12);
+        self.geometry_data.edge_lengths.insert(edge23, l23);
+        self.geometry_data.edge_lengths.insert(edge31, l31);
+
+        self.geometry_data.vertex_connections.entry(vertex1).or_insert_with(HashSet::new).insert(vertex2);
+        self.geometry_data.vertex_connections.entry(vertex2).or_insert_with(HashSet::new).insert(vertex1);
+        
+        self.geometry_data.vertex_connections.entry(vertex2).or_insert_with(HashSet::new).insert(vertex3);
+        self.geometry_data.vertex_connections.entry(vertex3).or_insert_with(HashSet::new).insert(vertex2);
+        
+        self.geometry_data.vertex_connections.entry(vertex3).or_insert_with(HashSet::new).insert(vertex1);
+        self.geometry_data.vertex_connections.entry(vertex1).or_insert_with(HashSet::new).insert(vertex3);
+
+    }
+
+    pub fn add_tetrahedron(&mut self, vertex1: usize, vertex2: usize, vertex3: usize, vertex4: usize) {
+
+        let p1 = self.tetrahedra.get_vertices()[vertex1];
+        let p2 = self.tetrahedra.get_vertices()[vertex2];
+        let p3 = self.tetrahedra.get_vertices()[vertex3];
+        let p4 = self.tetrahedra.get_vertices()[vertex4];
+
+        let n1 = Point3D::new((p1[0] as f32, p1[1] as f32, p1[2] as f32));
+        let n2 = Point3D::new((p2[0] as f32, p2[1] as f32, p2[2] as f32));
+        let n3 = Point3D::new((p3[0] as f32, p3[1] as f32, p3[2] as f32));
+        let n4 = Point3D::new((p4[0] as f32, p4[1] as f32, p4[2] as f32));
+
+        let l12 = &n1.distance(&n2);
+        let l23 = &n2.distance(&n3);
+        let l31 = &n3.distance(&n1);
+        let l41 = &n4.distance(&n1);
+        let l42 = &n4.distance(&n2);
+        let l43 = &n4.distance(&n3);
+
+        let edge12 = Edge(min(vertex1, vertex2), max(vertex1, vertex2));
+        let edge23 = Edge(min(vertex2, vertex3), max(vertex2, vertex3));
+        let edge31 = Edge(min(vertex3, vertex1), max(vertex3, vertex1));
+        let edge41 = Edge(min(vertex4, vertex1), max(vertex4, vertex1));
+        let edge42 = Edge(min(vertex4, vertex2), max(vertex4, vertex2));
+        let edge43 = Edge(min(vertex4, vertex3), max(vertex4, vertex3));
+
+        self.geometry_data.edge_lengths.insert(edge12, *l12);
+        self.geometry_data.edge_lengths.insert(edge23, *l23);
+        self.geometry_data.edge_lengths.insert(edge31, *l31);
+        self.geometry_data.edge_lengths.insert(edge41, *l41);
+        self.geometry_data.edge_lengths.insert(edge42, *l42);
+        self.geometry_data.edge_lengths.insert(edge43, *l43);
+
+        self.geometry_data.vertex_connections.entry(vertex1).or_insert_with(HashSet::new).insert(vertex2);
+        self.geometry_data.vertex_connections.entry(vertex2).or_insert_with(HashSet::new).insert(vertex1);
+        
+        self.geometry_data.vertex_connections.entry(vertex2).or_insert_with(HashSet::new).insert(vertex3);
+        self.geometry_data.vertex_connections.entry(vertex3).or_insert_with(HashSet::new).insert(vertex2);
+        
+        self.geometry_data.vertex_connections.entry(vertex3).or_insert_with(HashSet::new).insert(vertex1);
+        self.geometry_data.vertex_connections.entry(vertex1).or_insert_with(HashSet::new).insert(vertex3);
+        
+        self.geometry_data.vertex_connections.entry(vertex4).or_insert_with(HashSet::new).insert(vertex1);
+        self.geometry_data.vertex_connections.entry(vertex1).or_insert_with(HashSet::new).insert(vertex4);
+        
+        self.geometry_data.vertex_connections.entry(vertex4).or_insert_with(HashSet::new).insert(vertex2);
+        self.geometry_data.vertex_connections.entry(vertex2).or_insert_with(HashSet::new).insert(vertex4);
+        
+        self.geometry_data.vertex_connections.entry(vertex4).or_insert_with(HashSet::new).insert(vertex3);
+        self.geometry_data.vertex_connections.entry(vertex3).or_insert_with(HashSet::new).insert(vertex4);
+
+    }
+
+    pub fn preprocess_3d(&mut self) {
+        let structure = self.tetrahedra.get_simplicial();
+        let num_tetras = structure.get_nb_tetrahedra();
+    
+        for tetra_idx in 0..num_tetras {
+            let tetrahedron = match self.tetrahedra.get_simplicial().get_tetrahedron(tetra_idx) {
+                Ok(tetra) => tetra.nodes(),
+                Err(_) => {
+                    println!("Error getting tetrahedron at index {}", tetra_idx);
+                    continue; // Skip this iteration and proceed with the next one
+                }
+            };
+    
+            let [node1, node2, node3, node4] = tetrahedron;
+            match (node1, node2, node3, node4) {
+                (Node::Infinity, Node::Value(ind_v2), Node::Value(ind_v3), Node::Value(ind_v4)) => {
+                    self.add_triangle(ind_v2, ind_v3, ind_v4);
+                },
+                (Node::Value(ind_v1), Node::Infinity, Node::Value(ind_v3), Node::Value(ind_v4)) => {
+                    self.add_triangle(ind_v3, ind_v4, ind_v1);
+                },
+                (Node::Value(ind_v1), Node::Value(ind_v2), Node::Infinity, Node::Value(ind_v4)) => {
+                    self.add_triangle(ind_v4, ind_v1, ind_v2);
+                },
+                (Node::Value(ind_v1), Node::Value(ind_v2), Node::Value(ind_v3), Node::Infinity) => {
+                    self.add_triangle(ind_v1, ind_v2, ind_v3);
+                },
+                (Node::Value(ind_v1), Node::Value(ind_v2), Node::Value(ind_v3), Node::Value(ind_v4)) => {
+                    self.add_tetrahedron(ind_v1, ind_v2, ind_v3, ind_v4);
+                },
+                (_, _, _, _) => {
+                    println!("Encountered an unexpected pattern of nodes");
+                },
+            };
+        }
+    }
+    
+
+}
\ No newline at end of file
diff --git a/tests/threed.rs b/tests/threed.rs
new file mode 100644
index 0000000..aabf4ba
--- /dev/null
+++ b/tests/threed.rs
@@ -0,0 +1,12 @@
+use xenobalanus::Xenobalanus;
+
+#[test]
+fn threed() {
+    let mut geodesic = Xenobalanus::new();
+    geodesic.random_points_3d((0.0, 0.0, 0.0), 1.0, 1000);
+    geodesic.delaunay_3d();
+    geodesic.preprocess_3d();
+    let attractors = geodesic.dtscan(5, 0.085);
+
+    println!("{:#?}", attractors);
+}
\ No newline at end of file