huanld/Triangle.NET
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Some code reorganization

Browse Source 
git-svn-id: https://triangle.svn.codeplex.com/svn@74818 0e2699bc-83d4-4a8f-98e7-55e24ab8c7a5
This commit is contained in:
SND\wo80_cp
2014-04-06 15:20:07 +00:00
parent e80534ac05
commit 616f86c52b
9 changed files with 1004 additions and 948 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Triangle.NET/TestApp/FormMain.cs
+6 -3
View File
@@ -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();
}
Triangle.NET/TestApp/Generators/RingPolygon.cs
+2 -2
View File
@@ -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);
Triangle.NET/TestApp/IO/GeometryWriter.cs
+16 -7
View File
@@ -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;
/// <summary>
@@ -17,8 +16,18 @@ namespace MeshExplorer.IO
/// </summary>
public static class GeometryWriter
{
public static void Write(InputGeometry geometry, string filename)
private static int OFFSET = 0;
/// <summary>
/// Writes an InputGeometry to a Triangle format .poly file.
/// </summary>
/// <param name="geometry">The InputGeometry to write.</param>
/// <param name="filename">The filename.</param>
/// <param name="compatibleMode">If true, indices will start at 1 (compatible with original C code).</param>
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<Point> 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<Edge> 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<Point> holes)
{
int index = 0;
int index = OFFSET;
writer.WriteLine("{0}", holes.Count());
Triangle.NET/Triangle/Behavior.cs
+4 -2
View File
@@ -8,6 +8,7 @@
namespace TriangleNet
{
using System;
using TriangleNet.Geometry;
using TriangleNet.Log;
/// <summary>
@@ -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<ITriangle, double, bool> usertest;
int noBisect = 0;
int steiner = -1;
@@ -181,7 +183,7 @@ namespace TriangleNet
/// <summary>
/// Apply a user-defined triangle constraint.
/// </summary>
public bool Usertest
public Func<ITriangle, double, bool> UserTest
{
get { return usertest; }
set { usertest = value; }
Triangle.NET/Triangle/ConstraintMesher.cs
+776
View File
@@ -0,0 +1,776 @@
// -----------------------------------------------------------------------
// <copyright file="ConstraintMesher.cs">
// 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/
// </copyright>
// -----------------------------------------------------------------------
namespace TriangleNet
{
using System;
using TriangleNet.Data;
using TriangleNet.Geometry;
using TriangleNet.Log;
internal class ConstraintMesher
{
Mesh mesh;
Behavior behavior;
TriangleLocator locator;
ILog<SimpleLogItem> logger;
public ConstraintMesher(Mesh mesh)
{
this.mesh = mesh;
this.behavior = mesh.behavior;
this.locator = mesh.locator;
logger = SimpleLog.Instance;
}
/// <summary>
/// Create the segments of a triangulation, including PSLG segments and edges
/// on the convex hull.
/// </summary>
/// <param name="segmentlist"></param>
/// <param name="segmentmarkerlist"></param>
/// <param name="numberofsegments"></param>
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
/// <summary>
/// Find the first triangle on the path from one point to another.
/// </summary>
/// <param name="searchtri"></param>
/// <param name="searchpoint"></param>
/// <returns>
/// The return value notes whether the destination or apex of the found
/// triangle is collinear with the two points in question.</returns>
/// <remarks>
/// 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.
/// </remarks>
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;
}
}
/// <summary>
/// Find the intersection of an existing segment and a segment that is being
/// inserted. Insert a vertex at the intersection, splitting an existing subsegment.
/// </summary>
/// <param name="splittri"></param>
/// <param name="splitsubseg"></param>
/// <param name="endpoint2"></param>
/// <remarks>
/// 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.
/// </remarks>
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.
}
/// <summary>
/// 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.
/// </summary>
/// <param name="searchtri"></param>
/// <param name="endpoint2"></param>
/// <param name="newmark"></param>
/// <returns>Returns true if the entire segment is successfully inserted, and false
/// if the job must be finished by ConstrainedEdge().</returns>
/// <remarks>
/// 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.
/// </remarks>
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);
}
}
}
/// <summary>
/// Enforce the Delaunay condition at an edge, fanning out recursively from
/// an existing vertex. Pay special attention to stacking inverted triangles.
/// </summary>
/// <param name="fixuptri"></param>
/// <param name="leftside">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.)</param>
/// <remarks>
/// 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.
/// </remarks>
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);
}
/// <summary>
/// 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.
/// </summary>
/// <param name="starttri"></param>
/// <param name="endpoint2"></param>
/// <param name="newmark"></param>
/// <remarks>
/// 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.
/// </remarks>
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);
}
}
}
/// <summary>
/// Insert a PSLG segment into a triangulation.
/// </summary>
/// <param name="endpoint1"></param>
/// <param name="endpoint2"></param>
/// <param name="newmark"></param>
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);
}
/// <summary>
/// Cover the convex hull of a triangulation with subsegments.
/// </summary>
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
}
}
Triangle.NET/Triangle/Mesh.cs
+3 -756
View File
@@ -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
}
}
/// <summary>
/// Check mesh consistency and (constrained) Delaunay property.
/// </summary>
/// <param name="isConsistent">Value indicating if mesh topology is consistent.</param>
/// <param name="isDelaunay">Value indicating if mesh is Delaunay.</param>
public void Check(out bool isConsistent, out bool isDelaunay)
{
isConsistent = quality.CheckMesh();
isDelaunay = quality.CheckDelaunay();
}
#region Misc
/// <summary>
@@ -1982,750 +1973,6 @@ namespace TriangleNet
#endregion
#region Segment insertion
/// <summary>
/// Find the first triangle on the path from one point to another.
/// </summary>
/// <param name="searchtri"></param>
/// <param name="searchpoint"></param>
/// <returns>
/// The return value notes whether the destination or apex of the found
/// triangle is collinear with the two points in question.</returns>
/// <remarks>
/// 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.
/// </remarks>
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;
}
}
/// <summary>
/// Find the intersection of an existing segment and a segment that is being
/// inserted. Insert a vertex at the intersection, splitting an existing subsegment.
/// </summary>
/// <param name="splittri"></param>
/// <param name="splitsubseg"></param>
/// <param name="endpoint2"></param>
/// <remarks>
/// 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.
/// </remarks>
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.
}
/// <summary>
/// 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.
/// </summary>
/// <param name="searchtri"></param>
/// <param name="endpoint2"></param>
/// <param name="newmark"></param>
/// <returns>Returns true if the entire segment is successfully inserted, and false
/// if the job must be finished by ConstrainedEdge().</returns>
/// <remarks>
/// 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.
/// </remarks>
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);
}
}
}
/// <summary>
/// Enforce the Delaunay condition at an edge, fanning out recursively from
/// an existing vertex. Pay special attention to stacking inverted triangles.
/// </summary>
/// <param name="fixuptri"></param>
/// <param name="leftside">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.)</param>
/// <remarks>
/// 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.
/// </remarks>
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);
}
/// <summary>
/// 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.
/// </summary>
/// <param name="starttri"></param>
/// <param name="endpoint2"></param>
/// <param name="newmark"></param>
/// <remarks>
/// 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.
/// </remarks>
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);
}
}
}
/// <summary>
/// Insert a PSLG segment into a triangulation.
/// </summary>
/// <param name="endpoint1"></param>
/// <param name="endpoint2"></param>
/// <param name="newmark"></param>
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);
}
/// <summary>
/// Cover the convex hull of a triangulation with subsegments.
/// </summary>
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));
}
/// <summary>
/// Create the segments of a triangulation, including PSLG segments and edges
/// on the convex hull.
/// </summary>
/// <param name="segmentlist"></param>
/// <param name="segmentmarkerlist"></param>
/// <param name="numberofsegments"></param>
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
/// <summary>
Triangle.NET/Triangle/MeshValidator.cs
+191
View File
@@ -0,0 +1,191 @@
// -----------------------------------------------------------------------
// <copyright file="MeshValidator.cs">
// 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/
// </copyright>
// -----------------------------------------------------------------------
namespace TriangleNet
{
using System;
using TriangleNet.Data;
using TriangleNet.Log;
public static class MeshValidator
{
/// <summary>
/// Test the mesh for topological consistency.
/// </summary>
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);
}
/// <summary>
/// Ensure that the mesh is (constrained) Delaunay.
/// </summary>
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);
}
}
}
Triangle.NET/Triangle/Quality.cs
+4 -178
View File
@@ -25,9 +25,6 @@ namespace TriangleNet
NewLocation newLocation;
// Not used at the moment
Func<Point, Point, Point, double, bool> userTest;
ILog<SimpleLogItem> logger;
public Quality(Mesh mesh)
@@ -54,177 +51,6 @@ namespace TriangleNet
#region Check
/// <summary>
/// Test the mesh for topological consistency.
/// </summary>
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);
}
/// <summary>
/// Ensure that the mesh is (constrained) Delaunay.
/// </summary>
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);
}
/// <summary>
/// Check a subsegment to see if it is encroached; add it to the list if it is.
/// </summary>
@@ -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)
Triangle.NET/Triangle/Triangle.csproj
+2
View File
@@ -45,6 +45,7 @@
<Compile Include="BadTriQueue.cs" />
<Compile Include="Behavior.cs" />
<Compile Include="Carver.cs" />
<Compile Include="ConstraintMesher.cs" />
<Compile Include="Data\BadSubseg.cs" />
<Compile Include="Data\BadTriangle.cs" />
<Compile Include="Data\Osub.cs" />
@@ -73,6 +74,7 @@
<Compile Include="Log\ILogItem.cs" />
<Compile Include="Log\SimpleLog.cs" />
<Compile Include="Log\SimpleLogItem.cs" />
<Compile Include="MeshValidator.cs" />
<Compile Include="NewLocation.cs" />
<Compile Include="Quality.cs" />
<Compile Include="Enums.cs" />