Added QuadTree for point location

Added maximum angle to test app

git-svn-id: https://triangle.svn.codeplex.com/svn@73189 0e2699bc-83d4-4a8f-98e7-55e24ab8c7a5
This commit is contained in:
SND\wo80_cp
2013-06-18 09:38:34 +00:00
parent b8972bfcd3
commit e29212329d
6 changed files with 529 additions and 33 deletions
+17
View File
@@ -452,6 +452,13 @@ namespace MeshExplorer
mesh.Behavior.MinAngle = meshControlView.ParamMinAngleValue;
double maxAngle = meshControlView.ParamMaxAngleValue;
if (maxAngle < 180)
{
mesh.Behavior.MaxAngle = maxAngle;
}
// Ignore area constraints on initial triangulation.
//double area = slMaxArea.Value * 0.01;
@@ -507,6 +514,13 @@ namespace MeshExplorer
mesh.Behavior.MinAngle = meshControlView.ParamMinAngleValue;
double maxAngle = meshControlView.ParamMaxAngleValue;
if (maxAngle < 180)
{
mesh.Behavior.MaxAngle = maxAngle;
}
try
{
sw.Start();
@@ -715,7 +729,10 @@ namespace MeshExplorer
if (mesh != null)
{
bool isConsistent, isDelaunay;
Behavior.Verbose = true;
mesh.Check(out isConsistent, out isDelaunay);
Behavior.Verbose = false;
ShowLog();
}
+1 -1
View File
@@ -79,7 +79,7 @@
this.label19.Name = "label19";
this.label19.Size = new System.Drawing.Size(134, 40);
this.label19.TabIndex = 6;
this.label19.Text = "Beta 3 (2013-01-20)\r\nChristian Woltering\r\nMIT";
this.label19.Text = "Beta 3 (2013-06-18)\r\nChristian Woltering\r\nMIT";
//
// label18
//
+75 -32
View File
@@ -35,20 +35,23 @@
this.label9 = new System.Windows.Forms.Label();
this.label8 = new System.Windows.Forms.Label();
this.label5 = new System.Windows.Forms.Label();
this.label1 = new System.Windows.Forms.Label();
this.label2 = new System.Windows.Forms.Label();
this.lbMaxAngle = new System.Windows.Forms.Label();
this.cbSweepline = new MeshExplorer.Controls.DarkCheckBox();
this.slMaxAngle = new MeshExplorer.Controls.DarkSlider();
this.slMaxArea = new MeshExplorer.Controls.DarkSlider();
this.slMinAngle = new MeshExplorer.Controls.DarkSlider();
this.cbConformDel = new MeshExplorer.Controls.DarkCheckBox();
this.cbConvex = new MeshExplorer.Controls.DarkCheckBox();
this.cbQuality = new MeshExplorer.Controls.DarkCheckBox();
this.cbSweepline = new MeshExplorer.Controls.DarkCheckBox();
this.label1 = new System.Windows.Forms.Label();
this.SuspendLayout();
//
// lbMaxArea
//
this.lbMaxArea.AutoSize = true;
this.lbMaxArea.ForeColor = System.Drawing.Color.White;
this.lbMaxArea.Location = new System.Drawing.Point(227, 65);
this.lbMaxArea.Location = new System.Drawing.Point(227, 87);
this.lbMaxArea.Name = "lbMaxArea";
this.lbMaxArea.Size = new System.Drawing.Size(13, 13);
this.lbMaxArea.TabIndex = 20;
@@ -58,7 +61,7 @@
//
this.label6.AutoSize = true;
this.label6.ForeColor = System.Drawing.Color.White;
this.label6.Location = new System.Drawing.Point(8, 64);
this.label6.Location = new System.Drawing.Point(8, 86);
this.label6.Name = "label6";
this.label6.Size = new System.Drawing.Size(81, 13);
this.label6.TabIndex = 23;
@@ -78,7 +81,7 @@
//
this.label23.BackColor = System.Drawing.Color.DimGray;
this.label23.ForeColor = System.Drawing.Color.DarkGray;
this.label23.Location = new System.Drawing.Point(8, 151);
this.label23.Location = new System.Drawing.Point(8, 173);
this.label23.Name = "label23";
this.label23.Size = new System.Drawing.Size(251, 33);
this.label23.TabIndex = 24;
@@ -89,7 +92,7 @@
//
this.label9.BackColor = System.Drawing.Color.DimGray;
this.label9.ForeColor = System.Drawing.Color.DarkGray;
this.label9.Location = new System.Drawing.Point(8, 215);
this.label9.Location = new System.Drawing.Point(8, 237);
this.label9.Name = "label9";
this.label9.Size = new System.Drawing.Size(258, 33);
this.label9.TabIndex = 26;
@@ -99,7 +102,7 @@
//
this.label8.BackColor = System.Drawing.Color.DimGray;
this.label8.ForeColor = System.Drawing.Color.DarkGray;
this.label8.Location = new System.Drawing.Point(8, 88);
this.label8.Location = new System.Drawing.Point(8, 110);
this.label8.Name = "label8";
this.label8.Size = new System.Drawing.Size(258, 33);
this.label8.TabIndex = 25;
@@ -116,11 +119,66 @@
this.label5.TabIndex = 21;
this.label5.Text = "Minimum angle";
//
// label1
//
this.label1.BackColor = System.Drawing.Color.DimGray;
this.label1.ForeColor = System.Drawing.Color.DarkGray;
this.label1.Location = new System.Drawing.Point(8, 304);
this.label1.Name = "label1";
this.label1.Size = new System.Drawing.Size(258, 33);
this.label1.TabIndex = 26;
this.label1.Text = "Use Sweepline algorithm for triangulation instead of default Divide && Conquer.";
//
// label2
//
this.label2.AutoSize = true;
this.label2.ForeColor = System.Drawing.Color.White;
this.label2.Location = new System.Drawing.Point(8, 64);
this.label2.Name = "label2";
this.label2.Size = new System.Drawing.Size(88, 13);
this.label2.TabIndex = 21;
this.label2.Text = "Maximum angle";
//
// lbMaxAngle
//
this.lbMaxAngle.AutoSize = true;
this.lbMaxAngle.ForeColor = System.Drawing.Color.White;
this.lbMaxAngle.Location = new System.Drawing.Point(227, 65);
this.lbMaxAngle.Name = "lbMaxAngle";
this.lbMaxAngle.Size = new System.Drawing.Size(25, 13);
this.lbMaxAngle.TabIndex = 22;
this.lbMaxAngle.Text = "180";
//
// cbSweepline
//
this.cbSweepline.BackColor = System.Drawing.Color.DimGray;
this.cbSweepline.Checked = false;
this.cbSweepline.Location = new System.Drawing.Point(11, 278);
this.cbSweepline.Name = "cbSweepline";
this.cbSweepline.Size = new System.Drawing.Size(181, 23);
this.cbSweepline.TabIndex = 27;
this.cbSweepline.Text = "Use Sweepline algorithm";
this.cbSweepline.UseVisualStyleBackColor = false;
//
// slMaxAngle
//
this.slMaxAngle.BackColor = System.Drawing.Color.Transparent;
this.slMaxAngle.CriticalPercent = ((uint)(89u));
this.slMaxAngle.Location = new System.Drawing.Point(102, 61);
this.slMaxAngle.Maximum = 100;
this.slMaxAngle.Minimum = 0;
this.slMaxAngle.Name = "slMaxAngle";
this.slMaxAngle.Size = new System.Drawing.Size(119, 18);
this.slMaxAngle.TabIndex = 18;
this.slMaxAngle.Text = "darkSlider1";
this.slMaxAngle.Value = 0;
this.slMaxAngle.ValueChanging += new System.EventHandler(this.slMaxAngle_ValueChanging);
//
// slMaxArea
//
this.slMaxArea.BackColor = System.Drawing.Color.Transparent;
this.slMaxArea.CriticalPercent = ((uint)(0u));
this.slMaxArea.Location = new System.Drawing.Point(102, 61);
this.slMaxArea.Location = new System.Drawing.Point(102, 83);
this.slMaxArea.Maximum = 100;
this.slMaxArea.Minimum = 0;
this.slMaxArea.Name = "slMaxArea";
@@ -148,7 +206,7 @@
//
this.cbConformDel.BackColor = System.Drawing.Color.DimGray;
this.cbConformDel.Checked = false;
this.cbConformDel.Location = new System.Drawing.Point(11, 131);
this.cbConformDel.Location = new System.Drawing.Point(11, 153);
this.cbConformDel.Name = "cbConformDel";
this.cbConformDel.Size = new System.Drawing.Size(142, 17);
this.cbConformDel.TabIndex = 16;
@@ -159,7 +217,7 @@
//
this.cbConvex.BackColor = System.Drawing.Color.DimGray;
this.cbConvex.Checked = false;
this.cbConvex.Location = new System.Drawing.Point(11, 195);
this.cbConvex.Location = new System.Drawing.Point(11, 217);
this.cbConvex.Name = "cbConvex";
this.cbConvex.Size = new System.Drawing.Size(115, 17);
this.cbConvex.TabIndex = 15;
@@ -177,44 +235,26 @@
this.cbQuality.Text = "Quality mesh";
this.cbQuality.UseVisualStyleBackColor = false;
//
// cbSweepline
//
this.cbSweepline.BackColor = System.Drawing.Color.DimGray;
this.cbSweepline.Checked = false;
this.cbSweepline.Location = new System.Drawing.Point(11, 260);
this.cbSweepline.Name = "cbSweepline";
this.cbSweepline.Size = new System.Drawing.Size(187, 17);
this.cbSweepline.TabIndex = 15;
this.cbSweepline.Text = "Use Sweepline algorithm";
this.cbSweepline.UseVisualStyleBackColor = false;
//
// label1
//
this.label1.BackColor = System.Drawing.Color.DimGray;
this.label1.ForeColor = System.Drawing.Color.DarkGray;
this.label1.Location = new System.Drawing.Point(8, 280);
this.label1.Name = "label1";
this.label1.Size = new System.Drawing.Size(258, 33);
this.label1.TabIndex = 26;
this.label1.Text = "Use Sweepline algorithm for triangulation instead of default Divide && Conquer.";
//
// MeshControlView
//
this.AutoScaleDimensions = new System.Drawing.SizeF(6F, 13F);
this.AutoScaleMode = System.Windows.Forms.AutoScaleMode.Font;
this.BackColor = System.Drawing.Color.DimGray;
this.Controls.Add(this.cbSweepline);
this.Controls.Add(this.lbMaxArea);
this.Controls.Add(this.label6);
this.Controls.Add(this.lbMaxAngle);
this.Controls.Add(this.lbMinAngle);
this.Controls.Add(this.label23);
this.Controls.Add(this.label1);
this.Controls.Add(this.label9);
this.Controls.Add(this.label8);
this.Controls.Add(this.label2);
this.Controls.Add(this.label5);
this.Controls.Add(this.slMaxAngle);
this.Controls.Add(this.slMaxArea);
this.Controls.Add(this.slMinAngle);
this.Controls.Add(this.cbConformDel);
this.Controls.Add(this.cbSweepline);
this.Controls.Add(this.cbConvex);
this.Controls.Add(this.cbQuality);
this.Font = new System.Drawing.Font("Segoe UI", 8.25F, System.Drawing.FontStyle.Regular, System.Drawing.GraphicsUnit.Point, ((byte)(0)));
@@ -242,5 +282,8 @@
private Controls.DarkCheckBox cbQuality;
private Controls.DarkCheckBox cbSweepline;
private System.Windows.Forms.Label label1;
private Controls.DarkSlider slMaxAngle;
private System.Windows.Forms.Label label2;
private System.Windows.Forms.Label lbMaxAngle;
}
}
@@ -43,6 +43,11 @@ namespace MeshExplorer.Views
get { return (slMinAngle.Value * 40) / 100; }
}
public int ParamMaxAngleValue
{
get { return 80 + (100 - slMaxAngle.Value); }
}
public double ParamMaxAreaValue
{
get { return slMaxArea.Value * 0.01; }
@@ -55,6 +60,13 @@ namespace MeshExplorer.Views
lbMinAngle.Text = angle.ToString();
}
private void slMaxAngle_ValueChanging(object sender, EventArgs e)
{
// Between 180 and 80 (step 1)
int angle = 80 + (100 - slMaxAngle.Value);
lbMaxAngle.Text = angle.ToString();
}
private void slMaxArea_ValueChanging(object sender, EventArgs e)
{
// Between 0 and 1 (step 0.01)
+423
View File
@@ -0,0 +1,423 @@
// -----------------------------------------------------------------------
// <copyright file="QuadTree.cs" company="">
// Original code by Frank Dockhorn, http://sourceforge.net/projects/quadtreesim/
// Triangle.NET code by Christian Woltering, http://triangle.codeplex.com/
// </copyright>
// -----------------------------------------------------------------------
namespace TriangleNet.Tools
{
using System.Collections.Generic;
using System.Linq;
using TriangleNet.Geometry;
/// <summary>
/// A Quadtree implementation optimised for triangles.
/// </summary>
public class QuadTree
{
QuadNode root;
internal ITriangle[] triangles;
internal int sizeBound;
internal int maxDepth;
/// <summary>
/// Initializes a new instance of the <see cref="QuadTree" /> class.
/// </summary>
/// <param name="mesh">Mesh containing triangles.</param>
/// <param name="maxDepth">The maximum depth of the tree.</param>
/// <param name="sizeBound">The maximum number of triangles contained in a leaf.</param>
/// <remarks>
/// The quadtree does not track changes of the mesh. If a mesh is refined or
/// changed in any other way, a new quadtree has to be built to make the point
/// location work.
///
/// A node of the tree will be split, if its level if less than the max depth parameter
/// AND the number of triangles in the node is greater than the size bound.
/// </remarks>
public QuadTree(Mesh mesh, int maxDepth = 10, int sizeBound = 10)
{
this.maxDepth = maxDepth;
this.sizeBound = sizeBound;
triangles = mesh.Triangles.ToArray();
int currentDepth = 0;
root = new QuadNode(mesh.Bounds, this, true);
root.CreateSubRegion(++currentDepth);
}
public ITriangle Query(double x, double y)
{
var point = new Point(x, y);
var indices = root.FindTriangles(point);
var result = new List<ITriangle>();
foreach (var i in indices)
{
var tri = this.triangles[i];
if (IsPointInTriangle(point, tri.GetVertex(0), tri.GetVertex(1), tri.GetVertex(2)))
{
result.Add(tri);
}
}
return result.FirstOrDefault();
}
/// <summary>
/// Test, if a given point lies inside a triangle.
/// </summary>
/// <param name="p">Point to locate.</param>
/// <param name="t0">Corner point of triangle.</param>
/// <param name="t1">Corner point of triangle.</param>
/// <param name="t2">Corner point of triangle.</param>
/// <returns>True, if point is inside or on the edge of this triangle.</returns>
internal static bool IsPointInTriangle(Point p, Point t0, Point t1, Point t2)
{
// TODO: no need to create new Point instances here
Point d0 = new Point(t1.X - t0.X, t1.Y - t0.Y);
Point d1 = new Point(t2.X - t0.X, t2.Y - t0.Y);
Point d2 = new Point(p.X - t0.X, p.Y - t0.Y);
// crossproduct of (0, 0, 1) and d0
Point c0 = new Point(-d0.Y, d0.X);
// crossproduct of (0, 0, 1) and d1
Point c1 = new Point(-d1.Y, d1.X);
// Linear combination d2 = s * d0 + v * d1.
//
// Multiply both sides of the equation with c0 and c1
// and solve for s and v respectively
//
// s = d2 * c1 / d0 * c1
// v = d2 * c0 / d1 * c0
double s = DotProduct(d2, c1) / DotProduct(d0, c1);
double v = DotProduct(d2, c0) / DotProduct(d1, c0);
if (s >= 0 && v >= 0 && ((s + v) <= 1))
{
// Point is inside or on the edge of this triangle.
return true;
}
return false;
}
internal static double DotProduct(Point p, Point q)
{
return p.X * q.X + p.Y * q.Y;
}
}
#region QuadNode class
/// <summary>
/// A node of the quadtree.
/// </summary>
class QuadNode
{
const int SW = 0;
const int SE = 1;
const int NW = 2;
const int NE = 3;
const double EPS = 1e-6;
static readonly byte[] BITVECTOR = { 0x1, 0x2, 0x4, 0x8 };
BoundingBox bounds;
Point pivot;
QuadTree tree;
QuadNode[] regions;
List<int> triangles;
byte bitRegions;
public QuadNode(BoundingBox box, QuadTree tree)
: this(box, tree, false)
{
}
public QuadNode(BoundingBox box, QuadTree tree, bool init)
{
this.tree = tree;
this.bounds = new BoundingBox(box.Xmin, box.Ymin, box.Xmax, box.Ymax);
this.pivot = new Point((box.Xmin + box.Xmax) / 2, (box.Ymin + box.Ymax) / 2);
this.bitRegions = 0;
this.regions = new QuadNode[4];
this.triangles = new List<int>();
if (init)
{
// Allocate memory upfront
triangles.Capacity = tree.triangles.Length;
foreach (var tri in tree.triangles)
{
triangles.Add(tri.ID);
}
}
}
public List<int> FindTriangles(Point searchPoint)
{
int region = FindRegion(searchPoint);
if (regions[region] == null)
{
return triangles;
}
return regions[region].FindTriangles(searchPoint);
}
public void CreateSubRegion(int currentDepth)
{
// The four sub regions of the quad tree
// +--------------+
// | nw | ne |
// |------+pivot--|
// | sw | se |
// +--------------+
BoundingBox box;
// 1. region south west
box = new BoundingBox(bounds.Xmin, bounds.Ymin, pivot.X, pivot.Y);
regions[0] = new QuadNode(box, tree);
// 2. region south east
box = new BoundingBox(pivot.X, bounds.Ymin, bounds.Xmax, pivot.Y);
regions[1] = new QuadNode(box, tree);
// 3. region north west
box = new BoundingBox(bounds.Xmin, pivot.Y, pivot.X, bounds.Ymax);
regions[2] = new QuadNode(box, tree);
// 4. region north east
box = new BoundingBox(pivot.X, pivot.Y, bounds.Xmax, bounds.Ymax);
regions[3] = new QuadNode(box, tree);
Point[] triangle = new Point[3];
// Find region for every triangle vertex
foreach (var index in triangles)
{
ITriangle tri = tree.triangles[index];
triangle[0] = tri.GetVertex(0);
triangle[1] = tri.GetVertex(1);
triangle[2] = tri.GetVertex(2);
AddTriangleToRegion(triangle, tri.ID);
}
for (int i = 0; i < 4; i++)
{
if (regions[i].triangles.Count > tree.sizeBound && currentDepth < tree.maxDepth)
{
regions[i].CreateSubRegion(currentDepth + 1);
}
}
}
void AddTriangleToRegion(Point[] triangle, int index)
{
bitRegions = 0;
if (QuadTree.IsPointInTriangle(pivot, triangle[0], triangle[1], triangle[2]))
{
AddToRegion(index, SW);
AddToRegion(index, SE);
AddToRegion(index, NW);
AddToRegion(index, NE);
return;
}
FindTriangleIntersections(triangle, index);
if (bitRegions == 0)
{
// we didn't find any intersection so we add this triangle to a point's region
int region = FindRegion(triangle[0]);
regions[region].triangles.Add(index);
}
}
void FindTriangleIntersections(Point[] triangle, int index)
{
// PLEASE NOTE: Handling of component comparison is tightly associated with the implementation
// of the findRegion() function. That means when the point to be compared equals
// the pivot point the triangle must be put at least into region 2.
// Linear equations are in parametric form.
// m_pivot.dx = triangle[0].dx + t * (triangle[1].dx - triangle[0].dx)
// m_pivot.dy = triangle[0].dy + t * (triangle[1].dy - triangle[0].dy)
int k = 2;
double dx, dy;
// Iterate through all triangle laterals and find bounding box intersections
for (int i = 0; i < 3; k = i++)
{
dx = triangle[i].X - triangle[k].X;
dy = triangle[i].Y - triangle[k].Y;
if (dx != 0.0)
{
FindIntersectionsWithX(dx, dy, triangle, index, k);
}
if (dy != 0.0)
{
FindIntersectionsWithY(dx, dy, triangle, index, k);
}
}
}
void FindIntersectionsWithX(double dx, double dy, Point[] triangle, int index, int k)
{
// find intersection with plane x = m_pivot.dX
double t = (pivot.X - triangle[k].X) / dx;
if (t < (1 + EPS) && t > -EPS)
{
// we have an intersection
double yComponent = triangle[k].Y + t * dy;
if (yComponent < pivot.Y)
{
if (yComponent >= bounds.Ymin)
{
AddToRegion(index, SW);
AddToRegion(index, SE);
}
}
else if (yComponent <= bounds.Ymax)
{
AddToRegion(index, NW);
AddToRegion(index, NE);
}
}
// find intersection with plane x = m_boundingBox[0].dX
t = (bounds.Xmin - triangle[k].X) / dx;
if (t < (1 + EPS) && t > -EPS)
{
// we have an intersection
double yComponent = triangle[k].Y + t * dy;
if (yComponent <= pivot.Y && yComponent >= bounds.Ymin)
{
AddToRegion(index, SW);
}
else if (yComponent >= pivot.Y && yComponent <= bounds.Ymax)
{
AddToRegion(index, NW);
}
}
// find intersection with plane x = m_boundingBox[1].dX
t = (bounds.Xmax - triangle[k].X) / dx;
if (t < (1 + EPS) && t > -EPS)
{
// we have an intersection
double yComponent = triangle[k].Y + t * dy;
if (yComponent <= pivot.Y && yComponent >= bounds.Ymin)
{
AddToRegion(index, SE);
}
else if (yComponent >= pivot.Y && yComponent <= bounds.Ymax)
{
AddToRegion(index, NE);
}
}
}
void FindIntersectionsWithY(double dx, double dy, Point[] triangle, int index, int k)
{
// find intersection with plane y = m_pivot.dY
double t = (pivot.Y - triangle[k].Y) / (dy);
if (t < (1 + EPS) && t > -EPS)
{
// we have an intersection
double xComponent = triangle[k].X + t * (dy);
if (xComponent > pivot.X)
{
if (xComponent <= bounds.Xmax)
{
AddToRegion(index, SE);
AddToRegion(index, NE);
}
}
else if (xComponent >= bounds.Xmin)
{
AddToRegion(index, SW);
AddToRegion(index, NW);
}
}
// find intersection with plane y = m_boundingBox[0].dY
t = (bounds.Ymin - triangle[k].Y) / dy;
if (t < (1 + EPS) && t > -EPS)
{
// we have an intersection
double xComponent = triangle[k].X + t * dx;
if (xComponent <= pivot.X && xComponent >= bounds.Xmin)
{
AddToRegion(index, SW);
}
else if (xComponent >= pivot.X && xComponent <= bounds.Xmax)
{
AddToRegion(index, SE);
}
}
// find intersection with plane y = m_boundingBox[1].dY
t = (bounds.Ymax - triangle[k].Y) / dy;
if (t < (1 + EPS) && t > -EPS)
{
// we have an intersection
double xComponent = triangle[k].X + t * dx;
if (xComponent <= pivot.X && xComponent >= bounds.Xmin)
{
AddToRegion(index, NW);
}
else if (xComponent >= pivot.X && xComponent <= bounds.Xmax)
{
AddToRegion(index, NE);
}
}
}
int FindRegion(Point point)
{
int b = 2;
if (point.Y < pivot.Y)
{
b = 0;
}
if (point.X > pivot.X)
{
b++;
}
return b;
}
void AddToRegion(int index, int region)
{
//if (!(m_bitRegions & BITVECTOR[region]))
if ((bitRegions & BITVECTOR[region]) == 0)
{
regions[region].triangles.Add(index);
bitRegions |= BITVECTOR[region];
}
}
}
#endregion
}
+1
View File
@@ -87,6 +87,7 @@
<Compile Include="Tools\BoundedVoronoi.cs" />
<Compile Include="Tools\CuthillMcKee.cs" />
<Compile Include="Tools\IVoronoi.cs" />
<Compile Include="Tools\QuadTree.cs" />
<Compile Include="Tools\QualityMeasure.cs" />
<Compile Include="Tools\RegionIterator.cs" />
<Compile Include="Tools\Statistic.cs" />