Improved uniform data structure for Voronoi classes.

git-svn-id: https://triangle.svn.codeplex.com/svn@69976 0e2699bc-83d4-4a8f-98e7-55e24ab8c7a5
This commit is contained in:
SND\wo80_cp
2012-10-14 22:08:33 +00:00
parent 894c3e6eba
commit 04c45ed41a
9 changed files with 610 additions and 198 deletions
+49 -2
View File
@@ -8,9 +8,7 @@ namespace MeshRenderer.Core
{
using System.Collections.Generic;
using System.Linq;
using System.Drawing;
using TriangleNet;
using TriangleNet.Data;
using TriangleNet.Geometry;
using TriangleNet.Tools;
@@ -156,6 +154,55 @@ namespace MeshRenderer.Core
/// </summary>
public void SetVoronoi(IVoronoi voro, int infCount)
{
int i, n = voro.Points.Length;
// Copy points
this.VoronoiPoints = new float[2 * n + infCount];
foreach (var v in voro.Points)
{
if (v == null)
{
continue;
}
i = v.ID;
this.VoronoiPoints[2 * i] = (float)v.X;
this.VoronoiPoints[2 * i + 1] = (float)v.Y;
}
// Copy edges
Point first, last;
var edges = new List<uint>(voro.Regions.Count * 4);
foreach (var region in voro.Regions)
{
first = null;
last = null;
foreach (var pt in region.Vertices)
{
if (first == null)
{
first = pt;
last = pt;
}
else
{
edges.Add((uint)last.ID);
edges.Add((uint)pt.ID);
last = pt;
}
}
if (region.Bounded && first != null)
{
edges.Add((uint)last.ID);
edges.Add((uint)first.ID);
}
}
this.VoronoiEdges = edges.ToArray();
/*
int i, n = voro.VertexList.Count;
+24 -2
View File
@@ -234,7 +234,6 @@ namespace MeshExplorer
// Clear voronoi
menuViewVoronoi.Checked = false;
//renderManager.ShowVoronoi = false;
// Disable menu items
menuFileSave.Enabled = false;
@@ -290,6 +289,9 @@ namespace MeshExplorer
// Update Statistic view
statisticView.HandleMeshChange(mesh);
// TODO: Should the Voronoi diagram automaticaly update?
menuViewVoronoi.Checked = false;
// Enable menu items
menuFileSave.Enabled = true;
menuFileExport.Enabled = true;
@@ -632,8 +634,28 @@ namespace MeshExplorer
private void menuViewVoronoi_Click(object sender, EventArgs e)
{
if (mesh == null)
{
return;
}
IVoronoi voronoi;
if (mesh.IsPolygon)
{
voronoi = new BoundedVoronoi(mesh);
}
else
{
voronoi = new Voronoi(mesh);
}
voronoi.Generate();
renderData.SetVoronoi(voronoi);
renderManager.SetData(renderData);
menuViewVoronoi.Checked = !menuViewVoronoi.Checked;
//renderControl1.ShowVoronoi = menuViewVoronoi.Checked;
}
private void menuViewLog_Click(object sender, EventArgs e)
-11
View File
@@ -21,9 +21,6 @@ namespace TriangleNet.Data
// Hash for dictionary. Will be set by mesh instance.
internal int hash;
// The ID is only used for mesh output.
internal int id;
internal VertexType type;
internal Otri tri;
@@ -76,14 +73,6 @@ namespace TriangleNet.Data
#region Public properties
/// <summary>
/// Gets the vertex id.
/// </summary>
public int ID
{
get { return this.id; }
}
/// <summary>
/// Gets the vertex type.
/// </summary>
+9
View File
@@ -16,6 +16,7 @@ namespace TriangleNet.Geometry
/// </summary>
public class Point : IComparable<Point>, IEquatable<Point>
{
internal int id;
internal double x;
internal double y;
internal int mark;
@@ -40,6 +41,14 @@ namespace TriangleNet.Geometry
#region Public properties
/// <summary>
/// Gets the vertex id.
/// </summary>
public int ID
{
get { return this.id; }
}
/// <summary>
/// Gets the vertex x coordinate.
/// </summary>
+170 -97
View File
@@ -20,11 +20,16 @@ namespace TriangleNet.Tools
/// 2D Centroidal Voronoi Tessellations with Constraints, 2010,
/// Jane Tournois, Pierre Alliez and Olivier Devillers
/// </remarks>
public class BoundedVoronoi
public class BoundedVoronoi : IVoronoi
{
Mesh mesh;
Point[] points;
List<VoronoiRegion> regions;
int segIndex;
Dictionary<int, Segment> subsegMap;
List<Point[]> voronoi;
/// <summary>
/// Initializes a new instance of the <see cref="BoundedVoronoi" /> class.
@@ -33,19 +38,22 @@ namespace TriangleNet.Tools
public BoundedVoronoi(Mesh mesh)
{
this.mesh = mesh;
this.voronoi = new List<Point[]>();
// TODO: introduce isVertexMapValid in mesh and don't call
// MakeVertexMap() if not necessary.
this.mesh.MakeVertexMap();
}
/// <summary>
/// Gets the list of Voronoi cells.
/// Gets the list of Voronoi vertices.
/// </summary>
public List<Point[]> Cells
public Point[] Points
{
get { return voronoi; }
get { return points; }
}
/// <summary>
/// Gets the list of Voronoi regions.
/// </summary>
public List<VoronoiRegion> Regions
{
get { return regions; }
}
/// <summary>
@@ -53,7 +61,14 @@ namespace TriangleNet.Tools
/// </summary>
public void Generate()
{
voronoi.Clear();
mesh.Renumber();
mesh.MakeVertexMap();
// Allocate space for voronoi diagram
this.points = new Point[mesh.triangles.Count + mesh.subsegs.Count * 5]; // This is an upper bound.
this.regions = new List<VoronoiRegion>(mesh.vertices.Count);
ComputeCircumCenters();
TagBlindTriangles();
@@ -62,15 +77,33 @@ namespace TriangleNet.Tools
// TODO: Need a reliable way to check if a vertex is on a segment
if (v.type == VertexType.FreeVertex || v.Boundary == 0)
{
voronoi.Add(ConstructBvdCell(v));
ConstructBvdCell(v);
}
else
{
voronoi.Add(ConstructBoundaryBvdCell(v));
ConstructBoundaryBvdCell(v);
}
}
}
private void ComputeCircumCenters()
{
Otri tri = default(Otri);
double xi = 0, eta = 0;
Point pt;
// Compue triangle circumcenters
foreach (var item in mesh.triangles.Values)
{
tri.triangle = item;
pt = Primitives.FindCircumcenter(tri.Org(), tri.Dest(), tri.Apex(), ref xi, ref eta);
pt.id = item.id;
points[item.id] = pt;
}
}
/// <summary>
/// Tag all blind triangles.
/// </summary>
@@ -169,26 +202,28 @@ namespace TriangleNet.Tools
/// <returns>Returns true, if the triangle is blinded.</returns>
private bool TriangleIsBlinded(ref Otri tri, ref Osub seg)
{
Point cc, pt = new Point();
double xi = 0, eta = 0;
Point c, pt;
Vertex torg = tri.Org();
Vertex tdest = tri.Dest();
Vertex tapex = tri.Apex();
cc = Primitives.FindCircumcenter(torg, tdest, tapex, ref xi, ref eta);
Vertex sorg = seg.Org();
Vertex sdest = seg.Dest();
if (SegmentsIntersect(ref seg, cc, torg, ref pt, true))
c = this.points[tri.triangle.id];
if (SegmentsIntersect(sorg, sdest, c, torg, out pt, true))
{
return true;
}
if (SegmentsIntersect(ref seg, cc, tdest, ref pt, true))
if (SegmentsIntersect(sorg, sdest, c, tdest, out pt, true))
{
return true;
}
if (SegmentsIntersect(ref seg, cc, tapex, ref pt, true))
if (SegmentsIntersect(sorg, sdest, c, tapex, out pt, true))
{
return true;
}
@@ -196,25 +231,28 @@ namespace TriangleNet.Tools
return false;
}
private Point[] ConstructBvdCell(Vertex x)
private void ConstructBvdCell(Vertex vertex)
{
VoronoiRegion region = new VoronoiRegion(vertex);
regions.Add(region);
Otri f = default(Otri);
Otri f_init = default(Otri);
Otri f_next = default(Otri);
Osub sf = default(Osub);
Osub sfn = default(Osub);
Vertex torg, tdest, tapex;
Point cc_f, cc_f_next, p;
double xi = 0, eta = 0;
int n = mesh.triangles.Count;
// Call P the polygon (cell) in construction
List<Point> P = new List<Point>();
List<Point> vpoints = new List<Point>();
// Call f_init a triangle incident to x
x.tri.Copy(ref f_init);
vertex.tri.Copy(ref f_init);
if (f_init.Org() != x)
if (f_init.Org() != vertex)
{
throw new Exception("ConstructBvdCell: inconsistent topology.");
}
@@ -228,32 +266,28 @@ namespace TriangleNet.Tools
do
{
// Call Lffnext the line going through the circumcenters of f and f_next
torg = f.Org();
tdest = f.Dest();
tapex = f.Apex();
cc_f = Primitives.FindCircumcenter(torg, tdest, tapex, ref xi, ref eta);
torg = f_next.Org();
tdest = f_next.Dest();
tapex = f_next.Apex();
cc_f_next = Primitives.FindCircumcenter(torg, tdest, tapex, ref xi, ref eta);
cc_f = this.points[f.triangle.id];
cc_f_next = this.points[f_next.triangle.id];
// if f is tagged non-blind then
if (!f.triangle.infected)
{
// Insert the circumcenter of f into P
P.Add(cc_f);
vpoints.Add(cc_f);
if (f_next.triangle.infected)
{
// Call S_fnext the constrained edge blinding f_next
sfn.seg = subsegMap[f_next.triangle.hash];
p = new Point();
// Insert point Lf,f_next /\ Sf_next into P
if (SegmentsIntersect(ref sfn, cc_f, cc_f_next, ref p, true))
if (SegmentsIntersect(sfn.SegOrg(), sfn.SegDest(), cc_f, cc_f_next, out p, true))
{
P.Add(p);
p.id = n + segIndex;
points[n + segIndex] = p;
segIndex++;
vpoints.Add(p);
}
}
}
@@ -265,11 +299,14 @@ namespace TriangleNet.Tools
// if f_next is tagged non-blind then
if (!f_next.triangle.infected)
{
p = new Point();
// Insert point Lf,f_next /\ Sf into P
if (SegmentsIntersect(ref sf, cc_f, cc_f_next, ref p, true))
if (SegmentsIntersect(sf.SegOrg(), sf.SegDest(), cc_f, cc_f_next, out p, true))
{
P.Add(p);
p.id = n + segIndex;
points[n + segIndex] = p;
segIndex++;
vpoints.Add(p);
}
}
else
@@ -280,17 +317,23 @@ namespace TriangleNet.Tools
// if Sf != Sf_next then
if (!sf.Equal(sfn))
{
p = new Point();
// Insert Lf,fnext /\ Sf and Lf,fnext /\ Sfnext into P
if (SegmentsIntersect(ref sf, cc_f, cc_f_next, ref p, true))
if (SegmentsIntersect(sf.SegOrg(), sf.SegDest(), cc_f, cc_f_next, out p, true))
{
P.Add(p);
p.id = n + segIndex;
points[n + segIndex] = p;
segIndex++;
vpoints.Add(p);
}
p = new Point();
if (SegmentsIntersect(ref sfn, cc_f, cc_f_next, ref p, true))
if (SegmentsIntersect(sfn.SegOrg(), sfn.SegDest(), cc_f, cc_f_next, out p, true))
{
P.Add(p);
p.id = n + segIndex;
points[n + segIndex] = p;
segIndex++;
vpoints.Add(p);
}
}
}
@@ -305,11 +348,14 @@ namespace TriangleNet.Tools
while (!f.Equal(f_init));
// Output: Bounded Voronoi cell of x in counterclockwise order.
return P.ToArray();
region.Add(vpoints);
}
private Point[] ConstructBoundaryBvdCell(Vertex x)
private void ConstructBoundaryBvdCell(Vertex vertex)
{
VoronoiRegion region = new VoronoiRegion(vertex);
regions.Add(region);
Otri f = default(Otri);
Otri f_init = default(Otri);
Otri f_next = default(Otri);
@@ -317,17 +363,18 @@ namespace TriangleNet.Tools
Osub sf = default(Osub);
Osub sfn = default(Osub);
Vertex torg, tdest, tapex;
Vertex torg, tdest, tapex, sorg, sdest;
Point cc_f, cc_f_next, p;
double xi = 0, eta = 0;
int n = mesh.triangles.Count;
// Call P the polygon (cell) in construction
List<Point> P = new List<Point>();
List<Point> vpoints = new List<Point>();
// Call f_init a triangle incident to x
x.tri.Copy(ref f_init);
vertex.tri.Copy(ref f_init);
if (f_init.Org() != x)
if (f_init.Org() != vertex)
{
throw new Exception("ConstructBoundaryBvdCell: inconsistent topology.");
}
@@ -353,60 +400,72 @@ namespace TriangleNet.Tools
}
// Add vertex on border
P.Add(x);
p = new Point(vertex.x, vertex.y);
p.id = n + segIndex;
points[n + segIndex] = p;
segIndex++;
vpoints.Add(p);
// Add midpoint of start triangles' edge.
torg = f.Org();
tdest = f.Dest();
P.Add(new Point((torg.X + tdest.X) / 2, (torg.Y + tdest.Y) / 2));
p = new Point((torg.X + tdest.X) / 2, (torg.Y + tdest.Y) / 2);
p.id = n + segIndex;
points[n + segIndex] = p;
segIndex++;
vpoints.Add(p);
// repeat ... until f = f_init
do
{
// Call Lffnext the line going through the circumcenters of f and f_next
torg = f.Org();
tdest = f.Dest();
tapex = f.Apex();
cc_f = Primitives.FindCircumcenter(torg, tdest, tapex, ref xi, ref eta);
cc_f = this.points[f.triangle.id];
if (f_next.triangle == Mesh.dummytri)
{
if (!f.triangle.infected)
{
// Add last circumcenter
P.Add(cc_f);
vpoints.Add(cc_f);
}
// Add midpoint of last triangles' edge (chances are it has already
// been added, so post process cell to remove duplicates???)
torg = f.Org();
tapex = f.Apex();
P.Add(new Point((torg.X + tapex.X) / 2, (torg.Y + tapex.Y) / 2));
p = new Point((torg.X + tapex.X) / 2, (torg.Y + tapex.Y) / 2);
p.id = n + segIndex;
points[n + segIndex] = p;
segIndex++;
vpoints.Add(p);
break;
}
torg = f_next.Org();
tdest = f_next.Dest();
tapex = f_next.Apex();
cc_f_next = Primitives.FindCircumcenter(torg, tdest, tapex, ref xi, ref eta);
cc_f_next = this.points[f_next.triangle.id];
// if f is tagged non-blind then
if (!f.triangle.infected)
{
// Insert the circumcenter of f into P
P.Add(cc_f);
vpoints.Add(cc_f);
if (f_next.triangle.infected)
{
// Call S_fnext the constrained edge blinding f_next
sfn.seg = subsegMap[f_next.triangle.hash];
p = new Point();
// Insert point Lf,f_next /\ Sf_next into P
if (SegmentsIntersect(ref sfn, cc_f, cc_f_next, ref p, true))
if (SegmentsIntersect(sfn.SegOrg(), sfn.SegDest(), cc_f, cc_f_next, out p, true))
{
P.Add(p);
p.id = n + segIndex;
points[n + segIndex] = p;
segIndex++;
vpoints.Add(p);
}
}
}
@@ -415,6 +474,9 @@ namespace TriangleNet.Tools
// Call Sf the constrained edge blinding f
sf.seg = subsegMap[f.triangle.hash];
sorg = sf.SegOrg();
sdest = sf.SegDest();
// if f_next is tagged non-blind then
if (!f_next.triangle.infected)
{
@@ -423,24 +485,28 @@ namespace TriangleNet.Tools
// Both circumcenters lie on the blinded side, but we
// have to add the intersection with the segment.
p = new Point();
// Center of f edge dest->apex
Point bisec = new Point((tdest.X + tapex.X) / 2, (tdest.Y + tapex.Y) / 2);
// Find intersection of seg with line through f's bisector and circumcenter
if (SegmentsIntersect(ref sf, bisec, cc_f, ref p, false))
if (SegmentsIntersect(sorg, sdest, bisec, cc_f, out p, false))
{
P.Add(p);
p.id = n + segIndex;
points[n + segIndex] = p;
segIndex++;
vpoints.Add(p);
}
p = new Point();
// Insert point Lf,f_next /\ Sf into P
if (SegmentsIntersect(ref sf, cc_f, cc_f_next, ref p, true))
if (SegmentsIntersect(sorg, sdest, cc_f, cc_f_next, out p, true))
{
P.Add(p);
p.id = n + segIndex;
points[n + segIndex] = p;
segIndex++;
vpoints.Add(p);
}
}
else
@@ -451,32 +517,41 @@ namespace TriangleNet.Tools
// if Sf != Sf_next then
if (!sf.Equal(sfn))
{
p = new Point();
// Insert Lf,fnext /\ Sf and Lf,fnext /\ Sfnext into P
if (SegmentsIntersect(ref sf, cc_f, cc_f_next, ref p, true))
if (SegmentsIntersect(sorg, sdest, cc_f, cc_f_next, out p, true))
{
P.Add(p);
p.id = n + segIndex;
points[n + segIndex] = p;
segIndex++;
vpoints.Add(p);
}
p = new Point();
if (SegmentsIntersect(ref sfn, cc_f, cc_f_next, ref p, true))
if (SegmentsIntersect(sfn.SegOrg(), sfn.SegDest(), cc_f, cc_f_next, out p, true))
{
P.Add(p);
p.id = n + segIndex;
points[n + segIndex] = p;
segIndex++;
vpoints.Add(p);
}
}
else
{
// Both circumcenters lie on the blinded side, but we
// have to add the intersection with the segment.
p = new Point();
// Center of f_next edge org->dest
Point bisec = new Point((torg.X + tdest.X) / 2, (torg.Y + tdest.Y) / 2);
// Find intersection of seg with line through f_next's bisector and circumcenter
if (SegmentsIntersect(ref sf, bisec, cc_f_next, ref p, false))
if (SegmentsIntersect(sorg, sdest, bisec, cc_f_next, out p, false))
{
P.Add(p);
p.id = n + segIndex;
points[n + segIndex] = p;
segIndex++;
vpoints.Add(p);
}
}
}
@@ -491,7 +566,7 @@ namespace TriangleNet.Tools
while (!f.Equal(f_init));
// Output: Bounded Voronoi cell of x in counterclockwise order.
return P.ToArray();
region.Add(vpoints);
}
/// <summary>
@@ -507,15 +582,14 @@ namespace TriangleNet.Tools
/// Returns false if there is no determinable intersection point, in which case X,Y will
/// be unmodified.
/// </returns>
private bool SegmentsIntersect(ref Osub seg, Point pc, Point pd, ref Point p, bool strictIntersect)
private bool SegmentsIntersect(Point p1, Point p2, Point p3, Point p4, out Point p, bool strictIntersect)
{
Vertex s1 = seg.Org();
Vertex s2 = seg.Dest();
p = null;
double Ax = s1.X, Ay = s1.Y;
double Bx = s2.X, By = s2.Y;
double Cx = pc.X, Cy = pc.Y;
double Dx = pd.X, Dy = pd.Y;
double Ax = p1.X, Ay = p1.Y;
double Bx = p2.X, By = p2.Y;
double Cx = p3.X, Cy = p3.Y;
double Dx = p4.X, Dy = p4.Y;
double distAB, theCos, theSin, newX, ABpos;
@@ -555,8 +629,7 @@ namespace TriangleNet.Tools
if (ABpos < 0 || ABpos > distAB && strictIntersect) return false;
// (4) Apply the discovered position to line A-B in the original coordinate system.
p.x = Ax + ABpos * theCos;
p.y = Ay + ABpos * theSin;
p = new Point(Ax + ABpos * theCos, Ay + ABpos * theSin);
// Success.
return true;
+12 -3
View File
@@ -6,15 +6,24 @@
namespace TriangleNet.Tools
{
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using TriangleNet.Geometry;
/// <summary>
/// TODO: Update summary.
/// </summary>
public interface IVoronoi
{
/// <summary>
/// Gets the list of Voronoi vertices.
/// </summary>
Point[] Points { get; }
/// <summary>
/// Gets the directions for infinite Voronoi edges.
/// </summary>
List<VoronoiRegion> Regions { get; }
void Generate();
}
}
+263 -83
View File
@@ -7,21 +7,25 @@
using TriangleNet.Data;
using TriangleNet.Geometry;
using System.Collections.Generic;
namespace TriangleNet.Tools
{
/// <summary>
/// The Voronoi Diagram is the dual of a pointset triangulation.
/// </summary>
public class Voronoi
public class Voronoi : IVoronoi
{
Mesh mesh;
Point[] points;
Point[] directions; // Stores the direction for infinite voronoi edges
Edge[] edges;
List<VoronoiRegion> regions;
Dictionary<int, Point> rayPoints;
int rayIndex;
/// <summary>
/// Initializes a new instance of the <see cref="BoundedVoronoi" /> class.
/// Initializes a new instance of the <see cref="Voronoi" /> class.
/// </summary>
/// <param name="mesh"></param>
/// <remarks>
@@ -41,19 +45,11 @@ namespace TriangleNet.Tools
}
/// <summary>
/// Gets the directions for infinite Voronoi edges.
/// Gets the list of Voronoi regions.
/// </summary>
public Point[] Directions
public List<VoronoiRegion> Regions
{
get { return directions; }
}
/// <summary>
/// Gets the list of Voronoi edges.
/// </summary>
public Edge[] Edges
{
get { return edges; }
get { return regions; }
}
/// <summary>
@@ -69,81 +65,265 @@ namespace TriangleNet.Tools
///</remarks>
public void Generate()
{
Otri tri = default(Otri), trisym = default(Otri);
Vertex torg, tdest, tapex;
Point circumcenter;
double xi = 0, eta = 0;
int p1, p2;
mesh.Renumber();
mesh.MakeVertexMap();
// Allocate memory for Voronoi vertices.
this.points = new Point[mesh.triangles.Count];
// Allocate space for voronoi diagram
this.points = new Point[mesh.triangles.Count + mesh.hullsize];
this.regions = new List<VoronoiRegion>(mesh.vertices.Count);
tri.orient = 0;
ComputeCircumCenters();
int i = 0;
foreach (var item in mesh.triangles.Values)
rayPoints = new Dictionary<int, Point>();
rayIndex = 0;
// Loop over the mesh vertices (Voronoi generators).
foreach (var item in mesh.vertices.Values)
{
tri.triangle = item;
torg = tri.Org();
tdest = tri.Dest();
tapex = tri.Apex();
circumcenter = Primitives.FindCircumcenter(torg, tdest, tapex, ref xi, ref eta);
// X and y coordinates.
this.points[i] = new Point(circumcenter.x, circumcenter.y);
// Update element id
tri.triangle.id = i++;
}
// Allocate memory for output Voronoi edges.
this.edges = new Edge[mesh.edges];
// Allocate memory for output Voronoi norms.
this.directions = new Point[mesh.edges];
i = 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 pointer than its neighbor. This way, each edge is
// considered only once.
foreach (var item in mesh.triangles.Values)
{
tri.triangle = item;
for (tri.orient = 0; tri.orient < 3; tri.orient++)
//if (item.Boundary == 0)
{
tri.Sym(ref trisym);
if ((tri.triangle.id < trisym.triangle.id) || (trisym.triangle == Mesh.dummytri))
{
// Find the number of this triangle (and Voronoi vertex).
p1 = tri.triangle.id;
if (trisym.triangle == Mesh.dummytri)
{
torg = tri.Org();
tdest = tri.Dest();
// Copy an infinite ray. Index of one endpoint, and -1.
this.edges[i] = new Edge(p1, -1);
this.directions[i] = new Point(tdest.y - torg.y, torg.x - tdest.x);
}
else
{
// Find the number of the adjacent triangle (and Voronoi vertex).
p2 = trisym.triangle.id;
// Finite edge. Write indices of two endpoints.
this.edges[i] = new Edge(p1, p2);
this.directions[i] = new Point(0, 0);
}
i++;
}
ConstructVoronoiRegion(item);
}
}
}
private void ComputeCircumCenters()
{
Otri tri = default(Otri);
double xi = 0, eta = 0;
Point pt;
// Compue triangle circumcenters
foreach (var item in mesh.triangles.Values)
{
tri.triangle = item;
pt = Primitives.FindCircumcenter(tri.Org(), tri.Dest(), tri.Apex(), ref xi, ref eta);
pt.id = item.id;
points[item.id] = pt;
}
}
/// <summary>
///
/// </summary>
/// <param name="vertex"></param>
/// <returns>The circumcenter indices which make up the cell.</returns>
private void ConstructVoronoiRegion(Vertex vertex)
{
VoronoiRegion region = new VoronoiRegion(vertex);
regions.Add(region);
List<Point> vpoints = new List<Point>();
Otri f = default(Otri);
Otri f_init = default(Otri);
Otri f_next = default(Otri);
Otri f_prev = default(Otri);
Osub sub = default(Osub);
// Call f_init a triangle incident to x
vertex.tri.Copy(ref f_init);
f_init.Copy(ref f);
f_init.Onext(ref f_next);
if (f_next.triangle == Mesh.dummytri)
{
f_init.Oprev(ref f_prev);
if (f_prev.triangle != Mesh.dummytri)
{
f_init.Copy(ref f_next);
// Move one triangle clockwise
f_init.OprevSelf();
f_init.Copy(ref f);
}
}
// Go counterclockwise until we reach the border or the initial triangle.
while (f_next.triangle != Mesh.dummytri)
{
// Add circumcenter of current triangle
vpoints.Add(points[f.triangle.id]);
if (f_next.Equal(f_init))
{
// Voronoi cell is complete (bounded case).
region.Add(vpoints);
return;
}
f_next.Copy(ref f);
f_next.OnextSelf();
}
// Voronoi cell is unbounded
region.Bounded = false;
Vertex torg, tdest, tapex, intersection;
int sid, n = mesh.triangles.Count;
// Find the boundary segment id.
f.Lprev(ref f_next);
f_next.SegPivot(ref sub);
sid = sub.seg.hash;
// Last valid f lies at the boundary. Add the circumcenter.
vpoints.Add(points[f.triangle.id]);
// Check if the intersection with the bounding box has already been computed.
if (rayPoints.ContainsKey(sid))
{
vpoints.Add(rayPoints[sid]);
}
else
{
torg = f.Org();
tapex = f.Apex();
BoxRayIntersection(points[f.triangle.id], torg.y - tapex.y, tapex.x - torg.x, out intersection);
// Set the correct id for the vertex
intersection.id = n + rayIndex;
points[n + rayIndex] = intersection;
rayIndex++;
vpoints.Add(intersection);
rayPoints.Add(sid, intersection);
}
// Now walk from f_init clockwise till we reach the boundary.
vpoints.Reverse();
f_init.Copy(ref f);
f.Oprev(ref f_prev);
while (f_prev.triangle != Mesh.dummytri)
{
vpoints.Add(points[f_prev.triangle.id]);
f_prev.Copy(ref f);
f_prev.OprevSelf();
}
// Find the boundary segment id.
f.SegPivot(ref sub);
sid = sub.seg.hash;
// Last valid f lies at the boundary. Add the circumcenter.
vpoints.Add(points[f.triangle.id]);
if (rayPoints.ContainsKey(sid))
{
vpoints.Add(rayPoints[sid]);
}
else
{
// Intersection has not been computed yet.
torg = f.Org();
tdest = f.Dest();
BoxRayIntersection(points[f.triangle.id], tdest.y - torg.y, torg.x - tdest.x, out intersection);
// Set the correct id for the vertex
intersection.id = n + rayIndex;
points[n + rayIndex] = intersection;
rayIndex++;
vpoints.Add(intersection);
rayPoints.Add(sid, intersection);
}
// Add the new points to the region (in counter-clockwise order)
vpoints.Reverse();
region.Add(vpoints);
}
private bool BoxRayIntersection(Point pt, double dx, double dy, out Vertex intersect)
{
double x = pt.X;
double y = pt.Y;
double t1, x1, y1, t2, x2, y2;
// Bounding box (50% enlarged)
var box = mesh.Bounds;
double dw = box.Width * 0.5f;
double dh = box.Height * 0.5f;
double minX = box.Xmin - dw;
double maxX = box.Xmax + dw;
double minY = box.Ymin - dh;
double maxY = box.Ymax + dh;
// Check if point is inside the bounds
if (x < minX || x > maxX || y < minY || y > maxY)
{
intersect = null;
return false;
}
// Calculate the cut through the vertical boundaries
if (dx < 0)
{
// Line going to the left: intersect with x = minX
t1 = (minX - x) / dx;
x1 = minX;
y1 = y + t1 * dy;
}
else if (dx > 0)
{
// Line going to the right: intersect with x = maxX
t1 = (maxX - x) / dx;
x1 = maxX;
y1 = y + t1 * dy;
}
else
{
// Line going straight up or down: no intersection possible
t1 = double.MaxValue;
x1 = y1 = 0;
}
// Calculate the cut through upper and lower boundaries
if (dy < 0)
{
// Line going downwards: intersect with y = minY
t2 = (minY - y) / dy;
x2 = x + t2 * dx;
y2 = minY;
}
else if (dx > 0)
{
// Line going upwards: intersect with y = maxY
t2 = (maxY - y) / dy;
x2 = x + t2 * dx;
y2 = maxY;
}
else
{
// Horizontal line: no intersection possible
t2 = double.MaxValue;
x2 = y2 = 0;
}
if (t1 < t2)
{
intersect = new Vertex(x1, y1, -1);
}
else
{
intersect = new Vertex(x2, y2, -1);
}
return true;
}
}
}
@@ -0,0 +1,82 @@
// -----------------------------------------------------------------------
// <copyright file="VoronoiRegion.cs" company="">
// TODO: Update copyright text.
// </copyright>
// -----------------------------------------------------------------------
namespace TriangleNet.Tools
{
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using TriangleNet.Geometry;
using TriangleNet.Data;
/// <summary>
/// Represents a region in the Voronoi diagram.
/// </summary>
public class VoronoiRegion
{
int id;
Point generator;
List<Point> vertices;
bool bounded;
/// <summary>
/// Gets the Voronoi region id (which is the same as the generators vertex id).
/// </summary>
public int ID
{
get { return id; }
}
/// <summary>
/// Gets the Voronoi regions generator.
/// </summary>
public Point Generator
{
get { return generator; }
}
/// <summary>
/// Gets the Voronoi vertices on the regions boundary.
/// </summary>
public ICollection<Point> Vertices
{
get { return vertices; }
}
/// <summary>
/// Gets or sets whether the Voronoi region is bounded.
/// </summary>
public bool Bounded
{
get { return bounded; }
set { bounded = value; }
}
public VoronoiRegion(Vertex generator)
{
this.id = generator.id;
this.generator = generator;
this.vertices = new List<Point>();
this.bounded = true;
}
public void Add(Point point)
{
this.vertices.Add(point);
}
public void Add(List<Point> points)
{
this.vertices.AddRange(points);
}
public override string ToString()
{
return String.Format("R-ID {0}", id);
}
}
}
+1
View File
@@ -88,6 +88,7 @@
<Compile Include="Tools\Statistic.cs" />
<Compile Include="Algorithm\SweepLine.cs" />
<Compile Include="Tools\Voronoi.cs" />
<Compile Include="Tools\VoronoiRegion.cs" />
</ItemGroup>
<ItemGroup />
<Import Project="$(MSBuildToolsPath)\Microsoft.CSharp.targets" />