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

Update StandardVoronoi code,

Browse Source 
AdjacencyMatrix column pointers are now 0-based,
Rename QuadTree to TriangleQuadTree,
Minor fixes

git-svn-id: https://triangle.svn.codeplex.com/svn@77031 0e2699bc-83d4-4a8f-98e7-55e24ab8c7a5
This commit is contained in:
SND\wo80_cp
2015-07-09 17:05:07 +00:00
parent f00ac93176
commit 58755fa461
18 changed files with 316 additions and 354 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Triangle.NET/TestApp/FormMain.cs
+10 -1
View File
@@ -370,7 +370,16 @@ namespace MeshExplorer
"Do you want to import the mesh?", MessageBoxButtons.YesNo) == DialogResult.OK)
{
input = null;
mesh = FileProcessor.Import(filename);
try
{
mesh = FileProcessor.Import(filename);
}
catch (Exception e)
{
DarkMessageBox.Show("Import mesh error", e.Message, MessageBoxButtons.OK);
return false;
}
if (mesh != null)
{
Triangle.NET/TestApp/FormTopology.cs
+2 -2
View File
@@ -12,7 +12,7 @@ namespace MeshExplorer
public partial class FormTopology : Form
{
Mesh mesh;
QuadTree tree;
TriangleQuadTree tree;
Otri current;
public FormTopology()
@@ -63,7 +63,7 @@ namespace MeshExplorer
if (tree == null)
{
tree = new QuadTree(mesh, 5, 2);
tree = new TriangleQuadTree(mesh, 5, 2);
}
return tree.Query(p.X, p.Y);
Triangle.NET/TestApp/Generators/RandomPoints.cs
+3 -3
View File
@@ -23,7 +23,7 @@ namespace MeshExplorer.Generators
descriptions[1] = "Width:";
descriptions[2] = "Height:";
ranges[0] = new int[] { 5, 5000 };
ranges[0] = new int[] { 10, 5000 };
ranges[1] = new int[] { 10, 200 };
ranges[2] = new int[] { 10, 200 };
}
@@ -33,9 +33,9 @@ namespace MeshExplorer.Generators
int numPoints = GetParamValueInt(0, param0);
numPoints = (numPoints / 10) * 10;
if (numPoints < 5)
if (numPoints < ranges[0][0])
{
numPoints = 5;
numPoints = ranges[0][0];
}
var input = new Polygon(numPoints);
Triangle.NET/Triangle/Behavior.cs
+1 -1
View File
@@ -69,7 +69,7 @@ namespace TriangleNet
Log.Instance.Warning("Invalid quality option (minimum angle).", "Mesh.Behavior");
}
if ((this.maxAngle != 0.0) && this.maxAngle < 90 || this.maxAngle > 180)
if ((this.maxAngle != 0.0) && (this.maxAngle < 60 || this.maxAngle > 180))
{
this.maxAngle = 0;
this.quality = false;
Triangle.NET/Triangle/Mesh.cs
+1 -1
View File
@@ -1733,7 +1733,7 @@ namespace TriangleNet
// Count the degree of the vertex being deleted.
deltri.Onext(ref countingtri);
edgecount = 1;
while (!deltri.Equal(countingtri))
while (!deltri.Equals(countingtri))
{
edgecount++;
countingtri.Onext();
Triangle.NET/Triangle/Meshing/Algorithm/Dwyer.cs
+1 -1
View File
@@ -893,7 +893,7 @@ namespace TriangleNet.Meshing.Algorithm
// Delete the bounding triangle.
mesh.TriangleDealloc(deadtriangle.tri);
} while (!dissolveedge.Equal(startghost));
} while (!dissolveedge.Equals(startghost));
return hullsize;
}
Triangle.NET/Triangle/Meshing/Algorithm/Incremental.cs
+1 -1
View File
@@ -147,7 +147,7 @@ namespace TriangleNet.Meshing.Algorithm
mesh.dummytri.neighbors[0] = searchedge;
hullsize = -2;
while (!nextedge.Equal(finaledge))
while (!nextedge.Equals(finaledge))
{
hullsize++;
nextedge.Lprev(ref dissolveedge);
Triangle.NET/Triangle/Meshing/Algorithm/SweepLine.cs
+4 -4
View File
@@ -124,7 +124,7 @@ namespace TriangleNet.Meshing.Algorithm
fliptri.Onext(ref farrighttri);
Check4DeadEvent(ref farrighttri, eventheap, ref heapsize);
if (farlefttri.Equal(bottommost))
if (farlefttri.Equals(bottommost))
{
fliptri.Lprev(ref bottommost);
}
@@ -191,7 +191,7 @@ namespace TriangleNet.Meshing.Algorithm
righttri.Lprev();
lefttri.Bond(ref farlefttri);
righttri.Bond(ref farrighttri);
if (!farrightflag && farrighttri.Equal(bottommost))
if (!farrightflag && farrighttri.Equals(bottommost))
{
lefttri.Copy(ref bottommost);
}
@@ -585,7 +585,7 @@ namespace TriangleNet.Meshing.Algorithm
while (!farrightflag && RightOfHyperbola(ref searchtri, searchvertex))
{
searchtri.Onext();
farrightflag = searchtri.Equal(bottommost);
farrightflag = searchtri.Equals(bottommost);
}
farright = farrightflag;
return splayroot;
@@ -735,7 +735,7 @@ namespace TriangleNet.Meshing.Algorithm
// Delete the bounding triangle.
mesh.TriangleDealloc(deadtriangle.tri);
} while (!dissolveedge.Equal(startghost));
} while (!dissolveedge.Equals(startghost));
return hullsize;
}
Triangle.NET/Triangle/Meshing/ConstraintMesher.cs
+3 -3
View File
@@ -317,7 +317,7 @@ namespace TriangleNet.Meshing
hulltri.Oprev(ref nexttri);
}
} while (!hulltri.Equal(starttri));
} while (!hulltri.Equals(starttri));
}
/// <summary>
@@ -446,7 +446,7 @@ namespace TriangleNet.Meshing
testtri.Onext(ref neighbor);
// Stop upon reaching a boundary or the starting triangle.
while ((neighbor.tri.id != Mesh.DUMMY) &&
(!neighbor.Equal(testtri)))
(!neighbor.Equals(testtri)))
{
if (neighbor.IsInfected())
{
@@ -1179,7 +1179,7 @@ namespace TriangleNet.Meshing
nexttri.Copy(ref hulltri);
hulltri.Oprev(ref nexttri);
}
} while (!hulltri.Equal(starttri));
} while (!hulltri.Equals(starttri));
}
#endregion
Triangle.NET/Triangle/NewLocation.cs
+2 -2
View File
@@ -3429,8 +3429,8 @@ namespace TriangleNet
{
for (i = 0; i < numpolys; i++)
{
min = 99999999999999999;
max = -99999999999999999;
min = double.MaxValue;
max = double.MinValue;
// compute the minimum and maximum of the
// third coordinate of the cross product
for (j = 1; j <= 2 * polys[i][0] - 1; j = j + 2)
Triangle.NET/Triangle/Tools/AdjacencyMatrix.cs
+94 -221
View File
@@ -8,49 +8,74 @@
namespace TriangleNet.Tools
{
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
/// <summary>
/// The adjacency matrix of the mesh.
/// </summary>
public class AdjacencyMatrix
{
// Number of nodes in the mesh.
int node_num;
// Number of adjacency entries.
int adj_num;
int nnz;
// Pointers into the actual adjacency structure adj. Information about row k is
// stored in entries adj_row(k) through adj_row(k+1)-1 of adj. Size: node_num + 1
int[] adj_row;
// stored in entries pcol(k) through pcol(k+1)-1 of adj. Size: N + 1
int[] pcol;
// The adjacency structure. For each row, it contains the column indices
// of the nonzero entries. Size: adj_num
int[] adj;
// of the nonzero entries. Size: nnz
int[] irow;
public int[] AdjacencyRow
/// <summary>
/// Gets the number of columns (nodes of the mesh).
/// </summary>
public readonly int N;
/// <summary>
/// Gets the column pointers.
/// </summary>
public int[] ColumnPointers
{
get { return adj_row; }
get { return pcol; }
}
public int[] Adjacency
/// <summary>
/// Gets the row indices.
/// </summary>
public int[] RowIndices
{
get { return adj; }
get { return irow; }
}
public AdjacencyMatrix(Mesh mesh)
{
this.node_num = mesh.vertices.Count;
this.N = mesh.vertices.Count;
// Set up the adj_row adjacency pointer array.
this.adj_row = AdjacencyCount(mesh);
this.adj_num = adj_row[node_num] - 1;
this.pcol = AdjacencyCount(mesh);
this.nnz = pcol[N];
// Set up the adj adjacency array.
this.adj = AdjacencySet(mesh, this.adj_row);
this.irow = AdjacencySet(mesh, this.pcol);
}
public AdjacencyMatrix(int[] pcol, int[] irow)
{
this.N = pcol.Length - 1;
this.nnz = pcol[N];
this.pcol = pcol;
this.irow = irow;
if (pcol[0] != 0)
{
throw new ArgumentException("Expected 0-based indexing.", "pcol");
}
if (irow.Length < nnz)
{
throw new ArgumentException();
}
}
/// <summary>
@@ -67,11 +92,11 @@ namespace TriangleNet.Tools
band_lo = 0;
band_hi = 0;
for (i = 0; i < node_num; i++)
for (i = 0; i < N; i++)
{
for (j = adj_row[i]; j <= adj_row[i + 1] - 1; j++)
for (j = pcol[i]; j < pcol[i + 1]; j++)
{
col = adj[j - 1];
col = irow[j];
band_lo = Math.Max(band_lo, i - col);
band_hi = Math.Max(band_hi, col - i);
}
@@ -94,41 +119,25 @@ namespace TriangleNet.Tools
///
/// Two nodes are "adjacent" if they are both nodes in some triangle.
/// Also, a node is considered to be adjacent to itself.
///
/// Diagram:
///
/// 3
/// s |\
/// i | \
/// d | \
/// e | \ side 1
/// | \
/// 2 | \
/// | \
/// 1-------2
///
/// side 3
/// </remarks>
int[] AdjacencyCount(Mesh mesh)
{
int i;
int node;
int n = N;
int n1, n2, n3;
int tri_id;
int neigh_id;
int tid, nid;
int[] adj_rows = new int[node_num + 1];
int[] pcol = new int[n + 1];
// Set every node to be adjacent to itself.
for (node = 0; node < node_num; node++)
for (int i = 0; i < n; i++)
{
adj_rows[node] = 1;
pcol[i] = 1;
}
// Examine each triangle.
foreach (var tri in mesh.triangles.Values)
{
tri_id = tri.id;
tid = tri.id;
n1 = tri.vertices[0].id;
n2 = tri.vertices[1].id;
@@ -137,47 +146,47 @@ namespace TriangleNet.Tools
// Add edge (1,2) if this is the first occurrence, that is, if
// the edge (1,2) is on a boundary (nid <= 0) or if this triangle
// is the first of the pair in which the edge occurs (tid < nid).
neigh_id = tri.neighbors[2].tri.id;
nid = tri.neighbors[2].tri.id;
if (neigh_id < 0 || tri_id < neigh_id)
if (nid < 0 || tid < nid)
{
adj_rows[n1] += 1;
adj_rows[n2] += 1;
pcol[n1] += 1;
pcol[n2] += 1;
}
// Add edge (2,3).
neigh_id = tri.neighbors[0].tri.id;
nid = tri.neighbors[0].tri.id;
if (neigh_id < 0 || tri_id < neigh_id)
if (nid < 0 || tid < nid)
{
adj_rows[n2] += 1;
adj_rows[n3] += 1;
pcol[n2] += 1;
pcol[n3] += 1;
}
// Add edge (3,1).
neigh_id = tri.neighbors[1].tri.id;
nid = tri.neighbors[1].tri.id;
if (neigh_id < 0 || tri_id < neigh_id)
if (nid < 0 || tid < nid)
{
adj_rows[n3] += 1;
adj_rows[n1] += 1;
pcol[n3] += 1;
pcol[n1] += 1;
}
}
// We used ADJ_COL to count the number of entries in each column.
// We used PCOL to count the number of entries in each column.
// Convert it to pointers into the ADJ array.
for (node = node_num; 1 <= node; node--)
for (int i = n; i > 0; i--)
{
adj_rows[node] = adj_rows[node - 1];
pcol[i] = pcol[i - 1];
}
adj_rows[0] = 1;
for (i = 1; i <= node_num; i++)
pcol[0] = 0;
for (int i = 1; i <= n; i++)
{
adj_rows[i] = adj_rows[i - 1] + adj_rows[i];
pcol[i] = pcol[i - 1] + pcol[i];
}
return adj_rows;
return pcol;
}
/// <summary>
@@ -188,28 +197,25 @@ namespace TriangleNet.Tools
/// for a linear triangle finite element discretization of Poisson's
/// equation in two dimensions.
/// </remarks>
int[] AdjacencySet(Mesh mesh, int[] rows)
int[] AdjacencySet(Mesh mesh, int[] pcol)
{
// Output list, stores the actual adjacency information.
int[] list;
int n = this.N;
int[] col = new int[n];
// Copy of the adjacency rows input.
int[] rowsCopy = new int[node_num];
Array.Copy(rows, rowsCopy, node_num);
Array.Copy(pcol, col, n);
int i, n = rows[node_num] - 1;
int i, nnz = pcol[n];
list = new int[n];
for (i = 0; i < n; i++)
{
list[i] = -1;
}
// Output list, stores the actual adjacency information.
int[] list = new int[nnz];
// Set every node to be adjacent to itself.
for (i = 0; i < node_num; i++)
for (i = 0; i < n; i++)
{
list[rowsCopy[i] - 1] = i;
rowsCopy[i] += 1;
list[col[i]] = i;
col[i] += 1;
}
int n1, n2, n3; // Vertex numbers.
@@ -231,10 +237,8 @@ namespace TriangleNet.Tools
if (nid < 0 || tid < nid)
{
list[rowsCopy[n1] - 1] = n2;
rowsCopy[n1] += 1;
list[rowsCopy[n2] - 1] = n1;
rowsCopy[n2] += 1;
list[col[n1]++] = n2;
list[col[n2]++] = n1;
}
// Add edge (2,3).
@@ -242,10 +246,8 @@ namespace TriangleNet.Tools
if (nid < 0 || tid < nid)
{
list[rowsCopy[n2] - 1] = n3;
rowsCopy[n2] += 1;
list[rowsCopy[n3] - 1] = n2;
rowsCopy[n3] += 1;
list[col[n2]++] = n3;
list[col[n3]++] = n2;
}
// Add edge (3,1).
@@ -253,153 +255,24 @@ namespace TriangleNet.Tools
if (nid < 0 || tid < nid)
{
list[rowsCopy[n1] - 1] = n3;
rowsCopy[n1] += 1;
list[rowsCopy[n3] - 1] = n1;
rowsCopy[n3] += 1;
list[col[n1]++] = n3;
list[col[n3]++] = n1;
}
}
int k1, k2;
// Ascending sort the entries for each node.
for (i = 0; i < node_num; i++)
// Ascending sort the entries for each column.
for (i = 0; i < n; i++)
{
k1 = rows[i];
k2 = rows[i + 1] - 1;
HeapSort(list, k1 - 1, k2 + 1 - k1);
k1 = pcol[i];
k2 = pcol[i + 1];
Array.Sort(list, k1, k2 - k1);
}
return list;
}
#endregion
#region Heap sort
/// <summary>
/// Reorders an array of integers into a descending heap.
/// </summary>
/// <param name="size">the size of the input array.</param>
/// <param name="a">an unsorted array.</param>
/// <remarks>
/// A heap is an array A with the property that, for every index J,
/// A[J] >= A[2*J+1] and A[J] >= A[2*J+2], (as long as the indices
/// 2*J+1 and 2*J+2 are legal).
///
/// Diagram:
///
/// A(0)
/// / \
/// A(1) A(2)
/// / \ / \
/// A(3) A(4) A(5) A(6)
/// / \ / \
/// A(7) A(8) A(9) A(10)
/// </remarks>
private void CreateHeap(int[] a, int offset, int size)
{
int i;
int ifree;
int key;
int m;
// Only nodes (N/2)-1 down to 0 can be "parent" nodes.
for (i = (size / 2) - 1; 0 <= i; i--)
{
// Copy the value out of the parent node.
// Position IFREE is now "open".
key = a[offset + i];
ifree = i;
for (; ; )
{
// Positions 2*IFREE + 1 and 2*IFREE + 2 are the descendants of position
// IFREE. (One or both may not exist because they equal or exceed N.)
m = 2 * ifree + 1;
// Does the first position exist?
if (size <= m)
{
break;
}
else
{
// Does the second position exist?
if (m + 1 < size)
{
// If both positions exist, take the larger of the two values,
// and update M if necessary.
if (a[offset + m] < a[offset + m + 1])
{
m = m + 1;
}
}
// If the large descendant is larger than KEY, move it up,
// and update IFREE, the location of the free position, and
// consider the descendants of THIS position.
if (key < a[offset + m])
{
a[offset + ifree] = a[offset + m];
ifree = m;
}
else
{
break;
}
}
}
// When you have stopped shifting items up, return the item you
// pulled out back to the heap.
a[offset + ifree] = key;
}
return;
}
/// <summary>
/// ascending sorts an array of integers using heap sort.
/// </summary>
/// <param name="size">Number of entries in the array.</param>
/// <param name="a">Array to be sorted;</param>
private void HeapSort(int[] a, int offset, int size)
{
int n1;
int temp;
if (size <= 1)
{
return;
}
// 1: Put A into descending heap form.
CreateHeap(a, offset, size);
// 2: Sort A.
// The largest object in the heap is in A[0].
// Move it to position A[N-1].
temp = a[offset];
a[offset] = a[offset + size - 1];
a[offset + size - 1] = temp;
// Consider the diminished heap of size N1.
for (n1 = size - 1; 2 <= n1; n1--)
{
// Restore the heap structure of the initial N1 entries of A.
CreateHeap(a, offset, n1);
// Take the largest object from A[0] and move it to A[N1-1].
temp = a[offset];
a[offset] = a[offset + n1 - 1];
a[offset + n1 - 1] = temp;
}
return;
}
#endregion
}
}
Triangle.NET/Triangle/Tools/CuthillMcKee.cs
+119 -89
View File
@@ -8,10 +8,6 @@
namespace TriangleNet.Tools
{
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using TriangleNet.Logging;
/// <summary>
/// Applies the Cuthill and McKee renumbering algorithm to reduce the bandwidth of
@@ -19,9 +15,6 @@ namespace TriangleNet.Tools
/// </summary>
public class CuthillMcKee
{
// Number of nodes in the mesh.
int node_num;
// The adjacency matrix of the mesh.
AdjacencyMatrix matrix;
@@ -32,24 +25,36 @@ namespace TriangleNet.Tools
/// <returns>Permutation vector.</returns>
public int[] Renumber(Mesh mesh)
{
int bandwidth1, bandwidth2;
this.node_num = mesh.vertices.Count;
// Algorithm needs linear numbering of the nodes.
mesh.Renumber(NodeNumbering.Linear);
// Set up the adj_row adjacency pointer array.
matrix = new AdjacencyMatrix(mesh);
return Renumber(new AdjacencyMatrix(mesh));
}
bandwidth1 = matrix.Bandwidth();
/// <summary>
/// Gets the permutation vector for the Reverse Cuthill-McKee numbering.
/// </summary>
/// <param name="mesh">The mesh.</param>
/// <returns>Permutation vector.</returns>
public int[] Renumber(AdjacencyMatrix matrix)
{
this.matrix = matrix;
int bandwidth1 = matrix.Bandwidth();
var pcol = matrix.ColumnPointers;
// Adjust column pointers (1-based indexing).
Shift(pcol, true);
// TODO: Make RCM work with 0-based matrix.
// Compute the RCM permutation.
int[] perm = GenerateRcm();
int[] perm_inv = PermInverse(node_num, perm);
int[] perm_inv = PermInverse(perm);
bandwidth2 = PermBandwidth(perm, perm_inv);
int bandwidth2 = PermBandwidth(perm, perm_inv);
if (Log.Verbose)
{
@@ -57,42 +62,12 @@ namespace TriangleNet.Tools
bandwidth1, bandwidth2));
}
// Adjust column pointers (0-based indexing).
Shift(pcol, false);
return perm_inv;
}
/// <summary>
/// Computes the bandwidth of a permuted adjacency matrix.
/// </summary>
/// <param name="perm">The permutation.</param>
/// <param name="perm_inv">The inverse permutation.</param>
/// <returns>Bandwidth of the permuted adjacency matrix.</returns>
/// <remarks>
/// The matrix is defined by the adjacency information and a permutation.
/// The routine also computes the bandwidth and the size of the envelope.
/// </remarks>
int PermBandwidth(int[] perm, int[] perm_inv)
{
int[] adj_row = matrix.AdjacencyRow;
int[] adj = matrix.Adjacency;
int col, i, j;
int band_lo = 0;
int band_hi = 0;
for (i = 0; i < node_num; i++)
{
for (j = adj_row[perm[i]]; j <= adj_row[perm[i] + 1] - 1; j++)
{
col = perm_inv[adj[j - 1]];
band_lo = Math.Max(band_lo, i - col);
band_hi = Math.Max(band_hi, col - i);
}
}
return band_lo + 1 + band_hi;
}
#region RCM
/// <summary>
@@ -105,7 +80,10 @@ namespace TriangleNet.Tools
/// </remarks>
int[] GenerateRcm()
{
int[] perm = new int[node_num];
// Number of nodes in the mesh.
int n = matrix.N;
int[] perm = new int[n];
int i, num, root;
int iccsze = 0;
@@ -113,19 +91,19 @@ namespace TriangleNet.Tools
/// Index vector for a level structure. The level structure is stored in the
/// currently unused spaces in the permutation vector PERM.
int[] level_row = new int[node_num + 1];
int[] level_row = new int[n + 1];
/// Marks variables that have been numbered.
int[] mask = new int[node_num];
int[] mask = new int[n];
for (i = 0; i < node_num; i++)
for (i = 0; i < n; i++)
{
mask[i] = 1;
}
num = 1;
for (i = 0; i < node_num; i++)
for (i = 0; i < n; i++)
{
// For each masked connected component...
if (mask[i] != 0)
@@ -142,7 +120,7 @@ namespace TriangleNet.Tools
num += iccsze;
// We can stop once every node is in one of the connected components.
if (node_num < num)
if (n < num)
{
return perm;
}
@@ -179,8 +157,8 @@ namespace TriangleNet.Tools
/// </remarks>
void Rcm(int root, int[] mask, int[] perm, int offset, ref int iccsze)
{
int[] adj_row = matrix.AdjacencyRow;
int[] adj = matrix.Adjacency;
int[] pcol = matrix.ColumnPointers;
int[] irow = matrix.RowIndices;
int fnbr;
int i, j, k, l;
@@ -188,9 +166,12 @@ namespace TriangleNet.Tools
int lbegin, lnbr, lperm, lvlend;
int nbr, node;
// Number of nodes in the mesh.
int n = matrix.N;
/// Workspace, int DEG[NODE_NUM], a temporary vector used to hold
/// the degree of the nodes in the section graph specified by mask and root.
int[] deg = new int[node_num];
int[] deg = new int[n];
// Find the degrees of the nodes in the component specified by MASK and ROOT.
Degree(root, mask, deg, ref iccsze, perm, offset);
@@ -216,8 +197,8 @@ namespace TriangleNet.Tools
{
// For each node in the current level...
node = perm[offset + i - 1];
jstrt = adj_row[node];
jstop = adj_row[node + 1] - 1;
jstrt = pcol[node];
jstop = pcol[node + 1] - 1;
// Find the unnumbered neighbors of NODE.
@@ -227,7 +208,7 @@ namespace TriangleNet.Tools
for (j = jstrt; j <= jstop; j++)
{
nbr = adj[j - 1];
nbr = irow[j - 1];
if (mask[nbr] != 0)
{
@@ -330,11 +311,11 @@ namespace TriangleNet.Tools
/// ACM Transactions on Mathematical Software,
/// Volume 2, pages 378-387, 1976.
/// </remarks>
void FindRoot(ref int root, int[] mask, ref int level_num, int[] level_row,
void FindRoot(ref int root, int[] mask, ref int level_num, int[] level_row,
int[] level, int offset)
{
int[] adj_row = matrix.AdjacencyRow;
int[] adj = matrix.Adjacency;
int[] pcol = matrix.ColumnPointers;
int[] irow = matrix.RowIndices;
int iccsze;
int j, jstrt;
@@ -376,12 +357,12 @@ namespace TriangleNet.Tools
{
node = level[offset + j - 1];
ndeg = 0;
kstrt = adj_row[node - 1];
kstop = adj_row[node] - 1;
kstrt = pcol[node - 1];
kstop = pcol[node] - 1;
for (k = kstrt; k <= kstop; k++)
{
nghbor = adj[k - 1];
nghbor = irow[k - 1];
if (mask[nghbor] > 0)
{
ndeg += 1;
@@ -443,11 +424,11 @@ namespace TriangleNet.Tools
/// Computer Solution of Large Sparse Positive Definite Systems,
/// Prentice Hall, 1981.
/// </remarks>
void GetLevelSet(ref int root, int[] mask, ref int level_num, int[] level_row,
void GetLevelSet(ref int root, int[] mask, ref int level_num, int[] level_row,
int[] level, int offset)
{
int[] adj_row = matrix.AdjacencyRow;
int[] adj = matrix.Adjacency;
int[] pcol = matrix.ColumnPointers;
int[] irow = matrix.RowIndices;
int i, iccsze;
int j, jstop, jstrt;
@@ -475,12 +456,12 @@ namespace TriangleNet.Tools
for (i = lbegin; i <= lvlend; i++)
{
node = level[offset + i - 1];
jstrt = adj_row[node];
jstop = adj_row[node + 1] - 1;
jstrt = pcol[node];
jstop = pcol[node + 1] - 1;
for (j = jstrt; j <= jstop; j++)
{
nbr = adj[j - 1];
nbr = irow[j - 1];
if (mask[nbr] != 0)
{
@@ -533,18 +514,18 @@ namespace TriangleNet.Tools
/// </remarks>
void Degree(int root, int[] mask, int[] deg, ref int iccsze, int[] ls, int offset)
{
int[] adj_row = matrix.AdjacencyRow;
int[] adj = matrix.Adjacency;
int[] pcol = matrix.ColumnPointers;
int[] irow = matrix.RowIndices;
int i, ideg;
int j, jstop, jstrt;
int lbegin, lvlend;
int lvsize = 1;
int nghbr, node;
int nbr, node;
// The sign of ADJ_ROW(I) is used to indicate if node I has been considered.
ls[offset] = root;
adj_row[root] = -adj_row[root];
pcol[root] = -pcol[root];
lvlend = 0;
iccsze = 1;
@@ -561,23 +542,23 @@ namespace TriangleNet.Tools
for (i = lbegin; i <= lvlend; i++)
{
node = ls[offset + i - 1];
jstrt = -adj_row[node];
jstop = Math.Abs(adj_row[node + 1]) - 1;
jstrt = -pcol[node];
jstop = Math.Abs(pcol[node + 1]) - 1;
ideg = 0;
for (j = jstrt; j <= jstop; j++)
{
nghbr = adj[j - 1];
nbr = irow[j - 1];
if (mask[nghbr] != 0) // EDIT: [nbr - 1]
if (mask[nbr] != 0) // EDIT: [nbr - 1]
{
ideg = ideg + 1;
if (0 <= adj_row[nghbr]) // EDIT: [nbr - 1]
if (0 <= pcol[nbr]) // EDIT: [nbr - 1]
{
adj_row[nghbr] = -adj_row[nghbr]; // EDIT: [nbr - 1]
pcol[nbr] = -pcol[nbr]; // EDIT: [nbr - 1]
iccsze = iccsze + 1;
ls[offset + iccsze - 1] = nghbr;
ls[offset + iccsze - 1] = nbr;
}
}
}
@@ -592,7 +573,7 @@ namespace TriangleNet.Tools
for (i = 0; i < iccsze; i++)
{
node = ls[offset + i];
adj_row[node] = -adj_row[node];
pcol[node] = -pcol[node];
}
return;
@@ -602,19 +583,54 @@ namespace TriangleNet.Tools
#region Tools
/// <summary>
/// Computes the bandwidth of a permuted adjacency matrix.
/// </summary>
/// <param name="perm">The permutation.</param>
/// <param name="perm_inv">The inverse permutation.</param>
/// <returns>Bandwidth of the permuted adjacency matrix.</returns>
/// <remarks>
/// The matrix is defined by the adjacency information and a permutation.
/// The routine also computes the bandwidth and the size of the envelope.
/// </remarks>
int PermBandwidth(int[] perm, int[] perm_inv)
{
int[] pcol = matrix.ColumnPointers;
int[] irow = matrix.RowIndices;
int col, i, j;
int band_lo = 0;
int band_hi = 0;
int n = matrix.N;
for (i = 0; i < n; i++)
{
for (j = pcol[perm[i]]; j < pcol[perm[i] + 1]; j++)
{
col = perm_inv[irow[j - 1]];
band_lo = Math.Max(band_lo, i - col);
band_hi = Math.Max(band_hi, col - i);
}
}
return band_lo + 1 + band_hi;
}
/// <summary>
/// Produces the inverse of a given permutation.
/// </summary>
/// <param name="n">Number of items permuted.</param>
/// <param name="perm">PERM[N], a permutation.</param>
/// <returns>The inverse permutation.</returns>
int[] PermInverse(int n, int[] perm)
int[] PermInverse(int[] perm)
{
int[] perm_inv = new int[node_num];
int n = matrix.N;
int i;
int[] perm_inv = new int[n];
for (i = 0; i < n; i++)
for (int i = 0; i < n; i++)
{
perm_inv[perm[i]] = i;
}
@@ -650,6 +666,20 @@ namespace TriangleNet.Tools
return;
}
void Shift(int[] a, bool up)
{
int length = a.Length;
if (up)
{
for (int i = 0; i < length; a[i]++, i++) ;
}
else
{
for (int i = 0; i < length; a[i]--, i++) ;
}
}
#endregion
}
}
Triangle.NET/Triangle/Tools/QuadTree.cs → Triangle.NET/Triangle/Tools/TriangleQuadTree.cs
+9 -9
View File
@@ -1,6 +1,6 @@
// -----------------------------------------------------------------------
// <copyright file="QuadTree.cs" company="">
// Original code by Frank Dockhorn, http://sourceforge.net/projects/quadtreesim/
// <copyright file="TriangleQuadTree.cs" company="">
// Original code by Frank Dockhorn, [not available anymore: http://sourceforge.net/projects/quadtreesim/]
// Triangle.NET code by Christian Woltering, http://triangle.codeplex.com/
// </copyright>
// -----------------------------------------------------------------------
@@ -14,7 +14,7 @@ namespace TriangleNet.Tools
/// <summary>
/// A Quadtree implementation optimized for triangles.
/// </summary>
public class QuadTree
public class TriangleQuadTree
{
QuadNode root;
@@ -24,7 +24,7 @@ namespace TriangleNet.Tools
internal int maxDepth;
/// <summary>
/// Initializes a new instance of the <see cref="QuadTree" /> class.
/// Initializes a new instance of the <see cref="TriangleQuadTree" /> class.
/// </summary>
/// <param name="mesh">Mesh containing triangles.</param>
/// <param name="maxDepth">The maximum depth of the tree.</param>
@@ -37,7 +37,7 @@ namespace TriangleNet.Tools
/// 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)
public TriangleQuadTree(Mesh mesh, int maxDepth = 10, int sizeBound = 10)
{
this.maxDepth = maxDepth;
this.sizeBound = sizeBound;
@@ -130,18 +130,18 @@ namespace TriangleNet.Tools
Rectangle bounds;
Point pivot;
QuadTree tree;
TriangleQuadTree tree;
QuadNode[] regions;
List<int> triangles;
byte bitRegions;
public QuadNode(Rectangle box, QuadTree tree)
public QuadNode(Rectangle box, TriangleQuadTree tree)
: this(box, tree, false)
{
}
public QuadNode(Rectangle box, QuadTree tree, bool init)
public QuadNode(Rectangle box, TriangleQuadTree tree, bool init)
{
this.tree = tree;
@@ -232,7 +232,7 @@ namespace TriangleNet.Tools
void AddTriangleToRegion(Point[] triangle, int index)
{
bitRegions = 0;
if (QuadTree.IsPointInTriangle(pivot, triangle[0], triangle[1], triangle[2]))
if (TriangleQuadTree.IsPointInTriangle(pivot, triangle[0], triangle[1], triangle[2]))
{
AddToRegion(index, SW);
AddToRegion(index, SE);
Triangle.NET/Triangle/Topology/Otri.cs
+8 -8
View File
@@ -349,6 +349,14 @@ namespace TriangleNet.Topology
ot.orient = orient;
}
/// <summary>
/// Test for equality of oriented triangles.
/// </summary>
public bool Equals(Otri ot)
{
return ((tri == ot.tri) && (orient == ot.orient));
}
#endregion
#region Otri primitives (internal)
@@ -402,14 +410,6 @@ namespace TriangleNet.Topology
tri.neighbors[orient].orient = 0;
}
/// <summary>
/// Test for equality of oriented triangles.
/// </summary>
internal bool Equal(Otri ot)
{
return ((tri == ot.tri) && (orient == ot.orient));
}
/// <summary>
/// Infect a triangle with the virus.
/// </summary>
Triangle.NET/Triangle/Triangle.csproj
+1 -1
View File
@@ -102,7 +102,7 @@
<Compile Include="Voronoi\Legacy\BoundedVoronoiLegacy.cs" />
<Compile Include="Tools\CuthillMcKee.cs" />
<Compile Include="Voronoi\Legacy\IVoronoi.cs" />
<Compile Include="Tools\QuadTree.cs" />
<Compile Include="Tools\TriangleQuadTree.cs" />
<Compile Include="Tools\QualityMeasure.cs" />
<Compile Include="Meshing\Iterators\RegionIterator.cs" />
<Compile Include="Tools\Statistic.cs" />
Triangle.NET/Triangle/Voronoi/Legacy/BoundedVoronoiLegacy.cs
+3 -3
View File
@@ -372,7 +372,7 @@ namespace TriangleNet.Voronoi.Legacy
// Call f_next the next triangle counterclockwise around x
f_next.Onext();
}
while (!f.Equal(f_init));
while (!f.Equals(f_init));
// Output: Bounded Voronoi cell of x in counterclockwise order.
region.Add(vpoints);
@@ -416,7 +416,7 @@ namespace TriangleNet.Voronoi.Legacy
if (f_prev.tri.id != Mesh.DUMMY)
{
// Go clockwise until we reach the border (or the initial triangle)
while (f_prev.tri.id != Mesh.DUMMY && !f_prev.Equal(f_init))
while (f_prev.tri.id != Mesh.DUMMY && !f_prev.Equals(f_init))
{
f_prev.Copy(ref f);
f_prev.Oprev();
@@ -579,7 +579,7 @@ namespace TriangleNet.Voronoi.Legacy
// Call f_next the next triangle counterclockwise around x
f_next.Onext();
}
while (!f.Equal(f_init));
while (!f.Equals(f_init));
// Output: Bounded Voronoi cell of x in counterclockwise order.
region.Add(vpoints);
Triangle.NET/Triangle/Voronoi/Legacy/SimpleVoronoi.cs
+1 -1
View File
@@ -177,7 +177,7 @@ namespace TriangleNet.Voronoi.Legacy
region.AddNeighbor(f.tri.id, regions[f.Apex().id]);
if (f_next.Equal(f_init))
if (f_next.Equals(f_init))
{
// Voronoi cell is complete (bounded case).
region.Add(vpoints);
Triangle.NET/Triangle/Voronoi/StandardVoronoi.cs
+53 -3
View File
@@ -6,18 +6,68 @@
namespace TriangleNet.Voronoi
{
using System.Collections.Generic;
using TriangleNet.Geometry;
using TriangleNet.Tools;
using TriangleNet.Topology.DCEL;
public class StandardVoronoi : VoronoiBase
{
public StandardVoronoi(Mesh mesh)
: base(mesh, true)
: this(mesh, mesh.bounds)
{
}
private void Intersect(Rectangle box)
public StandardVoronoi(Mesh mesh, Rectangle box)
: base(mesh, true)
{
// TODO: compute edge intersections with bounding box.
// We assume the box to be at least as large as the mesh.
box.Expand(mesh.bounds);
// We explicitly told the base constructor to call the Generate method, so
// at this point the basic Voronoi diagram is already created.
PostProcess(box);
}
/// <summary>
/// Compute edge intersections with bounding box.
/// </summary>
private void PostProcess(Rectangle box)
{
var infEdges = new List<HalfEdge>();
// TODO: save the half-infinite boundary edge in base class
// so we don't have to process the complete list here.
foreach (var edge in base.edges)
{
if (edge.next == null)
{
infEdges.Add(edge);
}
}
foreach (var edge in infEdges)
{
// The vertices of the infinite edge.
var v1 = (Point)edge.origin;
var v2 = (Point)edge.twin.origin;
if (box.Contains(v1) || box.Contains(v2))
{
// Move infinite vertex v2 onto the box boundary.
IntersectionHelper.BoxRayIntersection(box, v1, v2, ref v2);
}
else
{
// There is actually no easy way to handle the second case. The two edges
// leaving v1, pointing towards the mesh, don't have to intersect the box
// (the could join with edges of other cells outside the box).
// A general intersection algorithm (DCEL <-> Rectangle) is needed, which
// computes intersections with all edges and discards objects outside the
// box.
}
}
}
}
}