diff --git a/Triangle.NET/Triangle/IO/TriangleWriter.cs b/Triangle.NET/Triangle/IO/TriangleWriter.cs
index 19a95af..99794a2 100644
--- a/Triangle.NET/Triangle/IO/TriangleWriter.cs
+++ b/Triangle.NET/Triangle/IO/TriangleWriter.cs
@@ -453,167 +453,5 @@ namespace TriangleNet.IO
}
}
}
-
- ///
- /// Write the Voronoi diagram to a .voro file.
- ///
- ///
- ///
- ///
- ///
- /// The Voronoi diagram is the geometric dual of the Delaunay triangulation.
- /// Hence, the Voronoi vertices are listed by traversing the Delaunay
- /// triangles, and the Voronoi edges are listed by traversing the Delaunay
- /// edges.
- ///
- /// WARNING: In order to assign numbers to the Voronoi vertices, this
- /// procedure messes up the subsegments or the extra nodes of every
- /// element. Hence, you should call this procedure last.
- public static void WriteVoronoi(Mesh mesh, string filename)
- {
- Otri tri = default(Otri), trisym = default(Otri);
- Vertex torg, tdest, tapex;
- Point circumcenter;
- double xi = 0, eta = 0;
-
- int p1, p2, index = 0;
- tri.orient = 0;
-
- using (StreamWriter writer = new StreamWriter(filename))
- {
- // Number of triangles, two dimensions, number of vertex attributes, no markers.
- writer.WriteLine("{0} 2 {1} 0", mesh.triangles.Count, mesh.nextras);
-
- foreach (var item in mesh.triangles.Values)
- {
- tri.tri = item;
- torg = tri.Org();
- tdest = tri.Dest();
- tapex = tri.Apex();
- circumcenter = RobustPredicates.FindCircumcenter(torg, tdest, tapex, ref xi, ref eta);
-
- // X and y coordinates.
- writer.Write("{0} {1} {2}", index, circumcenter.X.ToString(nfi),
- circumcenter.Y.ToString(nfi));
-
- for (int i = 0; i < mesh.nextras; i++)
- {
- writer.Write(" 0");
- // TODO
- // Interpolate the vertex attributes at the circumcenter.
- //writer.Write(" {0}", torg.attribs[i] + xi * (tdes.attribst[i] - torg.attribs[i]) +
- // eta * (tapex.attribs[i] - torg.attribs[i]));
- }
- writer.WriteLine();
-
- tri.tri.id = index++;
- }
-
-
- // Number of edges, zero boundary markers.
- writer.WriteLine("{0} 0", mesh.edges);
-
- index = 0;
- // To loop over the set of edges, loop over all triangles, and look at
- // the three edges of each triangle. If there isn't another triangle
- // adjacent to the edge, operate on the edge. If there is another
- // adjacent triangle, operate on the edge only if the current triangle
- // has a smaller pointer than its neighbor. This way, each edge is
- // considered only once.
- foreach (var item in mesh.triangles.Values)
- {
- tri.tri = item;
-
- for (tri.orient = 0; tri.orient < 3; tri.orient++)
- {
- tri.Sym(ref trisym);
- if ((tri.tri.id < trisym.tri.id) || (trisym.tri.id == Mesh.DUMMY))
- {
- // Find the number of this triangle (and Voronoi vertex).
- p1 = tri.tri.id;
-
- if (trisym.tri.id == Mesh.DUMMY)
- {
- torg = tri.Org();
- tdest = tri.Dest();
-
- // Write an infinite ray. Edge number, index of one endpoint,
- // -1, and x and y coordinates of a vector representing the
- // direction of the ray.
- writer.WriteLine("{0} {1} -1 {2} {3}", index, p1,
- (tdest[1] - torg[1]).ToString(nfi),
- (torg[0] - tdest[0]).ToString(nfi));
- }
- else
- {
- // Find the number of the adjacent triangle (and Voronoi vertex).
- p2 = trisym.tri.id;
- // Finite edge. Write indices of two endpoints.
- writer.WriteLine("{0} {1} {2}", index, p1, p2);
- }
-
- index++;
- }
- }
- }
- }
- }
-
- ///
- /// Write the triangulation to an .off file.
- ///
- ///
- ///
- ///
- /// OFF stands for the Object File Format, a format used by the Geometry
- /// Center's Geomview package.
- ///
- public static void WriteOffFile(Mesh mesh, string filename)
- {
- Otri tri;
- Vertex p1, p2, p3;
-
- long outvertices = mesh.vertices.Count;
-
- if (mesh.behavior.Jettison)
- {
- outvertices = mesh.vertices.Count - mesh.undeads;
- }
-
- int index = 0;
-
- using (StreamWriter writer = new StreamWriter(filename))
- {
- writer.WriteLine("OFF");
- writer.WriteLine("{0} {1} {2}", outvertices, mesh.triangles.Count, mesh.edges);
-
- foreach (var item in mesh.vertices.Values)
- {
- p1 = item;
-
- if (!mesh.behavior.Jettison || p1.type != VertexType.UndeadVertex)
- {
- // The "0.0" is here because the OFF format uses 3D coordinates.
- writer.WriteLine(" {0} {1} 0.0", p1[0].ToString(nfi), p1[1].ToString(nfi));
-
- p1.id = index++;
- }
- }
-
- // Write the triangles.
- tri.orient = 0;
- foreach (var item in mesh.triangles.Values)
- {
- tri.tri = item;
-
- p1 = tri.Org();
- p2 = tri.Dest();
- p3 = tri.Apex();
-
- // The "3" means a three-vertex polygon.
- writer.WriteLine(" 3 {0} {1} {2}", p1.id, p2.id, p3.id);
- }
- }
- }
}
}
diff --git a/Triangle.NET/Triangle/IPredicates.cs b/Triangle.NET/Triangle/IPredicates.cs
new file mode 100644
index 0000000..895690d
--- /dev/null
+++ b/Triangle.NET/Triangle/IPredicates.cs
@@ -0,0 +1,22 @@
+// -----------------------------------------------------------------------
+//
+// Triangle.NET code by Christian Woltering, http://triangle.codeplex.com/
+//
+// -----------------------------------------------------------------------
+
+namespace TriangleNet
+{
+ using TriangleNet.Geometry;
+
+ public interface IPredicates
+ {
+ double CounterClockwise(Point a, Point b, Point c);
+
+ double InCircle(Point a, Point b, Point c, Point p);
+
+ Point FindCircumcenter(Point org, Point dest, Point apex, ref double xi, ref double eta);
+
+ Point FindCircumcenter(Point org, Point dest, Point apex, ref double xi, ref double eta,
+ double offconstant);
+ }
+}
diff --git a/Triangle.NET/Triangle/Mesh.cs b/Triangle.NET/Triangle/Mesh.cs
index a3028df..0d15f61 100644
--- a/Triangle.NET/Triangle/Mesh.cs
+++ b/Triangle.NET/Triangle/Mesh.cs
@@ -24,6 +24,8 @@ namespace TriangleNet
{
#region Variables
+ IPredicates predicates;
+
ILog logger;
QualityMesher qualityMesher;
@@ -221,7 +223,7 @@ namespace TriangleNet
///
/// Initializes a new instance of the class.
///
- public Mesh()
+ public Mesh(IPredicates predicates)
{
Initialize();
@@ -238,11 +240,11 @@ namespace TriangleNet
holes = new List();
regions = new List();
- qualityMesher = new QualityMesher(this);
+ this.predicates = predicates;
- locator = new TriangleLocator(this);
+ this.qualityMesher = new QualityMesher(this, predicates);
- RobustPredicates.ExactInit();
+ this.locator = new TriangleLocator(this, predicates);
}
public void Refine(QualityOptions quality)
@@ -439,7 +441,7 @@ namespace TriangleNet
infvertex3 = null;
// Insert segments, carving holes.
- var mesher = new ConstraintMesher(this);
+ var mesher = new ConstraintMesher(this, predicates);
if (behavior.useSegments)
{
@@ -1125,7 +1127,7 @@ namespace TriangleNet
// the boundary of the triangulation. 'farvertex' might be
// infinite as well, but trust me, this same condition should
// be applied.
- doflip = RobustPredicates.CounterClockwise(newvertex, rightvertex, farvertex) > 0.0;
+ doflip = predicates.CounterClockwise(newvertex, rightvertex, farvertex) > 0.0;
}
else if ((rightvertex == infvertex1) ||
(rightvertex == infvertex2) ||
@@ -1135,7 +1137,7 @@ namespace TriangleNet
// the boundary of the triangulation. 'farvertex' might be
// infinite as well, but trust me, this same condition should
// be applied.
- doflip = RobustPredicates.CounterClockwise(farvertex, leftvertex, newvertex) > 0.0;
+ doflip = predicates.CounterClockwise(farvertex, leftvertex, newvertex) > 0.0;
}
else if ((farvertex == infvertex1) ||
(farvertex == infvertex2) ||
@@ -1148,7 +1150,7 @@ namespace TriangleNet
else
{
// Test whether the edge is locally Delaunay.
- doflip = RobustPredicates.InCircle(leftvertex, newvertex, rightvertex, farvertex) > 0.0;
+ doflip = predicates.InCircle(leftvertex, newvertex, rightvertex, farvertex) > 0.0;
}
if (doflip)
{
@@ -1675,7 +1677,7 @@ namespace TriangleNet
testtri.Onext();
testvertex = testtri.Dest();
// Is this a better vertex?
- if (RobustPredicates.InCircle(leftbasevertex, rightbasevertex, bestvertex, testvertex) > 0.0)
+ if (predicates.InCircle(leftbasevertex, rightbasevertex, bestvertex, testvertex) > 0.0)
{
testtri.Copy(ref besttri);
bestvertex = testvertex;
diff --git a/Triangle.NET/Triangle/MeshValidator.cs b/Triangle.NET/Triangle/MeshValidator.cs
index 59c8dd8..64bad8b 100644
--- a/Triangle.NET/Triangle/MeshValidator.cs
+++ b/Triangle.NET/Triangle/MeshValidator.cs
@@ -13,6 +13,8 @@ namespace TriangleNet
public static class MeshValidator
{
+ private static RobustPredicates predicates = RobustPredicates.Default;
+
///
/// Test the mesh for topological consistency.
///
@@ -46,7 +48,7 @@ namespace TriangleNet
// Only test for inversion once.
// Test if the triangle is flat or inverted.
apex = tri.Apex();
- if (RobustPredicates.CounterClockwise(org, dest, apex) <= 0.0)
+ if (predicates.CounterClockwise(org, dest, apex) <= 0.0)
{
if (Log.Verbose)
{
@@ -189,7 +191,7 @@ namespace TriangleNet
if (shouldbedelaunay)
{
- if (RobustPredicates.NonRegular(org, dest, apex, oppoapex) > 0.0)
+ if (predicates.NonRegular(org, dest, apex, oppoapex) > 0.0)
{
if (Log.Verbose)
{
diff --git a/Triangle.NET/Triangle/Meshing/Algorithm/Dwyer.cs b/Triangle.NET/Triangle/Meshing/Algorithm/Dwyer.cs
index 35cb035..f8b9950 100644
--- a/Triangle.NET/Triangle/Meshing/Algorithm/Dwyer.cs
+++ b/Triangle.NET/Triangle/Meshing/Algorithm/Dwyer.cs
@@ -45,13 +45,26 @@ namespace TriangleNet.Meshing.Algorithm
///
public class Dwyer : ITriangulator
{
- static Random rand = new Random(DateTime.Now.Millisecond);
+ // Random is not threadsafe, so don't make this static.
+ Random rand = new Random(DateTime.Now.Millisecond);
+
+ IPredicates predicates;
public bool UseDwyer = true;
Vertex[] sortarray;
Mesh mesh;
+ public Dwyer()
+ : this(RobustPredicates.Default)
+ {
+ }
+
+ public Dwyer(IPredicates predicates)
+ {
+ this.predicates = predicates;
+ }
+
///
/// Form a Delaunay triangulation by the divide-and-conquer method.
///
@@ -62,7 +75,7 @@ namespace TriangleNet.Meshing.Algorithm
///
public IMesh Triangulate(IList points)
{
- this.mesh = new Mesh();
+ this.mesh = new Mesh(predicates);
this.mesh.TransferNodes(points);
Otri hullleft = default(Otri), hullright = default(Otri);
@@ -438,7 +451,7 @@ namespace TriangleNet.Meshing.Algorithm
{
changemade = false;
// Make innerleftdest the "bottommost" vertex of the left hull.
- if (RobustPredicates.CounterClockwise(innerleftdest, innerleftapex, innerrightorg) > 0.0)
+ if (predicates.CounterClockwise(innerleftdest, innerleftapex, innerrightorg) > 0.0)
{
innerleft.Lprev();
innerleft.Sym();
@@ -447,7 +460,7 @@ namespace TriangleNet.Meshing.Algorithm
changemade = true;
}
// Make innerrightorg the "bottommost" vertex of the right hull.
- if (RobustPredicates.CounterClockwise(innerrightapex, innerrightorg, innerleftdest) > 0.0)
+ if (predicates.CounterClockwise(innerrightapex, innerrightorg, innerleftdest) > 0.0)
{
innerright.Lnext();
innerright.Sym();
@@ -495,8 +508,8 @@ namespace TriangleNet.Meshing.Algorithm
// because even though the left triangulation might seem finished now,
// moving up on the right triangulation might reveal a new vertex of
// the left triangulation. And vice-versa.)
- leftfinished = RobustPredicates.CounterClockwise(upperleft, lowerleft, lowerright) <= 0.0;
- rightfinished = RobustPredicates.CounterClockwise(upperright, lowerleft, lowerright) <= 0.0;
+ leftfinished = predicates.CounterClockwise(upperleft, lowerleft, lowerright) <= 0.0;
+ rightfinished = predicates.CounterClockwise(upperright, lowerleft, lowerright) <= 0.0;
if (leftfinished && rightfinished)
{
// Create the top new bounding triangle.
@@ -553,7 +566,7 @@ namespace TriangleNet.Meshing.Algorithm
if (nextapex != null)
{
// Check whether the edge is Delaunay.
- badedge = RobustPredicates.InCircle(lowerleft, lowerright, upperleft, nextapex) > 0.0;
+ badedge = predicates.InCircle(lowerleft, lowerright, upperleft, nextapex) > 0.0;
while (badedge)
{
// Eliminate the edge with an edge flip. As a result, the
@@ -583,7 +596,7 @@ namespace TriangleNet.Meshing.Algorithm
if (nextapex != null)
{
// Check whether the edge is Delaunay.
- badedge = RobustPredicates.InCircle(lowerleft, lowerright, upperleft, nextapex) > 0.0;
+ badedge = predicates.InCircle(lowerleft, lowerright, upperleft, nextapex) > 0.0;
}
else
{
@@ -605,7 +618,7 @@ namespace TriangleNet.Meshing.Algorithm
if (nextapex != null)
{
// Check whether the edge is Delaunay.
- badedge = RobustPredicates.InCircle(lowerleft, lowerright, upperright, nextapex) > 0.0;
+ badedge = predicates.InCircle(lowerleft, lowerright, upperright, nextapex) > 0.0;
while (badedge)
{
// Eliminate the edge with an edge flip. As a result, the
@@ -635,7 +648,7 @@ namespace TriangleNet.Meshing.Algorithm
if (nextapex != null)
{
// Check whether the edge is Delaunay.
- badedge = RobustPredicates.InCircle(lowerleft, lowerright, upperright, nextapex) > 0.0;
+ badedge = predicates.InCircle(lowerleft, lowerright, upperright, nextapex) > 0.0;
}
else
{
@@ -646,7 +659,7 @@ namespace TriangleNet.Meshing.Algorithm
}
}
if (leftfinished || (!rightfinished &&
- (RobustPredicates.InCircle(upperleft, lowerleft, lowerright, upperright) > 0.0)))
+ (predicates.InCircle(upperleft, lowerleft, lowerright, upperright) > 0.0)))
{
// Knit the triangulations, adding an edge from 'lowerleft'
// to 'upperright'.
@@ -735,7 +748,7 @@ namespace TriangleNet.Meshing.Algorithm
mesh.MakeTriangle(ref tri1);
mesh.MakeTriangle(ref tri2);
mesh.MakeTriangle(ref tri3);
- area = RobustPredicates.CounterClockwise(sortarray[left], sortarray[left + 1], sortarray[left + 2]);
+ area = predicates.CounterClockwise(sortarray[left], sortarray[left + 1], sortarray[left + 2]);
if (area == 0.0)
{
// Three collinear vertices; the triangulation is two edges.
diff --git a/Triangle.NET/Triangle/Meshing/Algorithm/Incremental.cs b/Triangle.NET/Triangle/Meshing/Algorithm/Incremental.cs
index abf73ab..ae420ea 100644
--- a/Triangle.NET/Triangle/Meshing/Algorithm/Incremental.cs
+++ b/Triangle.NET/Triangle/Meshing/Algorithm/Incremental.cs
@@ -16,8 +16,20 @@ namespace TriangleNet.Meshing.Algorithm
///
public class Incremental : ITriangulator
{
+ IPredicates predicates;
+
Mesh mesh;
+ public Incremental()
+ : this(RobustPredicates.Default)
+ {
+ }
+
+ public Incremental(IPredicates predicates)
+ {
+ this.predicates = predicates;
+ }
+
///
/// Form a Delaunay triangulation by incrementally inserting vertices.
///
@@ -25,7 +37,7 @@ namespace TriangleNet.Meshing.Algorithm
/// triangulation.
public IMesh Triangulate(IList points)
{
- this.mesh = new Mesh();
+ this.mesh = new Mesh(predicates);
this.mesh.TransferNodes(points);
Otri starttri = new Otri();
diff --git a/Triangle.NET/Triangle/Meshing/Algorithm/SweepLine.cs b/Triangle.NET/Triangle/Meshing/Algorithm/SweepLine.cs
index b53fd19..d5dff4e 100644
--- a/Triangle.NET/Triangle/Meshing/Algorithm/SweepLine.cs
+++ b/Triangle.NET/Triangle/Meshing/Algorithm/SweepLine.cs
@@ -27,13 +27,25 @@ namespace TriangleNet.Meshing.Algorithm
return randomseed / (714025 / choices + 1);
}
+ IPredicates predicates;
+
Mesh mesh;
double xminextreme; // Nonexistent x value used as a flag in sweepline.
List splaynodes;
+ public SweepLine()
+ : this(RobustPredicates.Default)
+ {
+ }
+
+ public SweepLine(IPredicates predicates)
+ {
+ this.predicates = predicates;
+ }
+
public IMesh Triangulate(IList points)
{
- this.mesh = new Mesh();
+ this.mesh = new Mesh(predicates);
this.mesh.TransferNodes(points);
// Nonexistent x value used as a flag to mark circle events in sweepline
@@ -213,7 +225,7 @@ namespace TriangleNet.Meshing.Algorithm
leftvertex = farlefttri.Apex();
midvertex = lefttri.Dest();
rightvertex = lefttri.Apex();
- lefttest = RobustPredicates.CounterClockwise(leftvertex, midvertex, rightvertex);
+ lefttest = predicates.CounterClockwise(leftvertex, midvertex, rightvertex);
if (lefttest > 0.0)
{
newevent = new SweepEvent();
@@ -228,7 +240,7 @@ namespace TriangleNet.Meshing.Algorithm
leftvertex = righttri.Apex();
midvertex = righttri.Org();
rightvertex = farrighttri.Apex();
- righttest = RobustPredicates.CounterClockwise(leftvertex, midvertex, rightvertex);
+ righttest = predicates.CounterClockwise(leftvertex, midvertex, rightvertex);
if (righttest > 0.0)
{
newevent = new SweepEvent();
@@ -600,7 +612,7 @@ namespace TriangleNet.Meshing.Algorithm
Point searchpoint = new Point(); // TODO: mesh.nextras
Otri dummytri = default(Otri);
- ccwabc = RobustPredicates.CounterClockwise(pa, pb, pc);
+ ccwabc = predicates.CounterClockwise(pa, pb, pc);
xac = pa.x - pc.x;
yac = pa.y - pc.y;
xbc = pb.x - pc.x;
diff --git a/Triangle.NET/Triangle/Meshing/ConstraintMesher.cs b/Triangle.NET/Triangle/Meshing/ConstraintMesher.cs
index 7fd3756..7e83a53 100644
--- a/Triangle.NET/Triangle/Meshing/ConstraintMesher.cs
+++ b/Triangle.NET/Triangle/Meshing/ConstraintMesher.cs
@@ -16,6 +16,8 @@ namespace TriangleNet.Meshing
internal class ConstraintMesher
{
+ IPredicates predicates;
+
Mesh mesh;
Behavior behavior;
TriangleLocator locator;
@@ -24,9 +26,11 @@ namespace TriangleNet.Meshing
ILog logger;
- public ConstraintMesher(Mesh mesh)
+ public ConstraintMesher(Mesh mesh, IPredicates predicates)
{
this.mesh = mesh;
+ this.predicates = predicates;
+
this.behavior = mesh.behavior;
this.locator = mesh.locator;
@@ -74,7 +78,7 @@ namespace TriangleNet.Meshing
// falls within the starting triangle.
searchorg = searchtri.Org();
searchdest = searchtri.Dest();
- if (RobustPredicates.CounterClockwise(searchorg, searchdest, hole) > 0.0)
+ if (predicates.CounterClockwise(searchorg, searchdest, hole) > 0.0)
{
// Find a triangle that contains the hole.
intersect = mesh.locator.Locate(hole, ref searchtri);
@@ -116,7 +120,7 @@ namespace TriangleNet.Meshing
// region point falls within the starting triangle.
searchorg = searchtri.Org();
searchdest = searchtri.Dest();
- if (RobustPredicates.CounterClockwise(searchorg, searchdest, region.point) > 0.0)
+ if (predicates.CounterClockwise(searchorg, searchdest, region.point) > 0.0)
{
// Find a triangle that contains the region point.
intersect = mesh.locator.Locate(region.point, ref searchtri);
@@ -534,10 +538,10 @@ namespace TriangleNet.Meshing
rightvertex = searchtri.Dest();
leftvertex = searchtri.Apex();
// Is 'searchpoint' to the left?
- leftccw = RobustPredicates.CounterClockwise(searchpoint, startvertex, leftvertex);
+ leftccw = predicates.CounterClockwise(searchpoint, startvertex, leftvertex);
leftflag = leftccw > 0.0;
// Is 'searchpoint' to the right?
- rightccw = RobustPredicates.CounterClockwise(startvertex, searchpoint, rightvertex);
+ rightccw = predicates.CounterClockwise(startvertex, searchpoint, rightvertex);
rightflag = rightccw > 0.0;
if (leftflag && rightflag)
{
@@ -564,7 +568,7 @@ namespace TriangleNet.Meshing
}
leftvertex = searchtri.Apex();
rightccw = leftccw;
- leftccw = RobustPredicates.CounterClockwise(searchpoint, startvertex, leftvertex);
+ leftccw = predicates.CounterClockwise(searchpoint, startvertex, leftvertex);
leftflag = leftccw > 0.0;
}
while (rightflag)
@@ -578,7 +582,7 @@ namespace TriangleNet.Meshing
}
rightvertex = searchtri.Dest();
leftccw = rightccw;
- rightccw = RobustPredicates.CounterClockwise(startvertex, searchpoint, rightvertex);
+ rightccw = predicates.CounterClockwise(startvertex, searchpoint, rightvertex);
rightflag = rightccw > 0.0;
}
if (leftccw == 0.0)
@@ -859,7 +863,7 @@ namespace TriangleNet.Meshing
// Check whether the previous polygon vertex is a reflex vertex.
if (leftside)
{
- if (RobustPredicates.CounterClockwise(nearvertex, leftvertex, farvertex) <= 0.0)
+ if (predicates.CounterClockwise(nearvertex, leftvertex, farvertex) <= 0.0)
{
// leftvertex is a reflex vertex too. Nothing can
// be done until a convex section is found.
@@ -868,20 +872,20 @@ namespace TriangleNet.Meshing
}
else
{
- if (RobustPredicates.CounterClockwise(farvertex, rightvertex, nearvertex) <= 0.0)
+ if (predicates.CounterClockwise(farvertex, rightvertex, nearvertex) <= 0.0)
{
// rightvertex is a reflex vertex too. Nothing can
// be done until a convex section is found.
return;
}
}
- if (RobustPredicates.CounterClockwise(rightvertex, leftvertex, farvertex) > 0.0)
+ if (predicates.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 (RobustPredicates.InCircle(leftvertex, farvertex, rightvertex, nearvertex) <= 0.0)
+ if (predicates.InCircle(leftvertex, farvertex, rightvertex, nearvertex) <= 0.0)
{
return;
}
@@ -983,7 +987,7 @@ namespace TriangleNet.Meshing
{
// 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 = RobustPredicates.CounterClockwise(endpoint1, endpoint2, farvertex);
+ area = predicates.CounterClockwise(endpoint1, endpoint2, farvertex);
if (area == 0.0)
{
// We've collided with a vertex between endpoint1 and endpoint2.
diff --git a/Triangle.NET/Triangle/Meshing/Converter.cs b/Triangle.NET/Triangle/Meshing/Converter.cs
index be6d9b9..d839983 100644
--- a/Triangle.NET/Triangle/Meshing/Converter.cs
+++ b/Triangle.NET/Triangle/Meshing/Converter.cs
@@ -44,7 +44,7 @@ namespace TriangleNet.Meshing
int elements = triangles == null ? 0 : triangles.Length;
int segments = polygon.Segments.Count;
- var mesh = new Mesh();
+ var mesh = new Mesh(RobustPredicates.Default);
mesh.TransferNodes(polygon.Points);
diff --git a/Triangle.NET/Triangle/Meshing/GenericMesher.cs b/Triangle.NET/Triangle/Meshing/GenericMesher.cs
index 4aeaf50..618f9cd 100644
--- a/Triangle.NET/Triangle/Meshing/GenericMesher.cs
+++ b/Triangle.NET/Triangle/Meshing/GenericMesher.cs
@@ -1,4 +1,9 @@
-
+// -----------------------------------------------------------------------
+//
+// Triangle.NET code by Christian Woltering, http://triangle.codeplex.com/
+//
+// -----------------------------------------------------------------------
+
namespace TriangleNet.Meshing
{
using System;
@@ -12,15 +17,24 @@ namespace TriangleNet.Meshing
///
public class GenericMesher : ITriangulator, IConstraintMesher, IQualityMesher
{
+ IPredicates predicates;
+
ITriangulator triangulator;
public GenericMesher()
- : this(new Dwyer())
{
+ this.predicates = RobustPredicates.Default;
+ this.triangulator = new Dwyer(this.predicates);
}
public GenericMesher(ITriangulator triangulator)
+ : this(triangulator, RobustPredicates.Default)
{
+ }
+
+ public GenericMesher(ITriangulator triangulator, IPredicates predicates)
+ {
+ this.predicates = predicates;
this.triangulator = triangulator;
}
diff --git a/Triangle.NET/Triangle/Meshing/QualityMesher.cs b/Triangle.NET/Triangle/Meshing/QualityMesher.cs
index 00625eb..cf6f29f 100644
--- a/Triangle.NET/Triangle/Meshing/QualityMesher.cs
+++ b/Triangle.NET/Triangle/Meshing/QualityMesher.cs
@@ -19,6 +19,8 @@ namespace TriangleNet.Meshing
///
class QualityMesher
{
+ IPredicates predicates;
+
Queue badsubsegs;
BadTriQueue queue;
Mesh mesh;
@@ -28,7 +30,7 @@ namespace TriangleNet.Meshing
ILog logger;
- public QualityMesher(Mesh mesh)
+ public QualityMesher(Mesh mesh, IPredicates predicates)
{
logger = Log.Instance;
@@ -36,9 +38,11 @@ namespace TriangleNet.Meshing
queue = new BadTriQueue();
this.mesh = mesh;
+ this.predicates = predicates;
+
this.behavior = mesh.behavior;
- newLocation = new NewLocation(mesh);
+ newLocation = new NewLocation(mesh, predicates);
}
///
@@ -570,7 +574,7 @@ namespace TriangleNet.Meshing
// Roundoff in the above calculation may yield a 'newvertex'
// that is not precisely collinear with 'eorg' and 'edest'.
// Improve collinearity by one step of iterative refinement.
- multiplier = RobustPredicates.CounterClockwise(eorg, edest, newvertex);
+ multiplier = predicates.CounterClockwise(eorg, edest, newvertex);
divisor = ((eorg.x - edest.x) * (eorg.x - edest.x) +
(eorg.y - edest.y) * (eorg.y - edest.y));
if ((multiplier != 0.0) && (divisor != 0.0))
@@ -671,7 +675,7 @@ namespace TriangleNet.Meshing
// reset VertexType?
if (behavior.fixedArea || behavior.VarArea)
{
- newloc = RobustPredicates.FindCircumcenter(borg, bdest, bapex, ref xi, ref eta, behavior.offconstant);
+ newloc = predicates.FindCircumcenter(borg, bdest, bapex, ref xi, ref eta, behavior.offconstant);
}
else
{
diff --git a/Triangle.NET/Triangle/NewLocation.cs b/Triangle.NET/Triangle/NewLocation.cs
index 8be0ee5..204248a 100644
--- a/Triangle.NET/Triangle/NewLocation.cs
+++ b/Triangle.NET/Triangle/NewLocation.cs
@@ -22,6 +22,8 @@ namespace TriangleNet
{
const double EPS = 1e-50;
+ IPredicates predicates;
+
Mesh mesh;
Behavior behavior;
@@ -42,9 +44,11 @@ namespace TriangleNet
double[] poly2 = new double[100];
double[][] polys = new double[3][];
- public NewLocation(Mesh mesh)
+ public NewLocation(Mesh mesh, IPredicates predicates)
{
this.mesh = mesh;
+ this.predicates = predicates;
+
this.behavior = mesh.behavior;
}
@@ -178,7 +182,7 @@ namespace TriangleNet
// Use the counterclockwise() routine to ensure a positive (and
// reasonably accurate) result, avoiding any possibility of
// division by zero.
- denominator = 0.5 / RobustPredicates.CounterClockwise(tdest, tapex, torg);
+ denominator = 0.5 / predicates.CounterClockwise(tdest, tapex, torg);
// Don't count the above as an orientation test.
Statistic.CounterClockwiseCount--;
}
@@ -473,7 +477,7 @@ namespace TriangleNet
neighborvertex_2 = neighborotri.Dest();
neighborvertex_3 = neighborotri.Apex();
// now calculate neighbor's circumcenter which is the voronoi site
- neighborCircumcenter = RobustPredicates.FindCircumcenter(neighborvertex_1, neighborvertex_2, neighborvertex_3,
+ neighborCircumcenter = predicates.FindCircumcenter(neighborvertex_1, neighborvertex_2, neighborvertex_3,
ref xi_tmp, ref eta_tmp);
/// compute petal and Voronoi edge intersection ///
@@ -604,7 +608,7 @@ namespace TriangleNet
neighborvertex_2 = neighborotri.Dest();
neighborvertex_3 = neighborotri.Apex();
// now calculate neighbor's circumcenter which is the voronoi site
- neighborCircumcenter = RobustPredicates.FindCircumcenter(neighborvertex_1, neighborvertex_2, neighborvertex_3,
+ neighborCircumcenter = predicates.FindCircumcenter(neighborvertex_1, neighborvertex_2, neighborvertex_3,
ref xi_tmp, ref eta_tmp);
/// compute petal and Voronoi edge intersection ///
@@ -891,7 +895,7 @@ namespace TriangleNet
// Use the counterclockwise() routine to ensure a positive (and
// reasonably accurate) result, avoiding any possibility of
// division by zero.
- denominator = 0.5 / RobustPredicates.CounterClockwise(tdest, tapex, torg);
+ denominator = 0.5 / predicates.CounterClockwise(tdest, tapex, torg);
// Don't count the above as an orientation test.
Statistic.CounterClockwiseCount--;
}
@@ -1234,7 +1238,7 @@ namespace TriangleNet
neighborvertex_2 = neighborotri.Dest();
neighborvertex_3 = neighborotri.Apex();
// now calculate neighbor's circumcenter which is the voronoi site
- neighborCircumcenter = RobustPredicates.FindCircumcenter(neighborvertex_1, neighborvertex_2, neighborvertex_3,
+ neighborCircumcenter = predicates.FindCircumcenter(neighborvertex_1, neighborvertex_2, neighborvertex_3,
ref xi_tmp, ref eta_tmp);
/// compute petal and Voronoi edge intersection ///
@@ -1516,7 +1520,7 @@ namespace TriangleNet
neighborvertex_2 = neighborotri.Dest();
neighborvertex_3 = neighborotri.Apex();
// now calculate neighbor's circumcenter which is the voronoi site
- neighborCircumcenter = RobustPredicates.FindCircumcenter(neighborvertex_1, neighborvertex_2, neighborvertex_3,
+ neighborCircumcenter = predicates.FindCircumcenter(neighborvertex_1, neighborvertex_2, neighborvertex_3,
ref xi_tmp, ref eta_tmp);
/// compute petal and Voronoi edge intersection ///
@@ -3419,7 +3423,7 @@ namespace TriangleNet
int intFound = 0;
dx = x2 - x1;
dy = y2 - y1;
- numpolys = SplitConvexPolygon(numvertices, convexPoly, x1, y1, x2, y2, ref polys);
+ numpolys = SplitConvexPolygon(numvertices, convexPoly, x1, y1, x2, y2, polys);
if (numpolys == 3)
{
@@ -3480,7 +3484,7 @@ namespace TriangleNet
///
/// http://www.mathematik.uni-ulm.de/stochastik/lehre/ws03_04/rt/Geometry2D.ps
///
- private int SplitConvexPolygon(int numvertices, double[] convexPoly, double x1, double y1, double x2, double y2, ref double[][] polys)
+ private int SplitConvexPolygon(int numvertices, double[] convexPoly, double x1, double y1, double x2, double y2, double[][] polys)
{
// state = 0: before the first intersection (with the line)
// state = 1: after the first intersection (with the line)
@@ -4099,7 +4103,7 @@ namespace TriangleNet
else
{
// Orient 'searchtri' to fit the preconditions of calling preciselocate().
- ahead = RobustPredicates.CounterClockwise(torg, tdest, newvertex);
+ ahead = predicates.CounterClockwise(torg, tdest, newvertex);
if (ahead < 0.0)
{
// Turn around so that 'searchpoint' is to the left of the
diff --git a/Triangle.NET/Triangle/RobustPredicates.cs b/Triangle.NET/Triangle/RobustPredicates.cs
index 7c93801..2f2810e 100644
--- a/Triangle.NET/Triangle/RobustPredicates.cs
+++ b/Triangle.NET/Triangle/RobustPredicates.cs
@@ -25,8 +25,39 @@ namespace TriangleNet
/// command prompt with the command "CL /P /C EXACT.C", see
/// http://msdn.microsoft.com/en-us/library/8z9z0bx6.aspx
///
- public static class RobustPredicates
+ public class RobustPredicates : IPredicates
{
+ #region Default predicates instance (Singleton)
+
+ private static readonly object creationLock = new object();
+ private static RobustPredicates _default;
+
+ ///
+ /// Gets the default configuration instance.
+ ///
+ public static RobustPredicates Default
+ {
+ get
+ {
+ if (_default == null)
+ {
+ lock (creationLock)
+ {
+ if (_default == null)
+ {
+ _default = new RobustPredicates();
+ }
+ }
+ }
+
+ return _default;
+ }
+ }
+
+ #endregion
+
+ #region Static initialization
+
private static double epsilon, splitter, resulterrbound;
private static double ccwerrboundA, ccwerrboundB, ccwerrboundC;
private static double iccerrboundA, iccerrboundB, iccerrboundC;
@@ -49,7 +80,7 @@ namespace TriangleNet
///
/// Don't change this routine unless you fully understand it.
///
- public static void ExactInit()
+ static RobustPredicates()
{
double half;
double check, lastcheck;
@@ -89,6 +120,13 @@ namespace TriangleNet
//o3derrboundC = (26.0 + 288.0 * epsilon) * epsilon * epsilon;
}
+ #endregion
+
+ public RobustPredicates()
+ {
+ AllocateWorkspace();
+ }
+
///
/// Check, if the three points appear in counterclockwise order. The result is
/// also a rough approximation of twice the signed area of the triangle defined
@@ -100,7 +138,7 @@ namespace TriangleNet
/// Return a positive value if the points pa, pb, and pc occur in
/// counterclockwise order; a negative value if they occur in clockwise order;
/// and zero if they are collinear.
- public static double CounterClockwise(Point pa, Point pb, Point pc)
+ public double CounterClockwise(Point pa, Point pb, Point pc)
{
double detleft, detright, det;
double detsum, errbound;
@@ -165,7 +203,7 @@ namespace TriangleNet
/// Return a positive value if the point pd lies inside the circle passing through
/// pa, pb, and pc; a negative value if it lies outside; and zero if the four points
/// are cocircular.
- public static double InCircle(Point pa, Point pb, Point pc, Point pd)
+ public double InCircle(Point pa, Point pb, Point pc, Point pd)
{
double adx, bdx, cdx, ady, bdy, cdy;
double bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
@@ -231,7 +269,7 @@ namespace TriangleNet
/// Return a positive value if the point pd lies inside the circle passing through
/// pa, pb, and pc; a negative value if it lies outside; and zero if the four points
/// are cocircular.
- public static double NonRegular(Point pa, Point pb, Point pc, Point pd)
+ public double NonRegular(Point pa, Point pb, Point pc, Point pd)
{
return InCircle(pa, pb, pc, pd);
}
@@ -246,7 +284,7 @@ namespace TriangleNet
/// Relative coordinate of new location.
/// Off-center constant.
/// Coordinates of the circumcenter (or off-center)
- public static Point FindCircumcenter(Point org, Point dest, Point apex,
+ public Point FindCircumcenter(Point org, Point dest, Point apex,
ref double xi, ref double eta, double offconstant)
{
double xdo, ydo, xao, yao;
@@ -365,7 +403,7 @@ namespace TriangleNet
/// This procedure also returns the square of the length of the triangle's
/// shortest edge.
///
- public static Point FindCircumcenter(Point org, Point dest, Point apex,
+ public Point FindCircumcenter(Point org, Point dest, Point apex,
ref double xi, ref double eta)
{
double xdo, ydo, xao, yao;
@@ -429,7 +467,7 @@ namespace TriangleNet
/// property. (That is, if e is strongly nonoverlapping, h will be also.) Does NOT
/// maintain the nonoverlapping or nonadjacent properties.
///
- static int FastExpansionSumZeroElim(int elen, double[] e, int flen, double[] f, double[] h)
+ private int FastExpansionSumZeroElim(int elen, double[] e, int flen, double[] f, double[] h)
{
double Q;
double Qnew;
@@ -555,7 +593,7 @@ namespace TriangleNet
/// maintains the strongly nonoverlapping and nonadjacent properties as well. (That is,
/// if e has one of these properties, so will h.)
///
- static int ScaleExpansionZeroElim(int elen, double[] e, double b, double[] h)
+ private int ScaleExpansionZeroElim(int elen, double[] e, double b, double[] h)
{
double Q, sum;
double hh;
@@ -605,7 +643,7 @@ namespace TriangleNet
///
///
///
- static double Estimate(int elen, double[] e)
+ private double Estimate(int elen, double[] e)
{
double Q;
int eindex;
@@ -636,7 +674,7 @@ namespace TriangleNet
/// the returned value has the correct sign. Hence, this function is usually quite fast,
/// but will run more slowly when the input points are collinear or nearly so.
///
- static double CounterClockwiseAdapt(Point pa, Point pb, Point pc, double detsum)
+ private double CounterClockwiseAdapt(Point pa, Point pb, Point pc, double detsum)
{
double acx, acy, bcx, bcy;
double acxtail, acytail, bcxtail, bcytail;
@@ -649,7 +687,7 @@ namespace TriangleNet
double[] C1 = new double[8], C2 = new double[12], D = new double[16];
double B3;
int C1length, C2length, Dlength;
-
+
double u3;
double s1, t1;
double s0, t0;
@@ -741,7 +779,7 @@ namespace TriangleNet
/// the returned value has the correct sign. Hence, this function is usually quite fast,
/// but will run more slowly when the input points are cocircular or nearly so.
///
- static double InCircleAdapt(Point pa, Point pb, Point pc, Point pd, double permanent)
+ private double InCircleAdapt(Point pa, Point pb, Point pc, Point pd, double permanent)
{
double adx, bdx, cdx, ady, bdy, cdy;
double det, errbound;
@@ -750,15 +788,10 @@ namespace TriangleNet
double bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0;
double[] bc = new double[4], ca = new double[4], ab = new double[4];
double bc3, ca3, ab3;
- double[] axbc = new double[8], axxbc = new double[16], aybc = new double[8], ayybc = new double[16], adet = new double[32];
int axbclen, axxbclen, aybclen, ayybclen, alen;
- double[] bxca = new double[8], bxxca = new double[16], byca = new double[8], byyca = new double[16], bdet = new double[32];
int bxcalen, bxxcalen, bycalen, byycalen, blen;
- double[] cxab = new double[8], cxxab = new double[16], cyab = new double[8], cyyab = new double[16], cdet = new double[32];
int cxablen, cxxablen, cyablen, cyyablen, clen;
- double[] abdet = new double[64];
int ablen;
- double[] fin1 = new double[1152], fin2 = new double[1152];
double[] finnow, finother, finswap;
int finlength;
@@ -773,8 +806,6 @@ namespace TriangleNet
// See unsafe indexing in FastExpansionSumZeroElim.
double[] u = new double[5], v = new double[5];
double u3, v3;
- double[] temp8 = new double[8], temp16a = new double[16], temp16b = new double[16], temp16c = new double[16];
- double[] temp32a = new double[32], temp32b = new double[32], temp48 = new double[48], temp64 = new double[64];
int temp8len, temp16alen, temp16blen, temp16clen;
int temp32alen, temp32blen, temp48len, temp64len;
double[] axtbb = new double[8], axtcc = new double[8], aytbb = new double[8], aytcc = new double[8];
@@ -887,24 +918,21 @@ namespace TriangleNet
finnow = fin1;
finother = fin2;
- if ((bdxtail != 0.0) || (bdytail != 0.0)
- || (cdxtail != 0.0) || (cdytail != 0.0))
+ if ((bdxtail != 0.0) || (bdytail != 0.0) || (cdxtail != 0.0) || (cdytail != 0.0))
{
adxadx1 = (double)(adx * adx); c = (double)(splitter * adx); abig = (double)(c - adx); ahi = c - abig; alo = adx - ahi; err1 = adxadx1 - (ahi * ahi); err3 = err1 - ((ahi + ahi) * alo); adxadx0 = (alo * alo) - err3;
adyady1 = (double)(ady * ady); c = (double)(splitter * ady); abig = (double)(c - ady); ahi = c - abig; alo = ady - ahi; err1 = adyady1 - (ahi * ahi); err3 = err1 - ((ahi + ahi) * alo); adyady0 = (alo * alo) - err3;
_i = (double)(adxadx0 + adyady0); bvirt = (double)(_i - adxadx0); avirt = _i - bvirt; bround = adyady0 - bvirt; around = adxadx0 - avirt; aa[0] = around + bround; _j = (double)(adxadx1 + _i); bvirt = (double)(_j - adxadx1); avirt = _j - bvirt; bround = _i - bvirt; around = adxadx1 - avirt; _0 = around + bround; _i = (double)(_0 + adyady1); bvirt = (double)(_i - _0); avirt = _i - bvirt; bround = adyady1 - bvirt; around = _0 - avirt; aa[1] = around + bround; aa3 = (double)(_j + _i); bvirt = (double)(aa3 - _j); avirt = aa3 - bvirt; bround = _i - bvirt; around = _j - avirt; aa[2] = around + bround;
aa[3] = aa3;
}
- if ((cdxtail != 0.0) || (cdytail != 0.0)
- || (adxtail != 0.0) || (adytail != 0.0))
+ if ((cdxtail != 0.0) || (cdytail != 0.0) || (adxtail != 0.0) || (adytail != 0.0))
{
bdxbdx1 = (double)(bdx * bdx); c = (double)(splitter * bdx); abig = (double)(c - bdx); ahi = c - abig; alo = bdx - ahi; err1 = bdxbdx1 - (ahi * ahi); err3 = err1 - ((ahi + ahi) * alo); bdxbdx0 = (alo * alo) - err3;
bdybdy1 = (double)(bdy * bdy); c = (double)(splitter * bdy); abig = (double)(c - bdy); ahi = c - abig; alo = bdy - ahi; err1 = bdybdy1 - (ahi * ahi); err3 = err1 - ((ahi + ahi) * alo); bdybdy0 = (alo * alo) - err3;
_i = (double)(bdxbdx0 + bdybdy0); bvirt = (double)(_i - bdxbdx0); avirt = _i - bvirt; bround = bdybdy0 - bvirt; around = bdxbdx0 - avirt; bb[0] = around + bround; _j = (double)(bdxbdx1 + _i); bvirt = (double)(_j - bdxbdx1); avirt = _j - bvirt; bround = _i - bvirt; around = bdxbdx1 - avirt; _0 = around + bround; _i = (double)(_0 + bdybdy1); bvirt = (double)(_i - _0); avirt = _i - bvirt; bround = bdybdy1 - bvirt; around = _0 - avirt; bb[1] = around + bround; bb3 = (double)(_j + _i); bvirt = (double)(bb3 - _j); avirt = bb3 - bvirt; bround = _i - bvirt; around = _j - avirt; bb[2] = around + bround;
bb[3] = bb3;
}
- if ((adxtail != 0.0) || (adytail != 0.0)
- || (bdxtail != 0.0) || (bdytail != 0.0))
+ if ((adxtail != 0.0) || (adytail != 0.0) || (bdxtail != 0.0) || (bdytail != 0.0))
{
cdxcdx1 = (double)(cdx * cdx); c = (double)(splitter * cdx); abig = (double)(c - cdx); ahi = c - abig; alo = cdx - ahi; err1 = cdxcdx1 - (ahi * ahi); err3 = err1 - ((ahi + ahi) * alo); cdxcdx0 = (alo * alo) - err3;
cdycdy1 = (double)(cdy * cdy); c = (double)(splitter * cdy); abig = (double)(c - cdy); ahi = c - abig; alo = cdy - ahi; err1 = cdycdy1 - (ahi * ahi); err3 = err1 - ((ahi + ahi) * alo); cdycdy0 = (alo * alo) - err3;
@@ -1261,6 +1289,59 @@ namespace TriangleNet
return finnow[finlength - 1];
}
+ #region Workspace
+
+ // InCircleAdapt workspace:
+ double[] fin1, fin2, abdet;
+
+ double[] axbc, axxbc, aybc, ayybc, adet;
+ double[] bxca, bxxca, byca, byyca, bdet;
+ double[] cxab, cxxab, cyab, cyyab, cdet;
+
+ double[] temp8, temp16a, temp16b, temp16c;
+ double[] temp32a, temp32b, temp48, temp64;
+
+ private void AllocateWorkspace()
+ {
+ fin1 = new double[1152];
+ fin2 = new double[1152];
+ abdet = new double[64];
+
+ axbc = new double[8];
+ axxbc = new double[16];
+ aybc = new double[8];
+ ayybc = new double[16];
+ adet = new double[32];
+
+ bxca = new double[8];
+ bxxca = new double[16];
+ byca = new double[8];
+ byyca = new double[16];
+ bdet = new double[32];
+
+ cxab = new double[8];
+ cxxab = new double[16];
+ cyab = new double[8];
+ cyyab = new double[16];
+ cdet = new double[32];
+
+ temp8 = new double[8];
+ temp16a = new double[16];
+ temp16b = new double[16];
+ temp16c = new double[16];
+
+ temp32a = new double[32];
+ temp32b = new double[32];
+ temp48 = new double[48];
+ temp64 = new double[64];
+ }
+
+ private void ClearWorkspace()
+ {
+ }
+
+ #endregion
+
#endregion
}
}
diff --git a/Triangle.NET/Triangle/Smoothing/SimpleSmoother.cs b/Triangle.NET/Triangle/Smoothing/SimpleSmoother.cs
index 09a4c6b..c47f4c5 100644
--- a/Triangle.NET/Triangle/Smoothing/SimpleSmoother.cs
+++ b/Triangle.NET/Triangle/Smoothing/SimpleSmoother.cs
@@ -20,13 +20,30 @@ namespace TriangleNet.Smoothing
///
public class SimpleSmoother : ISmoother
{
+ IPredicates predicates;
+
+ IVoronoiFactory factory;
+
ConstraintOptions options;
///
/// Initializes a new instance of the class.
///
public SimpleSmoother()
+ : this(new VoronoiFactory(), RobustPredicates.Default)
{
+ }
+
+ ///
+ /// Initializes a new instance of the class.
+ ///
+ /// Voronoi object factory.
+ /// Geometric predicates implementation.
+ public SimpleSmoother(IVoronoiFactory factory, IPredicates predicates)
+ {
+ this.factory = factory;
+ this.predicates = predicates;
+
this.options = new ConstraintOptions() { ConformingDelaunay = true };
}
@@ -45,20 +62,22 @@ namespace TriangleNet.Smoothing
// Take a few smoothing rounds (Lloyd's algorithm).
for (int i = 0; i < limit; i++)
{
- Step(smoothedMesh);
+ Step(smoothedMesh, factory);
// Actually, we only want to rebuild, if mesh is no longer
// Delaunay. Flipping edges could be the right choice instead
// of re-triangulating...
smoothedMesh = (Mesh)Rebuild(smoothedMesh).Triangulate(options);
+
+ factory.Reset();
}
smoothedMesh.CopyTo((Mesh)mesh);
}
- private void Step(Mesh mesh)
+ private void Step(Mesh mesh, IVoronoiFactory factory)
{
- var voronoi = new BoundedVoronoi(mesh);
+ var voronoi = new BoundedVoronoi(mesh, factory, predicates);
double x, y;
diff --git a/Triangle.NET/Triangle/Smoothing/VoronoiFactory.cs b/Triangle.NET/Triangle/Smoothing/VoronoiFactory.cs
new file mode 100644
index 0000000..8be31aa
--- /dev/null
+++ b/Triangle.NET/Triangle/Smoothing/VoronoiFactory.cs
@@ -0,0 +1,201 @@
+
+namespace TriangleNet.Smoothing
+{
+ using System;
+ using TriangleNet.Topology.DCEL;
+ using TriangleNet.Voronoi;
+
+ ///
+ /// Factory which re-uses objects in the smoothing loop to enhance performance.
+ ///
+ ///
+ /// See .
+ ///
+ class VoronoiFactory : IVoronoiFactory
+ {
+ ObjectPool vertices;
+ ObjectPool edges;
+ ObjectPool faces;
+
+ public VoronoiFactory()
+ {
+ vertices = new ObjectPool();
+ edges = new ObjectPool();
+ faces = new ObjectPool();
+ }
+
+ public void Initialize(int vertexCount, int edgeCount, int faceCount)
+ {
+ vertices.Capacity = vertexCount;
+ edges.Capacity = edgeCount;
+ faces.Capacity = faceCount;
+
+ for (int i = vertices.Count; i < vertexCount; i++)
+ {
+ vertices.Put(new Vertex(0, 0));
+ }
+
+
+ for (int i = edges.Count; i < edgeCount; i++)
+ {
+ edges.Put(new HalfEdge(null));
+ }
+
+ for (int i = faces.Count; i < faceCount; i++)
+ {
+ faces.Put(new Face(null));
+ }
+
+ Reset();
+ }
+
+ public void Reset()
+ {
+ vertices.Release();
+ edges.Release();
+ faces.Release();
+ }
+
+ public Vertex CreateVertex(double x, double y)
+ {
+ Vertex vertex;
+
+ if (vertices.TryGet(out vertex))
+ {
+ vertex.x = x;
+ vertex.y = y;
+ vertex.leaving = null;
+
+ return vertex;
+ }
+
+ vertex = new Vertex(x, y);
+
+ vertices.Put(vertex);
+
+ return vertex;
+ }
+
+ public HalfEdge CreateHalfEdge(Vertex origin, Face face)
+ {
+ HalfEdge edge;
+
+ if (edges.TryGet(out edge))
+ {
+ edge.origin = origin;
+ edge.face = face;
+ edge.next = null;
+ edge.twin = null;
+
+ if (face != null && face.edge == null)
+ {
+ face.edge = edge;
+ }
+
+ return edge;
+ }
+
+ edge = new HalfEdge(origin, face);
+
+ edges.Put(edge);
+
+ return edge;
+ }
+
+ public Face CreateFace(Geometry.Vertex vertex)
+ {
+ Face face;
+
+ if (faces.TryGet(out face))
+ {
+ face.id = vertex.id;
+ face.generator = vertex;
+ face.edge = null;
+
+ return face;
+ }
+
+ face = new Face(vertex);
+
+ faces.Put(face);
+
+ return face;
+ }
+
+ class ObjectPool where T : class
+ {
+ int index, count;
+
+ T[] pool;
+
+ public int Count
+ {
+ get { return count; }
+ }
+
+
+ public int Capacity
+ {
+ get { return this.pool.Length; }
+ set { Resize(value); }
+ }
+
+ public ObjectPool(int capacity = 3)
+ {
+ this.index = 0;
+ this.count = 0;
+
+ this.pool = new T[capacity];
+ }
+
+ public ObjectPool(T[] pool)
+ {
+ this.index = 0;
+ this.count = 0;
+
+ this.pool = pool;
+ }
+
+ public bool TryGet(out T obj)
+ {
+ if (this.index < this.count)
+ {
+ obj = this.pool[this.index++];
+
+ return true;
+ }
+
+ obj = null;
+
+ return false;
+ }
+
+ public void Put(T obj)
+ {
+ var capacity = this.pool.Length;
+
+ if (capacity <= this.count)
+ {
+ Resize(2 * capacity);
+ }
+
+ this.pool[this.count++] = obj;
+
+ this.index++;
+ }
+
+ public void Release()
+ {
+ this.index = 0;
+ }
+
+ private void Resize(int size)
+ {
+ if (size > this.count)
+ {
+ Array.Resize(ref this.pool, size);
+ }
+ }
+ }
+ }
+}
diff --git a/Triangle.NET/Triangle/Triangle.csproj b/Triangle.NET/Triangle/Triangle.csproj
index be22a09..60aa0cf 100644
--- a/Triangle.NET/Triangle/Triangle.csproj
+++ b/Triangle.NET/Triangle/Triangle.csproj
@@ -62,6 +62,7 @@
+
@@ -80,6 +81,7 @@
+
@@ -99,6 +101,8 @@
+
+
diff --git a/Triangle.NET/Triangle/TriangleLocator.cs b/Triangle.NET/Triangle/TriangleLocator.cs
index 28b848b..391521c 100644
--- a/Triangle.NET/Triangle/TriangleLocator.cs
+++ b/Triangle.NET/Triangle/TriangleLocator.cs
@@ -22,13 +22,16 @@ namespace TriangleNet
Sampler sampler;
Mesh mesh;
+ IPredicates predicates;
+
// Pointer to a recently visited triangle. Improves point location if
// proximate vertices are inserted sequentially.
internal Otri recenttri;
- public TriangleLocator(Mesh mesh)
+ public TriangleLocator(Mesh mesh, IPredicates predicates)
{
this.mesh = mesh;
+ this.predicates = predicates;
sampler = new Sampler();
}
@@ -133,10 +136,10 @@ namespace TriangleNet
}
// Does the point lie on the other side of the line defined by the
// triangle edge opposite the triangle's destination?
- destorient = RobustPredicates.CounterClockwise(forg, fapex, searchpoint);
+ destorient = predicates.CounterClockwise(forg, fapex, searchpoint);
// Does the point lie on the other side of the line defined by the
// triangle edge opposite the triangle's origin?
- orgorient = RobustPredicates.CounterClockwise(fapex, fdest, searchpoint);
+ orgorient = predicates.CounterClockwise(fapex, fdest, searchpoint);
if (destorient > 0.0)
{
if (orgorient > 0.0)
@@ -322,7 +325,7 @@ namespace TriangleNet
return LocateResult.OnVertex;
}
// Orient 'searchtri' to fit the preconditions of calling preciselocate().
- ahead = RobustPredicates.CounterClockwise(torg, tdest, searchpoint);
+ ahead = predicates.CounterClockwise(torg, tdest, searchpoint);
if (ahead < 0.0)
{
// Turn around so that 'searchpoint' is to the left of the
diff --git a/Triangle.NET/Triangle/Voronoi/BoundedVoronoi.cs b/Triangle.NET/Triangle/Voronoi/BoundedVoronoi.cs
index 2700b4f..dd7f6ca 100644
--- a/Triangle.NET/Triangle/Voronoi/BoundedVoronoi.cs
+++ b/Triangle.NET/Triangle/Voronoi/BoundedVoronoi.cs
@@ -19,7 +19,12 @@ namespace TriangleNet.Voronoi
int offset;
public BoundedVoronoi(Mesh mesh)
- : base(mesh, true)
+ : this(mesh, new DefaultVoronoiFactory(), RobustPredicates.Default)
+ {
+ }
+
+ public BoundedVoronoi(Mesh mesh, IVoronoiFactory factory, IPredicates predicates)
+ : base(mesh, factory, predicates, true)
{
// We explicitly told the base constructor to call the Generate method, so
// at this point the basic Voronoi diagram is already created.
@@ -46,7 +51,7 @@ namespace TriangleNet.Voronoi
var v1 = (TVertex)edge.face.generator;
var v2 = (TVertex)twin.face.generator;
- double dir = RobustPredicates.CounterClockwise(v1, v2, edge.origin);
+ double dir = predicates.CounterClockwise(v1, v2, edge.origin);
if (dir <= 0)
{
@@ -75,10 +80,10 @@ namespace TriangleNet.Voronoi
v.y = (v1.y + v2.y) / 2.0;
// Close the cell connected to edge.
- var gen = new HVertex(v1.x, v1.y);
+ var gen = factory.CreateVertex(v1.x, v1.y);
- var h1 = new HalfEdge(edge.twin.origin, edge.face);
- var h2 = new HalfEdge(gen, edge.face);
+ var h1 = factory.CreateHalfEdge(edge.twin.origin, edge.face);
+ var h2 = factory.CreateHalfEdge(gen, edge.face);
edge.next = h1;
h1.next = h2;
@@ -129,8 +134,8 @@ namespace TriangleNet.Voronoi
edge.twin = null;
// Close the cell.
- var gen = new HVertex(v1.x, v1.y);
- var he = new HalfEdge(gen, edge.face);
+ var gen = factory.CreateVertex(v1.x, v1.y);
+ var he = factory.CreateHalfEdge(gen, edge.face);
edge.next = he;
he.next = edge.face.edge;
diff --git a/Triangle.NET/Triangle/Voronoi/DefaultVoronoiFactory.cs b/Triangle.NET/Triangle/Voronoi/DefaultVoronoiFactory.cs
new file mode 100644
index 0000000..6702ff0
--- /dev/null
+++ b/Triangle.NET/Triangle/Voronoi/DefaultVoronoiFactory.cs
@@ -0,0 +1,35 @@
+
+namespace TriangleNet.Voronoi
+{
+ using System;
+ using TriangleNet.Topology.DCEL;
+
+ ///
+ /// Default factory for Voronoi / DCEL mesh objects.
+ ///
+ public class DefaultVoronoiFactory : IVoronoiFactory
+ {
+ public void Initialize(int vertexCount, int edgeCount, int faceCount)
+ {
+ }
+
+ public void Reset()
+ {
+ }
+
+ public Vertex CreateVertex(double x, double y)
+ {
+ return new Vertex(x, y);
+ }
+
+ public HalfEdge CreateHalfEdge(Vertex origin, Face face)
+ {
+ return new HalfEdge(origin, face);
+ }
+
+ public Face CreateFace(Geometry.Vertex vertex)
+ {
+ return new Face(vertex);
+ }
+ }
+}
diff --git a/Triangle.NET/Triangle/Voronoi/IVoronoiFactory.cs b/Triangle.NET/Triangle/Voronoi/IVoronoiFactory.cs
new file mode 100644
index 0000000..47276e0
--- /dev/null
+++ b/Triangle.NET/Triangle/Voronoi/IVoronoiFactory.cs
@@ -0,0 +1,18 @@
+
+namespace TriangleNet.Voronoi
+{
+ using TriangleNet.Topology.DCEL;
+
+ public interface IVoronoiFactory
+ {
+ void Initialize(int vertexCount, int edgeCount, int faceCount);
+
+ void Reset();
+
+ Vertex CreateVertex(double x, double y);
+
+ HalfEdge CreateHalfEdge(Vertex origin, Face face);
+
+ Face CreateFace(Geometry.Vertex vertex);
+ }
+}
diff --git a/Triangle.NET/Triangle/Voronoi/Legacy/BoundedVoronoiLegacy.cs b/Triangle.NET/Triangle/Voronoi/Legacy/BoundedVoronoiLegacy.cs
index ace56aa..47ccc39 100644
--- a/Triangle.NET/Triangle/Voronoi/Legacy/BoundedVoronoiLegacy.cs
+++ b/Triangle.NET/Triangle/Voronoi/Legacy/BoundedVoronoiLegacy.cs
@@ -21,6 +21,8 @@ namespace TriangleNet.Voronoi.Legacy
[Obsolete("Use TriangleNet.Voronoi.BoundedVoronoi class instead.")]
public class BoundedVoronoiLegacy : IVoronoi
{
+ IPredicates predicates = RobustPredicates.Default;
+
Mesh mesh;
Point[] points;
@@ -132,7 +134,7 @@ namespace TriangleNet.Voronoi.Legacy
{
tri.tri = item;
- pt = RobustPredicates.FindCircumcenter(tri.Org(), tri.Dest(), tri.Apex(), ref xi, ref eta);
+ pt = predicates.FindCircumcenter(tri.Org(), tri.Dest(), tri.Apex(), ref xi, ref eta);
pt.id = item.id;
points[item.id] = pt;
diff --git a/Triangle.NET/Triangle/Voronoi/Legacy/SimpleVoronoi.cs b/Triangle.NET/Triangle/Voronoi/Legacy/SimpleVoronoi.cs
index 7a54adf..17fd0ce 100644
--- a/Triangle.NET/Triangle/Voronoi/Legacy/SimpleVoronoi.cs
+++ b/Triangle.NET/Triangle/Voronoi/Legacy/SimpleVoronoi.cs
@@ -18,6 +18,8 @@ namespace TriangleNet.Voronoi.Legacy
[Obsolete("Use TriangleNet.Voronoi.StandardVoronoi class instead.")]
public class SimpleVoronoi : IVoronoi
{
+ IPredicates predicates = RobustPredicates.Default;
+
Mesh mesh;
Point[] points;
@@ -120,7 +122,7 @@ namespace TriangleNet.Voronoi.Legacy
{
tri.tri = item;
- pt = RobustPredicates.FindCircumcenter(tri.Org(), tri.Dest(), tri.Apex(), ref xi, ref eta);
+ pt = predicates.FindCircumcenter(tri.Org(), tri.Dest(), tri.Apex(), ref xi, ref eta);
pt.id = item.id;
points[item.id] = pt;
diff --git a/Triangle.NET/Triangle/Voronoi/StandardVoronoi.cs b/Triangle.NET/Triangle/Voronoi/StandardVoronoi.cs
index 54c69b3..b2ac312 100644
--- a/Triangle.NET/Triangle/Voronoi/StandardVoronoi.cs
+++ b/Triangle.NET/Triangle/Voronoi/StandardVoronoi.cs
@@ -14,12 +14,17 @@ namespace TriangleNet.Voronoi
public class StandardVoronoi : VoronoiBase
{
public StandardVoronoi(Mesh mesh)
- : this(mesh, mesh.bounds)
+ : this(mesh, mesh.bounds, new DefaultVoronoiFactory(), RobustPredicates.Default)
{
}
public StandardVoronoi(Mesh mesh, Rectangle box)
- : base(mesh, true)
+ : this(mesh, box, new DefaultVoronoiFactory(), RobustPredicates.Default)
+ {
+ }
+
+ public StandardVoronoi(Mesh mesh, Rectangle box, IVoronoiFactory factory, IPredicates predicates)
+ : base(mesh, factory, predicates, true)
{
// We assume the box to be at least as large as the mesh.
box.Expand(mesh.bounds);
diff --git a/Triangle.NET/Triangle/Voronoi/VoronoiBase.cs b/Triangle.NET/Triangle/Voronoi/VoronoiBase.cs
index 0741d4f..5a94332 100644
--- a/Triangle.NET/Triangle/Voronoi/VoronoiBase.cs
+++ b/Triangle.NET/Triangle/Voronoi/VoronoiBase.cs
@@ -20,6 +20,10 @@ namespace TriangleNet.Voronoi
///
public abstract class VoronoiBase : DcelMesh
{
+ protected IPredicates predicates;
+
+ protected IVoronoiFactory factory;
+
// List of infinite half-edges, i.e. half-edges that start at circumcenters of triangles
// which lie on the domain boundary.
protected List rays;
@@ -28,11 +32,16 @@ namespace TriangleNet.Voronoi
/// Initializes a new instance of the class.
///
/// Triangle mesh.
+ /// Voronoi object factory.
+ /// Geometric predicates implementation.
/// If set to true, the constuctor will call the Generate
/// method, which builds the Voronoi diagram.
- protected VoronoiBase(Mesh mesh, bool generate)
- : base(false)
+ protected VoronoiBase(Mesh mesh, IVoronoiFactory factory, IPredicates predicates,
+ bool generate) : base(false)
{
+ this.factory = factory;
+ this.predicates = predicates;
+
if (generate)
{
Generate(mesh);
@@ -55,13 +64,20 @@ namespace TriangleNet.Voronoi
var vertices = new Vertex[mesh.triangles.Count + mesh.hullsize];
var faces = new Face[mesh.vertices.Count];
+ if (factory == null)
+ {
+ factory = new DefaultVoronoiFactory();
+ }
+
+ factory.Initialize(vertices.Length, 2 * mesh.edges, faces.Length);
+
// Compute triangles circumcenters.
var map = ComputeVertices(mesh, vertices);
// Create all Voronoi faces.
foreach (var vertex in mesh.vertices.Values)
{
- faces[vertex.id] = new Face(vertex);
+ faces[vertex.id] = factory.CreateFace(vertex);
}
ComputeEdges(mesh, vertices, faces, map);
@@ -94,9 +110,9 @@ namespace TriangleNet.Voronoi
id = t.id;
tri.tri = t;
- pt = RobustPredicates.FindCircumcenter(tri.Org(), tri.Dest(), tri.Apex(), ref xi, ref eta);
+ pt = predicates.FindCircumcenter(tri.Org(), tri.Dest(), tri.Apex(), ref xi, ref eta);
- vertex = new Vertex(pt.x, pt.y);
+ vertex = factory.CreateVertex(pt.x, pt.y);
vertex.id = id;
vertices[id] = vertex;
@@ -169,13 +185,13 @@ namespace TriangleNet.Voronoi
px = dest.y - org.y;
py = org.x - dest.x;
- end = new Vertex(vertex.x + px, vertex.y + py);
+ end = factory.CreateVertex(vertex.x + px, vertex.y + py);
end.id = count + j++;
vertices[end.id] = end;
- edge = new HalfEdge(end, face);
- twin = new HalfEdge(vertex, neighborFace);
+ edge = factory.CreateHalfEdge(end, face);
+ twin = factory.CreateHalfEdge(vertex, neighborFace);
// Make (face.edge) always point to an edge that starts at an infinite
// vertex. This will allow traversing of unbounded faces.
@@ -191,8 +207,8 @@ namespace TriangleNet.Voronoi
end = vertices[nid];
// Create half-edges.
- edge = new HalfEdge(end, face);
- twin = new HalfEdge(vertex, neighborFace);
+ edge = factory.CreateHalfEdge(end, face);
+ twin = factory.CreateHalfEdge(vertex, neighborFace);
// Add to vertex map.
map[nid].Add(edge);