Polish collision box reshaping

common/src/collision.rs

@@ -1,3 +1,4 @@
+use crate::math;
 use crate::node::DualGraph;
 use crate::{dodeca::Vertex, graph::NodeId};
 use std::fmt;
@@ -164,46 +165,54 @@ impl ChunkBoundingBox {
         radius: f64,
         dimension: u8,
     ) -> Option<Self> {
-        let euclidean_position = {
-            let temp = chunk.node_to_chunk() * translated_position;
-            temp.xyz() / temp[3]
-        };
-
-        let mut min_xyz = na::Vector3::<u32>::new(0_u32, 0_u32, 0_u32);
-        let mut max_xyz = na::Vector3::<u32>::new(0_u32, 0_u32, 0_u32);
-
-        // It's important to note that euclidean_position is measured as chunk lengths, and radius is measured in absolute units.
-        // By coincidence, an absolute unit is approximately a chunk's diameter, and only because of that there is no unit conversion here.
-
-        // verify at least one box corner is within the chunk
-        if euclidean_position
-            .iter()
-            .all(|n| n + radius > 0_f64 && n - radius < 1_f64)
-        {
-            min_xyz.x =
-                ((euclidean_position.x - radius).max(0_f64) * dimension as f64).floor() as u32;
-            max_xyz.x =
-                ((euclidean_position.x + radius).min(1_f64) * dimension as f64).ceil() as u32;
-
-            min_xyz.y =
-                ((euclidean_position.y - radius).max(0_f64) * dimension as f64).floor() as u32;
-            max_xyz.y =
-                ((euclidean_position.y + radius).min(1_f64) * dimension as f64).ceil() as u32;
-
-            min_xyz.z =
-                ((euclidean_position.z - radius).max(0_f64) * dimension as f64).floor() as u32;
-            max_xyz.z =
-                ((euclidean_position.z + radius).min(1_f64) * dimension as f64).ceil() as u32;
-            Some(ChunkBoundingBox {
-                node,
-                chunk,
-                min_xyz,
-                max_xyz,
-                dimension,
-            })
-        } else {
-            None
+        // The following code computes the bounds of the minimum-size bounding box that contains a sphere of the given
+        // radius centered at translated_position. Computing this requires solving a quadratic equation, whose solution
+        // was used to write this code. Not all intermediate values used in the computation have any intrinsic meaning.
+
+        // Position of entity relative to chunk corner
+        let chunk_position = chunk.node_to_dual() * translated_position;
+
+        let sinh_radius = radius.sinh();
+        let common_factor = 1.0 / (chunk_position[3].powi(2) + sinh_radius.powi(2))
+            * Vertex::dual_to_chunk_factor();
+
+        let mut voxel_min = [0.0; 3];
+        let mut voxel_max = [0.0; 3];
+        let dimension_float = dimension as f64;
+        for i in 0..3 {
+            let aabb_center = chunk_position[i] * chunk_position[3];
+            let aabb_width = sinh_radius
+                * (chunk_position[3].powi(2) + sinh_radius.powi(2) - chunk_position[i].powi(2))
+                    .sqrt();
+
+            voxel_min[i] = ((1.0 - (aabb_center + aabb_width) * common_factor).max(0.0)
+                * dimension_float)
+                .floor();
+
+            voxel_max[i] = ((1.0 - (aabb_center - aabb_width) * common_factor).min(1.0)
+                * dimension_float)
+                .ceil();
+
+            if voxel_min[i] >= dimension_float || voxel_max[i] <= 0.0 {
+                return None;
+            }
         }
+
+        Some(ChunkBoundingBox {
+            node,
+            chunk,
+            min_xyz: na::Vector3::<u32>::new(
+                voxel_min[0] as u32,
+                voxel_min[1] as u32,
+                voxel_min[2] as u32,
+            ),
+            max_xyz: na::Vector3::<u32>::new(
+                voxel_max[0] as u32,
+                voxel_max[1] as u32,
+                voxel_max[2] as u32,
+            ),
+            dimension,
+        })
     }
 
     pub fn every_voxel(&self) -> impl Iterator<Item = u32> + '_ {
@@ -352,11 +361,11 @@ mod tests {
                 // Getting the correct estimation for the number of voxels can be tricky
                 let expected_voxel_count = (radius * 2.0 * CHUNK_SIZE_F).powf(3.0); // value to display
                 let minimum_expected_voxel_count =
-                    ((((radius * CHUNK_SIZE_F) - 1_f64).powf(3.0) / margin_of_error).floor() * 8_f64 )
-                        as i32;
+                    ((((radius * CHUNK_SIZE_F) - 1_f64).powf(3.0) / margin_of_error).floor()
+                        * 8_f64) as i32;
                 let maximum_expected_voxel_count =
-                    ((((radius * CHUNK_SIZE_F) + 1_f64).powf(3.0)  * margin_of_error).ceil() * 20_f64)
-                        as i32;
+                    ((((radius * CHUNK_SIZE_F) + 1_f64).powf(3.0) * margin_of_error).ceil()
+                        * 20_f64) as i32;
 
                 let position = central_chunk.chunk_to_node()
                     * na::Vector4::new(

common/src/dodeca.rs

@@ -49,6 +49,12 @@ impl Side {
         SIDE_VERTICES[self as usize]
     }
 
+    /// Outward normal vector of this side
+    #[inline]
+    pub fn normal(self) -> &'static na::Vector4<f64> {
+        &SIDE_NORMALS[self as usize]
+    }
+
     /// Reflection across this side
     #[inline]
     pub fn reflection(self) -> &'static na::Matrix4<f64> {
@@ -139,21 +145,38 @@ impl Vertex {
 
     /// Transform from euclidean chunk coordinates to hyperbolic node space
     pub fn chunk_to_node(self) -> na::Matrix4<f64> {
-        let origin = na::Vector4::new(0.0, 0.0, 0.0, 1.0);
         let [a, b, c] = self.canonical_sides();
         na::Matrix4::from_columns(&[
-            a.reflection().column(3) - origin,
-            b.reflection().column(3) - origin,
-            c.reflection().column(3) - origin,
-            origin,
-        ]) * na::Matrix4::new_scaling(0.5)
+            a.normal() * a.normal()[3],
+            b.normal() * b.normal()[3],
+            c.normal() * c.normal()[3],
+            math::origin(),
+        ])
     }
 
     /// Transform from hyperbolic node space to euclidean chunk coordinates
     pub fn node_to_chunk(self) -> na::Matrix4<f64> {
         self.chunk_to_node().try_inverse().unwrap()
     }
 
+    /// Transform from cube-centric coordinates to dodeca-centric coordinates
+    pub fn dual_to_node(self) -> &'static na::Matrix4<f64> {
+        &DUAL_TO_NODE[self as usize]
+    }
+
+    /// Transform from dodeca-centric coordinates to cube-centric coordinates
+    pub fn node_to_dual(self) -> &'static na::Matrix4<f64> {
+        &NODE_TO_DUAL[self as usize]
+    }
+
+    /// Scale factor used in conversion from cube-centric coordinates to euclidean chunk coordinates.
+    /// Scaling the x, y, and z components of a vector in cube-centric coordinates by this value
+    /// and dividing them by the w coordinate will yield chunk coordinates, but reflected such that 0
+    /// and 1 are swapped for each coordinate.
+    pub fn dual_to_chunk_factor() -> f64 {
+        (2.0 + 5.0f64.sqrt()).sqrt()
+    }
+
     /// Convenience method for `self.chunk_to_node().determinant() < 0`.
     pub fn parity(self) -> bool {
         CHUNK_TO_NODE_PARITY[self as usize]
@@ -180,33 +203,33 @@ lazy_static! {
         result
     };
 
-    /// Transform that moves from a neighbor to a reference node, for each side
-    static ref REFLECTIONS: [na::Matrix4<f64>; SIDE_COUNT] = {
+    /// Vector corresponding to the outer normal of each side
+    static ref SIDE_NORMALS: [na::Vector4<f64>; SIDE_COUNT] = {
         let phi = 1.25f64.sqrt() + 0.5; // golden ratio
-        let root_phi = phi.sqrt();
-        let f = math::lorentz_normalize(&na::Vector4::new(root_phi, phi * root_phi, 0.0, phi + 2.0));
+        let f = math::lorentz_normalize(&na::Vector4::new(1.0, phi, 0.0, phi.sqrt()));
 
-        let mut result = [na::zero(); SIDE_COUNT];
+        let mut result: [na::Vector4<f64>; SIDE_COUNT] = [na::zero(); SIDE_COUNT];
         let mut i = 0;
         for (x, y, z, w) in [
             (f.x, f.y, f.z, f.w),
             (-f.x, f.y, -f.z, f.w),
             (f.x, -f.y, -f.z, f.w),
             (-f.x, -f.y, f.z, f.w),
         ]
-            .iter()
-            .cloned()
         {
-            for (x, y, z, w) in [(x, y, z, w), (y, z, x, w), (z, x, y, w)].iter().cloned() {
-                result[i] = math::translate(&math::origin(), &na::Vector4::new(x, y, z, w))
-                    * math::euclidean_reflect(&na::Vector4::new(x, y, z, 0.0))
-                    * math::translate(&math::origin(), &na::Vector4::new(-x, -y, -z, w));
+            for (x, y, z, w) in [(x, y, z, w), (y, z, x, w), (z, x, y, w)] {
+                result[i] = na::Vector4::new(x, y, z, w);
                 i += 1;
             }
         }
         result
     };
 
+    /// Transform that moves from a neighbor to a reference node, for each side
+    static ref REFLECTIONS: [na::Matrix4<f64>; SIDE_COUNT] = {
+        SIDE_NORMALS.map(|r| math::reflect(&r))
+    };
+
     /// Sides incident to a vertex, in canonical order
     static ref VERTEX_SIDES: [[Side; 3]; VERTEX_COUNT] = {
         let mut result = [[Side::A; 3]; VERTEX_COUNT];
@@ -227,6 +250,25 @@ lazy_static! {
         result
     };
 
+    /// Transform that converts from cube-centric coordinates to dodeca-centric coordinates
+    static ref DUAL_TO_NODE: [na::Matrix4<f64>; VERTEX_COUNT] = {
+        let mip_origin_normal = math::mip(&math::origin(), &SIDE_NORMALS[0]); // This value is the same for every side
+        let mut result = [na::zero(); VERTEX_COUNT];
+        for i in 0..VERTEX_COUNT {
+            let [a, b, c] = VERTEX_SIDES[i];
+            let vertex_position = math::lorentz_normalize(
+                &(math::origin() - (a.normal() + b.normal() + c.normal()) * mip_origin_normal),
+            );
+            result[i] = na::Matrix4::from_columns(&[-a.normal(), -b.normal(), -c.normal(), vertex_position]);
+        }
+        result
+    };
+
+    /// Transform that converts from dodeca-centric coordinates to cube-centric coordinates
+    static ref NODE_TO_DUAL: [na::Matrix4<f64>; VERTEX_COUNT] = {
+        DUAL_TO_NODE.map(|m| m.try_inverse().unwrap())
+    };
+
     /// Vertex shared by 3 sides
     static ref SIDES_TO_VERTEX: [[[Option<Vertex>; SIDE_COUNT]; SIDE_COUNT]; SIDE_COUNT] = {
         let mut result = [[[None; SIDE_COUNT]; SIDE_COUNT]; SIDE_COUNT];