001//The MIT License (MIT) 002// 003//Copyright (c) 2015 Johannes Diemke 004// 005//Permission is hereby granted, free of charge, to any person obtaining a copy 006//of this software and associated documentation files (the "Software"), to deal 007//in the Software without restriction, including without limitation the rights 008//to use, copy, modify, merge, publish, distribute, sublicense, and/or sell 009//copies of the Software, and to permit persons to whom the Software is 010//furnished to do so, subject to the following conditions: 011// 012//The above copyright notice and this permission notice shall be included in all 013//copies or substantial portions of the Software. 014// 015//THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 016//IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 017//FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE 018//AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 019//LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, 020//OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE 021//SOFTWARE. 022 023package squidpony.squidmath; 024 025import java.io.Serializable; 026import java.util.ArrayList; 027import java.util.Collection; 028import java.util.Comparator; 029import java.util.ListIterator; 030 031/** 032 * A Java implementation of an incremental 2D Delaunay triangulation algorithm. 033 * This is a port of <a href="https://github.com/jdiemke/delaunay-triangulator">Johannes Diemke's code</a>, with 034 * some substantial non-algorithmic changes to better work in SquidLib and to reduce allocations. You should consider 035 * using the {@code com.badlogic.gdx.math.DelaunayTriangulator} class if you use libGDX, which allocates fewer objects, 036 * or {@link IndexedDelaunayTriangulator} for a version of that libGDX class ported to use doubles. 037 * @author Johannes Diemke 038 * @author Tommy Ettinger 039 */ 040public class DelaunayTriangulator implements Serializable { 041 private static final long serialVersionUID = 1L; 042 043 private CoordDouble[] points; 044 private ArrayList<Triangle> triangleSoup; 045 046 /** 047 * Constructs a triangulator instance but does not insert any points; you should add points to 048 * {@link #getPoints()}, which is an array that can hold 256 points, before running {@link #triangulate()}. 049 */ 050 public DelaunayTriangulator() { 051 this.points = new CoordDouble[256]; 052 this.triangleSoup = new ArrayList<>(256); 053 } 054 055 /** 056 * Constructs a new triangulator instance using the specified point set. 057 * 058 * @param points The point set to be triangulated 059 */ 060 public DelaunayTriangulator(Collection<CoordDouble> points) { 061 this.points = points.toArray(new CoordDouble[0]); 062 this.triangleSoup = new ArrayList<>(points.size()); 063 } 064 /** 065 * Returns the triangle from this triangle soup that contains the specified 066 * point or null if no triangle from the triangle soup contains the point. 067 * 068 * @param point 069 * The point 070 * @return Returns the triangle from this triangle soup that contains the 071 * specified point or null 072 */ 073 public Triangle findContainingTriangle(CoordDouble point) { 074 for (Triangle triangle : triangleSoup) { 075 if (triangle.contains(point)) { 076 return triangle; 077 } 078 } 079 return null; 080 } 081 /** 082 * Returns the neighbor triangle of the specified triangle sharing the same 083 * edge as specified. If no neighbor sharing the same edge exists null is 084 * returned. 085 * 086 * @param triangle 087 * The triangle 088 * @param edge 089 * The edge 090 * @return The triangles neighbor triangle sharing the same edge or null if 091 * no triangle exists 092 */ 093 public Triangle findNeighbor(Triangle triangle, Edge edge) { 094 for (Triangle triangleFromSoup : triangleSoup) { 095 if (triangleFromSoup.isNeighbor(edge) && triangleFromSoup != triangle) { 096 return triangleFromSoup; 097 } 098 } 099 return null; 100 } 101 public Triangle findNeighbor(Triangle triangle, CoordDouble ea, CoordDouble eb) { 102 for (Triangle triangleFromSoup : triangleSoup) { 103 if (triangleFromSoup.isNeighbor(ea, eb) && triangleFromSoup != triangle) { 104 return triangleFromSoup; 105 } 106 } 107 return null; 108 } 109 110 /** 111 * Returns one of the possible triangles sharing the specified edge. Based 112 * on the ordering of the triangles in this triangle soup the returned 113 * triangle may differ. To find the other triangle that shares this edge use 114 * the {@link #findNeighbor(Triangle, Edge)} method. 115 * 116 * @param edge 117 * The edge 118 * @return Returns one triangle that shares the specified edge 119 */ 120 public Triangle findOneTriangleSharing(Edge edge) { 121 for (Triangle triangle : triangleSoup) { 122 if (triangle.isNeighbor(edge)) { 123 return triangle; 124 } 125 } 126 return null; 127 } 128 129 /** 130 * Returns the edge from the triangle soup nearest to the specified point. 131 * 132 * @param point 133 * The point 134 * @return The edge from the triangle soup nearest to the specified point 135 */ 136 public Edge findNearestEdge(CoordDouble point) { 137// List<EdgeDistancePack> edgeList = new ArrayList<EdgeDistancePack>(); 138 139 OrderedMap<Edge, Double> edges = new OrderedMap<>(triangleSoup.size()); 140 for (Triangle triangle : triangleSoup) { 141 Edge ab = new Edge(triangle.a, triangle.b); 142 Edge bc = new Edge(triangle.b, triangle.c); 143 Edge ca = new Edge(triangle.c, triangle.a); 144 double abd = triangle.computeClosestPoint(ab, point).subtract(point).lengthSq(); 145 double bcd = triangle.computeClosestPoint(bc, point).subtract(point).lengthSq(); 146 double cad = triangle.computeClosestPoint(ca, point).subtract(point).lengthSq(); 147 148 if(abd <= bcd && abd <= cad) 149 edges.put(ab, abd); 150 else if(bcd <= abd && bcd <= cad) 151 edges.put(bc, bcd); 152 else 153 edges.put(ca, cad); 154 } 155 edges.sortByValue(doubleComparator); 156 return edges.keyAt(0); 157 } 158 159 /** 160 * Removes all triangles from this triangle soup that contain the specified 161 * vertex. 162 * 163 * @param vertex 164 * The vertex 165 */ 166 public void removeTrianglesUsing(CoordDouble vertex) { 167 ListIterator<Triangle> li = triangleSoup.listIterator(); 168 while (li.hasNext()) 169 { 170 Triangle triangle = li.next(); 171 if(triangle.hasVertex(vertex)) 172 li.remove(); 173 } 174 } 175 176 /** 177 * This method generates a Delaunay triangulation from the specified point 178 * set. 179 */ 180 public ArrayList<Triangle> triangulate() { 181 int numPoints; 182 if (points == null || (numPoints = points.length) < 3) { 183 throw new IllegalArgumentException("Less than three points in point set."); 184 } 185 triangleSoup.clear(); 186 187 /* 188 * In order for the in circumcircle test to not consider the vertices of 189 * the super triangle we have to start out with a big triangle 190 * containing the whole point set. We have to scale the super triangle 191 * to be very large. Otherwise the triangulation is not convex. 192 */ 193 double maxOfAnyCoordinate = 0.0; 194 195 for (int i = 0; i < numPoints; i++) { 196 final CoordDouble pt = points[i]; 197 maxOfAnyCoordinate = Math.max(Math.max(pt.x, pt.y), maxOfAnyCoordinate); 198 } 199 200 maxOfAnyCoordinate *= 48.0; 201 202 CoordDouble p1 = new CoordDouble(0.0, maxOfAnyCoordinate); 203 CoordDouble p2 = new CoordDouble(maxOfAnyCoordinate, 0.0); 204 CoordDouble p3 = new CoordDouble(-maxOfAnyCoordinate, -maxOfAnyCoordinate); 205 206 Triangle superTriangle = new Triangle(p1, p2, p3); 207 208 triangleSoup.add(superTriangle); 209 210 for (int i = 0; i < points.length; i++) { 211 Triangle triangle = findContainingTriangle(points[i]); 212 213 if (triangle == null) { 214 /* 215 If no containing triangle exists, then the vertex is not 216 inside a triangle (this can also happen due to numerical 217 errors) and lies on an edge. In order to find this edge we 218 search all edges of the triangle soup and select the one 219 which is nearest to the point we try to add. This edge is 220 removed and four new edges are added. 221 */ 222 Edge edge = findNearestEdge(points[i]); 223 224 Triangle first = findOneTriangleSharing(edge); 225 Triangle second = findNeighbor(first, edge); 226 227 CoordDouble firstNonEdgeVertex = first.getNonEdgeVertex(edge); 228 CoordDouble secondNonEdgeVertex = second.getNonEdgeVertex(edge); 229 CoordDouble p = points[i]; 230 231 triangleSoup.remove(first); 232 triangleSoup.remove(second); 233 234 Triangle triangle1 = new Triangle(edge.a, firstNonEdgeVertex, p); 235 Triangle triangle2 = new Triangle(edge.b, firstNonEdgeVertex, p); 236 Triangle triangle3 = new Triangle(edge.a, secondNonEdgeVertex, p); 237 Triangle triangle4 = new Triangle(edge.b, secondNonEdgeVertex, p); 238 239 triangleSoup.add(triangle1); 240 triangleSoup.add(triangle2); 241 triangleSoup.add(triangle3); 242 triangleSoup.add(triangle4); 243 244 legalizeEdge(triangle1, edge.a, firstNonEdgeVertex, p); 245 legalizeEdge(triangle2, edge.b, firstNonEdgeVertex, p); 246 legalizeEdge(triangle3, edge.a, secondNonEdgeVertex, p); 247 legalizeEdge(triangle4, edge.b, secondNonEdgeVertex, p); 248 } else { 249 /* 250 * The vertex is inside a triangle. 251 */ 252 CoordDouble a = triangle.a; 253 CoordDouble b = triangle.b; 254 CoordDouble c = triangle.c; 255 CoordDouble p = points[i]; 256 257 triangleSoup.remove(triangle); 258 259 Triangle first = new Triangle(a, b, p); 260 Triangle second = new Triangle(b, c, p); 261 Triangle third = new Triangle(c, a, p); 262 263 triangleSoup.add(first); 264 triangleSoup.add(second); 265 triangleSoup.add(third); 266 267 legalizeEdge(first, a, b, p); 268 legalizeEdge(second, b, c, p); 269 legalizeEdge(third, c, a, p); 270 } 271 } 272 273 /* 274 * Remove all triangles that contain vertices of the super triangle. 275 */ 276 removeTrianglesUsing(superTriangle.a); 277 removeTrianglesUsing(superTriangle.b); 278 removeTrianglesUsing(superTriangle.c); 279 return triangleSoup; 280 } 281 282 /** 283 * This method legalizes edges by recursively flipping all illegal edges. 284 * 285 * @param triangle 286 * The triangle 287 * @param ea 288 * The "a" point on the edge to be legalized 289 * @param eb 290 * The "b" point on the edge to be legalized 291 * @param newVertex 292 * The new vertex 293 */ 294 private void legalizeEdge(Triangle triangle, CoordDouble ea, CoordDouble eb, CoordDouble newVertex) { 295 Triangle neighborTriangle = findNeighbor(triangle, ea, eb); 296 297 /* 298 If the triangle has a neighbor, then legalize the edge 299 */ 300 if (neighborTriangle != null) { 301 if (neighborTriangle.isPointInCircumcircle(newVertex)) { 302 triangleSoup.remove(triangle); 303 triangleSoup.remove(neighborTriangle); 304 305 CoordDouble noneEdgeVertex = neighborTriangle.getNonEdgeVertex(ea, eb); 306 307 Triangle firstTriangle = new Triangle(noneEdgeVertex, ea, newVertex); 308 Triangle secondTriangle = new Triangle(noneEdgeVertex, eb, newVertex); 309 310 triangleSoup.add(firstTriangle); 311 triangleSoup.add(secondTriangle); 312 313 legalizeEdge(firstTriangle, noneEdgeVertex, ea, newVertex); 314 legalizeEdge(secondTriangle, noneEdgeVertex, eb, newVertex); 315 } 316 } 317 } 318 319 /** 320 * Creates a random permutation of the specified point set. Based on the 321 * implementation of the Delaunay algorithm this can speed up the 322 * computation. 323 */ 324 public void shuffle(IRNG rng) { 325 rng.shuffleInPlace(points); 326 } 327 328 /** 329 * Returns the point set in form of a vector of 2D vectors. 330 * 331 * @return Returns the points set. 332 */ 333 public CoordDouble[] getPoints() { 334 return points; 335 } 336 337 /** 338 * Returns the triangles of the triangulation in form of a list of 2D 339 * triangles. 340 * 341 * @return Returns the triangles of the triangulation. 342 */ 343 public ArrayList<Triangle> getTriangles() { 344 return triangleSoup; 345 } 346 347 public static class Edge implements Serializable 348 { 349 private static final long serialVersionUID = 1L; 350 351 public CoordDouble a; 352 public CoordDouble b; 353 354 public Edge() 355 { 356 a = new CoordDouble(0.0, 0.0); 357 b = new CoordDouble(0.0, 1.0); 358 } 359 public Edge(CoordDouble a, CoordDouble b) 360 { 361 this.a = a; 362 this.b = b; 363 } 364 public Edge(double ax, double ay, double bx, double by) 365 { 366 a = new CoordDouble(ax, ay); 367 b = new CoordDouble(bx, by); 368 } 369 } 370 371 public static class Triangle implements Serializable 372 { 373 private static final long serialVersionUID = 1L; 374 375 public CoordDouble a; 376 public CoordDouble b; 377 public CoordDouble c; 378 379 /** 380 * Constructor of the 2D triangle class used to create a new triangle 381 * instance from three 2D vectors describing the triangle's vertices. 382 * 383 * @param a 384 * The first vertex of the triangle 385 * @param b 386 * The second vertex of the triangle 387 * @param c 388 * The third vertex of the triangle 389 */ 390 public Triangle(CoordDouble a, CoordDouble b, CoordDouble c) { 391 this.a = a; 392 this.b = b; 393 this.c = c; 394 } 395 396 /** 397 * Tests if a 2D point lies inside this 2D triangle. See Real-Time Collision 398 * Detection, chap. 5, p. 206. 399 * 400 * @param point 401 * The point to be tested 402 * @return Returns true iff the point lies inside this 2D triangle 403 */ 404 public boolean contains(CoordDouble point) { 405 double px = point.x - a.x; 406 double py = point.y - a.y; 407 double sx = b.x - a.x; 408 double sy = b.y - a.y; 409 double pab = py * sx - px * sy;//point.sub(a).cross(b.sub(a)); 410 px = point.x - b.x; 411 py = point.y - b.y; 412 sx = c.x - b.x; 413 sy = c.y - b.y; 414 double pbc = py * sx - px * sy;//point.sub(b).cross(c.sub(b)); 415 416 if (!hasSameSign(pab, pbc)) { 417 return false; 418 } 419 420 px = point.x - c.x; 421 py = point.y - c.y; 422 sx = a.x - c.x; 423 sy = a.y - c.y; 424 double pca = py * sx - px * sy;//point.sub(c).cross(a.sub(c)); 425 426 return hasSameSign(pab, pca); 427 428 } 429 430 /** 431 * Tests if a given point lies in the circumcircle of this triangle. Let the 432 * triangle ABC appear in counterclockwise (CCW) order. Then when det > 433 * 0, the point lies inside the circumcircle through the three points a, b 434 * and c. If instead det < 0, the point lies outside the circumcircle. 435 * When det = 0, the four points are cocircular. If the triangle is oriented 436 * clockwise (CW) the result is reversed. See Real-Time Collision Detection, 437 * chap. 3, p. 34. 438 * 439 * @param point 440 * The point to be tested 441 * @return Returns true iff the point lies inside the circumcircle through 442 * the three points a, b, and c of the triangle 443 */ 444 public boolean isPointInCircumcircle(CoordDouble point) { 445 double a11 = a.x - point.x; 446 double a21 = b.x - point.x; 447 double a31 = c.x - point.x; 448 449 double a12 = a.y - point.y; 450 double a22 = b.y - point.y; 451 double a32 = c.y - point.y; 452 453 double a13 = (a.x - point.x) * (a.x - point.x) + (a.y - point.y) * (a.y - point.y); 454 double a23 = (b.x - point.x) * (b.x - point.x) + (b.y - point.y) * (b.y - point.y); 455 double a33 = (c.x - point.x) * (c.x - point.x) + (c.y - point.y) * (c.y - point.y); 456 457 double det = a11 * a22 * a33 + a12 * a23 * a31 + a13 * a21 * a32 - a13 * a22 * a31 - a12 * a21 * a33 458 - a11 * a23 * a32; 459 460 if (isOrientedCCW()) { 461 return det > 0.0; 462 } 463 464 return det < 0.0; 465 } 466 467 /** 468 * Test if this triangle is oriented counterclockwise (CCW). Let A, B and C 469 * be three 2D points. If det > 0, C lies to the left of the directed 470 * line AB. Equivalently the triangle ABC is oriented counterclockwise. When 471 * det < 0, C lies to the right of the directed line AB, and the triangle 472 * ABC is oriented clockwise. When det = 0, the three points are colinear. 473 * See Real-Time Collision Detection, chap. 3, p. 32 474 * 475 * @return Returns true iff the triangle ABC is oriented counterclockwise 476 * (CCW) 477 */ 478 public boolean isOrientedCCW() { 479 double a11 = a.x - c.x; 480 double a21 = b.x - c.x; 481 482 double a12 = a.y - c.y; 483 double a22 = b.y - c.y; 484 485 double det = a11 * a22 - a12 * a21; 486 487 return det > 0.0; 488 } 489 490 /** 491 * Returns true if this triangle contains the given edge. 492 * 493 * @param edge 494 * The edge to be tested 495 * @return Returns true if this triangle contains the edge 496 */ 497 public boolean isNeighbor(Edge edge) { 498 return (a == edge.a || b == edge.a || c == edge.a) && (a == edge.b || b == edge.b || c == edge.b); 499 } 500 public boolean isNeighbor(CoordDouble ea, CoordDouble eb) { 501 return (a == ea || b == ea || c == ea) && (a == eb || b == eb || c == eb); 502 } 503 504 /** 505 * Returns the vertex of this triangle that is not part of the given edge. 506 * 507 * @param edge 508 * The edge 509 * @return The vertex of this triangle that is not part of the edge 510 */ 511 public CoordDouble getNonEdgeVertex(Edge edge) { 512 if (a != edge.a && a != edge.b) { 513 return a; 514 } else if (b != edge.a && b != edge.b) { 515 return b; 516 } else if (c != edge.a && c != edge.b) { 517 return c; 518 } 519 520 return null; 521 } 522 public CoordDouble getNonEdgeVertex(CoordDouble ea, CoordDouble eb) { 523 if (a != ea && a != eb) { 524 return a; 525 } else if (b != ea && b != eb) { 526 return b; 527 } else if (c != ea && c != eb) { 528 return c; 529 } 530 531 return null; 532 } 533 534 /** 535 * Returns true if the given vertex is one of the vertices describing this 536 * triangle. 537 * 538 * @param vertex 539 * The vertex to be tested 540 * @return Returns true if the Vertex is one of the vertices describing this 541 * triangle 542 */ 543 public boolean hasVertex(CoordDouble vertex) { 544 return a == vertex || b == vertex || c == vertex; 545 } 546 547 /** 548 * Computes the closest point on the given edge to the specified point. 549 * 550 * @param edge 551 * The edge on which we search the closest point to the specified 552 * point 553 * @param point 554 * The point to which we search the closest point on the edge 555 * @return The closest point on the given edge to the specified point 556 */ 557 CoordDouble computeClosestPoint(Edge edge, CoordDouble point) { 558 //CoordDouble ab = edge.b.sub(edge.a); 559 double abx = edge.b.x = edge.a.x; 560 double aby = edge.b.y = edge.a.y; 561 double t = ((point.x - edge.a.x) * abx + (point.y - edge.a.y) * aby) / (abx * abx + aby * aby); 562 //double t = point.sub(edge.a).dot(ab) / ab.dot(ab); 563 564 if (t < 0.0d) { 565 t = 0.0d; 566 } else if (t > 1.0d) { 567 t = 1.0d; 568 } 569 return new CoordDouble(edge.a.x + abx * t, edge.a.y + aby * t); 570// return edge.a.add(ab.mult(t)); 571 } 572 573 /** 574 * Tests if the two arguments have the same sign. 575 * 576 * @param a 577 * The first floating point argument 578 * @param b 579 * The second floating point argument 580 * @return Returns true iff both arguments have the same sign 581 */ 582 private boolean hasSameSign(double a, double b) { 583 return (a == 0.0 && b == 0.0) || ((a > 0.0) == (b > 0.0)); 584 } 585 586 @Override 587 public String toString() { 588 return "Triangle[" + a + ", " + b + ", " + c + "]"; 589 } 590 } 591 private static final Comparator<Double> doubleComparator = new Comparator<Double>() { 592 @Override 593 public int compare(Double o1, Double o2) { 594 return o1.compareTo(o2); 595 } 596 }; 597 598}