1 /* 2 * Copyright (c) 2016 Vivid Solutions. 3 * 4 * All rights reserved. This program and the accompanying materials 5 * are made available under the terms of the Eclipse Public License 2.0 6 * and Eclipse Distribution License v. 1.0 which accompanies this distribution. 7 * The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html 8 * and the Eclipse Distribution License is available at 9 * 10 * http://www.eclipse.org/org/documents/edl-v10.php. 11 */ 12 13 package org.locationtech.jts.triangulate; 14 15 import java.util.ArrayList; 16 import java.util.Collection; 17 import java.util.Iterator; 18 import java.util.List; 19 20 import org.locationtech.jts.algorithm.ConvexHull; 21 import org.locationtech.jts.geom.Coordinate; 22 import org.locationtech.jts.geom.Envelope; 23 import org.locationtech.jts.geom.Geometry; 24 import org.locationtech.jts.geom.GeometryFactory; 25 import org.locationtech.jts.index.kdtree.KdNode; 26 import org.locationtech.jts.index.kdtree.KdTree; 27 import org.locationtech.jts.triangulate.quadedge.LastFoundQuadEdgeLocator; 28 import org.locationtech.jts.triangulate.quadedge.QuadEdgeSubdivision; 29 import org.locationtech.jts.triangulate.quadedge.Vertex; 30 import org.locationtech.jts.util.Debug; 31 32 33 /** 34 * Computes a Conforming Delaunay Triangulation over a set of sites and a set of 35 * linear constraints. 36 * <p> 37 * A conforming Delaunay triangulation is a true Delaunay triangulation. In it 38 * each constraint segment is present as a union of one or more triangulation 39 * edges. Constraint segments may be subdivided into two or more triangulation 40 * edges by the insertion of additional sites. The additional sites are called 41 * Steiner points, and are necessary to allow the segments to be faithfully 42 * reflected in the triangulation while maintaining the Delaunay property. 43 * Another way of stating this is that in a conforming Delaunay triangulation 44 * every constraint segment will be the union of a subset of the triangulation 45 * edges (up to tolerance). 46 * <p> 47 * A Conforming Delaunay triangulation is distinct from a Constrained Delaunay triangulation. 48 * A Constrained Delaunay triangulation is not necessarily fully Delaunay, 49 * and it contains the constraint segments exactly as edges of the triangulation. 50 * <p> 51 * A typical usage pattern for the triangulator is: 52 * <pre> 53 * ConformingDelaunayTriangulator cdt = new ConformingDelaunayTriangulator(sites, tolerance); 54 * 55 * // optional 56 * cdt.setSplitPointFinder(splitPointFinder); 57 * cdt.setVertexFactory(vertexFactory); 58 * 59 * cdt.setConstraints(segments, new ArrayList(vertexMap.values())); 60 * cdt.formInitialDelaunay(); 61 * cdt.enforceConstraints(); 62 * subdiv = cdt.getSubdivision(); 63 * </pre> 64 * 65 * @author David Skea 66 * @author Martin Davis 67 */ 68 public class ConformingDelaunayTriangulator 69 { computeVertexEnvelope(Collection vertices)70 private static Envelope computeVertexEnvelope(Collection vertices) { 71 Envelope env = new Envelope(); 72 for (Iterator i = vertices.iterator(); i.hasNext();) { 73 Vertex v = (Vertex) i.next(); 74 env.expandToInclude(v.getCoordinate()); 75 } 76 return env; 77 } 78 79 private List initialVertices; // List<Vertex> 80 private List segVertices; // List<Vertex> 81 82 // MD - using a Set doesn't seem to be much faster 83 // private Set segments = new HashSet(); 84 private List segments = new ArrayList(); // List<Segment> 85 private QuadEdgeSubdivision subdiv = null; 86 private IncrementalDelaunayTriangulator incDel; 87 private Geometry convexHull; 88 private ConstraintSplitPointFinder splitFinder = new NonEncroachingSplitPointFinder(); 89 private KdTree kdt = null; 90 private ConstraintVertexFactory vertexFactory = null; 91 92 // allPointsEnv expanded by a small buffer 93 private Envelope computeAreaEnv; 94 // records the last split point computed, for error reporting 95 private Coordinate splitPt = null; 96 97 private double tolerance; // defines if two sites are the same. 98 99 /** 100 * Creates a Conforming Delaunay Triangulation based on the given 101 * unconstrained initial vertices. The initial vertex set should not contain 102 * any vertices which appear in the constraint set. 103 * 104 * @param initialVertices 105 * a collection of {@link ConstraintVertex} 106 * @param tolerance 107 * the distance tolerance below which points are considered identical 108 */ ConformingDelaunayTriangulator(Collection initialVertices, double tolerance)109 public ConformingDelaunayTriangulator(Collection initialVertices, 110 double tolerance) { 111 this.initialVertices = new ArrayList(initialVertices); 112 this.tolerance = tolerance; 113 kdt = new KdTree(tolerance); 114 } 115 116 /** 117 * Sets the constraints to be conformed to by the computed triangulation. 118 * The constraints must not contain duplicate segments (up to orientation). 119 * The unique set of vertices (as {@link ConstraintVertex}es) 120 * forming the constraints must also be supplied. 121 * Supplying it explicitly allows the ConstraintVertexes to be initialized 122 * appropriately (e.g. with external data), and avoids re-computing the unique set 123 * if it is already available. 124 * 125 * @param segments a list of the constraint {@link Segment}s 126 * @param segVertices the set of unique {@link ConstraintVertex}es referenced by the segments 127 */ setConstraints(List segments, List segVertices)128 public void setConstraints(List segments, List segVertices) { 129 this.segments = segments; 130 this.segVertices = segVertices; 131 } 132 133 /** 134 * Sets the {@link ConstraintSplitPointFinder} to be 135 * used during constraint enforcement. 136 * Different splitting strategies may be appropriate 137 * for special situations. 138 * 139 * @param splitFinder the ConstraintSplitPointFinder to be used 140 */ setSplitPointFinder(ConstraintSplitPointFinder splitFinder)141 public void setSplitPointFinder(ConstraintSplitPointFinder splitFinder) { 142 this.splitFinder = splitFinder; 143 } 144 145 /** 146 * Gets the tolerance value used to construct the triangulation. 147 * 148 * @return a tolerance value 149 */ getTolerance()150 public double getTolerance() 151 { 152 return tolerance; 153 } 154 155 /** 156 * Gets the <tt>ConstraintVertexFactory</tt> used to create new constraint vertices at split points. 157 * 158 * @return a new constraint vertex 159 */ getVertexFactory()160 public ConstraintVertexFactory getVertexFactory() { 161 return vertexFactory; 162 } 163 164 /** 165 * Sets a custom {@link ConstraintVertexFactory} to be used 166 * to allow vertices carrying extra information to be created. 167 * 168 * @param vertexFactory the ConstraintVertexFactory to be used 169 */ setVertexFactory(ConstraintVertexFactory vertexFactory)170 public void setVertexFactory(ConstraintVertexFactory vertexFactory) { 171 this.vertexFactory = vertexFactory; 172 } 173 174 /** 175 * Gets the {@link QuadEdgeSubdivision} which represents the triangulation. 176 * 177 * @return a subdivision 178 */ getSubdivision()179 public QuadEdgeSubdivision getSubdivision() { 180 return subdiv; 181 } 182 183 /** 184 * Gets the {@link KdTree} which contains the vertices of the triangulation. 185 * 186 * @return a KdTree 187 */ getKDT()188 public KdTree getKDT() { 189 return kdt; 190 } 191 192 /** 193 * Gets the sites (vertices) used to initialize the triangulation. 194 * 195 * @return a List of Vertex 196 */ getInitialVertices()197 public List getInitialVertices() { 198 return initialVertices; 199 } 200 201 /** 202 * Gets the {@link Segment}s which represent the constraints. 203 * 204 * @return a collection of Segments 205 */ getConstraintSegments()206 public Collection getConstraintSegments() { 207 return segments; 208 } 209 210 /** 211 * Gets the convex hull of all the sites in the triangulation, 212 * including constraint vertices. 213 * Only valid after the constraints have been enforced. 214 * 215 * @return the convex hull of the sites 216 */ getConvexHull()217 public Geometry getConvexHull() { 218 return convexHull; 219 } 220 221 // ================================================================== 222 computeBoundingBox()223 private void computeBoundingBox() { 224 Envelope vertexEnv = computeVertexEnvelope(initialVertices); 225 Envelope segEnv = computeVertexEnvelope(segVertices); 226 227 Envelope allPointsEnv = new Envelope(vertexEnv); 228 allPointsEnv.expandToInclude(segEnv); 229 230 double deltaX = allPointsEnv.getWidth() * 0.2; 231 double deltaY = allPointsEnv.getHeight() * 0.2; 232 233 double delta = Math.max(deltaX, deltaY); 234 235 computeAreaEnv = new Envelope(allPointsEnv); 236 computeAreaEnv.expandBy(delta); 237 } 238 computeConvexHull()239 private void computeConvexHull() { 240 GeometryFactory fact = new GeometryFactory(); 241 Coordinate[] coords = getPointArray(); 242 ConvexHull hull = new ConvexHull(coords, fact); 243 convexHull = hull.getConvexHull(); 244 } 245 246 // /** 247 // * Adds the segments in the Convex Hull of all sites in the input data as 248 // linear constraints. 249 // * This is required if TIN Refinement is performed. The hull segments are 250 // flagged with a 251 // unique 252 // * data object to allow distinguishing them. 253 // * 254 // * @param convexHullSegmentData the data object to attach to each convex 255 // hull segment 256 // */ 257 // private void addConvexHullToConstraints(Object convexHullSegmentData) { 258 // Coordinate[] coords = convexHull.getCoordinates(); 259 // for (int i = 1; i < coords.length; i++) { 260 // Segment s = new Segment(coords[i - 1], coords[i], convexHullSegmentData); 261 // addConstraintIfUnique(s); 262 // } 263 // } 264 265 // private void addConstraintIfUnique(Segment r) { 266 // boolean exists = false; 267 // Iterator it = segments.iterator(); 268 // Segment s = null; 269 // while (it.hasNext()) { 270 // s = (Segment) it.next(); 271 // if (r.equalsTopo(s)) { 272 // exists = true; 273 // } 274 // } 275 // if (!exists) { 276 // segments.add((Object) r); 277 // } 278 // } 279 getPointArray()280 private Coordinate[] getPointArray() { 281 Coordinate[] pts = new Coordinate[initialVertices.size() 282 + segVertices.size()]; 283 int index = 0; 284 for (Iterator i = initialVertices.iterator(); i.hasNext();) { 285 Vertex v = (Vertex) i.next(); 286 pts[index++] = v.getCoordinate(); 287 } 288 for (Iterator i2 = segVertices.iterator(); i2.hasNext();) { 289 Vertex v = (Vertex) i2.next(); 290 pts[index++] = v.getCoordinate(); 291 } 292 return pts; 293 } 294 createVertex(Coordinate p)295 private ConstraintVertex createVertex(Coordinate p) { 296 ConstraintVertex v = null; 297 if (vertexFactory != null) 298 v = vertexFactory.createVertex(p, null); 299 else 300 v = new ConstraintVertex(p); 301 return v; 302 } 303 304 /** 305 * Creates a vertex on a constraint segment 306 * 307 * @param p the location of the vertex to create 308 * @param seg the constraint segment it lies on 309 * @return the new constraint vertex 310 */ createVertex(Coordinate p, Segment seg)311 private ConstraintVertex createVertex(Coordinate p, Segment seg) { 312 ConstraintVertex v = null; 313 if (vertexFactory != null) 314 v = vertexFactory.createVertex(p, seg); 315 else 316 v = new ConstraintVertex(p); 317 v.setOnConstraint(true); 318 return v; 319 } 320 321 /** 322 * Inserts all sites in a collection 323 * 324 * @param vertices a collection of ConstraintVertex 325 */ insertSites(Collection vertices)326 private void insertSites(Collection vertices) { 327 Debug.println("Adding sites: " + vertices.size()); 328 for (Iterator i = vertices.iterator(); i.hasNext();) { 329 ConstraintVertex v = (ConstraintVertex) i.next(); 330 insertSite(v); 331 } 332 } 333 insertSite(ConstraintVertex v)334 private ConstraintVertex insertSite(ConstraintVertex v) { 335 KdNode kdnode = kdt.insert(v.getCoordinate(), v); 336 if (!kdnode.isRepeated()) { 337 incDel.insertSite(v); 338 } else { 339 ConstraintVertex snappedV = (ConstraintVertex) kdnode.getData(); 340 snappedV.merge(v); 341 return snappedV; 342 // testing 343 // if ( v.isOnConstraint() && ! currV.isOnConstraint()) { 344 // System.out.println(v); 345 // } 346 } 347 return v; 348 } 349 350 /** 351 * Inserts a site into the triangulation, maintaining the conformal Delaunay property. 352 * This can be used to further refine the triangulation if required 353 * (e.g. to approximate the medial axis of the constraints, 354 * or to improve the grading of the triangulation). 355 * 356 * @param p the location of the site to insert 357 */ insertSite(Coordinate p)358 public void insertSite(Coordinate p) { 359 insertSite(createVertex(p)); 360 } 361 362 // ================================================================== 363 364 /** 365 * Computes the Delaunay triangulation of the initial sites. 366 */ formInitialDelaunay()367 public void formInitialDelaunay() { 368 computeBoundingBox(); 369 subdiv = new QuadEdgeSubdivision(computeAreaEnv, tolerance); 370 subdiv.setLocator(new LastFoundQuadEdgeLocator(subdiv)); 371 incDel = new IncrementalDelaunayTriangulator(subdiv); 372 insertSites(initialVertices); 373 } 374 375 // ================================================================== 376 377 private final static int MAX_SPLIT_ITER = 99; 378 379 /** 380 * Enforces the supplied constraints into the triangulation. 381 * 382 * @throws ConstraintEnforcementException 383 * if the constraints cannot be enforced 384 */ enforceConstraints()385 public void enforceConstraints() { 386 addConstraintVertices(); 387 // if (true) return; 388 389 int count = 0; 390 int splits = 0; 391 do { 392 splits = enforceGabriel(segments); 393 394 count++; 395 Debug.println("Iter: " + count + " Splits: " + splits 396 + " Current # segments = " + segments.size()); 397 } while (splits > 0 && count < MAX_SPLIT_ITER); 398 if (count == MAX_SPLIT_ITER) { 399 Debug.println("ABORTED! Too many iterations while enforcing constraints"); 400 if (!Debug.isDebugging()) 401 throw new ConstraintEnforcementException( 402 "Too many splitting iterations while enforcing constraints. Last split point was at: ", 403 splitPt); 404 } 405 } 406 addConstraintVertices()407 private void addConstraintVertices() { 408 computeConvexHull(); 409 // insert constraint vertices as sites 410 insertSites(segVertices); 411 } 412 413 /* 414 * private List findMissingConstraints() { List missingSegs = new ArrayList(); 415 * for (int i = 0; i < segments.size(); i++) { Segment s = (Segment) 416 * segments.get(i); QuadEdge q = subdiv.locate(s.getStart(), s.getEnd()); if 417 * (q == null) missingSegs.add(s); } return missingSegs; } 418 */ 419 enforceGabriel(Collection segsToInsert)420 private int enforceGabriel(Collection segsToInsert) { 421 List newSegments = new ArrayList(); 422 int splits = 0; 423 List segsToRemove = new ArrayList(); 424 425 /** 426 * On each iteration must always scan all constraint (sub)segments, since 427 * some constraints may be rebroken by Delaunay triangle flipping caused by 428 * insertion of another constraint. However, this process must converge 429 * eventually, with no splits remaining to find. 430 */ 431 for (Iterator i = segsToInsert.iterator(); i.hasNext();) { 432 Segment seg = (Segment) i.next(); 433 // System.out.println(seg); 434 435 Coordinate encroachPt = findNonGabrielPoint(seg); 436 // no encroachment found - segment must already be in subdivision 437 if (encroachPt == null) 438 continue; 439 440 // compute split point 441 splitPt = splitFinder.findSplitPoint(seg, encroachPt); 442 ConstraintVertex splitVertex = createVertex(splitPt, seg); 443 444 // DebugFeature.addLineSegment(DEBUG_SEG_SPLIT, encroachPt, splitPt, ""); 445 // Debug.println(WKTWriter.toLineString(encroachPt, splitPt)); 446 447 /** 448 * Check whether the inserted point still equals the split pt. This will 449 * not be the case if the split pt was too close to an existing site. If 450 * the point was snapped, the triangulation will not respect the inserted 451 * constraint - this is a failure. This can be caused by: 452 * <ul> 453 * <li>An initial site that lies very close to a constraint segment The 454 * cure for this is to remove any initial sites which are close to 455 * constraint segments in a preprocessing phase. 456 * <li>A narrow constraint angle which causing repeated splitting until 457 * the split segments are too small. The cure for this is to either choose 458 * better split points or "guard" narrow angles by cracking the segments 459 * equidistant from the corner. 460 * </ul> 461 */ 462 ConstraintVertex insertedVertex = insertSite(splitVertex); 463 if (!insertedVertex.getCoordinate().equals2D(splitPt)) { 464 Debug.println("Split pt snapped to: " + insertedVertex); 465 // throw new ConstraintEnforcementException("Split point snapped to 466 // existing point 467 // (tolerance too large or constraint interior narrow angle?)", 468 // splitPt); 469 } 470 471 // split segment and record the new halves 472 Segment s1 = new Segment(seg.getStartX(), seg.getStartY(), seg 473 .getStartZ(), splitVertex.getX(), splitVertex.getY(), splitVertex 474 .getZ(), seg.getData()); 475 Segment s2 = new Segment(splitVertex.getX(), splitVertex.getY(), 476 splitVertex.getZ(), seg.getEndX(), seg.getEndY(), seg.getEndZ(), seg 477 .getData()); 478 newSegments.add(s1); 479 newSegments.add(s2); 480 segsToRemove.add(seg); 481 482 splits = splits + 1; 483 } 484 segsToInsert.removeAll(segsToRemove); 485 segsToInsert.addAll(newSegments); 486 487 return splits; 488 } 489 490 // public static final String DEBUG_SEG_SPLIT = "C:\\proj\\CWB\\test\\segSplit.jml"; 491 492 /** 493 * Given a set of points stored in the kd-tree and a line segment defined by 494 * two points in this set, finds a {@link Coordinate} in the circumcircle of 495 * the line segment, if one exists. This is called the Gabriel point - if none 496 * exists then the segment is said to have the Gabriel condition. Uses the 497 * heuristic of finding the non-Gabriel point closest to the midpoint of the 498 * segment. 499 * 500 * @param p 501 * start of the line segment 502 * @param q 503 * end of the line segment 504 * @return a point which is non-Gabriel 505 * or null if no point is non-Gabriel 506 */ findNonGabrielPoint(Segment seg)507 private Coordinate findNonGabrielPoint(Segment seg) { 508 Coordinate p = seg.getStart(); 509 Coordinate q = seg.getEnd(); 510 // Find the mid point on the line and compute the radius of enclosing circle 511 Coordinate midPt = new Coordinate((p.x + q.x) / 2.0, (p.y + q.y) / 2.0); 512 double segRadius = p.distance(midPt); 513 514 // compute envelope of circumcircle 515 Envelope env = new Envelope(midPt); 516 env.expandBy(segRadius); 517 // Find all points in envelope 518 List result = kdt.query(env); 519 520 // For each point found, test if it falls strictly in the circle 521 // find closest point 522 Coordinate closestNonGabriel = null; 523 double minDist = Double.MAX_VALUE; 524 for (Iterator i = result.iterator(); i.hasNext();) { 525 KdNode nextNode = (KdNode) i.next(); 526 Coordinate testPt = nextNode.getCoordinate(); 527 // ignore segment endpoints 528 if (testPt.equals2D(p) || testPt.equals2D(q)) 529 continue; 530 531 double testRadius = midPt.distance(testPt); 532 if (testRadius < segRadius) { 533 // double testDist = seg.distance(testPt); 534 double testDist = testRadius; 535 if (closestNonGabriel == null || testDist < minDist) { 536 closestNonGabriel = testPt; 537 minDist = testDist; 538 } 539 } 540 } 541 return closestNonGabriel; 542 } 543 544 } 545