Doc Classic

Class ConformingDelaunayTriangulator

  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 {
 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      */
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      */
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      */
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      */
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      */
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      */
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      */
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      */
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      */
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      */
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      */
217     public Geometry getConvexHull() {
218         return convexHull;
219     }
220  
221     // ==================================================================
222  
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  
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  
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  
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      */
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      */
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  
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      */
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      */
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      */
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  
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  
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      */
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