I implemented a conforming (i.e. it creates Steiner points using Ruppert's algorithm) delaunay triangulator, which is working, but there is one step I am doing that I straight up don't understand and suggests something is broken.
After introducing a circumcenter of a triangle, I am having to check the entire mesh and restore the Delaunay property manually. I don't understand why I would need to do this?
I have tested all the subroutines the algorithm depends on, so I know the issue is coming from the logic of the algorithm itself.
This is my implementation in rust:
pub fn conforming_delaunay_triangulation_as_mesh(
points: &dyn Fn(usize) -> Vec3,
point_count: usize,
segments: &dyn Fn(usize) -> [usize; 2],
segment_count: usize,
angle_constraint: f32,
area_constraint: f32,
) -> HalfMesh<Vec3, (), ()>
{
// Compute a triangulation.
let mesh = constrained_delaunay_triangulation_hollow_as_mesh(
points,
point_count,
segments,
segment_count,
);
// Sort all *bad* faces from worst to best in apriority queue.
// === 1. Score all triangles and filter the bad ones ===
let mut bad_face_queue = PQueue::new();
for triangle in mesh
.iter_faces()
.filter(|t| triangle_is_bad(t, angle_constraint, area_constraint))
{
let weight = triangle_weight(&triangle);
bad_face_queue.push(triangle.id(), weight);
}
// Find all constrained and boundary segments.
let mut segment_set: BTreeSet<GeomId<1>> = (0..segment_count)
.map(|i| {
let [v1, v2] = segments(i);
mesh.vert_handle(GeomId::from(v1))
.shared_hedge(GeomId::from(v2))
.unwrap()
})
.map(|e| e.id().select_even())
.collect();
segment_set.extend(
mesh.iter_edges()
.filter(|e| e.is_in_boundary_edge())
.map(|e| e.id()),
);
// === Algorithm ===
while !bad_face_queue.is_empty()
{
// === 2. Pop current worst triangle ===
let (fid, _w) = bad_face_queue.pop().unwrap();
let points = mesh.face_handle(fid).vertices_tri();
let circumcenter = osculating_circle(&points[0], &points[1], &points[2]);
// === 3. Check if the circumcircle encroaches an edge ===
let mut segment_encroached = false;
let mut edge = mesh.hedge_handle(GeomId(0));
for segment in &segment_set
{
edge = mesh.hedge_handle(*segment);
if encroaches(&edge, &circumcenter)
{
segment_encroached = true;
break;
}
}
// === 4. If it encroaches the edge, split the edge ===
if segment_encroached
{
let new_edges;
let vid;
// Split the edge (check is needed because boundary edges must use special
// logic).
if !edge.is_in_boundary_edge()
{
(vid, _, new_edges) = edge.split::<f32>();
}
else
{
(vid, _, new_edges) = edge
.get_non_boundary_hedge()
.unwrap()
.split_boundary::<f32>();
}
// Restore Delaunay criterion.
lawsons_algorithm(vid, &mesh, |e| e.is_in_boundary_edge());
// We added new faces and changed existing ones, update the queue.
for face in mesh
.vert_handle(vid)
.iter_hedges()
.filter(|e| !e.is_boundary_hedge())
.map(|e| e.face())
{
if !triangle_is_bad(&face, angle_constraint, area_constraint)
{
continue;
}
let weight = triangle_weight(&face);
bad_face_queue.push(face.id(), weight);
}
segment_set.insert(new_edges[1]);
}
// === 5. If it does not encroach the edge, split the face ===
else
{
let face = find_containing_triangle(&circumcenter, &mesh);
if face.is_none()
{
continue;
}
let face = face.unwrap();
let vid = face.split_triangle_face::<f32>();
// Set face to circumcircle:
let v = mesh.vert_handle(vid);
v.mutate_data(&circumcenter);
// Restore Delaunay criterion.
lawsons_algorithm(vid, &mesh, |e| e.is_in_boundary_edge());
// TODO(low): why is this necessary?
for v in mesh.iter_verts()
{
lawsons_algorithm(v.id(), &mesh, |e| e.is_in_boundary_edge());
}
// We added new faces and changed existing ones, update the queue.
for face in mesh.vert_handle(vid).iter_hedges().map(|e| e.face())
{
if !triangle_is_bad(&face, angle_constraint, area_constraint)
{
continue;
}
let weight = triangle_weight(&face);
bad_face_queue.push(face.id(), weight);
}
}
}
mesh
}
/// Recursively flips triangles around a determined vertex until all pass the
/// delaunay test (or are constrained).
pub fn lawsons_algorithm<F>(
source: GeomId<0>,
mesh: &HalfMesh<Vec3, (), ()>,
is_constraint_edge: F,
) where
F: Fn(&HEdgeHandle<Vec3, (), ()>) -> bool,
{
let mut queue = VecDeque::from_iter(
mesh.vert_handle(source)
.iter_hedges()
.map(|e| e.next().id()),
);
while !queue.is_empty()
{
let edge_id = queue.pop_front().unwrap();
let edge_handle = mesh.hedge_handle(edge_id);
if delaunay_test(&edge_handle)
&& edge_handle.can_flip()
&& !is_constraint_edge(&edge_handle)
{
edge_handle.flip();
queue.push_back(edge_handle.next().id());
queue.push_back(edge_handle.pair().prev().id());
}
}
}
The summarizing pseudo code would be:
Compute a constrained Delaunay triangulation.
Create a priority queue of bad triangles from worst face to least bad
Create a set of segments
While there are bad triangles:
- Pop the worst triangle
- If its circumcenter encroaches a segment in the set, split that segment, add the new segment created by the split to the existing set, restore the Delaunay condition, update face queue.
- Else, find the face containing the circumcenter, split the face at the circumcenter, restore the Delaunay condition around the circumcenter.
- HACK!!!!! restore the Delaunay cirterion globally across the mesh
- Update the triangle queue with the faces surrounding the circumcenter.
I don't understand why I need to do the hack. Either my lawson algorithm is incorrect, or I have miss understood something key about the algorithm. I am hoping someone can lend me a hand.
