diff --git a/Triangle.NET/Triangle/Geometry/Vertex.cs b/Triangle.NET/Triangle/Geometry/Vertex.cs index bb817f9..49f1b58 100644 --- a/Triangle.NET/Triangle/Geometry/Vertex.cs +++ b/Triangle.NET/Triangle/Geometry/Vertex.cs @@ -22,7 +22,7 @@ namespace TriangleNet.Geometry internal double[] attributes; #endif internal VertexType type; - internal Otri tri; // TODO: change from Otri to Triangle + internal Otri tri; /// /// Initializes a new instance of the class. diff --git a/Triangle.NET/Triangle/Meshing/Algorithm/Dwyer.cs b/Triangle.NET/Triangle/Meshing/Algorithm/Dwyer.cs index a3e7cbb..a4c834d 100644 --- a/Triangle.NET/Triangle/Meshing/Algorithm/Dwyer.cs +++ b/Triangle.NET/Triangle/Meshing/Algorithm/Dwyer.cs @@ -9,8 +9,9 @@ namespace TriangleNet.Meshing.Algorithm { using System; using System.Collections.Generic; - using TriangleNet.Topology; using TriangleNet.Geometry; + using TriangleNet.Tools; + using TriangleNet.Topology; /// /// Builds a delaunay triangulation using the divide-and-conquer algorithm. @@ -71,22 +72,19 @@ namespace TriangleNet.Meshing.Algorithm this.mesh.TransferNodes(points); Otri hullleft = default(Otri), hullright = default(Otri); - int divider; int i, j, n = points.Count; - //DebugWriter.Session.Start("test-dbg"); - // Allocate an array of pointers to vertices for sorting. - // TODO: use ToArray this.sortarray = new Vertex[n]; i = 0; foreach (var v in points) { sortarray[i++] = v; } + // Sort the vertices. - //Array.Sort(sortarray); - VertexSort(0, n - 1); + VertexSorter.Sort(sortarray); + // Discard duplicate vertices, which can really mess up the algorithm. i = 0; for (j = 1; j < n; j++) @@ -112,225 +110,17 @@ namespace TriangleNet.Meshing.Algorithm if (UseDwyer) { // Re-sort the array of vertices to accommodate alternating cuts. - divider = i >> 1; - if (i - divider >= 2) - { - if (divider >= 2) - { - AlternateAxes(0, divider - 1, 1); - } - AlternateAxes(divider, i - 1, 1); - } + VertexSorter.Alternate(sortarray, i); } // Form the Delaunay triangulation. DivconqRecurse(0, i - 1, 0, ref hullleft, ref hullright); - //DebugWriter.Session.Write(mesh); - //DebugWriter.Session.Finish(); - this.mesh.hullsize = RemoveGhosts(ref hullleft); return this.mesh; } - /// - /// Sort an array of vertices by x-coordinate, using the y-coordinate as a secondary key. - /// - /// - /// - /// - /// Uses quicksort. Randomized O(n log n) time. No, I did not make any of - /// the usual quicksort mistakes. - /// - void VertexSort(int left, int right) - { - int oleft = left; - int oright = right; - int arraysize = right - left + 1; - int pivot; - double pivotx, pivoty; - Vertex temp; - - if (arraysize < 32) - { - // Insertion sort - for (int i = left + 1; i <= right; i++) - { - var a = sortarray[i]; - int j = i - 1; - while (j >= left && (sortarray[j].x > a.x || (sortarray[j].x == a.x && sortarray[j].y > a.y))) - { - sortarray[j + 1] = sortarray[j]; - j--; - } - sortarray[j + 1] = a; - } - - return; - } - - // Choose a random pivot to split the array. - pivot = rand.Next(left, right); - pivotx = sortarray[pivot].x; - pivoty = sortarray[pivot].y; - // Split the array. - left--; - right++; - while (left < right) - { - // Search for a vertex whose x-coordinate is too large for the left. - do - { - left++; - } - while ((left <= right) && ((sortarray[left].x < pivotx) || - ((sortarray[left].x == pivotx) && (sortarray[left].y < pivoty)))); - // Search for a vertex whose x-coordinate is too small for the right. - do - { - right--; - } - while ((left <= right) && ((sortarray[right].x > pivotx) || - ((sortarray[right].x == pivotx) && (sortarray[right].y > pivoty)))); - - if (left < right) - { - // Swap the left and right vertices. - temp = sortarray[left]; - sortarray[left] = sortarray[right]; - sortarray[right] = temp; - } - } - if (left > oleft) - { - // Recursively sort the left subset. - VertexSort(oleft, left); - } - if (oright > right + 1) - { - // Recursively sort the right subset. - VertexSort(right + 1, oright); - } - } - - /// - /// An order statistic algorithm, almost. Shuffles an array of vertices so that - /// the first 'median' vertices occur lexicographically before the remaining vertices. - /// - /// - /// - /// - /// - /// - /// Uses the x-coordinate as the primary key if axis == 0; the y-coordinate - /// if axis == 1. Very similar to the vertexsort() procedure, but runs in - /// randomized linear time. - /// - void VertexMedian(int left, int right, int median, int axis) - { - int arraysize = right - left + 1; - int oleft = left, oright = right; - int pivot; - double pivot1, pivot2; - Vertex temp; - - if (arraysize == 2) - { - // Recursive base case. - if ((sortarray[left][axis] > sortarray[right][axis]) || - ((sortarray[left][axis] == sortarray[right][axis]) && - (sortarray[left][1 - axis] > sortarray[right][1 - axis]))) - { - temp = sortarray[right]; - sortarray[right] = sortarray[left]; - sortarray[left] = temp; - } - return; - } - // Choose a random pivot to split the array. - pivot = rand.Next(left, right); //left + arraysize / 2; - pivot1 = sortarray[pivot][axis]; - pivot2 = sortarray[pivot][1 - axis]; - - left--; - right++; - while (left < right) - { - // Search for a vertex whose x-coordinate is too large for the left. - do - { - left++; - } - while ((left <= right) && ((sortarray[left][axis] < pivot1) || - ((sortarray[left][axis] == pivot1) && (sortarray[left][1 - axis] < pivot2)))); - // Search for a vertex whose x-coordinate is too small for the right. - do - { - right--; - } - while ((left <= right) && ((sortarray[right][axis] > pivot1) || - ((sortarray[right][axis] == pivot1) && (sortarray[right][1 - axis] > pivot2)))); - if (left < right) - { - // Swap the left and right vertices. - temp = sortarray[left]; - sortarray[left] = sortarray[right]; - sortarray[right] = temp; - } - } - - // Unlike in vertexsort(), at most one of the following conditionals is true. - if (left > median) - { - // Recursively shuffle the left subset. - VertexMedian(oleft, left - 1, median, axis); - } - if (right < median - 1) - { - // Recursively shuffle the right subset. - VertexMedian(right + 1, oright, median, axis); - } - } - - /// - /// Sorts the vertices as appropriate for the divide-and-conquer algorithm with - /// alternating cuts. - /// - /// - /// - /// - /// - /// Partitions by x-coordinate if axis == 0; by y-coordinate if axis == 1. - /// For the base case, subsets containing only two or three vertices are - /// always sorted by x-coordinate. - /// - void AlternateAxes(int left, int right, int axis) - { - int arraysize = right - left + 1; - int divider; - - divider = arraysize >> 1; - //divider += left; // TODO: check - if (arraysize <= 3) - { - // Recursive base case: subsets of two or three vertices will be - // handled specially, and should always be sorted by x-coordinate. - axis = 0; - } - // Partition with a horizontal or vertical cut. - VertexMedian(left, right, left + divider, axis); - // Recursively partition the subsets with a cross cut. - if (arraysize - divider >= 2) - { - if (divider >= 2) - { - AlternateAxes(left, left + divider - 1, 1 - axis); - } - AlternateAxes(left + divider, right, 1 - axis); - } - } - /// /// Merge two adjacent Delaunay triangulations into a single Delaunay triangulation. /// @@ -834,15 +624,13 @@ namespace TriangleNet.Meshing.Algorithm { // Split the vertices in half. divider = vertices >> 1; + // Recursively triangulate each half. DivconqRecurse(left, left + divider - 1, 1 - axis, ref farleft, ref innerleft); - //DebugWriter.Session.Write(mesh, true); DivconqRecurse(left + divider, right, 1 - axis, ref innerright, ref farright); - //DebugWriter.Session.Write(mesh, true); // Merge the two triangulations into one. MergeHulls(ref farleft, ref innerleft, ref innerright, ref farright, axis); - //DebugWriter.Session.Write(mesh, true); } } diff --git a/Triangle.NET/Triangle/Tools/PointSorter.cs b/Triangle.NET/Triangle/Tools/PointSorter.cs deleted file mode 100644 index d910933..0000000 --- a/Triangle.NET/Triangle/Tools/PointSorter.cs +++ /dev/null @@ -1,111 +0,0 @@ -// ----------------------------------------------------------------------- -// -// Original Triangle code by Jonathan Richard Shewchuk, http://www.cs.cmu.edu/~quake/triangle.html -// Triangle.NET code by Christian Woltering, http://triangle.codeplex.com/ -// -// ----------------------------------------------------------------------- - -namespace TriangleNet.Tools -{ - using System; - using TriangleNet.Geometry; - - /// - /// Sort an array of points using quicksort. - /// - public class PointSorter - { - static Random rand = new Random(57113); - - Point[] points; - - public void Sort(Point[] array) - { - this.points = array; - - VertexSort(0, array.Length - 1); - } - - /// - /// Sort an array of vertices by x-coordinate, using the y-coordinate as a secondary key. - /// - /// - /// - /// - /// Uses quicksort. Randomized O(n log n) time. No, I did not make any of - /// the usual quicksort mistakes. - /// - void VertexSort(int left, int right) - { - int oleft = left; - int oright = right; - int arraysize = right - left + 1; - int pivot; - double pivotx, pivoty; - Point temp; - - var array = this.points; - - if (arraysize < 32) - { - // Insertion sort - for (int i = left + 1; i <= right; i++) - { - var a = array[i]; - int j = i - 1; - while (j >= left && (array[j].x > a.x || (array[j].x == a.x && array[j].y > a.y))) - { - array[j + 1] = array[j]; - j--; - } - array[j + 1] = a; - } - - return; - } - - // Choose a random pivot to split the array. - pivot = rand.Next(left, right); - pivotx = array[pivot].x; - pivoty = array[pivot].y; - // Split the array. - left--; - right++; - while (left < right) - { - // Search for a vertex whose x-coordinate is too large for the left. - do - { - left++; - } - while ((left <= right) && ((array[left].x < pivotx) || - ((array[left].x == pivotx) && (array[left].y < pivoty)))); - // Search for a vertex whose x-coordinate is too small for the right. - do - { - right--; - } - while ((left <= right) && ((array[right].x > pivotx) || - ((array[right].x == pivotx) && (array[right].y > pivoty)))); - - if (left < right) - { - // Swap the left and right vertices. - temp = array[left]; - array[left] = array[right]; - array[right] = temp; - } - } - if (left > oleft) - { - // Recursively sort the left subset. - VertexSort(oleft, left); - } - if (oright > right + 1) - { - // Recursively sort the right subset. - VertexSort(right + 1, oright); - } - } - } -} diff --git a/Triangle.NET/Triangle/Tools/PolygonValidator.cs b/Triangle.NET/Triangle/Tools/PolygonValidator.cs index 56173bb..8ec9f59 100644 --- a/Triangle.NET/Triangle/Tools/PolygonValidator.cs +++ b/Triangle.NET/Triangle/Tools/PolygonValidator.cs @@ -102,7 +102,7 @@ namespace TriangleNet.Tools var points = poly.Points.ToArray(); - (new PointSorter()).Sort(points); + VertexSorter.Sort(points); for (int i = 1; i < points.Length; i++) { diff --git a/Triangle.NET/Triangle/Tools/VertexSorter.cs b/Triangle.NET/Triangle/Tools/VertexSorter.cs new file mode 100644 index 0000000..72e662e --- /dev/null +++ b/Triangle.NET/Triangle/Tools/VertexSorter.cs @@ -0,0 +1,371 @@ +// ----------------------------------------------------------------------- +// +// Original Triangle code by Jonathan Richard Shewchuk, http://www.cs.cmu.edu/~quake/triangle.html +// Triangle.NET code by Christian Woltering, http://triangle.codeplex.com/ +// +// ----------------------------------------------------------------------- + +namespace TriangleNet.Tools +{ + using System; + using TriangleNet.Geometry; + + /// + /// Sort an array of points using quicksort. + /// + public class VertexSorter + { + private const int RANDOM_SEED = 57113; + + Random rand; + + Vertex[] points; + + VertexSorter(Vertex[] points, int seed) + { + this.points = points; + this.rand = new Random(seed); + } + + /// + /// Sorts the given vertex array by x-coordinate. + /// + /// The vertex array. + /// Random seed used for pivoting. + public static void Sort(Vertex[] array, int seed = RANDOM_SEED) + { + var qs = new VertexSorter(array, seed); + + qs.QuickSort(0, array.Length - 1); + } + + /// + /// Impose alternating cuts on given vertex array. + /// + /// The vertex array. + /// The number of vertices to sort. + /// Random seed used for pivoting. + public static void Alternate(Vertex[] array, int length, int seed = RANDOM_SEED) + { + var qs = new VertexSorter(array, seed); + + int divider = length >> 1; + + // Re-sort the array of vertices to accommodate alternating cuts. + if (length - divider >= 2) + { + if (divider >= 2) + { + qs.AlternateAxes(0, divider - 1, 1); + } + + qs.AlternateAxes(divider, length - 1, 1); + } + } + + #region Quicksort + + /// + /// Sort an array of vertices by x-coordinate, using the y-coordinate as a secondary key. + /// + /// + /// + /// + /// Uses quicksort. Randomized O(n log n) time. No, I did not make any of + /// the usual quicksort mistakes. + /// + private void QuickSort(int left, int right) + { + int oleft = left; + int oright = right; + int arraysize = right - left + 1; + int pivot; + double pivotx, pivoty; + Vertex temp; + + var array = this.points; + + if (arraysize < 32) + { + // Insertion sort + for (int i = left + 1; i <= right; i++) + { + var a = array[i]; + int j = i - 1; + while (j >= left && (array[j].x > a.x || (array[j].x == a.x && array[j].y > a.y))) + { + array[j + 1] = array[j]; + j--; + } + array[j + 1] = a; + } + + return; + } + + // Choose a random pivot to split the array. + pivot = rand.Next(left, right); + pivotx = array[pivot].x; + pivoty = array[pivot].y; + // Split the array. + left--; + right++; + while (left < right) + { + // Search for a vertex whose x-coordinate is too large for the left. + do + { + left++; + } + while ((left <= right) && ((array[left].x < pivotx) || + ((array[left].x == pivotx) && (array[left].y < pivoty)))); + + // Search for a vertex whose x-coordinate is too small for the right. + do + { + right--; + } + while ((left <= right) && ((array[right].x > pivotx) || + ((array[right].x == pivotx) && (array[right].y > pivoty)))); + + if (left < right) + { + // Swap the left and right vertices. + temp = array[left]; + array[left] = array[right]; + array[right] = temp; + } + } + + if (left > oleft) + { + // Recursively sort the left subset. + QuickSort(oleft, left); + } + + if (oright > right + 1) + { + // Recursively sort the right subset. + QuickSort(right + 1, oright); + } + } + + #endregion + + #region Alternate axes + + /// + /// Sorts the vertices as appropriate for the divide-and-conquer algorithm with + /// alternating cuts. + /// + /// + /// + /// + /// + /// Partitions by x-coordinate if axis == 0; by y-coordinate if axis == 1. + /// For the base case, subsets containing only two or three vertices are + /// always sorted by x-coordinate. + /// + private void AlternateAxes(int left, int right, int axis) + { + int size = right - left + 1; + int divider = size >> 1; + + if (size <= 3) + { + // Recursive base case: subsets of two or three vertices will be + // handled specially, and should always be sorted by x-coordinate. + axis = 0; + } + + // Partition with a horizontal or vertical cut. + if (axis == 0) + { + VertexMedianX(left, right, left + divider); + } + else + { + VertexMedianY(left, right, left + divider); + } + + // Recursively partition the subsets with a cross cut. + if (size - divider >= 2) + { + if (divider >= 2) + { + AlternateAxes(left, left + divider - 1, 1 - axis); + } + + AlternateAxes(left + divider, right, 1 - axis); + } + } + + /// + /// An order statistic algorithm, almost. Shuffles an array of vertices so that the + /// first 'median' vertices occur lexicographically before the remaining vertices. + /// + /// + /// + /// + /// + /// Uses the x-coordinate as the primary key. Very similar to the QuickSort() + /// procedure, but runs in randomized linear time. + /// + private void VertexMedianX(int left, int right, int median) + { + int arraysize = right - left + 1; + int oleft = left, oright = right; + int pivot; + double pivot1, pivot2; + Vertex temp; + + var array = this.points; + + if (arraysize == 2) + { + // Recursive base case. + if ((array[left].x > array[right].x) || + ((array[left].x == array[right].x) && + (array[left].y > array[right].y))) + { + temp = array[right]; + array[right] = array[left]; + array[left] = temp; + } + return; + } + + // Choose a random pivot to split the array. + pivot = rand.Next(left, right); + pivot1 = array[pivot].x; + pivot2 = array[pivot].y; + + left--; + right++; + while (left < right) + { + // Search for a vertex whose x-coordinate is too large for the left. + do + { + left++; + } + while ((left <= right) && ((array[left].x < pivot1) || + ((array[left].x == pivot1) && (array[left].y < pivot2)))); + + // Search for a vertex whose x-coordinate is too small for the right. + do + { + right--; + } + while ((left <= right) && ((array[right].x > pivot1) || + ((array[right].x == pivot1) && (array[right].y > pivot2)))); + + if (left < right) + { + // Swap the left and right vertices. + temp = array[left]; + array[left] = array[right]; + array[right] = temp; + } + } + + // Unlike in vertexsort(), at most one of the following conditionals is true. + if (left > median) + { + // Recursively shuffle the left subset. + VertexMedianX(oleft, left - 1, median); + } + + if (right < median - 1) + { + // Recursively shuffle the right subset. + VertexMedianX(right + 1, oright, median); + } + } + + /// + /// An order statistic algorithm, almost. Shuffles an array of vertices so that + /// the first 'median' vertices occur lexicographically before the remaining vertices. + /// + /// + /// + /// + /// + /// Uses the y-coordinate as the primary key. Very similar to the QuickSort() + /// procedure, but runs in randomized linear time. + /// + private void VertexMedianY(int left, int right, int median) + { + int arraysize = right - left + 1; + int oleft = left, oright = right; + int pivot; + double pivot1, pivot2; + Vertex temp; + + var array = this.points; + + if (arraysize == 2) + { + // Recursive base case. + if ((array[left].y > array[right].y) || + ((array[left].y == array[right].y) && + (array[left].x > array[right].x))) + { + temp = array[right]; + array[right] = array[left]; + array[left] = temp; + } + return; + } + + // Choose a random pivot to split the array. + pivot = rand.Next(left, right); + pivot1 = array[pivot].y; + pivot2 = array[pivot].x; + + left--; + right++; + while (left < right) + { + // Search for a vertex whose x-coordinate is too large for the left. + do + { + left++; + } + while ((left <= right) && ((array[left].y < pivot1) || + ((array[left].y == pivot1) && (array[left].x < pivot2)))); + + // Search for a vertex whose x-coordinate is too small for the right. + do + { + right--; + } + while ((left <= right) && ((array[right].y > pivot1) || + ((array[right].y == pivot1) && (array[right].x > pivot2)))); + + if (left < right) + { + // Swap the left and right vertices. + temp = array[left]; + array[left] = array[right]; + array[right] = temp; + } + } + + // Unlike in QuickSort(), at most one of the following conditionals is true. + if (left > median) + { + // Recursively shuffle the left subset. + VertexMedianY(oleft, left - 1, median); + } + + if (right < median - 1) + { + // Recursively shuffle the right subset. + VertexMedianY(right + 1, oright, median); + } + } + + #endregion + } +} diff --git a/Triangle.NET/Triangle/Triangle.csproj b/Triangle.NET/Triangle/Triangle.csproj index 79bdc0b..d0a8638 100644 --- a/Triangle.NET/Triangle/Triangle.csproj +++ b/Triangle.NET/Triangle/Triangle.csproj @@ -88,7 +88,7 @@ - +