/*
* Copyright (c) 2016 Vivid Solutions.
*
* All rights reserved. This program and the accompanying materials
* are made available under the terms of the Eclipse Public License v1.0
* and Eclipse Distribution License v. 1.0 which accompanies this distribution.
* The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v10.html
* and the Eclipse Distribution License is available at
*
* http://www.eclipse.org/org/documents/edl-v10.php.
*/
package org.locationtech.jts.triangulate;
import java.util.ArrayList;
import java.util.Collection;
import java.util.Iterator;
import java.util.List;
import org.locationtech.jts.algorithm.ConvexHull;
import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.Envelope;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.index.kdtree.KdNode;
import org.locationtech.jts.index.kdtree.KdTree;
import org.locationtech.jts.triangulate.quadedge.LastFoundQuadEdgeLocator;
import org.locationtech.jts.triangulate.quadedge.QuadEdgeSubdivision;
import org.locationtech.jts.triangulate.quadedge.Vertex;
import org.locationtech.jts.util.Debug;
/**
* Computes a Conforming Delaunay Triangulation over a set of sites and a set of
* linear constraints.
* <p>
* A conforming Delaunay triangulation is a true Delaunay triangulation. In it
* each constraint segment is present as a union of one or more triangulation
* edges. Constraint segments may be subdivided into two or more triangulation
* edges by the insertion of additional sites. The additional sites are called
* Steiner points, and are necessary to allow the segments to be faithfully
* reflected in the triangulation while maintaining the Delaunay property.
* Another way of stating this is that in a conforming Delaunay triangulation
* every constraint segment will be the union of a subset of the triangulation
* edges (up to tolerance).
* <p>
* A Conforming Delaunay triangulation is distinct from a Constrained Delaunay triangulation.
* A Constrained Delaunay triangulation is not necessarily fully Delaunay,
* and it contains the constraint segments exactly as edges of the triangulation.
* <p>
* A typical usage pattern for the triangulator is:
* <pre>
* ConformingDelaunayTriangulator cdt = new ConformingDelaunayTriangulator(sites, tolerance);
*
* // optional
* cdt.setSplitPointFinder(splitPointFinder);
* cdt.setVertexFactory(vertexFactory);
*
* cdt.setConstraints(segments, new ArrayList(vertexMap.values()));
* cdt.formInitialDelaunay();
* cdt.enforceConstraints();
* subdiv = cdt.getSubdivision();
* </pre>
*
* @author David Skea
* @author Martin Davis
*/
public class ConformingDelaunayTriangulator
{
private static Envelope computeVertexEnvelope(Collection vertices) {
Envelope env = new Envelope();
for (Iterator i = vertices.iterator(); i.hasNext();) {
Vertex v = (Vertex) i.next();
env.expandToInclude(v.getCoordinate());
}
return env;
}
private List initialVertices; // List<Vertex>
private List segVertices; // List<Vertex>
// MD - using a Set doesn't seem to be much faster
// private Set segments = new HashSet();
private List segments = new ArrayList(); // List<Segment>
private QuadEdgeSubdivision subdiv = null;
private IncrementalDelaunayTriangulator incDel;
private Geometry convexHull;
private ConstraintSplitPointFinder splitFinder = new NonEncroachingSplitPointFinder();
private KdTree kdt = null;
private ConstraintVertexFactory vertexFactory = null;
// allPointsEnv expanded by a small buffer
private Envelope computeAreaEnv;
// records the last split point computed, for error reporting
private Coordinate splitPt = null;
private double tolerance; // defines if two sites are the same.
/**
* Creates a Conforming Delaunay Triangulation based on the given
* unconstrained initial vertices. The initial vertex set should not contain
* any vertices which appear in the constraint set.
*
* @param initialVertices
* a collection of {@link ConstraintVertex}
* @param tolerance
* the distance tolerance below which points are considered identical
*/
public ConformingDelaunayTriangulator(Collection initialVertices,
double tolerance) {
this.initialVertices = new ArrayList(initialVertices);
this.tolerance = tolerance;
kdt = new KdTree(tolerance);
}
/**
* Sets the constraints to be conformed to by the computed triangulation.
* The constraints must not contain duplicate segments (up to orientation).
* The unique set of vertices (as {@link ConstraintVertex}es)
* forming the constraints must also be supplied.
* Supplying it explicitly allows the ConstraintVertexes to be initialized
* appropriately (e.g. with external data), and avoids re-computing the unique set
* if it is already available.
*
* @param segments a list of the constraint {@link Segment}s
* @param segVertices the set of unique {@link ConstraintVertex}es referenced by the segments
*/
public void setConstraints(List segments, List segVertices) {
this.segments = segments;
this.segVertices = segVertices;
}
/**
* Sets the {@link ConstraintSplitPointFinder} to be
* used during constraint enforcement.
* Different splitting strategies may be appropriate
* for special situations.
*
* @param splitFinder the ConstraintSplitPointFinder to be used
*/
public void setSplitPointFinder(ConstraintSplitPointFinder splitFinder) {
this.splitFinder = splitFinder;
}
/**
* Gets the tolerance value used to construct the triangulation.
*
* @return a tolerance value
*/
public double getTolerance()
{
return tolerance;
}
/**
* Gets the <tt>ConstraintVertexFactory</tt> used to create new constraint vertices at split points.
*
* @return a new constraint vertex
*/
public ConstraintVertexFactory getVertexFactory() {
return vertexFactory;
}
/**
* Sets a custom {@link ConstraintVertexFactory} to be used
* to allow vertices carrying extra information to be created.
*
* @param vertexFactory the ConstraintVertexFactory to be used
*/
public void setVertexFactory(ConstraintVertexFactory vertexFactory) {
this.vertexFactory = vertexFactory;
}
/**
* Gets the {@link QuadEdgeSubdivision} which represents the triangulation.
*
* @return a subdivision
*/
public QuadEdgeSubdivision getSubdivision() {
return subdiv;
}
/**
* Gets the {@link KdTree} which contains the vertices of the triangulation.
*
* @return a KdTree
*/
public KdTree getKDT() {
return kdt;
}
/**
* Gets the sites (vertices) used to initialize the triangulation.
*
* @return a List of Vertex
*/
public List getInitialVertices() {
return initialVertices;
}
/**
* Gets the {@link Segment}s which represent the constraints.
*
* @return a collection of Segments
*/
public Collection getConstraintSegments() {
return segments;
}
/**
* Gets the convex hull of all the sites in the triangulation,
* including constraint vertices.
* Only valid after the constraints have been enforced.
*
* @return the convex hull of the sites
*/
public Geometry getConvexHull() {
return convexHull;
}
// ==================================================================
private void computeBoundingBox() {
Envelope vertexEnv = computeVertexEnvelope(initialVertices);
Envelope segEnv = computeVertexEnvelope(segVertices);
Envelope allPointsEnv = new Envelope(vertexEnv);
allPointsEnv.expandToInclude(segEnv);
double deltaX = allPointsEnv.getWidth() * 0.2;
double deltaY = allPointsEnv.getHeight() * 0.2;
double delta = Math.max(deltaX, deltaY);
computeAreaEnv = new Envelope(allPointsEnv);
computeAreaEnv.expandBy(delta);
}
private void computeConvexHull() {
GeometryFactory fact = new GeometryFactory();
Coordinate[] coords = getPointArray();
ConvexHull hull = new ConvexHull(coords, fact);
convexHull = hull.getConvexHull();
}
// /**
// * Adds the segments in the Convex Hull of all sites in the input data as
// linear constraints.
// * This is required if TIN Refinement is performed. The hull segments are
// flagged with a
// unique
// * data object to allow distinguishing them.
// *
// * @param convexHullSegmentData the data object to attach to each convex
// hull segment
// */
// private void addConvexHullToConstraints(Object convexHullSegmentData) {
// Coordinate[] coords = convexHull.getCoordinates();
// for (int i = 1; i < coords.length; i++) {
// Segment s = new Segment(coords[i - 1], coords[i], convexHullSegmentData);
// addConstraintIfUnique(s);
// }
// }
// private void addConstraintIfUnique(Segment r) {
// boolean exists = false;
// Iterator it = segments.iterator();
// Segment s = null;
// while (it.hasNext()) {
// s = (Segment) it.next();
// if (r.equalsTopo(s)) {
// exists = true;
// }
// }
// if (!exists) {
// segments.add((Object) r);
// }
// }
private Coordinate[] getPointArray() {
Coordinate[] pts = new Coordinate[initialVertices.size()
+ segVertices.size()];
int index = 0;
for (Iterator i = initialVertices.iterator(); i.hasNext();) {
Vertex v = (Vertex) i.next();
pts[index++] = v.getCoordinate();
}
for (Iterator i2 = segVertices.iterator(); i2.hasNext();) {
Vertex v = (Vertex) i2.next();
pts[index++] = v.getCoordinate();
}
return pts;
}
private ConstraintVertex createVertex(Coordinate p) {
ConstraintVertex v = null;
if (vertexFactory != null)
v = vertexFactory.createVertex(p, null);
else
v = new ConstraintVertex(p);
return v;
}
/**
* Creates a vertex on a constraint segment
*
* @param p the location of the vertex to create
* @param seg the constraint segment it lies on
* @return the new constraint vertex
*/
private ConstraintVertex createVertex(Coordinate p, Segment seg) {
ConstraintVertex v = null;
if (vertexFactory != null)
v = vertexFactory.createVertex(p, seg);
else
v = new ConstraintVertex(p);
v.setOnConstraint(true);
return v;
}
/**
* Inserts all sites in a collection
*
* @param vertices a collection of ConstraintVertex
*/
private void insertSites(Collection vertices) {
Debug.println("Adding sites: " + vertices.size());
for (Iterator i = vertices.iterator(); i.hasNext();) {
ConstraintVertex v = (ConstraintVertex) i.next();
insertSite(v);
}
}
private ConstraintVertex insertSite(ConstraintVertex v) {
KdNode kdnode = kdt.insert(v.getCoordinate(), v);
if (!kdnode.isRepeated()) {
incDel.insertSite(v);
} else {
ConstraintVertex snappedV = (ConstraintVertex) kdnode.getData();
snappedV.merge(v);
return snappedV;
// testing
// if ( v.isOnConstraint() && ! currV.isOnConstraint()) {
// System.out.println(v);
// }
}
return v;
}
/**
* Inserts a site into the triangulation, maintaining the conformal Delaunay property.
* This can be used to further refine the triangulation if required
* (e.g. to approximate the medial axis of the constraints,
* or to improve the grading of the triangulation).
*
* @param p the location of the site to insert
*/
public void insertSite(Coordinate p) {
insertSite(createVertex(p));
}
// ==================================================================
/**
* Computes the Delaunay triangulation of the initial sites.
*/
public void formInitialDelaunay() {
computeBoundingBox();
subdiv = new QuadEdgeSubdivision(computeAreaEnv, tolerance);
subdiv.setLocator(new LastFoundQuadEdgeLocator(subdiv));
incDel = new IncrementalDelaunayTriangulator(subdiv);
insertSites(initialVertices);
}
// ==================================================================
private final static int MAX_SPLIT_ITER = 99;
/**
* Enforces the supplied constraints into the triangulation.
*
* @throws ConstraintEnforcementException
* if the constraints cannot be enforced
*/
public void enforceConstraints() {
addConstraintVertices();
// if (true) return;
int count = 0;
int splits = 0;
do {
splits = enforceGabriel(segments);
count++;
Debug.println("Iter: " + count + " Splits: " + splits
+ " Current # segments = " + segments.size());
} while (splits > 0 && count < MAX_SPLIT_ITER);
if (count == MAX_SPLIT_ITER) {
Debug.println("ABORTED! Too many iterations while enforcing constraints");
if (!Debug.isDebugging())
throw new ConstraintEnforcementException(
"Too many splitting iterations while enforcing constraints. Last split point was at: ",
splitPt);
}
}
private void addConstraintVertices() {
computeConvexHull();
// insert constraint vertices as sites
insertSites(segVertices);
}
/*
* private List findMissingConstraints() { List missingSegs = new ArrayList();
* for (int i = 0; i < segments.size(); i++) { Segment s = (Segment)
* segments.get(i); QuadEdge q = subdiv.locate(s.getStart(), s.getEnd()); if
* (q == null) missingSegs.add(s); } return missingSegs; }
*/
private int enforceGabriel(Collection segsToInsert) {
List newSegments = new ArrayList();
int splits = 0;
List segsToRemove = new ArrayList();
/**
* On each iteration must always scan all constraint (sub)segments, since
* some constraints may be rebroken by Delaunay triangle flipping caused by
* insertion of another constraint. However, this process must converge
* eventually, with no splits remaining to find.
*/
for (Iterator i = segsToInsert.iterator(); i.hasNext();) {
Segment seg = (Segment) i.next();
// System.out.println(seg);
Coordinate encroachPt = findNonGabrielPoint(seg);
// no encroachment found - segment must already be in subdivision
if (encroachPt == null)
continue;
// compute split point
splitPt = splitFinder.findSplitPoint(seg, encroachPt);
ConstraintVertex splitVertex = createVertex(splitPt, seg);
// DebugFeature.addLineSegment(DEBUG_SEG_SPLIT, encroachPt, splitPt, "");
// Debug.println(WKTWriter.toLineString(encroachPt, splitPt));
/**
* Check whether the inserted point still equals the split pt. This will
* not be the case if the split pt was too close to an existing site. If
* the point was snapped, the triangulation will not respect the inserted
* constraint - this is a failure. This can be caused by:
* <ul>
* <li>An initial site that lies very close to a constraint segment The
* cure for this is to remove any initial sites which are close to
* constraint segments in a preprocessing phase.
* <li>A narrow constraint angle which causing repeated splitting until
* the split segments are too small. The cure for this is to either choose
* better split points or "guard" narrow angles by cracking the segments
* equidistant from the corner.
* </ul>
*/
ConstraintVertex insertedVertex = insertSite(splitVertex);
if (!insertedVertex.getCoordinate().equals2D(splitPt)) {
Debug.println("Split pt snapped to: " + insertedVertex);
// throw new ConstraintEnforcementException("Split point snapped to
// existing point
// (tolerance too large or constraint interior narrow angle?)",
// splitPt);
}
// split segment and record the new halves
Segment s1 = new Segment(seg.getStartX(), seg.getStartY(), seg
.getStartZ(), splitVertex.getX(), splitVertex.getY(), splitVertex
.getZ(), seg.getData());
Segment s2 = new Segment(splitVertex.getX(), splitVertex.getY(),
splitVertex.getZ(), seg.getEndX(), seg.getEndY(), seg.getEndZ(), seg
.getData());
newSegments.add(s1);
newSegments.add(s2);
segsToRemove.add(seg);
splits = splits + 1;
}
segsToInsert.removeAll(segsToRemove);
segsToInsert.addAll(newSegments);
return splits;
}
// public static final String DEBUG_SEG_SPLIT = "C:\\proj\\CWB\\test\\segSplit.jml";
/**
* Given a set of points stored in the kd-tree and a line segment defined by
* two points in this set, finds a {@link Coordinate} in the circumcircle of
* the line segment, if one exists. This is called the Gabriel point - if none
* exists then the segment is said to have the Gabriel condition. Uses the
* heuristic of finding the non-Gabriel point closest to the midpoint of the
* segment.
*
* @param p
* start of the line segment
* @param q
* end of the line segment
* @return a point which is non-Gabriel
* or null if no point is non-Gabriel
*/
private Coordinate findNonGabrielPoint(Segment seg) {
Coordinate p = seg.getStart();
Coordinate q = seg.getEnd();
// Find the mid point on the line and compute the radius of enclosing circle
Coordinate midPt = new Coordinate((p.x + q.x) / 2.0, (p.y + q.y) / 2.0);
double segRadius = p.distance(midPt);
// compute envelope of circumcircle
Envelope env = new Envelope(midPt);
env.expandBy(segRadius);
// Find all points in envelope
List result = kdt.query(env);
// For each point found, test if it falls strictly in the circle
// find closest point
Coordinate closestNonGabriel = null;
double minDist = Double.MAX_VALUE;
for (Iterator i = result.iterator(); i.hasNext();) {
KdNode nextNode = (KdNode) i.next();
Coordinate testPt = nextNode.getCoordinate();
// ignore segment endpoints
if (testPt.equals2D(p) || testPt.equals2D(q))
continue;
double testRadius = midPt.distance(testPt);
if (testRadius < segRadius) {
// double testDist = seg.distance(testPt);
double testDist = testRadius;
if (closestNonGabriel == null || testDist < minDist) {
closestNonGabriel = testPt;
minDist = testDist;
}
}
}
return closestNonGabriel;
}
}