diff --git a/benches/benchmark_utilities.rs b/benches/benchmark_utilities.rs index 23029a0..7cf2817 100644 --- a/benches/benchmark_utilities.rs +++ b/benches/benchmark_utilities.rs @@ -1,7 +1,7 @@ use core::fmt::{self, Display, Formatter}; use criterion::{measurement::WallTime, BenchmarkGroup, BenchmarkId, Throughput}; -use rand::{distributions::uniform::SampleUniform, Rng, SeedableRng}; +use rand::{distr::uniform::SampleUniform, Rng, SeedableRng}; use spade::{ DelaunayTriangulation, HierarchyHintGenerator, LastUsedVertexHintGenerator, Point2, SpadeNum, Triangulation, @@ -27,7 +27,7 @@ pub fn uniform_distribution( where S::Sampler: Copy, { - let range = rand::distributions::Uniform::new_inclusive(-range, range); + let range = rand::distr::Uniform::new_inclusive(-range, range).unwrap(); let mut rng = rand::rngs::StdRng::from_seed(seed); core::iter::from_fn(move || Some(Point2::new(rng.sample(range), rng.sample(range)))) } @@ -43,7 +43,7 @@ pub fn random_walk_distribution( where S::Sampler: Copy, { - let range = rand::distributions::Uniform::new_inclusive(-step_size, step_size); + let range = rand::distr::Uniform::new_inclusive(-step_size, step_size).unwrap(); let mut last_x = S::zero(); let mut last_y = S::one(); diff --git a/benches/locate_benchmark.rs b/benches/locate_benchmark.rs index 59b71b2..83ab9a7 100644 --- a/benches/locate_benchmark.rs +++ b/benches/locate_benchmark.rs @@ -1,7 +1,7 @@ use std::time::Duration; use criterion::{measurement::WallTime, BenchmarkGroup, Criterion}; -use rand::distributions::uniform::SampleUniform; +use rand::distr::uniform::SampleUniform; use spade::{ DelaunayTriangulation, HierarchyHintGeneratorWithBranchFactor, HintGenerator, Point2, SpadeNum, Triangulation, diff --git a/delaunay_compare/examples/real_data_benchmark.rs b/delaunay_compare/examples/real_data_benchmark.rs index 21b174f..7bcc289 100644 --- a/delaunay_compare/examples/real_data_benchmark.rs +++ b/delaunay_compare/examples/real_data_benchmark.rs @@ -37,15 +37,14 @@ fn main() -> anyhow::Result<()> { vertices.shrink_to_fit(); edges.shrink_to_fit(); - load_with_spade(&vertices, &edges)?; - println!(); - load_with_cdt_crate(&vertices, &edges)?; + println!(); + load_with_spade(vertices, edges)?; Ok(()) } -fn load_with_spade(vertices: &Vec>, edges: &Vec<[usize; 2]>) -> anyhow::Result<()> { +fn load_with_spade(vertices: Vec>, edges: Vec<[usize; 2]>) -> anyhow::Result<()> { let vertices_clone = vertices.clone(); let edges_clone = edges.clone(); @@ -172,8 +171,10 @@ fn draw_to_pixmap(cdt: ConstrainedDelaunayTriangulation>) -> anyhow: .post_scale(scale_x as f32, -scale_y as f32) .post_translate(0.0, res_y as f32); - let mut stroke = Stroke::default(); - stroke.width = 0.1 / scale_x as f32; + let mut stroke = Stroke { + width: 0.1 / scale_x as f32, + ..Default::default() + }; let mut paint = Paint::default(); diff --git a/examples/interpolation/main.rs b/examples/interpolation/main.rs index d38c67c..681d42d 100644 --- a/examples/interpolation/main.rs +++ b/examples/interpolation/main.rs @@ -44,7 +44,6 @@ const VERTICES: &[PointWithHeight] = &[ PointWithHeight::new(-10.0, 10.0, 0.0), PointWithHeight::new(-10.0, -10.0, 0.1), PointWithHeight::new(-15.0, 20.0, 0.5), - PointWithHeight::new(-20.0, 20.0, 0.0), PointWithHeight::new(-5.0, 0.0, 0.25), PointWithHeight::new(5.0, 7.0, 0.75), PointWithHeight::new(12.0, -10.0, 0.4), @@ -71,7 +70,9 @@ pub fn main() -> Result<()> { let dimensions = 512; let mut nn_c0_pixmap = Pixmap::new(dimensions, dimensions).context("Failed to allocate image")?; - let mut nn_c1_pixmap = + let mut nn_c1_pixmap_1 = + Pixmap::new(dimensions, dimensions).context("Failed to allocate image")?; + let mut nn_c1_pixmap_05 = Pixmap::new(dimensions, dimensions).context("Failed to allocate image")?; let mut barycentric_pixmap = Pixmap::new(dimensions, dimensions).context("Failed to allocate image")?; @@ -92,15 +93,19 @@ pub fn main() -> Result<()> { let nn = t.natural_neighbor(); let barycentric = t.barycentric(); + let grads = nn.estimate_gradients(|v| v.data().height); + for y in 0..dimensions { for x in 0..dimensions { let coords = px_to_coords(x, y); let value_nn_c0 = nn.interpolate(|v| v.data().height, coords); set_pixel(&mut nn_c0_pixmap, x, y, value_nn_c0); - let value_nn_c1 = - nn.interpolate_gradient(|v| v.data().height, |_| [0.0, 0.0], 1.0, coords); - set_pixel(&mut nn_c1_pixmap, x, y, value_nn_c1); + let value_nn_c1_1 = nn.interpolate_gradient(|v| v.data().height, &grads, 1.0, coords); + set_pixel(&mut nn_c1_pixmap_1, x, y, value_nn_c1_1); + + let value_nn_c1_05 = nn.interpolate_gradient(|v| v.data().height, &grads, 0.5, coords); + set_pixel(&mut nn_c1_pixmap_05, x, y, value_nn_c1_05); let value_barycentric = barycentric.interpolate(|v| v.data().height, coords); set_pixel(&mut barycentric_pixmap, x, y, value_barycentric); @@ -133,7 +138,8 @@ pub fn main() -> Result<()> { for pixmap in [ &mut nn_c0_pixmap, - &mut nn_c1_pixmap, + &mut nn_c1_pixmap_1, + &mut nn_c1_pixmap_05, &mut barycentric_pixmap, &mut nearest_neighbor_pixmap, ] { @@ -170,7 +176,8 @@ pub fn main() -> Result<()> { } save_pixmap(nn_c0_pixmap, "interpolation_nn_c0")?; - save_pixmap(nn_c1_pixmap, "interpolation_nn_c1")?; + save_pixmap(nn_c1_pixmap_1, "interpolation_nn_c1_flatness_1")?; + save_pixmap(nn_c1_pixmap_05, "interpolation_nn_c1_flatness_05")?; save_pixmap(barycentric_pixmap, "interpolation_barycentric")?; save_pixmap(nearest_neighbor_pixmap, "interpolation_nearest_neighbor")?; diff --git a/images/interpolation_barycentric.img b/images/interpolation_barycentric.img index 91a7dab..a3cd4e0 100644 --- a/images/interpolation_barycentric.img +++ b/images/interpolation_barycentric.img @@ -1 +1 @@ - \ No newline at end of file + \ No newline at end of file diff --git a/images/interpolation_barycentric.jpg b/images/interpolation_barycentric.jpg index ea36e85..baba804 100644 Binary files a/images/interpolation_barycentric.jpg and b/images/interpolation_barycentric.jpg differ diff --git a/images/interpolation_nearest_neighbor.img b/images/interpolation_nearest_neighbor.img index 621d3c7..47f13bd 100644 --- a/images/interpolation_nearest_neighbor.img +++ b/images/interpolation_nearest_neighbor.img @@ -1 +1 @@ - \ No newline at end of file + \ No newline at end of file diff --git a/images/interpolation_nearest_neighbor.jpg b/images/interpolation_nearest_neighbor.jpg index aa457bc..e1c9c50 100644 Binary files a/images/interpolation_nearest_neighbor.jpg and b/images/interpolation_nearest_neighbor.jpg differ diff --git a/images/interpolation_nn_c0.img b/images/interpolation_nn_c0.img index 9412ea4..3ed6f76 100644 --- a/images/interpolation_nn_c0.img +++ b/images/interpolation_nn_c0.img @@ -1 +1 @@ - \ No newline at end of file + \ No newline at end of file diff --git a/images/interpolation_nn_c0.jpg b/images/interpolation_nn_c0.jpg index 1541960..15ca142 100644 Binary files a/images/interpolation_nn_c0.jpg and b/images/interpolation_nn_c0.jpg differ diff --git a/images/interpolation_nn_c1.img b/images/interpolation_nn_c1.img deleted file mode 100644 index dcc8e10..0000000 --- a/images/interpolation_nn_c1.img +++ /dev/null @@ -1 +0,0 @@ - \ No newline at end of file diff --git a/images/interpolation_nn_c1.jpg b/images/interpolation_nn_c1.jpg deleted file mode 100644 index 8a4591b..0000000 Binary files a/images/interpolation_nn_c1.jpg and /dev/null differ diff --git a/images/interpolation_nn_c1_flatness_05.img b/images/interpolation_nn_c1_flatness_05.img new file mode 100644 index 0000000..6767d63 --- /dev/null +++ b/images/interpolation_nn_c1_flatness_05.img @@ -0,0 +1 @@ + \ No newline at end of file diff --git a/images/interpolation_nn_c1_flatness_05.jpg b/images/interpolation_nn_c1_flatness_05.jpg new file mode 100644 index 0000000..fe4f2e2 Binary files /dev/null and b/images/interpolation_nn_c1_flatness_05.jpg differ diff --git a/images/interpolation_nn_c1_flatness_1.img b/images/interpolation_nn_c1_flatness_1.img new file mode 100644 index 0000000..cc44a7d --- /dev/null +++ b/images/interpolation_nn_c1_flatness_1.img @@ -0,0 +1 @@ + \ No newline at end of file diff --git a/images/interpolation_nn_c1_flatness_1.jpg b/images/interpolation_nn_c1_flatness_1.jpg new file mode 100644 index 0000000..f32750e Binary files /dev/null and b/images/interpolation_nn_c1_flatness_1.jpg differ diff --git a/src/cdt.rs b/src/cdt.rs index 3eda667..b49aa60 100644 --- a/src/cdt.rs +++ b/src/cdt.rs @@ -2003,7 +2003,7 @@ mod test { assert_eq!(cdt.num_vertices(), initial_num_vertices + 1); assert_eq!(edges.len(), 2); - check_returned_edges(&mut cdt, &edges, from, to); + check_returned_edges(&cdt, &edges, from, to); Ok(()) } @@ -2021,7 +2021,7 @@ mod test { // 3 new points should be added as the constraint intersects all 3 existing edges assert_eq!(cdt.num_vertices(), initial_num_vertices + 3); assert_eq!(edges.len(), 4); - check_returned_edges(&mut cdt, &edges, from, to); + check_returned_edges(&cdt, &edges, from, to); Ok(()) } diff --git a/src/delaunay_core/bulk_load.rs b/src/delaunay_core/bulk_load.rs index 73bc145..9d51f0e 100644 --- a/src/delaunay_core/bulk_load.rs +++ b/src/delaunay_core/bulk_load.rs @@ -788,7 +788,7 @@ impl Hull { const INVALID: usize = usize::MAX; self.buckets - .extend(core::iter::repeat(INVALID).take(target_size)); + .extend(core::iter::repeat_n(INVALID, target_size)); let (first_index, current_node) = self .data diff --git a/src/delaunay_core/hint_generator.rs b/src/delaunay_core/hint_generator.rs index b0340cd..e710b02 100644 --- a/src/delaunay_core/hint_generator.rs +++ b/src/delaunay_core/hint_generator.rs @@ -26,7 +26,7 @@ use alloc::vec::Vec; /// fulfill most needs: /// - A heuristic that uses the last inserted vertex as hint ([LastUsedVertexHintGenerator]) /// - A hint generator based on a hierarchy of triangulations that improves walk time to `O(log(n))` -/// ([HierarchyHintGenerator]) +/// ([HierarchyHintGenerator]) pub trait HintGenerator: Default { /// Returns a vertex handle that should be close to a given position. /// diff --git a/src/delaunay_core/interpolation.rs b/src/delaunay_core/interpolation.rs index 4c699c9..01e2017 100644 --- a/src/delaunay_core/interpolation.rs +++ b/src/delaunay_core/interpolation.rs @@ -109,10 +109,11 @@ use super::VertexHandle; #[doc = include_str!("../../images/interpolation_nn_c0.img")] /// /// With a gradient, the sharp peaks are gone - the surface will smoothly approximate a linear function as defined -/// by the gradient in the vicinity of each vertex. In the image below, a gradient of `(0.0, 0.0)` is used which -/// leads to a small "disc" around each vertex with values close to the vertex value. -/// -#[doc = include_str!("../../images/interpolation_nn_c1.img")] +/// by the gradient in the vicinity of each vertex. The gradients are estimated with the estimate_gradients function. +/// With flatness = 0.5 +#[doc = include_str!("../../images/interpolation_nn_c1_flatness_05.img")] +/// With flatness = 1.0 +#[doc = include_str!("../../images/interpolation_nn_c1_flatness_1.img")] /// #[doc(alias = "Interpolation")] pub struct NaturalNeighbor<'a, T> @@ -324,7 +325,39 @@ where /// /// Refer to [NaturalNeighbor] for more information and a visual example. /// - /// # Example + /// # Example 1 + /// + /// ``` + /// use spade::{DelaunayTriangulation, HasPosition, Point2}; + /// # use spade::Triangulation; + /// + /// struct PointWithData { + /// position: Point2, + /// height: f64, + /// grad: [f64; 2], + /// } + /// + /// impl HasPosition for PointWithData { + /// type Scalar = f64; + /// fn position(&self) -> Point2 { self.position } + /// } + /// + /// let mut triangulation: DelaunayTriangulation = Default::default(); + /// // Insert some points into the triangulation + /// triangulation.insert(PointWithData { position: Point2::new(10.0, 10.0), height: 0.0, grad: [-1., -1.] }); + /// triangulation.insert(PointWithData { position: Point2::new(10.0, -10.0), height: 0.0, grad: [-1., 1.] }); + /// triangulation.insert(PointWithData { position: Point2::new(-10.0, 10.0), height: 0.0, grad: [1., -1.] }); + /// triangulation.insert(PointWithData { position: Point2::new(-10.0, -10.0), height: 0.0, grad: [1., 1.] }); + /// + /// let nn = triangulation.natural_neighbor(); + /// + /// // Interpolate point at coordinates (1.0, 2.0). This example uses a known gradients at each vertex. + /// // The gradient is the direction of steepest ascent for the function at that point + /// let query_point = Point2::new(1.0, 2.0); + /// let value: f64 = nn.interpolate_gradient(|v| v.data().height, |v| v.data().grad, 1.0, query_point).unwrap(); + /// ``` + /// + /// /// # Example 2 /// /// ``` /// use spade::{DelaunayTriangulation, HasPosition, Point2}; @@ -349,11 +382,14 @@ where /// /// let nn = triangulation.natural_neighbor(); /// - /// // Interpolate point at coordinates (1.0, 2.0). This example uses a fixed gradient of (0.0, 0.0) which - /// // means that the interpolation will have normal vector parallel to the z-axis at each input point. - /// // Realistically, the gradient might be stored as an additional property of `PointWithHeight`. + /// // Interpolate point at coordinates (1.0, 2.0). /// let query_point = Point2::new(1.0, 2.0); - /// let value: f64 = nn.interpolate_gradient(|v| v.data().height, |_| [0.0, 0.0], 1.0, query_point).unwrap(); + /// // Let's interpolate the value with estimated gradients. The gradient is the direction of steepest ascent for the function at that point + /// let value: f64 = nn.interpolate_gradient(|v| v.data().height, |v| nn.estimate_gradient(v, |h| h.data().height), 1.0, query_point).unwrap(); + /// + /// // could also have gone and calculated all grads upfront like this, notice the s in gradients + /// let grads = nn.estimate_gradients(|v| v.data().height); + /// let value: f64 = nn.interpolate_gradient(|v| v.data().height, &grads, 1.0, query_point).unwrap(); /// ``` /// /// # References @@ -398,6 +434,11 @@ where let h_i = i(handle); let diff = position.sub(pos_i); let r_i2 = diff.length2(); + + if r_i2 == zero() { + return Some(h_i); + } + let r_i = r_i2.powf(flatness); let c1_weight_i = *weight / r_i; let grad_i = g(handle); @@ -573,6 +614,91 @@ where tuple.1 = tuple.1 / total_area; } } + + /// Estimates and returns the gradient for a every vertex in this triangulation. + /// + /// This assumes that the triangulation models some kind of height field. + /// The gradient is calculated from the weighted and normalized average of the normals of all triangles + /// adjacent to the given vertex (weighted by triangle size). + /// + /// NOTE: This is not Sibson's gradient estimator + pub fn estimate_gradients( + &self, + i: I, + ) -> impl Fn(VertexHandle) -> [::Scalar; 2] + where + I: Fn(VertexHandle) -> ::Scalar, + { + let grads = self + .triangulation + .vertices() + .map(|v| self.estimate_gradient(v, &i)) + .collect::>(); + + move |v: VertexHandle| grads[v.index()] + } + + /// Estimates and returns the gradient for a single vertex in this triangulation. + /// + /// This assumes that the triangulation models some kind of height field. + /// The gradient is calculated from the weighted and normalized average of the normals of all triangles + /// adjacent to the given vertex (weighted by triangle size). + /// + /// NOTE: This is not Sibson's gradient estimator + pub fn estimate_gradient( + &self, + v: VertexHandle, + i: I, + ) -> [::Scalar; 2] + where + I: Fn(VertexHandle) -> ::Scalar, + { + let v_2d = v.position(); + let v_pos = [v_2d.x, v_2d.y, i(v)]; + + let neighbor_positions = { + v.out_edges() + .map(|e| { + let pos = e.to().position(); + [pos.x, pos.y, i(e.to())] + }) + .collect::>() + }; + let mut final_normal: [::Scalar; 3] = [zero(); 3]; + for index in 0..neighbor_positions.len() { + let p0 = neighbor_positions[index]; + let p1 = neighbor_positions[(index + 1) % neighbor_positions.len()]; + + let d0 = [p0[0] - v_pos[0], p0[1] - v_pos[1], p0[2] - v_pos[2]]; + let d1 = [p1[0] - v_pos[0], p1[1] - v_pos[1], p1[2] - v_pos[2]]; + + let normal = [ + d0[1] * d1[2] - d0[2] * d1[1], + d0[2] * d1[0] - d0[0] * d1[2], + d0[0] * d1[1] - d0[1] * d1[0], + ]; + + // this should be true for every normal + // given that the height field assumption holds + if normal[2] > zero() { + final_normal = [ + final_normal[0] + normal[0], + final_normal[1] + normal[1], + final_normal[2] + normal[2], + ]; + } + } + + // Calculate gradient from normal + if final_normal[2] != zero() { + [ + -final_normal[0] / final_normal[2], + -final_normal[1] / final_normal[2], + ] + } else { + [zero(), zero()] + } + } } fn get_natural_neighbor_edges( @@ -970,4 +1096,54 @@ mod test { assert!(!result.is_nan()); Ok(()) } + + #[test] + fn test_gradient_estimation_planar() -> Result<(), InsertionError> { + let points = vec![ + PointWithHeight::new(Point2::new(0., 0.), 0.), + PointWithHeight::new(Point2::new(1., 0.), 0.), + PointWithHeight::new(Point2::new(1., 1.), 1.), + PointWithHeight::new(Point2::new(0., 1.), 1.), + PointWithHeight::new(Point2::new(0.5, 0.5), 0.5), + ]; + + let t = DelaunayTriangulation::<_>::bulk_load_stable(points).unwrap(); + let nn = t.natural_neighbor(); + + let v = nn + .triangulation + .locate_vertex(Point2::new(0.5, 0.5)) + .unwrap(); + + let grad = nn.estimate_gradient(v, |v| v.data().height); + + assert!(grad == [0., 1.0]); + + Ok(()) + } + + #[test] + fn test_gradient_estimation_flat() -> Result<(), InsertionError> { + let points = vec![ + PointWithHeight::new(Point2::new(0., 0.), 1.), + PointWithHeight::new(Point2::new(1., 0.), 1.), + PointWithHeight::new(Point2::new(1., 1.), 1.), + PointWithHeight::new(Point2::new(0., 1.), 1.), + PointWithHeight::new(Point2::new(0.5, 0.5), 0.), + ]; + + let t = DelaunayTriangulation::<_>::bulk_load_stable(points).unwrap(); + let nn = t.natural_neighbor(); + + let v = nn + .triangulation + .locate_vertex(Point2::new(0.5, 0.5)) + .unwrap(); + + let grad = nn.estimate_gradient(v, |v| v.data().height); + + assert!(grad == [0., 0.]); + + Ok(()) + } } diff --git a/src/lib.rs b/src/lib.rs index 472a081..ff0fc65 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -84,7 +84,7 @@ mod test_utilities; /// Reference handles come in one of four variants: /// * [FaceHandle](handles::FaceHandle)s refer to a single face (triangle) of the triangulation. /// They are used get the triangle's adjacent vertices and edges. They also may refer to -/// the single outer face. +/// the single outer face. /// * [VertexHandle](handles::VertexHandle)s refer to a single vertex of the triangulation. They /// allow to retrieve the vertex position and its outgoing edges. /// * [DirectedEdgeHandle](handles::DirectedEdgeHandle)s refer to a single directed edge. They