diff --git a/Triangle.NET/TestApp/Controls/RendererControl.cs b/Triangle.NET/TestApp/Controls/RendererControl.cs index 955b9a9..c15a8b6 100644 --- a/Triangle.NET/TestApp/Controls/RendererControl.cs +++ b/Triangle.NET/TestApp/Controls/RendererControl.cs @@ -73,7 +73,8 @@ namespace MeshExplorer.Controls timer = new Timer(); timer.Interval = 3000; - timer.Tick += (sender, e) => { + timer.Tick += (sender, e) => + { timer.Stop(); coordinate = String.Empty; this.Invalidate(); @@ -90,13 +91,6 @@ namespace MeshExplorer.Controls this.Invalidate(); } - public void ShowQuality(int measure) - { - //Tuple[] q = TriangleQuality.Measure(data, measure); - - //this.RenderQualities(q); - } - public void SetData(InputGeometry mesh) { data.SetData(mesh); @@ -202,66 +196,20 @@ namespace MeshExplorer.Controls g.SmoothingMode = SmoothingMode.AntiAlias; - if (meshRenderer != null) + if (voronoiRenderer != null && this.showVoronoi) + { + meshRenderer.RenderMesh(g, zoom, renderColors); + voronoiRenderer.Render(g, zoom, renderColors); + meshRenderer.RenderGeometry(g, zoom, renderColors); + } + else if (meshRenderer != null) { meshRenderer.Render(g, zoom, renderColors); } - if (voronoiRenderer != null && this.showVoronoi) - { - voronoiRenderer.Render(g, zoom, renderColors); - } - this.Invalidate(); } - /* - private void RenderQualitiesX(Tuple[] q) - { - PointF[] p = new PointF[3]; - PointF[] pts = data.Points; - - int[] tri; - - Brush q1 = new SolidBrush(Color.FromArgb(50, Color.Orange)); - Brush q2 = new SolidBrush(Color.FromArgb(50, Color.Red)); - - Graphics g = buffer.Graphics; - - g.SmoothingMode = SmoothingMode.AntiAlias; - - int n = q.Length; - - for (int i = 0; i < n; i++) - { - tri = data.Triangles[q[i].Item1]; - - if (zoom.ViewportContains(pts[tri[0]]) || - zoom.ViewportContains(pts[tri[1]]) || - zoom.ViewportContains(pts[tri[2]])) - { - p[0] = zoom.WorldToScreen(pts[tri[0]]); - p[1] = zoom.WorldToScreen(pts[tri[1]]); - p[2] = zoom.WorldToScreen(pts[tri[2]]); - - // Fill - g.FillPolygon(q[i].Item2 > 1 ? q2 : q1, p); - - // Outline - g.DrawLine(lines, p[0], p[1]); - g.DrawLine(lines, p[1], p[2]); - g.DrawLine(lines, p[2], p[0]); - - // Points - g.FillEllipse(Brushes.Green, p[0].X - 1.5f, p[0].Y - 1.5f, 3, 3); - g.FillEllipse(Brushes.Green, p[1].X - 1.5f, p[1].Y - 1.5f, 3, 3); - g.FillEllipse(Brushes.Green, p[2].X - 1.5f, p[2].Y - 1.5f, 3, 3); - } - } - - this.Invalidate(); - }*/ - #region Control overrides protected override void OnPaint(PaintEventArgs pe) @@ -297,8 +245,8 @@ namespace MeshExplorer.Controls timer.Stop(); PointF c = zoom.ScreenToWorld((float)e.X / this.Width, (float)e.Y / this.Height); - coordinate = String.Format("X:{0} Y:{1}", - c.X.ToString(Util.Nfi), + coordinate = String.Format("X:{0} Y:{1}", + c.X.ToString(Util.Nfi), c.Y.ToString(Util.Nfi)); this.Invalidate(); diff --git a/Triangle.NET/TestApp/Rendering/MeshRenderer.cs b/Triangle.NET/TestApp/Rendering/MeshRenderer.cs index 4b2dd93..217a2c0 100644 --- a/Triangle.NET/TestApp/Rendering/MeshRenderer.cs +++ b/Triangle.NET/TestApp/Rendering/MeshRenderer.cs @@ -146,5 +146,42 @@ namespace MeshExplorer.Rendering this.RenderPoints(g); } } + + /// + /// Renders only the mesh edges (no points or segments). + /// + public void RenderMesh(Graphics g, Zoom zoom, RenderColors renderColors) + { + this.renderColors = renderColors; + this.zoom = zoom; + + if (data.Edges != null) + { + this.RenderEdges(g); + } + else if (data.Triangles != null) + { + this.RenderTriangles(g); + } + } + + /// + /// Renders only points and segments (no mesh triangles). + /// + public void RenderGeometry(Graphics g, Zoom zoom, RenderColors renderColors) + { + this.renderColors = renderColors; + this.zoom = zoom; + + if (data.Segments != null) + { + this.RenderSegments(g); + } + + if (data.Points != null) + { + this.RenderPoints(g); + } + } } } diff --git a/Triangle.NET/Triangle/Algorithm/Dwyer.cs b/Triangle.NET/Triangle/Algorithm/Dwyer.cs index fb2415e..ecd57c4 100644 --- a/Triangle.NET/Triangle/Algorithm/Dwyer.cs +++ b/Triangle.NET/Triangle/Algorithm/Dwyer.cs @@ -781,8 +781,11 @@ namespace TriangleNet.Algorithm Otri dissolveedge = default(Otri); Otri deadtriangle = default(Otri); Vertex markorg; + int hullsize; + bool noPoly = !mesh.behavior.Poly; + // Find an edge on the convex hull to start point location from. startghost.Lprev(ref searchedge); searchedge.SymSelf(); @@ -799,7 +802,7 @@ namespace TriangleNet.Algorithm // If no PSLG is involved, set the boundary markers of all the vertices // on the convex hull. If a PSLG is used, this step is done later. - if (!Behavior.Poly) + if (noPoly) { // Watch out for the case where all the input vertices are collinear. if (dissolveedge.triangle != Mesh.dummytri) diff --git a/Triangle.NET/Triangle/Algorithm/Incremental.cs b/Triangle.NET/Triangle/Algorithm/Incremental.cs index d1c4850..65f3f48 100644 --- a/Triangle.NET/Triangle/Algorithm/Incremental.cs +++ b/Triangle.NET/Triangle/Algorithm/Incremental.cs @@ -75,6 +75,8 @@ namespace TriangleNet.Algorithm Vertex markorg; int hullsize; + bool noPoly = !mesh.behavior.Poly; + // Find a boundary triangle. nextedge.triangle = Mesh.dummytri; nextedge.orient = 0; @@ -110,7 +112,7 @@ namespace TriangleNet.Algorithm dissolveedge.SymSelf(); // If not using a PSLG, the vertices should be marked now. // (If using a PSLG, markhull() will do the job.) - if (!Behavior.Poly) + if (noPoly) { // Be careful! One must check for the case where all the input // vertices are collinear, and thus all the triangles are part of diff --git a/Triangle.NET/Triangle/Algorithm/SweepLine.cs b/Triangle.NET/Triangle/Algorithm/SweepLine.cs index 5252085..c106712 100644 --- a/Triangle.NET/Triangle/Algorithm/SweepLine.cs +++ b/Triangle.NET/Triangle/Algorithm/SweepLine.cs @@ -480,6 +480,8 @@ namespace TriangleNet.Algorithm Vertex markorg; int hullsize; + bool noPoly = !mesh.behavior.Poly; + // Find an edge on the convex hull to start point location from. startghost.Lprev(ref searchedge); searchedge.SymSelf(); @@ -496,7 +498,7 @@ namespace TriangleNet.Algorithm // If no PSLG is involved, set the boundary markers of all the vertices // on the convex hull. If a PSLG is used, this step is done later. - if (!Behavior.Poly) + if (noPoly) { // Watch out for the case where all the input vertices are collinear. if (dissolveedge.triangle != Mesh.dummytri) diff --git a/Triangle.NET/Triangle/Behavior.cs b/Triangle.NET/Triangle/Behavior.cs index df358e5..24a1cb8 100644 --- a/Triangle.NET/Triangle/Behavior.cs +++ b/Triangle.NET/Triangle/Behavior.cs @@ -5,115 +5,19 @@ // // ----------------------------------------------------------------------- -using System; namespace TriangleNet { - // TODO: Make Behavior non-static and an instance member of mesh class. + using System; /// /// Controls the behavior of the meshing software. /// - static class Behavior + class Behavior { - /// - /// Input is a Planar Straight Line Graph. - /// - public static bool Poly { get; set; } - /// - /// Quality mesh generation. - /// - public static bool Quality { get; set; } - /// - /// Apply a maximum triangle area constraint. - /// - public static bool VarArea { get; set; } - /// - /// Apply a maximum triangle area constraint. - /// - public static bool FixedArea { get; set; } - /// - /// Apply a user-defined triangle constraint. - /// - public static bool Usertest { get; set; } - /// - /// Apply attributes to identify triangles in certain regions. - /// - public static bool RegionAttrib { get; set; } - /// - /// Enclose the convex hull with segments. - /// - public static bool Convex { get; set; } - /// - /// Jettison unused vertices from output. - /// - public static bool Jettison { get; set; } - /// - /// Compute boundary information. - /// - public static bool UseBoundaryMarkers { get; set; } - /// - /// Ignores holes in polygons. - /// - public static bool NoHoles { get; set; } - /// - /// No exact arithmetic. - /// - public static bool NoExact { get; set; } - /// - /// Conforming Delaunay (all triangles are truly Delaunay). - /// - public static bool ConformDel { get; set; } - /// - /// Algorithm to use for triangulation. - /// - public static TriangulationAlgorithm Algorithm { get; set; } - /// - /// Log detailed information. - /// - public static bool Verbose { get; set; } - /// - /// Use segments (should not be set manually) - /// - public static bool UseSegments { get; set; } // TODO: internal set - - /// - /// Suppresses boundary segment splitting. - /// - public static int NoBisect { get; set; } // <- int ! - /// - /// Use maximum number of added Steiner points. - /// - public static int Steiner { get; set; } - - /// - /// Minimum angle constraint. - /// - public static double MinAngle { get; set; } - /// - /// (should not be set manually) - /// - public static double GoodAngle { get; set; } - /// - /// (should not be set manually) - /// - public static double Offconstant { get; set; } - /// - /// Maximum area constraint. - /// - public static double MaxArea { get; set; } - /// - /// Maximum angle constraint. - /// - public static double MaxAngle { get; set; } - /// - /// (should not be set manually) - /// - public static double MaxGoodAngle { get; set; } - /// /// Load behavior defaults. /// - public static void Init() + public Behavior() { Poly = false; Quality = false; @@ -125,10 +29,8 @@ namespace TriangleNet Jettison = false; UseBoundaryMarkers = true; NoHoles = false; - NoExact = false; ConformDel = false; Algorithm = TriangulationAlgorithm.Dwyer; - Verbose = true; UseSegments = true; NoBisect = 0; @@ -140,6 +42,114 @@ namespace TriangleNet MaxGoodAngle = 0.0; MaxArea = -1.0; Offconstant = 0.0; + + Verbose = true; + NoExact = false; } + + #region Static properties + + /// + /// No exact arithmetic. + /// + public static bool NoExact { get; set; } + + /// + /// Log detailed information. + /// + public static bool Verbose { get; set; } + + #endregion + + #region Public properties + + /// + /// Input is a Planar Straight Line Graph. + /// + public bool Poly { get; set; } + /// + /// Quality mesh generation. + /// + public bool Quality { get; set; } + /// + /// Apply a maximum triangle area constraint. + /// + public bool VarArea { get; set; } + /// + /// Apply a maximum triangle area constraint. + /// + public bool FixedArea { get; set; } + /// + /// Apply a user-defined triangle constraint. + /// + public bool Usertest { get; set; } + /// + /// Apply attributes to identify triangles in certain regions. + /// + public bool RegionAttrib { get; set; } + /// + /// Enclose the convex hull with segments. + /// + public bool Convex { get; set; } + /// + /// Jettison unused vertices from output. + /// + public bool Jettison { get; set; } + /// + /// Compute boundary information. + /// + public bool UseBoundaryMarkers { get; set; } + /// + /// Ignores holes in polygons. + /// + public bool NoHoles { get; set; } + /// + /// Conforming Delaunay (all triangles are truly Delaunay). + /// + public bool ConformDel { get; set; } + /// + /// Algorithm to use for triangulation. + /// + public TriangulationAlgorithm Algorithm { get; set; } + /// + /// Use segments (should not be set manually) + /// + public bool UseSegments { get; set; } // TODO: internal set + + /// + /// Suppresses boundary segment splitting. + /// + public int NoBisect { get; set; } // <- int ! + /// + /// Use maximum number of added Steiner points. + /// + public int Steiner { get; set; } + + /// + /// Minimum angle constraint. + /// + public double MinAngle { get; set; } + /// + /// (should not be set manually) + /// + public double GoodAngle { get; set; } + /// + /// (should not be set manually) + /// + public double Offconstant { get; set; } + /// + /// Maximum area constraint. + /// + public double MaxArea { get; set; } + /// + /// Maximum angle constraint. + /// + public double MaxAngle { get; set; } + /// + /// (should not be set manually) + /// + public double MaxGoodAngle { get; set; } + + #endregion } } diff --git a/Triangle.NET/Triangle/Carver.cs b/Triangle.NET/Triangle/Carver.cs index 5883194..4496f3e 100644 --- a/Triangle.NET/Triangle/Carver.cs +++ b/Triangle.NET/Triangle/Carver.cs @@ -19,7 +19,6 @@ namespace TriangleNet { Mesh mesh; - public Carver(Mesh mesh) { this.mesh = mesh; @@ -311,6 +310,8 @@ namespace TriangleNet Otri neighbor = default(Otri); Osub neighborsubseg = default(Osub); + Behavior behavior = mesh.behavior; + // Loop through all the infected triangles, spreading the attribute // and/or area constraint to their neighbors, then to their neighbors' // neighbors. @@ -325,12 +326,12 @@ namespace TriangleNet // adjacent subsegments. // TODO: Not true in the C# version (so we could skip this). testtri.Uninfect(); - if (Behavior.RegionAttrib) + if (behavior.RegionAttrib) { // Set an attribute (Note: the attributes array was resized before). testtri.triangle.attributes[mesh.eextras] = attribute; } - if (Behavior.VarArea) + if (behavior.VarArea) { // Set an area constraint. testtri.triangle.area = area; @@ -383,14 +384,14 @@ namespace TriangleNet Otri[] regionTris = null; - if (!Behavior.Convex) + if (!mesh.behavior.Convex) { // Mark as infected any unprotected triangles on the boundary. // This is one way by which concavities are created. InfectHull(); } - if (!Behavior.NoHoles) + if (!mesh.behavior.NoHoles) { // Infect each triangle in which a hole lies. foreach (var hole in mesh.holes) @@ -475,7 +476,7 @@ namespace TriangleNet if (regionTris != null) { - if (Behavior.RegionAttrib) + if (mesh.behavior.RegionAttrib) { // Make the triangle's attributes larger. double[] attributes = new double[mesh.eextras + 1]; @@ -508,7 +509,7 @@ namespace TriangleNet } } - if (Behavior.RegionAttrib) + if (mesh.behavior.RegionAttrib) { // Note the fact that each triangle has an additional attribute. mesh.eextras++; diff --git a/Triangle.NET/Triangle/Data/Triangle.cs b/Triangle.NET/Triangle/Data/Triangle.cs index 4143590..ca74619 100644 --- a/Triangle.NET/Triangle/Data/Triangle.cs +++ b/Triangle.NET/Triangle/Data/Triangle.cs @@ -48,7 +48,7 @@ namespace TriangleNet.Data // Three NULL vertices. vertices = new Vertex[3]; - if (Behavior.UseSegments) + // TODO: if (Behavior.UseSegments) { // Initialize the three adjoining subsegments to be the // omnipresent subsegment. @@ -63,10 +63,11 @@ namespace TriangleNet.Data attributes = new double[numAttributes]; } - if (Behavior.VarArea) - { - area = -1.0; - } + // TODO: + //if (Behavior.VarArea) + //{ + // area = -1.0; + //} } diff --git a/Triangle.NET/Triangle/IO/DataReader.cs b/Triangle.NET/Triangle/IO/DataReader.cs index 1ba541f..596cabf 100644 --- a/Triangle.NET/Triangle/IO/DataReader.cs +++ b/Triangle.NET/Triangle/IO/DataReader.cs @@ -85,7 +85,7 @@ namespace TriangleNet.IO //tri.triangle.neighbors[0].triangle = tri.triangle; } - if (Behavior.Poly) + if (mesh.behavior.Poly) { mesh.insegments = numberofsegments; @@ -134,7 +134,7 @@ namespace TriangleNet.IO tri.triangle.attributes = triangles[i].Attributes; // TODO - if (Behavior.VarArea) + if (mesh.behavior.VarArea) { tri.triangle.area = triangles[i].Area; } @@ -197,7 +197,7 @@ namespace TriangleNet.IO // Prepare to count the boundary edges. hullsize = 0; - if (Behavior.Poly) + if (mesh.behavior.Poly) { // Read the segments from the .poly file, and link them // to their neighboring triangles. diff --git a/Triangle.NET/Triangle/IO/FileWriter.cs b/Triangle.NET/Triangle/IO/FileWriter.cs index a0a9cd2..0ded508 100644 --- a/Triangle.NET/Triangle/IO/FileWriter.cs +++ b/Triangle.NET/Triangle/IO/FileWriter.cs @@ -43,7 +43,9 @@ namespace TriangleNet.IO Vertex vertex; long outvertices = mesh.vertices.Count; - if (Behavior.Jettison) + Behavior behavior = mesh.behavior; + + if (behavior.Jettison) { outvertices = mesh.vertices.Count - mesh.undeads; } @@ -55,13 +57,13 @@ namespace TriangleNet.IO // Number of vertices, number of dimensions, number of vertex attributes, // and number of boundary markers (zero or one). writer.WriteLine("{0} {1} {2} {3}", outvertices, mesh.mesh_dim, mesh.nextras, - Behavior.UseBoundaryMarkers ? "1" : "0"); + behavior.UseBoundaryMarkers ? "1" : "0"); foreach (var item in mesh.vertices.Values) { vertex = item; - if (!Behavior.Jettison || vertex.type != VertexType.UndeadVertex) + if (!behavior.Jettison || vertex.type != VertexType.UndeadVertex) { // Vertex number, x and y coordinates. writer.Write("{0} {1} {2}", index, vertex.x.ToString(nfi), vertex.y.ToString(nfi)); @@ -72,7 +74,7 @@ namespace TriangleNet.IO writer.Write(" {0}", vertex.attributes[j].ToString(nfi)); } - if (Behavior.UseBoundaryMarkers) + if (behavior.UseBoundaryMarkers) { // Write the boundary marker. writer.Write(" {0}", vertex.mark); @@ -154,6 +156,8 @@ namespace TriangleNet.IO Osub subseg = default(Osub); Vertex pt1, pt2; + bool useBoundaryMarkers = mesh.behavior.UseBoundaryMarkers; + using (StreamWriter writer = new StreamWriter(filename)) { if (writeNodes) @@ -167,12 +171,12 @@ namespace TriangleNet.IO // Followed by number of dimensions, number of vertex attributes, // and number of boundary markers (zero or one). writer.WriteLine("0 {0} {1} {2}", mesh.mesh_dim, mesh.nextras, - Behavior.UseBoundaryMarkers ? "1" : "0"); + useBoundaryMarkers ? "1" : "0"); } // Number of segments, number of boundary markers (zero or one). writer.WriteLine("{0} {1}", mesh.subsegs.Count, - Behavior.UseBoundaryMarkers ? "1" : "0"); + useBoundaryMarkers ? "1" : "0"); subseg.orient = 0; @@ -185,7 +189,7 @@ namespace TriangleNet.IO pt2 = subseg.Dest(); // Segment number, indices of its two endpoints, and possibly a marker. - if (Behavior.UseBoundaryMarkers) + if (useBoundaryMarkers) { writer.WriteLine("{0} {1} {2} {3}", j, pt1.id, pt2.id, subseg.seg.boundary); } @@ -232,10 +236,12 @@ namespace TriangleNet.IO Osub checkmark = default(Osub); Vertex p1, p2; + Behavior behavior = mesh.behavior; + using (StreamWriter writer = new StreamWriter(filename)) { // Number of edges, number of boundary markers (zero or one). - writer.WriteLine("{0} {1}", mesh.edges, Behavior.UseBoundaryMarkers ? "1" : "0"); + writer.WriteLine("{0} {1}", mesh.edges, behavior.UseBoundaryMarkers ? "1" : "0"); long index = 0; // To loop over the set of edges, loop over all triangles, and look at @@ -256,11 +262,11 @@ namespace TriangleNet.IO p1 = tri.Org(); p2 = tri.Dest(); - if (Behavior.UseBoundaryMarkers) + if (behavior.UseBoundaryMarkers) { // Edge number, indices of two endpoints, and a boundary marker. // If there's no subsegment, the boundary marker is zero. - if (Behavior.UseSegments) + if (behavior.UseSegments) { tri.SegPivot(ref checkmark); @@ -371,7 +377,7 @@ namespace TriangleNet.IO torg = tri.Org(); tdest = tri.Dest(); tapex = tri.Apex(); - circumcenter = Primitives.FindCircumcenter(torg, tdest, tapex, ref xi, ref eta, false); + circumcenter = Primitives.FindCircumcenter(torg, tdest, tapex, ref xi, ref eta); // X and y coordinates. writer.Write("{0} {1} {2}", index, circumcenter.X.ToString(nfi), @@ -456,7 +462,7 @@ namespace TriangleNet.IO long outvertices = mesh.vertices.Count; - if (Behavior.Jettison) + if (mesh.behavior.Jettison) { outvertices = mesh.vertices.Count - mesh.undeads; } @@ -472,7 +478,7 @@ namespace TriangleNet.IO { p1 = item; - if (!Behavior.Jettison || p1.type != VertexType.UndeadVertex) + if (!mesh.behavior.Jettison || p1.type != VertexType.UndeadVertex) { // The "0.0" is here because the OFF format uses 3D coordinates. writer.WriteLine(" {0} {1} 0.0", p1[0].ToString(nfi), p1[1].ToString(nfi)); diff --git a/Triangle.NET/Triangle/Mesh.cs b/Triangle.NET/Triangle/Mesh.cs index 61bfdd2..da9d25f 100644 --- a/Triangle.NET/Triangle/Mesh.cs +++ b/Triangle.NET/Triangle/Mesh.cs @@ -81,7 +81,8 @@ namespace TriangleNet // proximate vertices are inserted sequentially. internal Otri recenttri; - //static FlipStacker dummyflip; + // Controls the behavior of the mesh instance. + internal Behavior behavior; #endregion @@ -166,7 +167,7 @@ namespace TriangleNet { logger = SimpleLog.Instance; - Behavior.Init(); + behavior = new Behavior(); vertices = new Dictionary(); triangles = new Dictionary(); @@ -222,28 +223,28 @@ namespace TriangleNet if (input.HasSegments) { - Behavior.Poly = true; + behavior.Poly = true; } //if (input.EdgeMarkers != null) //{ - // Behavior.UseBoundaryMarkers = true; + // behavior.UseBoundaryMarkers = true; //} //if (input.TriangleAreas != null) //{ - // Behavior.VarArea = true; + // behavior.VarArea = true; //} - if (!Behavior.Poly) + if (!behavior.Poly) { // Be careful not to allocate space for element area constraints that // will never be assigned any value (other than the default -1.0). - Behavior.VarArea = false; + behavior.VarArea = false; // Be careful not to add an extra attribute to each element unless the // input supports it (PSLG in, but not refining a preexisting mesh). - Behavior.RegionAttrib = false; + behavior.RegionAttrib = false; } TransferNodes(input); @@ -277,26 +278,26 @@ namespace TriangleNet if (input.HasSegments) { - Behavior.Poly = true; + behavior.Poly = true; } //if (input.EdgeMarkers != null) //{ - // Behavior.UseBoundaryMarkers = true; + // behavior.UseBoundaryMarkers = true; //} - if (!Behavior.Poly) + if (!behavior.Poly) { // Be careful not to allocate space for element area constraints that // will never be assigned any value (other than the default -1.0). - Behavior.VarArea = false; + behavior.VarArea = false; // Be careful not to add an extra attribute to each element unless the // input supports it (PSLG in, but not refining a preexisting mesh). - Behavior.RegionAttrib = false; + behavior.RegionAttrib = false; } - steinerleft = Behavior.Steiner; + steinerleft = behavior.Steiner; TransferNodes(input); @@ -308,7 +309,7 @@ namespace TriangleNet infvertex2 = null; infvertex3 = null; - if (Behavior.UseSegments) + if (behavior.UseSegments) { // Segments will be introduced next. checksegments = true; @@ -317,7 +318,7 @@ namespace TriangleNet FormSkeleton(input); } - if (Behavior.Poly && (triangles.Count > 0)) + if (behavior.Poly && (triangles.Count > 0)) { // Copy holes foreach (var item in input.holes) @@ -348,7 +349,7 @@ namespace TriangleNet regions.Clear(); } - if (Behavior.Quality && (triangles.Count > 0)) + if (behavior.Quality && (triangles.Count > 0)) { quality.EnforceQuality(); // Enforce angle and area constraints. } @@ -394,14 +395,14 @@ namespace TriangleNet /// Global area constraint. public void Refine(double areaConstraint) { - Behavior.FixedArea = true; - Behavior.MaxArea = areaConstraint; + behavior.FixedArea = true; + behavior.MaxArea = areaConstraint; this.Refine(); // Reset option for sanity - Behavior.FixedArea = false; - Behavior.MaxArea = -1.0; + behavior.FixedArea = false; + behavior.MaxArea = -1.0; } /// @@ -414,9 +415,9 @@ namespace TriangleNet // TODO: Set all vertex types to input (i.e. NOT free)? - if (Behavior.Poly) + if (behavior.Poly) { - if (Behavior.UseSegments) + if (behavior.UseSegments) { insegments = subsegs.Count; } @@ -428,7 +429,7 @@ namespace TriangleNet Reset(); - steinerleft = Behavior.Steiner; + steinerleft = behavior.Steiner; // Ensure that no vertex can be mistaken for a triangular bounding // box vertex in insertvertex(). @@ -436,7 +437,7 @@ namespace TriangleNet infvertex2 = null; infvertex3 = null; - if (Behavior.UseSegments) + if (behavior.UseSegments) { checksegments = true; } @@ -502,23 +503,23 @@ namespace TriangleNet { if (option == Options.ConformingDelaunay) { - Behavior.ConformDel = value; - Behavior.Quality = value; // TODO: ok? + behavior.ConformDel = value; + behavior.Quality = value; // TODO: ok? return; } else if (option == Options.BoundaryMarkers) { - Behavior.UseBoundaryMarkers = value; + behavior.UseBoundaryMarkers = value; return; } else if (option == Options.Quality) { - Behavior.Quality = value; + behavior.Quality = value; if (value) { - Behavior.MinAngle = 20.0; - Behavior.MaxAngle = 140.0; + behavior.MinAngle = 20.0; + behavior.MaxAngle = 140.0; UpdateOptions(); } @@ -526,7 +527,7 @@ namespace TriangleNet } else if (option == Options.Convex) { - Behavior.Convex = value; + behavior.Convex = value; return; } @@ -543,23 +544,23 @@ namespace TriangleNet if (option == Options.MinAngle) { - Behavior.MinAngle = value; - Behavior.Quality = (value >= 0); // TODO: ok? + behavior.MinAngle = value; + behavior.Quality = (value >= 0); // TODO: ok? UpdateOptions(); return; } else if (option == Options.MaxAngle) { - Behavior.MaxAngle = value; - Behavior.Quality = (value >= 0); // TODO: ok? + behavior.MaxAngle = value; + behavior.Quality = (value >= 0); // TODO: ok? UpdateOptions(); return; } else if (option == Options.MaxArea) { - Behavior.MaxArea = value; - Behavior.FixedArea = true; - Behavior.Quality = (value >= 0); // TODO: ok? + behavior.MaxArea = value; + behavior.FixedArea = true; + behavior.Quality = (value >= 0); // TODO: ok? return; } @@ -575,33 +576,33 @@ namespace TriangleNet { if (option == Options.NoBisect) { - Behavior.NoBisect = value; + behavior.NoBisect = value; return; } else if (option == Options.SteinerPoints) { - Behavior.Steiner = value; + behavior.Steiner = value; return; } else if (option == Options.MinAngle) { - Behavior.MinAngle = value; - Behavior.Quality = (value >= 0); // TODO: ok? + behavior.MinAngle = value; + behavior.Quality = (value >= 0); // TODO: ok? UpdateOptions(); return; } else if (option == Options.MaxAngle) { - Behavior.MaxAngle = value; - Behavior.Quality = (value >= 0); // TODO: ok? + behavior.MaxAngle = value; + behavior.Quality = (value >= 0); // TODO: ok? UpdateOptions(); return; } else if (option == Options.MaxArea) { - Behavior.MaxArea = value; - Behavior.Quality = (value >= 0); // TODO: ok? - Behavior.FixedArea = true; + behavior.MaxArea = value; + behavior.Quality = (value >= 0); // TODO: ok? + behavior.FixedArea = true; return; } @@ -613,20 +614,20 @@ namespace TriangleNet /// private void UpdateOptions() { - Behavior.UseSegments = Behavior.Poly || Behavior.Quality || Behavior.Convex; - Behavior.GoodAngle = Math.Cos(Behavior.MinAngle * Math.PI / 180.0); - Behavior.MaxGoodAngle = Math.Cos(Behavior.MaxAngle * Math.PI / 180.0); + behavior.UseSegments = behavior.Poly || behavior.Quality || behavior.Convex; + behavior.GoodAngle = Math.Cos(behavior.MinAngle * Math.PI / 180.0); + behavior.MaxGoodAngle = Math.Cos(behavior.MaxAngle * Math.PI / 180.0); - if (Behavior.GoodAngle == 1.0) + if (behavior.GoodAngle == 1.0) { - Behavior.Offconstant = 0.0; + behavior.Offconstant = 0.0; } else { - Behavior.Offconstant = 0.475 * Math.Sqrt((1.0 + Behavior.GoodAngle) / (1.0 - Behavior.GoodAngle)); + behavior.Offconstant = 0.475 * Math.Sqrt((1.0 + behavior.GoodAngle) / (1.0 - behavior.GoodAngle)); } - Behavior.GoodAngle *= Behavior.GoodAngle; + behavior.GoodAngle *= behavior.GoodAngle; } #endregion @@ -643,12 +644,12 @@ namespace TriangleNet eextras = 0; - if (Behavior.Algorithm == TriangulationAlgorithm.Dwyer) + if (behavior.Algorithm == TriangulationAlgorithm.Dwyer) { Dwyer alg = new Dwyer(); hulledges = alg.Triangulate(this); } - else if (Behavior.Algorithm == TriangulationAlgorithm.SweepLine) + else if (behavior.Algorithm == TriangulationAlgorithm.SweepLine) { SweepLine alg = new SweepLine(); hulledges = alg.Triangulate(this); @@ -750,7 +751,7 @@ namespace TriangleNet dummytri.neighbors[1].triangle = dummytri; dummytri.neighbors[2].triangle = dummytri; - if (Behavior.UseSegments) + if (behavior.UseSegments) { // Set up 'dummysub', the omnipresent subsegment pointed to by any // triangle side or subsegment end that isn't attached to a real @@ -973,8 +974,8 @@ namespace TriangleNet // The vertex falls on a subsegment, and hence will not be inserted. if (segmentflaws) { - enq = Behavior.NoBisect != 2; - if (enq && (Behavior.NoBisect == 1)) + enq = behavior.NoBisect != 2; + if (enq && (behavior.NoBisect == 1)) { // This subsegment may be split only if it is an // internal boundary. @@ -1035,7 +1036,7 @@ namespace TriangleNet newbotright.triangle.attributes[i] = botright.triangle.attributes[i]; } - if (Behavior.VarArea) + if (behavior.VarArea) { // Set the area constraint of a new triangle. newbotright.triangle.area = botright.triangle.area; @@ -1055,7 +1056,7 @@ namespace TriangleNet newtopright.triangle.attributes[i] = topright.triangle.attributes[i]; } - if (Behavior.VarArea) + if (behavior.VarArea) { // Set the area constraint of another new triangle. newtopright.triangle.area = topright.triangle.area; @@ -1166,7 +1167,7 @@ namespace TriangleNet newbotright.triangle.attributes[i] = attrib; } - if (Behavior.VarArea) + if (behavior.VarArea) { // Set the area constraint of the new triangles. area = horiz.triangle.area; @@ -1369,7 +1370,7 @@ namespace TriangleNet horiz.triangle.attributes[i] = attrib; } - if (Behavior.VarArea) + if (behavior.VarArea) { if ((top.triangle.area <= 0.0) || (horiz.triangle.area <= 0.0)) { @@ -1904,7 +1905,7 @@ namespace TriangleNet // the resulting triangles. deltri.Onext(ref firstedge); deltri.Oprev(ref lastedge); - TriangulatePolygon(firstedge, lastedge, edgecount, false, Behavior.NoBisect == 0); + TriangulatePolygon(firstedge, lastedge, edgecount, false, behavior.NoBisect == 0); } // Splice out two triangles. deltri.Lprev(ref deltriright); @@ -1928,7 +1929,7 @@ namespace TriangleNet // Set the new origin of 'deltri' and check its quality. neworg = lefttri.Org(); deltri.SetOrg(neworg); - if (Behavior.NoBisect == 0) + if (behavior.NoBisect == 0) { quality.TestTriangle(ref deltri); } @@ -3028,7 +3029,7 @@ namespace TriangleNet this.insegments = 0; - if (Behavior.Poly) + if (behavior.Poly) { // If the input vertices are collinear, there is no triangulation, // so don't try to insert segments. @@ -3092,7 +3093,7 @@ namespace TriangleNet } } - if (Behavior.Convex || !Behavior.Poly) + if (behavior.Convex || !behavior.Poly) { // Enclose the convex hull with subsegments. MarkHull(); diff --git a/Triangle.NET/Triangle/NewLocation.cs b/Triangle.NET/Triangle/NewLocation.cs index c4cf09f..fa1f7ca 100644 --- a/Triangle.NET/Triangle/NewLocation.cs +++ b/Triangle.NET/Triangle/NewLocation.cs @@ -18,14 +18,23 @@ namespace TriangleNet /// /// http://www.cise.ufl.edu/~ungor/aCute/index.html /// - static class NewLocation + class NewLocation { const double EPS = 0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000001; + Mesh mesh; + Behavior behavior; + + public NewLocation(Mesh mesh) + { + this.mesh = mesh; + this.behavior = mesh.behavior; + } + /// /// Find a new location for a Steiner point. /// - /// + /// /// /// /// @@ -34,17 +43,17 @@ namespace TriangleNet /// /// /// - public static void FindLocation(Mesh m, Vertex torg, Vertex tdest, Vertex tapex, + public void FindLocation(Vertex torg, Vertex tdest, Vertex tapex, Vertex circumcenter, ref double xi, ref double eta, bool offcenter, Otri badotri) // TODO: ref circumcenter??? { // Based on using -U switch, call the corresponding function - if (Behavior.MaxAngle == 0.0) + if (behavior.MaxAngle == 0.0) { - FindNewLocationWithoutMaxAngle(m, torg, tdest, tapex, circumcenter, ref xi, ref eta, true, badotri); + FindNewLocationWithoutMaxAngle(torg, tdest, tapex, circumcenter, ref xi, ref eta, true, badotri); } else { - FindNewLocation(m, torg, tdest, tapex, circumcenter, ref xi, ref eta, true, badotri); + FindNewLocation(torg, tdest, tapex, circumcenter, ref xi, ref eta, true, badotri); } } @@ -52,7 +61,6 @@ namespace TriangleNet /// /// Find a new location for a Steiner point. /// - /// /// /// /// @@ -61,9 +69,10 @@ namespace TriangleNet /// /// /// - static void FindNewLocationWithoutMaxAngle(Mesh m, Vertex torg, Vertex tdest, Vertex tapex, + private void FindNewLocationWithoutMaxAngle(Vertex torg, Vertex tdest, Vertex tapex, Vertex circumcenter, ref double xi, ref double eta, bool offcenter, Otri badotri) { + double offconstant = behavior.Offconstant; // for calculating the distances of the edges double xdo, ydo, xao, yao, xda, yda; @@ -284,14 +293,14 @@ namespace TriangleNet }// end of switch // check for offcenter condition - if (offcenter && (Behavior.Offconstant > 0.0)) + if (offcenter && (offconstant > 0.0)) { // origin has the smallest angle if (orientation == 213 || orientation == 231) { // Find the position of the off-center, as described by Alper Ungor. - dxoff = 0.5 * xShortestEdge - Behavior.Offconstant * yShortestEdge; - dyoff = 0.5 * yShortestEdge + Behavior.Offconstant * xShortestEdge; + dxoff = 0.5 * xShortestEdge - offconstant * yShortestEdge; + dyoff = 0.5 * yShortestEdge + offconstant * xShortestEdge; // If the off-center is closer to destination than the // circumcenter, use the off-center instead. /// doubleLY BAD CASE /// @@ -311,8 +320,8 @@ namespace TriangleNet else if (orientation == 123 || orientation == 132) { // Find the position of the off-center, as described by Alper Ungor. - dxoff = 0.5 * xShortestEdge + Behavior.Offconstant * yShortestEdge; - dyoff = 0.5 * yShortestEdge - Behavior.Offconstant * xShortestEdge; + dxoff = 0.5 * xShortestEdge + offconstant * yShortestEdge; + dyoff = 0.5 * yShortestEdge - offconstant * xShortestEdge; // If the off-center is closer to the origin than the // circumcenter, use the off-center instead. /// doubleLY BAD CASE /// @@ -331,8 +340,8 @@ namespace TriangleNet else {//orientation == 312 || orientation == 321 // Find the position of the off-center, as described by Alper Ungor. - dxoff = 0.5 * xShortestEdge - Behavior.Offconstant * yShortestEdge; - dyoff = 0.5 * yShortestEdge + Behavior.Offconstant * xShortestEdge; + dxoff = 0.5 * xShortestEdge - offconstant * yShortestEdge; + dyoff = 0.5 * yShortestEdge + offconstant * xShortestEdge; // If the off-center is closer to the origin than the // circumcenter, use the off-center instead. /// doubleLY BAD CASE /// @@ -371,7 +380,7 @@ namespace TriangleNet } /// RELOCATION (LOCAL SMOOTHING) /// /// check for possible relocation of one of triangle's points /// - relocated = DoSmoothing(m, delotri, torg, tdest, tapex, ref newloc); + relocated = DoSmoothing(delotri, torg, tdest, tapex, ref newloc); /// if relocation is possible, delete that vertex and insert a vertex at the new location /// if (relocated > 0) { @@ -385,17 +394,17 @@ namespace TriangleNet { case 1: //printf("Relocate: (%f,%f)\n", torg[0],torg[1]); - m.DeleteVertex(ref delotri); + mesh.DeleteVertex(ref delotri); break; case 2: //printf("Relocate: (%f,%f)\n", tdest[0],tdest[1]); delotri.LnextSelf(); - m.DeleteVertex(ref delotri); + mesh.DeleteVertex(ref delotri); break; case 3: //printf("Relocate: (%f,%f)\n", tapex[0],tapex[1]); delotri.LprevSelf(); - m.DeleteVertex(ref delotri); + mesh.DeleteVertex(ref delotri); break; } @@ -405,7 +414,7 @@ namespace TriangleNet // calculate radius of the petal according to angle constraint // first find the visible region, PETAL // find the center of the circle and radius - petalRadius = Math.Sqrt(shortestEdgeDist) / (2 * Math.Sin(Behavior.MinAngle * Math.PI / 180.0)); + petalRadius = Math.Sqrt(shortestEdgeDist) / (2 * Math.Sin(behavior.MinAngle * Math.PI / 180.0)); /// compute two possible centers of the petal /// // finding the center // first find the middle point of smallest edge @@ -439,7 +448,7 @@ namespace TriangleNet } /// find the third point of the neighbor triangle /// - neighborNotFound = GetNeighborsVertex(m, badotri, middleAngleCorner.x, middleAngleCorner.y, + neighborNotFound = GetNeighborsVertex(badotri, middleAngleCorner.x, middleAngleCorner.y, smallestAngleCorner.x, smallestAngleCorner.y, ref thirdPoint, ref neighborotri); /// find the circumcenter of the neighbor triangle /// dxFirstSuggestion = dx; // if we cannot find any appropriate suggestion, we use circumcenter @@ -451,7 +460,9 @@ namespace TriangleNet neighborvertex_2 = neighborotri.Dest(); neighborvertex_3 = neighborotri.Apex(); // now calculate neighbor's circumcenter which is the voronoi site - neighborCircumcenter = Primitives.FindCircumcenter(neighborvertex_1, neighborvertex_2, neighborvertex_3, ref xi_tmp, ref eta_tmp, false); + neighborCircumcenter = Primitives.FindCircumcenter(neighborvertex_1, neighborvertex_2, neighborvertex_3, + ref xi_tmp, ref eta_tmp); + /// compute petal and Voronoi edge intersection /// // in order to avoid degenerate cases, we need to do a vector based calculation for line vector_x = (middleAngleCorner.y - smallestAngleCorner.y);//(-y, x) @@ -568,7 +579,7 @@ namespace TriangleNet /// DO THE SAME THING FOR THE OTHER DIRECTION /// /// find the third point of the neighbor triangle /// - neighborNotFound = GetNeighborsVertex(m, badotri, largestAngleCorner.x, largestAngleCorner.y, + neighborNotFound = GetNeighborsVertex(badotri, largestAngleCorner.x, largestAngleCorner.y, smallestAngleCorner.x, smallestAngleCorner.y, ref thirdPoint, ref neighborotri); /// find the circumcenter of the neighbor triangle /// dxSecondSuggestion = dx; // if we cannot find any appropriate suggestion, we use circumcenter @@ -580,7 +591,8 @@ namespace TriangleNet neighborvertex_2 = neighborotri.Dest(); neighborvertex_3 = neighborotri.Apex(); // now calculate neighbor's circumcenter which is the voronoi site - neighborCircumcenter = Primitives.FindCircumcenter(neighborvertex_1, neighborvertex_2, neighborvertex_3, ref xi_tmp, ref eta_tmp, false); + neighborCircumcenter = Primitives.FindCircumcenter(neighborvertex_1, neighborvertex_2, neighborvertex_3, + ref xi_tmp, ref eta_tmp); /// compute petal and Voronoi edge intersection /// // in order to avoid degenerate cases, we need to do a vector based calculation for line @@ -747,7 +759,6 @@ namespace TriangleNet /// /// Find a new location for a Steiner point. /// - /// /// /// /// @@ -756,10 +767,10 @@ namespace TriangleNet /// /// /// - static void FindNewLocation(Mesh m, - Vertex torg, Vertex tdest, Vertex tapex, + private void FindNewLocation(Vertex torg, Vertex tdest, Vertex tapex, Vertex circumcenter, ref double xi, ref double eta, bool offcenter, Otri badotri) { + double offconstant = behavior.Offconstant; // for calculating the distances of the edges double xdo, ydo, xao, yao, xda, yda; @@ -992,14 +1003,14 @@ namespace TriangleNet }// end of switch // check for offcenter condition - if (offcenter && (Behavior.Offconstant > 0.0)) + if (offcenter && (offconstant > 0.0)) { // origin has the smallest angle if (orientation == 213 || orientation == 231) { // Find the position of the off-center, as described by Alper Ungor. - dxoff = 0.5 * xShortestEdge - Behavior.Offconstant * yShortestEdge; - dyoff = 0.5 * yShortestEdge + Behavior.Offconstant * xShortestEdge; + dxoff = 0.5 * xShortestEdge - offconstant * yShortestEdge; + dyoff = 0.5 * yShortestEdge + offconstant * xShortestEdge; // If the off-center is closer to destination than the // circumcenter, use the off-center instead. /// doubleLY BAD CASE /// @@ -1019,8 +1030,8 @@ namespace TriangleNet else if (orientation == 123 || orientation == 132) { // Find the position of the off-center, as described by Alper Ungor. - dxoff = 0.5 * xShortestEdge + Behavior.Offconstant * yShortestEdge; - dyoff = 0.5 * yShortestEdge - Behavior.Offconstant * xShortestEdge; + dxoff = 0.5 * xShortestEdge + offconstant * yShortestEdge; + dyoff = 0.5 * yShortestEdge - offconstant * xShortestEdge; // If the off-center is closer to the origin than the // circumcenter, use the off-center instead. /// doubleLY BAD CASE /// @@ -1039,8 +1050,8 @@ namespace TriangleNet else {//orientation == 312 || orientation == 321 // Find the position of the off-center, as described by Alper Ungor. - dxoff = 0.5 * xShortestEdge - Behavior.Offconstant * yShortestEdge; - dyoff = 0.5 * yShortestEdge + Behavior.Offconstant * xShortestEdge; + dxoff = 0.5 * xShortestEdge - offconstant * yShortestEdge; + dyoff = 0.5 * yShortestEdge + offconstant * xShortestEdge; // If the off-center is closer to the origin than the // circumcenter, use the off-center instead. /// doubleLY BAD CASE /// @@ -1079,7 +1090,7 @@ namespace TriangleNet } /// RELOCATION (LOCAL SMOOTHING) /// /// check for possible relocation of one of triangle's points /// - relocated = DoSmoothing(m, delotri, torg, tdest, tapex, ref newloc); + relocated = DoSmoothing(delotri, torg, tdest, tapex, ref newloc); /// if relocation is possible, delete that vertex and insert a vertex at the new location /// if (relocated > 0) { @@ -1093,17 +1104,17 @@ namespace TriangleNet { case 1: //printf("Relocate: (%f,%f)\n", torg[0],torg[1]); - m.DeleteVertex(ref delotri); + mesh.DeleteVertex(ref delotri); break; case 2: //printf("Relocate: (%f,%f)\n", tdest[0],tdest[1]); delotri.LnextSelf(); - m.DeleteVertex(ref delotri); + mesh.DeleteVertex(ref delotri); break; case 3: //printf("Relocate: (%f,%f)\n", tapex[0],tapex[1]); delotri.LprevSelf(); - m.DeleteVertex(ref delotri); + mesh.DeleteVertex(ref delotri); break; } } @@ -1114,9 +1125,9 @@ namespace TriangleNet // find the center of the circle and radius // choose minimum angle as the maximum of quality angle and the minimum angle of the bad triangle minangle = Math.Acos((middleEdgeDist + longestEdgeDist - shortestEdgeDist) / (2 * Math.Sqrt(middleEdgeDist) * Math.Sqrt(longestEdgeDist))) * 180.0 / Math.PI; - if (Behavior.MinAngle > minangle) + if (behavior.MinAngle > minangle) { - minangle = Behavior.MinAngle; + minangle = behavior.MinAngle; } else { @@ -1155,7 +1166,7 @@ namespace TriangleNet xPetalCtr = xPetalCtr_2; yPetalCtr = yPetalCtr_2; } /// find the third point of the neighbor triangle /// - neighborNotFound_first = GetNeighborsVertex(m, badotri, middleAngleCorner.x, middleAngleCorner.y, + neighborNotFound_first = GetNeighborsVertex(badotri, middleAngleCorner.x, middleAngleCorner.y, smallestAngleCorner.x, smallestAngleCorner.y, ref thirdPoint, ref neighborotri); /// find the circumcenter of the neighbor triangle /// dxFirstSuggestion = dx; // if we cannot find any appropriate suggestion, we use circumcenter @@ -1171,7 +1182,7 @@ namespace TriangleNet // find the third point other than p and q petal_bisector_x = xPetalCtr + line_vector_x * petalRadius; petal_bisector_y = yPetalCtr + line_vector_y * petalRadius; - alpha = (2.0 * Behavior.MaxAngle + minangle - 180.0) * Math.PI / 180.0; + alpha = (2.0 * behavior.MaxAngle + minangle - 180.0) * Math.PI / 180.0; // rotate the vector cw around the petal center x_1 = petal_bisector_x * Math.Cos(alpha) + petal_bisector_y * Math.Sin(alpha) + xPetalCtr - xPetalCtr * Math.Cos(alpha) - yPetalCtr * Math.Sin(alpha); y_1 = -petal_bisector_x * Math.Sin(alpha) + petal_bisector_y * Math.Cos(alpha) + yPetalCtr + xPetalCtr * Math.Sin(alpha) - yPetalCtr * Math.Cos(alpha); @@ -1207,7 +1218,9 @@ namespace TriangleNet neighborvertex_2 = neighborotri.Dest(); neighborvertex_3 = neighborotri.Apex(); // now calculate neighbor's circumcenter which is the voronoi site - neighborCircumcenter = Primitives.FindCircumcenter(neighborvertex_1, neighborvertex_2, neighborvertex_3, ref xi_tmp, ref eta_tmp, false); + neighborCircumcenter = Primitives.FindCircumcenter(neighborvertex_1, neighborvertex_2, neighborvertex_3, + ref xi_tmp, ref eta_tmp); + /// compute petal and Voronoi edge intersection /// // in order to avoid degenerate cases, we need to do a vector based calculation for line vector_x = (middleAngleCorner.y - smallestAngleCorner.y);//(-y, x) @@ -1285,7 +1298,7 @@ namespace TriangleNet (smallestAngleCorner.y - line_inter_y) * (smallestAngleCorner.y - line_inter_y))) && (IsBadTriangleAngle(middleAngleCorner.x, middleAngleCorner.y, largestAngleCorner.x, largestAngleCorner.y, petal_slab_inter_x_first, petal_slab_inter_y_first)) - && MinDistanceToNeighbor(m, petal_slab_inter_x_first, petal_slab_inter_y_first, ref neighborotri) > MinDistanceToNeighbor(m, line_inter_x, line_inter_y, ref neighborotri)) + && MinDistanceToNeighbor(petal_slab_inter_x_first, petal_slab_inter_y_first, ref neighborotri) > MinDistanceToNeighbor(line_inter_x, line_inter_y, ref neighborotri)) { // /// check the neighbor's vertices also, which one if better @@ -1379,7 +1392,7 @@ namespace TriangleNet (smallestAngleCorner.y - line_inter_y) * (smallestAngleCorner.y - line_inter_y))) && (IsBadTriangleAngle(middleAngleCorner.x, middleAngleCorner.y, largestAngleCorner.x, largestAngleCorner.y, petal_slab_inter_x_first, petal_slab_inter_y_first)) - && MinDistanceToNeighbor(m, petal_slab_inter_x_first, petal_slab_inter_y_first, ref neighborotri) > MinDistanceToNeighbor(m, line_inter_x, line_inter_y, ref neighborotri)) + && MinDistanceToNeighbor(petal_slab_inter_x_first, petal_slab_inter_y_first, ref neighborotri) > MinDistanceToNeighbor(line_inter_x, line_inter_y, ref neighborotri)) { //slab and petal intersection is advised dxFirstSuggestion = petal_slab_inter_x_first - torg.x; @@ -1497,7 +1510,7 @@ namespace TriangleNet /// DO THE SAME THING FOR THE OTHER DIRECTION /// /// find the third point of the neighbor triangle /// - neighborNotFound_second = GetNeighborsVertex(m, badotri, largestAngleCorner.x, largestAngleCorner.y, + neighborNotFound_second = GetNeighborsVertex(badotri, largestAngleCorner.x, largestAngleCorner.y, smallestAngleCorner.x, smallestAngleCorner.y, ref thirdPoint, ref neighborotri); /// find the circumcenter of the neighbor triangle /// dxSecondSuggestion = dx; // if we cannot find any appropriate suggestion, we use circumcenter @@ -1514,7 +1527,8 @@ namespace TriangleNet neighborvertex_2 = neighborotri.Dest(); neighborvertex_3 = neighborotri.Apex(); // now calculate neighbor's circumcenter which is the voronoi site - neighborCircumcenter = Primitives.FindCircumcenter(neighborvertex_1, neighborvertex_2, neighborvertex_3, ref xi_tmp, ref eta_tmp, false); + neighborCircumcenter = Primitives.FindCircumcenter(neighborvertex_1, neighborvertex_2, neighborvertex_3, + ref xi_tmp, ref eta_tmp); /// compute petal and Voronoi edge intersection /// // in order to avoid degenerate cases, we need to do a vector based calculation for line @@ -1595,7 +1609,7 @@ namespace TriangleNet (smallestAngleCorner.y - line_inter_y) * (smallestAngleCorner.y - line_inter_y))) && (IsBadTriangleAngle(middleAngleCorner.x, middleAngleCorner.y, largestAngleCorner.x, largestAngleCorner.y, petal_slab_inter_x_second, petal_slab_inter_y_second)) - && MinDistanceToNeighbor(m, petal_slab_inter_x_second, petal_slab_inter_y_second, ref neighborotri) > MinDistanceToNeighbor(m, line_inter_x, line_inter_y, ref neighborotri)) + && MinDistanceToNeighbor(petal_slab_inter_x_second, petal_slab_inter_y_second, ref neighborotri) > MinDistanceToNeighbor(line_inter_x, line_inter_y, ref neighborotri)) { // slab and petal intersection is advised dxSecondSuggestion = petal_slab_inter_x_second - torg.x; @@ -1684,7 +1698,7 @@ namespace TriangleNet (smallestAngleCorner.y - line_inter_y) * (smallestAngleCorner.y - line_inter_y))) && (IsBadTriangleAngle(middleAngleCorner.x, middleAngleCorner.y, largestAngleCorner.x, largestAngleCorner.y, petal_slab_inter_x_second, petal_slab_inter_y_second)) - && MinDistanceToNeighbor(m, petal_slab_inter_x_second, petal_slab_inter_y_second, ref neighborotri) > MinDistanceToNeighbor(m, line_inter_x, line_inter_y, ref neighborotri)) + && MinDistanceToNeighbor(petal_slab_inter_x_second, petal_slab_inter_y_second, ref neighborotri) > MinDistanceToNeighbor(line_inter_x, line_inter_y, ref neighborotri)) { // slab and petal intersection is advised dxSecondSuggestion = petal_slab_inter_x_second - torg.x; @@ -2015,7 +2029,7 @@ namespace TriangleNet /// /// /// Returns a number indicating an orientation. - static int LongestShortestEdge(double aodist, double dadist, double dodist) + private int LongestShortestEdge(double aodist, double dadist, double dodist) { // 123: shortest: aodist // 213: shortest: dadist // 312: shortest: dodist // middle: dadist // middle: aodist // middle: aodist @@ -2075,15 +2089,14 @@ namespace TriangleNet /// /// Checks if smothing is possible for a given bad triangle. /// - /// /// /// /// /// /// The new location for the point, if somothing is possible. /// Returns 1, 2 or 3 if smoothing will work, 0 otherwise. - static int DoSmoothing(Mesh m, Otri badotri, - Vertex torg, Vertex tdest, Vertex tapex, ref double[] newloc) + private int DoSmoothing(Otri badotri, Vertex torg, Vertex tdest, Vertex tapex, + ref double[] newloc) { int numpoints_p = 0;// keeps the number of points in a star of point p, q, r @@ -2106,7 +2119,7 @@ namespace TriangleNet //********************* TRY TO RELOCATE POINT "p" *************** // get the surrounding points of p, so this gives us the triangles - numpoints_p = GetStarPoints(m, badotri, torg, tdest, tapex, 1, ref points_p); + numpoints_p = GetStarPoints(badotri, torg, tdest, tapex, 1, ref points_p); // check if the points in counterclockwise order // p1[0] = points_p[0]; p1[1] = points_p[1]; // p2[0] = points_p[2]; p2[1] = points_p[3]; @@ -2130,13 +2143,13 @@ namespace TriangleNet { //newLocFound = getPetalIntersection(m, b, numpoints_p, points_p, newloc); //newLocFound = getPetalIntersectionBruteForce(m, b,numpoints_p, points_p, newloc,torg[0],torg[1]); - if (Behavior.MaxAngle == 0.0) + if (behavior.MaxAngle == 0.0) { - newLocFound = GetWedgeIntersectionWithoutMaxAngle(m, numpoints_p, points_p, ref newloc); + newLocFound = GetWedgeIntersectionWithoutMaxAngle(numpoints_p, points_p, ref newloc); } else { - newLocFound = GetWedgeIntersection(m, numpoints_p, points_p, ref newloc); + newLocFound = GetWedgeIntersection(numpoints_p, points_p, ref newloc); } //printf("call petal intersection for p\n"); // make sure the relocated point is a free vertex @@ -2152,7 +2165,7 @@ namespace TriangleNet //********************* TRY TO RELOCATE POINT "q" *************** // get the surrounding points of q, so this gives us the triangles - numpoints_q = GetStarPoints(m, badotri, torg, tdest, tapex, 2, ref points_q); + numpoints_q = GetStarPoints(badotri, torg, tdest, tapex, 2, ref points_q); // // check if the points in counterclockwise order // v1[0] = points_q[0]; v1[1] = points_q[1]; // v2[0] = points_q[2]; v2[1] = points_q[3]; @@ -2175,13 +2188,13 @@ namespace TriangleNet { //newLocFound = getPetalIntersection(m, b,numpoints_q, points_q, newloc); //newLocFound = getPetalIntersectionBruteForce(m, b,numpoints_q, points_q, newloc,tapex[0],tapex[1]); - if (Behavior.MaxAngle == 0.0) + if (behavior.MaxAngle == 0.0) { - newLocFound = GetWedgeIntersectionWithoutMaxAngle(m, numpoints_q, points_q, ref newloc); + newLocFound = GetWedgeIntersectionWithoutMaxAngle(numpoints_q, points_q, ref newloc); } else { - newLocFound = GetWedgeIntersection(m, numpoints_q, points_q, ref newloc); + newLocFound = GetWedgeIntersection(numpoints_q, points_q, ref newloc); } //printf("call petal intersection for q\n"); @@ -2198,7 +2211,7 @@ namespace TriangleNet //********************* TRY TO RELOCATE POINT "q" *************** // get the surrounding points of r, so this gives us the triangles - numpoints_r = GetStarPoints(m, badotri, torg, tdest, tapex, 3, ref points_r); + numpoints_r = GetStarPoints(badotri, torg, tdest, tapex, 3, ref points_r); // check if the points in counterclockwise order // v1[0] = points_r[0]; v1[1] = points_r[1]; // v2[0] = points_r[2]; v2[1] = points_r[3]; @@ -2221,13 +2234,13 @@ namespace TriangleNet { //newLocFound = getPetalIntersection(m, b,numpoints_r, points_r, newloc); //newLocFound = getPetalIntersectionBruteForce(m, b,numpoints_r, points_r, newloc,tdest[0],tdest[1]); - if (Behavior.MaxAngle == 0.0) + if (behavior.MaxAngle == 0.0) { - newLocFound = GetWedgeIntersectionWithoutMaxAngle(m, numpoints_r, points_r, ref newloc); + newLocFound = GetWedgeIntersectionWithoutMaxAngle(numpoints_r, points_r, ref newloc); } else { - newLocFound = GetWedgeIntersection(m, numpoints_r, points_r, ref newloc); + newLocFound = GetWedgeIntersection(numpoints_r, points_r, ref newloc); } //printf("call petal intersection for r\n"); @@ -2281,7 +2294,6 @@ namespace TriangleNet /// /// Finds the star of a given point. /// - /// /// /// /// @@ -2289,12 +2301,8 @@ namespace TriangleNet /// /// List of points on the star of the given point. /// Number of points on the star of the given point. - static int GetStarPoints(Mesh m, Otri badotri, - Vertex p, - Vertex q, - Vertex r, - int whichPoint, - ref double[] points) + private int GetStarPoints(Otri badotri, Vertex p, Vertex q, Vertex r, + int whichPoint, ref double[] points) { Otri neighotri = default(Otri); // for return value of the function @@ -2345,7 +2353,7 @@ namespace TriangleNet do { // find the neighbor's third point where it is incident to given edge - if (!GetNeighborsVertex(m, tempotri, first_x, first_y, second_x, second_y, ref returnPoint, ref neighotri)) + if (!GetNeighborsVertex(tempotri, first_x, first_y, second_x, second_y, ref returnPoint, ref neighotri)) { // go to next triangle tempotri = neighotri; @@ -2373,7 +2381,6 @@ namespace TriangleNet /// /// Gets a neighbours vertex. /// - /// /// /// /// @@ -2382,7 +2389,7 @@ namespace TriangleNet /// Neighbor's third vertex incident to given edge. /// Pointer for the neighbor triangle. /// Returns true if vertex was found. - static bool GetNeighborsVertex(Mesh m, Otri badotri, + private bool GetNeighborsVertex(Otri badotri, double first_x, double first_y, double second_x, double second_y, ref double[] thirdpoint, ref Otri neighotri) @@ -2534,13 +2541,12 @@ namespace TriangleNet /// /// Find a new point location by wedge intersection. /// - /// /// /// /// A new location for the point according to surrounding points. /// Returns true if new location found - static bool GetWedgeIntersectionWithoutMaxAngle(Mesh m, - int numpoints, double[] points, ref double[] newloc) + private bool GetWedgeIntersectionWithoutMaxAngle(int numpoints, + double[] points, ref double[] newloc) { //double total_x = 0; //double total_y = 0; @@ -2586,9 +2592,9 @@ namespace TriangleNet y1 = points[2 * numpoints - 1]; // minimum angle - alpha = Behavior.MinAngle * Math.PI / 180.0; + alpha = behavior.MinAngle * Math.PI / 180.0; // initialize the constants - if (Behavior.GoodAngle == 1.0) + if (behavior.GoodAngle == 1.0) { petalcenterconstant = 0; petalradiusconstant = 0; @@ -2759,7 +2765,7 @@ namespace TriangleNet /// IF THERE IS A FEASIBLE INTERSECTION POLYGON, FIND ITS CENTROID AS THE NEW LOCATION FindPolyCentroid(numpolypoints, initialConvexPoly, ref newloc); - if (Behavior.FixedArea) + if (behavior.FixedArea) { // numBadTriangle = 0; // for(j= 0; j < numpoints *2-2; j = j+2){ @@ -2794,12 +2800,11 @@ namespace TriangleNet /// /// Find a new point location by wedge intersection. /// - /// /// /// /// A new location for the point according to surrounding points. /// Returns true if new location found - static bool GetWedgeIntersection(Mesh m, int numpoints, double[] points, ref double[] newloc) + private bool GetWedgeIntersection(int numpoints, double[] points, ref double[] newloc) { //double total_x = 0; //double total_y = 0; @@ -2851,15 +2856,15 @@ namespace TriangleNet // minimum / maximum angle double alpha, sinAlpha, cosAlpha, beta, sinBeta, cosBeta; - alpha = Behavior.MinAngle * Math.PI / 180.0; + alpha = behavior.MinAngle * Math.PI / 180.0; sinAlpha = Math.Sin(alpha); cosAlpha = Math.Cos(alpha); - beta = Behavior.MaxAngle * Math.PI / 180.0; + beta = behavior.MaxAngle * Math.PI / 180.0; sinBeta = Math.Sin(beta); cosBeta = Math.Cos(beta); // initialize the constants - if (Behavior.GoodAngle == 1.0) + if (behavior.GoodAngle == 1.0) { petalcenterconstant = 0; petalradiusconstant = 0; @@ -2941,7 +2946,7 @@ namespace TriangleNet /// DETERMINE HOW MANY POINTS TO USE ACCORDING TO THE MINANGLE-MAXANGLE COMBINATION // petal center angle - alpha = (2.0 * Behavior.MaxAngle + Behavior.MinAngle - 180.0); + alpha = (2.0 * behavior.MaxAngle + behavior.MinAngle - 180.0); if (alpha <= 0.0) {// when only angle lines needed // 4 point case @@ -3199,7 +3204,7 @@ namespace TriangleNet /// IF THERE IS A FEASIBLE INTERSECTION POLYGON, FIND ITS CENTROID AS THE NEW LOCATION FindPolyCentroid(numpolypoints, initialConvexPoly, ref newloc); - if (Behavior.MaxAngle != 0.0) + if (behavior.MaxAngle != 0.0) { numBadTriangle = 0; for (j = 0; j < numpoints * 2 - 2; j = j + 2) @@ -3284,7 +3289,7 @@ namespace TriangleNet /// /// /// Returns true if the polygon has angles greater than 2*minangle. - static bool ValidPolygonAngles(int numpoints, double[] points) + private bool ValidPolygonAngles(int numpoints, double[] points) { int i;//,j for (i = 0; i < numpoints; i++) @@ -3326,7 +3331,7 @@ namespace TriangleNet /// /// Returns true, if it is a BAD polygon corner, returns false if it is a GOOD /// polygon corner - static bool IsBadPolygonAngle(double x1, double y1, + private bool IsBadPolygonAngle(double x1, double y1, double x2, double y2, double x3, double y3) { // variables keeping the distance values for the edges @@ -3352,7 +3357,7 @@ namespace TriangleNet cosAngle = (dist12 + dist23 - dist31) / (2 * Math.Sqrt(dist12) * Math.Sqrt(dist23)); // Check whether the angle is smaller than permitted which is 2*minangle!!! //printf("angle: %f 2*minangle = %f\n",acos(cosAngle)*180/PI, 2*acos(Math.Sqrt(b.goodangle))*180/PI); - if (Math.Acos(cosAngle) < 2 * Math.Acos(Math.Sqrt(Behavior.GoodAngle))) + if (Math.Acos(cosAngle) < 2 * Math.Acos(Math.Sqrt(behavior.GoodAngle))) { return true;// it is a BAD triangle } @@ -3375,7 +3380,7 @@ namespace TriangleNet /// // referenced to: http://local.wasp.uwa.edu.au/~pbourke/geometry/ /// - static void LineLineIntersection( + private void LineLineIntersection( double x1, double y1, double x2, double y2, double x3, double y3, @@ -3431,7 +3436,7 @@ namespace TriangleNet /// /// http://www.mathematik.uni-ulm.de/stochastik/lehre/ws03_04/rt/Geometry2D.ps /// - static int HalfPlaneIntersection(int numvertices, ref double[] convexPoly, double x1, double y1, double x2, double y2) + private int HalfPlaneIntersection(int numvertices, ref double[] convexPoly, double x1, double y1, double x2, double y2) { double dx, dy; // direction of the line double z, min, max; @@ -3509,7 +3514,7 @@ namespace TriangleNet /// /// http://www.mathematik.uni-ulm.de/stochastik/lehre/ws03_04/rt/Geometry2D.ps /// - static int SplitConvexPolygon(int numvertices, double[] convexPoly, double x1, double y1, double x2, double y2, ref double[][] polys) + private int SplitConvexPolygon(int numvertices, double[] convexPoly, double x1, double y1, double x2, double y2, ref double[][] polys) { // state = 0: before the first intersection (with the line) // state = 1: after the first intersection (with the line) @@ -3699,7 +3704,7 @@ namespace TriangleNet /// /// http://www.mathematik.uni-ulm.de/stochastik/lehre/ws03_04/rt/Geometry2D.ps /// - static int LinePointLocation(double x1, double y1, double x2, double y2, double x, double y) + private int LinePointLocation(double x1, double y1, double x2, double y2, double x, double y) { double z; if (Math.Atan((y2 - y1) / (x2 - x1)) * 180.0 / Math.PI == 90.0) @@ -3743,7 +3748,7 @@ namespace TriangleNet /// /// referenced to: http://local.wasp.uwa.edu.au/~pbourke/geometry/ /// - static void LineLineSegmentIntersection( + private void LineLineSegmentIntersection( double x1, double y1, double x2, double y2, double x3, double y3, @@ -3806,7 +3811,7 @@ namespace TriangleNet /// /// /// Centroid of a given polygon - static void FindPolyCentroid(int numpoints, double[] points, ref double[] centroid) + private void FindPolyCentroid(int numpoints, double[] points, ref double[] centroid) { int i; //double area = 0.0;//, temp @@ -3838,7 +3843,7 @@ namespace TriangleNet /// /// referenced to: http://local.wasp.uwa.edu.au/~pbourke/geometry/sphereline/ /// - static void CircleLineIntersection( + private void CircleLineIntersection( double x1, double y1, double x2, double y2, double x3, double y3, double r, ref double[] p) @@ -3906,7 +3911,7 @@ namespace TriangleNet /// P3 coordinates (circumcenter point) /// /// Returns true, if given point is the correct one otherwise return false. - static bool ChooseCorrectPoint( + private bool ChooseCorrectPoint( double x1, double y1, double x2, double y2, double x3, double y3, bool isObtuse) @@ -3959,7 +3964,7 @@ namespace TriangleNet /// P3 coordinates [point to be compared] (neighbor's circumcenter) /// P3 coordinates [point to be compared] (neighbor's circumcenter) /// - static void PointBetweenPoints(double x1, double y1, double x2, double y2, double x, double y, ref double[] p) + private void PointBetweenPoints(double x1, double y1, double x2, double y2, double x, double y, ref double[] p) { // now check whether the point is close to circumcenter than intersection point // BETWEEN THE POINTS @@ -3991,7 +3996,7 @@ namespace TriangleNet /// /// /// Returns true, if it is a BAD triangle, returns false if it is a GOOD triangle. - static bool IsBadTriangleAngle(double x1, double y1, double x2, double y2, double x3, double y3) + private bool IsBadTriangleAngle(double x1, double y1, double x2, double y2, double x3, double y3) { // variables keeping the distance values for the edges double dxod, dyod, dxda, dyda, dxao, dyao; @@ -4077,7 +4082,7 @@ namespace TriangleNet // Check whether the angle is smaller than permitted. - if ((angle > Behavior.GoodAngle) || (Behavior.MaxAngle != 0.00 && maxangle < Behavior.MaxGoodAngle)) + if ((angle > behavior.GoodAngle) || (behavior.MaxAngle != 0.00 && maxangle < behavior.MaxGoodAngle)) { return true;// it is a bad triangle } @@ -4089,12 +4094,11 @@ namespace TriangleNet /// Given the triangulation, and a vertex returns the minimum distance to the /// vertices of the triangle where the given vertex located. /// - /// /// /// /// /// - static double MinDistanceToNeighbor(Mesh m, double newlocX, double newlocY, ref Otri searchtri) + private double MinDistanceToNeighbor(double newlocX, double newlocY, ref Otri searchtri) { Otri horiz = default(Otri); // for search operation LocateResult intersect = LocateResult.Outside; @@ -4139,7 +4143,7 @@ namespace TriangleNet // edge specified by 'searchtri'. searchtri.SymSelf(); searchtri.Copy(ref horiz); - intersect = m.PreciseLocate(newvertex, ref horiz, false); + intersect = mesh.PreciseLocate(newvertex, ref horiz, false); } else if (ahead == 0.0) { @@ -4155,7 +4159,7 @@ namespace TriangleNet else { searchtri.Copy(ref horiz); - intersect = m.PreciseLocate(newvertex, ref horiz, false); + intersect = mesh.PreciseLocate(newvertex, ref horiz, false); } } if (intersect == LocateResult.OnVertex || intersect == LocateResult.Outside) diff --git a/Triangle.NET/Triangle/Primitives.cs b/Triangle.NET/Triangle/Primitives.cs index 3f22673..709fe97 100644 --- a/Triangle.NET/Triangle/Primitives.cs +++ b/Triangle.NET/Triangle/Primitives.cs @@ -261,19 +261,10 @@ namespace TriangleNet /// Triangle point. /// Relative coordinate of new location. /// Relative coordinate of new location. - /// Use off-center for new location. + /// Off-center constant. /// Coordinates of the circumcenter (or off-center) - /// - /// The result is returned both in terms of x-y coordinates and xi-eta - /// (barycentric) coordinates. The xi-eta coordinate system is defined in - /// terms of the triangle: the origin of the triangle is the origin of the - /// coordinate system; the destination of the triangle is one unit along the - /// xi axis; and the apex of the triangle is one unit along the eta axis. - /// This procedure also returns the square of the length of the triangle's - /// shortest edge. - /// public static Point FindCircumcenter(Point torg, Point tdest, Point tapex, - ref double xi, ref double eta, bool offcenter) + ref double xi, ref double eta, double offconstant) { double xdo, ydo, xao, yao; double dodist, aodist, dadist; @@ -291,6 +282,7 @@ namespace TriangleNet aodist = xao * xao + yao * yao; dadist = (tdest.x - tapex.x) * (tdest.x - tapex.x) + (tdest.y - tapex.y) * (tdest.y - tapex.y); + if (Behavior.NoExact) { denominator = 0.5 / (xdo * yao - xao * ydo); @@ -304,6 +296,7 @@ namespace TriangleNet // Don't count the above as an orientation test. Statistic.CounterClockwiseCount--; } + dx = (yao * dodist - ydo * aodist) * denominator; dy = (xdo * aodist - xao * dodist) * denominator; @@ -314,11 +307,11 @@ namespace TriangleNet // the input PSLG. if ((dodist < aodist) && (dodist < dadist)) { - if (offcenter && (Behavior.Offconstant > 0.0)) + if (offconstant > 0.0) { // Find the position of the off-center, as described by Alper Ungor. - dxoff = 0.5 * xdo - Behavior.Offconstant * ydo; - dyoff = 0.5 * ydo + Behavior.Offconstant * xdo; + dxoff = 0.5 * xdo - offconstant * ydo; + dyoff = 0.5 * ydo + offconstant * xdo; // If the off-center is closer to the origin than the // circumcenter, use the off-center instead. if (dxoff * dxoff + dyoff * dyoff < dx * dx + dy * dy) @@ -330,10 +323,10 @@ namespace TriangleNet } else if (aodist < dadist) { - if (offcenter && (Behavior.Offconstant > 0.0)) + if (offconstant > 0.0) { - dxoff = 0.5 * xao + Behavior.Offconstant * yao; - dyoff = 0.5 * yao - Behavior.Offconstant * xao; + dxoff = 0.5 * xao + offconstant * yao; + dyoff = 0.5 * yao - offconstant * xao; // If the off-center is closer to the origin than the // circumcenter, use the off-center instead. if (dxoff * dxoff + dyoff * dyoff < dx * dx + dy * dy) @@ -345,12 +338,10 @@ namespace TriangleNet } else { - if (offcenter && (Behavior.Offconstant > 0.0)) + if (offconstant > 0.0) { - dxoff = 0.5 * (tapex.x - tdest.x) - - Behavior.Offconstant * (tapex.y - tdest.y); - dyoff = 0.5 * (tapex.y - tdest.y) + - Behavior.Offconstant * (tapex.x - tdest.x); + dxoff = 0.5 * (tapex.x - tdest.x) - offconstant * (tapex.y - tdest.y); + dyoff = 0.5 * (tapex.y - tdest.y) + offconstant * (tapex.x - tdest.x); // If the off-center is closer to the destination than the // circumcenter, use the off-center instead. if (dxoff * dxoff + dyoff * dyoff < @@ -362,7 +353,71 @@ namespace TriangleNet } } - Point circumcenter = new Point(torg.x + dx, torg.y + dy); + // To interpolate vertex attributes for the new vertex inserted at + // the circumcenter, define a coordinate system with a xi-axis, + // directed from the triangle's origin to its destination, and + // an eta-axis, directed from its origin to its apex. + // Calculate the xi and eta coordinates of the circumcenter. + xi = (yao * dx - xao * dy) * (2.0 * denominator); + eta = (xdo * dy - ydo * dx) * (2.0 * denominator); + + return new Point(torg.x + dx, torg.y + dy); + } + + /// + /// Find the circumcenter of a triangle. + /// + /// Triangle point. + /// Triangle point. + /// Triangle point. + /// Relative coordinate of new location. + /// Relative coordinate of new location. + /// Coordinates of the circumcenter + /// + /// The result is returned both in terms of x-y coordinates and xi-eta + /// (barycentric) coordinates. The xi-eta coordinate system is defined in + /// terms of the triangle: the origin of the triangle is the origin of the + /// coordinate system; the destination of the triangle is one unit along the + /// xi axis; and the apex of the triangle is one unit along the eta axis. + /// This procedure also returns the square of the length of the triangle's + /// shortest edge. + /// + public static Point FindCircumcenter(Point torg, Point tdest, Point tapex, + ref double xi, ref double eta) + { + double xdo, ydo, xao, yao; + double dodist, aodist, dadist; + double denominator; + double dx, dy; + + Statistic.CircumcenterCount++; + + // Compute the circumcenter of the triangle. + xdo = tdest.x - torg.x; + ydo = tdest.y - torg.y; + xao = tapex.x - torg.x; + yao = tapex.y - torg.y; + dodist = xdo * xdo + ydo * ydo; + aodist = xao * xao + yao * yao; + dadist = (tdest.x - tapex.x) * (tdest.x - tapex.x) + + (tdest.y - tapex.y) * (tdest.y - tapex.y); + + if (Behavior.NoExact) + { + denominator = 0.5 / (xdo * yao - xao * ydo); + } + else + { + // Use the counterclockwise() routine to ensure a positive (and + // reasonably accurate) result, avoiding any possibility of + // division by zero. + denominator = 0.5 / CounterClockwise(tdest, tapex, torg); + // Don't count the above as an orientation test. + Statistic.CounterClockwiseCount--; + } + + dx = (yao * dodist - ydo * aodist) * denominator; + dy = (xdo * aodist - xao * dodist) * denominator; // To interpolate vertex attributes for the new vertex inserted at // the circumcenter, define a coordinate system with a xi-axis, @@ -372,7 +427,7 @@ namespace TriangleNet xi = (yao * dx - xao * dy) * (2.0 * denominator); eta = (xdo * dy - ydo * dx) * (2.0 * denominator); - return circumcenter; + return new Point(torg.x + dx, torg.y + dy); } } } diff --git a/Triangle.NET/Triangle/Quality.cs b/Triangle.NET/Triangle/Quality.cs index dab80cf..1ae3be1 100644 --- a/Triangle.NET/Triangle/Quality.cs +++ b/Triangle.NET/Triangle/Quality.cs @@ -21,19 +21,26 @@ namespace TriangleNet Queue badsubsegs; BadTriQueue queue; Mesh mesh; + Behavior behavior; + + NewLocation newLocation; // Not used at the moment Func userTest; ILog logger; - public Quality(Mesh m) + public Quality(Mesh mesh) { logger = SimpleLog.Instance; badsubsegs = new Queue(); queue = new BadTriQueue(); - mesh = m; + + this.mesh = mesh; + this.behavior = mesh.behavior; + + newLocation = new NewLocation(mesh); } /// @@ -278,9 +285,9 @@ namespace TriangleNet (eorg.y - eapex.y) * (edest.y - eapex.y); if (dotproduct < 0.0) { - if (Behavior.ConformDel || + if (behavior.ConformDel || (dotproduct * dotproduct >= - (2.0 * Behavior.GoodAngle - 1.0) * (2.0 * Behavior.GoodAngle - 1.0) * + (2.0 * behavior.GoodAngle - 1.0) * (2.0 * behavior.GoodAngle - 1.0) * ((eorg.x - eapex.x) * (eorg.x - eapex.x) + (eorg.y - eapex.y) * (eorg.y - eapex.y)) * ((edest.x - eapex.x) * (edest.x - eapex.x) + @@ -305,9 +312,9 @@ namespace TriangleNet (eorg.y - eapex.y) * (edest.y - eapex.y); if (dotproduct < 0.0) { - if (Behavior.ConformDel || + if (behavior.ConformDel || (dotproduct * dotproduct >= - (2.0 * Behavior.GoodAngle - 1.0) * (2.0 * Behavior.GoodAngle - 1.0) * + (2.0 * behavior.GoodAngle - 1.0) * (2.0 * behavior.GoodAngle - 1.0) * ((eorg.x - eapex.x) * (eorg.x - eapex.x) + (eorg.y - eapex.y) * (eorg.y - eapex.y)) * ((edest.x - eapex.x) * (edest.x - eapex.x) + @@ -318,7 +325,7 @@ namespace TriangleNet } } - if (encroached > 0 && (Behavior.NoBisect == 0 || ((Behavior.NoBisect == 1) && (sides == 2)))) + if (encroached > 0 && (behavior.NoBisect == 0 || ((behavior.NoBisect == 1) && (sides == 2)))) { // Add the subsegment to the list of encroached subsegments. // Be sure to get the orientation right. @@ -422,11 +429,11 @@ namespace TriangleNet testtri.Lprev(ref tri1); } - if (Behavior.VarArea || Behavior.FixedArea || Behavior.Usertest) + if (behavior.VarArea || behavior.FixedArea || behavior.Usertest) { // Check whether the area is larger than permitted. area = 0.5 * (dxod * dyda - dyod * dxda); - if (Behavior.FixedArea && (area > Behavior.MaxArea)) + if (behavior.FixedArea && (area > behavior.MaxArea)) { // Add this triangle to the list of bad triangles. queue.Enqueue(ref testtri, minedge, tapex, torg, tdest); @@ -434,7 +441,7 @@ namespace TriangleNet } // Nonpositive area constraints are treated as unconstrained. - if ((Behavior.VarArea) && (area > testtri.triangle.area) && (testtri.triangle.area > 0.0)) + if ((behavior.VarArea) && (area > testtri.triangle.area) && (testtri.triangle.area > 0.0)) { // Add this triangle to the list of bad triangles. queue.Enqueue(ref testtri, minedge, tapex, torg, tdest); @@ -442,7 +449,7 @@ namespace TriangleNet } // Check whether the user thinks this triangle is too large. - if (Behavior.Usertest && userTest != null) + if (behavior.Usertest && userTest != null) { if (userTest(torg, tdest, tapex, area)) { @@ -476,7 +483,7 @@ namespace TriangleNet } // Check whether the angle is smaller than permitted. - if ((angle > Behavior.GoodAngle) || (maxangle < Behavior.MaxGoodAngle && Behavior.MaxAngle != 0.0)) + if ((angle > behavior.GoodAngle) || (maxangle < behavior.MaxGoodAngle && behavior.MaxAngle != 0.0)) { // Use the rules of Miller, Pav, and Walkington to decide that certain // triangles should not be split, even if they have bad angles. @@ -643,7 +650,7 @@ namespace TriangleNet // If we're using Chew's algorithm (rather than Ruppert's) // to define encroachment, delete free vertices from the // subsegment's diametral circle. - if (!Behavior.ConformDel && !acuteorg && !acutedest) + if (!behavior.ConformDel && !acuteorg && !acutedest) { eapex = enctri.Apex(); while ((eapex.type == VertexType.FreeVertex) && @@ -673,7 +680,7 @@ namespace TriangleNet acuteorg = acuteorg || acuteorg2; // Delete free vertices from the subsegment's diametral circle. - if (!Behavior.ConformDel && !acuteorg2 && !acutedest2) + if (!behavior.ConformDel && !acuteorg2 && !acutedest2) { eapex = testtri.Org(); while ((eapex.type == VertexType.FreeVertex) && @@ -845,15 +852,15 @@ namespace TriangleNet // for mesh refinement. // TODO: NewLocation doesn't work for refinement. Why? Maybe // reset VertexType? - if (Behavior.FixedArea || Behavior.VarArea) + if (behavior.FixedArea || behavior.VarArea) { - Point tmp = Primitives.FindCircumcenter(borg, bdest, bapex, ref xi, ref eta, true); + Point tmp = Primitives.FindCircumcenter(borg, bdest, bapex, ref xi, ref eta, behavior.Offconstant); newvertex.x = tmp.x; newvertex.y = tmp.y; } else { - NewLocation.FindLocation(mesh, borg, bdest, bapex, newvertex, ref xi, ref eta, true, badotri); + newLocation.FindLocation(borg, bdest, bapex, newvertex, ref xi, ref eta, true, badotri); } // Check whether the new vertex lies on a triangle vertex. @@ -959,7 +966,7 @@ namespace TriangleNet // triangulation should be (conforming) Delaunay. // Next, we worry about enforcing triangle quality. - if ((Behavior.MinAngle > 0.0) || Behavior.VarArea || Behavior.FixedArea || Behavior.Usertest) + if ((behavior.MinAngle > 0.0) || behavior.VarArea || behavior.FixedArea || behavior.Usertest) { // TODO: Reset queue? (Or is it always empty at this point) @@ -989,7 +996,7 @@ namespace TriangleNet // and have no low-quality triangles. // Might we have run out of Steiner points too soon? - if (Behavior.Verbose && Behavior.ConformDel && (badsubsegs.Count > 0) && (mesh.steinerleft == 0)) + if (Behavior.Verbose && behavior.ConformDel && (badsubsegs.Count > 0) && (mesh.steinerleft == 0)) { logger.Warning("I ran out of Steiner points, but the mesh has encroached subsegments, " diff --git a/Triangle.NET/Triangle/Tools/BoundedVoronoi.cs b/Triangle.NET/Triangle/Tools/BoundedVoronoi.cs index 7301e91..2e46eec 100644 --- a/Triangle.NET/Triangle/Tools/BoundedVoronoi.cs +++ b/Triangle.NET/Triangle/Tools/BoundedVoronoi.cs @@ -162,7 +162,7 @@ namespace TriangleNet.Tools Vertex tdest = tri.Dest(); Vertex tapex = tri.Apex(); - cc = Primitives.FindCircumcenter(torg, tdest, tapex, ref xi, ref eta, false); + cc = Primitives.FindCircumcenter(torg, tdest, tapex, ref xi, ref eta); if (SegmentsIntersect(ref seg, cc, torg, ref pt, true)) { @@ -217,12 +217,12 @@ namespace TriangleNet.Tools torg = f.Org(); tdest = f.Dest(); tapex = f.Apex(); - cc_f = Primitives.FindCircumcenter(torg, tdest, tapex, ref xi, ref eta, false); + 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, false); + cc_f_next = Primitives.FindCircumcenter(torg, tdest, tapex, ref xi, ref eta); // if f is tagged non-blind then if (!f.triangle.infected) @@ -353,7 +353,7 @@ namespace TriangleNet.Tools torg = f.Org(); tdest = f.Dest(); tapex = f.Apex(); - cc_f = Primitives.FindCircumcenter(torg, tdest, tapex, ref xi, ref eta, false); + cc_f = Primitives.FindCircumcenter(torg, tdest, tapex, ref xi, ref eta); if (f_next.triangle == Mesh.dummytri) { @@ -375,7 +375,7 @@ namespace TriangleNet.Tools 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); + cc_f_next = Primitives.FindCircumcenter(torg, tdest, tapex, ref xi, ref eta); // if f is tagged non-blind then if (!f.triangle.infected) diff --git a/Triangle.NET/Triangle/Tools/Voronoi.cs b/Triangle.NET/Triangle/Tools/Voronoi.cs index 46fd671..435c992 100644 --- a/Triangle.NET/Triangle/Tools/Voronoi.cs +++ b/Triangle.NET/Triangle/Tools/Voronoi.cs @@ -87,7 +87,7 @@ namespace TriangleNet.Tools torg = tri.Org(); tdest = tri.Dest(); tapex = tri.Apex(); - circumcenter = Primitives.FindCircumcenter(torg, tdest, tapex, ref xi, ref eta, false); + circumcenter = Primitives.FindCircumcenter(torg, tdest, tapex, ref xi, ref eta); // X and y coordinates. this.points[i] = new Point(circumcenter.x, circumcenter.y);