Add license and readme

- Add custom error type
- Update more mesh based calculations to f32

Cargo.toml

@@ -5,3 +5,5 @@ edition = "2021"
 
 [dependencies]
 anyhow = "*"
+vecmath = "*"
+itertools = "*"

LICENSE

@@ -0,0 +1,22 @@
+MIT License
+
+Copyright (c) [year] [fullname]
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+

README.md

@@ -0,0 +1,4 @@
+# Roctree
+
+This is planned to be an Octree rust crate enabling adaptive refinement for finite element analysis of voxel meshes.
+The plan is to also ship a tool that can convert an STL into a voxel volume mesh for further processing by finite element analysis codes.

src/bvh.rs

@@ -0,0 +1,263 @@
+use core::f64;
+
+use crate::{errors::RoctreeError, primitives::Box3, stl::parser::StlSolid};
+
+/// Number of potential split planes to check
+const NBINS: usize = 8;
+
+/// Struct to hold bin information, this includes an axis-aligned bounding box as well as the count
+/// of triangles with the
+#[derive(Default, Clone)]
+struct Bin {
+    aabb: Box3,
+    tri_count: usize,
+}
+
+/// A node representing a node single in the bounding volume hierarchy
+#[derive(Default)]
+struct BvhNode {
+    aabb: Box3,
+    left_child: usize,
+    right_child: usize,
+    prim_idx: usize,
+    prim_count: usize,
+}
+
+impl BvhNode {
+    /// This computes the computational cost associated with the specific node.
+    /// The cost is an aggregate that considers the size of the bounding box (how likely it will
+    /// intersect a ray) * the number of triangles (the number of ray-triangle checks we have to
+    /// perform if an intersect occurs)
+    fn calculate_node_cost(&self) -> f32 {
+        self.aabb.surface_area() * self.prim_count as f32
+    }
+}
+
+/// Struct representing the total bounding volume hierarchy
+pub struct Bvh {
+    nodes: Vec<BvhNode>,
+    tri_idx: Vec<usize>,
+}
+
+/// Return type for get_subdivision_position_and_axis
+struct BestSplitPlane {
+    axis: usize,
+    pos: f64,
+    cost: f64,
+}
+
+impl Bvh {
+    /// Create a new Bounding volume hierarchy from an STL solid.
+    pub fn new(solid: &StlSolid) -> Result<Bvh, RoctreeError> {
+        // Initialize the root node
+        let mut nodes = vec![BvhNode {
+            prim_idx: 0,
+            prim_count: solid.triangles.len(),
+            ..Default::default()
+        }];
+        nodes[0]
+            .aabb
+            .shrink_wrap_primitives(&solid.vertices, &solid.triangles);
+
+        let mut bvh = Bvh {
+            nodes,
+            tri_idx: Vec::from_iter(0..solid.triangles.len()),
+        };
+
+        bvh.subdivide(0, solid)?;
+        Ok(bvh)
+    }
+
+    /// Subdivide a node, subdivision doesn't occur if the node has only 2 primitives
+    pub fn subdivide(&mut self, node_idx: usize, solid: &StlSolid) -> Result<(), RoctreeError> {
+        let BestSplitPlane {
+            cost,
+            axis: split_axis,
+            pos: split_position,
+        } = self.get_best_split(node_idx, solid)?;
+
+        // If the split cost is greater than keeping the node intact just return
+        if cost > self.nodes[node_idx].calculate_node_cost() as f64 {
+            return Ok(());
+        }
+
+        let prim_count = self.nodes[node_idx].prim_count;
+        let prim_idx = self.nodes[node_idx].prim_idx;
+        // Reorganize the triangles so the triangles are ordered, this allows triangles belonging
+        // to a specific node to be adjacent
+        let mut start_idx = prim_idx;
+        let mut end_idx = prim_idx + prim_count - 1;
+        while start_idx <= end_idx {
+            if solid.triangle_centroids[self.tri_idx[start_idx]][split_axis] < split_position as f32
+            {
+                start_idx += 1;
+            } else {
+                self.tri_idx.swap(start_idx, end_idx);
+                end_idx -= 1;
+            }
+        }
+
+        // Now split the primitives into two nodes
+        let left_count = start_idx - prim_idx;
+        let right_count = prim_count - left_count;
+
+        // If all primitives belong to either the left or right half exit because the parent node
+        // is a leaf
+        if left_count == 0 || left_count == prim_count {
+            return Ok(());
+        }
+
+        // Otherwise split the nodes into left and right child
+        let left_child = self.nodes.len();
+        let right_child = self.nodes.len() + 1;
+
+        // FIX: Having to copy the [usize; 3] that defines each triangle just to create the
+        // axis-aligned bounding box seems wasteful although it makes the code interface a lot
+        // cleaner. Maybe spend some time to redesign this part to make it cleaner once we really
+        // start focusing on performance tuning
+        let left_prim_triangles: Vec<_> = (prim_idx..left_count)
+            .map(|idx| solid.triangles[idx])
+            .collect();
+        let mut aabb = Box3::default();
+        aabb.shrink_wrap_primitives(&solid.vertices, &left_prim_triangles);
+
+        self.nodes.push(BvhNode {
+            aabb,
+            left_child: 0,
+            right_child: 0,
+            prim_idx,
+            prim_count: left_count,
+        });
+
+        let right_prim_triangles: Vec<_> = ((prim_idx + left_count)..prim_count)
+            .map(|idx| solid.triangles[idx])
+            .collect();
+        let mut aabb = Box3::default();
+        aabb.shrink_wrap_primitives(&solid.vertices, &right_prim_triangles);
+
+        self.nodes.push(BvhNode {
+            aabb,
+            left_child: 0,
+            right_child: 0,
+            prim_idx: prim_idx + left_count,
+            prim_count: right_count,
+        });
+        self.nodes[node_idx].prim_count = 0;
+        self.nodes[node_idx].left_child = left_child;
+        self.nodes[node_idx].right_child = right_child;
+
+        // Recurse to divide children
+        self.subdivide(left_child, solid)?;
+        self.subdivide(right_child, solid)?;
+
+        Ok(())
+    }
+
+    /// Get the best split plane based on the surface area heuristic
+    /// This procedure works by binning
+    fn get_best_split(
+        &self,
+        node_idx: usize,
+        solid: &StlSolid,
+    ) -> Result<BestSplitPlane, RoctreeError> {
+        // Initialize the best split plane information with bad values so we know if no plane gets
+        // assigned
+        let mut best_cost = f64::INFINITY;
+        let mut best_pos = f64::NAN;
+
+        // Pick the longest axis as the splitting axis
+        let axis = self.nodes[node_idx]
+            .aabb
+            .get_extents()
+            .iter()
+            .enumerate()
+            .fold((0_usize, &f32::NEG_INFINITY), |acc, a| {
+                if acc.1 > a.1 {
+                    acc
+                } else {
+                    a
+                }
+            })
+            .0;
+
+        let node = &self.nodes[node_idx];
+        // Get minimum and maximum bounds over the centroids off the triangles. This improves
+        // the binning compared to doing the binning over the vertices
+        let (min_bd, max_bd) = self.tri_idx[node.prim_idx..node.prim_idx + node.prim_count]
+            .iter()
+            .fold((f64::INFINITY, f64::NEG_INFINITY), |bounds, tri_idx| {
+                (
+                    bounds
+                        .0
+                        .min(solid.triangle_centroids[*tri_idx][axis].into()),
+                    bounds
+                        .1
+                        .max(solid.triangle_centroids[*tri_idx][axis].into()),
+                )
+            });
+        if min_bd == max_bd {
+            return Err(RoctreeError::SplitPlaneError);
+        };
+
+        let mut bins = vec![Bin::default(); NBINS];
+        let scale = NBINS as f64 / (max_bd - min_bd);
+
+        // Assign the triangles to each of the bins
+        for tri_idx in self.tri_idx[node.prim_idx..node.prim_idx + node.prim_count].iter() {
+            let bin_index = usize::min(
+                NBINS - 1,
+                ((solid.triangle_centroids[*tri_idx][axis] as f64 - min_bd) * scale) as usize,
+            );
+
+            bins[bin_index].tri_count += 1;
+            bins[bin_index]
+                .aabb
+                .grow(solid.vertices[solid.triangles[*tri_idx][0]]);
+            bins[bin_index]
+                .aabb
+                .grow(solid.vertices[solid.triangles[*tri_idx][1]]);
+            bins[bin_index]
+                .aabb
+                .grow(solid.vertices[solid.triangles[*tri_idx][2]]);
+        }
+
+        // Now get the left and right counts for each split plane (based on the bins)
+        let mut left_area = [0.0; NBINS - 1];
+        let mut right_area = [0.0; NBINS - 1];
+        let mut left_count = [0; NBINS - 1];
+        let mut right_count = [0; NBINS - 1];
+        let mut left_box = Box3::default();
+        let mut right_box = Box3::default();
+        let mut left_sum = 0;
+        let mut right_sum = 0;
+        for plane_idx in 0..NBINS - 1 {
+            // Calculate the left side of the plane
+            left_sum += bins[plane_idx].tri_count;
+            left_count[plane_idx] = left_sum;
+            left_box.grow_by_other_box(&bins[plane_idx].aabb);
+            left_area[plane_idx] = left_box.surface_area();
+
+            // Calculate the right side of the plane
+            right_sum += bins[NBINS - 1 - plane_idx].tri_count;
+            right_count[NBINS - 2 - plane_idx] = right_sum;
+            right_box.grow_by_other_box(&bins[NBINS - 1 - plane_idx].aabb);
+            right_area[NBINS - 2 - plane_idx] = right_box.surface_area();
+        }
+        let scale = 1.0 / scale;
+        for plane_idx in 0..NBINS - 1 {
+            let plane_cost: f64 = (left_count[plane_idx] as f32 * left_area[plane_idx]
+                + right_count[plane_idx] as f32 * right_area[plane_idx])
+                .into();
+            if plane_cost < best_cost {
+                best_cost = plane_cost;
+                best_pos = scale * plane_idx as f64;
+            }
+        }
+
+        Ok(BestSplitPlane {
+            axis,
+            pos: best_pos,
+            cost: best_cost,
+        })
+    }
+}

src/errors.rs

@@ -0,0 +1,19 @@
+use std::{error::Error, fmt::Display};
+
+#[derive(Debug)]
+pub enum RoctreeError {
+    SplitPlaneError,
+}
+
+impl Error for RoctreeError {}
+
+impl Display for RoctreeError {
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        match self {
+            RoctreeError::SplitPlaneError => write!(
+                f,
+                "Failed to get split plane, centroids overlap along longest axis"
+            ),
+        }
+    }
+}

src/lib.rs

@@ -1,2 +1,6 @@
+pub mod bvh;
+pub mod errors;
+pub mod octree;
+pub mod primitives;
 pub mod stl;
 pub mod utilities;

src/octree/base.rs

@@ -0,0 +1,9 @@
+pub struct Octree {
+    _octants: Vec<Octant>,
+    _min_size: f64,
+}
+
+pub struct Octant {
+    _morton_id: usize,
+    _level: usize,
+}

src/octree/mod.rs

@@ -0,0 +1,2 @@
+pub mod base;
+pub mod morton_ids;

src/octree/morton_ids.rs

@@ -0,0 +1,17 @@
+/// Calculate the morton id in 3d space based on integer position
+pub fn encode(mut x: usize, mut y: usize, mut z: usize) {
+    let mut id = 0;
+    let mut level = 0;
+    while x > 0 || y > 0 || z > 0 {
+        let xbit = x % 2;
+        let ybit = x % 2;
+        let zbit = x % 2;
+        let curr_code = zbit + 2 * ybit + 4 * xbit;
+
+        id += curr_code << (3 * level);
+        x /= 2;
+        y /= 2;
+        z /= 2;
+        level += 1;
+    }
+}

src/primitives.rs

@@ -0,0 +1,91 @@
+use vecmath::vec3_sub;
+
+#[derive(Default, Clone)]
+pub struct Box3 {
+    pub min_bd: [f32; 3],
+    pub max_bd: [f32; 3],
+}
+
+impl Box3 {
+    /// Create a new box from the bounding values
+    pub fn new_from_bounds(min_bd: [f32; 3], max_bd: [f32; 3]) -> Self {
+        Box3 { min_bd, max_bd }
+    }
+
+    /// This function returns the 8 vertices of an axis-aligned bounding box.
+    /// These are ordered
+    ///     7----6
+    /// |  /|   /|
+    /// | 4----5
+    /// | | 3--| 2
+    /// | |/   |/
+    /// | 0----1
+    pub fn get_vertices(&self) -> Vec<[f32; 3]> {
+        vec![
+            self.min_bd,
+            [self.max_bd[0], self.min_bd[1], self.min_bd[2]],
+            [self.max_bd[0], self.max_bd[1], self.min_bd[2]],
+            [self.min_bd[0], self.max_bd[1], self.min_bd[2]],
+            [self.max_bd[0], self.min_bd[1], self.max_bd[2]],
+            [self.max_bd[0], self.max_bd[1], self.max_bd[2]],
+            self.max_bd,
+            [self.min_bd[0], self.max_bd[1], self.max_bd[2]],
+        ]
+    }
+
+    /// Expand a box to fit all primitive inside of it.
+    /// NOTE: vertices does not need to include only vertices being considered. Instead only
+    /// vertices included as indices in the primitives list are considered
+    pub fn shrink_wrap_primitives<T>(&mut self, vertices: &[[f32; 3]], primitives: &[T])
+    where
+        for<'a> &'a T: IntoIterator<Item = &'a usize>,
+    {
+        // If the current primitive list is empty, zero bounding box and exit
+        if primitives.is_empty() {
+            self.min_bd = [0.0; 3];
+            self.max_bd = [0.0; 3];
+            return;
+        }
+
+        for vertex_idx in primitives.iter().flatten() {
+            self.grow(vertices[*vertex_idx]);
+        }
+    }
+
+    /// Compute the surface area of a box
+    pub fn surface_area(&self) -> f32 {
+        let extents = vec3_sub(self.max_bd, self.min_bd);
+        extents[0] * extents[1] + extents[1] * extents[2] + extents[2] * extents[0]
+    }
+
+    /// Get the extents of the current box
+    #[inline(always)]
+    pub fn get_extents(&self) -> [f32; 3] {
+        vec3_sub(self.max_bd, self.min_bd)
+    }
+
+    /// Grow the axis aligned bounding box to include a specific vertex
+    #[inline(always)]
+    pub fn grow(&mut self, vertex: [f32; 3]) {
+        for (axis, pos) in vertex.iter().enumerate() {
+            self.min_bd[axis] = self.min_bd[axis].min(*pos);
+            self.max_bd[axis] = self.max_bd[axis].min(*pos);
+        }
+    }
+
+    /// Grow the axis aligned bounding box to include a min and max boundary
+    /// This is useful when trying to combine two bounding boxes
+    #[inline(always)]
+    pub fn grow_by_other_box(&mut self, other: &Self) {
+        for axis in 0..3 {
+            self.min_bd[axis] = self.min_bd[axis].min(other.min_bd[axis]);
+            self.max_bd[axis] = self.max_bd[axis].max(other.max_bd[axis]);
+        }
+    }
+}
+
+pub struct Triangle {
+    pub v1: [f32; 3],
+    pub v2: [f32; 3],
+    pub v3: [f32; 3],
+}

src/stl/parser.rs

@@ -4,15 +4,17 @@ use std::collections::HashMap;
 use std::fs::File;
 use std::io::{BufRead, BufReader, Read};
 use std::path::Path;
+use vecmath::vec3_add;
 
 /// Representation of an STL mesh
 /// Vertices are stored in vec with repeat vertices removed. Each vertex is represented as an
 /// [f32; 3]. Triangles are represented as a [usize; 3] where each entry is an index in the
 /// vertices Vec
 pub struct StlSolid {
-    vertices: Vec<[f32; 3]>,
+    pub vertices: Vec<[f32; 3]>,
-    triangles: Vec<[usize; 3]>,
+    pub triangles: Vec<[usize; 3]>,
-    normals: Vec<[f32; 3]>,
+    pub triangle_centroids: Vec<[f32; 3]>,
+    pub normals: Vec<[f32; 3]>,
 }
 
 pub fn parse_ascii_stl(file_path: &impl AsRef<Path>) -> Result<StlSolid> {
@@ -110,9 +112,11 @@ pub fn parse_ascii_stl(file_path: &impl AsRef<Path>) -> Result<StlSolid> {
         ];
     }
 
+    let triangle_centroids = calc_triangle_centroids(&vertices, &triangles);
     Ok(StlSolid {
         vertices,
         triangles,
+        triangle_centroids,
         normals,
     })
 }
@@ -183,13 +187,26 @@ pub fn parse_binary_stl(file_path: &impl AsRef<Path>) -> Result<StlSolid> {
         triangles.push(indices);
     }
 
+    let triangle_centroids = calc_triangle_centroids(&vertices, &triangles);
     Ok(StlSolid {
         vertices,
         triangles,
+        triangle_centroids,
         normals,
     })
 }
 
+/// Calculate the centroids for each triangle in a list
+fn calc_triangle_centroids(vertices: &[[f32; 3]], triangles: &[[usize; 3]]) -> Vec<[f32; 3]> {
+    triangles
+        .iter()
+        .map(|tri| {
+            tri.iter()
+                .fold([0.0; 3], |acc, idx| vec3_add(vertices[*idx], acc))
+        })
+        .collect()
+}
+
 #[cfg(test)]
 pub mod test {
     use super::*;

src/utilities/intersection_checks.rs

@@ -0,0 +1,48 @@
+use vecmath::{vec3_cross, vec3_dot, vec3_sub};
+
+use crate::primitives::{Box3, Triangle};
+
+/// Determine whether two shapes interesect using the separating axis theorem
+pub fn tri_aabb_intersect(bx: &Box3, tri: &Triangle) -> bool {
+    // Get edge vectors for triangle
+    let e1 = vec3_sub(tri.v2, tri.v1);
+    let e2 = vec3_sub(tri.v3, tri.v2);
+    let e3 = vec3_sub(tri.v1, tri.v3);
+
+    // Get face and edge normals
+    let tri_face = vec3_cross(e1, e2);
+    let n1 = vec3_cross(e1, e1);
+    let n2 = vec3_cross(e2, e1);
+    let n3 = vec3_cross(e3, e1);
+
+    // Face normals AABB
+    let u0 = [1.0, 0.0, 0.0];
+    let u1 = [0.0, 1.0, 0.0];
+    let u2 = [0.0, 0.0, 1.0];
+
+    // Now check all of the axes
+    let box_vertices = bx.get_vertices();
+    for axis in [u0, u1, u2, tri_face, n1, n2, n3] {
+        let box_projection = project_pnts(&box_vertices, axis);
+        let tri_projection = project_pnts(&[tri.v1, tri.v2, tri.v3], axis);
+
+        let start = f32::max(box_projection.0, tri_projection.0);
+        let end = f32::min(box_projection.1, tri_projection.1);
+        if start < end {
+            return false;
+        }
+    }
+    true
+}
+
+/// Project points unto an Axis and return the max and min
+fn project_pnts(pnts: &[[f32; 3]], axis: [f32; 3]) -> (f32, f32) {
+    let mut min_proj = f32::INFINITY;
+    let mut max_proj = f32::NEG_INFINITY;
+    for pnt in pnts {
+        let projection = vec3_dot(axis, *pnt);
+        min_proj = min_proj.min(projection);
+        max_proj = max_proj.max(projection);
+    }
+    (min_proj, max_proj)
+}

src/utilities/mod.rs

@@ -1,2 +1,3 @@
+pub mod intersection_checks;
 #[cfg(test)]
 pub mod test_macros;