Commit 889495e0 authored by Philip Trettner's avatar Philip Trettner
Browse files

delaunay triangulation

parent 428ea3ad
......@@ -117,6 +117,11 @@ Scalar angle_defect(vertex_handle v, vertex_attribute<Vec3> const& position);
template <class Vec3, class Scalar = typename field_3d<Vec3>::Scalar>
vertex_attribute<Scalar> vertex_voronoi_areas(Mesh const& m, vertex_attribute<Vec3> const& position);
/// returns true if the edge satisfies the delaunay property
/// NOTE: only works on triangles
template <class Vec3>
bool is_delaunay(edge_handle e, vertex_attribute<Vec3> const& position);
/// ======== IMPLEMENTATION ========
inline bool is_boundary(vertex_handle v) { return v.is_boundary(); }
......@@ -359,4 +364,27 @@ vertex_attribute<Scalar> vertex_voronoi_areas(Mesh const& m, vertex_attribute<Ve
return areas;
}
template <class Vec3>
bool is_delaunay(edge_handle e, vertex_attribute<Vec3> const& position)
{
auto h0 = e.halfedgeA();
auto h1 = e.halfedgeB();
auto pi = position[h0.vertex_to()];
auto pj = position[h1.vertex_to()];
auto pa = position[h0.next().vertex_to()];
auto pb = position[h1.next().vertex_to()];
auto e_ia = pi - pa;
auto e_ja = pj - pa;
auto e_ib = pi - pb;
auto e_jb = pj - pb;
auto cot_a = field_3d<Vec3>::dot(e_ia, e_ja) / field_3d<Vec3>::length(field_3d<Vec3>::cross(e_ia, e_ja));
auto cot_b = field_3d<Vec3>::dot(e_ib, e_jb) / field_3d<Vec3>::length(field_3d<Vec3>::cross(e_ib, e_jb));
return cot_a + cot_b >= 0;
}
}
#pragma once
#include <polymesh/Mesh.hh>
#include <polymesh/algorithms/properties.hh>
namespace polymesh
{
/// Given a triangular mesh, performs edge flips until all flippable edges are delaunay
/// (extrinsic delaunay triangulation)
template <class Vec3>
void make_delaunay(Mesh& m, vertex_attribute<Vec3> const& position);
/// ======== IMPLEMENTATION ========
template <class Vec3>
void make_delaunay(Mesh& m, vertex_attribute<Vec3> const& position)
{
std::vector<edge_index> queue;
for (auto e : m.edges())
queue.push_back(e);
while (!queue.empty())
{
auto e = queue.back().of(m);
queue.pop_back();
if (e.is_boundary())
continue;
if (is_delaunay(e, position))
continue;
if (valence(e.vertexA()) <= 2 || valence(e.vertexB()) <= 2)
continue;
queue.push_back(e.halfedgeA().next().edge());
queue.push_back(e.halfedgeA().prev().edge());
queue.push_back(e.halfedgeB().next().edge());
queue.push_back(e.halfedgeB().prev().edge());
m.edges().rotate_next(e);
}
}
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment