From 913386dafb28b00b6ebc7e3a90235f0dc8ee12a1 Mon Sep 17 00:00:00 2001 From: Daniel Oom Date: Fri, 9 Aug 2024 11:22:23 +0200 Subject: [PATCH] Split flat kdtree cells --- kdtree/src/build.rs | 147 +++++++++++++++++++++++++++++++------ kdtree/src/partitioning.rs | 3 +- kdtree/src/sah.rs | 120 +++++++++++++++++++++--------- 3 files changed, 213 insertions(+), 57 deletions(-) diff --git a/kdtree/src/build.rs b/kdtree/src/build.rs index 7f8931e5..d71f378f 100644 --- a/kdtree/src/build.rs +++ b/kdtree/src/build.rs @@ -17,7 +17,7 @@ fn starting_box(geometries: &[Geometry]) -> KdCell { fn build_helper( geometries: &[Geometry], - cost: &SahCost, + sah: &SahCost, max_depth: u32, depth: u32, cell: KdCell, @@ -26,28 +26,28 @@ fn build_helper( return KdNode::new_leaf(cell.indices); } - match find_best_split(geometries, cost, &cell) { + match find_best_split(geometries, sah, &cell) { None => KdNode::new_leaf(cell.indices), Some(split) => { - if should_terminate(cost, &cell, &split) { + if should_terminate(sah, &cell, &split) { KdNode::new_leaf(cell.indices) } else { - let left = build_helper(geometries, cost, max_depth, depth + 1, split.left); - let right = build_helper(geometries, cost, max_depth, depth + 1, split.right); + let left = build_helper(geometries, sah, max_depth, depth + 1, split.left); + let right = build_helper(geometries, sah, max_depth, depth + 1, split.right); KdNode::new_node(split.plane, left, right) } } } } -pub fn build_kdtree(geometries: &[Geometry], max_depth: u32, cost: &SahCost) -> KdNode { +pub fn build_kdtree(geometries: &[Geometry], max_depth: u32, sah: &SahCost) -> KdNode { if max_depth as usize > MAX_DEPTH { panic!( "Max depth ({}) must be smaller than hard coded value ({}).", max_depth, MAX_DEPTH ); } - *build_helper(geometries, cost, max_depth, 1, starting_box(geometries)) + *build_helper(geometries, sah, max_depth, 1, starting_box(geometries)) } #[cfg(test)] @@ -60,7 +60,7 @@ mod tests { use super::*; #[test] - fn non_axially_aligned_triangle() { + fn two_oriented_triangles() { let triangle1 = Triangle { v0: Vec3::new(0.0, 0.0, 0.0), v1: Vec3::new(1.0, 0.0, 0.0), @@ -72,12 +72,12 @@ mod tests { v2: Vec3::new(2.0, 1.0, 1.0), }; let geometries = [triangle1, triangle2].map(Geometry::from); - let cost = SahCost { + let sah = SahCost { traverse_cost: 0.1, intersect_cost: 1.0, empty_factor: 0.8, }; - let actual = build_kdtree(&geometries, 7, &cost); + let actual = build_kdtree(&geometries, 7, &sah); let expected = KdNode::new_node( Aap::new_x(1.0), @@ -92,7 +92,7 @@ mod tests { } #[test] - fn axially_aligned_triangle() { + fn two_axially_aligned_triangles() { let triangle1 = Triangle { v0: Vec3::new(0.0, 0.0, 0.0), v1: Vec3::new(1.0, 0.0, 0.0), @@ -103,23 +103,128 @@ mod tests { v1: Vec3::new(1.0, 0.0, 1.0), v2: Vec3::new(1.0, 1.0, 1.0), }; - let triangle3 = Triangle { - v0: Vec3::new(0.0, 0.0, 2.0), - v1: Vec3::new(1.0, 0.0, 2.0), - v2: Vec3::new(1.0, 1.0, 2.0), + let geometries = [triangle1, triangle2].map(Geometry::from); + let sah = SahCost { + traverse_cost: 0.0, + intersect_cost: 1.0, + empty_factor: 1.0, }; - let geometries = [triangle1, triangle2, triangle3].map(Geometry::from); - let cost = SahCost { + let actual = build_kdtree(&geometries, 6, &sah); + + let expected = KdNode::new_node( + Aap::new_z(0.0), + KdNode::new_leaf(vec![0]), + KdNode::new_node(Aap::new_z(1.0), KdNode::empty(), KdNode::new_leaf(vec![1])), + ); + assert_eq!( + actual, *expected, + "\n actual: {}\n expected: {}", + actual, *expected + ); + } + + #[test] + fn one_cube() { + let triangles = [ + // Front + Triangle { + v0: Vec3::new(0.0, 0.0, 0.0), + v1: Vec3::new(1.0, 0.0, 0.0), + v2: Vec3::new(1.0, 1.0, 0.0), + }, + Triangle { + v0: Vec3::new(0.0, 0.0, 0.0), + v1: Vec3::new(0.0, 1.0, 0.0), + v2: Vec3::new(1.0, 1.0, 0.0), + }, + // Back + Triangle { + v0: Vec3::new(0.0, 0.0, 1.0), + v1: Vec3::new(1.0, 0.0, 1.0), + v2: Vec3::new(1.0, 1.0, 1.0), + }, + Triangle { + v0: Vec3::new(0.0, 0.0, 1.0), + v1: Vec3::new(0.0, 1.0, 1.0), + v2: Vec3::new(1.0, 1.0, 1.0), + }, + // Bottom + Triangle { + v0: Vec3::new(0.0, 0.0, 0.0), + v1: Vec3::new(1.0, 0.0, 0.0), + v2: Vec3::new(1.0, 0.0, 1.0), + }, + Triangle { + v0: Vec3::new(0.0, 0.0, 0.0), + v1: Vec3::new(0.0, 0.0, 1.0), + v2: Vec3::new(1.0, 0.0, 1.0), + }, + // Top + Triangle { + v0: Vec3::new(0.0, 1.0, 0.0), + v1: Vec3::new(1.0, 1.0, 0.0), + v2: Vec3::new(1.0, 1.0, 1.0), + }, + Triangle { + v0: Vec3::new(0.0, 1.0, 0.0), + v1: Vec3::new(0.0, 1.0, 1.0), + v2: Vec3::new(1.0, 1.0, 1.0), + }, + // Right + Triangle { + v0: Vec3::new(0.0, 0.0, 0.0), + v1: Vec3::new(0.0, 0.0, 1.0), + v2: Vec3::new(0.0, 1.0, 1.0), + }, + Triangle { + v0: Vec3::new(0.0, 0.0, 0.0), + v1: Vec3::new(0.0, 1.0, 0.0), + v2: Vec3::new(0.0, 1.0, 1.0), + }, + // Left + Triangle { + v0: Vec3::new(1.0, 0.0, 0.0), + v1: Vec3::new(1.0, 1.0, 0.0), + v2: Vec3::new(1.0, 1.0, 1.0), + }, + Triangle { + v0: Vec3::new(1.0, 0.0, 0.0), + v1: Vec3::new(1.0, 0.0, 1.0), + v2: Vec3::new(1.0, 1.0, 1.0), + }, + ]; + let geometries = triangles.map(Geometry::from); + let sah = SahCost { traverse_cost: 0.0, intersect_cost: 1.0, empty_factor: 1.0, }; - let actual = build_kdtree(&geometries, 6, &cost); + let actual = build_kdtree(&geometries, 10, &sah); let expected = KdNode::new_node( - Aap::new_z(1.0), - KdNode::new_leaf(vec![0, 1]), - KdNode::new_leaf(vec![2, 1]), + Aap::new_x(0.0), + KdNode::new_leaf(vec![8, 9]), + KdNode::new_node( + Aap::new_x(1.0), + KdNode::new_node( + Aap::new_y(0.0), + KdNode::new_leaf(vec![4, 5]), + KdNode::new_node( + Aap::new_y(1.0), + KdNode::new_node( + Aap::new_z(0.0), + KdNode::new_leaf(vec![0, 1]), + KdNode::new_node( + Aap::new_z(1.0), + KdNode::empty(), + KdNode::new_leaf(vec![2, 3]), + ), + ), + KdNode::new_leaf(vec![6, 7]), + ), + ), + KdNode::new_leaf(vec![10, 11]), + ), ); assert_eq!( actual, *expected, diff --git a/kdtree/src/partitioning.rs b/kdtree/src/partitioning.rs index fcda1841..2577dd9d 100644 --- a/kdtree/src/partitioning.rs +++ b/kdtree/src/partitioning.rs @@ -44,7 +44,8 @@ pub(crate) fn split_and_partition_clipped_geometries( plane: Aap, ) -> KdPartitioning { let (left_aabb, right_aabb) = parent_aabb.split(&plane); - let (left_indices, middle_indices, right_indices) = partition_clipped_geometries(clipped, &plane); + let (left_indices, middle_indices, right_indices) = + partition_clipped_geometries(clipped, &plane); KdPartitioning { plane, parent_aabb, diff --git a/kdtree/src/sah.rs b/kdtree/src/sah.rs index ebef180d..fad9349f 100644 --- a/kdtree/src/sah.rs +++ b/kdtree/src/sah.rs @@ -12,10 +12,16 @@ pub struct SahCost { } impl SahCost { - fn calculate(&self, probability: (f32, f32), counts: (usize, usize)) -> f32 { + fn calculate( + &self, + volume: (f32, f32), + probability: (f32, f32), + counts: (usize, usize), + ) -> f32 { debug_assert!((0.0..=1.0).contains(&probability.0) && (0.0..=1.0).contains(&probability.1)); debug_assert!(probability.0 > 0.0 || probability.1 > 0.0); - let empty_factor = if counts.0 == 0 || counts.1 == 0 { + // TODO: if empty side is flat, apply no empty factor + let empty_factor = if counts.0 == 0 && volume.0 > 0.01 || counts.1 == 0 && volume.1 > 0.01 { self.empty_factor } else { 1.0 @@ -24,25 +30,13 @@ impl SahCost { * (probability.0 * counts.0 as f32 + probability.1 * counts.1 as f32); empty_factor * (self.traverse_cost + intersect_cost) } - - fn calculate_for_split(&self, partitioning: &KdPartitioning) -> f32 { - let surface_area = partitioning.parent_aabb.surface_area(); - let probability_left = partitioning.left_aabb.surface_area() / surface_area; - let probability_right = partitioning.right_aabb.surface_area() / surface_area; - let probability = (probability_left, probability_right); - let counts = ( - partitioning.left_indices.len() + partitioning.middle_indices.len(), - partitioning.right_indices.len() + partitioning.middle_indices.len(), - ); - self.calculate(probability, counts) - } } impl Default for SahCost { fn default() -> Self { SahCost { - traverse_cost: 2.0, - intersect_cost: 1.0, + traverse_cost: 1.0, + intersect_cost: 2.0, empty_factor: 0.8, } } @@ -55,27 +49,82 @@ pub(crate) struct KdSplit { pub(crate) right: KdCell, } -fn calculate_cost(cost: &SahCost, partitioning: KdPartitioning) -> Option<(KdSplit, f32)> { - // TODO: Place planes to the left or to the right depending on what gives best cost. - if (partitioning.left_aabb.volume() == 0.0 || partitioning.right_aabb.volume() == 0.0) - && partitioning.middle_indices.is_empty() - { - return None; - } - let cost = cost.calculate_for_split(&partitioning); - let mut left_indices = partitioning.left_indices; - let mut right_indices = partitioning.right_indices; - left_indices.extend(&partitioning.middle_indices); - right_indices.extend(partitioning.middle_indices); +fn select_best_partitioning(sah: &SahCost, partitioning: KdPartitioning) -> Option<(KdSplit, f32)> { + let volume = ( + partitioning.left_aabb.volume(), + partitioning.right_aabb.volume(), + ); + let parent_surface_area = partitioning.parent_aabb.surface_area(); + let probability = ( + partitioning.left_aabb.surface_area() / parent_surface_area, + partitioning.right_aabb.surface_area() / parent_surface_area, + ); + let (cost, left_indices, right_indices) = if partitioning.middle_indices.is_empty() { + let counts = ( + partitioning.left_indices.len(), + partitioning.right_indices.len(), + ); + if counts.0 == 0 && volume.0 <= 0.01 || counts.1 == 0 && volume.1 <= 0.01 { + return None; + } + let cost = sah.calculate(volume, probability, counts); + (cost, partitioning.left_indices, partitioning.right_indices) + } else if volume.0 <= 0.01 { + let counts = ( + partitioning.left_indices.len() + partitioning.middle_indices.len(), + partitioning.right_indices.len(), + ); + if counts.0 == 0 && volume.0 <= 0.01 || counts.1 == 0 && volume.1 <= 0.01 { + return None; + } + let cost = sah.calculate(volume, probability, counts); + let mut left_indices = partitioning.left_indices; + left_indices.extend(partitioning.middle_indices); + (cost, left_indices, partitioning.right_indices) + } else if volume.1 <= 0.01 { + let counts = ( + partitioning.left_indices.len(), + partitioning.right_indices.len() + partitioning.middle_indices.len(), + ); + if counts.0 == 0 && volume.0 <= 0.01 || counts.1 == 0 && volume.1 <= 0.01 { + return None; + } + let cost = sah.calculate(volume, probability, counts); + let mut right_indices = partitioning.right_indices; + right_indices.extend(partitioning.middle_indices); + (cost, partitioning.left_indices, right_indices) + } else { + let counts_middle_left = ( + partitioning.left_indices.len() + partitioning.middle_indices.len(), + partitioning.right_indices.len(), + ); + let counts_middle_right = ( + partitioning.left_indices.len(), + partitioning.right_indices.len() + partitioning.middle_indices.len(), + ); + let cost = ( + sah.calculate(volume, probability, counts_middle_left), + sah.calculate(volume, probability, counts_middle_right), + ); + if cost.0 <= cost.1 { + let mut left_indices = partitioning.left_indices; + left_indices.extend(partitioning.middle_indices); + (cost.0, left_indices, partitioning.right_indices) + } else { + let mut right_indices = partitioning.right_indices; + right_indices.extend(partitioning.middle_indices); + (cost.1, partitioning.left_indices, right_indices) + } + }; + let plane = partitioning.plane; let left = KdCell::new(partitioning.left_aabb, left_indices); let right = KdCell::new(partitioning.right_aabb, right_indices); - let plane = partitioning.plane; Some((KdSplit { plane, left, right }, cost)) } pub(crate) fn find_best_split( geometries: &[Geometry], - cost: &SahCost, + sah: &SahCost, cell: &KdCell, ) -> Option { debug_assert!( @@ -96,8 +145,8 @@ pub(crate) fn find_best_split( splits .into_iter() .filter_map(|plane| { - calculate_cost( - cost, + select_best_partitioning( + sah, split_and_partition_clipped_geometries(&clipped, cell.boundary.clone(), plane), ) }) @@ -107,8 +156,8 @@ pub(crate) fn find_best_split( splits .into_par_iter() .filter_map(|plane| { - calculate_cost( - cost, + select_best_partitioning( + sah, split_and_partition_clipped_geometries(&clipped, cell.boundary.clone(), plane), ) }) @@ -118,12 +167,13 @@ pub(crate) fn find_best_split( } pub(crate) fn should_terminate(cost: &SahCost, cell: &KdCell, split: &KdSplit) -> bool { + let volume = (split.left.boundary.volume(), split.right.boundary.volume()); let surface_area = cell.boundary.surface_area(); let probability_left = split.left.boundary.surface_area() / surface_area; let probability_right = split.right.boundary.surface_area() / surface_area; let probability = (probability_left, probability_right); let counts = (split.left.indices.len(), split.right.indices.len()); - let split_cost = cost.calculate(probability, counts); + let split_cost = cost.calculate(volume, probability, counts); let intersect_cost = cost.intersect_cost * cell.indices.len() as f32; split_cost >= intersect_cost }