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

Add interface for geometric predicates

Browse Source 
Add factory for Voronoi mesh entities (reduces allocations in SimpleSmoother)
Fix issue #10959

git-svn-id: https://triangle.svn.codeplex.com/svn@77246 0e2699bc-83d4-4a8f-98e7-55e24ab8c7a5
This commit is contained in:
SND\wo80_cp
2015-09-10 23:05:25 +00:00
parent da90ca5d0e
commit a2a6d8c839
24 changed files with 591 additions and 273 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Triangle.NET/Triangle/IO/TriangleWriter.cs
-162
View File
@@ -453,167 +453,5 @@ namespace TriangleNet.IO
}
}
}
/// <summary>
/// Write the Voronoi diagram to a .voro file.
/// </summary>
/// <param name="mesh"></param>
/// <param name="filename"></param>
/// <returns></returns>
/// <remarks>
/// 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.</remarks>
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++;
}
}
}
}
}
/// <summary>
/// Write the triangulation to an .off file.
/// </summary>
/// <param name="mesh"></param>
/// <param name="filename"></param>
/// <remarks>
/// OFF stands for the Object File Format, a format used by the Geometry
/// Center's Geomview package.
/// </remarks>
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);
}
}
}
}
}
Triangle.NET/Triangle/IPredicates.cs
+22
View File
@@ -0,0 +1,22 @@
// -----------------------------------------------------------------------
// <copyright file="IPredicates.cs">
// Triangle.NET code by Christian Woltering, http://triangle.codeplex.com/
// </copyright>
// -----------------------------------------------------------------------
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);
}
}
Triangle.NET/Triangle/Mesh.cs
+11 -9
View File
@@ -24,6 +24,8 @@ namespace TriangleNet
{
#region Variables
IPredicates predicates;
ILog<LogItem> logger;
QualityMesher qualityMesher;
@@ -221,7 +223,7 @@ namespace TriangleNet
/// <summary>
/// Initializes a new instance of the <see cref="Mesh" /> class.
/// </summary>
public Mesh()
public Mesh(IPredicates predicates)
{
Initialize();
@@ -238,11 +240,11 @@ namespace TriangleNet
holes = new List<Point>();
regions = new List<RegionPointer>();
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;
Triangle.NET/Triangle/MeshValidator.cs
+4 -2
View File
@@ -13,6 +13,8 @@ namespace TriangleNet
public static class MeshValidator
{
private static RobustPredicates predicates = RobustPredicates.Default;
/// <summary>
/// Test the mesh for topological consistency.
/// </summary>
@@ -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)
{
Triangle.NET/Triangle/Meshing/Algorithm/Dwyer.cs
+25 -12
View File
@@ -45,13 +45,26 @@ namespace TriangleNet.Meshing.Algorithm
/// </remarks>
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;
}
/// <summary>
/// Form a Delaunay triangulation by the divide-and-conquer method.
/// </summary>
@@ -62,7 +75,7 @@ namespace TriangleNet.Meshing.Algorithm
/// </remarks>
public IMesh Triangulate(IList<Vertex> 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.
Triangle.NET/Triangle/Meshing/Algorithm/Incremental.cs
+13 -1
View File
@@ -16,8 +16,20 @@ namespace TriangleNet.Meshing.Algorithm
/// </summary>
public class Incremental : ITriangulator
{
IPredicates predicates;
Mesh mesh;
public Incremental()
: this(RobustPredicates.Default)
{
}
public Incremental(IPredicates predicates)
{
this.predicates = predicates;
}
/// <summary>
/// Form a Delaunay triangulation by incrementally inserting vertices.
/// </summary>
@@ -25,7 +37,7 @@ namespace TriangleNet.Meshing.Algorithm
/// triangulation.</returns>
public IMesh Triangulate(IList<Vertex> points)
{
this.mesh = new Mesh();
this.mesh = new Mesh(predicates);
this.mesh.TransferNodes(points);
Otri starttri = new Otri();
Triangle.NET/Triangle/Meshing/Algorithm/SweepLine.cs
+16 -4
View File
@@ -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<SplayNode> splaynodes;
public SweepLine()
: this(RobustPredicates.Default)
{
}
public SweepLine(IPredicates predicates)
{
this.predicates = predicates;
}
public IMesh Triangulate(IList<Vertex> 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;
Triangle.NET/Triangle/Meshing/ConstraintMesher.cs
+16 -12
View File
@@ -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<LogItem> 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.
Triangle.NET/Triangle/Meshing/Converter.cs
+1 -1
View File
@@ -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);
Triangle.NET/Triangle/Meshing/GenericMesher.cs
+16 -2
View File
@@ -1,4 +1,9 @@

// -----------------------------------------------------------------------
// <copyright file="GenericMesher.cs">
// Triangle.NET code by Christian Woltering, http://triangle.codeplex.com/
// </copyright>
// -----------------------------------------------------------------------
namespace TriangleNet.Meshing
{
using System;
@@ -12,15 +17,24 @@ namespace TriangleNet.Meshing
/// </summary>
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;
}
Triangle.NET/Triangle/Meshing/QualityMesher.cs
+8 -4
View File
@@ -19,6 +19,8 @@ namespace TriangleNet.Meshing
/// </summary>
class QualityMesher
{
IPredicates predicates;
Queue<BadSubseg> badsubsegs;
BadTriQueue queue;
Mesh mesh;
@@ -28,7 +30,7 @@ namespace TriangleNet.Meshing
ILog<LogItem> 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);
}
/// <summary>
@@ -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
{
Triangle.NET/Triangle/NewLocation.cs
+14 -10
View File
@@ -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
/// <remarks>
/// http://www.mathematik.uni-ulm.de/stochastik/lehre/ws03_04/rt/Geometry2D.ps
/// </remarks>
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
Triangle.NET/Triangle/RobustPredicates.cs
+107 -26
View File
@@ -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
/// </remarks>
public static class RobustPredicates
public class RobustPredicates : IPredicates
{
#region Default predicates instance (Singleton)
private static readonly object creationLock = new object();
private static RobustPredicates _default;
/// <summary>
/// Gets the default configuration instance.
/// </summary>
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.
/// </remarks>
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();
}
/// <summary>
/// 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
/// <returns>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.</returns>
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
/// <returns>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.</returns>
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
/// <returns>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.</returns>
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
/// <param name="eta">Relative coordinate of new location.</param>
/// <param name="offconstant">Off-center constant.</param>
/// <returns>Coordinates of the circumcenter (or off-center)</returns>
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.
/// </remarks>
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.
/// </remarks>
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.)
/// </remarks>
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
/// <param name="elen"></param>
/// <param name="e"></param>
/// <returns></returns>
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.
/// </remarks>
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.
/// </remarks>
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
}
}
Triangle.NET/Triangle/Smoothing/SimpleSmoother.cs
+22 -3
View File
@@ -20,13 +20,30 @@ namespace TriangleNet.Smoothing
/// </remarks>
public class SimpleSmoother : ISmoother
{
IPredicates predicates;
IVoronoiFactory factory;
ConstraintOptions options;
/// <summary>
/// Initializes a new instance of the <see cref="SimpleSmoother" /> class.
/// </summary>
public SimpleSmoother()
: this(new VoronoiFactory(), RobustPredicates.Default)
{
}
/// <summary>
/// Initializes a new instance of the <see cref="SimpleSmoother" /> class.
/// </summary>
/// <param name="factory">Voronoi object factory.</param>
/// <param name="predicates">Geometric predicates implementation.</param>
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;
Triangle.NET/Triangle/Smoothing/VoronoiFactory.cs
+201
View File
@@ -0,0 +1,201 @@

namespace TriangleNet.Smoothing
{
using System;
using TriangleNet.Topology.DCEL;
using TriangleNet.Voronoi;
/// <summary>
/// Factory which re-uses objects in the smoothing loop to enhance performance.
/// </summary>
/// <remarks>
/// See <see cref="SimpleSmoother"/>.
/// </remarks>
class VoronoiFactory : IVoronoiFactory
{
ObjectPool<Vertex> vertices;
ObjectPool<HalfEdge> edges;
ObjectPool<Face> faces;
public VoronoiFactory()
{
vertices = new ObjectPool<Vertex>();
edges = new ObjectPool<HalfEdge>();
faces = new ObjectPool<Face>();
}
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<T> 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);
}
}
}
}
}
Triangle.NET/Triangle/Triangle.csproj
+4
View File
@@ -62,6 +62,7 @@
<Compile Include="IO\TriangleFormat.cs" />
<Compile Include="IO\TriangleReader.cs" />
<Compile Include="IO\TriangleWriter.cs" />
<Compile Include="IPredicates.cs" />
<Compile Include="Meshing\ConstraintMesher.cs" />
<Compile Include="Meshing\ConstraintOptions.cs" />
<Compile Include="Meshing\Converter.cs" />
@@ -80,6 +81,7 @@
<Compile Include="Meshing\Data\BadTriQueue.cs" />
<Compile Include="Meshing\Iterators\EdgeIterator.cs" />
<Compile Include="Meshing\Iterators\RegionIterator.cs" />
<Compile Include="Smoothing\VoronoiFactory.cs" />
<Compile Include="Tools\AdjacencyMatrix.cs" />
<Compile Include="Tools\CuthillMcKee.cs" />
<Compile Include="Tools\IntersectionHelper.cs" />
@@ -99,6 +101,8 @@
<Compile Include="Smoothing\ISmoother.cs" />
<Compile Include="Smoothing\SimpleSmoother.cs" />
<Compile Include="Voronoi\BoundedVoronoi.cs" />
<Compile Include="Voronoi\DefaultVoronoiFactory.cs" />
<Compile Include="Voronoi\IVoronoiFactory.cs" />
<Compile Include="Voronoi\StandardVoronoi.cs" />
<Compile Include="Voronoi\VoronoiBase.cs" />
<Compile Include="Voronoi\Legacy\BoundedVoronoiLegacy.cs" />
Triangle.NET/Triangle/TriangleLocator.cs
+7 -4
View File
@@ -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
Triangle.NET/Triangle/Voronoi/BoundedVoronoi.cs
+12 -7
View File
@@ -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;
Triangle.NET/Triangle/Voronoi/DefaultVoronoiFactory.cs
+35
View File
@@ -0,0 +1,35 @@

namespace TriangleNet.Voronoi
{
using System;
using TriangleNet.Topology.DCEL;
/// <summary>
/// Default factory for Voronoi / DCEL mesh objects.
/// </summary>
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);
}
}
}
Triangle.NET/Triangle/Voronoi/IVoronoiFactory.cs
+18
View File
@@ -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);
}
}
Triangle.NET/Triangle/Voronoi/Legacy/BoundedVoronoiLegacy.cs
+3 -1
View File
@@ -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;
Triangle.NET/Triangle/Voronoi/Legacy/SimpleVoronoi.cs
+3 -1
View File
@@ -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;
Triangle.NET/Triangle/Voronoi/StandardVoronoi.cs
+7 -2
View File
@@ -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);
Triangle.NET/Triangle/Voronoi/VoronoiBase.cs
+26 -10
View File
@@ -20,6 +20,10 @@ namespace TriangleNet.Voronoi
/// </summary>
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<HalfEdge> rays;
@@ -28,11 +32,16 @@ namespace TriangleNet.Voronoi
/// Initializes a new instance of the <see cref="VoronoiBase" /> class.
/// </summary>
/// <param name="mesh">Triangle mesh.</param>
/// <param name="factory">Voronoi object factory.</param>
/// <param name="predicates">Geometric predicates implementation.</param>
/// <param name="generate">If set to true, the constuctor will call the Generate
/// method, which builds the Voronoi diagram.</param>
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);