001package squidpony.squidmath; 002/****************************************************************************** 003 Copyright 2011 See AUTHORS file. 004 005 Licensed under the Apache License, Version 2.0 (the "License"); 006 you may not use this file except in compliance with the License. 007 You may obtain a copy of the License at 008 009 http://www.apache.org/licenses/LICENSE-2.0 010 011 Unless required by applicable law or agreed to in writing, software 012 distributed under the License is distributed on an "AS IS" BASIS, 013 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 See the License for the specific language governing permissions and 015 limitations under the License. 016 */ 017 018/** Delaunay triangulation. Adapted from Paul Bourke's triangulate: http://paulbourke.net/papers/triangulate/ 019 * @author Nathan Sweet */ 020public class IndexedDelaunayTriangulator { 021 static private final double EPSILON = 0.000001; 022 static private final int INSIDE = 0; 023 static private final int COMPLETE = 1; 024 static private final int INCOMPLETE = 2; 025 026 private final IntVLA quicksortStack = new IntVLA(); 027 private double[] sortedPoints; 028 private final IntVLA triangles = new IntVLA(16); 029 private final IntVLA originalIndices = new IntVLA(0); 030 private final IntVLA edges = new IntVLA(); 031 private final ShortVLA complete = new ShortVLA(false, 16); // only a ShortVLA because we don't have BooleanArray 032 private final double[] superTriangle = new double[6]; 033 034 /** @see #computeTriangles(double[], int, int, boolean) */ 035 public IntVLA computeTriangles (double[] polygon, boolean sorted) { 036 return computeTriangles(polygon, 0, polygon.length, sorted); 037 } 038 039 /** Triangulates the given point cloud to a list of triangle indices that make up the Delaunay triangulation. 040 * @param points x,y pairs describing points. Duplicate points will result in undefined behavior. 041 * @param sorted If false, the points will be sorted by the x coordinate, which is required by the triangulation 042 * algorithm. In that case, the input array is not modified, the returned indices are for the input 043 * array, and count*2 additional working memory is needed. 044 * @return triples of indices into the points that describe the triangles in clockwise order. Note the returned array is reused 045 * for later calls to the same method. */ 046 public IntVLA computeTriangles (double[] points, int offset, int count, boolean sorted) { 047 IntVLA triangles = this.triangles; 048 triangles.clear(); 049 if (count < 6) return triangles; 050 triangles.ensureCapacity(count); 051 052 if (!sorted) { 053 if (sortedPoints == null || sortedPoints.length < count) sortedPoints = new double[count]; 054 System.arraycopy(points, offset, sortedPoints, 0, count); 055 points = sortedPoints; 056 offset = 0; 057 sort(points, count); 058 } 059 060 int end = offset + count; 061 062 // Determine bounds for super triangle. 063 double xmin = points[0], ymin = points[1]; 064 double xmax = xmin, ymax = ymin; 065 for (int i = offset + 2; i < end; i++) { 066 double value = points[i]; 067 if (value < xmin) xmin = value; 068 if (value > xmax) xmax = value; 069 i++; 070 value = points[i]; 071 if (value < ymin) ymin = value; 072 if (value > ymax) ymax = value; 073 } 074 double dx = xmax - xmin, dy = ymax - ymin; 075 double dmax = (Math.max(dx, dy)) * 20f; 076 double xmid = (xmax + xmin) / 2f, ymid = (ymax + ymin) / 2f; 077 078 // Setup the super triangle, which contains all points. 079 double[] superTriangle = this.superTriangle; 080 superTriangle[0] = xmid - dmax; 081 superTriangle[1] = ymid - dmax; 082 superTriangle[2] = xmid; 083 superTriangle[3] = ymid + dmax; 084 superTriangle[4] = xmid + dmax; 085 superTriangle[5] = ymid - dmax; 086 087 IntVLA edges = this.edges; 088 edges.ensureCapacity(count / 2); 089 090 ShortVLA complete = this.complete; 091 complete.clear(); 092 complete.ensureCapacity(count); 093 094 // Add super triangle. 095 triangles.add(end); 096 triangles.add(end + 2); 097 triangles.add(end + 4); 098 complete.add((short) 0); 099 100 // Include each point one at a time into the existing mesh. 101 for (int pointIndex = offset; pointIndex < end; pointIndex += 2) { 102 double x = points[pointIndex], y = points[pointIndex + 1]; 103 104 // If x,y lies inside the circumcircle of a triangle, the edges are stored and the triangle removed. 105 int[] trianglesArray = triangles.items; 106 short[] completeArray = complete.items; 107 for (int triangleIndex = triangles.size - 1; triangleIndex >= 0; triangleIndex -= 3) { 108 int completeIndex = triangleIndex / 3; 109 if (completeArray[completeIndex] != 0) continue; 110 int p1 = trianglesArray[triangleIndex - 2]; 111 int p2 = trianglesArray[triangleIndex - 1]; 112 int p3 = trianglesArray[triangleIndex]; 113 double x1, y1, x2, y2, x3, y3; 114 if (p1 >= end) { 115 int i = p1 - end; 116 x1 = superTriangle[i]; 117 y1 = superTriangle[i + 1]; 118 } else { 119 x1 = points[p1]; 120 y1 = points[p1 + 1]; 121 } 122 if (p2 >= end) { 123 int i = p2 - end; 124 x2 = superTriangle[i]; 125 y2 = superTriangle[i + 1]; 126 } else { 127 x2 = points[p2]; 128 y2 = points[p2 + 1]; 129 } 130 if (p3 >= end) { 131 int i = p3 - end; 132 x3 = superTriangle[i]; 133 y3 = superTriangle[i + 1]; 134 } else { 135 x3 = points[p3]; 136 y3 = points[p3 + 1]; 137 } 138 switch (circumCircle(x, y, x1, y1, x2, y2, x3, y3)) { 139 case COMPLETE: 140 completeArray[completeIndex] = 1; 141 break; 142 case INSIDE: 143 edges.add(p1); 144 edges.add(p2); 145 edges.add(p2); 146 edges.add(p3); 147 edges.add(p3); 148 edges.add(p1); 149 150 trianglesArray[triangleIndex] = trianglesArray[--triangles.size]; 151 trianglesArray[triangleIndex - 1] = trianglesArray[--triangles.size]; 152 trianglesArray[triangleIndex - 2] = trianglesArray[--triangles.size]; 153 completeArray[completeIndex] = completeArray[--complete.size]; 154 break; 155 } 156 } 157 158 int[] edgesArray = edges.items; 159 for (int i = 0, n = edges.size; i < n; i += 2) { 160 // Skip multiple edges. If all triangles are anticlockwise then all interior edges are opposite pointing in direction. 161 int p1 = edgesArray[i]; 162 if (p1 == -1) continue; 163 int p2 = edgesArray[i + 1]; 164 boolean skip = false; 165 for (int ii = i + 2; ii < n; ii += 2) { 166 if (p1 == edgesArray[ii + 1] && p2 == edgesArray[ii]) { 167 skip = true; 168 edgesArray[ii] = -1; 169 } 170 } 171 if (skip) continue; 172 173 // Form new triangles for the current point. Edges are arranged in clockwise order. 174 triangles.add(p1); 175 triangles.add(edgesArray[i + 1]); 176 triangles.add(pointIndex); 177 complete.add((short)0); 178 } 179 edges.clear(); 180 } 181 182 // Remove triangles with super triangle vertices. 183 int[] trianglesArray = triangles.items; 184 for (int i = triangles.size - 1; i >= 0; i -= 3) { 185 if (trianglesArray[i] >= end || trianglesArray[i - 1] >= end || trianglesArray[i - 2] >= end) { 186 trianglesArray[i] = trianglesArray[--triangles.size]; 187 trianglesArray[i - 1] = trianglesArray[--triangles.size]; 188 trianglesArray[i - 2] = trianglesArray[--triangles.size]; 189 } 190 } 191 192 // Convert sorted to unsorted indices. 193 if (!sorted) { 194 int[] originalIndicesArray = originalIndices.items; 195 for (int i = 0, n = triangles.size; i < n; i++) 196 trianglesArray[i] = (originalIndicesArray[trianglesArray[i] / 2] * 2); 197 } 198 199 // Adjust triangles to start from zero and count by 1, not by vertex x,y coordinate pairs. 200 if (offset == 0) { 201 for (int i = 0, n = triangles.size; i < n; i++) 202 trianglesArray[i] = (trianglesArray[i] / 2); 203 } else { 204 for (int i = 0, n = triangles.size; i < n; i++) 205 trianglesArray[i] = ((trianglesArray[i] - offset) / 2); 206 } 207 208 return triangles; 209 } 210 211 /** Returns INSIDE if point xp,yp is inside the circumcircle made up of the points x1,y1, x2,y2, x3,y3. Returns COMPLETE if xp 212 * is to the right of the entire circumcircle. Otherwise returns INCOMPLETE. Note: a point on the circumcircle edge is 213 * considered inside. */ 214 private int circumCircle (double xp, double yp, double x1, double y1, double x2, double y2, double x3, double y3) { 215 double xc, yc; 216 double y1y2 = Math.abs(y1 - y2); 217 double y2y3 = Math.abs(y2 - y3); 218 if (y1y2 < EPSILON) { 219 if (y2y3 < EPSILON) return INCOMPLETE; 220 double m2 = -(x3 - x2) / (y3 - y2); 221 double mx2 = (x2 + x3) / 2f; 222 double my2 = (y2 + y3) / 2f; 223 xc = (x2 + x1) / 2f; 224 yc = m2 * (xc - mx2) + my2; 225 } else { 226 double m1 = -(x2 - x1) / (y2 - y1); 227 double mx1 = (x1 + x2) / 2f; 228 double my1 = (y1 + y2) / 2f; 229 if (y2y3 < EPSILON) { 230 xc = (x3 + x2) / 2f; 231 yc = m1 * (xc - mx1) + my1; 232 } else { 233 double m2 = -(x3 - x2) / (y3 - y2); 234 double mx2 = (x2 + x3) / 2f; 235 double my2 = (y2 + y3) / 2f; 236 xc = (m1 * mx1 - m2 * mx2 + my2 - my1) / (m1 - m2); 237 yc = m1 * (xc - mx1) + my1; 238 } 239 } 240 241 double dx = x2 - xc; 242 double dy = y2 - yc; 243 double rsqr = dx * dx + dy * dy; 244 245 dx = xp - xc; 246 dx *= dx; 247 dy = yp - yc; 248 if (dx + dy * dy - rsqr <= EPSILON) return INSIDE; 249 return xp > xc && dx > rsqr ? COMPLETE : INCOMPLETE; 250 } 251 252// public void sortPairs(double[] points) 253// { 254// final int pointCount = points.length / 2; 255// if (sortedPoints == null || sortedPoints.length < points.length) sortedPoints = new double[points.length]; 256// System.arraycopy(points, 0, sortedPoints, 0, points.length); 257//// points = sortedPoints; 258// sort(sortedPoints, sortedPoints.length); 259// int[] originalIndicesArray = originalIndices.items; 260// int p; 261// for (int i = 0; i < pointCount; i++) { 262// p = originalIndicesArray[i]; 263// points[p<<1] = sortedPoints[i<<1]; 264// points[p<<1|1] = sortedPoints[i<<1|1]; 265// } 266// System.out.println(sortedPoints.length); 267// } 268 269 /** Sorts x,y pairs of values by the x value. 270 * @param count Number of indices, must be even. */ 271 private void sort (double[] values, int count) { 272 int pointCount = count / 2; 273 originalIndices.clear(); 274 originalIndices.ensureCapacity(pointCount); 275 int[] originalIndicesArray = originalIndices.items; 276 for (int i = 0; i < pointCount; i++) 277 originalIndicesArray[i] = i; 278 279 int lower = 0; 280 int upper = count - 1; 281 IntVLA stack = quicksortStack; 282 stack.add(lower); 283 stack.add(upper - 1); 284 while (stack.size > 0) { 285 upper = stack.pop(); 286 lower = stack.pop(); 287 if (upper <= lower) continue; 288 int i = quicksortPartition(values, lower, upper, originalIndicesArray); 289 if (i - lower > upper - i) { 290 stack.add(lower); 291 stack.add(i - 2); 292 } 293 stack.add(i + 2); 294 stack.add(upper); 295 if (upper - i >= i - lower) { 296 stack.add(lower); 297 stack.add(i - 2); 298 } 299 } 300 } 301 302 private int quicksortPartition (final double[] values, int lower, int upper, int[] originalIndices) { 303 double value = values[lower]; 304 int up = upper; 305 int down = lower + 2; 306 double tempValue; 307 int tempIndex; 308 while (down < up) { 309 while (down < up && values[down] <= value) 310 down = down + 2; 311 while (values[up] > value) 312 up = up - 2; 313 if (down < up) { 314 tempValue = values[down]; 315 values[down] = values[up]; 316 values[up] = tempValue; 317 318 tempValue = values[down + 1]; 319 values[down + 1] = values[up + 1]; 320 values[up + 1] = tempValue; 321 322 tempIndex = originalIndices[down / 2]; 323 originalIndices[down / 2] = originalIndices[up / 2]; 324 originalIndices[up / 2] = tempIndex; 325 } 326 } 327 values[lower] = values[up]; 328 values[up] = value; 329 330 tempValue = values[lower + 1]; 331 values[lower + 1] = values[up + 1]; 332 values[up + 1] = tempValue; 333 334 tempIndex = originalIndices[lower / 2]; 335 originalIndices[lower / 2] = originalIndices[up / 2]; 336 originalIndices[up / 2] = tempIndex; 337 return up; 338 } 339 340 /** Removes all triangles with a centroid outside the specified hull, which may be concave. Note some triangulations may have 341 * triangles whose centroid is inside the hull but a portion is outside. */ 342 public void trim (IntVLA triangles, double[] points, double[] hull, int offset, int count) { 343 int[] trianglesArray = triangles.items; 344 for (int i = triangles.size - 1; i >= 0; i -= 3) { 345 int p1 = trianglesArray[i - 2] * 2; 346 int p2 = trianglesArray[i - 1] * 2; 347 int p3 = trianglesArray[i] * 2; 348 double centroidX = (points[p1] + points[p2] + points[p3]) / 3.0; 349 double centroidY = (points[p1+1] + points[p2+1] + points[p3+1]) / 3.0; 350 if (!isPointInPolygon(hull, offset, count, centroidX, centroidY)) { 351 trianglesArray[i] = trianglesArray[--triangles.size]; 352 trianglesArray[i - 1] = trianglesArray[--triangles.size]; 353 trianglesArray[i - 2] = trianglesArray[--triangles.size]; 354 } 355 } 356 } 357 /** Returns true if the specified point is in the polygon. 358 * @param offset Starting polygon index. 359 * @param count Number of array indices to use after offset. */ 360 public static boolean isPointInPolygon (double[] polygon, int offset, int count, double x, double y) { 361 boolean oddNodes = false; 362 int j = offset + count - 2; 363 for (int i = offset, n = j; i <= n; i += 2) { 364 double yi = polygon[i + 1]; 365 double yj = polygon[j + 1]; 366 if ((yi < y && yj >= y) || (yj < y && yi >= y)) { 367 double xi = polygon[i]; 368 if (xi + (y - yi) / (yj - yi) * (polygon[j] - xi) < x) oddNodes = !oddNodes; 369 } 370 j = i; 371 } 372 return oddNodes; 373 } 374} 375