From c0d89aebebe3127a72bfe063e72ccc581880dcdb Mon Sep 17 00:00:00 2001 From: "SND\\wo80_cp" Date: Fri, 8 Jun 2012 20:21:45 +0000 Subject: [PATCH] Added (bounded) Voronoi diagram class git-svn-id: https://triangle.svn.codeplex.com/svn@67930 0e2699bc-83d4-4a8f-98e7-55e24ab8c7a5 --- .../TestApp/Controls/DarkToolStripRenderer.cs | 19 + Triangle.NET/TestApp/Controls/MeshRenderer.cs | 25 + Triangle.NET/TestApp/Examples.cs | 1 + Triangle.NET/TestApp/{ => IO}/ImageWriter.cs | 30 +- Triangle.NET/TestApp/Mesh Explorer.csproj | 3 +- .../TestApp/Properties/AssemblyInfo.cs | 4 +- .../TestApp/Rendering/VoronoiRenderer.cs | 273 +++++++++ Triangle.NET/TestApp/Rendering/Zoom.cs | 2 +- Triangle.NET/Triangle/IO/VoronoiData.cs | 24 - Triangle.NET/Triangle/Mesh.cs | 7 +- Triangle.NET/Triangle/Quality.cs | 20 +- Triangle.NET/Triangle/Tools/BoundedVoronoi.cs | 551 ++++++++++++++++++ Triangle.NET/Triangle/Tools/QualityMeasure.cs | 541 +++++++++++++++++ .../{IO/DataWriter.cs => Tools/Voronoi.cs} | 112 ++-- Triangle.NET/Triangle/Triangle.csproj | 5 +- 15 files changed, 1511 insertions(+), 106 deletions(-) rename Triangle.NET/TestApp/{ => IO}/ImageWriter.cs (93%) create mode 100644 Triangle.NET/TestApp/Rendering/VoronoiRenderer.cs delete mode 100644 Triangle.NET/Triangle/IO/VoronoiData.cs create mode 100644 Triangle.NET/Triangle/Tools/BoundedVoronoi.cs create mode 100644 Triangle.NET/Triangle/Tools/QualityMeasure.cs rename Triangle.NET/Triangle/{IO/DataWriter.cs => Tools/Voronoi.cs} (62%) diff --git a/Triangle.NET/TestApp/Controls/DarkToolStripRenderer.cs b/Triangle.NET/TestApp/Controls/DarkToolStripRenderer.cs index d55cae6..2d09d35 100644 --- a/Triangle.NET/TestApp/Controls/DarkToolStripRenderer.cs +++ b/Triangle.NET/TestApp/Controls/DarkToolStripRenderer.cs @@ -12,6 +12,7 @@ namespace MeshExplorer.Controls using System.Text; using System.Windows.Forms; using System.Drawing; + using System.Drawing.Drawing2D; /// /// TODO: Update summary. @@ -25,6 +26,24 @@ namespace MeshExplorer.Controls base.OnRenderItemText(e); } + protected override void OnRenderItemCheck(ToolStripItemImageRenderEventArgs e) + { + int h = e.Item.Height; + int size = h / 2 - 2; // Box size + Rectangle rect = new Rectangle(10, (h - size) / 2, size - 2, size); + + var mode = e.Graphics.SmoothingMode; + + using (Pen pen = new Pen(Color.White, 1.6f)) + { + e.Graphics.SmoothingMode = SmoothingMode.AntiAlias; + e.Graphics.DrawLine(pen, rect.Left, rect.Bottom - size / 2, rect.Left + size / 2.5f, rect.Bottom - 2); + e.Graphics.DrawLine(pen, rect.Left + size / 2.6f, rect.Bottom - 2, rect.Right, rect.Top); + } + + e.Graphics.SmoothingMode = mode; + } + protected override void OnRenderArrow(ToolStripArrowRenderEventArgs e) { if (!e.Item.Selected && !e.Item.Pressed) diff --git a/Triangle.NET/TestApp/Controls/MeshRenderer.cs b/Triangle.NET/TestApp/Controls/MeshRenderer.cs index 98b72b8..f5a7d72 100644 --- a/Triangle.NET/TestApp/Controls/MeshRenderer.cs +++ b/Triangle.NET/TestApp/Controls/MeshRenderer.cs @@ -32,6 +32,8 @@ namespace MeshExplorer.Controls Zoom zoom; RenderData data; bool initialized = false; + VoronoiRenderer voronoi; + bool showVoronoi = false; string coordinate = String.Empty; @@ -39,6 +41,21 @@ namespace MeshExplorer.Controls public long RenderTime { get; private set; } public RenderData Data { get { return data; } } + public bool ShowVoronoi + { + get { return showVoronoi; } + set + { + showVoronoi = value; + + if (voronoi != null && showVoronoi) + { + voronoi.Update(); + } + + this.Render(); + } + } public MeshRenderer() { @@ -83,6 +100,9 @@ namespace MeshExplorer.Controls public void SetData(Mesh mesh) { + voronoi = new VoronoiRenderer(mesh); + voronoi.Update(); + data.SetData(mesh); initialized = true; @@ -259,6 +279,11 @@ namespace MeshExplorer.Controls this.RenderTriangles(g); } + if (voronoi != null && this.showVoronoi) + { + voronoi.Render(g, zoom); + } + if (data.Segments != null) { this.RenderSegments(g); diff --git a/Triangle.NET/TestApp/Examples.cs b/Triangle.NET/TestApp/Examples.cs index 9ce16a4..d2de5f5 100644 --- a/Triangle.NET/TestApp/Examples.cs +++ b/Triangle.NET/TestApp/Examples.cs @@ -13,6 +13,7 @@ namespace MeshExplorer using TriangleNet; using TriangleNet.IO; using TriangleNet.Geometry; + using MeshExplorer.IO; /// /// Code of the online examples. diff --git a/Triangle.NET/TestApp/ImageWriter.cs b/Triangle.NET/TestApp/IO/ImageWriter.cs similarity index 93% rename from Triangle.NET/TestApp/ImageWriter.cs rename to Triangle.NET/TestApp/IO/ImageWriter.cs index b3037cb..8101c45 100644 --- a/Triangle.NET/TestApp/ImageWriter.cs +++ b/Triangle.NET/TestApp/IO/ImageWriter.cs @@ -4,7 +4,7 @@ // // ----------------------------------------------------------------------- -namespace MeshExplorer +namespace MeshExplorer.IO { using System; using System.Drawing; @@ -14,6 +14,7 @@ namespace MeshExplorer using TriangleNet; using TriangleNet.Data; using TriangleNet.IO; + using TriangleNet.Tools; /// /// Writes an image of the mesh to disk. @@ -187,7 +188,7 @@ namespace MeshExplorer filename = String.Format("mesh-{0}.png", DateTime.Now.ToString("yyyy-M-d-hh-mm-ss")); } - VoronoiData data = DataWriter.WriteVoronoi(mesh); + Voronoi data = null; // DataWriter.WriteVoronoi(mesh); int n = data.Points.Length; @@ -332,7 +333,7 @@ namespace MeshExplorer /// /// Draw mesh to the graphics object. /// - private static void DrawVoronoi(Graphics g, Mesh mesh, VoronoiData voronoi, float scale) + private static void DrawVoronoi(Graphics g, Mesh mesh, Voronoi voronoi, float scale) { g.SmoothingMode = SmoothingMode.AntiAlias; // Colors @@ -399,16 +400,16 @@ namespace MeshExplorer } // Draw input points - n = voronoi.InputPoints.Length; + //n = voronoi.InputPoints.Length; - for (int i = 0; i < n; i++) - { - x = (float)voronoi.InputPoints[i].X; - y = (float)voronoi.InputPoints[i].Y; + //for (int i = 0; i < n; i++) + //{ + // x = (float)voronoi.InputPoints[i].X; + // y = (float)voronoi.InputPoints[i].Y; - g.FillEllipse(spBrush, x - radius, y - radius, 2 * radius, 2 * radius); - //g.DrawEllipse(spBrush, x - radius, y - radius, 2 * radius, 2 * radius); - } + // g.FillEllipse(spBrush, x - radius, y - radius, 2 * radius, 2 * radius); + // //g.DrawEllipse(spBrush, x - radius, y - radius, 2 * radius, 2 * radius); + //} bgBrush.Dispose(); ptBrush.Dispose(); @@ -418,12 +419,13 @@ namespace MeshExplorer trBrush.Dispose(); } - private static PointF VoronoiBoxIntersection(BBox bounds, TriangleNet.Geometry.Point pt, double[] direction) + private static PointF VoronoiBoxIntersection(BBox bounds, TriangleNet.Geometry.Point pt, + TriangleNet.Geometry.Point direction) { double x = pt.X; double y = pt.Y; - double dx = direction[0]; - double dy = direction[1]; + double dx = direction.X; + double dy = direction.Y; double t1, x1, y1, t2, x2, y2; diff --git a/Triangle.NET/TestApp/Mesh Explorer.csproj b/Triangle.NET/TestApp/Mesh Explorer.csproj index 49a561a..f5d38be 100644 --- a/Triangle.NET/TestApp/Mesh Explorer.csproj +++ b/Triangle.NET/TestApp/Mesh Explorer.csproj @@ -94,7 +94,7 @@ FormQuality.cs - + @@ -105,6 +105,7 @@ + diff --git a/Triangle.NET/TestApp/Properties/AssemblyInfo.cs b/Triangle.NET/TestApp/Properties/AssemblyInfo.cs index 8b41fd1..d240e61 100644 --- a/Triangle.NET/TestApp/Properties/AssemblyInfo.cs +++ b/Triangle.NET/TestApp/Properties/AssemblyInfo.cs @@ -5,11 +5,11 @@ using System.Runtime.InteropServices; // General Information about an assembly is controlled through the following // set of attributes. Change these attribute values to modify the information // associated with an assembly. -[assembly: AssemblyTitle("Example1")] +[assembly: AssemblyTitle("Mesh Explorer")] [assembly: AssemblyDescription("")] [assembly: AssemblyConfiguration("")] [assembly: AssemblyCompany("")] -[assembly: AssemblyProduct("Example1")] +[assembly: AssemblyProduct("Mesh Explorer")] [assembly: AssemblyCopyright("Copyright © 2012")] [assembly: AssemblyTrademark("")] [assembly: AssemblyCulture("")] diff --git a/Triangle.NET/TestApp/Rendering/VoronoiRenderer.cs b/Triangle.NET/TestApp/Rendering/VoronoiRenderer.cs new file mode 100644 index 0000000..fad33b4 --- /dev/null +++ b/Triangle.NET/TestApp/Rendering/VoronoiRenderer.cs @@ -0,0 +1,273 @@ +// ----------------------------------------------------------------------- +// +// Christian Woltering, http://home.edo.tu-dortmund.de/~woltering/triangle/ +// +// ----------------------------------------------------------------------- + +namespace MeshExplorer.Rendering +{ + using System; + using System.Collections.Generic; + using System.Drawing; + using System.Linq; + using System.Text; + using TriangleNet; + using TriangleNet.Tools; + + /// + /// Renders a (bounded) Voronoi diagram. + /// + public class VoronoiRenderer + { + Mesh mesh; + Voronoi simpleVoro; + BoundedVoronoi boundedVoro; + + Pen lines = new Pen(Color.FromArgb(40, 50, 60)); + + public VoronoiRenderer(Mesh mesh) + { + this.mesh = mesh; + + //if (mesh.NumberOfSegments > 0) + if (mesh.IsPolygon) + { + boundedVoro = new BoundedVoronoi(mesh); + } + else + { + simpleVoro = new Voronoi(mesh); + } + } + + public void Update() + { + if (simpleVoro != null) + { + simpleVoro.Generate(); + } + + if (boundedVoro != null) + { + boundedVoro.Generate(); + } + } + + public void Reset() + { + simpleVoro = null; + boundedVoro = null; + } + + internal void Render(Graphics g, Zoom zoom) + { + if (simpleVoro != null) + { + RenderSimple(g, zoom); + return; + } + + if (boundedVoro != null) + { + RenderBounded(g, zoom); + return; + } + } + + private void RenderSimple(Graphics g, Zoom zoom) + { + PointF p0, p1; + + BBox bounds = new BBox(mesh.Bounds); + + TriangleNet.Geometry.Point[] points = simpleVoro.Points; + + // Enlarge 50% + bounds.Extend(0.5f); + + // Draw edges + int n = simpleVoro.Edges == null ? 0 : simpleVoro.Edges.Length; + + for (int i = 0; i < n; i++) + { + var seg = simpleVoro.Edges[i]; + + if (seg.P1 == -1) + { + // Infinite voronoi edge + p0 = new PointF((float)points[seg.P0].X, (float)points[seg.P0].Y); + + if (!BoxRayIntersection(bounds, points[seg.P0], simpleVoro.Directions[i], out p1)) + { + continue; + } + } + else + { + p0 = new PointF((float)points[seg.P0].X, (float)points[seg.P0].Y); + p1 = new PointF((float)points[seg.P1].X, (float)points[seg.P1].Y); + } + + // Draw line + RenderEdge(g, zoom, p0, p1); + } + + // Shrink 50% + bounds.Extend(-0.5f); + + // Scale the points radius to 2 pixel. + //float radius = 1.5f / scale, x, y; + + // Draw points + //n = voro.Points.Length; + + //for (int i = 0; i < n; i++) + //{ + // x = (float)voro.Points[i].X; + // y = (float)voro.Points[i].Y; + + // g.FillEllipse(Brushes.Blue, x - radius, y - radius, 2 * radius, 2 * radius); + //} + } + + private void RenderBounded(Graphics g, Zoom zoom) + { + PointF p0, p1; + int n; + + foreach (var cell in boundedVoro.Cells) + { + n = cell.Length; + + for (int i = 1; i < n; i++) + { + p0 = new PointF((float)cell[i - 1].X, (float)cell[i - 1].Y); + p1 = new PointF((float)cell[i].X, (float)cell[i].Y); + + RenderEdge(g, zoom, p0, p1); + } + } + } + + private void RenderEdge(Graphics g, Zoom zoom, PointF p0, PointF p1) + { + if (zoom.ViewportContains(p0) || + zoom.ViewportContains(p1)) + { + p0 = zoom.WorldToScreen(p0); + p1 = zoom.WorldToScreen(p1); + + g.DrawLine(lines, p0, p1); + } + } + + private bool BoxRayIntersection(BBox bounds, TriangleNet.Geometry.Point pt, + TriangleNet.Geometry.Point direction, out PointF intersect) + { + double x = pt.X; + double y = pt.Y; + double dx = direction.X; + double dy = direction.Y; + + double t1, x1, y1, t2, x2, y2; + + intersect = new PointF(); + + // Check if point is inside the bounds + if (x < bounds.MinX || x > bounds.MaxX || y < bounds.MinY || y > bounds.MaxY) + { + return false; + } + + // Calculate the cut through the vertical boundaries + if (dx < 0) + { + // Line going to the left: intersect with x = bounds.MinX + t1 = (bounds.MinX - x) / dx; + x1 = bounds.MinX; + y1 = y + t1 * dy; + } + else if (dx > 0) + { + // Line going to the right: intersect with x = bounds.MaxX + t1 = (bounds.MaxX - x) / dx; + x1 = bounds.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 = bounds.MinY + t2 = (bounds.MinY - y) / dy; + x2 = x + t2 * dx; + y2 = bounds.MinY; + } + else if (dx > 0) + { + // Line going upwards: intersect with y = bounds.MaxY + t2 = (bounds.MaxY - y) / dy; + x2 = x + t2 * dx; + y2 = bounds.MaxY; + } + else + { + // Horizontal line: no intersection possible + t2 = double.MaxValue; + x2 = y2 = 0; + } + + if (t1 < t2) + { + intersect.X = (float)x1; + intersect.Y = (float)y1; + } + else + { + intersect.X = (float)x2; + intersect.Y = (float)y2; + } + + return true; + } + + /// + /// Bounding box. + /// + struct BBox + { + public float MinX; + public float MaxX; + public float MinY; + public float MaxY; + + public float Width { get { return MaxX - MinX; } } + public float Height { get { return MaxY - MinY; } } + + public BBox(TriangleNet.Geometry.BoundingBox box) + { + MinX = (float)box.Xmin; + MaxX = (float)box.Xmax; + MinY = (float)box.Ymin; + MaxY = (float)box.Ymax; + } + + public void Extend(float amount) + { + float dx = amount * this.Width; + float dy = amount * this.Height; + + MinX -= dx; + MaxX += dx; + MinY -= dy; + MaxY += dy; + } + } + } +} diff --git a/Triangle.NET/TestApp/Rendering/Zoom.cs b/Triangle.NET/TestApp/Rendering/Zoom.cs index 5512a5e..7333a9e 100644 --- a/Triangle.NET/TestApp/Rendering/Zoom.cs +++ b/Triangle.NET/TestApp/Rendering/Zoom.cs @@ -15,7 +15,7 @@ namespace MeshExplorer.Rendering /// /// Manages the current world to screen transformation /// - class Zoom + internal class Zoom { // The complete mesh Rectangle Screen { get; set; } diff --git a/Triangle.NET/Triangle/IO/VoronoiData.cs b/Triangle.NET/Triangle/IO/VoronoiData.cs deleted file mode 100644 index 41dc354..0000000 --- a/Triangle.NET/Triangle/IO/VoronoiData.cs +++ /dev/null @@ -1,24 +0,0 @@ -// ----------------------------------------------------------------------- -// -// Original Triangle code by Jonathan Richard Shewchuk, http://www.cs.cmu.edu/~quake/triangle.html -// Triangle.NET code by Christian Woltering, http://home.edo.tu-dortmund.de/~woltering/index.html -// -// ----------------------------------------------------------------------- - -namespace TriangleNet.IO -{ - using TriangleNet.Geometry; - - /// - /// Stores the voronoi data (output only). - /// - public class VoronoiData - { - public Point[] InputPoints; - public Point[] Points; - public Edge[] Edges; - - // Stores the direction for infinite voronoi edges - public double[][] Directions; - } -} \ No newline at end of file diff --git a/Triangle.NET/Triangle/Mesh.cs b/Triangle.NET/Triangle/Mesh.cs index 74563fd..25ff5b6 100644 --- a/Triangle.NET/Triangle/Mesh.cs +++ b/Triangle.NET/Triangle/Mesh.cs @@ -152,6 +152,11 @@ namespace TriangleNet /// public int NumberOfEdges { get { return this.edges; } } + /// + /// Indicates whether the input is a PSLG or a point set. + /// + public bool IsPolygon { get { return this.insegments > 0; } } + #endregion /// @@ -2040,7 +2045,7 @@ namespace TriangleNet /// triangles, but in the end every vertex will point to some triangle /// that contains it. /// - private void MakeVertexMap() + internal void MakeVertexMap() { Otri tri = default(Otri); Vertex triorg; diff --git a/Triangle.NET/Triangle/Quality.cs b/Triangle.NET/Triangle/Quality.cs index f1db6a6..ed1c273 100644 --- a/Triangle.NET/Triangle/Quality.cs +++ b/Triangle.NET/Triangle/Quality.cs @@ -142,7 +142,7 @@ namespace TriangleNet /// public bool CheckDelaunay() { - Otri triangleloop = default(Otri); + Otri loop = default(Otri); Otri oppotri = default(Otri); Osub opposubseg = default(Osub); Vertex triorg, tridest, triapex; @@ -157,18 +157,18 @@ namespace TriangleNet horrors = 0; // Run through the list of triangles, checking each one. - foreach (var t in mesh.triangles.Values) + foreach (var tri in mesh.triangles.Values) { - triangleloop.triangle = t; + loop.triangle = tri; // Check all three edges of the triangle. - for (triangleloop.orient = 0; triangleloop.orient < 3; - triangleloop.orient++) + for (loop.orient = 0; loop.orient < 3; + loop.orient++) { - triorg = triangleloop.Org(); - tridest = triangleloop.Dest(); - triapex = triangleloop.Apex(); - triangleloop.Sym(ref oppotri); + triorg = loop.Org(); + tridest = loop.Dest(); + triapex = loop.Apex(); + loop.Sym(ref oppotri); oppoapex = oppotri.Apex(); // Only test that the edge is locally Delaunay if there is an // adjoining triangle whose pointer is larger (to ensure that @@ -187,7 +187,7 @@ namespace TriangleNet { // If a subsegment separates the triangles, then the edge is // constrained, so no local Delaunay test should be done. - triangleloop.SegPivot(ref opposubseg); + loop.SegPivot(ref opposubseg); if (opposubseg.seg != Mesh.dummysub) { shouldbedelaunay = false; diff --git a/Triangle.NET/Triangle/Tools/BoundedVoronoi.cs b/Triangle.NET/Triangle/Tools/BoundedVoronoi.cs new file mode 100644 index 0000000..7301e91 --- /dev/null +++ b/Triangle.NET/Triangle/Tools/BoundedVoronoi.cs @@ -0,0 +1,551 @@ +// ----------------------------------------------------------------------- +// +// Triangle.NET code by Christian Woltering, http://home.edo.tu-dortmund.de/~woltering/index.html +// +// ----------------------------------------------------------------------- + +namespace TriangleNet.Tools +{ + using System; + using System.Collections.Generic; + using System.Linq; + using System.Text; + using TriangleNet.Data; + using TriangleNet.Geometry; + + /// + /// The Bounded Voronoi Diagram is the dual of a PSLG triangulation. + /// + /// + /// 2D Centroidal Voronoi Tessellations with Constraints, 2010, + /// Jane Tournois, Pierre Alliez and Olivier Devillers + /// + public class BoundedVoronoi + { + Mesh mesh; + Dictionary subsegMap; + List voronoi; + + /// + /// Initializes a new instance of the class. + /// + /// Mesh instance. + 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. + /// + public List Cells + { + get { return voronoi; } + } + + /// + /// Computes the bounded voronoi diagram. + /// + public void Generate() + { + voronoi.Clear(); + + TagBlindTriangles(); + + foreach (var v in mesh.vertices.Values) + { + if (v.type == VertexType.FreeVertex) + { + voronoi.Add(ConstructBvdCell(v)); + } + else + { + voronoi.Add(ConstructBoundaryBvdCell(v)); + } + } + } + + private void TagBlindTriangles() + { + int blinded = 0; + + Stack triangles; + subsegMap = new Dictionary(); + + Otri f = default(Otri); + Otri f0 = default(Otri); + Osub e = default(Osub); + Osub sub1 = default(Osub); + + // Tag all triangles non-blind + foreach (var t in mesh.triangles.Values) + { + // Use the infected flag for 'blinded' attribute. + t.infected = false; + } + + // for each constrained edge e of cdt do + foreach (var ss in mesh.subsegs.Values) + { + // Create a stack: triangles + triangles = new Stack(); + + // for both adjacent triangles fe to e tagged non-blind do + // Push fe into triangles + e.seg = ss; + e.orient = 0; + e.TriPivot(ref f); + + if (f.triangle != Mesh.dummytri && !f.triangle.infected) + { + triangles.Push(f.triangle); + } + + e.SymSelf(); + e.TriPivot(ref f); + + if (f.triangle != Mesh.dummytri && !f.triangle.infected) + { + triangles.Push(f.triangle); + } + + // while triangles is non-empty + while (triangles.Count > 0) + { + // Pop f from stack triangles + f.triangle = triangles.Pop(); + f.orient = 0; + + // if f is blinded by e (use P) then + if (TriangleIsBlinded(ref f, ref e)) + { + // Tag f as blinded by e + f.triangle.infected = true; + blinded++; + + // Store association triangle -> subseg + subsegMap.Add(f.triangle.hash, e.seg); + + // for each adjacent triangle f0 to f do + for (f.orient = 0; f.orient < 3; f.orient++) + { + f.Sym(ref f0); + + f0.SegPivot(ref sub1); + + // if f0 is finite and tagged non-blind & the common edge + // between f and f0 is unconstrained then + if (f0.triangle != Mesh.dummytri && !f0.triangle.infected && sub1.seg == Mesh.dummysub) + { + // Push f0 into triangles. + triangles.Push(f0.triangle); + } + } + } + } + } + + blinded = 0; + } + + private bool TriangleIsBlinded(ref Otri tri, ref Osub seg) + { + Point cc, pt = new Point(); + double xi = 0, eta = 0; + + Vertex torg = tri.Org(); + Vertex tdest = tri.Dest(); + Vertex tapex = tri.Apex(); + + cc = Primitives.FindCircumcenter(torg, tdest, tapex, ref xi, ref eta, false); + + if (SegmentsIntersect(ref seg, cc, torg, ref pt, true)) + { + return true; + } + + if (SegmentsIntersect(ref seg, cc, tdest, ref pt, true)) + { + return true; + } + + if (SegmentsIntersect(ref seg, cc, tapex, ref pt, true)) + { + return true; + } + + return false; + } + + private Point[] ConstructBvdCell(Vertex x) + { + 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; + + // Call P the polygon (cell) in construction + List P = new List(); + + // Call f_init a triangle incident to x + x.tri.Copy(ref f_init); + + if (f_init.Org() != x) + { + throw new Exception("ConstructBvdCell: inconsistent topology."); + } + + // Let f be initialized to f_init + f_init.Copy(ref f); + // Call f_next the next triangle counterclockwise around x + f_init.Onext(ref f_next); + + // 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, false); + + torg = f_next.Org(); + tdest = f_next.Dest(); + tapex = f_next.Apex(); + cc_f_next = Primitives.FindCircumcenter(torg, tdest, tapex, ref xi, ref eta, false); + + // if f is tagged non-blind then + if (!f.triangle.infected) + { + // Insert the circumcenter of f into P + P.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)) + { + P.Add(p); + } + } + } + else + { + // Call Sf the constrained edge blinding f + sf.seg = subsegMap[f.triangle.hash]; + + // 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)) + { + P.Add(p); + } + } + else + { + // Call Sf_next the constrained edge blinding f_next + sfn.seg = subsegMap[f_next.triangle.hash]; + + // 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)) + { + P.Add(p); + } + + p = new Point(); + if (SegmentsIntersect(ref sfn, cc_f, cc_f_next, ref p, true)) + { + P.Add(p); + } + } + } + } + + // f <- f_next + f_next.Copy(ref f); + + // Call f_next the next triangle counterclockwise around x + f_next.OnextSelf(); + } + while (!f.Equal(f_init)); + + // Output: Bounded Voronoi cell of x in counterclockwise order. + return P.ToArray(); + } + + private Point[] ConstructBoundaryBvdCell(Vertex x) + { + Otri f = default(Otri); + Otri f_init = default(Otri); + Otri f_next = default(Otri); + Otri f_prev = 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; + + // Call P the polygon (cell) in construction + List P = new List(); + + // Call f_init a triangle incident to x + x.tri.Copy(ref f_init); + + if (f_init.Org() != x) + { + throw new Exception("ConstructBoundaryBvdCell: inconsistent topology."); + } + // Let f be initialized to f_init + f_init.Copy(ref f); + // Call f_next the next triangle counterclockwise around x + f_init.Onext(ref f_next); + + f_init.Oprev(ref f_prev); + + // Is the border to the left? + if (f_prev.triangle != Mesh.dummytri) + { + // Go clockwise until we reach the border (or the initial triangle) + while (f_prev.triangle != Mesh.dummytri && !f_prev.Equal(f_init)) + { + f_prev.Copy(ref f); + f_prev.OprevSelf(); + } + + f.Copy(ref f_init); + f.Onext(ref f_next); + } + + // Add vertex on border + P.Add(x); + + // 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)); + + // 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, false); + + if (f_next.triangle == Mesh.dummytri) + { + if (!f.triangle.infected) + { + // Add last circumcenter + P.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)); + + 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, false); + + // if f is tagged non-blind then + if (!f.triangle.infected) + { + // Insert the circumcenter of f into P + P.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)) + { + P.Add(p); + } + } + } + else + { + // Call Sf the constrained edge blinding f + sf.seg = subsegMap[f.triangle.hash]; + + // if f_next is tagged non-blind then + if (!f_next.triangle.infected) + { + tdest = f.Dest(); + tapex = f.Apex(); + + // 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)) + { + P.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)) + { + P.Add(p); + } + } + else + { + // Call Sf_next the constrained edge blinding f_next + sfn.seg = subsegMap[f_next.triangle.hash]; + + // 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)) + { + P.Add(p); + } + + p = new Point(); + if (SegmentsIntersect(ref sfn, cc_f, cc_f_next, ref p, true)) + { + P.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)) + { + P.Add(p); + } + } + } + } + + // f <- f_next + f_next.Copy(ref f); + + // Call f_next the next triangle counterclockwise around x + f_next.OnextSelf(); + } + while (!f.Equal(f_init)); + + // Output: Bounded Voronoi cell of x in counterclockwise order. + return P.ToArray(); + } + + /// + /// Determines the intersection point of the line segment defined by points A and B with the + /// line segment defined by points C and D. + /// + /// The first segment AB. + /// Endpoint C of second segment. + /// Endpoint D of second segment. + /// Reference to the intersection point. + /// If false, pa and pb represent a line. + /// Returns true if the intersection point was found, and stores that point in X,Y. + /// 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) + { + Vertex s1 = seg.Org(); + Vertex s2 = seg.Dest(); + + 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 distAB, theCos, theSin, newX, ABpos; + + // Fail if either line segment is zero-length. + if (Ax == Bx && Ay == By || Cx == Dx && Cy == Dy) return false; + + // Fail if the segments share an end-point. + if (Ax == Cx && Ay == Cy || Bx == Cx && By == Cy + || Ax == Dx && Ay == Dy || Bx == Dx && By == Dy) + { + return false; + } + + // (1) Translate the system so that point A is on the origin. + Bx -= Ax; By -= Ay; + Cx -= Ax; Cy -= Ay; + Dx -= Ax; Dy -= Ay; + + // Discover the length of segment A-B. + distAB = Math.Sqrt(Bx * Bx + By * By); + + // (2) Rotate the system so that point B is on the positive X axis. + theCos = Bx / distAB; + theSin = By / distAB; + newX = Cx * theCos + Cy * theSin; + Cy = Cy * theCos - Cx * theSin; Cx = newX; + newX = Dx * theCos + Dy * theSin; + Dy = Dy * theCos - Dx * theSin; Dx = newX; + + // Fail if segment C-D doesn't cross line A-B. + if (Cy < 0 && Dy < 0 || Cy >= 0 && Dy >= 0 && strictIntersect) return false; + + // (3) Discover the position of the intersection point along line A-B. + ABpos = Dx + (Cx - Dx) * Dy / (Dy - Cy); + + // Fail if segment C-D crosses line A-B outside of segment A-B. + 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; + + // Success. + return true; + } + } +} diff --git a/Triangle.NET/Triangle/Tools/QualityMeasure.cs b/Triangle.NET/Triangle/Tools/QualityMeasure.cs new file mode 100644 index 0000000..40b664d --- /dev/null +++ b/Triangle.NET/Triangle/Tools/QualityMeasure.cs @@ -0,0 +1,541 @@ +// ----------------------------------------------------------------------- +// +// TODO: Update copyright text. +// +// ----------------------------------------------------------------------- + +namespace TriangleNet.Tools +{ + using System; + using System.Collections.Generic; + using System.Linq; + using System.Text; + using TriangleNet.Geometry; + + /// + /// Provides mesh quality information. + /// + /// + /// Given a triangle abc with points A (ax, ay), B (bx, by), C (cx, cy). + /// + /// The side lengths are given as + /// a = sqrt((cx - bx)^2 + (cy - by)^2) -- side BC opposite of A + /// b = sqrt((cx - ax)^2 + (cy - ay)^2) -- side CA opposite of B + /// c = sqrt((ax - bx)^2 + (ay - by)^2) -- side AB opposite of C + /// + /// The angles are given as + /// ang_a = acos((b^2 + c^2 - a^2) / (2 * b * c)) -- angle at A + /// ang_b = acos((c^2 + a^2 - b^2) / (2 * c * a)) -- angle at B + /// ang_c = acos((a^2 + b^2 - c^2) / (2 * a * b)) -- angle at C + /// + /// The semiperimeter is given as + /// s = (a + b + c) / 2 + /// + /// The area is given as + /// D = abs(ax * (by - cy) + bx * (cy - ay) + cx * (ay - by)) / 2 + /// = sqrt(s * (s - a) * (s - b) * (s - c)) + /// + /// The inradius is given as + /// r = D / s + /// + /// The circumradius is given as + /// R = a * b * c / (4 * D) + /// + /// The altitudes are given as + /// alt_a = 2 * D / a -- altitude above side a + /// alt_b = 2 * D / b -- altitude above side b + /// alt_c = 2 * D / c -- altitude above side c + /// + /// The aspect ratio may be given as the ratio of the longest to the + /// shortest edge or, more commonly as the ratio of the circumradius + /// to twice the inradius + /// ar = R / (2 * r) + /// = a * b * c / (8 * (s - a) * (s - b) * (s - c)) + /// = a * b * c / ((b + c - a) * (c + a - b) * (a + b - c)) + /// + public class QualityMeasure + { + AreaMeasure areaMeasure; + AlphaMeasure alphaMeasure; + Q_Measure qMeasure; + + Mesh mesh; + + public QualityMeasure() + { + areaMeasure = new AreaMeasure(); + alphaMeasure = new AlphaMeasure(); + qMeasure = new Q_Measure(); + } + + #region Public properties + + /// + /// Minimum triangle area. + /// + public double AreaMinimum + { + get { return areaMeasure.area_min; } + } + + /// + /// Maximum triangle area. + /// + public double AreaMaximum + { + get { return areaMeasure.area_max; } + } + + /// + /// Ratio of maximum and minimum triangle area. + /// + public double AreaRatio + { + get { return areaMeasure.area_max / areaMeasure.area_min; } + } + + /// + /// Smallest angle. + /// + public double AlphaMinimum + { + get { return alphaMeasure.alpha_min; } + } + + /// + /// Maximum smallest angle. + /// + public double AlphaMaximum + { + get { return alphaMeasure.alpha_max; } + } + + /// + /// Average angle. + /// + public double AlphaAverage + { + get { return alphaMeasure.alpha_ave; } + } + + /// + /// Average angle weighted by area. + /// + public double AlphaArea + { + get { return alphaMeasure.alpha_area; } + } + + /// + /// Smallest aspect ratio. + /// + public double Q_Minimum + { + get { return qMeasure.q_min; } + } + + /// + /// Largest aspect ratio. + /// + public double Q_Maximum + { + get { return qMeasure.q_max; } + } + + /// + /// Average aspect ratio. + /// + public double Q_Average + { + get { return qMeasure.q_ave; } + } + + /// + /// Average aspect ratio weighted by area. + /// + public double Q_Area + { + get { return qMeasure.q_area; } + } + + #endregion + + public void Update(Mesh mesh) + { + this.mesh = mesh; + + // Reset all measures. + areaMeasure.Reset(); + alphaMeasure.Reset(); + qMeasure.Reset(); + + Compute(); + } + + private void Compute() + { + Point a, b, c; + double ab, bc, ca; + double lx, ly; + double area; + + int n = 0; + + foreach (var tri in mesh.triangles.Values) + { + n++; + + a = tri.vertices[0]; + b = tri.vertices[1]; + c = tri.vertices[2]; + + lx = a.x - b.x; + ly = a.y - b.y; + ab = Math.Sqrt(lx * lx + ly * ly); + lx = b.x - c.x; + ly = b.y - c.y; + bc = Math.Sqrt(lx * lx + ly * ly); + lx = c.x - a.x; + ly = c.y - a.y; + ca = Math.Sqrt(lx * lx + ly * ly); + + area = areaMeasure.Measure(a, b, c); + alphaMeasure.Measure(ab, bc, ca, area); + qMeasure.Measure(ab, bc, ca, area); + } + + // Normalize measures + alphaMeasure.Normalize(n, areaMeasure.area_total); + qMeasure.Normalize(n, areaMeasure.area_total); + } + + /// + /// Determines the bandwidth of the coefficient matrix. + /// + /// Bandwidth of the coefficient matrix. + /// + /// The quantity computed here is the "geometric" bandwidth determined + /// by the finite element mesh alone. + /// + /// If a single finite element variable is associated with each node + /// of the mesh, and if the nodes and variables are numbered in the + /// same way, then the geometric bandwidth is the same as the bandwidth + /// of a typical finite element matrix. + /// + /// The bandwidth M is defined in terms of the lower and upper bandwidths: + /// + /// M = ML + 1 + MU + /// + /// where + /// + /// ML = maximum distance from any diagonal entry to a nonzero + /// entry in the same row, but earlier column, + /// + /// MU = maximum distance from any diagonal entry to a nonzero + /// entry in the same row, but later column. + /// + /// Because the finite element node adjacency relationship is symmetric, + /// we are guaranteed that ML = MU. + /// + public int Bandwidth() + { + // Lower and upper bandwidth of the matrix + int ml = 0, mu = 0; + + int gi, gj; + + foreach (var tri in mesh.triangles.Values) + { + for (int j = 0; j < 3; j++) + { + gi = tri[j]; + + for (int k = 0; k < 3; k++) + { + gj = tri[k]; + + mu = Math.Max(mu, gj - gi); + ml = Math.Max(ml, gi - gj); + } + } + } + + return ml + 1 + mu; + } + + class AreaMeasure + { + // Minimum area + public double area_min = double.MaxValue; + // Maximum area + public double area_max = -double.MaxValue; + // Total area of geometry + public double area_total = 0; + // Nmber of triangles with zero area + public int area_zero = 0; + + /// + /// Reset all values. + /// + public void Reset() + { + area_min = double.MaxValue; + area_max = -double.MaxValue; + area_total = 0; + area_zero = 0; + } + + /// + /// Compute the area of given triangle. + /// + /// Triangle corner a. + /// Triangle corner b. + /// Triangle corner c. + /// Triangle area. + public double Measure(Point a, Point b, Point c) + { + double area = 0.5 * Math.Abs(a.x * (b.y - c.y) + b.x * (c.y - a.y) + c.x * (a.y - b.y)); + + area_min = Math.Min(area_min, area); + area_max = Math.Max(area_max, area); + area_total += area; + + if (area == 0.0) + { + area_zero = area_zero + 1; + } + + return area; + } + } + + /// + /// The alpha measure determines the triangulated pointset quality. + /// + /// + /// The alpha measure evaluates the uniformity of the shapes of the triangles + /// defined by a triangulated pointset. + /// + /// We compute the minimum angle among all the triangles in the triangulated + /// dataset and divide by the maximum possible value (which, in degrees, + /// is 60). The best possible value is 1, and the worst 0. A good + /// triangulation should have an alpha score close to 1. + /// + class AlphaMeasure + { + // Minimum value over all triangles + public double alpha_min; + // Maximum value over all triangles + public double alpha_max; + // Value averaged over all triangles + public double alpha_ave; + // Value averaged over all triangles and weighted by area + public double alpha_area; + + /// + /// Reset all values. + /// + public void Reset() + { + alpha_min = double.MaxValue; + alpha_max = -double.MaxValue; + alpha_ave = 0; + alpha_area = 0; + } + + double acos(double c) + { + if (c <= -1.0) + { + return Math.PI; + } + else if (1.0 <= c) + { + return 0.0; + } + else + { + return Math.Acos(c); + } + } + + /// + /// Compute q value of given triangle. + /// + /// Side length ab. + /// Side length bc. + /// Side length ca. + /// Triangle area. + /// + public double Measure(double ab, double bc, double ca, double area) + { + double alpha = double.MaxValue; + + double ab2 = ab * ab; + double bc2 = bc * bc; + double ca2 = ca * ca; + + double a_angle; + double b_angle; + double c_angle; + + // Take care of a ridiculous special case. + if (ab == 0.0 && bc == 0.0 && ca == 0.0) + { + a_angle = 2.0 * Math.PI / 3.0; + b_angle = 2.0 * Math.PI / 3.0; + c_angle = 2.0 * Math.PI / 3.0; + } + else + { + if (ca == 0.0 || ab == 0.0) + { + a_angle = Math.PI; + } + else + { + a_angle = acos((ca2 + ab2 - bc2) / (2.0 * ca * ab)); + } + + if (ab == 0.0 || bc == 0.0) + { + b_angle = Math.PI; + } + else + { + b_angle = acos((ab2 + bc2 - ca2) / (2.0 * ab * bc)); + } + + if (bc == 0.0 || ca == 0.0) + { + c_angle = Math.PI; + } + else + { + c_angle = acos((bc2 + ca2 - ab2) / (2.0 * bc * ca)); + } + } + + alpha = Math.Min(alpha, a_angle); + alpha = Math.Min(alpha, b_angle); + alpha = Math.Min(alpha, c_angle); + + // Normalize angle from [0,pi/3] radians into qualities in [0,1]. + alpha = alpha * 3.0 / Math.PI; + + alpha_ave += alpha; + alpha_area += area * alpha; + + alpha_min = Math.Min(alpha, alpha_min); + alpha_max = Math.Max(alpha, alpha_max); + + return alpha; + } + + /// + /// Normalize values. + /// + public void Normalize(int n, double area_total) + { + if (n > 0) + { + alpha_ave /= n; + } + else + { + alpha_ave = 0.0; + } + + if (0.0 < area_total) + { + alpha_area /= area_total; + } + else + { + alpha_area = 0.0; + } + } + } + + /// + /// The Q measure determines the triangulated pointset quality. + /// + /// + /// The Q measure evaluates the uniformity of the shapes of the triangles + /// defined by a triangulated pointset. It uses the aspect ratio + /// + /// 2 * (incircle radius) / (circumcircle radius) + /// + /// In an ideally regular mesh, all triangles would have the same + /// equilateral shape, for which Q = 1. A good mesh would have + /// 0.5 < Q. + /// + class Q_Measure + { + // Minimum value over all triangles + public double q_min; + // Maximum value over all triangles + public double q_max; + // Average value + public double q_ave; + // Average value weighted by the area of each triangle + public double q_area; + + /// + /// Reset all values. + /// + public void Reset() + { + q_min = double.MaxValue; + q_max = -double.MaxValue; + q_ave = 0; + q_area = 0; + } + + /// + /// Compute q value of given triangle. + /// + /// Side length ab. + /// Side length bc. + /// Side length ca. + /// Triangle area. + /// + public double Measure(double ab, double bc, double ca, double area) + { + double q = (bc + ca - ab) * (ca + ab - bc) * (ab + bc - ca) / (ab * bc * ca); + + q_min = Math.Min(q_min, q); + q_max = Math.Max(q_max, q); + + q_ave += q; + q_area += q * area; + + return q; + } + + /// + /// Normalize values. + /// + public void Normalize(int n, double area_total) + { + if (n > 0) + { + q_ave /= n; + } + else + { + q_ave = 0.0; + } + + if (area_total > 0.0) + { + q_area /= area_total; + } + else + { + q_area = 0.0; + } + } + } + } +} diff --git a/Triangle.NET/Triangle/IO/DataWriter.cs b/Triangle.NET/Triangle/Tools/Voronoi.cs similarity index 62% rename from Triangle.NET/Triangle/IO/DataWriter.cs rename to Triangle.NET/Triangle/Tools/Voronoi.cs index 1d5dcca..46fd671 100644 --- a/Triangle.NET/Triangle/IO/DataWriter.cs +++ b/Triangle.NET/Triangle/Tools/Voronoi.cs @@ -1,25 +1,60 @@ // ----------------------------------------------------------------------- -// +// // Original Triangle code by Jonathan Richard Shewchuk, http://www.cs.cmu.edu/~quake/triangle.html -// Triangle.NET code by Christian Woltering, http://home.edo.tu-dortmund.de/~woltering/triangle/ +// Triangle.NET code by Christian Woltering, http://home.edo.tu-dortmund.de/~woltering/index.html // // ----------------------------------------------------------------------- -namespace TriangleNet.IO -{ - using System; - using System.IO; - using System.Collections.Generic; - using System.Globalization; - using TriangleNet.Data; - using TriangleNet.Geometry; +using TriangleNet.Data; +using TriangleNet.Geometry; +namespace TriangleNet.Tools +{ /// - /// Generates a mesh representaion using arrays. + /// The Voronoi Diagram is the dual of a pointset triangulation. /// - public static class DataWriter + public class Voronoi { - #region Library + Mesh mesh; + Point[] points; + Point[] directions; // Stores the direction for infinite voronoi edges + Edge[] edges; + + /// + /// Initializes a new instance of the class. + /// + /// + /// + /// Be sure MakeVertexMap has been called (should always be the case). + /// + public Voronoi(Mesh mesh) + { + this.mesh = mesh; + } + + /// + /// Gets the list of Voronoi vertices. + /// + public Point[] Points + { + get { return points; } + } + + /// + /// Gets the directions for infinite Voronoi edges. + /// + public Point[] Directions + { + get { return directions; } + } + + /// + /// Gets the list of Voronoi edges. + /// + public Edge[] Edges + { + get { return edges; } + } /// /// Gets the Voronoi diagram as raw output data. @@ -31,42 +66,21 @@ namespace TriangleNet.IO /// Hence, the Voronoi vertices are listed by traversing the Delaunay /// triangles, and the Voronoi edges are listed by traversing the Delaunay /// edges. - /// - /// WARNING: In order to assign numbers to the Voronoi vertices, this - /// procedure messes up the subsegments or the extra nodes of every - /// element. Hence, you should call this procedure last. - public static VoronoiData WriteVoronoi(Mesh mesh) + /// + public void Generate() { - VoronoiData data = new VoronoiData(); - Otri tri = default(Otri), trisym = default(Otri); Vertex torg, tdest, tapex; Point circumcenter; double xi = 0, eta = 0; int p1, p2; - int i = 0; - - // Copy input points (actually not part of the voronoi diagram) - data.InputPoints = new Vertex[mesh.vertices.Count]; - - foreach (var item in mesh.vertices.Values) - { - if (item.type != VertexType.UndeadVertex) - { - data.InputPoints[i] = item; - i++; - } - } - // Allocate memory for Voronoi vertices. - data.Points = new Vertex[mesh.triangles.Count]; - - int index = 0; + this.points = new Point[mesh.triangles.Count]; tri.orient = 0; - i = 0; + int i = 0; foreach (var item in mesh.triangles.Values) { tri.triangle = item; @@ -76,19 +90,19 @@ namespace TriangleNet.IO circumcenter = Primitives.FindCircumcenter(torg, tdest, tapex, ref xi, ref eta, false); // X and y coordinates. - data.Points[i] = new Point(circumcenter.x, circumcenter.y); + this.points[i] = new Point(circumcenter.x, circumcenter.y); // Update element id tri.triangle.id = i++; } // Allocate memory for output Voronoi edges. - data.Edges = new Edge[mesh.edges]; + this.edges = new Edge[mesh.edges]; // Allocate memory for output Voronoi norms. - data.Directions = new double[mesh.edges][]; + this.directions = new Point[mesh.edges]; - index = 0; + 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 @@ -113,8 +127,8 @@ namespace TriangleNet.IO tdest = tri.Dest(); // Copy an infinite ray. Index of one endpoint, and -1. - data.Edges[index] = new Edge(p1, -1); - data.Directions[index] = new double[] { tdest[1] - torg[1], torg[0] - tdest[0] }; + this.edges[i] = new Edge(p1, -1); + this.directions[i] = new Point(tdest.y - torg.y, torg.x - tdest.x); } else { @@ -122,18 +136,14 @@ namespace TriangleNet.IO p2 = trisym.triangle.id; // Finite edge. Write indices of two endpoints. - data.Edges[index] = new Edge(p1, p2); - data.Directions[index] = new double[] { 0, 0 }; + this.edges[i] = new Edge(p1, p2); + this.directions[i] = new Point(0, 0); } - index++; + i++; } } } - - return data; } - - #endregion } } diff --git a/Triangle.NET/Triangle/Triangle.csproj b/Triangle.NET/Triangle/Triangle.csproj index fe17888..0a09499 100644 --- a/Triangle.NET/Triangle/Triangle.csproj +++ b/Triangle.NET/Triangle/Triangle.csproj @@ -64,8 +64,6 @@ - - @@ -80,8 +78,11 @@ + + +