From 04c45ed41ad33402f8ff95a2c9acc8972221433c Mon Sep 17 00:00:00 2001 From: "SND\\wo80_cp" Date: Sun, 14 Oct 2012 22:08:33 +0000 Subject: [PATCH] Improved uniform data structure for Voronoi classes. git-svn-id: https://triangle.svn.codeplex.com/svn@69976 0e2699bc-83d4-4a8f-98e7-55e24ab8c7a5 --- Triangle.NET/MeshRenderer.Core/RenderData.cs | 51 ++- Triangle.NET/TestApp/FormMain.cs | 26 +- Triangle.NET/Triangle/Data/Vertex.cs | 11 - Triangle.NET/Triangle/Geometry/Point.cs | 9 + Triangle.NET/Triangle/Tools/BoundedVoronoi.cs | 267 +++++++++----- Triangle.NET/Triangle/Tools/IVoronoi.cs | 15 +- Triangle.NET/Triangle/Tools/Voronoi.cs | 346 +++++++++++++----- Triangle.NET/Triangle/Tools/VoronoiRegion.cs | 82 +++++ Triangle.NET/Triangle/Triangle.csproj | 1 + 9 files changed, 610 insertions(+), 198 deletions(-) create mode 100644 Triangle.NET/Triangle/Tools/VoronoiRegion.cs diff --git a/Triangle.NET/MeshRenderer.Core/RenderData.cs b/Triangle.NET/MeshRenderer.Core/RenderData.cs index be99868..c9a7084 100644 --- a/Triangle.NET/MeshRenderer.Core/RenderData.cs +++ b/Triangle.NET/MeshRenderer.Core/RenderData.cs @@ -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 /// 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(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; diff --git a/Triangle.NET/TestApp/FormMain.cs b/Triangle.NET/TestApp/FormMain.cs index 4fdc11b..00c6f47 100644 --- a/Triangle.NET/TestApp/FormMain.cs +++ b/Triangle.NET/TestApp/FormMain.cs @@ -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) diff --git a/Triangle.NET/Triangle/Data/Vertex.cs b/Triangle.NET/Triangle/Data/Vertex.cs index 4b25ffd..b09193c 100644 --- a/Triangle.NET/Triangle/Data/Vertex.cs +++ b/Triangle.NET/Triangle/Data/Vertex.cs @@ -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 - /// - /// Gets the vertex id. - /// - public int ID - { - get { return this.id; } - } - /// /// Gets the vertex type. /// diff --git a/Triangle.NET/Triangle/Geometry/Point.cs b/Triangle.NET/Triangle/Geometry/Point.cs index 8dabcc2..21a2332 100644 --- a/Triangle.NET/Triangle/Geometry/Point.cs +++ b/Triangle.NET/Triangle/Geometry/Point.cs @@ -16,6 +16,7 @@ namespace TriangleNet.Geometry /// public class Point : IComparable, IEquatable { + internal int id; internal double x; internal double y; internal int mark; @@ -40,6 +41,14 @@ namespace TriangleNet.Geometry #region Public properties + /// + /// Gets the vertex id. + /// + public int ID + { + get { return this.id; } + } + /// /// Gets the vertex x coordinate. /// diff --git a/Triangle.NET/Triangle/Tools/BoundedVoronoi.cs b/Triangle.NET/Triangle/Tools/BoundedVoronoi.cs index ad31964..3bce64f 100644 --- a/Triangle.NET/Triangle/Tools/BoundedVoronoi.cs +++ b/Triangle.NET/Triangle/Tools/BoundedVoronoi.cs @@ -20,11 +20,16 @@ namespace TriangleNet.Tools /// 2D Centroidal Voronoi Tessellations with Constraints, 2010, /// Jane Tournois, Pierre Alliez and Olivier Devillers /// - public class BoundedVoronoi + public class BoundedVoronoi : IVoronoi { Mesh mesh; + + Point[] points; + List regions; + + int segIndex; + Dictionary subsegMap; - List voronoi; /// /// Initializes a new instance of the class. @@ -33,19 +38,22 @@ namespace TriangleNet.Tools public BoundedVoronoi(Mesh mesh) { this.mesh = mesh; - this.voronoi = new List(); - - // TODO: introduce isVertexMapValid in mesh and don't call - // MakeVertexMap() if not necessary. - this.mesh.MakeVertexMap(); } /// - /// Gets the list of Voronoi cells. + /// Gets the list of Voronoi vertices. /// - public List Cells + public Point[] Points { - get { return voronoi; } + get { return points; } + } + + /// + /// Gets the list of Voronoi regions. + /// + public List Regions + { + get { return regions; } } /// @@ -53,7 +61,14 @@ namespace TriangleNet.Tools /// 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(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; + } + } + /// /// Tag all blind triangles. /// @@ -169,26 +202,28 @@ namespace TriangleNet.Tools /// Returns true, if the triangle is blinded. 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 P = new List(); + List vpoints = new List(); // 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 P = new List(); + List vpoints = new List(); // 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); } /// @@ -507,15 +582,14 @@ namespace TriangleNet.Tools /// Returns false if there is no determinable intersection point, in which case X,Y will /// be unmodified. /// - 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; diff --git a/Triangle.NET/Triangle/Tools/IVoronoi.cs b/Triangle.NET/Triangle/Tools/IVoronoi.cs index f534589..ca63f50 100644 --- a/Triangle.NET/Triangle/Tools/IVoronoi.cs +++ b/Triangle.NET/Triangle/Tools/IVoronoi.cs @@ -6,15 +6,24 @@ namespace TriangleNet.Tools { - using System; using System.Collections.Generic; - using System.Linq; - using System.Text; + using TriangleNet.Geometry; /// /// TODO: Update summary. /// public interface IVoronoi { + /// + /// Gets the list of Voronoi vertices. + /// + Point[] Points { get; } + + /// + /// Gets the directions for infinite Voronoi edges. + /// + List Regions { get; } + + void Generate(); } } diff --git a/Triangle.NET/Triangle/Tools/Voronoi.cs b/Triangle.NET/Triangle/Tools/Voronoi.cs index cfc3629..370cd0f 100644 --- a/Triangle.NET/Triangle/Tools/Voronoi.cs +++ b/Triangle.NET/Triangle/Tools/Voronoi.cs @@ -7,21 +7,25 @@ using TriangleNet.Data; using TriangleNet.Geometry; +using System.Collections.Generic; namespace TriangleNet.Tools { /// /// The Voronoi Diagram is the dual of a pointset triangulation. /// - public class Voronoi + public class Voronoi : IVoronoi { Mesh mesh; + Point[] points; - Point[] directions; // Stores the direction for infinite voronoi edges - Edge[] edges; + List regions; + + Dictionary rayPoints; + int rayIndex; /// - /// Initializes a new instance of the class. + /// Initializes a new instance of the class. /// /// /// @@ -41,19 +45,11 @@ namespace TriangleNet.Tools } /// - /// Gets the directions for infinite Voronoi edges. + /// Gets the list of Voronoi regions. /// - public Point[] Directions + public List Regions { - get { return directions; } - } - - /// - /// Gets the list of Voronoi edges. - /// - public Edge[] Edges - { - get { return edges; } + get { return regions; } } /// @@ -69,81 +65,265 @@ namespace TriangleNet.Tools /// 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(mesh.vertices.Count); - tri.orient = 0; + ComputeCircumCenters(); - int i = 0; - foreach (var item in mesh.triangles.Values) + rayPoints = new Dictionary(); + 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; + } + } + + /// + /// + /// + /// + /// The circumcenter indices which make up the cell. + private void ConstructVoronoiRegion(Vertex vertex) + { + VoronoiRegion region = new VoronoiRegion(vertex); + regions.Add(region); + + List vpoints = new List(); + + 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; + } } } diff --git a/Triangle.NET/Triangle/Tools/VoronoiRegion.cs b/Triangle.NET/Triangle/Tools/VoronoiRegion.cs new file mode 100644 index 0000000..0d00471 --- /dev/null +++ b/Triangle.NET/Triangle/Tools/VoronoiRegion.cs @@ -0,0 +1,82 @@ +// ----------------------------------------------------------------------- +// +// TODO: Update copyright text. +// +// ----------------------------------------------------------------------- + +namespace TriangleNet.Tools +{ + using System; + using System.Collections.Generic; + using System.Linq; + using System.Text; + using TriangleNet.Geometry; + using TriangleNet.Data; + + /// + /// Represents a region in the Voronoi diagram. + /// + public class VoronoiRegion + { + int id; + Point generator; + List vertices; + bool bounded; + + /// + /// Gets the Voronoi region id (which is the same as the generators vertex id). + /// + public int ID + { + get { return id; } + } + + /// + /// Gets the Voronoi regions generator. + /// + public Point Generator + { + get { return generator; } + } + + /// + /// Gets the Voronoi vertices on the regions boundary. + /// + public ICollection Vertices + { + get { return vertices; } + } + + /// + /// Gets or sets whether the Voronoi region is bounded. + /// + public bool Bounded + { + get { return bounded; } + set { bounded = value; } + } + + public VoronoiRegion(Vertex generator) + { + this.id = generator.id; + this.generator = generator; + this.vertices = new List(); + this.bounded = true; + } + + public void Add(Point point) + { + this.vertices.Add(point); + } + + public void Add(List points) + { + this.vertices.AddRange(points); + } + + public override string ToString() + { + return String.Format("R-ID {0}", id); + } + } +} diff --git a/Triangle.NET/Triangle/Triangle.csproj b/Triangle.NET/Triangle/Triangle.csproj index c4a7680..7b9aad5 100644 --- a/Triangle.NET/Triangle/Triangle.csproj +++ b/Triangle.NET/Triangle/Triangle.csproj @@ -88,6 +88,7 @@ +