diff --git a/Triangle.NET/TestApp/FormMain.cs b/Triangle.NET/TestApp/FormMain.cs index 0ab4a8c..2681602 100644 --- a/Triangle.NET/TestApp/FormMain.cs +++ b/Triangle.NET/TestApp/FormMain.cs @@ -728,11 +728,14 @@ namespace MeshExplorer { if (mesh != null) { - bool isConsistent, isDelaunay; + bool save = Behavior.Verbose; Behavior.Verbose = true; - mesh.Check(out isConsistent, out isDelaunay); - Behavior.Verbose = false; + + bool isConsistent = MeshValidator.IsConsistent(mesh); + bool isDelaunay = MeshValidator.IsDelaunay(mesh); + + Behavior.Verbose = save; ShowLog(); } diff --git a/Triangle.NET/TestApp/Generators/RingPolygon.cs b/Triangle.NET/TestApp/Generators/RingPolygon.cs index b508452..dea4d7f 100644 --- a/Triangle.NET/TestApp/Generators/RingPolygon.cs +++ b/Triangle.NET/TestApp/Generators/RingPolygon.cs @@ -61,7 +61,7 @@ namespace MeshExplorer.Generators for (int i = 0; i < m; i++) { input.AddPoint(r * Math.Cos(i * step), r * Math.Sin(i * step)); - input.AddSegment(i, (i + 1) % m); + input.AddSegment(i, (i + 1) % m, 1); } r = 1.5 * r; @@ -81,7 +81,7 @@ namespace MeshExplorer.Generators } input.AddPoint(ro * Math.Cos(i * step + offset), ro * Math.Sin(i * step + offset)); - input.AddSegment(m + i, m + ((i + 1) % n)); + input.AddSegment(m + i, m + ((i + 1) % n), 2); } input.AddHole(0, 0); diff --git a/Triangle.NET/TestApp/IO/GeometryWriter.cs b/Triangle.NET/TestApp/IO/GeometryWriter.cs index c416866..79746e1 100644 --- a/Triangle.NET/TestApp/IO/GeometryWriter.cs +++ b/Triangle.NET/TestApp/IO/GeometryWriter.cs @@ -6,10 +6,9 @@ namespace MeshExplorer.IO { - using System; using System.Collections.Generic; - using System.Linq; using System.IO; + using System.Linq; using TriangleNet.Geometry; /// @@ -17,8 +16,18 @@ namespace MeshExplorer.IO /// public static class GeometryWriter { - public static void Write(InputGeometry geometry, string filename) + private static int OFFSET = 0; + + /// + /// Writes an InputGeometry to a Triangle format .poly file. + /// + /// The InputGeometry to write. + /// The filename. + /// If true, indices will start at 1 (compatible with original C code). + public static void Write(InputGeometry geometry, string filename, bool compatibleMode = false) { + OFFSET = compatibleMode ? 1 : 0; + using (StreamWriter writer = new StreamWriter(filename)) { WritePoints(writer, geometry.Points, geometry.Count); @@ -29,7 +38,7 @@ namespace MeshExplorer.IO private static void WritePoints(StreamWriter writer, IEnumerable points, int count) { - int attributes = 0, index = 0; + int attributes = 0, index = OFFSET; var first = points.FirstOrDefault(); @@ -60,13 +69,13 @@ namespace MeshExplorer.IO private static void WriteSegments(StreamWriter writer, IEnumerable edges) { - int index = 0; + int index = OFFSET; writer.WriteLine("{0} {1}", edges.Count(), 1); foreach (var item in edges) { - writer.WriteLine("{0} {1} {2} {3}", index, item.P0, item.P1, item.Boundary); + writer.WriteLine("{0} {1} {2} {3}", index, item.P0 + OFFSET, item.P1 + OFFSET, item.Boundary); index++; } @@ -74,7 +83,7 @@ namespace MeshExplorer.IO private static void WriteHoles(StreamWriter writer, IEnumerable holes) { - int index = 0; + int index = OFFSET; writer.WriteLine("{0}", holes.Count()); diff --git a/Triangle.NET/Triangle/Behavior.cs b/Triangle.NET/Triangle/Behavior.cs index 8d4c87a..cae622f 100644 --- a/Triangle.NET/Triangle/Behavior.cs +++ b/Triangle.NET/Triangle/Behavior.cs @@ -8,6 +8,7 @@ namespace TriangleNet { using System; + using TriangleNet.Geometry; using TriangleNet.Log; /// @@ -20,7 +21,6 @@ namespace TriangleNet bool poly = false; bool quality = false; bool varArea = false; - bool usertest = false; bool convex = false; bool jettison = false; bool boundaryMarkers = true; @@ -28,6 +28,8 @@ namespace TriangleNet bool conformDel = false; TriangulationAlgorithm algorithm = TriangulationAlgorithm.Dwyer; + Func usertest; + int noBisect = 0; int steiner = -1; @@ -181,7 +183,7 @@ namespace TriangleNet /// /// Apply a user-defined triangle constraint. /// - public bool Usertest + public Func UserTest { get { return usertest; } set { usertest = value; } diff --git a/Triangle.NET/Triangle/ConstraintMesher.cs b/Triangle.NET/Triangle/ConstraintMesher.cs new file mode 100644 index 0000000..13ec654 --- /dev/null +++ b/Triangle.NET/Triangle/ConstraintMesher.cs @@ -0,0 +1,776 @@ +// ----------------------------------------------------------------------- +// +// Original Triangle code by Jonathan Richard Shewchuk, http://www.cs.cmu.edu/~quake/triangle.html +// Triangle.NET code by Christian Woltering, http://triangle.codeplex.com/ +// +// ----------------------------------------------------------------------- + +namespace TriangleNet +{ + using System; + using TriangleNet.Data; + using TriangleNet.Geometry; + using TriangleNet.Log; + + internal class ConstraintMesher + { + Mesh mesh; + Behavior behavior; + TriangleLocator locator; + + ILog logger; + + public ConstraintMesher(Mesh mesh) + { + this.mesh = mesh; + this.behavior = mesh.behavior; + this.locator = mesh.locator; + + logger = SimpleLog.Instance; + } + + /// + /// Create the segments of a triangulation, including PSLG segments and edges + /// on the convex hull. + /// + /// + /// + /// + public void FormSkeleton(InputGeometry input) + { + Vertex endpoint1, endpoint2; + int end1, end2; + int boundmarker; + + mesh.insegments = 0; + + if (behavior.Poly) + { + // If the input vertices are collinear, there is no triangulation, + // so don't try to insert segments. + if (mesh.triangles.Count == 0) + { + return; + } + + // If segments are to be inserted, compute a mapping + // from vertices to triangles. + if (input.HasSegments) + { + mesh.MakeVertexMap(); + } + + boundmarker = 0; + + // Read and insert the segments. + foreach (var seg in input.segments) + { + mesh.insegments++; + + end1 = seg.P0; + end2 = seg.P1; + boundmarker = seg.Boundary; + + if ((end1 < 0) || (end1 >= mesh.invertices)) + { + if (Behavior.Verbose) + { + logger.Warning("Invalid first endpoint of segment.", "Mesh.FormSkeleton().1"); + } + } + else if ((end2 < 0) || (end2 >= mesh.invertices)) + { + if (Behavior.Verbose) + { + logger.Warning("Invalid second endpoint of segment.", "Mesh.FormSkeleton().2"); + } + } + else + { + // TODO: Is using the vertex ID reliable??? + // It should be. The ID gets appropriately set in TransferNodes(). + + // Find the vertices numbered 'end1' and 'end2'. + endpoint1 = mesh.vertices[end1]; + endpoint2 = mesh.vertices[end2]; + if ((endpoint1.x == endpoint2.x) && (endpoint1.y == endpoint2.y)) + { + if (Behavior.Verbose) + { + logger.Warning("Endpoints of segment (IDs " + end1 + "/" + end2 + ") are coincident.", + "Mesh.FormSkeleton()"); + } + } + else + { + InsertSegment(endpoint1, endpoint2, boundmarker); + } + } + } + } + + if (behavior.Convex || !behavior.Poly) + { + // Enclose the convex hull with subsegments. + MarkHull(); + } + } + + #region Segment insertion + + /// + /// Find the first triangle on the path from one point to another. + /// + /// + /// + /// + /// The return value notes whether the destination or apex of the found + /// triangle is collinear with the two points in question. + /// + /// Finds the triangle that intersects a line segment drawn from the + /// origin of 'searchtri' to the point 'searchpoint', and returns the result + /// in 'searchtri'. The origin of 'searchtri' does not change, even though + /// the triangle returned may differ from the one passed in. This routine + /// is used to find the direction to move in to get from one point to + /// another. + /// + private FindDirectionResult FindDirection(ref Otri searchtri, Vertex searchpoint) + { + Otri checktri = default(Otri); + Vertex startvertex; + Vertex leftvertex, rightvertex; + double leftccw, rightccw; + bool leftflag, rightflag; + + startvertex = searchtri.Org(); + rightvertex = searchtri.Dest(); + leftvertex = searchtri.Apex(); + // Is 'searchpoint' to the left? + leftccw = Primitives.CounterClockwise(searchpoint, startvertex, leftvertex); + leftflag = leftccw > 0.0; + // Is 'searchpoint' to the right? + rightccw = Primitives.CounterClockwise(startvertex, searchpoint, rightvertex); + rightflag = rightccw > 0.0; + if (leftflag && rightflag) + { + // 'searchtri' faces directly away from 'searchpoint'. We could go left + // or right. Ask whether it's a triangle or a boundary on the left. + searchtri.Onext(ref checktri); + if (checktri.triangle == Mesh.dummytri) + { + leftflag = false; + } + else + { + rightflag = false; + } + } + while (leftflag) + { + // Turn left until satisfied. + searchtri.OnextSelf(); + if (searchtri.triangle == Mesh.dummytri) + { + logger.Error("Unable to find a triangle on path.", "Mesh.FindDirection().1"); + throw new Exception("Unable to find a triangle on path."); + } + leftvertex = searchtri.Apex(); + rightccw = leftccw; + leftccw = Primitives.CounterClockwise(searchpoint, startvertex, leftvertex); + leftflag = leftccw > 0.0; + } + while (rightflag) + { + // Turn right until satisfied. + searchtri.OprevSelf(); + if (searchtri.triangle == Mesh.dummytri) + { + logger.Error("Unable to find a triangle on path.", "Mesh.FindDirection().2"); + throw new Exception("Unable to find a triangle on path."); + } + rightvertex = searchtri.Dest(); + leftccw = rightccw; + rightccw = Primitives.CounterClockwise(startvertex, searchpoint, rightvertex); + rightflag = rightccw > 0.0; + } + if (leftccw == 0.0) + { + return FindDirectionResult.Leftcollinear; + } + else if (rightccw == 0.0) + { + return FindDirectionResult.Rightcollinear; + } + else + { + return FindDirectionResult.Within; + } + } + + /// + /// Find the intersection of an existing segment and a segment that is being + /// inserted. Insert a vertex at the intersection, splitting an existing subsegment. + /// + /// + /// + /// + /// + /// The segment being inserted connects the apex of splittri to endpoint2. + /// splitsubseg is the subsegment being split, and MUST adjoin splittri. + /// Hence, endpoints of the subsegment being split are the origin and + /// destination of splittri. + /// + /// On completion, splittri is a handle having the newly inserted + /// intersection point as its origin, and endpoint1 as its destination. + /// + private void SegmentIntersection(ref Otri splittri, ref Osub splitsubseg, Vertex endpoint2) + { + Osub opposubseg = default(Osub); + Vertex endpoint1; + Vertex torg, tdest; + Vertex leftvertex, rightvertex; + Vertex newvertex; + InsertVertexResult success; + + double ex, ey; + double tx, ty; + double etx, ety; + double split, denom; + + // Find the other three segment endpoints. + endpoint1 = splittri.Apex(); + torg = splittri.Org(); + tdest = splittri.Dest(); + // Segment intersection formulae; see the Antonio reference. + tx = tdest.x - torg.x; + ty = tdest.y - torg.y; + ex = endpoint2.x - endpoint1.x; + ey = endpoint2.y - endpoint1.y; + etx = torg.x - endpoint2.x; + ety = torg.y - endpoint2.y; + denom = ty * ex - tx * ey; + if (denom == 0.0) + { + logger.Error("Attempt to find intersection of parallel segments.", + "Mesh.SegmentIntersection()"); + throw new Exception("Attempt to find intersection of parallel segments."); + } + split = (ey * etx - ex * ety) / denom; + + // Create the new vertex. + newvertex = new Vertex( + torg.x + split * (tdest.x - torg.x), + torg.y + split * (tdest.y - torg.y), + splitsubseg.seg.boundary, + mesh.nextras); + + newvertex.hash = mesh.hash_vtx++; + newvertex.id = newvertex.hash; + + // Interpolate its attributes. + for (int i = 0; i < mesh.nextras; i++) + { + newvertex.attributes[i] = torg.attributes[i] + split * (tdest.attributes[i] - torg.attributes[i]); + } + + mesh.vertices.Add(newvertex.hash, newvertex); + + // Insert the intersection vertex. This should always succeed. + success = mesh.InsertVertex(newvertex, ref splittri, ref splitsubseg, false, false); + if (success != InsertVertexResult.Successful) + { + logger.Error("Failure to split a segment.", "Mesh.SegmentIntersection()"); + throw new Exception("Failure to split a segment."); + } + // Record a triangle whose origin is the new vertex. + newvertex.tri = splittri; + if (mesh.steinerleft > 0) + { + mesh.steinerleft--; + } + + // Divide the segment into two, and correct the segment endpoints. + splitsubseg.SymSelf(); + splitsubseg.Pivot(ref opposubseg); + splitsubseg.Dissolve(); + opposubseg.Dissolve(); + do + { + splitsubseg.SetSegOrg(newvertex); + splitsubseg.NextSelf(); + } while (splitsubseg.seg != Mesh.dummysub); + do + { + opposubseg.SetSegOrg(newvertex); + opposubseg.NextSelf(); + } while (opposubseg.seg != Mesh.dummysub); + + // Inserting the vertex may have caused edge flips. We wish to rediscover + // the edge connecting endpoint1 to the new intersection vertex. + FindDirection(ref splittri, endpoint1); + + rightvertex = splittri.Dest(); + leftvertex = splittri.Apex(); + if ((leftvertex.x == endpoint1.x) && (leftvertex.y == endpoint1.y)) + { + splittri.OnextSelf(); + } + else if ((rightvertex.x != endpoint1.x) || (rightvertex.y != endpoint1.y)) + { + logger.Error("Topological inconsistency after splitting a segment.", "Mesh.SegmentIntersection()"); + throw new Exception("Topological inconsistency after splitting a segment."); + } + // 'splittri' should have destination endpoint1. + } + + /// + /// Scout the first triangle on the path from one endpoint to another, and check + /// for completion (reaching the second endpoint), a collinear vertex, or the + /// intersection of two segments. + /// + /// + /// + /// + /// Returns true if the entire segment is successfully inserted, and false + /// if the job must be finished by ConstrainedEdge(). + /// + /// If the first triangle on the path has the second endpoint as its + /// destination or apex, a subsegment is inserted and the job is done. + /// + /// If the first triangle on the path has a destination or apex that lies on + /// the segment, a subsegment is inserted connecting the first endpoint to + /// the collinear vertex, and the search is continued from the collinear + /// vertex. + /// + /// If the first triangle on the path has a subsegment opposite its origin, + /// then there is a segment that intersects the segment being inserted. + /// Their intersection vertex is inserted, splitting the subsegment. + /// + private bool ScoutSegment(ref Otri searchtri, Vertex endpoint2, int newmark) + { + Otri crosstri = default(Otri); + Osub crosssubseg = default(Osub); + Vertex leftvertex, rightvertex; + FindDirectionResult collinear; + + collinear = FindDirection(ref searchtri, endpoint2); + rightvertex = searchtri.Dest(); + leftvertex = searchtri.Apex(); + if (((leftvertex.x == endpoint2.x) && (leftvertex.y == endpoint2.y)) || + ((rightvertex.x == endpoint2.x) && (rightvertex.y == endpoint2.y))) + { + // The segment is already an edge in the mesh. + if ((leftvertex.x == endpoint2.x) && (leftvertex.y == endpoint2.y)) + { + searchtri.LprevSelf(); + } + // Insert a subsegment, if there isn't already one there. + mesh.InsertSubseg(ref searchtri, newmark); + return true; + } + else if (collinear == FindDirectionResult.Leftcollinear) + { + // We've collided with a vertex between the segment's endpoints. + // Make the collinear vertex be the triangle's origin. + searchtri.LprevSelf(); + mesh.InsertSubseg(ref searchtri, newmark); + // Insert the remainder of the segment. + return ScoutSegment(ref searchtri, endpoint2, newmark); + } + else if (collinear == FindDirectionResult.Rightcollinear) + { + // We've collided with a vertex between the segment's endpoints. + mesh.InsertSubseg(ref searchtri, newmark); + // Make the collinear vertex be the triangle's origin. + searchtri.LnextSelf(); + // Insert the remainder of the segment. + return ScoutSegment(ref searchtri, endpoint2, newmark); + } + else + { + searchtri.Lnext(ref crosstri); + crosstri.SegPivot(ref crosssubseg); + // Check for a crossing segment. + if (crosssubseg.seg == Mesh.dummysub) + { + return false; + } + else + { + // Insert a vertex at the intersection. + SegmentIntersection(ref crosstri, ref crosssubseg, endpoint2); + crosstri.Copy(ref searchtri); + mesh.InsertSubseg(ref searchtri, newmark); + // Insert the remainder of the segment. + return ScoutSegment(ref searchtri, endpoint2, newmark); + } + } + } + + /// + /// Enforce the Delaunay condition at an edge, fanning out recursively from + /// an existing vertex. Pay special attention to stacking inverted triangles. + /// + /// + /// Indicates whether or not fixuptri is to the left of + /// the segment being inserted. (Imagine that the segment is pointing up from + /// endpoint1 to endpoint2.) + /// + /// This is a support routine for inserting segments into a constrained + /// Delaunay triangulation. + /// + /// The origin of fixuptri is treated as if it has just been inserted, and + /// the local Delaunay condition needs to be enforced. It is only enforced + /// in one sector, however, that being the angular range defined by + /// fixuptri. + /// + /// This routine also needs to make decisions regarding the "stacking" of + /// triangles. (Read the description of ConstrainedEdge() below before + /// reading on here, so you understand the algorithm.) If the position of + /// the new vertex (the origin of fixuptri) indicates that the vertex before + /// it on the polygon is a reflex vertex, then "stack" the triangle by + /// doing nothing. (fixuptri is an inverted triangle, which is how stacked + /// triangles are identified.) + /// + /// Otherwise, check whether the vertex before that was a reflex vertex. + /// If so, perform an edge flip, thereby eliminating an inverted triangle + /// (popping it off the stack). The edge flip may result in the creation + /// of a new inverted triangle, depending on whether or not the new vertex + /// is visible to the vertex three edges behind on the polygon. + /// + /// If neither of the two vertices behind the new vertex are reflex + /// vertices, fixuptri and fartri, the triangle opposite it, are not + /// inverted; hence, ensure that the edge between them is locally Delaunay. + /// + private void DelaunayFixup(ref Otri fixuptri, bool leftside) + { + Otri neartri = default(Otri); + Otri fartri = default(Otri); + Osub faredge = default(Osub); + Vertex nearvertex, leftvertex, rightvertex, farvertex; + + fixuptri.Lnext(ref neartri); + neartri.Sym(ref fartri); + // Check if the edge opposite the origin of fixuptri can be flipped. + if (fartri.triangle == Mesh.dummytri) + { + return; + } + neartri.SegPivot(ref faredge); + if (faredge.seg != Mesh.dummysub) + { + return; + } + // Find all the relevant vertices. + nearvertex = neartri.Apex(); + leftvertex = neartri.Org(); + rightvertex = neartri.Dest(); + farvertex = fartri.Apex(); + // Check whether the previous polygon vertex is a reflex vertex. + if (leftside) + { + if (Primitives.CounterClockwise(nearvertex, leftvertex, farvertex) <= 0.0) + { + // leftvertex is a reflex vertex too. Nothing can + // be done until a convex section is found. + return; + } + } + else + { + if (Primitives.CounterClockwise(farvertex, rightvertex, nearvertex) <= 0.0) + { + // rightvertex is a reflex vertex too. Nothing can + // be done until a convex section is found. + return; + } + } + if (Primitives.CounterClockwise(rightvertex, leftvertex, farvertex) > 0.0) + { + // fartri is not an inverted triangle, and farvertex is not a reflex + // vertex. As there are no reflex vertices, fixuptri isn't an + // inverted triangle, either. Hence, test the edge between the + // triangles to ensure it is locally Delaunay. + if (Primitives.InCircle(leftvertex, farvertex, rightvertex, nearvertex) <= 0.0) + { + return; + } + // Not locally Delaunay; go on to an edge flip. + } + // else fartri is inverted; remove it from the stack by flipping. + mesh.Flip(ref neartri); + fixuptri.LprevSelf(); // Restore the origin of fixuptri after the flip. + // Recursively process the two triangles that result from the flip. + DelaunayFixup(ref fixuptri, leftside); + DelaunayFixup(ref fartri, leftside); + } + + /// + /// Force a segment into a constrained Delaunay triangulation by deleting the + /// triangles it intersects, and triangulating the polygons that form on each + /// side of it. + /// + /// + /// + /// + /// + /// Generates a single subsegment connecting 'endpoint1' to 'endpoint2'. + /// The triangle 'starttri' has 'endpoint1' as its origin. 'newmark' is the + /// boundary marker of the segment. + /// + /// To insert a segment, every triangle whose interior intersects the + /// segment is deleted. The union of these deleted triangles is a polygon + /// (which is not necessarily monotone, but is close enough), which is + /// divided into two polygons by the new segment. This routine's task is + /// to generate the Delaunay triangulation of these two polygons. + /// + /// You might think of this routine's behavior as a two-step process. The + /// first step is to walk from endpoint1 to endpoint2, flipping each edge + /// encountered. This step creates a fan of edges connected to endpoint1, + /// including the desired edge to endpoint2. The second step enforces the + /// Delaunay condition on each side of the segment in an incremental manner: + /// proceeding along the polygon from endpoint1 to endpoint2 (this is done + /// independently on each side of the segment), each vertex is "enforced" + /// as if it had just been inserted, but affecting only the previous + /// vertices. The result is the same as if the vertices had been inserted + /// in the order they appear on the polygon, so the result is Delaunay. + /// + /// In truth, ConstrainedEdge() interleaves these two steps. The procedure + /// walks from endpoint1 to endpoint2, and each time an edge is encountered + /// and flipped, the newly exposed vertex (at the far end of the flipped + /// edge) is "enforced" upon the previously flipped edges, usually affecting + /// only one side of the polygon (depending upon which side of the segment + /// the vertex falls on). + /// + /// The algorithm is complicated by the need to handle polygons that are not + /// convex. Although the polygon is not necessarily monotone, it can be + /// triangulated in a manner similar to the stack-based algorithms for + /// monotone polygons. For each reflex vertex (local concavity) of the + /// polygon, there will be an inverted triangle formed by one of the edge + /// flips. (An inverted triangle is one with negative area - that is, its + /// vertices are arranged in clockwise order - and is best thought of as a + /// wrinkle in the fabric of the mesh.) Each inverted triangle can be + /// thought of as a reflex vertex pushed on the stack, waiting to be fixed + /// later. + /// + /// A reflex vertex is popped from the stack when a vertex is inserted that + /// is visible to the reflex vertex. (However, if the vertex behind the + /// reflex vertex is not visible to the reflex vertex, a new inverted + /// triangle will take its place on the stack.) These details are handled + /// by the DelaunayFixup() routine above. + /// + private void ConstrainedEdge(ref Otri starttri, Vertex endpoint2, int newmark) + { + Otri fixuptri = default(Otri), fixuptri2 = default(Otri); + Osub crosssubseg = default(Osub); + Vertex endpoint1; + Vertex farvertex; + double area; + bool collision; + bool done; + + endpoint1 = starttri.Org(); + starttri.Lnext(ref fixuptri); + mesh.Flip(ref fixuptri); + // 'collision' indicates whether we have found a vertex directly + // between endpoint1 and endpoint2. + collision = false; + done = false; + do + { + farvertex = fixuptri.Org(); + // 'farvertex' is the extreme point of the polygon we are "digging" + // to get from endpoint1 to endpoint2. + if ((farvertex.x == endpoint2.x) && (farvertex.y == endpoint2.y)) + { + fixuptri.Oprev(ref fixuptri2); + // Enforce the Delaunay condition around endpoint2. + DelaunayFixup(ref fixuptri, false); + DelaunayFixup(ref fixuptri2, true); + done = true; + } + else + { + // Check whether farvertex is to the left or right of the segment being + // inserted, to decide which edge of fixuptri to dig through next. + area = Primitives.CounterClockwise(endpoint1, endpoint2, farvertex); + if (area == 0.0) + { + // We've collided with a vertex between endpoint1 and endpoint2. + collision = true; + fixuptri.Oprev(ref fixuptri2); + // Enforce the Delaunay condition around farvertex. + DelaunayFixup(ref fixuptri, false); + DelaunayFixup(ref fixuptri2, true); + done = true; + } + else + { + if (area > 0.0) + { + // farvertex is to the left of the segment. + fixuptri.Oprev(ref fixuptri2); + // Enforce the Delaunay condition around farvertex, on the + // left side of the segment only. + DelaunayFixup(ref fixuptri2, true); + // Flip the edge that crosses the segment. After the edge is + // flipped, one of its endpoints is the fan vertex, and the + // destination of fixuptri is the fan vertex. + fixuptri.LprevSelf(); + } + else + { + // farvertex is to the right of the segment. + DelaunayFixup(ref fixuptri, false); + // Flip the edge that crosses the segment. After the edge is + // flipped, one of its endpoints is the fan vertex, and the + // destination of fixuptri is the fan vertex. + fixuptri.OprevSelf(); + } + // Check for two intersecting segments. + fixuptri.SegPivot(ref crosssubseg); + if (crosssubseg.seg == Mesh.dummysub) + { + mesh.Flip(ref fixuptri); // May create inverted triangle at left. + } + else + { + // We've collided with a segment between endpoint1 and endpoint2. + collision = true; + // Insert a vertex at the intersection. + SegmentIntersection(ref fixuptri, ref crosssubseg, endpoint2); + done = true; + } + } + } + } while (!done); + // Insert a subsegment to make the segment permanent. + mesh.InsertSubseg(ref fixuptri, newmark); + // If there was a collision with an interceding vertex, install another + // segment connecting that vertex with endpoint2. + if (collision) + { + // Insert the remainder of the segment. + if (!ScoutSegment(ref fixuptri, endpoint2, newmark)) + { + ConstrainedEdge(ref fixuptri, endpoint2, newmark); + } + } + } + + /// + /// Insert a PSLG segment into a triangulation. + /// + /// + /// + /// + private void InsertSegment(Vertex endpoint1, Vertex endpoint2, int newmark) + { + Otri searchtri1 = default(Otri), searchtri2 = default(Otri); + Vertex checkvertex = null; + + // Find a triangle whose origin is the segment's first endpoint. + searchtri1 = endpoint1.tri; + if (searchtri1.triangle != null) + { + checkvertex = searchtri1.Org(); + } + + if (checkvertex != endpoint1) + { + // Find a boundary triangle to search from. + searchtri1.triangle = Mesh.dummytri; + searchtri1.orient = 0; + searchtri1.SymSelf(); + // Search for the segment's first endpoint by point location. + if (locator.Locate(endpoint1, ref searchtri1) != LocateResult.OnVertex) + { + logger.Error("Unable to locate PSLG vertex in triangulation.", "Mesh.InsertSegment().1"); + throw new Exception("Unable to locate PSLG vertex in triangulation."); + } + } + // Remember this triangle to improve subsequent point location. + locator.Update(ref searchtri1); + + // Scout the beginnings of a path from the first endpoint + // toward the second. + if (ScoutSegment(ref searchtri1, endpoint2, newmark)) + { + // The segment was easily inserted. + return; + } + // The first endpoint may have changed if a collision with an intervening + // vertex on the segment occurred. + endpoint1 = searchtri1.Org(); + + // Find a triangle whose origin is the segment's second endpoint. + checkvertex = null; + searchtri2 = endpoint2.tri; + if (searchtri2.triangle != null) + { + checkvertex = searchtri2.Org(); + } + if (checkvertex != endpoint2) + { + // Find a boundary triangle to search from. + searchtri2.triangle = Mesh.dummytri; + searchtri2.orient = 0; + searchtri2.SymSelf(); + // Search for the segment's second endpoint by point location. + if (locator.Locate(endpoint2, ref searchtri2) != LocateResult.OnVertex) + { + logger.Error("Unable to locate PSLG vertex in triangulation.", "Mesh.InsertSegment().2"); + throw new Exception("Unable to locate PSLG vertex in triangulation."); + } + } + // Remember this triangle to improve subsequent point location. + locator.Update(ref searchtri2); + // Scout the beginnings of a path from the second endpoint + // toward the first. + if (ScoutSegment(ref searchtri2, endpoint1, newmark)) + { + // The segment was easily inserted. + return; + } + // The second endpoint may have changed if a collision with an intervening + // vertex on the segment occurred. + endpoint2 = searchtri2.Org(); + + // Insert the segment directly into the triangulation. + ConstrainedEdge(ref searchtri1, endpoint2, newmark); + } + + /// + /// Cover the convex hull of a triangulation with subsegments. + /// + private void MarkHull() + { + Otri hulltri = default(Otri); + Otri nexttri = default(Otri); + Otri starttri = default(Otri); + + // Find a triangle handle on the hull. + hulltri.triangle = Mesh.dummytri; + hulltri.orient = 0; + hulltri.SymSelf(); + // Remember where we started so we know when to stop. + hulltri.Copy(ref starttri); + // Go once counterclockwise around the convex hull. + do + { + // Create a subsegment if there isn't already one here. + mesh.InsertSubseg(ref hulltri, 1); + // To find the next hull edge, go clockwise around the next vertex. + hulltri.LnextSelf(); + hulltri.Oprev(ref nexttri); + while (nexttri.triangle != Mesh.dummytri) + { + nexttri.Copy(ref hulltri); + hulltri.Oprev(ref nexttri); + } + } while (!hulltri.Equal(starttri)); + } + + #endregion + } +} diff --git a/Triangle.NET/Triangle/Mesh.cs b/Triangle.NET/Triangle/Mesh.cs index daf00a6..c0eaa7f 100644 --- a/Triangle.NET/Triangle/Mesh.cs +++ b/Triangle.NET/Triangle/Mesh.cs @@ -340,8 +340,10 @@ namespace TriangleNet // Segments will be introduced next. checksegments = true; + var mesher = new ConstraintMesher(this); + // Insert PSLG segments and/or convex hull segments. - FormSkeleton(input); + mesher.FormSkeleton(input); } if (behavior.Poly && (triangles.Count > 0)) @@ -545,17 +547,6 @@ namespace TriangleNet } } - /// - /// Check mesh consistency and (constrained) Delaunay property. - /// - /// Value indicating if mesh topology is consistent. - /// Value indicating if mesh is Delaunay. - public void Check(out bool isConsistent, out bool isDelaunay) - { - isConsistent = quality.CheckMesh(); - isDelaunay = quality.CheckDelaunay(); - } - #region Misc /// @@ -1982,750 +1973,6 @@ namespace TriangleNet #endregion - #region Segment insertion - - /// - /// Find the first triangle on the path from one point to another. - /// - /// - /// - /// - /// The return value notes whether the destination or apex of the found - /// triangle is collinear with the two points in question. - /// - /// Finds the triangle that intersects a line segment drawn from the - /// origin of 'searchtri' to the point 'searchpoint', and returns the result - /// in 'searchtri'. The origin of 'searchtri' does not change, even though - /// the triangle returned may differ from the one passed in. This routine - /// is used to find the direction to move in to get from one point to - /// another. - /// - private FindDirectionResult FindDirection(ref Otri searchtri, Vertex searchpoint) - { - Otri checktri = default(Otri); - Vertex startvertex; - Vertex leftvertex, rightvertex; - double leftccw, rightccw; - bool leftflag, rightflag; - - startvertex = searchtri.Org(); - rightvertex = searchtri.Dest(); - leftvertex = searchtri.Apex(); - // Is 'searchpoint' to the left? - leftccw = Primitives.CounterClockwise(searchpoint, startvertex, leftvertex); - leftflag = leftccw > 0.0; - // Is 'searchpoint' to the right? - rightccw = Primitives.CounterClockwise(startvertex, searchpoint, rightvertex); - rightflag = rightccw > 0.0; - if (leftflag && rightflag) - { - // 'searchtri' faces directly away from 'searchpoint'. We could go left - // or right. Ask whether it's a triangle or a boundary on the left. - searchtri.Onext(ref checktri); - if (checktri.triangle == Mesh.dummytri) - { - leftflag = false; - } - else - { - rightflag = false; - } - } - while (leftflag) - { - // Turn left until satisfied. - searchtri.OnextSelf(); - if (searchtri.triangle == Mesh.dummytri) - { - logger.Error("Unable to find a triangle on path.", "Mesh.FindDirection().1"); - throw new Exception("Unable to find a triangle on path."); - } - leftvertex = searchtri.Apex(); - rightccw = leftccw; - leftccw = Primitives.CounterClockwise(searchpoint, startvertex, leftvertex); - leftflag = leftccw > 0.0; - } - while (rightflag) - { - // Turn right until satisfied. - searchtri.OprevSelf(); - if (searchtri.triangle == Mesh.dummytri) - { - logger.Error("Unable to find a triangle on path.", "Mesh.FindDirection().2"); - throw new Exception("Unable to find a triangle on path."); - } - rightvertex = searchtri.Dest(); - leftccw = rightccw; - rightccw = Primitives.CounterClockwise(startvertex, searchpoint, rightvertex); - rightflag = rightccw > 0.0; - } - if (leftccw == 0.0) - { - return FindDirectionResult.Leftcollinear; - } - else if (rightccw == 0.0) - { - return FindDirectionResult.Rightcollinear; - } - else - { - return FindDirectionResult.Within; - } - } - - /// - /// Find the intersection of an existing segment and a segment that is being - /// inserted. Insert a vertex at the intersection, splitting an existing subsegment. - /// - /// - /// - /// - /// - /// The segment being inserted connects the apex of splittri to endpoint2. - /// splitsubseg is the subsegment being split, and MUST adjoin splittri. - /// Hence, endpoints of the subsegment being split are the origin and - /// destination of splittri. - /// - /// On completion, splittri is a handle having the newly inserted - /// intersection point as its origin, and endpoint1 as its destination. - /// - private void SegmentIntersection(ref Otri splittri, ref Osub splitsubseg, Vertex endpoint2) - { - Osub opposubseg = default(Osub); - Vertex endpoint1; - Vertex torg, tdest; - Vertex leftvertex, rightvertex; - Vertex newvertex; - InsertVertexResult success; - - double ex, ey; - double tx, ty; - double etx, ety; - double split, denom; - - // Find the other three segment endpoints. - endpoint1 = splittri.Apex(); - torg = splittri.Org(); - tdest = splittri.Dest(); - // Segment intersection formulae; see the Antonio reference. - tx = tdest.x - torg.x; - ty = tdest.y - torg.y; - ex = endpoint2.x - endpoint1.x; - ey = endpoint2.y - endpoint1.y; - etx = torg.x - endpoint2.x; - ety = torg.y - endpoint2.y; - denom = ty * ex - tx * ey; - if (denom == 0.0) - { - logger.Error("Attempt to find intersection of parallel segments.", - "Mesh.SegmentIntersection()"); - throw new Exception("Attempt to find intersection of parallel segments."); - } - split = (ey * etx - ex * ety) / denom; - - // Create the new vertex. - newvertex = new Vertex( - torg.x + split * (tdest.x - torg.x), - torg.y + split * (tdest.y - torg.y), - splitsubseg.seg.boundary, - this.nextras); - - newvertex.hash = this.hash_vtx++; - newvertex.id = newvertex.hash; - - // Interpolate its attributes. - for (int i = 0; i < nextras; i++) - { - newvertex.attributes[i] = torg.attributes[i] + split * (tdest.attributes[i] - torg.attributes[i]); - } - - vertices.Add(newvertex.hash, newvertex); - - // Insert the intersection vertex. This should always succeed. - success = InsertVertex(newvertex, ref splittri, ref splitsubseg, false, false); - if (success != InsertVertexResult.Successful) - { - logger.Error("Failure to split a segment.", "Mesh.SegmentIntersection()"); - throw new Exception("Failure to split a segment."); - } - // Record a triangle whose origin is the new vertex. - newvertex.tri = splittri; - if (steinerleft > 0) - { - steinerleft--; - } - - // Divide the segment into two, and correct the segment endpoints. - splitsubseg.SymSelf(); - splitsubseg.Pivot(ref opposubseg); - splitsubseg.Dissolve(); - opposubseg.Dissolve(); - do - { - splitsubseg.SetSegOrg(newvertex); - splitsubseg.NextSelf(); - } while (splitsubseg.seg != Mesh.dummysub); - do - { - opposubseg.SetSegOrg(newvertex); - opposubseg.NextSelf(); - } while (opposubseg.seg != Mesh.dummysub); - - // Inserting the vertex may have caused edge flips. We wish to rediscover - // the edge connecting endpoint1 to the new intersection vertex. - FindDirection(ref splittri, endpoint1); - - rightvertex = splittri.Dest(); - leftvertex = splittri.Apex(); - if ((leftvertex.x == endpoint1.x) && (leftvertex.y == endpoint1.y)) - { - splittri.OnextSelf(); - } - else if ((rightvertex.x != endpoint1.x) || (rightvertex.y != endpoint1.y)) - { - logger.Error("Topological inconsistency after splitting a segment.", "Mesh.SegmentIntersection()"); - throw new Exception("Topological inconsistency after splitting a segment."); - } - // 'splittri' should have destination endpoint1. - } - - /// - /// Scout the first triangle on the path from one endpoint to another, and check - /// for completion (reaching the second endpoint), a collinear vertex, or the - /// intersection of two segments. - /// - /// - /// - /// - /// Returns true if the entire segment is successfully inserted, and false - /// if the job must be finished by ConstrainedEdge(). - /// - /// If the first triangle on the path has the second endpoint as its - /// destination or apex, a subsegment is inserted and the job is done. - /// - /// If the first triangle on the path has a destination or apex that lies on - /// the segment, a subsegment is inserted connecting the first endpoint to - /// the collinear vertex, and the search is continued from the collinear - /// vertex. - /// - /// If the first triangle on the path has a subsegment opposite its origin, - /// then there is a segment that intersects the segment being inserted. - /// Their intersection vertex is inserted, splitting the subsegment. - /// - private bool ScoutSegment(ref Otri searchtri, Vertex endpoint2, int newmark) - { - Otri crosstri = default(Otri); - Osub crosssubseg = default(Osub); - Vertex leftvertex, rightvertex; - FindDirectionResult collinear; - - collinear = FindDirection(ref searchtri, endpoint2); - rightvertex = searchtri.Dest(); - leftvertex = searchtri.Apex(); - if (((leftvertex.x == endpoint2.x) && (leftvertex.y == endpoint2.y)) || - ((rightvertex.x == endpoint2.x) && (rightvertex.y == endpoint2.y))) - { - // The segment is already an edge in the mesh. - if ((leftvertex.x == endpoint2.x) && (leftvertex.y == endpoint2.y)) - { - searchtri.LprevSelf(); - } - // Insert a subsegment, if there isn't already one there. - InsertSubseg(ref searchtri, newmark); - return true; - } - else if (collinear == FindDirectionResult.Leftcollinear) - { - // We've collided with a vertex between the segment's endpoints. - // Make the collinear vertex be the triangle's origin. - searchtri.LprevSelf(); - InsertSubseg(ref searchtri, newmark); - // Insert the remainder of the segment. - return ScoutSegment(ref searchtri, endpoint2, newmark); - } - else if (collinear == FindDirectionResult.Rightcollinear) - { - // We've collided with a vertex between the segment's endpoints. - InsertSubseg(ref searchtri, newmark); - // Make the collinear vertex be the triangle's origin. - searchtri.LnextSelf(); - // Insert the remainder of the segment. - return ScoutSegment(ref searchtri, endpoint2, newmark); - } - else - { - searchtri.Lnext(ref crosstri); - crosstri.SegPivot(ref crosssubseg); - // Check for a crossing segment. - if (crosssubseg.seg == Mesh.dummysub) - { - return false; - } - else - { - // Insert a vertex at the intersection. - SegmentIntersection(ref crosstri, ref crosssubseg, endpoint2); - crosstri.Copy(ref searchtri); - InsertSubseg(ref searchtri, newmark); - // Insert the remainder of the segment. - return ScoutSegment(ref searchtri, endpoint2, newmark); - } - } - } - - /// - /// Enforce the Delaunay condition at an edge, fanning out recursively from - /// an existing vertex. Pay special attention to stacking inverted triangles. - /// - /// - /// Indicates whether or not fixuptri is to the left of - /// the segment being inserted. (Imagine that the segment is pointing up from - /// endpoint1 to endpoint2.) - /// - /// This is a support routine for inserting segments into a constrained - /// Delaunay triangulation. - /// - /// The origin of fixuptri is treated as if it has just been inserted, and - /// the local Delaunay condition needs to be enforced. It is only enforced - /// in one sector, however, that being the angular range defined by - /// fixuptri. - /// - /// This routine also needs to make decisions regarding the "stacking" of - /// triangles. (Read the description of ConstrainedEdge() below before - /// reading on here, so you understand the algorithm.) If the position of - /// the new vertex (the origin of fixuptri) indicates that the vertex before - /// it on the polygon is a reflex vertex, then "stack" the triangle by - /// doing nothing. (fixuptri is an inverted triangle, which is how stacked - /// triangles are identified.) - /// - /// Otherwise, check whether the vertex before that was a reflex vertex. - /// If so, perform an edge flip, thereby eliminating an inverted triangle - /// (popping it off the stack). The edge flip may result in the creation - /// of a new inverted triangle, depending on whether or not the new vertex - /// is visible to the vertex three edges behind on the polygon. - /// - /// If neither of the two vertices behind the new vertex are reflex - /// vertices, fixuptri and fartri, the triangle opposite it, are not - /// inverted; hence, ensure that the edge between them is locally Delaunay. - /// - private void DelaunayFixup(ref Otri fixuptri, bool leftside) - { - Otri neartri = default(Otri); - Otri fartri = default(Otri); - Osub faredge = default(Osub); - Vertex nearvertex, leftvertex, rightvertex, farvertex; - - fixuptri.Lnext(ref neartri); - neartri.Sym(ref fartri); - // Check if the edge opposite the origin of fixuptri can be flipped. - if (fartri.triangle == Mesh.dummytri) - { - return; - } - neartri.SegPivot(ref faredge); - if (faredge.seg != Mesh.dummysub) - { - return; - } - // Find all the relevant vertices. - nearvertex = neartri.Apex(); - leftvertex = neartri.Org(); - rightvertex = neartri.Dest(); - farvertex = fartri.Apex(); - // Check whether the previous polygon vertex is a reflex vertex. - if (leftside) - { - if (Primitives.CounterClockwise(nearvertex, leftvertex, farvertex) <= 0.0) - { - // leftvertex is a reflex vertex too. Nothing can - // be done until a convex section is found. - return; - } - } - else - { - if (Primitives.CounterClockwise(farvertex, rightvertex, nearvertex) <= 0.0) - { - // rightvertex is a reflex vertex too. Nothing can - // be done until a convex section is found. - return; - } - } - if (Primitives.CounterClockwise(rightvertex, leftvertex, farvertex) > 0.0) - { - // fartri is not an inverted triangle, and farvertex is not a reflex - // vertex. As there are no reflex vertices, fixuptri isn't an - // inverted triangle, either. Hence, test the edge between the - // triangles to ensure it is locally Delaunay. - if (Primitives.InCircle(leftvertex, farvertex, rightvertex, nearvertex) <= 0.0) - { - return; - } - // Not locally Delaunay; go on to an edge flip. - } - // else fartri is inverted; remove it from the stack by flipping. - Flip(ref neartri); - fixuptri.LprevSelf(); // Restore the origin of fixuptri after the flip. - // Recursively process the two triangles that result from the flip. - DelaunayFixup(ref fixuptri, leftside); - DelaunayFixup(ref fartri, leftside); - } - - /// - /// Force a segment into a constrained Delaunay triangulation by deleting the - /// triangles it intersects, and triangulating the polygons that form on each - /// side of it. - /// - /// - /// - /// - /// - /// Generates a single subsegment connecting 'endpoint1' to 'endpoint2'. - /// The triangle 'starttri' has 'endpoint1' as its origin. 'newmark' is the - /// boundary marker of the segment. - /// - /// To insert a segment, every triangle whose interior intersects the - /// segment is deleted. The union of these deleted triangles is a polygon - /// (which is not necessarily monotone, but is close enough), which is - /// divided into two polygons by the new segment. This routine's task is - /// to generate the Delaunay triangulation of these two polygons. - /// - /// You might think of this routine's behavior as a two-step process. The - /// first step is to walk from endpoint1 to endpoint2, flipping each edge - /// encountered. This step creates a fan of edges connected to endpoint1, - /// including the desired edge to endpoint2. The second step enforces the - /// Delaunay condition on each side of the segment in an incremental manner: - /// proceeding along the polygon from endpoint1 to endpoint2 (this is done - /// independently on each side of the segment), each vertex is "enforced" - /// as if it had just been inserted, but affecting only the previous - /// vertices. The result is the same as if the vertices had been inserted - /// in the order they appear on the polygon, so the result is Delaunay. - /// - /// In truth, ConstrainedEdge() interleaves these two steps. The procedure - /// walks from endpoint1 to endpoint2, and each time an edge is encountered - /// and flipped, the newly exposed vertex (at the far end of the flipped - /// edge) is "enforced" upon the previously flipped edges, usually affecting - /// only one side of the polygon (depending upon which side of the segment - /// the vertex falls on). - /// - /// The algorithm is complicated by the need to handle polygons that are not - /// convex. Although the polygon is not necessarily monotone, it can be - /// triangulated in a manner similar to the stack-based algorithms for - /// monotone polygons. For each reflex vertex (local concavity) of the - /// polygon, there will be an inverted triangle formed by one of the edge - /// flips. (An inverted triangle is one with negative area - that is, its - /// vertices are arranged in clockwise order - and is best thought of as a - /// wrinkle in the fabric of the mesh.) Each inverted triangle can be - /// thought of as a reflex vertex pushed on the stack, waiting to be fixed - /// later. - /// - /// A reflex vertex is popped from the stack when a vertex is inserted that - /// is visible to the reflex vertex. (However, if the vertex behind the - /// reflex vertex is not visible to the reflex vertex, a new inverted - /// triangle will take its place on the stack.) These details are handled - /// by the DelaunayFixup() routine above. - /// - private void ConstrainedEdge(ref Otri starttri, Vertex endpoint2, int newmark) - { - Otri fixuptri = default(Otri), fixuptri2 = default(Otri); - Osub crosssubseg = default(Osub); - Vertex endpoint1; - Vertex farvertex; - double area; - bool collision; - bool done; - - endpoint1 = starttri.Org(); - starttri.Lnext(ref fixuptri); - Flip(ref fixuptri); - // 'collision' indicates whether we have found a vertex directly - // between endpoint1 and endpoint2. - collision = false; - done = false; - do - { - farvertex = fixuptri.Org(); - // 'farvertex' is the extreme point of the polygon we are "digging" - // to get from endpoint1 to endpoint2. - if ((farvertex.x == endpoint2.x) && (farvertex.y == endpoint2.y)) - { - fixuptri.Oprev(ref fixuptri2); - // Enforce the Delaunay condition around endpoint2. - DelaunayFixup(ref fixuptri, false); - DelaunayFixup(ref fixuptri2, true); - done = true; - } - else - { - // Check whether farvertex is to the left or right of the segment being - // inserted, to decide which edge of fixuptri to dig through next. - area = Primitives.CounterClockwise(endpoint1, endpoint2, farvertex); - if (area == 0.0) - { - // We've collided with a vertex between endpoint1 and endpoint2. - collision = true; - fixuptri.Oprev(ref fixuptri2); - // Enforce the Delaunay condition around farvertex. - DelaunayFixup(ref fixuptri, false); - DelaunayFixup(ref fixuptri2, true); - done = true; - } - else - { - if (area > 0.0) - { - // farvertex is to the left of the segment. - fixuptri.Oprev(ref fixuptri2); - // Enforce the Delaunay condition around farvertex, on the - // left side of the segment only. - DelaunayFixup(ref fixuptri2, true); - // Flip the edge that crosses the segment. After the edge is - // flipped, one of its endpoints is the fan vertex, and the - // destination of fixuptri is the fan vertex. - fixuptri.LprevSelf(); - } - else - { - // farvertex is to the right of the segment. - DelaunayFixup(ref fixuptri, false); - // Flip the edge that crosses the segment. After the edge is - // flipped, one of its endpoints is the fan vertex, and the - // destination of fixuptri is the fan vertex. - fixuptri.OprevSelf(); - } - // Check for two intersecting segments. - fixuptri.SegPivot(ref crosssubseg); - if (crosssubseg.seg == Mesh.dummysub) - { - Flip(ref fixuptri); // May create inverted triangle at left. - } - else - { - // We've collided with a segment between endpoint1 and endpoint2. - collision = true; - // Insert a vertex at the intersection. - SegmentIntersection(ref fixuptri, ref crosssubseg, endpoint2); - done = true; - } - } - } - } while (!done); - // Insert a subsegment to make the segment permanent. - InsertSubseg(ref fixuptri, newmark); - // If there was a collision with an interceding vertex, install another - // segment connecting that vertex with endpoint2. - if (collision) - { - // Insert the remainder of the segment. - if (!ScoutSegment(ref fixuptri, endpoint2, newmark)) - { - ConstrainedEdge(ref fixuptri, endpoint2, newmark); - } - } - } - - /// - /// Insert a PSLG segment into a triangulation. - /// - /// - /// - /// - private void InsertSegment(Vertex endpoint1, Vertex endpoint2, int newmark) - { - Otri searchtri1 = default(Otri), searchtri2 = default(Otri); - Vertex checkvertex = null; - - // Find a triangle whose origin is the segment's first endpoint. - searchtri1 = endpoint1.tri; - if (searchtri1.triangle != null) - { - checkvertex = searchtri1.Org(); - } - - if (checkvertex != endpoint1) - { - // Find a boundary triangle to search from. - searchtri1.triangle = Mesh.dummytri; - searchtri1.orient = 0; - searchtri1.SymSelf(); - // Search for the segment's first endpoint by point location. - if (locator.Locate(endpoint1, ref searchtri1) != LocateResult.OnVertex) - { - logger.Error("Unable to locate PSLG vertex in triangulation.", "Mesh.InsertSegment().1"); - throw new Exception("Unable to locate PSLG vertex in triangulation."); - } - } - // Remember this triangle to improve subsequent point location. - locator.Update(ref searchtri1); - - // Scout the beginnings of a path from the first endpoint - // toward the second. - if (ScoutSegment(ref searchtri1, endpoint2, newmark)) - { - // The segment was easily inserted. - return; - } - // The first endpoint may have changed if a collision with an intervening - // vertex on the segment occurred. - endpoint1 = searchtri1.Org(); - - // Find a triangle whose origin is the segment's second endpoint. - checkvertex = null; - searchtri2 = endpoint2.tri; - if (searchtri2.triangle != null) - { - checkvertex = searchtri2.Org(); - } - if (checkvertex != endpoint2) - { - // Find a boundary triangle to search from. - searchtri2.triangle = Mesh.dummytri; - searchtri2.orient = 0; - searchtri2.SymSelf(); - // Search for the segment's second endpoint by point location. - if (locator.Locate(endpoint2, ref searchtri2) != LocateResult.OnVertex) - { - logger.Error("Unable to locate PSLG vertex in triangulation.", "Mesh.InsertSegment().2"); - throw new Exception("Unable to locate PSLG vertex in triangulation."); - } - } - // Remember this triangle to improve subsequent point location. - locator.Update(ref searchtri2); - // Scout the beginnings of a path from the second endpoint - // toward the first. - if (ScoutSegment(ref searchtri2, endpoint1, newmark)) - { - // The segment was easily inserted. - return; - } - // The second endpoint may have changed if a collision with an intervening - // vertex on the segment occurred. - endpoint2 = searchtri2.Org(); - - // Insert the segment directly into the triangulation. - ConstrainedEdge(ref searchtri1, endpoint2, newmark); - } - - /// - /// Cover the convex hull of a triangulation with subsegments. - /// - private void MarkHull() - { - Otri hulltri = default(Otri); - Otri nexttri = default(Otri); - Otri starttri = default(Otri); - - // Find a triangle handle on the hull. - hulltri.triangle = Mesh.dummytri; - hulltri.orient = 0; - hulltri.SymSelf(); - // Remember where we started so we know when to stop. - hulltri.Copy(ref starttri); - // Go once counterclockwise around the convex hull. - do - { - // Create a subsegment if there isn't already one here. - InsertSubseg(ref hulltri, 1); - // To find the next hull edge, go clockwise around the next vertex. - hulltri.LnextSelf(); - hulltri.Oprev(ref nexttri); - while (nexttri.triangle != Mesh.dummytri) - { - nexttri.Copy(ref hulltri); - hulltri.Oprev(ref nexttri); - } - } while (!hulltri.Equal(starttri)); - } - - /// - /// Create the segments of a triangulation, including PSLG segments and edges - /// on the convex hull. - /// - /// - /// - /// - private void FormSkeleton(InputGeometry input) - { - Vertex endpoint1, endpoint2; - int end1, end2; - int boundmarker; - - this.insegments = 0; - - if (behavior.Poly) - { - // If the input vertices are collinear, there is no triangulation, - // so don't try to insert segments. - if (triangles.Count == 0) - { - return; - } - - // If segments are to be inserted, compute a mapping - // from vertices to triangles. - if (input.HasSegments) - { - MakeVertexMap(); - } - - boundmarker = 0; - - // Read and insert the segments. - foreach (var seg in input.segments) - { - this.insegments++; - - end1 = seg.P0; - end2 = seg.P1; - boundmarker = seg.Boundary; - - if ((end1 < 0) || (end1 >= invertices)) - { - if (Behavior.Verbose) - { - logger.Warning("Invalid first endpoint of segment.", "Mesh.FormSkeleton().1"); - } - } - else if ((end2 < 0) || (end2 >= invertices)) - { - if (Behavior.Verbose) - { - logger.Warning("Invalid second endpoint of segment.", "Mesh.FormSkeleton().2"); - } - } - else - { - // TODO: Is using the vertex ID reliable??? - // It should be. The ID gets appropriately set in TransferNodes(). - - // Find the vertices numbered 'end1' and 'end2'. - endpoint1 = vertices[end1]; - endpoint2 = vertices[end2]; - if ((endpoint1.x == endpoint2.x) && (endpoint1.y == endpoint2.y)) - { - if (Behavior.Verbose) - { - logger.Warning("Endpoints of segment (IDs " + end1 + "/" + end2 + ") are coincident.", - "Mesh.FormSkeleton()"); - } - } - else - { - InsertSegment(endpoint1, endpoint2, boundmarker); - } - } - } - } - - if (behavior.Convex || !behavior.Poly) - { - // Enclose the convex hull with subsegments. - MarkHull(); - } - } - - #endregion - #region Dealloc /// diff --git a/Triangle.NET/Triangle/MeshValidator.cs b/Triangle.NET/Triangle/MeshValidator.cs new file mode 100644 index 0000000..eda3ead --- /dev/null +++ b/Triangle.NET/Triangle/MeshValidator.cs @@ -0,0 +1,191 @@ +// ----------------------------------------------------------------------- +// +// Original Triangle code by Jonathan Richard Shewchuk, http://www.cs.cmu.edu/~quake/triangle.html +// Triangle.NET code by Christian Woltering, http://triangle.codeplex.com/ +// +// ----------------------------------------------------------------------- + +namespace TriangleNet +{ + using System; + using TriangleNet.Data; + using TriangleNet.Log; + + public static class MeshValidator + { + /// + /// Test the mesh for topological consistency. + /// + public static bool IsConsistent(Mesh mesh) + { + Otri tri = default(Otri); + Otri oppotri = default(Otri), oppooppotri = default(Otri); + Vertex triorg, tridest, triapex; + Vertex oppoorg, oppodest; + int horrors; + bool saveexact; + + var logger = SimpleLog.Instance; + + // Temporarily turn on exact arithmetic if it's off. + saveexact = Behavior.NoExact; + Behavior.NoExact = false; + horrors = 0; + + // Run through the list of triangles, checking each one. + foreach (var t in mesh.triangles.Values) + { + tri.triangle = t; + + // Check all three edges of the triangle. + for (tri.orient = 0; tri.orient < 3; tri.orient++) + { + triorg = tri.Org(); + tridest = tri.Dest(); + if (tri.orient == 0) + { // Only test for inversion once. + // Test if the triangle is flat or inverted. + triapex = tri.Apex(); + if (Primitives.CounterClockwise(triorg, tridest, triapex) <= 0.0) + { + logger.Warning("Triangle is flat or inverted.", + "Quality.CheckMesh()"); + horrors++; + } + } + // Find the neighboring triangle on this edge. + tri.Sym(ref oppotri); + if (oppotri.triangle != Mesh.dummytri) + { + // Check that the triangle's neighbor knows it's a neighbor. + oppotri.Sym(ref oppooppotri); + if ((tri.triangle != oppooppotri.triangle) || (tri.orient != oppooppotri.orient)) + { + if (tri.triangle == oppooppotri.triangle) + { + logger.Warning("Asymmetric triangle-triangle bond: (Right triangle, wrong orientation)", + "Quality.CheckMesh()"); + } + + horrors++; + } + // Check that both triangles agree on the identities + // of their shared vertices. + oppoorg = oppotri.Org(); + oppodest = oppotri.Dest(); + if ((triorg != oppodest) || (tridest != oppoorg)) + { + logger.Warning("Mismatched edge coordinates between two triangles.", + "Quality.CheckMesh()"); + + horrors++; + } + } + } + } + + // Check for unconnected vertices + mesh.MakeVertexMap(); + foreach (var v in mesh.vertices.Values) + { + if (v.tri.triangle == null) + { + logger.Warning("Vertex (ID " + v.id + ") not connected to mesh (duplicate input vertex?)", + "Quality.CheckMesh()"); + } + } + + if (horrors == 0) // && Behavior.Verbose + { + logger.Info("Mesh topology appears to be consistent."); + } + + // Restore the status of exact arithmetic. + Behavior.NoExact = saveexact; + + return (horrors == 0); + } + + /// + /// Ensure that the mesh is (constrained) Delaunay. + /// + public static bool IsDelaunay(Mesh mesh) + { + Otri loop = default(Otri); + Otri oppotri = default(Otri); + Osub opposubseg = default(Osub); + Vertex triorg, tridest, triapex; + Vertex oppoapex; + bool shouldbedelaunay; + int horrors; + bool saveexact; + + var logger = SimpleLog.Instance; + + // Temporarily turn on exact arithmetic if it's off. + saveexact = Behavior.NoExact; + Behavior.NoExact = false; + horrors = 0; + + // Run through the list of triangles, checking each one. + foreach (var tri in mesh.triangles.Values) + { + loop.triangle = tri; + + // Check all three edges of the triangle. + for (loop.orient = 0; loop.orient < 3; + loop.orient++) + { + 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 + // each pair isn't tested twice). + shouldbedelaunay = (oppotri.triangle != Mesh.dummytri) && + !Otri.IsDead(oppotri.triangle) && loop.triangle.id < oppotri.triangle.id && + (triorg != mesh.infvertex1) && (triorg != mesh.infvertex2) && + (triorg != mesh.infvertex3) && + (tridest != mesh.infvertex1) && (tridest != mesh.infvertex2) && + (tridest != mesh.infvertex3) && + (triapex != mesh.infvertex1) && (triapex != mesh.infvertex2) && + (triapex != mesh.infvertex3) && + (oppoapex != mesh.infvertex1) && (oppoapex != mesh.infvertex2) && + (oppoapex != mesh.infvertex3); + if (mesh.checksegments && shouldbedelaunay) + { + // If a subsegment separates the triangles, then the edge is + // constrained, so no local Delaunay test should be done. + loop.SegPivot(ref opposubseg); + if (opposubseg.seg != Mesh.dummysub) + { + shouldbedelaunay = false; + } + } + if (shouldbedelaunay) + { + if (Primitives.NonRegular(triorg, tridest, triapex, oppoapex) > 0.0) + { + logger.Warning(String.Format("Non-regular pair of triangles found (IDs {0}/{1}).", + loop.triangle.id, oppotri.triangle.id), "Quality.CheckDelaunay()"); + horrors++; + } + } + } + + } + + if (horrors == 0) // && Behavior.Verbose + { + logger.Info("Mesh is Delaunay."); + } + + // Restore the status of exact arithmetic. + Behavior.NoExact = saveexact; + + return (horrors == 0); + } + } +} diff --git a/Triangle.NET/Triangle/Quality.cs b/Triangle.NET/Triangle/Quality.cs index bc41819..522192d 100644 --- a/Triangle.NET/Triangle/Quality.cs +++ b/Triangle.NET/Triangle/Quality.cs @@ -25,9 +25,6 @@ namespace TriangleNet NewLocation newLocation; - // Not used at the moment - Func userTest; - ILog logger; public Quality(Mesh mesh) @@ -54,177 +51,6 @@ namespace TriangleNet #region Check - /// - /// Test the mesh for topological consistency. - /// - public bool CheckMesh() - { - Otri tri = default(Otri); - Otri oppotri = default(Otri), oppooppotri = default(Otri); - Vertex triorg, tridest, triapex; - Vertex oppoorg, oppodest; - int horrors; - bool saveexact; - - // Temporarily turn on exact arithmetic if it's off. - saveexact = Behavior.NoExact; - Behavior.NoExact = false; - horrors = 0; - - // Run through the list of triangles, checking each one. - foreach (var t in mesh.triangles.Values) - { - tri.triangle = t; - - // Check all three edges of the triangle. - for (tri.orient = 0; tri.orient < 3; tri.orient++) - { - triorg = tri.Org(); - tridest = tri.Dest(); - if (tri.orient == 0) - { // Only test for inversion once. - // Test if the triangle is flat or inverted. - triapex = tri.Apex(); - if (Primitives.CounterClockwise(triorg, tridest, triapex) <= 0.0) - { - logger.Warning("Triangle is flat or inverted.", - "Quality.CheckMesh()"); - horrors++; - } - } - // Find the neighboring triangle on this edge. - tri.Sym(ref oppotri); - if (oppotri.triangle != Mesh.dummytri) - { - // Check that the triangle's neighbor knows it's a neighbor. - oppotri.Sym(ref oppooppotri); - if ((tri.triangle != oppooppotri.triangle) || (tri.orient != oppooppotri.orient)) - { - if (tri.triangle == oppooppotri.triangle) - { - logger.Warning("Asymmetric triangle-triangle bond: (Right triangle, wrong orientation)", - "Quality.CheckMesh()"); - } - - horrors++; - } - // Check that both triangles agree on the identities - // of their shared vertices. - oppoorg = oppotri.Org(); - oppodest = oppotri.Dest(); - if ((triorg != oppodest) || (tridest != oppoorg)) - { - logger.Warning("Mismatched edge coordinates between two triangles.", - "Quality.CheckMesh()"); - - horrors++; - } - } - } - } - - // Check for unconnected vertices - mesh.MakeVertexMap(); - foreach (var v in mesh.vertices.Values) - { - if (v.tri.triangle == null) - { - logger.Warning("Vertex (ID " + v.id + ") not connected to mesh (duplicate input vertex?)", - "Quality.CheckMesh()"); - } - } - - if (horrors == 0) // && Behavior.Verbose - { - logger.Info("Mesh topology appears to be consistent."); - } - - // Restore the status of exact arithmetic. - Behavior.NoExact = saveexact; - - return (horrors == 0); - } - - /// - /// Ensure that the mesh is (constrained) Delaunay. - /// - public bool CheckDelaunay() - { - Otri loop = default(Otri); - Otri oppotri = default(Otri); - Osub opposubseg = default(Osub); - Vertex triorg, tridest, triapex; - Vertex oppoapex; - bool shouldbedelaunay; - int horrors; - bool saveexact; - - // Temporarily turn on exact arithmetic if it's off. - saveexact = Behavior.NoExact; - Behavior.NoExact = false; - horrors = 0; - - // Run through the list of triangles, checking each one. - foreach (var tri in mesh.triangles.Values) - { - loop.triangle = tri; - - // Check all three edges of the triangle. - for (loop.orient = 0; loop.orient < 3; - loop.orient++) - { - 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 - // each pair isn't tested twice). - shouldbedelaunay = (oppotri.triangle != Mesh.dummytri) && - !Otri.IsDead(oppotri.triangle) && loop.triangle.id < oppotri.triangle.id && - (triorg != mesh.infvertex1) && (triorg != mesh.infvertex2) && - (triorg != mesh.infvertex3) && - (tridest != mesh.infvertex1) && (tridest != mesh.infvertex2) && - (tridest != mesh.infvertex3) && - (triapex != mesh.infvertex1) && (triapex != mesh.infvertex2) && - (triapex != mesh.infvertex3) && - (oppoapex != mesh.infvertex1) && (oppoapex != mesh.infvertex2) && - (oppoapex != mesh.infvertex3); - if (mesh.checksegments && shouldbedelaunay) - { - // If a subsegment separates the triangles, then the edge is - // constrained, so no local Delaunay test should be done. - loop.SegPivot(ref opposubseg); - if (opposubseg.seg != Mesh.dummysub) - { - shouldbedelaunay = false; - } - } - if (shouldbedelaunay) - { - if (Primitives.NonRegular(triorg, tridest, triapex, oppoapex) > 0.0) - { - logger.Warning(String.Format("Non-regular pair of triangles found (IDs {0}/{1}).", - loop.triangle.id, oppotri.triangle.id), "Quality.CheckDelaunay()"); - horrors++; - } - } - } - - } - - if (horrors == 0) // && Behavior.Verbose - { - logger.Info("Mesh is Delaunay."); - } - - // Restore the status of exact arithmetic. - Behavior.NoExact = saveexact; - - return (horrors == 0); - } - /// /// Check a subsegment to see if it is encroached; add it to the list if it is. /// @@ -421,7 +247,7 @@ namespace TriangleNet testtri.Lprev(ref tri1); } - if (behavior.VarArea || behavior.fixedArea || behavior.Usertest) + if (behavior.VarArea || behavior.fixedArea || (behavior.UserTest != null)) { // Check whether the area is larger than permitted. area = 0.5 * (dxod * dyda - dyod * dxda); @@ -441,9 +267,9 @@ namespace TriangleNet } // Check whether the user thinks this triangle is too large. - if (behavior.Usertest && userTest != null) + if (behavior.UserTest != null) { - if (userTest(torg, tdest, tapex, area)) + if (behavior.UserTest(testtri.triangle, area)) { queue.Enqueue(ref testtri, minedge, tapex, torg, tdest); return; @@ -955,7 +781,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 != null) { // TODO: Reset queue? (Or is it always empty at this point) diff --git a/Triangle.NET/Triangle/Triangle.csproj b/Triangle.NET/Triangle/Triangle.csproj index c73b25b..488da34 100644 --- a/Triangle.NET/Triangle/Triangle.csproj +++ b/Triangle.NET/Triangle/Triangle.csproj @@ -45,6 +45,7 @@ + @@ -73,6 +74,7 @@ +