262 lines
8.2 KiB
Rust
262 lines
8.2 KiB
Rust
//! 3D Delaunay tetrahedralization using Bowyer-Watson algorithm.
|
|
|
|
use super::predicates::{in_circumsphere, orient3d};
|
|
use crate::error::{MeshGenError, Result};
|
|
use crate::mesh::{TetrahedralMesh, Tetrahedron, Vertex};
|
|
use nalgebra::Point3;
|
|
use std::collections::{HashMap, HashSet};
|
|
|
|
/// Compute 3D Delaunay tetrahedralization of a point set.
|
|
///
|
|
/// Uses the incremental Bowyer-Watson algorithm.
|
|
///
|
|
/// # Arguments
|
|
/// * `points` - Points to triangulate
|
|
///
|
|
/// # Returns
|
|
/// A tetrahedral mesh
|
|
pub fn delaunay_3d(points: &[Point3<f64>]) -> Result<TetrahedralMesh> {
|
|
if points.len() < 4 {
|
|
return Err(MeshGenError::InvalidInput(
|
|
"At least 4 points required for tetrahedralization".to_string(),
|
|
));
|
|
}
|
|
|
|
// Compute bounding box
|
|
let (min_pt, max_pt) = bounding_box(points);
|
|
let span = max_pt - min_pt;
|
|
let max_span = span.x.max(span.y).max(span.z);
|
|
|
|
// Create super-tetrahedron that contains all points
|
|
// Make it large enough with some margin
|
|
let center = Point3::from((min_pt.coords + max_pt.coords) / 2.0);
|
|
let size = max_span * 10.0; // Large enough to contain all points
|
|
|
|
let super_tet = create_super_tetrahedron(¢er, size);
|
|
|
|
// Initialize triangulation with super-tetrahedron
|
|
let mut vertices: Vec<Point3<f64>> = super_tet.to_vec();
|
|
let mut tetrahedra: Vec<[usize; 4]> = vec![[0, 1, 2, 3]];
|
|
|
|
// Insert points one by one
|
|
for (pi, &point) in points.iter().enumerate() {
|
|
let vi = vertices.len();
|
|
vertices.push(point);
|
|
|
|
// Find all tetrahedra whose circumsphere contains the new point
|
|
let mut bad_tets: HashSet<usize> = HashSet::new();
|
|
for (ti, tet) in tetrahedra.iter().enumerate() {
|
|
let a = &vertices[tet[0]];
|
|
let b = &vertices[tet[1]];
|
|
let c = &vertices[tet[2]];
|
|
let d = &vertices[tet[3]];
|
|
|
|
// Check orientation and adjust if needed
|
|
let orient = orient3d(a, b, c, d);
|
|
if orient.abs() < 1e-15 {
|
|
continue; // Degenerate, skip
|
|
}
|
|
|
|
let test = if orient > 0.0 {
|
|
in_circumsphere(a, b, c, d, &point)
|
|
} else {
|
|
in_circumsphere(a, b, d, c, &point)
|
|
};
|
|
|
|
if test > 0.0 {
|
|
bad_tets.insert(ti);
|
|
}
|
|
}
|
|
|
|
if bad_tets.is_empty() {
|
|
// Point might be outside all circumspheres (degenerate case)
|
|
// Skip this point
|
|
vertices.pop();
|
|
continue;
|
|
}
|
|
|
|
// Find the boundary of the cavity (faces not shared by bad tets)
|
|
let boundary_faces = find_cavity_boundary(&tetrahedra, &bad_tets);
|
|
|
|
// Remove bad tetrahedra (mark for removal, will compact later)
|
|
let mut new_tets: Vec<[usize; 4]> = Vec::new();
|
|
for (ti, tet) in tetrahedra.iter().enumerate() {
|
|
if !bad_tets.contains(&ti) {
|
|
new_tets.push(*tet);
|
|
}
|
|
}
|
|
|
|
// Create new tetrahedra by connecting boundary faces to new point
|
|
for face in boundary_faces {
|
|
// Ensure correct orientation
|
|
let a = &vertices[face[0]];
|
|
let b = &vertices[face[1]];
|
|
let c = &vertices[face[2]];
|
|
let p = &vertices[vi];
|
|
|
|
let orient = orient3d(a, b, c, p);
|
|
let new_tet = if orient > 0.0 {
|
|
[face[0], face[1], face[2], vi]
|
|
} else {
|
|
[face[0], face[2], face[1], vi]
|
|
};
|
|
|
|
new_tets.push(new_tet);
|
|
}
|
|
|
|
tetrahedra = new_tets;
|
|
|
|
// Progress reporting for large point sets
|
|
if points.len() > 1000 && pi % 1000 == 0 {
|
|
// Could add callback here for progress
|
|
}
|
|
}
|
|
|
|
// Remove tetrahedra connected to super-tetrahedron vertices (indices 0-3)
|
|
let super_indices: HashSet<usize> = [0, 1, 2, 3].into_iter().collect();
|
|
tetrahedra.retain(|tet| !tet.iter().any(|&vi| super_indices.contains(&vi)));
|
|
|
|
// Compact vertices (remove super-tet vertices and reindex)
|
|
let mut mesh = TetrahedralMesh::new();
|
|
let mut old_to_new: HashMap<usize, usize> = HashMap::new();
|
|
|
|
for (old_idx, &point) in vertices.iter().enumerate().skip(4) {
|
|
let new_idx = mesh.add_vertex(Vertex::new(point.x, point.y, point.z));
|
|
old_to_new.insert(old_idx, new_idx);
|
|
}
|
|
|
|
for tet in tetrahedra {
|
|
if let (Some(&v0), Some(&v1), Some(&v2), Some(&v3)) = (
|
|
old_to_new.get(&tet[0]),
|
|
old_to_new.get(&tet[1]),
|
|
old_to_new.get(&tet[2]),
|
|
old_to_new.get(&tet[3]),
|
|
) {
|
|
mesh.add_tetrahedron(Tetrahedron::new(v0, v1, v2, v3));
|
|
}
|
|
}
|
|
|
|
Ok(mesh)
|
|
}
|
|
|
|
/// Create a super-tetrahedron that contains all points.
|
|
fn create_super_tetrahedron(center: &Point3<f64>, size: f64) -> [Point3<f64>; 4] {
|
|
// Create a large regular tetrahedron centered at `center`
|
|
let h = size * 2.0; // Height from center to vertex
|
|
|
|
[
|
|
Point3::new(center.x, center.y + h, center.z),
|
|
Point3::new(center.x - h * 0.866, center.y - h * 0.5, center.z - h * 0.5),
|
|
Point3::new(center.x + h * 0.866, center.y - h * 0.5, center.z - h * 0.5),
|
|
Point3::new(center.x, center.y - h * 0.5, center.z + h),
|
|
]
|
|
}
|
|
|
|
/// Find boundary faces of the cavity formed by bad tetrahedra.
|
|
fn find_cavity_boundary(tetrahedra: &[[usize; 4]], bad_tets: &HashSet<usize>) -> Vec<[usize; 3]> {
|
|
let mut face_count: HashMap<[usize; 3], usize> = HashMap::new();
|
|
let mut face_original: HashMap<[usize; 3], [usize; 3]> = HashMap::new();
|
|
|
|
for &ti in bad_tets {
|
|
let tet = tetrahedra[ti];
|
|
let faces = [
|
|
[tet[1], tet[2], tet[3]],
|
|
[tet[0], tet[3], tet[2]],
|
|
[tet[0], tet[1], tet[3]],
|
|
[tet[0], tet[2], tet[1]],
|
|
];
|
|
|
|
for face in faces {
|
|
let mut sorted = face;
|
|
sorted.sort_unstable();
|
|
|
|
*face_count.entry(sorted).or_insert(0) += 1;
|
|
face_original.entry(sorted).or_insert(face);
|
|
}
|
|
}
|
|
|
|
// Boundary faces appear exactly once
|
|
face_count
|
|
.into_iter()
|
|
.filter(|(_, count)| *count == 1)
|
|
.map(|(sorted, _)| *face_original.get(&sorted).unwrap())
|
|
.collect()
|
|
}
|
|
|
|
/// Compute bounding box of points.
|
|
fn bounding_box(points: &[Point3<f64>]) -> (Point3<f64>, Point3<f64>) {
|
|
let mut min_pt = points[0];
|
|
let mut max_pt = points[0];
|
|
|
|
for p in points {
|
|
min_pt.x = min_pt.x.min(p.x);
|
|
min_pt.y = min_pt.y.min(p.y);
|
|
min_pt.z = min_pt.z.min(p.z);
|
|
max_pt.x = max_pt.x.max(p.x);
|
|
max_pt.y = max_pt.y.max(p.y);
|
|
max_pt.z = max_pt.z.max(p.z);
|
|
}
|
|
|
|
(min_pt, max_pt)
|
|
}
|
|
|
|
#[cfg(test)]
|
|
mod tests {
|
|
use super::*;
|
|
|
|
#[test]
|
|
fn test_delaunay_simple() {
|
|
// Simple cube vertices
|
|
let points = vec![
|
|
Point3::new(0.0, 0.0, 0.0),
|
|
Point3::new(1.0, 0.0, 0.0),
|
|
Point3::new(0.0, 1.0, 0.0),
|
|
Point3::new(1.0, 1.0, 0.0),
|
|
Point3::new(0.0, 0.0, 1.0),
|
|
Point3::new(1.0, 0.0, 1.0),
|
|
Point3::new(0.0, 1.0, 1.0),
|
|
Point3::new(1.0, 1.0, 1.0),
|
|
];
|
|
|
|
let mesh = delaunay_3d(&points).unwrap();
|
|
|
|
assert_eq!(mesh.num_vertices(), 8);
|
|
// A cube should be divided into 5 or 6 tetrahedra
|
|
assert!(mesh.num_tetrahedra() >= 5);
|
|
assert!(mesh.num_tetrahedra() <= 6);
|
|
}
|
|
|
|
#[test]
|
|
fn test_delaunay_random() {
|
|
use rand::Rng;
|
|
let mut rng = rand::thread_rng();
|
|
|
|
let points: Vec<Point3<f64>> = (0..20)
|
|
.map(|_| Point3::new(rng.r#gen::<f64>(), rng.r#gen::<f64>(), rng.r#gen::<f64>()))
|
|
.collect();
|
|
|
|
let mesh = delaunay_3d(&points).unwrap();
|
|
|
|
assert_eq!(mesh.num_vertices(), 20);
|
|
assert!(mesh.num_tetrahedra() > 0);
|
|
|
|
// All tetrahedra should have positive volume
|
|
for i in 0..mesh.num_tetrahedra() {
|
|
let vol = mesh.tetrahedron_volume(i);
|
|
assert!(vol > 0.0, "Tetrahedron {} has non-positive volume", i);
|
|
}
|
|
}
|
|
|
|
#[test]
|
|
fn test_too_few_points() {
|
|
let points = vec![
|
|
Point3::new(0.0, 0.0, 0.0),
|
|
Point3::new(1.0, 0.0, 0.0),
|
|
Point3::new(0.0, 1.0, 0.0),
|
|
];
|
|
|
|
let result = delaunay_3d(&points);
|
|
assert!(result.is_err());
|
|
}
|
|
}
|