292 lines
10 KiB
C#
292 lines
10 KiB
C#
// -----------------------------------------------------------------------
|
|
// <copyright file="VoronoiBase.cs">
|
|
// Original Triangle code by Jonathan Richard Shewchuk, http://www.cs.cmu.edu/~quake/triangle.html
|
|
// Triangle.NET code by Christian Woltering, http://triangle.codeplex.com/
|
|
// </copyright>
|
|
// -----------------------------------------------------------------------
|
|
|
|
namespace TriangleNet.Voronoi
|
|
{
|
|
using System.Collections.Generic;
|
|
|
|
using TriangleNet.Topology;
|
|
using TriangleNet.Geometry;
|
|
using TriangleNet.Topology.DCEL;
|
|
|
|
using Vertex = TriangleNet.Topology.DCEL.Vertex;
|
|
|
|
/// <summary>
|
|
/// The Voronoi diagram is the dual of a pointset triangulation.
|
|
/// </summary>
|
|
public abstract class VoronoiBase : DcelMesh
|
|
{
|
|
protected IPredicates predicates;
|
|
|
|
protected IVoronoiFactory factory;
|
|
|
|
// List of infinite half-edges, i.e. half-edges that start at circumcenters of triangles
|
|
// which lie on the domain boundary.
|
|
protected List<HalfEdge> rays;
|
|
|
|
/// <summary>
|
|
/// Initializes a new instance of the <see cref="VoronoiBase" /> class.
|
|
/// </summary>
|
|
/// <param name="mesh">Triangle mesh.</param>
|
|
/// <param name="factory">Voronoi object factory.</param>
|
|
/// <param name="predicates">Geometric predicates implementation.</param>
|
|
/// <param name="generate">If set to true, the constuctor will call the Generate
|
|
/// method, which builds the Voronoi diagram.</param>
|
|
protected VoronoiBase(Mesh mesh, IVoronoiFactory factory, IPredicates predicates,
|
|
bool generate)
|
|
: base(false)
|
|
{
|
|
this.factory = factory;
|
|
this.predicates = predicates;
|
|
|
|
if (generate)
|
|
{
|
|
Generate(mesh);
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Generate the Voronoi diagram from given triangle mesh..
|
|
/// </summary>
|
|
/// <param name="mesh"></param>
|
|
/// <param name="bounded"></param>
|
|
protected void Generate(Mesh mesh)
|
|
{
|
|
mesh.Renumber();
|
|
|
|
base.edges = new List<HalfEdge>();
|
|
this.rays = new List<HalfEdge>();
|
|
|
|
// Allocate space for Voronoi diagram.
|
|
var vertices = new Vertex[mesh.triangles.Count + mesh.hullsize];
|
|
var faces = new Face[mesh.vertices.Count];
|
|
|
|
if (factory == null)
|
|
{
|
|
factory = new DefaultVoronoiFactory();
|
|
}
|
|
|
|
factory.Initialize(vertices.Length, 2 * mesh.NumberOfEdges, faces.Length);
|
|
|
|
// Compute triangles circumcenters.
|
|
var map = ComputeVertices(mesh, vertices);
|
|
|
|
// Create all Voronoi faces.
|
|
foreach (var vertex in mesh.vertices.Values)
|
|
{
|
|
faces[vertex.id] = factory.CreateFace(vertex);
|
|
}
|
|
|
|
ComputeEdges(mesh, vertices, faces, map);
|
|
|
|
// At this point all edges are computed, but the (edge.next) pointers aren't set.
|
|
ConnectEdges(map);
|
|
|
|
base.vertices = new List<Vertex>(vertices);
|
|
base.faces = new List<Face>(faces);
|
|
}
|
|
|
|
/// <summary>
|
|
/// Compute the Voronoi vertices (the circumcenters of the triangles).
|
|
/// </summary>
|
|
/// <returns>An empty map, which will map all vertices to a list of leaving edges.</returns>
|
|
protected List<HalfEdge>[] ComputeVertices(Mesh mesh, Vertex[] vertices)
|
|
{
|
|
Otri tri = default(Otri);
|
|
double xi = 0, eta = 0;
|
|
Vertex vertex;
|
|
Point pt;
|
|
int id;
|
|
|
|
// Maps all vertices to a list of leaving edges.
|
|
var map = new List<HalfEdge>[mesh.triangles.Count];
|
|
|
|
// Compue triangle circumcenters
|
|
foreach (var t in mesh.triangles)
|
|
{
|
|
id = t.id;
|
|
tri.tri = t;
|
|
|
|
pt = predicates.FindCircumcenter(tri.Org(), tri.Dest(), tri.Apex(), ref xi, ref eta);
|
|
|
|
vertex = factory.CreateVertex(pt.x, pt.y);
|
|
vertex.id = id;
|
|
|
|
vertices[id] = vertex;
|
|
map[id] = new List<HalfEdge>();
|
|
}
|
|
|
|
return map;
|
|
}
|
|
|
|
/// <summary>
|
|
/// Compute the edges of the Voronoi diagram.
|
|
/// </summary>
|
|
/// <param name="mesh"></param>
|
|
/// <param name="vertices"></param>
|
|
/// <param name="faces"></param>
|
|
/// <param name="map">Empty vertex map.</param>
|
|
protected void ComputeEdges(Mesh mesh, Vertex[] vertices, Face[] faces, List<HalfEdge>[] map)
|
|
{
|
|
Otri tri, neighbor = default(Otri);
|
|
TriangleNet.Geometry.Vertex org, dest;
|
|
|
|
double px, py;
|
|
int id, nid, count = mesh.triangles.Count;
|
|
|
|
Face face, neighborFace;
|
|
HalfEdge edge, twin;
|
|
Vertex vertex, end;
|
|
|
|
// Count infinte edges (vertex id for their endpoints).
|
|
int j = 0;
|
|
|
|
// Count half-edges (edge ids).
|
|
int k = 0;
|
|
|
|
// To loop over the set of edges, loop over all triangles, and look at the
|
|
// three edges of each triangle. If there isn't another triangle adjacent
|
|
// to the edge, operate on the edge. If there is another adjacent triangle,
|
|
// operate on the edge only if the current triangle has a smaller id than
|
|
// its neighbor. This way, each edge is considered only once.
|
|
foreach (var t in mesh.triangles)
|
|
{
|
|
id = t.id;
|
|
|
|
tri.tri = t;
|
|
|
|
for (int i = 0; i < 3; i++)
|
|
{
|
|
tri.orient = i;
|
|
tri.Sym(ref neighbor);
|
|
|
|
nid = neighbor.tri.id;
|
|
|
|
if (id < nid || nid < 0)
|
|
{
|
|
// Get the endpoints of the current triangle edge.
|
|
org = tri.Org();
|
|
dest = tri.Dest();
|
|
|
|
face = faces[org.id];
|
|
neighborFace = faces[dest.id];
|
|
|
|
vertex = vertices[id];
|
|
|
|
// For each edge in the triangle mesh, there's a corresponding edge
|
|
// in the Voronoi diagram, i.e. two half-edges will be created.
|
|
if (nid < 0)
|
|
{
|
|
// Unbounded edge, direction perpendicular to the boundary edge,
|
|
// pointing outwards.
|
|
px = dest.y - org.y;
|
|
py = org.x - dest.x;
|
|
|
|
end = factory.CreateVertex(vertex.x + px, vertex.y + py);
|
|
end.id = count + j++;
|
|
|
|
vertices[end.id] = end;
|
|
|
|
edge = factory.CreateHalfEdge(end, face);
|
|
twin = factory.CreateHalfEdge(vertex, neighborFace);
|
|
|
|
// Make (face.edge) always point to an edge that starts at an infinite
|
|
// vertex. This will allow traversing of unbounded faces.
|
|
face.edge = edge;
|
|
face.bounded = false;
|
|
|
|
map[id].Add(twin);
|
|
|
|
rays.Add(twin);
|
|
}
|
|
else
|
|
{
|
|
end = vertices[nid];
|
|
|
|
// Create half-edges.
|
|
edge = factory.CreateHalfEdge(end, face);
|
|
twin = factory.CreateHalfEdge(vertex, neighborFace);
|
|
|
|
// Add to vertex map.
|
|
map[nid].Add(edge);
|
|
map[id].Add(twin);
|
|
}
|
|
|
|
vertex.leaving = twin;
|
|
end.leaving = edge;
|
|
|
|
edge.twin = twin;
|
|
twin.twin = edge;
|
|
|
|
edge.id = k++;
|
|
twin.id = k++;
|
|
|
|
this.edges.Add(edge);
|
|
this.edges.Add(twin);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Connect all edges of the Voronoi diagram.
|
|
/// </summary>
|
|
/// <param name="map">Maps all vertices to a list of leaving edges.</param>
|
|
protected virtual void ConnectEdges(List<HalfEdge>[] map)
|
|
{
|
|
int length = map.Length;
|
|
|
|
// For each half-edge, find its successor in the connected face.
|
|
foreach (var edge in this.edges)
|
|
{
|
|
var face = edge.face.generator.id;
|
|
|
|
// The id of the dest vertex of current edge.
|
|
int id = edge.twin.origin.id;
|
|
|
|
// The edge origin can also be an infinite vertex. Sort them out
|
|
// by checking the id.
|
|
if (id < length)
|
|
{
|
|
// Look for the edge that is connected to the current face. Each
|
|
// Voronoi vertex has degree 3, so this loop is actually O(1).
|
|
foreach (var next in map[id])
|
|
{
|
|
if (next.face.generator.id == face)
|
|
{
|
|
edge.next = next;
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
protected override IEnumerable<IEdge> EnumerateEdges()
|
|
{
|
|
var edges = new List<IEdge>(this.edges.Count / 2);
|
|
|
|
foreach (var edge in this.edges)
|
|
{
|
|
var twin = edge.twin;
|
|
|
|
// Report edge only once.
|
|
if (twin == null)
|
|
{
|
|
edges.Add(new Edge(edge.origin.id, edge.next.origin.id));
|
|
}
|
|
else if (edge.id < twin.id)
|
|
{
|
|
edges.Add(new Edge(edge.origin.id, twin.origin.id));
|
|
}
|
|
}
|
|
|
|
return edges;
|
|
}
|
|
}
|
|
}
|