diff --git a/Cargo.lock b/Cargo.lock index f051aa0..6d4db0d 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -139,13 +139,6 @@ dependencies = [ "libc", ] -[[package]] -name = "hilbert" -version = "0.1.0" -dependencies = [ - "image", -] - [[package]] name = "image" version = "0.23.14" @@ -324,6 +317,13 @@ dependencies = [ "weezl", ] +[[package]] +name = "voronoi" +version = "0.1.0" +dependencies = [ + "image", +] + [[package]] name = "weezl" version = "0.1.6" diff --git a/Cargo.toml b/Cargo.toml index beda60c..be0bf43 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,5 +1,5 @@ [package] -name = "hilbert" +name = "voronoi" version = "0.1.0" edition = "2021" diff --git a/src/main.rs b/src/main.rs index b34ac53..a03a0fa 100644 --- a/src/main.rs +++ b/src/main.rs @@ -27,7 +27,12 @@ impl Distance for Manhattan { } impl Distance for HilbertPoly { + fn dist(&self, a: [f64; N], b: [f64; N]) -> Result { + if a == b { + return Ok(0.0); + } + // ("#0: {:?} {:?}", a, b); // 2d Hilbert Convex Polygon, 3d is slightly different as it will be to a side (edge of // convex set ) @@ -57,6 +62,7 @@ impl Distance for HilbertPoly { } // Check if points are both inside the set + let (mut i, mut j, mut c) = (0, self.vertices.len() - 1, false); loop { if i >= self.vertices.len() { @@ -94,6 +100,37 @@ impl Distance for HilbertPoly { return Err("Second point not in set."); } + // Check if points are on the boundary + // + let c = edges.clone().into_iter().map(|e| [[e[0][0] - a[0], e[0][1] - a[1]], [e[1][0] - a[0], e[1][1] - a[1]]]).map(|s| { + let ab = f64::sqrt(s[0][0] * s[0][0] + s[0][1] * s[0][1]); + let bc = f64::sqrt(s[1][0] * s[1][0] + s[1][1] * s[1][1]); + let ac = f64::sqrt((s[1][0] - s[0][0]) * (s[1][0] - s[0][0]) + (s[1][1] - s[0][1]) * (s[1][1] - s[0][1])); + // ("#1: boundary check: {} {}", f64::round(ab + bc), f64::round(ac)); + if f64::round(ab + bc) != f64::round(ac) { + return true; + } + false + }).all(|x| x); + + if !c { + return Err("First point not in set."); + } + + let c = edges.clone().into_iter().map(|e| [[e[0][0] - b[0], e[0][1] - b[1]], [e[1][0] - b[0], e[1][1] - b[1]]]).map(|s| { + let ab = f64::sqrt(s[0][0] * s[0][0] + s[0][1] * s[0][1]); + let bc = f64::sqrt(s[1][0] * s[1][0] + s[1][1] * s[1][1]); + let ac = f64::sqrt((s[1][0] - s[0][0]) * (s[1][0] - s[0][0]) + (s[1][1] - s[0][1]) * (s[1][1] - s[0][1])); + if f64::round(ab + bc) != f64::round(ac) { + return true; + } + false + }).all(|x| x); + + if !c { + return Err("Second point not in set."); + } + // Intersect the line with each polygon side let slope = (b[1] - a[1]) / (b[0] - a[0]); let p_intersect = edges.into_iter().map(|p| { @@ -101,14 +138,17 @@ impl Distance for HilbertPoly { let d = edge_slope - slope; if d == 0.0 { None + } else if d == f64::INFINITY || d == f64::NEG_INFINITY { + Some([a[0], edge_slope * a[0] - edge_slope * p[0][0] + p[0][1]]) } else { + // ("#2: slopes: {:?} {:?} {} {} {:?}", a, b, slope, edge_slope, d); Some([((edge_slope * p[0][0] - p[0][1]) - (slope * a[0] - a[1])) / d, ((slope * (edge_slope * p[0][0] - p[0][1])) - (edge_slope * (slope * a[0] - a[1]))) / d]) } }).filter(|x| x.is_some()).map(|x| x.map_or([0.0, 0.0], |x| x)).collect::>(); // Find the closest distance let a_closest = p_intersect.clone().into_iter().map(|p| (p, (p[0] - a[0]) * (p[0] - a[0]) + (p[1] - a[1]) * (p[1] - a[1]))).min_by(|a, b| a.1.partial_cmp(&b.1).unwrap()).unwrap(); - let b_closest = p_intersect.into_iter().map(|p| (p, (p[0] - b[0]) * (p[0] - b[0]) + (p[1] - b[1]) * (p[1] - b[1]))).min_by(|a, b| a.1.partial_cmp(&b.1).unwrap()).unwrap(); + let b_closest = p_intersect.into_iter().map(|p| (p, (p[0] - b[0]) * (p[0] - b[0]) + (p[1] - b[1]) * (p[1] - b[1]))).filter(|p| p.0 != a_closest.0).min_by(|a, b| a.1.partial_cmp(&b.1).unwrap()).unwrap(); let pb = f64::sqrt((a_closest.0[0] - b[0]) * (a_closest.0[0] - b[0]) + (a_closest.0[1] - b[1]) * (a_closest.0[1] - b[1])); let pa = f64::sqrt(a_closest.1); @@ -119,23 +159,45 @@ impl Distance for HilbertPoly { } } -fn naive_voronoi_2d(metric: &mut dyn Distance<2>, w: u32, h: u32, points: Vec<(u32, u32)>) { +fn naive_voronoi_2d(metric: &mut dyn Distance<2>, w: u32, h: u32, points: Vec<(u32, u32, Rgb)>) -> RgbImage { let mut image: RgbImage = ImageBuffer::new(w, h); // All points must be in the width / height - for p in points { + for p in points.clone() { assert!(p.0 <= w); assert!(p.1 <= h); } + + for i in 0..w { for j in 0..h { - + let c = points.clone().into_iter().map(|p| metric.dist([i as f64, j as f64], [p.0 as f64, p.1 as f64]) /* ("#3: points for dist: ({} {}) ({} {}) {:?}", i, j, p.0, p.1, u)*/); + if c.clone().map(|d| d.is_ok()).all(|x| x) { + let closest = c.clone().map(|d| d.unwrap()).enumerate().min_by(|a, b| a.1.partial_cmp(&b.1).unwrap()).unwrap(); + if closest.1 == 0.0 { + *image.get_pixel_mut(i, j) = Rgb([0, 0, 0]); + } else { + *image.get_pixel_mut(i, j) = points[closest.0].2; + } + } else { + *image.get_pixel_mut(i, j) = Rgb([255, 255, 255]); + } } } + image } fn main() { - let triangle = HilbertPoly { vertices: vec![[-8.0, 0.0], [8.0, 0.0], [0.0, 8.0]] }; - println!("{}", triangle.dist([-2.0, 2.0], [2.0, 2.0]).unwrap()); + let mut triangle = HilbertPoly { vertices: vec![[300.0, 0.0], [0.0, 600.0], [600.0, 600.0]] }; + let mut euclidean: Euclidean<2> = Euclidean(); + let mut manhattan: Manhattan<2> = Manhattan(); + let points = vec![(360, 340, Rgb([255, 0, 0])), (340, 340, Rgb([0, 255, 0])), (400, 400, Rgb([0, 0, 255]))]; + let h = naive_voronoi_2d(&mut triangle, 600, 600, points.clone()); + h.save("hilbert.png").unwrap(); + let e = naive_voronoi_2d(&mut euclidean, 600, 600, points.clone()); + e.save("euclidean.png").unwrap(); + let m = naive_voronoi_2d(&mut manhattan, 600, 600, points); + m.save("manhattan.png").unwrap(); + }