From 2a42402bfa9c2a8b109e6f35908a09f464fd2f27 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Crozet?= Date: Fri, 3 Oct 2025 17:43:28 +0200 Subject: [PATCH 1/3] =?UTF-8?q?feat:=E2=80=AFmake=20the=20voxels=20shape?= =?UTF-8?q?=20use=20a=20sparse=20internal=20storage?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- crates/parry2d-f64/Cargo.toml | 4 +- crates/parry2d/Cargo.toml | 4 +- crates/parry3d-f64/Cargo.toml | 4 +- crates/parry3d/Cargo.toml | 4 +- src/bounding_volume/aabb_voxels.rs | 11 +- src/bounding_volume/bounding_sphere_voxels.rs | 11 +- src/partitioning/bvh/bvh_queries.rs | 25 + ...ontact_manifolds_voxels_composite_shape.rs | 2 +- .../contact_manifolds_voxels_shape.rs | 20 +- .../contact_manifolds_voxels_voxels.rs | 4 +- src/query/point/point_voxels.rs | 70 +- src/query/ray/ray_voxels.rs | 49 +- src/shape/mod.rs | 2 +- src/shape/voxels.rs | 1555 ----------------- src/shape/voxels/mod.rs | 11 + src/shape/voxels/voxels.rs | 535 ++++++ src/shape/voxels/voxels_chunk.rs | 274 +++ src/shape/voxels/voxels_consts.rs | 639 +++++++ src/shape/voxels/voxels_edition.rs | 145 ++ src/shape/voxels/voxels_neighborhood.rs | 197 +++ src/shape/voxels/voxels_sizeof.rs | 29 + .../voxelization/voxelized_volume.rs | 6 +- 22 files changed, 1982 insertions(+), 1619 deletions(-) delete mode 100644 src/shape/voxels.rs create mode 100644 src/shape/voxels/mod.rs create mode 100644 src/shape/voxels/voxels.rs create mode 100644 src/shape/voxels/voxels_chunk.rs create mode 100644 src/shape/voxels/voxels_consts.rs create mode 100644 src/shape/voxels/voxels_edition.rs create mode 100644 src/shape/voxels/voxels_neighborhood.rs create mode 100644 src/shape/voxels/voxels_sizeof.rs diff --git a/crates/parry2d-f64/Cargo.toml b/crates/parry2d-f64/Cargo.toml index bfa30daa..c26b27fa 100644 --- a/crates/parry2d-f64/Cargo.toml +++ b/crates/parry2d-f64/Cargo.toml @@ -36,6 +36,7 @@ dim2 = [] f64 = [] serde-serialize = [ "serde", + "serde_arrays", "nalgebra/serde-serialize", "arrayvec/serde", "bitflags/serde", @@ -73,9 +74,10 @@ num-traits = { version = "0.2", default-features = false } slab = { version = "0.4", optional = true } arrayvec = { version = "0.7", default-features = false } simba = { version = "0.9", default-features = false } -nalgebra = { version = "0.34", default-features = false, features = ["libm"] } +nalgebra = { version = "0.34", default-features = false, features = ["libm", "macros"] } approx = { version = "0.5", default-features = false } serde = { version = "1.0", optional = true, features = ["derive"] } +serde_arrays = { version = "0.2", optional = true } rkyv = { version = "0.7.41", optional = true } num-derive = "0.4" indexmap = { version = "2", features = ["serde"], optional = true } diff --git a/crates/parry2d/Cargo.toml b/crates/parry2d/Cargo.toml index 56d5cef0..e489bff7 100644 --- a/crates/parry2d/Cargo.toml +++ b/crates/parry2d/Cargo.toml @@ -36,6 +36,7 @@ dim2 = [] f32 = [] serde-serialize = [ "serde", + "serde_arrays", "nalgebra/serde-serialize", "arrayvec/serde", "bitflags/serde", @@ -73,9 +74,10 @@ num-traits = { version = "0.2", default-features = false } slab = { version = "0.4", optional = true } arrayvec = { version = "0.7", default-features = false } simba = { version = "0.9", default-features = false } -nalgebra = { version = "0.34", default-features = false, features = ["libm"] } +nalgebra = { version = "0.34", default-features = false, features = ["libm", "macros"] } approx = { version = "0.5", default-features = false } serde = { version = "1.0", optional = true, features = ["derive"] } +serde_arrays = { version = "0.2", optional = true } rkyv = { version = "0.7.41", optional = true } num-derive = "0.4" indexmap = { version = "2", features = ["serde"], optional = true } diff --git a/crates/parry3d-f64/Cargo.toml b/crates/parry3d-f64/Cargo.toml index ad33fa85..036d328f 100644 --- a/crates/parry3d-f64/Cargo.toml +++ b/crates/parry3d-f64/Cargo.toml @@ -35,6 +35,7 @@ dim3 = [] f64 = [] serde-serialize = [ "serde", + "serde_arrays", "nalgebra/serde-serialize", "bitflags/serde", "hashbrown?/serde", @@ -73,9 +74,10 @@ num-traits = { version = "0.2", default-features = false } slab = { version = "0.4", optional = true } arrayvec = { version = "0.7", default-features = false } simba = { version = "0.9", default-features = false } -nalgebra = { version = "0.34", default-features = false, features = ["libm"] } +nalgebra = { version = "0.34", default-features = false, features = ["libm", "macros"] } approx = { version = "0.5", default-features = false } serde = { version = "1.0", optional = true, features = ["derive", "rc"] } +serde_arrays = { version = "0.2", optional = true } rkyv = { version = "0.7.41", optional = true } num-derive = "0.4" indexmap = { version = "2", features = ["serde"], optional = true } diff --git a/crates/parry3d/Cargo.toml b/crates/parry3d/Cargo.toml index fa787efa..c31b231d 100644 --- a/crates/parry3d/Cargo.toml +++ b/crates/parry3d/Cargo.toml @@ -35,6 +35,7 @@ dim3 = [] f32 = [] serde-serialize = [ "serde", + "serde_arrays", "nalgebra/serde-serialize", "bitflags/serde", "hashbrown?/serde", @@ -74,9 +75,10 @@ num-traits = { version = "0.2", default-features = false } slab = { version = "0.4", optional = true } arrayvec = { version = "0.7", default-features = false } simba = { version = "0.9", default-features = false } -nalgebra = { version = "0.34", default-features = false, features = ["libm"] } +nalgebra = { version = "0.34", default-features = false, features = ["libm", "macros"] } approx = { version = "0.5", default-features = false } serde = { version = "1.0", optional = true, features = ["derive", "rc"] } +serde_arrays = { version = "0.2", optional = true } rkyv = { version = "0.7.41", optional = true } num-derive = "0.4" indexmap = { version = "2", features = ["serde"], optional = true } diff --git a/src/bounding_volume/aabb_voxels.rs b/src/bounding_volume/aabb_voxels.rs index e28efc7d..f59e26fa 100644 --- a/src/bounding_volume/aabb_voxels.rs +++ b/src/bounding_volume/aabb_voxels.rs @@ -1,20 +1,17 @@ use crate::bounding_volume::Aabb; -use crate::math::{Isometry, Real, Translation}; -use crate::shape::{Cuboid, Voxels}; +use crate::math::{Isometry, Real}; +use crate::shape::Voxels; impl Voxels { /// Computes the world-space Aabb of this set of voxels, transformed by `pos`. #[inline] pub fn aabb(&self, pos: &Isometry) -> Aabb { - let shift = Translation::from(self.domain_center()); - Cuboid::new(self.extents() / 2.0).aabb(&(pos * shift)) + self.chunk_bvh().root_aabb().transform_by(pos) } /// Computes the local-space Aabb of this set of voxels. #[inline] pub fn local_aabb(&self) -> Aabb { - Cuboid::new(self.extents() / 2.0) - .local_aabb() - .translated(&self.domain_center().coords) + self.chunk_bvh().root_aabb() } } diff --git a/src/bounding_volume/bounding_sphere_voxels.rs b/src/bounding_volume/bounding_sphere_voxels.rs index c4d5bf77..54e54b9a 100644 --- a/src/bounding_volume/bounding_sphere_voxels.rs +++ b/src/bounding_volume/bounding_sphere_voxels.rs @@ -1,20 +1,17 @@ use crate::bounding_volume::BoundingSphere; -use crate::math::{Isometry, Real, Translation}; -use crate::shape::{Cuboid, Voxels}; +use crate::math::{Isometry, Real}; +use crate::shape::Voxels; impl Voxels { /// Computes the world-space bounding sphere of this set of voxels, transformed by `pos`. #[inline] pub fn bounding_sphere(&self, pos: &Isometry) -> BoundingSphere { - let shift = Translation::from(self.domain_center().coords); - Cuboid::new(self.extents() / 2.0).bounding_sphere(&(pos * shift)) + self.local_aabb().bounding_sphere().transform_by(pos) } /// Computes the local-space bounding sphere of this set of voxels. #[inline] pub fn local_bounding_sphere(&self) -> BoundingSphere { - Cuboid::new(self.extents() / 2.0) - .local_bounding_sphere() - .translated(&(self.domain_center().coords)) + self.local_aabb().bounding_sphere() } } diff --git a/src/partitioning/bvh/bvh_queries.rs b/src/partitioning/bvh/bvh_queries.rs index 754cf0e3..f97883cc 100644 --- a/src/partitioning/bvh/bvh_queries.rs +++ b/src/partitioning/bvh/bvh_queries.rs @@ -4,6 +4,7 @@ use crate::math::Point; use crate::math::Real; use crate::query::PointProjection; use crate::query::{PointQuery, Ray}; +use crate::shape::FeatureId; #[cfg(all(feature = "simd-is-enabled", feature = "dim3", feature = "f32"))] pub(super) struct SimdInvRay { @@ -59,6 +60,30 @@ impl Bvh { ) } + /// Projects a point on this BVH using the provided leaf projection function. + /// + /// Also returns the feature the point was projected on. + /// + /// The `primitive_check` delegates the point-projection task to an external function that + /// is assumed to map a leaf index to an actual geometry to project on. The `Real` argument + /// given to that closure is the distance to the closest point found so far (or is equal to + /// `max_distance` if no projection was found so far). + pub fn project_point_and_get_feature( + &self, + point: &Point, + max_distance: Real, + primitive_check: impl Fn(u32, Real) -> Option<(PointProjection, FeatureId)>, + ) -> Option<(u32, (Real, (PointProjection, FeatureId)))> { + self.find_best( + max_distance, + |node: &BvhNode, _| node.aabb().distance_to_local_point(point, true), + |primitive, _| { + let proj = primitive_check(primitive, max_distance)?; + Some((na::distance(&proj.0.point, point), proj)) + }, + ) + } + /// Casts a ray on this BVH using the provided leaf ray-cast function. /// /// The `primitive_check` delegates the ray-casting task to an external function that diff --git a/src/query/contact_manifolds/contact_manifolds_voxels_composite_shape.rs b/src/query/contact_manifolds/contact_manifolds_voxels_composite_shape.rs index 002eea9a..ed6118b2 100644 --- a/src/query/contact_manifolds/contact_manifolds_voxels_composite_shape.rs +++ b/src/query/contact_manifolds/contact_manifolds_voxels_composite_shape.rs @@ -135,7 +135,7 @@ pub fn contact_manifolds_voxels_composite_shape( timestamp: new_timestamp, }; - let vox_id = vox1.linear_id; + let vox_id = vox1.linear_id.flat_id() as u32; let (id1, id2) = if flipped { (leaf2, vox_id) } else { diff --git a/src/query/contact_manifolds/contact_manifolds_voxels_shape.rs b/src/query/contact_manifolds/contact_manifolds_voxels_shape.rs index 43fafaaf..6aaa57e1 100644 --- a/src/query/contact_manifolds/contact_manifolds_voxels_shape.rs +++ b/src/query/contact_manifolds/contact_manifolds_voxels_shape.rs @@ -171,7 +171,7 @@ pub fn contact_manifolds_voxels_shape( timestamp: new_timestamp, }; - let vid = vox1.linear_id; + let vid = vox1.linear_id.flat_id() as u32; let (id1, id2) = if flipped { (0, vid) } else { (vid, 0) }; manifolds.push(ContactManifold::with_data( id1, @@ -329,6 +329,7 @@ impl CanonicalVoxelShape { // detection). let mins = voxels.domain()[0] - Vector::repeat(1); let maxs = voxels.domain()[1]; + let counts = maxs - mins; let mask1 = vox.state.free_faces(); let adjust_canon = |axis: AxisMask, i: usize, key: &mut Point, val: i32| { @@ -348,12 +349,21 @@ impl CanonicalVoxelShape { adjust_canon(AxisMask::Z_NEG, 2, &mut key_low, mins[2]); } + #[cfg(feature = "dim2")] + let workspace_id = |vox: Point| { + let local = vox - mins; + (local.x + local.y * counts.x) as u32 + }; + + #[cfg(feature = "dim3")] + let workspace_id = |vox: Point| { + let local = vox - mins; + (local.x + local.y * counts.x + local.z * counts.x * counts.y) as u32 + }; + Self { range: [key_low, key_high], - workspace_key: Vector2::new( - voxels.linear_index(key_low), - voxels.linear_index(key_high), - ), + workspace_key: Vector2::new(workspace_id(key_low), workspace_id(key_high)), } } diff --git a/src/query/contact_manifolds/contact_manifolds_voxels_voxels.rs b/src/query/contact_manifolds/contact_manifolds_voxels_voxels.rs index 5a9d2da3..ae12eeb3 100644 --- a/src/query/contact_manifolds/contact_manifolds_voxels_voxels.rs +++ b/src/query/contact_manifolds/contact_manifolds_voxels_voxels.rs @@ -132,8 +132,8 @@ pub fn contact_manifolds_voxels_voxels<'a, ManifoldData, ContactData>( }; manifolds.push(ContactManifold::with_data( - vox1.linear_id, - vox2.linear_id, + vox1.linear_id.flat_id() as u32, + vox2.linear_id.flat_id() as u32, ManifoldData::default(), )); diff --git a/src/query/point/point_voxels.rs b/src/query/point/point_voxels.rs index 3fe39418..1c0d4f5f 100644 --- a/src/query/point/point_voxels.rs +++ b/src/query/point/point_voxels.rs @@ -1,38 +1,64 @@ use crate::math::{Point, Real}; use crate::query::{PointProjection, PointQuery}; -use crate::shape::{Cuboid, FeatureId, VoxelType, Voxels}; +use crate::shape::{Cuboid, FeatureId, Voxels, VoxelsChunkRef}; impl PointQuery for Voxels { #[inline] fn project_local_point(&self, pt: &Point, solid: bool) -> PointProjection { - // TODO: optimize this very naive implementation. - let base_cuboid = Cuboid::new(self.voxel_size() / 2.0); + self.chunk_bvh() + .project_point(pt, Real::MAX, |chunk_id, _| { + let chunk = self.chunk_ref(chunk_id); + Some(chunk.project_local_point_and_get_vox_id(pt, solid).0) + }) + .unwrap() + .1 + .1 + } + + #[inline] + fn project_local_point_and_get_feature( + &self, + pt: &Point, + ) -> (PointProjection, FeatureId) { + self.chunk_bvh() + .project_point_and_get_feature(pt, Real::MAX, |chunk_id, _| { + let chunk = self.chunk_ref(chunk_id); + let (proj, vox) = chunk.project_local_point_and_get_vox_id(pt, false); + // TODO: we need a way to return both the voxel id, and the feature on the voxel. + Some((proj, FeatureId::Face(vox))) + }) + .unwrap() + .1 + .1 + } +} + +impl<'a> VoxelsChunkRef<'a> { + #[inline] + fn project_local_point_and_get_vox_id( + &self, + pt: &Point, + solid: bool, + ) -> (PointProjection, u32) { + // TODO: optimize this naive implementation that just iterates on all the voxels + // from this chunk. + let base_cuboid = Cuboid::new(self.parent.voxel_size() / 2.0); let mut smallest_dist = Real::MAX; let mut result = PointProjection::new(false, *pt); + let mut result_vox_id = 0; for vox in self.voxels() { - if vox.state.voxel_type() != VoxelType::Empty { - let mut candidate = - base_cuboid.project_local_point(&(pt - vox.center.coords), solid); - candidate.point += vox.center.coords; + let mut candidate = base_cuboid.project_local_point(&(pt - vox.center.coords), solid); + candidate.point += vox.center.coords; - let candidate_dist = (candidate.point - pt).norm(); - if candidate_dist < smallest_dist { - result = candidate; - smallest_dist = candidate_dist; - } + let candidate_dist = (candidate.point - pt).norm(); + if candidate_dist < smallest_dist { + result = candidate; + result_vox_id = vox.linear_id.flat_id(); + smallest_dist = candidate_dist; } } - result - } - - #[inline] - fn project_local_point_and_get_feature( - &self, - pt: &Point, - ) -> (PointProjection, FeatureId) { - // TODO: get the actual feature. - (self.project_local_point(pt, false), FeatureId::Unknown) + (result, result_vox_id as u32) } } diff --git a/src/query/ray/ray_voxels.rs b/src/query/ray/ray_voxels.rs index aa6f9612..5077b30b 100644 --- a/src/query/ray/ray_voxels.rs +++ b/src/query/ray/ray_voxels.rs @@ -1,8 +1,30 @@ use crate::math::{Real, Vector}; +use crate::partitioning::BvhNode; use crate::query::{Ray, RayCast, RayIntersection}; -use crate::shape::{FeatureId, Voxels}; +use crate::shape::{FeatureId, Voxels, VoxelsChunkRef}; impl RayCast for Voxels { + #[inline] + fn cast_local_ray_and_get_normal( + &self, + ray: &Ray, + max_time_of_impact: Real, + solid: bool, + ) -> Option { + self.chunk_bvh() + .find_best( + max_time_of_impact, + |node: &BvhNode, best_so_far| node.cast_ray(ray, best_so_far), + |primitive, best_so_far| { + let chunk = self.chunk_ref(primitive); + chunk.cast_local_ray_and_get_normal(ray, best_so_far, solid) + }, + ) + .map(|(_chunk_id, hit)| hit) + } +} + +impl<'a> RayCast for VoxelsChunkRef<'a> { #[inline] fn cast_local_ray_and_get_normal( &self, @@ -31,19 +53,22 @@ impl RayCast for Voxels { let [domain_mins, domain_maxs] = self.domain(); loop { - let voxel = self.voxel_state(voxel_key); - let aabb = self.voxel_aabb(voxel_key); + let aabb = self.voxel_aabb_unchecked(voxel_key); - if !voxel.is_empty() { - // We hit a voxel! - // TODO: if `solid` is false, and we started hitting from the first iteration, - // then we should continue the ray propagation until we reach empty space again. - let hit = aabb.cast_local_ray_and_get_normal(ray, max_t, solid); + if let Some(voxel) = self.voxel_state(voxel_key) { + if !voxel.is_empty() { + // We hit a voxel! + // TODO: if `solid` is false, and we started hitting from the first iteration, + // then we should continue the ray propagation until we reach empty space again. + let hit = aabb.cast_local_ray_and_get_normal(ray, max_t, solid); - if let Some(mut hit) = hit { - // TODO: have the feature id be based on the voxel type? - hit.feature = FeatureId::Face(self.linear_index(voxel_key)); - return Some(hit); + if let Some(mut hit) = hit { + // TODO: have the feature id be based on the voxel type? + hit.feature = FeatureId::Face( + self.flat_id(voxel_key).unwrap_or_else(|| unreachable!()), + ); + return Some(hit); + } } } diff --git a/src/shape/mod.rs b/src/shape/mod.rs index 93241e55..156a3458 100644 --- a/src/shape/mod.rs +++ b/src/shape/mod.rs @@ -22,7 +22,7 @@ pub use self::{ compound::Compound, polyline::Polyline, shared_shape::SharedShape, - voxels::{AxisMask, OctantPattern, VoxelData, VoxelState, VoxelType, Voxels}, + voxels::{AxisMask, OctantPattern, VoxelData, VoxelState, VoxelType, Voxels, VoxelsChunkRef}, }; #[cfg(feature = "dim2")] diff --git a/src/shape/voxels.rs b/src/shape/voxels.rs deleted file mode 100644 index 7b7fc1fc..00000000 --- a/src/shape/voxels.rs +++ /dev/null @@ -1,1555 +0,0 @@ -use crate::bounding_volume::Aabb; -use crate::math::{Point, Real, Vector, DIM}; -use alloc::{vec, vec::Vec}; -#[cfg(not(feature = "std"))] -use na::ComplexField; - -/// Categorization of a voxel based on its neighbors. -#[derive(Copy, Clone, Debug, PartialEq, Eq)] -pub enum VoxelType { - /// The voxel is empty. - Empty, - /// The voxel is a vertex if all three coordinate axis directions have at - /// least one empty neighbor. - Vertex, - /// The voxel is on an edge if it has non-empty neighbors in both directions of - /// a single coordinate axis. - #[cfg(feature = "dim3")] - Edge, - /// The voxel is on an edge if it has non-empty neighbors in both directions of - /// two coordinate axes. - Face, - /// The voxel is on an edge if it has non-empty neighbors in both directions of - /// all coordinate axes. - Interior, -} - -#[derive(Clone, Copy, Debug, Default, Eq, Hash, Ord, PartialEq, PartialOrd)] -/// The status of the cell of an heightfield. -pub struct AxisMask(u8); - -bitflags::bitflags! { - /// Flags for identifying signed directions along coordinate axes, or faces of a voxel. - impl AxisMask: u8 { - /// The direction or face along the `+x` coordinate axis. - const X_POS = 1 << 0; - /// The direction or face along the `-x` coordinate axis. - const X_NEG = 1 << 1; - /// The direction or face along the `+y` coordinate axis. - const Y_POS = 1 << 2; - /// The direction or face along the `-y` coordinate axis. - const Y_NEG = 1 << 3; - /// The direction or face along the `+z` coordinate axis. - #[cfg(feature= "dim3")] - const Z_POS = 1 << 4; - /// The direction or face along the `-z` coordinate axis. - #[cfg(feature= "dim3")] - const Z_NEG = 1 << 5; - } -} - -/// Indicates the local shape of a voxel on each octant. -/// -/// This provides geometric information of the shape’s exposed features on each octant. -// This is an alternative to `FACES_TO_FEATURE_MASKS` that can be more convenient for some -// collision-detection algorithms. -#[derive(Copy, Clone, Debug, PartialEq, Eq)] -pub struct OctantPattern; - -// NOTE: it is important that the max value of any OctantPattern variant -// is 7 because we don’t allocate more than 3 bits to store it in -// `FACES_TO_OCTANT_MASKS`. -/// Indicates the local shape of a voxel on each octant. -/// -/// This provides geometric information of the shape’s exposed features on each octant. -// This is an alternative to `FACES_TO_FEATURE_MASKS` that can be more convenient for some -// collision-detection algorithms. -#[cfg(feature = "dim3")] -impl OctantPattern { - /// The voxel doesn't have any exposed feature on the octant with this mask. - pub const INTERIOR: u32 = 0; - /// The voxel has an exposed vertex on the octant with this mask. - pub const VERTEX: u32 = 1; - /// The voxel has an exposed edges with direction X on the octant with this mask. - pub const EDGE_X: u32 = 2; - /// The voxel has an exposed edges with direction Y on the octant with this mask. - pub const EDGE_Y: u32 = 3; - /// The voxel has an exposed edges with direction Z on the octant with this mask. - pub const EDGE_Z: u32 = 4; - /// The voxel has an exposed face with normal X on the octant with this mask. - pub const FACE_X: u32 = 5; - /// The voxel has an exposed face with normal Y on the octant with this mask. - pub const FACE_Y: u32 = 6; - /// The voxel has an exposed face with normal Z on the octant with this mask. - pub const FACE_Z: u32 = 7; -} - -// NOTE: it is important that the max value of any OctantPattern variant -// is 7 because we don’t allocate more than 3 bits to store it in -// `FACES_TO_OCTANT_MASKS`. -/// Indicates the local shape of a voxel on each octant. -/// -/// This provides geometric information of the shape’s exposed features on each octant. -/// This is an alternative to `FACES_TO_FEATURE_MASKS` that can be more convenient for some -/// collision-detection algorithms. -#[cfg(feature = "dim2")] -impl OctantPattern { - /// The voxel doesn't have any exposed feature on the octant with this mask. - pub const INTERIOR: u32 = 0; - /// The voxel has an exposed vertex on the octant with this mask. - pub const VERTEX: u32 = 1; - /// The voxel has an exposed face with normal X on the octant with this mask. - pub const FACE_X: u32 = 2; - /// The voxel has an exposed face with normal Y on the octant with this mask. - pub const FACE_Y: u32 = 3; -} - -// The local neighborhood information is encoded in a 8-bits number in groups of two bits -// per coordinate axis: `0bwwzzyyxx`. In each group of two bits, e.g. `xx`, the rightmost (resp. -// leftmost) bit set to 1 means that the neighbor voxel with coordinate `+1` (resp `-1`) relative -// to the current voxel along the `x` axis is filled. If the bit is 0, then the corresponding -// neighbor is empty. See the `AxisMask` bitflags. -// For example, in 2D, the mask `0b00_00_10_01` matches the following configuration (assuming +y goes -// up, and +x goes right): -// -// ```txt -// 0 0 0 -// 0 x 1 -// 0 1 0 -// ``` -// -// The special value `0b01000000` indicates that the voxel is empty. -// And the value `0b00111111` (`0b00001111` in 2D) indicates that the voxel is an interior voxel (its whole neighborhood -// is filled). -/// A description of the local neighborhood of a voxel. -#[derive(Copy, Clone, Debug, PartialEq, Eq)] -#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))] -pub struct VoxelState(u8); - -impl VoxelState { - /// The value of empty voxels. - pub const EMPTY: VoxelState = VoxelState(EMPTY_FACE_MASK); - /// The value of a voxel with non-empty neighbors in all directions. - pub const INTERIOR: VoxelState = VoxelState(INTERIOR_FACE_MASK); - - /// Is this voxel empty? - pub const fn is_empty(self) -> bool { - self.0 == EMPTY_FACE_MASK - } - - /// A bit mask indicating which faces of the voxel don’t have any - /// adjacent non-empty voxel. - pub const fn free_faces(self) -> AxisMask { - if self.0 == INTERIOR_FACE_MASK || self.0 == EMPTY_FACE_MASK { - AxisMask::empty() - } else { - AxisMask::from_bits_truncate((!self.0) & INTERIOR_FACE_MASK) - } - } - - /// The [`VoxelType`] of this voxel. - pub const fn voxel_type(self) -> VoxelType { - FACES_TO_VOXEL_TYPES[self.0 as usize] - } - - // Bitmask indicating what vertices, edges, or faces of the voxel are "free". - pub(crate) const fn feature_mask(self) -> u16 { - FACES_TO_FEATURE_MASKS[self.0 as usize] - } - - pub(crate) const fn octant_mask(self) -> u32 { - FACES_TO_OCTANT_MASKS[self.0 as usize] - } -} - -/// Information associated to a voxel. -#[derive(Copy, Clone, Debug, PartialEq)] -pub struct VoxelData { - /// The temporary index in the internal voxels’ storage. - /// - /// This index can be invalidated after a call to [`Voxels::set_voxel`], or - /// [`Voxels::resize_domain`]. - pub linear_id: u32, - /// The voxel’s integer grid coordinates. - pub grid_coords: Point, - /// The voxel’s center position in the local-space of the [`Voxels`] shape it is part of. - pub center: Point, - /// The voxel’s state, indicating if it’s empty or full. - pub state: VoxelState, -} - -/// A shape made of axis-aligned, uniformly sized, cubes (aka. voxels). -/// -/// This shape is specialized to handle voxel worlds and voxelized obojects efficiently why ensuring -/// that collision-detection isn’t affected by the so-called "internal edges problem" that can create -/// artifacts when another objects rolls or slides against a flat voxelized surface. -/// -/// The internal storage is compact (but not sparse at the moment), storing only one byte per voxel -/// in the allowed domain. This has a generally smaller memory footprint than a mesh representation -/// of the voxels. -#[derive(Clone, Debug)] -#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))] -pub struct Voxels { - domain_mins: Point, - domain_maxs: Point, - states: Vec, // Somehow switch to a sparse representation? - voxel_size: Vector, -} - -impl Voxels { - /// Initializes a voxel shapes from the voxels grid coordinates. - /// - /// Each voxel will have its bottom-left-back corner located at - /// `grid_coordinates * voxel_size`; and its center at `(grid_coordinates + 0.5) * voxel_size`. - pub fn new(voxel_size: Vector, grid_coordinates: &[Point]) -> Self { - let mut domain_mins = grid_coordinates[0]; - let mut domain_maxs = grid_coordinates[0]; - - for vox in grid_coordinates { - domain_mins = domain_mins.inf(vox); - domain_maxs = domain_maxs.sup(vox); - } - - domain_maxs += Vector::repeat(1); - let dimensions = domain_maxs - domain_mins; - let voxels_count = dimensions.product(); - let mut result = Self { - domain_mins, - domain_maxs, - states: vec![VoxelState::EMPTY; voxels_count as usize], - voxel_size, - }; - - for vox in grid_coordinates { - let index = result.linear_index(*vox); - result.states[index as usize] = VoxelState::INTERIOR; - } - - result.recompute_voxels_data(); - result - } - - /// Computes a voxels shape from the set of `points`. - /// - /// The points are mapped to a regular grid centered at the provided point with smallest - /// coordinates, and with grid cell size equal to `scale`. It is OK if multiple points - /// fall into the same grid cell. - pub fn from_points(voxel_size: Vector, points: &[Point]) -> Self { - let voxels: Vec<_> = points - .iter() - .map(|pt| { - Point::from( - pt.coords - .component_div(&voxel_size) - .map(|x| x.floor() as i32), - ) - }) - .collect(); - Self::new(voxel_size, &voxels) - } - - // TODO: support a crate like get_size2 (will require support on nalgebra too)? - /// An approximation of the memory usage (in bytes) for this struct plus - /// the memory it allocates dynamically. - pub fn total_memory_size(&self) -> usize { - size_of::() + self.heap_memory_size() - } - - /// An approximation of the memory dynamically-allocated by this struct. - pub fn heap_memory_size(&self) -> usize { - // NOTE: if a new field is added to `Self`, adjust this function result. - let Self { - domain_mins: _, - domain_maxs: _, - states: data, - voxel_size: _, - } = self; - data.capacity() * size_of::() - } - - /// The extents of the total axis-aligned volume covered by this [`Voxels`] shape. - /// - /// This accounts for all the voxels reserved in the internal buffer of `self`, including empty - /// ones. - pub fn extents(&self) -> Vector { - self.dimensions() - .cast::() - .component_mul(&self.voxel_size) - } - - /// The center of this shape’s domain (accounting for both empty and filled voxels). - pub fn domain_center(&self) -> Point { - (self - .domain_mins - .coords - .cast::() - .component_mul(&self.voxel_size) - + self.extents() / 2.0) - .into() - } - - /// Sets the size of each voxel along each local coordinate axis. - pub fn set_voxel_size(&mut self, size: Vector) { - self.voxel_size = size; - } - - /// The valid semi-open range of voxel grid indices. - /// - /// With `let [mins, maxs] = voxels.domain();` the valid indices along the dimension `i` are - /// all the indices in the range `mins[i]..maxs[i]` (i.e. `maxs[i]` is excluded). - pub fn domain(&self) -> [&Point; 2] { - [&self.domain_mins, &self.domain_maxs] - } - - /// The number of voxels along each coordinate axis. - pub fn dimensions(&self) -> Vector { - (self.domain_maxs - self.domain_mins).map(|e| e as u32) - } - - /// The size of each voxel part this [`Voxels`] shape. - pub fn voxel_size(&self) -> Vector { - self.voxel_size - } - - /// Merges voxel state (neighborhood) information of a given voxel (and all its neighbors) - /// from `self` and `other`, to account for a recent change to the given `voxel` in `self`. - /// - /// This is designed to be called after `self` was modified with [`Voxels::set_voxel`] or - /// [`Voxels::try_set_voxel`]. - /// - /// This is the same as [`Voxels::combine_voxel_states`] but localized to a single voxel and its - /// neighbors. - pub fn propagate_voxel_change( - &mut self, - other: &mut Self, - voxel: Point, - origin_shift: Vector, - ) { - let center_is_empty = self - .get_voxel_state(voxel) - .map(|vox| vox.is_empty()) - .unwrap_or(true); - let center_state_delta = - other.update_neighbors_state(voxel - origin_shift, center_is_empty); - - if let Some(state_id) = self.get_linear_index(voxel) { - self.states[state_id as usize].0 |= center_state_delta.0; - } - } - - /// Merges voxel state (neighborhood) information of each voxel from `self` and `other`. - /// - /// This allows each voxel from one shape to be aware of the presence of neighbors belonging to - /// the other so that collision detection is capable of transitioning between the boundaries of - /// one shape to the other without hitting an internal edge. - /// - /// Both voxels shapes are assumed to have the same [`Self::voxel_size`]. - /// If `other` lives in a coordinate space with a different origin than `self`, then - /// `origin_shift` represents the distance (as a multiple of the `voxel_size`) from the origin - /// of `self` to the origin of `other`. Therefore, a voxel with coordinates `key` on `other` - /// will have coordinates `key + origin_shift` on `self`. - pub fn combine_voxel_states(&mut self, other: &mut Self, origin_shift: Vector) { - let one = Vector::repeat(1); - - // Intersect the domains + 1 cell. - let d0 = [self.domain_mins - one, self.domain_maxs + one * 2]; - let d1 = [ - other.domain_mins - one + origin_shift, - other.domain_maxs + one * 2 + origin_shift, - ]; - - let d01 = [d0[0].sup(&d1[0]), d0[1].inf(&d1[1])]; - // Iterate on the domain intersection. If the voxel exists (and is non-empty) on both shapes, we - // simply need to combine their bitmasks. If it doesn’t exist on both shapes, we need to - // actually check the neighbors. - // - // The `domain` is expressed in the grid coordinate space of `self`. - for i in d01[0].x..d01[1].x { - for j in d01[0].y..d01[1].y { - #[cfg(feature = "dim2")] - let k_range = 0..1; - #[cfg(feature = "dim3")] - let k_range = d01[0].z..d01[1].z; - for _k in k_range { - #[cfg(feature = "dim2")] - let key0 = Point::new(i, j); - #[cfg(feature = "dim3")] - let key0 = Point::new(i, j, _k); - let key1 = key0 - origin_shift; - let id0 = self - .get_linear_index(key0) - .filter(|id| !self.states[*id as usize].is_empty()); - let id1 = other - .get_linear_index(key1) - .filter(|id| !other.states[*id as usize].is_empty()); - - match (id0, id1) { - (Some(id0), Some(id1)) => { - self.states[id0 as usize].0 |= other.states[id1 as usize].0; - other.states[id1 as usize].0 |= self.states[id0 as usize].0; - } - (Some(id0), None) => { - self.states[id0 as usize].0 |= - other.compute_voxel_neighborhood_bits(key1).0; - } - (None, Some(id1)) => { - other.states[id1 as usize].0 |= - self.compute_voxel_neighborhood_bits(key0).0; - } - (None, None) => { /* Nothing to adjust. */ } - } - } - } - } - } - - fn recompute_voxels_data(&mut self) { - for i in 0..self.states.len() { - let key = self.voxel_at_id(i as u32); - self.states[i] = self.compute_voxel_state(key); - } - } - - /// Scale this shape. - pub fn scaled(mut self, scale: &Vector) -> Self { - self.voxel_size.component_mul_assign(scale); - self - } - - /// Sets the voxel at the given grid coordinates, returning `None` if it lies outside [`Self::domain`]. - /// - /// See [`Self::set_voxel`] for a method that automatically resizes the internal - /// storage of `self` if the key is out of the valid bounds. - pub fn try_set_voxel(&mut self, key: Point, is_filled: bool) -> Option { - let id = self.get_linear_index(key)?; - let prev = self.states[id as usize]; - let new_is_empty = !is_filled; - - if prev.is_empty() ^ new_is_empty { - self.states[id as usize] = self.update_neighbors_state(key, new_is_empty); - } - - Some(prev) - } - - /// Inserts a voxel at the given key, even if it is out of the bounds of this shape. - /// - /// If `is_filed` is `true` and the key lies out of the bounds on this shape, the internal - /// voxels storage will be resized automatically. If a resize happens, the cost of the insertion - /// is `O(n)` where `n` is the capacity of `self`. If no resize happens, then the cost of - /// insertion is `O(1)`. - /// - /// Use [`Self::try_set_voxel`] instead for a version that will be a no-op if the provided - /// coordinates are outside the [`Self::domain`], avoiding potential internal reallocations. - pub fn set_voxel(&mut self, key: Point, is_filled: bool) -> Option { - if !self.is_voxel_in_bounds(key) && is_filled { - let dims = self.dimensions(); - - // Add 10% extra padding. - let extra = dims.map(|k| k * 10 / 100); - let mut new_domain_mins = self.domain_mins; - let mut new_domain_maxs = self.domain_maxs; - - for k in 0..DIM { - if key[k] < self.domain_mins[k] { - new_domain_mins[k] = key[k] - extra[k] as i32; - } - - if key[k] >= self.domain_maxs[k] { - new_domain_maxs[k] = key[k] + extra[k] as i32 + 1; - } - } - - self.resize_domain(new_domain_mins, new_domain_maxs); - - self.try_set_voxel(key, is_filled) - } else { - self.try_set_voxel(key, is_filled) - } - } - - /// Set the model domain. - /// - /// The domain can be smaller or larger than the current one. Voxels in any overlap between - /// the current and new domain will be preserved. - /// - /// If for any index `i`, `domain_maxs[i] <= domain_mins[i]`, then the new domain is invalid - /// and this operation will result in a no-op. - pub fn resize_domain(&mut self, domain_mins: Point, domain_maxs: Point) { - if self.domain_mins == domain_mins && self.domain_maxs == domain_maxs { - // Nothing to change. - return; - } - - if let Some(new_shape) = self.with_resized_domain(domain_mins, domain_maxs) { - *self = new_shape; - } - } - - /// Clone this voxels shape with a new size. - /// - /// The domain can be smaller or larger than the current one. Voxels in any overlap between - /// the current and new domain will be preserved. - /// - /// If for any index `i`, `domain_maxs[i] <= domain_mins[i]`, then the new domain is invalid - /// and this operation returns `None`. - #[must_use] - pub fn with_resized_domain( - &self, - domain_mins: Point, - domain_maxs: Point, - ) -> Option { - if self.domain_mins == domain_mins && self.domain_maxs == domain_maxs { - // Nothing to change, just clone as-is. - return Some(self.clone()); - } - - let new_dim = domain_maxs - domain_mins; - if new_dim.iter().any(|d| *d <= 0) { - log::error!("Invalid domain provided for resizing a voxels shape. New domain: {domain_mins:?} - {domain_maxs:?}; new domain size: {new_dim:?}"); - return None; - } - - let new_len = new_dim.iter().map(|x| *x as usize).product(); - - let mut new_shape = Self { - domain_mins, - domain_maxs, - states: vec![VoxelState::EMPTY; new_len], - voxel_size: self.voxel_size, - }; - - for i in 0..self.states.len() { - let key = self.voxel_at_id(i as u32); - let new_i = new_shape.linear_index(key); - new_shape.states[new_i as usize] = self.states[i]; - } - - Some(new_shape) - } - - /// Checks if the given key is within [`Self::domain`]. - #[cfg(feature = "dim2")] - pub fn is_voxel_in_bounds(&self, key: Point) -> bool { - key[0] >= self.domain_mins[1] - && key[0] < self.domain_maxs[0] - && key[1] >= self.domain_mins[1] - && key[1] < self.domain_maxs[1] - } - - /// Checks if the given key is within [`Self::domain`]. - #[cfg(feature = "dim3")] - pub fn is_voxel_in_bounds(&self, key: Point) -> bool { - key[0] >= self.domain_mins[0] - && key[0] < self.domain_maxs[0] - && key[1] >= self.domain_mins[1] - && key[1] < self.domain_maxs[1] - && key[2] >= self.domain_mins[2] - && key[2] < self.domain_maxs[2] - } - - /// Updates the state of the neighbors of the voxel `key`. - /// - /// Modifies the state of the neighbors of `key` to account for it being empty or full. - /// Returns (but doesn’t modify) the new state of the voxel specified by `key`. - #[must_use] - fn update_neighbors_state(&mut self, key: Point, center_is_empty: bool) -> VoxelState { - let mut key_data = 0; - - for k in 0..DIM { - let mut left = key; - let mut right = key; - left[k] -= 1; - right[k] += 1; - - if let Some(left_id) = self.get_linear_index(left) { - if !self.states[left_id as usize].is_empty() { - if center_is_empty { - self.states[left_id as usize].0 &= !(1 << (k * 2)); - } else { - self.states[left_id as usize].0 |= 1 << (k * 2); - key_data |= 1 << (k * 2 + 1); - } - } - } - - if let Some(right_id) = self.get_linear_index(right) { - if !self.states[right_id as usize].is_empty() { - if center_is_empty { - self.states[right_id as usize].0 &= !(1 << (k * 2 + 1)); - } else { - self.states[right_id as usize].0 |= 1 << (k * 2 + 1); - key_data |= 1 << (k * 2); - } - } - } - } - - if center_is_empty { - VoxelState::EMPTY - } else { - VoxelState(key_data) - } - } - - /// The AABB of the voxel with the given quantized `key`. - pub fn voxel_aabb(&self, key: Point) -> Aabb { - let center = self.voxel_center(key); - let hext = self.voxel_size / 2.0; - Aabb::from_half_extents(center, hext) - } - - /// Returns the state of a given voxel. - /// - /// Panics if the key is out of the bounds defined by [`Self::domain`]. - pub fn voxel_state(&self, key: Point) -> VoxelState { - self.states[self.linear_index(key) as usize] - } - - /// Returns the state of a given voxel, if it is within the bounds of [`Self::domain`]. - pub fn get_voxel_state(&self, key: Point) -> Option { - Some(self.states[self.get_linear_index(key)? as usize]) - } - - /// Calculates the grid coordinates of the voxel containing the given `point`, regardless - /// of [`Self::domain`]. - pub fn voxel_at_point_unchecked(&self, point: Point) -> Point { - point - .coords - .component_div(&self.voxel_size) - .map(|x| x.floor() as i32) - .into() - } - - /// Gets the key of the voxel containing the given `pt`. - /// - /// Note that the returned key might address a voxel that is empty. - /// `None` is returned if the point is out of the domain of `self`. - pub fn voxel_at_point(&self, pt: Point) -> Option> { - let quant = self.voxel_at_point_unchecked(pt); - if quant[0] < self.domain_mins[0] - || quant[1] < self.domain_mins[1] - || quant[0] >= self.domain_maxs[0] - || quant[1] >= self.domain_maxs[1] - { - return None; - } - - #[cfg(feature = "dim3")] - if quant[2] < self.domain_mins[2] || quant[2] >= self.domain_maxs[2] { - return None; - } - - Some(quant) - } - - /// Clamps an arbitrary voxel into the valid domain of `self`. - pub fn clamp_voxel(&self, key: Point) -> Point { - key.coords - .zip_zip_map( - &self.domain_mins.coords, - &self.domain_maxs.coords, - |k, min, max| k.clamp(min, max - 1), - ) - .into() - } - - /// The range of grid coordinates of voxels intersecting the given AABB. - /// - /// The returned range covers both empty and non-empty voxels, and is not limited to the - /// bounds defined by [`Self::domain`]. - /// The range is semi, open, i.e., the range along each dimension `i` is understood as - /// the semi-open interval: `range[0][i]..range[1][i]`. - pub fn voxel_range_intersecting_local_aabb(&self, aabb: &Aabb) -> [Point; 2] { - let mins = aabb - .mins - .coords - .component_div(&self.voxel_size) - .map(|x| x.floor() as i32); - let maxs = aabb - .maxs - .coords - .component_div(&self.voxel_size) - .map(|x| x.ceil() as i32); - [mins.into(), maxs.into()] - } - - /// The AABB of a given range of voxels. - /// - /// The AABB is computed independently of [`Self::domain`] and independently of whether - /// the voxels contained within are empty or not. - pub fn voxel_range_aabb(&self, mins: Point, maxs: Point) -> Aabb { - Aabb { - mins: mins - .cast::() - .coords - .component_mul(&self.voxel_size) - .into(), - maxs: maxs - .cast::() - .coords - .component_mul(&self.voxel_size) - .into(), - } - } - - /// Aligns the given AABB with the voxelized grid. - /// - /// The aligned is calculated such that the returned AABB has corners lying at the grid - /// intersections (i.e. matches voxel corners) and fully contains the input `aabb`. - pub fn align_aabb_to_grid(&self, aabb: &Aabb) -> Aabb { - let mins = aabb - .mins - .coords - .zip_map(&self.voxel_size, |m, sz| (m / sz).floor() * m) - .into(); - let maxs = aabb - .maxs - .coords - .zip_map(&self.voxel_size, |m, sz| (m / sz).ceil() * m) - .into(); - Aabb { mins, maxs } - } - - /// Iterates through every voxel intersecting the given aabb. - /// - /// Returns the voxel’s linearized id, center, and state. - pub fn voxels_intersecting_local_aabb( - &self, - aabb: &Aabb, - ) -> impl Iterator + '_ { - let [mins, maxs] = self.voxel_range_intersecting_local_aabb(aabb); - self.voxels_in_range(mins, maxs) - } - - /// The center point of all the voxels in this shape (including empty ones). - /// - /// The voxel data associated to each center is provided to determine what kind of voxel - /// it is (and, in particular, if it is empty or full). - pub fn voxels(&self) -> impl Iterator + '_ { - self.voxels_in_range(self.domain_mins, self.domain_maxs) - } - - /// Splits this voxels shape into two subshapes. - /// - /// The first subshape contains all the voxels which centers are inside the `aabb`. - /// The second subshape contains all the remaining voxels. - pub fn split_with_box(&self, aabb: &Aabb) -> (Option, Option) { - // TODO: optimize this? - let mut in_box = vec![]; - let mut rest = vec![]; - for vox in self.voxels() { - if !vox.state.is_empty() { - if aabb.contains_local_point(&vox.center) { - in_box.push(vox.center); - } else { - rest.push(vox.center); - } - } - } - - let in_box = if !in_box.is_empty() { - Some(Voxels::from_points(self.voxel_size, &in_box)) - } else { - None - }; - - let rest = if !rest.is_empty() { - Some(Voxels::from_points(self.voxel_size, &rest)) - } else { - None - }; - - (in_box, rest) - } - - /// Iterate through the data of all the voxels within the given (semi-open) voxel grid indices. - /// - /// Note that this yields both empty and non-empty voxels within the range. This does not - /// include any voxel that falls outside [`Self::domain`]. - #[cfg(feature = "dim2")] - pub fn voxels_in_range( - &self, - mins: Point, - maxs: Point, - ) -> impl Iterator + '_ { - let mins = mins.coords.sup(&self.domain_mins.coords); - let maxs = maxs.coords.inf(&self.domain_maxs.coords); - - (mins[0]..maxs[0]).flat_map(move |ix| { - (mins[1]..maxs[1]).map(move |iy| { - let grid_coords = Point::new(ix, iy); - let vid = self.linear_index(grid_coords); - let center = - Vector::new(ix as Real + 0.5, iy as Real + 0.5).component_mul(&self.voxel_size); - VoxelData { - linear_id: vid, - grid_coords, - center: center.into(), - state: self.states[vid as usize], - } - }) - }) - } - - /// Iterate through the data of all the voxels within the given (semi-open) voxel grid indices. - /// - /// Note that this yields both empty and non-empty voxels within the range. This does not - /// include any voxel that falls outside [`Self::domain`]. - #[cfg(feature = "dim3")] - pub fn voxels_in_range( - &self, - mins: Point, - maxs: Point, - ) -> impl Iterator + '_ { - let mins = mins.coords.sup(&self.domain_mins.coords); - let maxs = maxs.coords.inf(&self.domain_maxs.coords); - - (mins[0]..maxs[0]).flat_map(move |ix| { - (mins[1]..maxs[1]).flat_map(move |iy| { - (mins[2]..maxs[2]).map(move |iz| { - let grid_coords = Point::new(ix, iy, iz); - let vid = self.linear_index(grid_coords); - let center = Vector::new(ix as Real + 0.5, iy as Real + 0.5, iz as Real + 0.5) - .component_mul(&self.voxel_size) - .into(); - VoxelData { - linear_id: vid, - grid_coords, - center, - state: self.states[vid as usize], - } - }) - }) - }) - } - - /// Gets linearized index associated to the given voxel key if it lies within [`Self::domain`]. - pub fn get_linear_index(&self, key: Point) -> Option { - if key[0] < self.domain_mins[0] - || key[0] >= self.domain_maxs[0] - || key[1] < self.domain_mins[1] - || key[1] >= self.domain_maxs[1] - { - return None; - } - - #[cfg(feature = "dim3")] - if key[2] < self.domain_mins[2] || key[2] >= self.domain_maxs[2] { - return None; - } - - Some(self.linear_index(key)) - } - - /// The linearized index associated to the given voxel key. - #[cfg(feature = "dim2")] - pub fn linear_index(&self, voxel_key: Point) -> u32 { - let dims = self.dimensions(); - let rel_key = voxel_key - self.domain_mins; - (rel_key.x + rel_key.y * dims[0] as i32) as u32 - } - - /// The linearized index associated to the given voxel key. - #[cfg(feature = "dim3")] - pub fn linear_index(&self, voxel_key: Point) -> u32 { - let dims = self.dimensions(); - let rel_key = voxel_key - self.domain_mins; - rel_key.x as u32 + rel_key.y as u32 * dims[0] + rel_key.z as u32 * dims[0] * dims[1] - } - - /// The key of the voxel at the given linearized index. - #[cfg(feature = "dim2")] - pub fn voxel_at_id(&self, linear_index: u32) -> Point { - let dim0 = self.domain_maxs[0] - self.domain_mins[0]; - let y = linear_index as i32 / dim0; - let x = linear_index as i32 % dim0; - self.domain_mins + Vector::new(x, y) - } - - /// The key of the voxel at the given linearized index. - #[cfg(feature = "dim3")] - pub fn voxel_at_id(&self, linear_index: u32) -> Point { - let dims = self.dimensions(); - - let d0d1 = dims[0] * dims[1]; - let z = linear_index / d0d1; - let y = (linear_index - z * d0d1) / dims[0]; - let x = linear_index % dims[0]; - self.domain_mins + Vector::new(x as i32, y as i32, z as i32) - } - - /// The center of the voxel with the given key. - pub fn voxel_center(&self, key: Point) -> Point { - (key.cast::() + Vector::repeat(0.5)) - .coords - .component_mul(&self.voxel_size) - .into() - } - - fn compute_voxel_state(&self, key: Point) -> VoxelState { - if self.states[self.linear_index(key) as usize].is_empty() { - return VoxelState::EMPTY; - } - - self.compute_voxel_neighborhood_bits(key) - } - - fn compute_voxel_neighborhood_bits(&self, key: Point) -> VoxelState { - let mut occupied_faces = 0; - - for k in 0..DIM { - let (mut prev, mut next) = (key, key); - prev[k] -= 1; - next[k] += 1; - - if let Some(next_id) = self.get_linear_index(next) { - if !self.states[next_id as usize].is_empty() { - occupied_faces |= 1 << (k * 2); - } - } - if let Some(prev_id) = self.get_linear_index(prev) { - if !self.states[prev_id as usize].is_empty() { - occupied_faces |= 1 << (k * 2 + 1); - } - } - } - - VoxelState(occupied_faces) - } -} - -// NOTE: this code is used to generate the constant tables -// FACES_TO_VOXEL_TYPES, FACES_TO_FEATURE_MASKS, FACES_TO_OCTANT_MASKS. -#[allow(dead_code)] -#[cfg(feature = "dim2")] -#[cfg(test)] -fn gen_const_tables() { - // The `j-th` bit of `faces_adj_to_vtx[i]` is set to 1, if the j-th face of the AABB (based on - // the face order depicted in `AABB::FACES_VERTEX_IDS`) is adjacent to the `i` vertex of the AABB - // (vertices are indexed as per the diagram depicted in the `FACES_VERTEX_IDS` doc. - // Each entry of this will always have exactly 3 bits set. - let mut faces_adj_to_vtx = [0usize; 4]; - - for fid in 0..4 { - let vids = Aabb::FACES_VERTEX_IDS[fid]; - let key = 1 << fid; - faces_adj_to_vtx[vids.0] |= key; - faces_adj_to_vtx[vids.1] |= key; - } - - /* - * FACES_TO_VOXEL_TYPES - */ - std::println!("const FACES_TO_VOXEL_TYPES: [VoxelType; 17] = ["); - 'outer: for i in 0usize..16 { - // If any vertex of the voxel has three faces with no adjacent voxels, - // then the voxel type is Vertex. - for adjs in faces_adj_to_vtx.iter() { - if (*adjs & i) == 0 { - std::println!("VoxelType::Vertex,"); - continue 'outer; - } - } - - // If one face doesn’t have any adjacent voxel, - // then the voxel type is Face. - for fid in 0..4 { - if ((1 << fid) & i) == 0 { - std::println!("VoxelType::Face,"); - continue 'outer; - } - } - } - - // Add final entries for special values. - std::println!("VoxelType::Interior,"); - std::println!("VoxelType::Empty,"); - std::println!("];"); - - /* - * FACES_TO_FEATURE_MASKS - */ - std::println!("const FACES_TO_FEATURE_MASKS: [u16; 17] = ["); - for i in 0usize..16 { - // Each bit set indicates a convex vertex that can lead to collisions. - // The result will be nonzero only for `VoxelType::Vertex` voxels. - let mut vtx_key = 0; - for (vid, adjs) in faces_adj_to_vtx.iter().enumerate() { - if (*adjs & i) == 0 { - vtx_key |= 1 << vid; - } - } - - if vtx_key != 0 { - std::println!("0b{:b},", vtx_key as u16); - continue; - } - - // Each bit set indicates an exposed face that can lead to collisions. - // The result will be nonzero only for `VoxelType::Face` voxels. - let mut face_key = 0; - for fid in 0..4 { - if ((1 << fid) & i) == 0 { - face_key |= 1 << fid; - } - } - - if face_key != 0 { - std::println!("0b{:b},", face_key as u16); - continue; - } - } - - std::println!("0b{:b},", u16::MAX); - std::println!("0,"); - std::println!("];"); - - /* - * Faces to octant masks. - */ - std::println!("const FACES_TO_OCTANT_MASKS: [u32; 17] = ["); - for i in 0usize..16 { - // First test if we have vertices. - let mut octant_mask = 0; - let mut set_mask = |mask, octant| { - // NOTE: we don’t overwrite any mask already set for the octant. - if (octant_mask >> (octant * 3)) & 0b0111 == 0 { - octant_mask |= mask << (octant * 3); - } - }; - - for (vid, adjs) in faces_adj_to_vtx.iter().enumerate() { - if (*adjs & i) == 0 { - set_mask(1, vid); - } - } - - // This is the index of the normal of the faces given by - // Aabb::FACES_VERTEX_IDS. - const FX: u32 = OctantPattern::FACE_X; - const FY: u32 = OctantPattern::FACE_Y; - const FACE_NORMALS: [u32; 4] = [FX, FX, FY, FY]; - - #[allow(clippy::needless_range_loop)] - for fid in 0..4 { - if ((1 << fid) & i) == 0 { - let vid = Aabb::FACES_VERTEX_IDS[fid]; - let mask = FACE_NORMALS[fid]; - - set_mask(mask, vid.0); - set_mask(mask, vid.1); - } - } - std::println!("0b{:b},", octant_mask); - } - std::println!("0,"); - std::println!("];"); -} - -// NOTE: this code is used to generate the constant tables -// FACES_TO_VOXEL_TYPES, FACES_TO_FEATURE_MASKS, FACES_TO_OCTANT_MASKS. -#[allow(dead_code)] -#[cfg(feature = "dim3")] -#[cfg(test)] -fn gen_const_tables() { - // The `j-th` bit of `faces_adj_to_vtx[i]` is set to 1, if the j-th face of the AABB (based on - // the face order depicted in `AABB::FACES_VERTEX_IDS`) is adjacent to the `i` vertex of the AABB - // (vertices are indexed as per the diagram depicted in the `FACES_VERTEX_IDS` doc. - // Each entry of this will always have exactly 3 bits set. - let mut faces_adj_to_vtx = [0usize; 8]; - - // The `j-th` bit of `faces_adj_to_vtx[i]` is set to 1, if the j-th edge of the AABB (based on - // the edge order depicted in `AABB::EDGES_VERTEX_IDS`) is adjacent to the `i` vertex of the AABB - // (vertices are indexed as per the diagram depicted in the `FACES_VERTEX_IDS` doc. - // Each entry of this will always have exactly 2 bits set. - let mut faces_adj_to_edge = [0usize; 12]; - - for fid in 0..6 { - let vids = Aabb::FACES_VERTEX_IDS[fid]; - let key = 1 << fid; - faces_adj_to_vtx[vids.0] |= key; - faces_adj_to_vtx[vids.1] |= key; - faces_adj_to_vtx[vids.2] |= key; - faces_adj_to_vtx[vids.3] |= key; - } - - #[allow(clippy::needless_range_loop)] - for eid in 0..12 { - let evids = Aabb::EDGES_VERTEX_IDS[eid]; - for fid in 0..6 { - let fvids = Aabb::FACES_VERTEX_IDS[fid]; - if (fvids.0 == evids.0 - || fvids.1 == evids.0 - || fvids.2 == evids.0 - || fvids.3 == evids.0) - && (fvids.0 == evids.1 - || fvids.1 == evids.1 - || fvids.2 == evids.1 - || fvids.3 == evids.1) - { - let key = 1 << fid; - faces_adj_to_edge[eid] |= key; - } - } - } - - /* - * FACES_TO_VOXEL_TYPES - */ - std::println!("const FACES_TO_VOXEL_TYPES: [VoxelType; 65] = ["); - 'outer: for i in 0usize..64 { - // If any vertex of the voxel has three faces with no adjacent voxels, - // then the voxel type is Vertex. - for adjs in faces_adj_to_vtx.iter() { - if (*adjs & i) == 0 { - std::println!("VoxelType::Vertex,"); - continue 'outer; - } - } - - // If any vertex of the voxel has three faces with no adjacent voxels, - // then the voxel type is Edge. - for adjs in faces_adj_to_edge.iter() { - if (*adjs & i) == 0 { - std::println!("VoxelType::Edge,"); - continue 'outer; - } - } - - // If one face doesn’t have any adjacent voxel, - // then the voxel type is Face. - for fid in 0..6 { - if ((1 << fid) & i) == 0 { - std::println!("VoxelType::Face,"); - continue 'outer; - } - } - } - - // Add final entries for special values. - std::println!("VoxelType::Interior,"); - std::println!("VoxelType::Empty,"); - std::println!("];"); - - /* - * FACES_TO_FEATURE_MASKS - */ - std::println!("const FACES_TO_FEATURE_MASKS: [u16; 65] = ["); - for i in 0usize..64 { - // Each bit set indicates a convex vertex that can lead to collisions. - // The result will be nonzero only for `VoxelType::Vertex` voxels. - let mut vtx_key = 0; - for (vid, adjs) in faces_adj_to_vtx.iter().enumerate() { - if (*adjs & i) == 0 { - vtx_key |= 1 << vid; - } - } - - if vtx_key != 0 { - std::println!("0b{:b},", vtx_key as u16); - continue; - } - - // Each bit set indicates a convex edge that can lead to collisions. - // The result will be nonzero only for `VoxelType::Edge` voxels. - let mut edge_key = 0; - for (eid, adjs) in faces_adj_to_edge.iter().enumerate() { - if (*adjs & i) == 0 { - edge_key |= 1 << eid; - } - } - - if edge_key != 0 { - std::println!("0b{:b},", edge_key as u16); - continue; - } - - // Each bit set indicates an exposed face that can lead to collisions. - // The result will be nonzero only for `VoxelType::Face` voxels. - let mut face_key = 0; - for fid in 0..6 { - if ((1 << fid) & i) == 0 { - face_key |= 1 << fid; - } - } - - if face_key != 0 { - std::println!("0b{:b},", face_key as u16); - continue; - } - } - - std::println!("0b{:b},", u16::MAX); - std::println!("0,"); - std::println!("];"); - - /* - * Faces to octant masks. - */ - std::println!("const FACES_TO_OCTANT_MASKS: [u32; 65] = ["); - for i in 0usize..64 { - // First test if we have vertices. - let mut octant_mask = 0; - let mut set_mask = |mask, octant| { - // NOTE: we don’t overwrite any mask already set for the octant. - if (octant_mask >> (octant * 3)) & 0b0111 == 0 { - octant_mask |= mask << (octant * 3); - } - }; - - for (vid, adjs) in faces_adj_to_vtx.iter().enumerate() { - if (*adjs & i) == 0 { - set_mask(1, vid); - } - } - - // This is the index of the axis porting the edges given by - // Aabb::EDGES_VERTEX_IDS. - const EX: u32 = OctantPattern::EDGE_X; - const EY: u32 = OctantPattern::EDGE_Y; - const EZ: u32 = OctantPattern::EDGE_Z; - const EDGE_AXIS: [u32; 12] = [EX, EY, EX, EY, EX, EY, EX, EY, EZ, EZ, EZ, EZ]; - for (eid, adjs) in faces_adj_to_edge.iter().enumerate() { - if (*adjs & i) == 0 { - let vid = Aabb::EDGES_VERTEX_IDS[eid]; - let mask = EDGE_AXIS[eid]; - - set_mask(mask, vid.0); - set_mask(mask, vid.1); - } - } - - // This is the index of the normal of the faces given by - // Aabb::FACES_VERTEX_IDS. - const FX: u32 = OctantPattern::FACE_X; - const FY: u32 = OctantPattern::FACE_Y; - const FZ: u32 = OctantPattern::FACE_Z; - const FACE_NORMALS: [u32; 6] = [FX, FX, FY, FY, FZ, FZ]; - - #[allow(clippy::needless_range_loop)] - for fid in 0..6 { - if ((1 << fid) & i) == 0 { - let vid = Aabb::FACES_VERTEX_IDS[fid]; - let mask = FACE_NORMALS[fid]; - - set_mask(mask, vid.0); - set_mask(mask, vid.1); - set_mask(mask, vid.2); - set_mask(mask, vid.3); - } - } - std::println!("0b{:b},", octant_mask); - } - std::println!("0,"); - std::println!("];"); -} - -// Index to the item of FACES_TO_VOXEL_TYPES which identifies interior voxels. -#[cfg(feature = "dim2")] -const INTERIOR_FACE_MASK: u8 = 0b0000_1111; -#[cfg(feature = "dim3")] -const INTERIOR_FACE_MASK: u8 = 0b0011_1111; -// Index to the item of FACES_TO_VOXEL_TYPES which identifies empty voxels. - -#[cfg(feature = "dim2")] -const EMPTY_FACE_MASK: u8 = 0b0001_0000; -#[cfg(feature = "dim3")] -const EMPTY_FACE_MASK: u8 = 0b0100_0000; - -/// The voxel type deduced from adjacency information. -/// -/// See the documentation of [`VoxelType`] for additional information on what each enum variant -/// means. -/// -/// In 3D there are 6 neighbor faces => 64 cases + 1 empty case. -#[cfg(feature = "dim3")] -const FACES_TO_VOXEL_TYPES: [VoxelType; 65] = [ - VoxelType::Vertex, - VoxelType::Vertex, - VoxelType::Vertex, - VoxelType::Edge, - VoxelType::Vertex, - VoxelType::Vertex, - VoxelType::Vertex, - VoxelType::Edge, - VoxelType::Vertex, - VoxelType::Vertex, - VoxelType::Vertex, - VoxelType::Edge, - VoxelType::Edge, - VoxelType::Edge, - VoxelType::Edge, - VoxelType::Face, - VoxelType::Vertex, - VoxelType::Vertex, - VoxelType::Vertex, - VoxelType::Edge, - VoxelType::Vertex, - VoxelType::Vertex, - VoxelType::Vertex, - VoxelType::Edge, - VoxelType::Vertex, - VoxelType::Vertex, - VoxelType::Vertex, - VoxelType::Edge, - VoxelType::Edge, - VoxelType::Edge, - VoxelType::Edge, - VoxelType::Face, - VoxelType::Vertex, - VoxelType::Vertex, - VoxelType::Vertex, - VoxelType::Edge, - VoxelType::Vertex, - VoxelType::Vertex, - VoxelType::Vertex, - VoxelType::Edge, - VoxelType::Vertex, - VoxelType::Vertex, - VoxelType::Vertex, - VoxelType::Edge, - VoxelType::Edge, - VoxelType::Edge, - VoxelType::Edge, - VoxelType::Face, - VoxelType::Edge, - VoxelType::Edge, - VoxelType::Edge, - VoxelType::Face, - VoxelType::Edge, - VoxelType::Edge, - VoxelType::Edge, - VoxelType::Face, - VoxelType::Edge, - VoxelType::Edge, - VoxelType::Edge, - VoxelType::Face, - VoxelType::Face, - VoxelType::Face, - VoxelType::Face, - VoxelType::Interior, - VoxelType::Empty, -]; - -/// Indicates the convex features of each voxel that can lead to collisions. -/// -/// The interpretation of each bit differs depending on the corresponding voxel type in -/// `FACES_TO_VOXEL_TYPES`: -/// - For `VoxelType::Vertex`: the i-th bit set to `1` indicates that the i-th AABB vertex is convex -/// and might lead to collisions. -/// - For `VoxelType::Edge`: the i-th bit set to `1` indicates that the i-th edge from `Aabb::EDGES_VERTEX_IDS` -/// is convex and might lead to collisions. -/// - For `VoxelType::Face`: the i-th bit set to `1` indicates that the i-th face from `Aabb::FACES_VERTEX_IDS` -/// is exposed and might lead to collisions. -#[cfg(feature = "dim3")] -const FACES_TO_FEATURE_MASKS: [u16; 65] = [ - 0b11111111, - 0b10011001, - 0b1100110, - 0b1010101, - 0b110011, - 0b10001, - 0b100010, - 0b10001, - 0b11001100, - 0b10001000, - 0b1000100, - 0b1000100, - 0b10101010, - 0b10001000, - 0b100010, - 0b110000, - 0b1111, - 0b1001, - 0b110, - 0b101, - 0b11, - 0b1, - 0b10, - 0b1, - 0b1100, - 0b1000, - 0b100, - 0b100, - 0b1010, - 0b1000, - 0b10, - 0b100000, - 0b11110000, - 0b10010000, - 0b1100000, - 0b1010000, - 0b110000, - 0b10000, - 0b100000, - 0b10000, - 0b11000000, - 0b10000000, - 0b1000000, - 0b1000000, - 0b10100000, - 0b10000000, - 0b100000, - 0b10000, - 0b111100000000, - 0b100100000000, - 0b11000000000, - 0b1100, - 0b1100000000, - 0b100000000, - 0b1000000000, - 0b1000, - 0b110000000000, - 0b100000000000, - 0b10000000000, - 0b100, - 0b11, - 0b10, - 0b1, - 0b1111111111111111, - 0, -]; - -/// Each octant is assigned three contiguous bits. -#[cfg(feature = "dim3")] -const FACES_TO_OCTANT_MASKS: [u32; 65] = [ - 0b1001001001001001001001, - 0b1010010001001010010001, - 0b10001001010010001001010, - 0b10010010010010010010010, - 0b11011001001011011001001, - 0b11111010001011111010001, - 0b111011001010111011001010, - 0b111111010010111111010010, - 0b1001011011001001011011, - 0b1010111011001010111011, - 0b10001011111010001011111, - 0b10010111111010010111111, - 0b11011011011011011011011, - 0b11111111011011111111011, - 0b111011011111111011011111, - 0b111111111111111111111111, - 0b100100100100001001001001, - 0b100110110100001010010001, - 0b110100100110010001001010, - 0b110110110110010010010010, - 0b101101100100011011001001, - 0b101000110100011111010001, - 0b101100110111011001010, - 0b110110111111010010, - 0b100100101101001001011011, - 0b100110000101001010111011, - 0b110100101000010001011111, - 0b110110000000010010111111, - 0b101101101101011011011011, - 0b101000000101011111111011, - 0b101101000111011011111, - 0b111111111111, - 0b1001001001100100100100, - 0b1010010001100110110100, - 0b10001001010110100100110, - 0b10010010010110110110110, - 0b11011001001101101100100, - 0b11111010001101000110100, - 0b111011001010000101100110, - 0b111111010010000000110110, - 0b1001011011100100101101, - 0b1010111011100110000101, - 0b10001011111110100101000, - 0b10010111111110110000000, - 0b11011011011101101101101, - 0b11111111011101000000101, - 0b111011011111000101101000, - 0b111111111111000000000000, - 0b100100100100100100100100, - 0b100110110100100110110100, - 0b110100100110110100100110, - 0b110110110110110110110110, - 0b101101100100101101100100, - 0b101000110100101000110100, - 0b101100110000101100110, - 0b110110000000110110, - 0b100100101101100100101101, - 0b100110000101100110000101, - 0b110100101000110100101000, - 0b110110000000110110000000, - 0b101101101101101101101101, - 0b101000000101101000000101, - 0b101101000000101101000, - 0b0, - 0, -]; - -#[cfg(feature = "dim2")] -const FACES_TO_VOXEL_TYPES: [VoxelType; 17] = [ - VoxelType::Vertex, - VoxelType::Vertex, - VoxelType::Vertex, - VoxelType::Face, - VoxelType::Vertex, - VoxelType::Vertex, - VoxelType::Vertex, - VoxelType::Face, - VoxelType::Vertex, - VoxelType::Vertex, - VoxelType::Vertex, - VoxelType::Face, - VoxelType::Face, - VoxelType::Face, - VoxelType::Face, - VoxelType::Interior, - VoxelType::Empty, -]; - -#[cfg(feature = "dim2")] -const FACES_TO_FEATURE_MASKS: [u16; 17] = [ - 0b1111, - 0b1001, - 0b110, - 0b1100, - 0b11, - 0b1, - 0b10, - 0b1000, - 0b1100, - 0b1000, - 0b100, - 0b100, - 0b11, - 0b10, - 0b1, - 0b1111111111111111, - 0, -]; - -// NOTE: in 2D we are also using 3 bits per octant even though we technically only need two. -// This keeps some collision-detection easier by avoiding some special-casing. -#[cfg(feature = "dim2")] -const FACES_TO_OCTANT_MASKS: [u32; 17] = [ - 0b1001001001, - 0b1011011001, - 0b11001001011, - 0b11011011011, - 0b10010001001, - 0b10000011001, - 0b10001011, - 0b11011, - 0b1001010010, - 0b1011000010, - 0b11001010000, - 0b11011000000, - 0b10010010010, - 0b10000000010, - 0b10010000, - 0b0, - 0, -]; - -#[cfg(test)] -mod test { - #[test] - fn gen_const_tables() { - super::gen_const_tables(); - } -} diff --git a/src/shape/voxels/mod.rs b/src/shape/voxels/mod.rs new file mode 100644 index 00000000..ea8b7c0a --- /dev/null +++ b/src/shape/voxels/mod.rs @@ -0,0 +1,11 @@ +pub use voxels::*; +pub use voxels_chunk::*; + +use voxels_consts::*; + +mod voxels; +mod voxels_chunk; +mod voxels_consts; +mod voxels_edition; +mod voxels_neighborhood; +mod voxels_sizeof; diff --git a/src/shape/voxels/voxels.rs b/src/shape/voxels/voxels.rs new file mode 100644 index 00000000..68718757 --- /dev/null +++ b/src/shape/voxels/voxels.rs @@ -0,0 +1,535 @@ +use super::{ + VoxelIndex, EMPTY_FACE_MASK, FACES_TO_FEATURE_MASKS, FACES_TO_OCTANT_MASKS, + FACES_TO_VOXEL_TYPES, INTERIOR_FACE_MASK, +}; +use crate::bounding_volume::{Aabb, BoundingVolume}; +use crate::math::{Point, Real, Vector}; +use crate::partitioning::{Bvh, BvhBuildStrategy, BvhNode}; +use crate::shape::voxels::voxels_chunk::{VoxelsChunk, VoxelsChunkHeader}; +use crate::shape::VoxelsChunkRef; +use crate::utils::hashmap::HashMap; +use alloc::{vec, vec::Vec}; +#[cfg(not(feature = "std"))] +use na::ComplexField; + +/// Categorization of a voxel based on its neighbors. +#[derive(Copy, Clone, Debug, PartialEq, Eq)] +pub enum VoxelType { + /// The voxel is empty. + Empty, + /// The voxel is a vertex if all three coordinate axis directions have at + /// least one empty neighbor. + Vertex, + /// The voxel is on an edge if it has non-empty neighbors in both directions of + /// a single coordinate axis. + #[cfg(feature = "dim3")] + Edge, + /// The voxel is on an edge if it has non-empty neighbors in both directions of + /// two coordinate axes. + Face, + /// The voxel is on an edge if it has non-empty neighbors in both directions of + /// all coordinate axes. + Interior, +} + +#[derive(Clone, Copy, Debug, Default, Eq, Hash, Ord, PartialEq, PartialOrd)] +/// The status of the cell of an heightfield. +pub struct AxisMask(u8); + +bitflags::bitflags! { + /// Flags for identifying signed directions along coordinate axes, or faces of a voxel. + impl AxisMask: u8 { + /// The direction or face along the `+x` coordinate axis. + const X_POS = 1 << 0; + /// The direction or face along the `-x` coordinate axis. + const X_NEG = 1 << 1; + /// The direction or face along the `+y` coordinate axis. + const Y_POS = 1 << 2; + /// The direction or face along the `-y` coordinate axis. + const Y_NEG = 1 << 3; + /// The direction or face along the `+z` coordinate axis. + #[cfg(feature= "dim3")] + const Z_POS = 1 << 4; + /// The direction or face along the `-z` coordinate axis. + #[cfg(feature= "dim3")] + const Z_NEG = 1 << 5; + } +} + +/// Indicates the local shape of a voxel on each octant. +/// +/// This provides geometric information of the shape’s exposed features on each octant. +// This is an alternative to `FACES_TO_FEATURE_MASKS` that can be more convenient for some +// collision-detection algorithms. +#[derive(Copy, Clone, Debug, PartialEq, Eq)] +pub struct OctantPattern; + +// NOTE: it is important that the max value of any OctantPattern variant +// is 7 because we don’t allocate more than 3 bits to store it in +// `FACES_TO_OCTANT_MASKS`. +/// Indicates the local shape of a voxel on each octant. +/// +/// This provides geometric information of the shape’s exposed features on each octant. +// This is an alternative to `FACES_TO_FEATURE_MASKS` that can be more convenient for some +// collision-detection algorithms. +#[cfg(feature = "dim3")] +impl OctantPattern { + /// The voxel doesn't have any exposed feature on the octant with this mask. + pub const INTERIOR: u32 = 0; + /// The voxel has an exposed vertex on the octant with this mask. + pub const VERTEX: u32 = 1; + /// The voxel has an exposed edges with direction X on the octant with this mask. + pub const EDGE_X: u32 = 2; + /// The voxel has an exposed edges with direction Y on the octant with this mask. + pub const EDGE_Y: u32 = 3; + /// The voxel has an exposed edges with direction Z on the octant with this mask. + pub const EDGE_Z: u32 = 4; + /// The voxel has an exposed face with normal X on the octant with this mask. + pub const FACE_X: u32 = 5; + /// The voxel has an exposed face with normal Y on the octant with this mask. + pub const FACE_Y: u32 = 6; + /// The voxel has an exposed face with normal Z on the octant with this mask. + pub const FACE_Z: u32 = 7; +} + +// NOTE: it is important that the max value of any OctantPattern variant +// is 7 because we don’t allocate more than 3 bits to store it in +// `FACES_TO_OCTANT_MASKS`. +/// Indicates the local shape of a voxel on each octant. +/// +/// This provides geometric information of the shape’s exposed features on each octant. +/// This is an alternative to `FACES_TO_FEATURE_MASKS` that can be more convenient for some +/// collision-detection algorithms. +#[cfg(feature = "dim2")] +impl OctantPattern { + /// The voxel doesn't have any exposed feature on the octant with this mask. + pub const INTERIOR: u32 = 0; + /// The voxel has an exposed vertex on the octant with this mask. + pub const VERTEX: u32 = 1; + /// The voxel has an exposed face with normal X on the octant with this mask. + pub const FACE_X: u32 = 2; + /// The voxel has an exposed face with normal Y on the octant with this mask. + pub const FACE_Y: u32 = 3; +} + +// The local neighborhood information is encoded in a 8-bits number in groups of two bits +// per coordinate axis: `0bwwzzyyxx`. In each group of two bits, e.g. `xx`, the rightmost (resp. +// leftmost) bit set to 1 means that the neighbor voxel with coordinate `+1` (resp `-1`) relative +// to the current voxel along the `x` axis is filled. If the bit is 0, then the corresponding +// neighbor is empty. See the `AxisMask` bitflags. +// For example, in 2D, the mask `0b00_00_10_01` matches the following configuration (assuming +y goes +// up, and +x goes right): +// +// ```txt +// 0 0 0 +// 0 x 1 +// 0 1 0 +// ``` +// +// The special value `0b01000000` indicates that the voxel is empty. +// And the value `0b00111111` (`0b00001111` in 2D) indicates that the voxel is an interior voxel (its whole neighborhood +// is filled). +/// A description of the local neighborhood of a voxel. +#[derive(Copy, Clone, Debug, PartialEq, Eq)] +#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))] +pub struct VoxelState(pub(super) u8); + +impl VoxelState { + /// The value of empty voxels. + pub const EMPTY: VoxelState = VoxelState(EMPTY_FACE_MASK); + /// The value of a voxel with non-empty neighbors in all directions. + pub const INTERIOR: VoxelState = VoxelState(INTERIOR_FACE_MASK); + + pub(crate) const fn new(state: u8) -> Self { + Self(state) + } + + /// Is this voxel empty? + pub const fn is_empty(self) -> bool { + self.0 == EMPTY_FACE_MASK + } + + /// A bit mask indicating which faces of the voxel don’t have any + /// adjacent non-empty voxel. + pub const fn free_faces(self) -> AxisMask { + if self.0 == INTERIOR_FACE_MASK || self.0 == EMPTY_FACE_MASK { + AxisMask::empty() + } else { + AxisMask::from_bits_truncate((!self.0) & INTERIOR_FACE_MASK) + } + } + + /// The [`VoxelType`] of this voxel. + pub const fn voxel_type(self) -> VoxelType { + FACES_TO_VOXEL_TYPES[self.0 as usize] + } + + // Bitmask indicating what vertices, edges, or faces of the voxel are "free". + pub(crate) const fn feature_mask(self) -> u16 { + FACES_TO_FEATURE_MASKS[self.0 as usize] + } + + pub(crate) const fn octant_mask(self) -> u32 { + FACES_TO_OCTANT_MASKS[self.0 as usize] + } +} + +/// Information associated to a voxel. +#[derive(Copy, Clone, Debug, PartialEq)] +pub struct VoxelData { + /// The temporary index in the internal voxels’ storage. + /// + /// This index can be invalidated after a call to [`Voxels::set_voxel`], or + /// [`Voxels::crop`]. + pub linear_id: VoxelIndex, + /// The voxel’s integer grid coordinates. + pub grid_coords: Point, + /// The voxel’s center position in the local-space of the [`Voxels`] shape it is part of. + pub center: Point, + /// The voxel’s state, indicating if it’s empty or full. + pub state: VoxelState, +} + +/// A shape made of axis-aligned, uniformly sized, cubes (aka. voxels). +/// +/// This shape is specialized to handle voxel worlds and voxelized obojects efficiently why ensuring +/// that collision-detection isn’t affected by the so-called "internal edges problem" that can create +/// artifacts when another object rolls or slides against a flat voxelized surface. +/// +/// The internal storage is compact (but not sparse at the moment), storing only one byte per voxel +/// in the allowed domain. This has a generally smaller memory footprint than a mesh representation +/// of the voxels. +#[derive(Clone, Debug)] +#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))] +pub struct Voxels { + /// A BVH of chunk keys. + /// + /// The bounding boxes are the ones of the chunk’s voxels **keys**. This is equivalent to a bvh + /// of the chunks with a uniform voxel size of 1. + pub(super) chunk_bvh: Bvh, + pub(super) chunk_headers: HashMap, VoxelsChunkHeader>, + pub(super) chunk_keys: Vec>, + pub(super) chunks: Vec, + pub(super) free_chunks: Vec, + pub(super) voxel_size: Vector, +} + +impl Voxels { + /// Initializes a voxel shapes from the voxels grid coordinates. + /// + /// Each voxel will have its bottom-left-back corner located at + /// `grid_coordinates * voxel_size`; and its center at `(grid_coordinates + 0.5) * voxel_size`. + pub fn new(voxel_size: Vector, grid_coordinates: &[Point]) -> Self { + let mut result = Self { + chunk_bvh: Bvh::new(), + chunk_headers: HashMap::default(), + chunk_keys: vec![], + chunks: vec![], + free_chunks: vec![], + voxel_size, + }; + + for vox in grid_coordinates { + let (chunk_key, id_in_chunk) = Self::chunk_key_and_id_in_chunk(*vox); + let chunk_header = result.chunk_headers.entry(chunk_key).or_insert_with(|| { + let id = result.chunks.len(); + result.chunks.push(VoxelsChunk::default()); + result.chunk_keys.push(chunk_key); + VoxelsChunkHeader { id, len: 0 } + }); + chunk_header.len += 1; + result.chunks[chunk_header.id].states[id_in_chunk] = VoxelState::INTERIOR; + } + + result.chunk_bvh = Bvh::from_iter( + BvhBuildStrategy::Ploc, + result.chunk_headers.iter().map(|(chunk_key, chunk_id)| { + ( + chunk_id.id, + VoxelsChunk::aabb(chunk_key, &result.voxel_size), + ) + }), + ); + + result.recompute_all_voxels_states(); + result + } + + /// Computes a voxels shape from the set of `points`. + /// + /// The points are mapped to a regular grid centered at the provided point with smallest + /// coordinates, and with grid cell size equal to `scale`. It is OK if multiple points + /// fall into the same grid cell. + pub fn from_points(voxel_size: Vector, points: &[Point]) -> Self { + let voxels: Vec<_> = points + .iter() + .map(|pt| { + Point::from( + pt.coords + .component_div(&voxel_size) + .map(|x| x.floor() as i32), + ) + }) + .collect(); + Self::new(voxel_size, &voxels) + } + + pub(crate) fn chunk_bvh(&self) -> &Bvh { + &self.chunk_bvh + } + + // /// The extents of the total axis-aligned volume covered by this [`Voxels`] shape. + // /// + // /// This accounts for all the voxels reserved in the internal buffer of `self`, including empty + // /// ones. + // pub fn extents(&self) -> Vector { + // self.dimensions() + // .cast::() + // .component_mul(&self.voxel_size) + // } + + // /// The center of this shape’s domain (accounting for both empty and filled voxels). + // pub fn domain_center(&self) -> Point { + // (self + // .domain_mins + // .coords + // .cast::() + // .component_mul(&self.voxel_size) + // + self.extents() / 2.0) + // .into() + // } + + /// The semi-open range of voxels in shape. + /// + /// This provides conservative bounds on the range of voxel indices that might be set to filled. + pub fn domain(&self) -> [Point; 2] { + let aabb = self.chunk_bvh.root_aabb(); + + // NOTE that we shift the AABB’s bounds so the endpoint matches a voxel center + // to avoid rounding errors. + let half_sz = self.voxel_size() / 2.0; + let mins = self.voxel_at_point(aabb.mins + half_sz); + // + 1 because the range is semi-open. + let maxs = self.voxel_at_point(aabb.maxs - half_sz) + Vector::repeat(1); + [mins, maxs] + } + + // /// The number of voxels along each coordinate axis. + // pub fn dimensions(&self) -> Vector { + // (self.domain_maxs - self.domain_mins).map(|e| e as u32) + // } + + /// The size of each voxel part this [`Voxels`] shape. + pub fn voxel_size(&self) -> Vector { + self.voxel_size + } + + /// Scale this shape. + pub fn scaled(mut self, scale: &Vector) -> Self { + self.voxel_size.component_mul_assign(scale); + self + } + + /// A reference to the chunk with id `chunk_id`. + /// + /// Panics if the chunk doesn’t exist. + pub fn chunk_ref(&self, chunk_id: u32) -> VoxelsChunkRef<'_> { + VoxelsChunkRef { + my_id: chunk_id as usize, + parent: self, + states: &self.chunks[chunk_id as usize].states, + key: &self.chunk_keys[chunk_id as usize], + } + } + + /// The AABB of the voxel with the given quantized `key`. + pub fn voxel_aabb(&self, key: Point) -> Aabb { + let center = self.voxel_center(key); + let hext = self.voxel_size / 2.0; + Aabb::from_half_extents(center, hext) + } + + /// Returns the state of a given voxel. + pub fn voxel_state(&self, key: Point) -> Option { + let vid = self.linear_index(key)?; + Some(self.chunks[vid.chunk_id].states[vid.id_in_chunk]) + } + + /// Calculates the grid coordinates of the voxel containing the given `point`, regardless + /// of whether this voxel is filled oor empty. + pub fn voxel_at_point(&self, point: Point) -> Point { + point + .coords + .component_div(&self.voxel_size) + .map(|x| x.floor() as i32) + .into() + } + + /// Gets the voxel at the given flat voxel index. + pub fn voxel_at_flat_id(&self, id: u32) -> Option> { + let vid = VoxelIndex::from_flat_id(id as usize); + let chunk_key = self.chunk_keys.get(vid.chunk_id)?; + if *chunk_key == VoxelsChunk::INVALID_CHUNK_KEY { + return None; + } + + Some(VoxelsChunk::voxel_key_at_id( + *chunk_key, + vid.id_in_chunk as u32, + )) + } + + /// The range of grid coordinates of voxels intersecting the given AABB. + /// + /// The returned range covers both empty and non-empty voxels, and is not limited to the + /// bounds defined by [`Self::domain`]. + /// The range is semi, open, i.e., the range along each dimension `i` is understood as + /// the semi-open interval: `range[0][i]..range[1][i]`. + pub fn voxel_range_intersecting_local_aabb(&self, aabb: &Aabb) -> [Point; 2] { + let mins = aabb + .mins + .coords + .component_div(&self.voxel_size) + .map(|x| x.floor() as i32); + let maxs = aabb + .maxs + .coords + .component_div(&self.voxel_size) + .map(|x| x.ceil() as i32); + [mins.into(), maxs.into()] + } + + /// The AABB of a given range of voxels. + /// + /// The AABB is computed independently of [`Self::domain`] and independently of whether + /// the voxels contained within are empty or not. + pub fn voxel_range_aabb(&self, mins: Point, maxs: Point) -> Aabb { + Aabb { + mins: mins + .cast::() + .coords + .component_mul(&self.voxel_size) + .into(), + maxs: maxs + .cast::() + .coords + .component_mul(&self.voxel_size) + .into(), + } + } + + /// Aligns the given AABB with the voxelized grid. + /// + /// The aligned is calculated such that the returned AABB has corners lying at the grid + /// intersections (i.e. matches voxel corners) and fully contains the input `aabb`. + pub fn align_aabb_to_grid(&self, aabb: &Aabb) -> Aabb { + let mins = aabb + .mins + .coords + .zip_map(&self.voxel_size, |m, sz| (m / sz).floor() * m) + .into(); + let maxs = aabb + .maxs + .coords + .zip_map(&self.voxel_size, |m, sz| (m / sz).ceil() * m) + .into(); + Aabb { mins, maxs } + } + + /// Iterates through every voxel intersecting the given aabb. + /// + /// Returns the voxel’s linearized id, center, and state. + pub fn voxels_intersecting_local_aabb( + &self, + aabb: &Aabb, + ) -> impl Iterator + '_ { + let [mins, maxs] = self.voxel_range_intersecting_local_aabb(aabb); + self.voxels_in_range(mins, maxs) + } + + /// The center point of all the voxels in this shape (including empty ones). + /// + /// The voxel data associated to each center is provided to determine what kind of voxel + /// it is (and, in particular, if it is empty or full). + pub fn voxels(&self) -> impl Iterator + '_ { + let aabb = self.chunk_bvh.root_aabb(); + self.voxels_in_range( + self.voxel_at_point(aabb.mins), + self.voxel_at_point(aabb.maxs), + ) + } + + /// Iterate through the data of all the voxels within the given (semi-open) voxel grid indices. + /// + /// Note that this yields both empty and non-empty voxels within the range. This does not + /// include any voxel that falls outside [`Self::domain`]. + pub fn voxels_in_range( + &self, + mins: Point, + maxs: Point, + ) -> impl Iterator + '_ { + let range_aabb = Aabb::new(self.voxel_center(mins), self.voxel_center(maxs)); + + self.chunk_bvh + .leaves(move |node: &BvhNode| node.aabb().intersects(&range_aabb)) + .flat_map(move |chunk_id| { + let chunk = self.chunk_ref(chunk_id); + chunk.voxels_in_range(mins, maxs) + }) + } + + fn voxel_to_chunk_key(voxel_key: Point) -> Point { + fn div_floor(a: i32, b: usize) -> i32 { + let sign = (a < 0) as i32; + (a + sign) / b as i32 - sign + } + + voxel_key.map(|e| div_floor(e, VoxelsChunk::VOXELS_PER_CHUNK_DIM)) + } + + /// Given a voxel key, returns the key of the voxel chunk that contains it, as well as the + /// linear index of the voxel within that chunk. + #[cfg(feature = "dim2")] + pub(super) fn chunk_key_and_id_in_chunk(voxel_key: Point) -> (Point, usize) { + let chunk_key = Self::voxel_to_chunk_key(voxel_key); + // NOTE: always positive since we subtracted the smallest possible key on that chunk. + let voxel_key_in_chunk = voxel_key - chunk_key * VoxelsChunk::VOXELS_PER_CHUNK_DIM as i32; + let id_in_chunk = (voxel_key_in_chunk.x + + voxel_key_in_chunk.y * VoxelsChunk::VOXELS_PER_CHUNK_DIM as i32) + as usize; + (chunk_key, id_in_chunk) + } + + /// Given a voxel key, returns the key of the voxel chunk that contains it, as well as the + /// linear index of the voxel within that chunk. + #[cfg(feature = "dim3")] + pub(super) fn chunk_key_and_id_in_chunk(voxel_key: Point) -> (Point, usize) { + let chunk_key = Self::voxel_to_chunk_key(voxel_key); + // NOTE: always positive since we subtracted the smallest possible key on that chunk. + let voxel_key_in_chunk = voxel_key - chunk_key * VoxelsChunk::VOXELS_PER_CHUNK_DIM as i32; + let id_in_chunk = (voxel_key_in_chunk.x + + voxel_key_in_chunk.y * VoxelsChunk::VOXELS_PER_CHUNK_DIM as i32 + + voxel_key_in_chunk.z + * VoxelsChunk::VOXELS_PER_CHUNK_DIM as i32 + * VoxelsChunk::VOXELS_PER_CHUNK_DIM as i32) as usize; + (chunk_key, id_in_chunk) + } + + /// The linearized index associated to the given voxel key. + pub fn linear_index(&self, voxel_key: Point) -> Option { + let (chunk_key, id_in_chunk) = Self::chunk_key_and_id_in_chunk(voxel_key); + let chunk_id = self.chunk_headers.get(&chunk_key)?.id; + Some(VoxelIndex { + chunk_id, + id_in_chunk, + }) + } + + /// The center of the voxel with the given key. + pub fn voxel_center(&self, key: Point) -> Point { + (key.cast::() + Vector::repeat(0.5)) + .coords + .component_mul(&self.voxel_size) + .into() + } +} diff --git a/src/shape/voxels/voxels_chunk.rs b/src/shape/voxels/voxels_chunk.rs new file mode 100644 index 00000000..604c903e --- /dev/null +++ b/src/shape/voxels/voxels_chunk.rs @@ -0,0 +1,274 @@ +use crate::bounding_volume::Aabb; +use crate::math::{Point, Real, Vector}; +use crate::shape::{VoxelData, VoxelState, Voxels}; +use na::point; + +#[derive(Clone, Debug)] +#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))] +pub(super) struct VoxelsChunkHeader { + pub(super) id: usize, + // The number of non-empty voxels in the chunk. + // This is used for detecting when a chunk can be removed + // if it becomes fully empty. + pub(super) len: usize, +} + +#[derive(Clone, Debug)] +#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))] +#[repr(C)] +#[repr(align(64))] +pub(super) struct VoxelsChunk { + #[cfg_attr(feature = "serde-serialize", serde(with = "serde_arrays"))] + pub(super) states: [VoxelState; VoxelsChunk::VOXELS_PER_CHUNK], +} + +impl Default for VoxelsChunk { + fn default() -> Self { + Self { + states: [VoxelState::EMPTY; VoxelsChunk::VOXELS_PER_CHUNK], + } + } +} + +#[derive(Copy, Clone, Debug, PartialEq)] +pub struct VoxelIndex { + pub(super) chunk_id: usize, + pub(super) id_in_chunk: usize, +} + +impl VoxelIndex { + pub fn flat_id(&self) -> usize { + self.chunk_id * VoxelsChunk::VOXELS_PER_CHUNK + self.id_in_chunk + } + + pub fn from_flat_id(id: usize) -> Self { + Self { + chunk_id: id / VoxelsChunk::VOXELS_PER_CHUNK, + id_in_chunk: id % VoxelsChunk::VOXELS_PER_CHUNK, + } + } +} + +impl VoxelsChunk { + // TODO: find the ideal number. We want a good balance between cache locality + // and number of BVH leaf nodes. + #[cfg(feature = "dim2")] + pub(super) const VOXELS_PER_CHUNK_DIM: usize = 16; + #[cfg(feature = "dim3")] + pub(super) const VOXELS_PER_CHUNK_DIM: usize = 8; + #[cfg(feature = "dim2")] + pub(super) const VOXELS_PER_CHUNK: usize = + Self::VOXELS_PER_CHUNK_DIM * Self::VOXELS_PER_CHUNK_DIM; + #[cfg(feature = "dim3")] + pub(super) const VOXELS_PER_CHUNK: usize = + Self::VOXELS_PER_CHUNK_DIM * Self::VOXELS_PER_CHUNK_DIM * Self::VOXELS_PER_CHUNK_DIM; + + #[cfg(feature = "dim2")] + pub(super) const INVALID_CHUNK_KEY: Point = point![i32::MAX, i32::MAX]; + #[cfg(feature = "dim3")] + pub(super) const INVALID_CHUNK_KEY: Point = point![i32::MAX, i32::MAX, i32::MAX]; + + /// The key of the voxel at the given linearized index within this chunk. + #[cfg(feature = "dim2")] + pub(super) fn voxel_key_at_id(chunk_key: Point, id_in_chunk: u32) -> Point { + let y = id_in_chunk as i32 / Self::VOXELS_PER_CHUNK_DIM as i32; + let x = id_in_chunk as i32 % Self::VOXELS_PER_CHUNK_DIM as i32; + chunk_key * (Self::VOXELS_PER_CHUNK_DIM as i32) + Vector::new(x, y) + } + + /// The key of the voxel at the given linearized index. + #[cfg(feature = "dim3")] + pub(super) fn voxel_key_at_id(chunk_key: Point, id_in_chunk: u32) -> Point { + let d0d1 = (Self::VOXELS_PER_CHUNK_DIM * Self::VOXELS_PER_CHUNK_DIM) as u32; + let z = id_in_chunk / d0d1; + let y = (id_in_chunk - z * d0d1) / Self::VOXELS_PER_CHUNK_DIM as u32; + let x = id_in_chunk % Self::VOXELS_PER_CHUNK_DIM as u32; + chunk_key * (Self::VOXELS_PER_CHUNK_DIM as i32) + Vector::new(x as i32, y as i32, z as i32) + } + + /// The semi-open range of valid voxel keys for this chunk. + pub(super) fn keys_bounds(chunk_key: &Point) -> [Point; 2] { + let imins = chunk_key * Self::VOXELS_PER_CHUNK_DIM as i32; + let imaxs = imins + Vector::repeat(Self::VOXELS_PER_CHUNK_DIM as i32); + [imins, imaxs] + } + + pub(super) fn aabb(chunk_key: &Point, voxel_size: &Vector) -> Aabb { + let [imins, imaxs] = Self::keys_bounds(chunk_key); + let mut aabb = Aabb::new(imins.cast(), imaxs.cast()); + aabb.mins.coords.component_mul_assign(voxel_size); + aabb.maxs.coords.component_mul_assign(voxel_size); + aabb + } +} + +/// A reference to a voxel chunk. +#[derive(Copy, Clone)] +pub struct VoxelsChunkRef<'a> { + /// The linear index of this chunk within the `Voxels` shape. + pub my_id: usize, + /// The voxel shape this chunk is part of. + pub parent: &'a Voxels, + /// The fill status of each voxel from this chunk. + pub states: &'a [VoxelState; VoxelsChunk::VOXELS_PER_CHUNK], + /// The spatial index of this chunk. + pub key: &'a Point, +} + +impl<'a> VoxelsChunkRef<'a> { + /// The AABB of this chunk of voxels. + /// + /// Note that this return the AABB of the whole chunk, without taking into account the fact + /// that some voxels are empty. + pub fn local_aabb(&self) -> Aabb { + VoxelsChunk::aabb(self.key, &self.parent.voxel_size) + } + + /// The domain of this chunk of voxels. + pub fn domain(&self) -> [Point; 2] { + VoxelsChunk::keys_bounds(self.key) + } + + /// Returns the spatial index of the voxel containing the given point. + pub fn voxel_at_point_unchecked(&self, pt: Point) -> Point { + self.parent.voxel_at_point(pt) + } + + /// The state of the voxel with key `voxel_key`. + pub fn voxel_state(&self, voxel_key: Point) -> Option { + let (chunk_key, id_in_chunk) = Voxels::chunk_key_and_id_in_chunk(voxel_key); + if &chunk_key != self.key { + return None; + } + Some(self.states[id_in_chunk]) + } + + /// Clamps the `voxel_key` so it is within the bounds of `self`. + pub fn clamp_voxel(&self, voxel_key: Point) -> Point { + let [mins, maxs] = self.domain(); + voxel_key + .coords + .zip_zip_map(&mins.coords, &maxs.coords, |k, min, max| k.clamp(min, max)) + .into() + } + + /// The AABB of the voxel with this key. + /// + /// Returns a result even if the voxel doesn’t belong to this chunk. + pub fn voxel_aabb_unchecked(&self, voxel_key: Point) -> Aabb { + self.parent.voxel_aabb(voxel_key) + } + + /// Convert a voxel index (expressed relative to the main `Voxels` shape, not relative to the + /// chunk alone) into a flat index within a voxel chunk. + /// + /// Returns `None` if the voxel isn’t part of this chunk. + pub fn flat_id(&self, voxel_key: Point) -> Option { + let (chunk_key, id_in_chunk) = Voxels::chunk_key_and_id_in_chunk(voxel_key); + if &chunk_key != self.key { + return None; + } + + Some( + VoxelIndex { + chunk_id: self.my_id, + id_in_chunk, + } + .flat_id() as u32, + ) + } + + /// Iterates through all the voxels in this chunk. + /// + /// Note that this yields both empty and non-empty voxels within the range. This does not + /// include any voxel that falls outside [`Self::domain`]. + pub fn voxels(&self) -> impl Iterator + '_ { + let range = self.domain(); + self.voxels_in_range(range[0], range[1]) + } + + /// Iterate through the data of all the voxels within the given (semi-open) voxel grid indices. + /// + /// Note that this yields both empty and non-empty voxels within the range. This does not + /// include any voxel that falls outside [`Self::domain`]. + #[cfg(feature = "dim2")] + pub fn voxels_in_range( + self, + mins: Point, + maxs: Point, + ) -> impl Iterator + use<'a> { + let [chunk_mins, chunk_maxs] = VoxelsChunk::keys_bounds(self.key); + let mins = mins.coords.sup(&chunk_mins.coords); + let maxs = maxs.coords.inf(&chunk_maxs.coords); + + (mins[0]..maxs[0]).flat_map(move |ix| { + (mins[1]..maxs[1]).flat_map(move |iy| { + let id_in_chunk = (ix - chunk_mins[0]) as usize + + (iy - chunk_mins[1]) as usize * VoxelsChunk::VOXELS_PER_CHUNK_DIM; + let state = self.states[id_in_chunk]; + + if state.is_empty() { + return None; + } + + let grid_coords = Point::new(ix, iy); + let center = Vector::new(ix as Real + 0.5, iy as Real + 0.5) + .component_mul(&self.parent.voxel_size); + Some(VoxelData { + linear_id: VoxelIndex { + chunk_id: self.my_id, + id_in_chunk, + }, + grid_coords, + center: center.into(), + state, + }) + }) + }) + } + + /// Iterate through the data of all the voxels within the given (semi-open) voxel grid indices. + /// + /// Note that this yields both empty and non-empty voxels within the range. This does not + /// include any voxel that falls outside [`Self::domain`]. + #[cfg(feature = "dim3")] + pub fn voxels_in_range( + self, + mins: Point, + maxs: Point, + ) -> impl Iterator + use<'a> { + let [chunk_mins, chunk_maxs] = VoxelsChunk::keys_bounds(self.key); + let mins = mins.coords.sup(&chunk_mins.coords); + let maxs = maxs.coords.inf(&chunk_maxs.coords); + + (mins[0]..maxs[0]).flat_map(move |ix| { + (mins[1]..maxs[1]).flat_map(move |iy| { + (mins[2]..maxs[2]).filter_map(move |iz| { + let id_in_chunk = (ix - chunk_mins[0]) as usize + + (iy - chunk_mins[1]) as usize * VoxelsChunk::VOXELS_PER_CHUNK_DIM + + (iz - chunk_mins[2]) as usize + * VoxelsChunk::VOXELS_PER_CHUNK_DIM + * VoxelsChunk::VOXELS_PER_CHUNK_DIM; + let state = self.states[id_in_chunk]; + + if state.is_empty() { + return None; + } + + let grid_coords = Point::new(ix, iy, iz); + let center = Vector::new(ix as Real + 0.5, iy as Real + 0.5, iz as Real + 0.5) + .component_mul(&self.parent.voxel_size); + Some(VoxelData { + linear_id: VoxelIndex { + chunk_id: self.my_id, + id_in_chunk, + }, + grid_coords, + center: center.into(), + state, + }) + }) + }) + }) + } +} diff --git a/src/shape/voxels/voxels_consts.rs b/src/shape/voxels/voxels_consts.rs new file mode 100644 index 00000000..9d8ad7dc --- /dev/null +++ b/src/shape/voxels/voxels_consts.rs @@ -0,0 +1,639 @@ +use super::VoxelType; + +#[cfg(test)] +use {crate::bounding_volume::Aabb, OctantPattern}; + +// Index to the item of FACES_TO_VOXEL_TYPES which identifies interior voxels. +#[cfg(feature = "dim2")] +pub(super) const INTERIOR_FACE_MASK: u8 = 0b0000_1111; +#[cfg(feature = "dim3")] +pub(super) const INTERIOR_FACE_MASK: u8 = 0b0011_1111; +// Index to the item of FACES_TO_VOXEL_TYPES which identifies empty voxels. + +#[cfg(feature = "dim2")] +pub(super) const EMPTY_FACE_MASK: u8 = 0b0001_0000; +#[cfg(feature = "dim3")] +pub(super) const EMPTY_FACE_MASK: u8 = 0b0100_0000; + +/// The voxel type deduced from adjacency information. +/// +/// See the documentation of [`VoxelType`] for additional information on what each enum variant +/// means. +/// +/// In 3D there are 6 neighbor faces => 64 cases + 1 empty case. +#[cfg(feature = "dim3")] +pub(super) const FACES_TO_VOXEL_TYPES: [VoxelType; 65] = [ + VoxelType::Vertex, + VoxelType::Vertex, + VoxelType::Vertex, + VoxelType::Edge, + VoxelType::Vertex, + VoxelType::Vertex, + VoxelType::Vertex, + VoxelType::Edge, + VoxelType::Vertex, + VoxelType::Vertex, + VoxelType::Vertex, + VoxelType::Edge, + VoxelType::Edge, + VoxelType::Edge, + VoxelType::Edge, + VoxelType::Face, + VoxelType::Vertex, + VoxelType::Vertex, + VoxelType::Vertex, + VoxelType::Edge, + VoxelType::Vertex, + VoxelType::Vertex, + VoxelType::Vertex, + VoxelType::Edge, + VoxelType::Vertex, + VoxelType::Vertex, + VoxelType::Vertex, + VoxelType::Edge, + VoxelType::Edge, + VoxelType::Edge, + VoxelType::Edge, + VoxelType::Face, + VoxelType::Vertex, + VoxelType::Vertex, + VoxelType::Vertex, + VoxelType::Edge, + VoxelType::Vertex, + VoxelType::Vertex, + VoxelType::Vertex, + VoxelType::Edge, + VoxelType::Vertex, + VoxelType::Vertex, + VoxelType::Vertex, + VoxelType::Edge, + VoxelType::Edge, + VoxelType::Edge, + VoxelType::Edge, + VoxelType::Face, + VoxelType::Edge, + VoxelType::Edge, + VoxelType::Edge, + VoxelType::Face, + VoxelType::Edge, + VoxelType::Edge, + VoxelType::Edge, + VoxelType::Face, + VoxelType::Edge, + VoxelType::Edge, + VoxelType::Edge, + VoxelType::Face, + VoxelType::Face, + VoxelType::Face, + VoxelType::Face, + VoxelType::Interior, + VoxelType::Empty, +]; + +/// Indicates the convex features of each voxel that can lead to collisions. +/// +/// The interpretation of each bit differs depending on the corresponding voxel type in +/// `FACES_TO_VOXEL_TYPES`: +/// - For `VoxelType::Vertex`: the i-th bit set to `1` indicates that the i-th AABB vertex is convex +/// and might lead to collisions. +/// - For `VoxelType::Edge`: the i-th bit set to `1` indicates that the i-th edge from `Aabb::EDGES_VERTEX_IDS` +/// is convex and might lead to collisions. +/// - For `VoxelType::Face`: the i-th bit set to `1` indicates that the i-th face from `Aabb::FACES_VERTEX_IDS` +/// is exposed and might lead to collisions. +#[cfg(feature = "dim3")] +pub(super) const FACES_TO_FEATURE_MASKS: [u16; 65] = [ + 0b11111111, + 0b10011001, + 0b1100110, + 0b1010101, + 0b110011, + 0b10001, + 0b100010, + 0b10001, + 0b11001100, + 0b10001000, + 0b1000100, + 0b1000100, + 0b10101010, + 0b10001000, + 0b100010, + 0b110000, + 0b1111, + 0b1001, + 0b110, + 0b101, + 0b11, + 0b1, + 0b10, + 0b1, + 0b1100, + 0b1000, + 0b100, + 0b100, + 0b1010, + 0b1000, + 0b10, + 0b100000, + 0b11110000, + 0b10010000, + 0b1100000, + 0b1010000, + 0b110000, + 0b10000, + 0b100000, + 0b10000, + 0b11000000, + 0b10000000, + 0b1000000, + 0b1000000, + 0b10100000, + 0b10000000, + 0b100000, + 0b10000, + 0b111100000000, + 0b100100000000, + 0b11000000000, + 0b1100, + 0b1100000000, + 0b100000000, + 0b1000000000, + 0b1000, + 0b110000000000, + 0b100000000000, + 0b10000000000, + 0b100, + 0b11, + 0b10, + 0b1, + 0b1111111111111111, + 0, +]; + +/// Each octant is assigned three contiguous bits. +#[cfg(feature = "dim3")] +pub(super) const FACES_TO_OCTANT_MASKS: [u32; 65] = [ + 0b1001001001001001001001, + 0b1010010001001010010001, + 0b10001001010010001001010, + 0b10010010010010010010010, + 0b11011001001011011001001, + 0b11111010001011111010001, + 0b111011001010111011001010, + 0b111111010010111111010010, + 0b1001011011001001011011, + 0b1010111011001010111011, + 0b10001011111010001011111, + 0b10010111111010010111111, + 0b11011011011011011011011, + 0b11111111011011111111011, + 0b111011011111111011011111, + 0b111111111111111111111111, + 0b100100100100001001001001, + 0b100110110100001010010001, + 0b110100100110010001001010, + 0b110110110110010010010010, + 0b101101100100011011001001, + 0b101000110100011111010001, + 0b101100110111011001010, + 0b110110111111010010, + 0b100100101101001001011011, + 0b100110000101001010111011, + 0b110100101000010001011111, + 0b110110000000010010111111, + 0b101101101101011011011011, + 0b101000000101011111111011, + 0b101101000111011011111, + 0b111111111111, + 0b1001001001100100100100, + 0b1010010001100110110100, + 0b10001001010110100100110, + 0b10010010010110110110110, + 0b11011001001101101100100, + 0b11111010001101000110100, + 0b111011001010000101100110, + 0b111111010010000000110110, + 0b1001011011100100101101, + 0b1010111011100110000101, + 0b10001011111110100101000, + 0b10010111111110110000000, + 0b11011011011101101101101, + 0b11111111011101000000101, + 0b111011011111000101101000, + 0b111111111111000000000000, + 0b100100100100100100100100, + 0b100110110100100110110100, + 0b110100100110110100100110, + 0b110110110110110110110110, + 0b101101100100101101100100, + 0b101000110100101000110100, + 0b101100110000101100110, + 0b110110000000110110, + 0b100100101101100100101101, + 0b100110000101100110000101, + 0b110100101000110100101000, + 0b110110000000110110000000, + 0b101101101101101101101101, + 0b101000000101101000000101, + 0b101101000000101101000, + 0b0, + 0, +]; + +#[cfg(feature = "dim2")] +pub(super) const FACES_TO_VOXEL_TYPES: [VoxelType; 17] = [ + VoxelType::Vertex, + VoxelType::Vertex, + VoxelType::Vertex, + VoxelType::Face, + VoxelType::Vertex, + VoxelType::Vertex, + VoxelType::Vertex, + VoxelType::Face, + VoxelType::Vertex, + VoxelType::Vertex, + VoxelType::Vertex, + VoxelType::Face, + VoxelType::Face, + VoxelType::Face, + VoxelType::Face, + VoxelType::Interior, + VoxelType::Empty, +]; + +#[cfg(feature = "dim2")] +pub(super) const FACES_TO_FEATURE_MASKS: [u16; 17] = [ + 0b1111, + 0b1001, + 0b110, + 0b1100, + 0b11, + 0b1, + 0b10, + 0b1000, + 0b1100, + 0b1000, + 0b100, + 0b100, + 0b11, + 0b10, + 0b1, + 0b1111111111111111, + 0, +]; + +// NOTE: in 2D we are also using 3 bits per octant even though we technically only need two. +// This keeps some collision-detection easier by avoiding some special-casing. +#[cfg(feature = "dim2")] +pub(super) const FACES_TO_OCTANT_MASKS: [u32; 17] = [ + 0b1001001001, + 0b1011011001, + 0b11001001011, + 0b11011011011, + 0b10010001001, + 0b10000011001, + 0b10001011, + 0b11011, + 0b1001010010, + 0b1011000010, + 0b11001010000, + 0b11011000000, + 0b10010010010, + 0b10000000010, + 0b10010000, + 0b0, + 0, +]; + +// NOTE: this code is used to generate the constant tables +// FACES_TO_VOXEL_TYPES, FACES_TO_FEATURE_MASKS, FACES_TO_OCTANT_MASKS. +#[allow(dead_code)] +#[cfg(feature = "dim2")] +#[cfg(test)] +fn gen_const_tables() { + // The `j-th` bit of `faces_adj_to_vtx[i]` is set to 1, if the j-th face of the AABB (based on + // the face order depicted in `AABB::FACES_VERTEX_IDS`) is adjacent to the `i` vertex of the AABB + // (vertices are indexed as per the diagram depicted in the `FACES_VERTEX_IDS` doc. + // Each entry of this will always have exactly 3 bits set. + let mut faces_adj_to_vtx = [0usize; 4]; + + for fid in 0..4 { + let vids = Aabb::FACES_VERTEX_IDS[fid]; + let key = 1 << fid; + faces_adj_to_vtx[vids.0] |= key; + faces_adj_to_vtx[vids.1] |= key; + } + + /* + * FACES_TO_VOXEL_TYPES + */ + std::println!("const FACES_TO_VOXEL_TYPES: [VoxelType; 17] = ["); + 'outer: for i in 0usize..16 { + // If any vertex of the voxel has three faces with no adjacent voxels, + // then the voxel type is Vertex. + for adjs in faces_adj_to_vtx.iter() { + if (*adjs & i) == 0 { + std::println!("VoxelType::Vertex,"); + continue 'outer; + } + } + + // If one face doesn’t have any adjacent voxel, + // then the voxel type is Face. + for fid in 0..4 { + if ((1 << fid) & i) == 0 { + std::println!("VoxelType::Face,"); + continue 'outer; + } + } + } + + // Add final entries for special values. + std::println!("VoxelType::Interior,"); + std::println!("VoxelType::Empty,"); + std::println!("];"); + + /* + * FACES_TO_FEATURE_MASKS + */ + std::println!("const FACES_TO_FEATURE_MASKS: [u16; 17] = ["); + for i in 0usize..16 { + // Each bit set indicates a convex vertex that can lead to collisions. + // The result will be nonzero only for `VoxelType::Vertex` voxels. + let mut vtx_key = 0; + for (vid, adjs) in faces_adj_to_vtx.iter().enumerate() { + if (*adjs & i) == 0 { + vtx_key |= 1 << vid; + } + } + + if vtx_key != 0 { + std::println!("0b{:b},", vtx_key as u16); + continue; + } + + // Each bit set indicates an exposed face that can lead to collisions. + // The result will be nonzero only for `VoxelType::Face` voxels. + let mut face_key = 0; + for fid in 0..4 { + if ((1 << fid) & i) == 0 { + face_key |= 1 << fid; + } + } + + if face_key != 0 { + std::println!("0b{:b},", face_key as u16); + continue; + } + } + + std::println!("0b{:b},", u16::MAX); + std::println!("0,"); + std::println!("];"); + + /* + * Faces to octant masks. + */ + std::println!("const FACES_TO_OCTANT_MASKS: [u32; 17] = ["); + for i in 0usize..16 { + // First test if we have vertices. + let mut octant_mask = 0; + let mut set_mask = |mask, octant| { + // NOTE: we don’t overwrite any mask already set for the octant. + if (octant_mask >> (octant * 3)) & 0b0111 == 0 { + octant_mask |= mask << (octant * 3); + } + }; + + for (vid, adjs) in faces_adj_to_vtx.iter().enumerate() { + if (*adjs & i) == 0 { + set_mask(1, vid); + } + } + + // This is the index of the normal of the faces given by + // Aabb::FACES_VERTEX_IDS. + const FX: u32 = OctantPattern::FACE_X; + const FY: u32 = OctantPattern::FACE_Y; + const FACE_NORMALS: [u32; 4] = [FX, FX, FY, FY]; + + #[allow(clippy::needless_range_loop)] + for fid in 0..4 { + if ((1 << fid) & i) == 0 { + let vid = Aabb::FACES_VERTEX_IDS[fid]; + let mask = FACE_NORMALS[fid]; + + set_mask(mask, vid.0); + set_mask(mask, vid.1); + } + } + std::println!("0b{:b},", octant_mask); + } + std::println!("0,"); + std::println!("];"); +} + +// NOTE: this code is used to generate the constant tables +// FACES_TO_VOXEL_TYPES, FACES_TO_FEATURE_MASKS, FACES_TO_OCTANT_MASKS. +#[allow(dead_code)] +#[cfg(feature = "dim3")] +#[cfg(test)] +fn gen_const_tables() { + // The `j-th` bit of `faces_adj_to_vtx[i]` is set to 1, if the j-th face of the AABB (based on + // the face order depicted in `AABB::FACES_VERTEX_IDS`) is adjacent to the `i` vertex of the AABB + // (vertices are indexed as per the diagram depicted in the `FACES_VERTEX_IDS` doc. + // Each entry of this will always have exactly 3 bits set. + let mut faces_adj_to_vtx = [0usize; 8]; + + // The `j-th` bit of `faces_adj_to_vtx[i]` is set to 1, if the j-th edge of the AABB (based on + // the edge order depicted in `AABB::EDGES_VERTEX_IDS`) is adjacent to the `i` vertex of the AABB + // (vertices are indexed as per the diagram depicted in the `FACES_VERTEX_IDS` doc. + // Each entry of this will always have exactly 2 bits set. + let mut faces_adj_to_edge = [0usize; 12]; + + for fid in 0..6 { + let vids = Aabb::FACES_VERTEX_IDS[fid]; + let key = 1 << fid; + faces_adj_to_vtx[vids.0] |= key; + faces_adj_to_vtx[vids.1] |= key; + faces_adj_to_vtx[vids.2] |= key; + faces_adj_to_vtx[vids.3] |= key; + } + + #[allow(clippy::needless_range_loop)] + for eid in 0..12 { + let evids = Aabb::EDGES_VERTEX_IDS[eid]; + for fid in 0..6 { + let fvids = Aabb::FACES_VERTEX_IDS[fid]; + if (fvids.0 == evids.0 + || fvids.1 == evids.0 + || fvids.2 == evids.0 + || fvids.3 == evids.0) + && (fvids.0 == evids.1 + || fvids.1 == evids.1 + || fvids.2 == evids.1 + || fvids.3 == evids.1) + { + let key = 1 << fid; + faces_adj_to_edge[eid] |= key; + } + } + } + + /* + * FACES_TO_VOXEL_TYPES + */ + std::println!("const FACES_TO_VOXEL_TYPES: [VoxelType; 65] = ["); + 'outer: for i in 0usize..64 { + // If any vertex of the voxel has three faces with no adjacent voxels, + // then the voxel type is Vertex. + for adjs in faces_adj_to_vtx.iter() { + if (*adjs & i) == 0 { + std::println!("VoxelType::Vertex,"); + continue 'outer; + } + } + + // If any vertex of the voxel has three faces with no adjacent voxels, + // then the voxel type is Edge. + for adjs in faces_adj_to_edge.iter() { + if (*adjs & i) == 0 { + std::println!("VoxelType::Edge,"); + continue 'outer; + } + } + + // If one face doesn’t have any adjacent voxel, + // then the voxel type is Face. + for fid in 0..6 { + if ((1 << fid) & i) == 0 { + std::println!("VoxelType::Face,"); + continue 'outer; + } + } + } + + // Add final entries for special values. + std::println!("VoxelType::Interior,"); + std::println!("VoxelType::Empty,"); + std::println!("];"); + + /* + * FACES_TO_FEATURE_MASKS + */ + std::println!("const FACES_TO_FEATURE_MASKS: [u16; 65] = ["); + for i in 0usize..64 { + // Each bit set indicates a convex vertex that can lead to collisions. + // The result will be nonzero only for `VoxelType::Vertex` voxels. + let mut vtx_key = 0; + for (vid, adjs) in faces_adj_to_vtx.iter().enumerate() { + if (*adjs & i) == 0 { + vtx_key |= 1 << vid; + } + } + + if vtx_key != 0 { + std::println!("0b{:b},", vtx_key as u16); + continue; + } + + // Each bit set indicates a convex edge that can lead to collisions. + // The result will be nonzero only for `VoxelType::Edge` voxels. + let mut edge_key = 0; + for (eid, adjs) in faces_adj_to_edge.iter().enumerate() { + if (*adjs & i) == 0 { + edge_key |= 1 << eid; + } + } + + if edge_key != 0 { + std::println!("0b{:b},", edge_key as u16); + continue; + } + + // Each bit set indicates an exposed face that can lead to collisions. + // The result will be nonzero only for `VoxelType::Face` voxels. + let mut face_key = 0; + for fid in 0..6 { + if ((1 << fid) & i) == 0 { + face_key |= 1 << fid; + } + } + + if face_key != 0 { + std::println!("0b{:b},", face_key as u16); + continue; + } + } + + std::println!("0b{:b},", u16::MAX); + std::println!("0,"); + std::println!("];"); + + /* + * Faces to octant masks. + */ + std::println!("const FACES_TO_OCTANT_MASKS: [u32; 65] = ["); + for i in 0usize..64 { + // First test if we have vertices. + let mut octant_mask = 0; + let mut set_mask = |mask, octant| { + // NOTE: we don’t overwrite any mask already set for the octant. + if (octant_mask >> (octant * 3)) & 0b0111 == 0 { + octant_mask |= mask << (octant * 3); + } + }; + + for (vid, adjs) in faces_adj_to_vtx.iter().enumerate() { + if (*adjs & i) == 0 { + set_mask(1, vid); + } + } + + // This is the index of the axis porting the edges given by + // Aabb::EDGES_VERTEX_IDS. + const EX: u32 = OctantPattern::EDGE_X; + const EY: u32 = OctantPattern::EDGE_Y; + const EZ: u32 = OctantPattern::EDGE_Z; + const EDGE_AXIS: [u32; 12] = [EX, EY, EX, EY, EX, EY, EX, EY, EZ, EZ, EZ, EZ]; + for (eid, adjs) in faces_adj_to_edge.iter().enumerate() { + if (*adjs & i) == 0 { + let vid = Aabb::EDGES_VERTEX_IDS[eid]; + let mask = EDGE_AXIS[eid]; + + set_mask(mask, vid.0); + set_mask(mask, vid.1); + } + } + + // This is the index of the normal of the faces given by + // Aabb::FACES_VERTEX_IDS. + const FX: u32 = OctantPattern::FACE_X; + const FY: u32 = OctantPattern::FACE_Y; + const FZ: u32 = OctantPattern::FACE_Z; + const FACE_NORMALS: [u32; 6] = [FX, FX, FY, FY, FZ, FZ]; + + #[allow(clippy::needless_range_loop)] + for fid in 0..6 { + if ((1 << fid) & i) == 0 { + let vid = Aabb::FACES_VERTEX_IDS[fid]; + let mask = FACE_NORMALS[fid]; + + set_mask(mask, vid.0); + set_mask(mask, vid.1); + set_mask(mask, vid.2); + set_mask(mask, vid.3); + } + } + std::println!("0b{:b},", octant_mask); + } + std::println!("0,"); + std::println!("];"); +} + +#[cfg(test)] +mod test { + #[test] + fn gen_const_tables() { + super::gen_const_tables(); + } +} diff --git a/src/shape/voxels/voxels_edition.rs b/src/shape/voxels/voxels_edition.rs new file mode 100644 index 00000000..100b5204 --- /dev/null +++ b/src/shape/voxels/voxels_edition.rs @@ -0,0 +1,145 @@ +use crate::bounding_volume::Aabb; +use crate::math::{Point, Real, Vector, DIM}; +use crate::shape::voxels::voxels_chunk::{VoxelsChunk, VoxelsChunkHeader}; +use crate::shape::{VoxelState, Voxels}; +use crate::utils::hashmap::Entry; +use alloc::vec; + +impl Voxels { + /// Sets the size of each voxel along each local coordinate axis. + /// + /// Since the internal spatial acceleration structure needs to be updated, this + /// operation runs in `O(n)` time, where `n` is the number of voxels. + pub fn set_voxel_size(&mut self, new_size: Vector) { + let scale = new_size.component_div(&self.voxel_size); + self.chunk_bvh.scale(&scale); + self.voxel_size = new_size; + } + + /// Inserts or remove a voxel from this shape. + /// + /// Return the previous `VoxelState` of this voxel. + pub fn set_voxel(&mut self, key: Point, is_filled: bool) -> VoxelState { + let (chunk_key, id_in_chunk) = Self::chunk_key_and_id_in_chunk(key); + let header_entry = self.chunk_headers.entry(chunk_key); + + if !is_filled && matches!(header_entry, Entry::Vacant(_)) { + // The voxel is already empty (it doesn’t exist at all). + // Nothing more to do. + return VoxelState::EMPTY; + } + + let chunk_header = header_entry.or_insert_with(|| { + let id = self.free_chunks.pop().unwrap_or_else(|| { + self.chunks.push(VoxelsChunk::default()); + self.chunk_keys.push(chunk_key); + self.chunks.len() - 1 + }); + + self.chunk_keys[id] = chunk_key; + self.chunk_bvh + .insert(VoxelsChunk::aabb(&chunk_key, &self.voxel_size), id as u32); + VoxelsChunkHeader { id, len: 0 } + }); + let chunk_id = chunk_header.id; + + let prev = self.chunks[chunk_id].states[id_in_chunk]; + let new_is_empty = !is_filled; + + if prev.is_empty() ^ new_is_empty { + let can_remove_chunk = if new_is_empty { + chunk_header.len -= 1; + chunk_header.len == 0 + } else { + chunk_header.len += 1; + false + }; + + self.chunks[chunk_id].states[id_in_chunk] = + self.update_neighbors_state(key, new_is_empty); + + if can_remove_chunk { + self.chunk_bvh.remove(chunk_id as u32); + let _ = self.chunk_headers.remove(&chunk_key); + self.free_chunks.push(chunk_id); + self.chunk_keys[chunk_id] = VoxelsChunk::INVALID_CHUNK_KEY; + } + } + + prev + } + + /// Crops in-place the voxel shape with a rectangular domain. + /// + /// This removes every voxels out of the `[domain_mins, domain_maxs]` bounds. + pub fn crop(&mut self, domain_mins: Point, domain_maxs: Point) { + // TODO PERF: this could be done more efficiently. + if let Some(new_shape) = self.cropped(domain_mins, domain_maxs) { + *self = new_shape; + } + } + + /// Returns a cropped version of this voxel shape with a rectangular domain. + /// + /// This removes every voxels out of the `[domain_mins, domain_maxs]` bounds. + pub fn cropped(&self, domain_mins: Point, domain_maxs: Point) -> Option { + // TODO PERF: can be optimized significantly. + let mut in_box = vec![]; + for vox in self.voxels() { + if !vox.state.is_empty() + && grid_aabb_contains_point(&domain_mins, &domain_maxs, &vox.grid_coords) + { + in_box.push(vox.grid_coords); + } + } + + if !in_box.is_empty() { + Some(Voxels::new(self.voxel_size, &in_box)) + } else { + None + } + } + + /// Splits this voxels shape into two subshapes. + /// + /// The first subshape contains all the voxels which centers are inside the `aabb`. + /// The second subshape contains all the remaining voxels. + pub fn split_with_box(&self, aabb: &Aabb) -> (Option, Option) { + // TODO PERF: can be optimized significantly. + let mut in_box = vec![]; + let mut rest = vec![]; + for vox in self.voxels() { + if !vox.state.is_empty() { + if aabb.contains_local_point(&vox.center) { + in_box.push(vox.grid_coords); + } else { + rest.push(vox.grid_coords); + } + } + } + + let in_box = if !in_box.is_empty() { + Some(Voxels::new(self.voxel_size, &in_box)) + } else { + None + }; + + let rest = if !rest.is_empty() { + Some(Voxels::new(self.voxel_size, &rest)) + } else { + None + }; + + (in_box, rest) + } +} + +fn grid_aabb_contains_point(mins: &Point, maxs: &Point, point: &Point) -> bool { + for i in 0..DIM { + if point[i] < mins[i] || point[i] > maxs[i] { + return false; + } + } + + true +} diff --git a/src/shape/voxels/voxels_neighborhood.rs b/src/shape/voxels/voxels_neighborhood.rs new file mode 100644 index 00000000..2fde0e73 --- /dev/null +++ b/src/shape/voxels/voxels_neighborhood.rs @@ -0,0 +1,197 @@ +use super::VoxelsChunk; +use crate::math::{Point, Vector, DIM}; +use crate::shape::{VoxelState, Voxels}; + +impl Voxels { + /// Updates the state of the neighbors of the voxel `key`. + /// + /// Modifies the state of the neighbors of `key` to account for it being empty or full. + /// Returns (but doesn’t modify) the new state of the voxel specified by `key`. + #[must_use] + pub(super) fn update_neighbors_state( + &mut self, + key: Point, + center_is_empty: bool, + ) -> VoxelState { + let mut key_data = 0; + + for k in 0..DIM { + let mut left = key; + let mut right = key; + left[k] -= 1; + right[k] += 1; + + // TODO PERF: all the calls to `linear_index` result in a hashmap lookup each time. + // We should instead be smarter and detect if left/right are in the same chunk + // to only look it up once. + if let Some(left_id) = self.linear_index(left) { + let left_state = &mut self.chunks[left_id.chunk_id].states[left_id.id_in_chunk]; + if !left_state.is_empty() { + if center_is_empty { + left_state.0 &= !(1 << (k * 2)); + } else { + left_state.0 |= 1 << (k * 2); + key_data |= 1 << (k * 2 + 1); + } + } + } + + if let Some(right_id) = self.linear_index(right) { + let right_state = &mut self.chunks[right_id.chunk_id].states[right_id.id_in_chunk]; + if !right_state.is_empty() { + if center_is_empty { + right_state.0 &= !(1 << (k * 2 + 1)); + } else { + right_state.0 |= 1 << (k * 2 + 1); + key_data |= 1 << (k * 2); + } + } + } + } + + if center_is_empty { + VoxelState::EMPTY + } else { + VoxelState::new(key_data) + } + } + + pub(super) fn recompute_all_voxels_states(&mut self) { + for (chunk_key, chunk_header) in self.chunk_headers.iter() { + for id_in_chunk in 0..VoxelsChunk::VOXELS_PER_CHUNK { + let voxel_key = VoxelsChunk::voxel_key_at_id(*chunk_key, id_in_chunk as u32); + self.chunks[chunk_header.id].states[id_in_chunk] = + self.compute_voxel_state(voxel_key); + } + } + } + + fn compute_voxel_state(&self, key: Point) -> VoxelState { + let Some(id) = self.linear_index(key) else { + return VoxelState::EMPTY; + }; + + if self.chunks[id.chunk_id].states[id.id_in_chunk].is_empty() { + return VoxelState::EMPTY; + } + + self.compute_voxel_neighborhood_bits(key) + } + + pub(super) fn compute_voxel_neighborhood_bits(&self, key: Point) -> VoxelState { + let mut occupied_faces = 0; + + for k in 0..DIM { + let (mut prev, mut next) = (key, key); + prev[k] -= 1; + next[k] += 1; + + if let Some(next_id) = self.linear_index(next) { + if !self.chunks[next_id.chunk_id].states[next_id.id_in_chunk].is_empty() { + occupied_faces |= 1 << (k * 2); + } + } + if let Some(prev_id) = self.linear_index(prev) { + if !self.chunks[prev_id.chunk_id].states[prev_id.id_in_chunk].is_empty() { + occupied_faces |= 1 << (k * 2 + 1); + } + } + } + + VoxelState::new(occupied_faces) + } + + /// Merges voxel state (neighborhood) information of a given voxel (and all its neighbors) + /// from `self` and `other`, to account for a recent change to the given `voxel` in `self`. + /// + /// This is designed to be called after `self` was modified with [`Voxels::set_voxel`] or + /// [`Voxels::try_set_voxel`]. + /// + /// This is the same as [`Voxels::combine_voxel_states`] but localized to a single voxel and its + /// neighbors. + pub fn propagate_voxel_change( + &mut self, + other: &mut Self, + voxel: Point, + origin_shift: Vector, + ) { + let center_is_empty = self + .voxel_state(voxel) + .map(|vox| vox.is_empty()) + .unwrap_or(true); + let center_state_delta = + other.update_neighbors_state(voxel - origin_shift, center_is_empty); + + if let Some(vid) = self.linear_index(voxel) { + self.chunks[vid.chunk_id].states[vid.id_in_chunk].0 |= center_state_delta.0; + } + } + + /// Merges voxel state (neighborhood) information of each voxel from `self` and `other`. + /// + /// This allows each voxel from one shape to be aware of the presence of neighbors belonging to + /// the other so that collision detection is capable of transitioning between the boundaries of + /// one shape to the other without hitting an internal edge. + /// + /// Both voxels shapes are assumed to have the same [`Self::voxel_size`]. + /// If `other` lives in a coordinate space with a different origin than `self`, then + /// `origin_shift` represents the distance (as a multiple of the `voxel_size`) from the origin + /// of `self` to the origin of `other`. Therefore, a voxel with coordinates `key` on `other` + /// will have coordinates `key + origin_shift` on `self`. + pub fn combine_voxel_states(&mut self, other: &mut Self, origin_shift: Vector) { + let one = Vector::repeat(1); + + // Intersect the domains + 1 cell. + let [my_domain_mins, my_domain_maxs] = self.domain(); + let [other_domain_mins, other_domain_maxs] = other.domain(); + let d0 = [my_domain_mins - one, my_domain_maxs + one * 2]; + let d1 = [ + other_domain_mins - one + origin_shift, + other_domain_maxs + one * 2 + origin_shift, + ]; + + let d01 = [d0[0].sup(&d1[0]), d0[1].inf(&d1[1])]; + // Iterate on the domain intersection. If the voxel exists (and is non-empty) on both shapes, we + // simply need to combine their bitmasks. If it doesn’t exist on both shapes, we need to + // actually check the neighbors. + // + // The `domain` is expressed in the grid coordinate space of `self`. + for i in d01[0].x..d01[1].x { + for j in d01[0].y..d01[1].y { + #[cfg(feature = "dim2")] + let k_range = 0..1; + #[cfg(feature = "dim3")] + let k_range = d01[0].z..d01[1].z; + for _k in k_range { + #[cfg(feature = "dim2")] + let key0 = Point::new(i, j); + #[cfg(feature = "dim3")] + let key0 = Point::new(i, j, _k); + let key1 = key0 - origin_shift; + let vox0 = self + .linear_index(key0) + .map(|id| &mut self.chunks[id.chunk_id].states[id.id_in_chunk]) + .filter(|state| !state.is_empty()); + let vox1 = other + .linear_index(key1) + .map(|id| &mut other.chunks[id.chunk_id].states[id.id_in_chunk]) + .filter(|state| !state.is_empty()); + + match (vox0, vox1) { + (Some(vox0), Some(vox1)) => { + vox0.0 |= vox1.0; + vox1.0 |= vox0.0; + } + (Some(vox0), None) => { + vox0.0 |= other.compute_voxel_neighborhood_bits(key1).0; + } + (None, Some(vox1)) => { + vox1.0 |= self.compute_voxel_neighborhood_bits(key0).0; + } + (None, None) => { /* Nothing to adjust. */ } + } + } + } + } + } +} diff --git a/src/shape/voxels/voxels_sizeof.rs b/src/shape/voxels/voxels_sizeof.rs new file mode 100644 index 00000000..08a9aa3f --- /dev/null +++ b/src/shape/voxels/voxels_sizeof.rs @@ -0,0 +1,29 @@ +use super::{Voxels, VoxelsChunk, VoxelsChunkHeader}; +use crate::math::Point; + +impl Voxels { + // TODO: support a crate like get_size2 (will require support on nalgebra too)? + /// An approximation of the memory usage (in bytes) for this struct plus + /// the memory it allocates dynamically. + pub fn total_memory_size(&self) -> usize { + size_of::() + self.heap_memory_size() + } + + /// An approximation of the memory dynamically-allocated by this struct. + pub fn heap_memory_size(&self) -> usize { + // NOTE: if a new field is added to `Self`, adjust this function result. + let Self { + chunk_bvh, + chunk_headers, + chunks, + free_chunks, + chunk_keys, + voxel_size: _, + } = self; + chunks.capacity() * size_of::() + + free_chunks.capacity() * size_of::() + + chunk_keys.capacity() * size_of::>() + + chunk_headers.capacity() * (size_of::() + size_of::>()) + + chunk_bvh.heap_memory_size() + } +} diff --git a/src/transformation/voxelization/voxelized_volume.rs b/src/transformation/voxelization/voxelized_volume.rs index b8745869..f57dea5a 100644 --- a/src/transformation/voxelization/voxelized_volume.rs +++ b/src/transformation/voxelization/voxelized_volume.rs @@ -112,7 +112,7 @@ pub enum VoxelValue { } #[derive(Copy, Clone, Debug, PartialEq, Eq)] -struct VoxelState { +struct VoxelData { #[cfg(feature = "dim2")] multiplicity: u32, num_primitive_intersections: u32, @@ -124,7 +124,7 @@ pub struct VoxelizedVolume { scale: Real, resolution: [u32; DIM], values: Vec, - data: Vec, + data: Vec, primitive_intersections: Vec<(u32, u32)>, } @@ -514,7 +514,7 @@ impl VoxelizedVolume { .resize(len as usize, VoxelValue::PrimitiveUndefined); self.data.resize( len as usize, - VoxelState { + VoxelData { #[cfg(feature = "dim2")] multiplicity: 0, num_primitive_intersections: 0, From fe7b12e0d838316ea44dc1ef0df6b9776f2bd76f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Crozet?= Date: Fri, 3 Oct 2025 17:57:03 +0200 Subject: [PATCH 2/3] fix compilation of docs and examples --- src/shape/voxels/voxels_consts.rs | 2 +- src/shape/voxels/voxels_edition.rs | 5 +++++ src/shape/voxels/voxels_neighborhood.rs | 3 +-- 3 files changed, 7 insertions(+), 3 deletions(-) diff --git a/src/shape/voxels/voxels_consts.rs b/src/shape/voxels/voxels_consts.rs index 9d8ad7dc..b2adf885 100644 --- a/src/shape/voxels/voxels_consts.rs +++ b/src/shape/voxels/voxels_consts.rs @@ -1,7 +1,7 @@ use super::VoxelType; #[cfg(test)] -use {crate::bounding_volume::Aabb, OctantPattern}; +use {crate::bounding_volume::Aabb, super::OctantPattern}; // Index to the item of FACES_TO_VOXEL_TYPES which identifies interior voxels. #[cfg(feature = "dim2")] diff --git a/src/shape/voxels/voxels_edition.rs b/src/shape/voxels/voxels_edition.rs index 100b5204..c4538a1a 100644 --- a/src/shape/voxels/voxels_edition.rs +++ b/src/shape/voxels/voxels_edition.rs @@ -60,7 +60,12 @@ impl Voxels { if can_remove_chunk { self.chunk_bvh.remove(chunk_id as u32); + + #[cfg(feature = "enhanced-determinism")] + let _ = self.chunk_headers.swap_remove(&chunk_key); + #[cfg(not(feature = "enhanced-determinism"))] let _ = self.chunk_headers.remove(&chunk_key); + self.free_chunks.push(chunk_id); self.chunk_keys[chunk_id] = VoxelsChunk::INVALID_CHUNK_KEY; } diff --git a/src/shape/voxels/voxels_neighborhood.rs b/src/shape/voxels/voxels_neighborhood.rs index 2fde0e73..b18eb589 100644 --- a/src/shape/voxels/voxels_neighborhood.rs +++ b/src/shape/voxels/voxels_neighborhood.rs @@ -104,8 +104,7 @@ impl Voxels { /// Merges voxel state (neighborhood) information of a given voxel (and all its neighbors) /// from `self` and `other`, to account for a recent change to the given `voxel` in `self`. /// - /// This is designed to be called after `self` was modified with [`Voxels::set_voxel`] or - /// [`Voxels::try_set_voxel`]. + /// This is designed to be called after `self` was modified with [`Voxels::set_voxel`]. /// /// This is the same as [`Voxels::combine_voxel_states`] but localized to a single voxel and its /// neighbors. From cdca8cd588575f1e298dc60187d70c19ea0149bb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Crozet?= Date: Fri, 3 Oct 2025 17:57:47 +0200 Subject: [PATCH 3/3] =?UTF-8?q?chore:=E2=80=AFcargo=20fmt?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/shape/voxels/voxels_consts.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/shape/voxels/voxels_consts.rs b/src/shape/voxels/voxels_consts.rs index b2adf885..a0b75560 100644 --- a/src/shape/voxels/voxels_consts.rs +++ b/src/shape/voxels/voxels_consts.rs @@ -1,7 +1,7 @@ use super::VoxelType; #[cfg(test)] -use {crate::bounding_volume::Aabb, super::OctantPattern}; +use {super::OctantPattern, crate::bounding_volume::Aabb}; // Index to the item of FACES_TO_VOXEL_TYPES which identifies interior voxels. #[cfg(feature = "dim2")]