From 58755fa461bac46b35430c2314bcda88c0b9fc3c Mon Sep 17 00:00:00 2001 From: "SND\\wo80_cp" Date: Thu, 9 Jul 2015 17:05:07 +0000 Subject: [PATCH] Update StandardVoronoi code, 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 --- Triangle.NET/TestApp/FormMain.cs | 11 +- Triangle.NET/TestApp/FormTopology.cs | 4 +- .../TestApp/Generators/RandomPoints.cs | 6 +- Triangle.NET/Triangle/Behavior.cs | 2 +- Triangle.NET/Triangle/Mesh.cs | 2 +- .../Triangle/Meshing/Algorithm/Dwyer.cs | 2 +- .../Triangle/Meshing/Algorithm/Incremental.cs | 2 +- .../Triangle/Meshing/Algorithm/SweepLine.cs | 8 +- .../Triangle/Meshing/ConstraintMesher.cs | 6 +- Triangle.NET/Triangle/NewLocation.cs | 4 +- .../Triangle/Tools/AdjacencyMatrix.cs | 315 ++++++------------ Triangle.NET/Triangle/Tools/CuthillMcKee.cs | 208 +++++++----- .../{QuadTree.cs => TriangleQuadTree.cs} | 18 +- Triangle.NET/Triangle/Topology/Otri.cs | 16 +- Triangle.NET/Triangle/Triangle.csproj | 2 +- .../Voronoi/Legacy/BoundedVoronoiLegacy.cs | 6 +- .../Triangle/Voronoi/Legacy/SimpleVoronoi.cs | 2 +- .../Triangle/Voronoi/StandardVoronoi.cs | 56 +++- 18 files changed, 316 insertions(+), 354 deletions(-) rename Triangle.NET/Triangle/Tools/{QuadTree.cs => TriangleQuadTree.cs} (93%) diff --git a/Triangle.NET/TestApp/FormMain.cs b/Triangle.NET/TestApp/FormMain.cs index 3a3509c..e2f5bb3 100644 --- a/Triangle.NET/TestApp/FormMain.cs +++ b/Triangle.NET/TestApp/FormMain.cs @@ -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) { diff --git a/Triangle.NET/TestApp/FormTopology.cs b/Triangle.NET/TestApp/FormTopology.cs index 8e235c2..c2889be 100644 --- a/Triangle.NET/TestApp/FormTopology.cs +++ b/Triangle.NET/TestApp/FormTopology.cs @@ -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); diff --git a/Triangle.NET/TestApp/Generators/RandomPoints.cs b/Triangle.NET/TestApp/Generators/RandomPoints.cs index ded256d..f0441f5 100644 --- a/Triangle.NET/TestApp/Generators/RandomPoints.cs +++ b/Triangle.NET/TestApp/Generators/RandomPoints.cs @@ -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); diff --git a/Triangle.NET/Triangle/Behavior.cs b/Triangle.NET/Triangle/Behavior.cs index b0d9d91..50c4b89 100644 --- a/Triangle.NET/Triangle/Behavior.cs +++ b/Triangle.NET/Triangle/Behavior.cs @@ -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; diff --git a/Triangle.NET/Triangle/Mesh.cs b/Triangle.NET/Triangle/Mesh.cs index 021208b..9b3b21f 100644 --- a/Triangle.NET/Triangle/Mesh.cs +++ b/Triangle.NET/Triangle/Mesh.cs @@ -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(); diff --git a/Triangle.NET/Triangle/Meshing/Algorithm/Dwyer.cs b/Triangle.NET/Triangle/Meshing/Algorithm/Dwyer.cs index 1a0da2a..91905ea 100644 --- a/Triangle.NET/Triangle/Meshing/Algorithm/Dwyer.cs +++ b/Triangle.NET/Triangle/Meshing/Algorithm/Dwyer.cs @@ -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; } diff --git a/Triangle.NET/Triangle/Meshing/Algorithm/Incremental.cs b/Triangle.NET/Triangle/Meshing/Algorithm/Incremental.cs index d03527f..7fd40f9 100644 --- a/Triangle.NET/Triangle/Meshing/Algorithm/Incremental.cs +++ b/Triangle.NET/Triangle/Meshing/Algorithm/Incremental.cs @@ -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); diff --git a/Triangle.NET/Triangle/Meshing/Algorithm/SweepLine.cs b/Triangle.NET/Triangle/Meshing/Algorithm/SweepLine.cs index 1ff2c2f..97175fc 100644 --- a/Triangle.NET/Triangle/Meshing/Algorithm/SweepLine.cs +++ b/Triangle.NET/Triangle/Meshing/Algorithm/SweepLine.cs @@ -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; } diff --git a/Triangle.NET/Triangle/Meshing/ConstraintMesher.cs b/Triangle.NET/Triangle/Meshing/ConstraintMesher.cs index 9d2dc92..49a54aa 100644 --- a/Triangle.NET/Triangle/Meshing/ConstraintMesher.cs +++ b/Triangle.NET/Triangle/Meshing/ConstraintMesher.cs @@ -317,7 +317,7 @@ namespace TriangleNet.Meshing hulltri.Oprev(ref nexttri); } - } while (!hulltri.Equal(starttri)); + } while (!hulltri.Equals(starttri)); } /// @@ -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 diff --git a/Triangle.NET/Triangle/NewLocation.cs b/Triangle.NET/Triangle/NewLocation.cs index 68138a0..8be0ee5 100644 --- a/Triangle.NET/Triangle/NewLocation.cs +++ b/Triangle.NET/Triangle/NewLocation.cs @@ -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) diff --git a/Triangle.NET/Triangle/Tools/AdjacencyMatrix.cs b/Triangle.NET/Triangle/Tools/AdjacencyMatrix.cs index 90b5a19..3031b1a 100644 --- a/Triangle.NET/Triangle/Tools/AdjacencyMatrix.cs +++ b/Triangle.NET/Triangle/Tools/AdjacencyMatrix.cs @@ -8,49 +8,74 @@ namespace TriangleNet.Tools { using System; - using System.Collections.Generic; - using System.Linq; - using System.Text; /// /// The adjacency matrix of the mesh. /// 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 + /// + /// Gets the number of columns (nodes of the mesh). + /// + public readonly int N; + + /// + /// Gets the column pointers. + /// + public int[] ColumnPointers { - get { return adj_row; } + get { return pcol; } } - public int[] Adjacency + /// + /// Gets the row indices. + /// + 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(); + } } /// @@ -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 /// 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; } /// @@ -188,28 +197,25 @@ namespace TriangleNet.Tools /// for a linear triangle finite element discretization of Poisson's /// equation in two dimensions. /// - 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 - - /// - /// Reorders an array of integers into a descending heap. - /// - /// the size of the input array. - /// an unsorted array. - /// - /// 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) - /// - 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; - } - - - /// - /// ascending sorts an array of integers using heap sort. - /// - /// Number of entries in the array. - /// Array to be sorted; - 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 } } diff --git a/Triangle.NET/Triangle/Tools/CuthillMcKee.cs b/Triangle.NET/Triangle/Tools/CuthillMcKee.cs index cd924c0..722ef49 100644 --- a/Triangle.NET/Triangle/Tools/CuthillMcKee.cs +++ b/Triangle.NET/Triangle/Tools/CuthillMcKee.cs @@ -8,10 +8,6 @@ namespace TriangleNet.Tools { using System; - using System.Collections.Generic; - using System.Linq; - using System.Text; - using TriangleNet.Logging; /// /// Applies the Cuthill and McKee renumbering algorithm to reduce the bandwidth of @@ -19,9 +15,6 @@ namespace TriangleNet.Tools /// 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 /// Permutation vector. 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(); + /// + /// Gets the permutation vector for the Reverse Cuthill-McKee numbering. + /// + /// The mesh. + /// Permutation vector. + 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; } - /// - /// Computes the bandwidth of a permuted adjacency matrix. - /// - /// The permutation. - /// The inverse permutation. - /// Bandwidth of the permuted adjacency matrix. - /// - /// The matrix is defined by the adjacency information and a permutation. - /// The routine also computes the bandwidth and the size of the envelope. - /// - 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 /// @@ -105,7 +80,10 @@ namespace TriangleNet.Tools /// 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 /// 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. /// - 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. /// - 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 /// 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 + /// + /// Computes the bandwidth of a permuted adjacency matrix. + /// + /// The permutation. + /// The inverse permutation. + /// Bandwidth of the permuted adjacency matrix. + /// + /// The matrix is defined by the adjacency information and a permutation. + /// The routine also computes the bandwidth and the size of the envelope. + /// + 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; + } + /// /// Produces the inverse of a given permutation. /// /// Number of items permuted. /// PERM[N], a permutation. /// The inverse permutation. - 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 } } diff --git a/Triangle.NET/Triangle/Tools/QuadTree.cs b/Triangle.NET/Triangle/Tools/TriangleQuadTree.cs similarity index 93% rename from Triangle.NET/Triangle/Tools/QuadTree.cs rename to Triangle.NET/Triangle/Tools/TriangleQuadTree.cs index 8caefbb..5f7e717 100644 --- a/Triangle.NET/Triangle/Tools/QuadTree.cs +++ b/Triangle.NET/Triangle/Tools/TriangleQuadTree.cs @@ -1,6 +1,6 @@ // ----------------------------------------------------------------------- -// -// Original code by Frank Dockhorn, http://sourceforge.net/projects/quadtreesim/ +// +// Original code by Frank Dockhorn, [not available anymore: http://sourceforge.net/projects/quadtreesim/] // Triangle.NET code by Christian Woltering, http://triangle.codeplex.com/ // // ----------------------------------------------------------------------- @@ -14,7 +14,7 @@ namespace TriangleNet.Tools /// /// A Quadtree implementation optimized for triangles. /// - public class QuadTree + public class TriangleQuadTree { QuadNode root; @@ -24,7 +24,7 @@ namespace TriangleNet.Tools internal int maxDepth; /// - /// Initializes a new instance of the class. + /// Initializes a new instance of the class. /// /// Mesh containing triangles. /// The maximum depth of the tree. @@ -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. /// - 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 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); diff --git a/Triangle.NET/Triangle/Topology/Otri.cs b/Triangle.NET/Triangle/Topology/Otri.cs index 51544f8..cc8da35 100644 --- a/Triangle.NET/Triangle/Topology/Otri.cs +++ b/Triangle.NET/Triangle/Topology/Otri.cs @@ -349,6 +349,14 @@ namespace TriangleNet.Topology ot.orient = orient; } + /// + /// Test for equality of oriented triangles. + /// + 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; } - /// - /// Test for equality of oriented triangles. - /// - internal bool Equal(Otri ot) - { - return ((tri == ot.tri) && (orient == ot.orient)); - } - /// /// Infect a triangle with the virus. /// diff --git a/Triangle.NET/Triangle/Triangle.csproj b/Triangle.NET/Triangle/Triangle.csproj index 9485261..a5c1417 100644 --- a/Triangle.NET/Triangle/Triangle.csproj +++ b/Triangle.NET/Triangle/Triangle.csproj @@ -102,7 +102,7 @@ - + diff --git a/Triangle.NET/Triangle/Voronoi/Legacy/BoundedVoronoiLegacy.cs b/Triangle.NET/Triangle/Voronoi/Legacy/BoundedVoronoiLegacy.cs index 5196c25..c9993a3 100644 --- a/Triangle.NET/Triangle/Voronoi/Legacy/BoundedVoronoiLegacy.cs +++ b/Triangle.NET/Triangle/Voronoi/Legacy/BoundedVoronoiLegacy.cs @@ -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); diff --git a/Triangle.NET/Triangle/Voronoi/Legacy/SimpleVoronoi.cs b/Triangle.NET/Triangle/Voronoi/Legacy/SimpleVoronoi.cs index 8f0333b..7a54adf 100644 --- a/Triangle.NET/Triangle/Voronoi/Legacy/SimpleVoronoi.cs +++ b/Triangle.NET/Triangle/Voronoi/Legacy/SimpleVoronoi.cs @@ -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); diff --git a/Triangle.NET/Triangle/Voronoi/StandardVoronoi.cs b/Triangle.NET/Triangle/Voronoi/StandardVoronoi.cs index 0270e2a..d67e0c2 100644 --- a/Triangle.NET/Triangle/Voronoi/StandardVoronoi.cs +++ b/Triangle.NET/Triangle/Voronoi/StandardVoronoi.cs @@ -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); + } + + /// + /// Compute edge intersections with bounding box. + /// + private void PostProcess(Rectangle box) + { + var infEdges = new List(); + + // 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. + } + } } } }