# How to avoid global delaunay check in conforming triangulation?

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 ===
for triangle in mesh
.iter_faces()
{
let weight = triangle_weight(&triangle);
}

// 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 ===
{
// === 2. Pop current worst triangle ===
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())
{
{
continue;
}
let weight = triangle_weight(&face);
}
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())
{
{
continue;
}
let weight = triangle_weight(&face);
}
}
}

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.

I figured it out. The implementation of lawson's algorithm I had was incorrect, this is the right one:


/// 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.face().id()),
);

let find_opposite = |face: &FaceHandle<Vec3, (), ()>| {
face.hedge()
.iter_loop()
.find(|e| e.source().id() != source && e.dest().id() != source)
.unwrap()
};
while !queue.is_empty()
{
let face_id = queue.pop_front().unwrap();
if face_id == ABSENT()
{
continue;
}
let edge_handle = find_opposite(&mesh.face_handle(face_id));
if delaunay_test(&edge_handle)
&& edge_handle.can_flip()
&& !is_constraint_edge(&edge_handle)
{
edge_handle.flip();
queue.push_back(edge_handle.face().id());
queue.push_back(edge_handle.pair().face().id());
}
}
}
$$$$
`