org.apache.lucene.geo.Tessellator Maven / Gradle / Ivy
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.lucene.geo;
import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import org.apache.lucene.geo.GeoUtils.WindingOrder;
import org.apache.lucene.util.BitUtil;
import static org.apache.lucene.geo.GeoEncodingUtils.encodeLatitude;
import static org.apache.lucene.geo.GeoEncodingUtils.encodeLongitude;
import static org.apache.lucene.geo.GeoUtils.orient;
/**
* Computes a triangular mesh tessellation for a given polygon.
*
* This is inspired by mapbox's earcut algorithm (https://github.com/mapbox/earcut)
* which is a modification to FIST (https://www.cosy.sbg.ac.at/~held/projects/triang/triang.html)
* written by Martin Held, and ear clipping (https://www.geometrictools.com/Documentation/TriangulationByEarClipping.pdf)
* written by David Eberly.
*
* Notes:
*
* - Requires valid polygons:
*
* - No self intersections
*
- Holes may only touch at one vertex
*
- Polygon must have an area (e.g., no "line" boxes)
*
- sensitive to overflow (e.g, subatomic values such as E-200 can cause unexpected behavior)
*
*
*
* The code is a modified version of the javascript implementation provided by MapBox
* under the following license:
*
* ISC License
*
* Copyright (c) 2016, Mapbox
*
* Permission to use, copy, modify, and/or distribute this software for any purpose
* with or without fee is hereby granted, provided that the above copyright notice
* and this permission notice appear in all copies.
*
* THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES WITH'
* REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
* FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT,
* INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS
* OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
* TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
* THIS SOFTWARE.
*
* @lucene.internal
*/
final public class Tessellator {
// this is a dumb heuristic to control whether we cut over to sorted morton values
private static final int VERTEX_THRESHOLD = 80;
/** state of the tessellated split - avoids recursion */
private enum State {
INIT, CURE, SPLIT
}
// No Instance:
private Tessellator() {}
public static final List tessellate(final Polygon polygon) {
// Attempt to establish a doubly-linked list of the provided shell points (should be CCW, but this will correct);
// then filter instances of intersections.
Node outerNode = createDoublyLinkedList(polygon.getPolyLons(), polygon.getPolyLats(),polygon.getWindingOrder(), true,
0, WindingOrder.CW);
// If an outer node hasn't been detected, the shape is malformed. (must comply with OGC SFA specification)
if(outerNode == null) {
throw new IllegalArgumentException("Malformed shape detected in Tessellator!");
}
// Determine if the specified list of points contains holes
if (polygon.numHoles() > 0) {
// Eliminate the hole triangulation.
outerNode = eliminateHoles(polygon, outerNode);
}
// If the shape crosses VERTEX_THRESHOLD, use z-order curve hashing:
final boolean mortonOptimized;
{
int threshold = VERTEX_THRESHOLD - polygon.numPoints();
for (int i = 0; threshold >= 0 && i < polygon.numHoles(); ++i) {
threshold -= polygon.getHole(i).numPoints();
}
// Link polygon nodes in Z-Order
mortonOptimized = threshold < 0;
if (mortonOptimized == true) {
sortByMorton(outerNode);
}
}
// Calculate the tessellation using the doubly LinkedList.
List result = earcutLinkedList(polygon, outerNode, new ArrayList<>(), State.INIT, mortonOptimized);
if (result.size() == 0) {
throw new IllegalArgumentException("Unable to Tessellate shape [" + polygon + "]. Possible malformed shape detected.");
}
return result;
}
public static final List tessellate(final XYPolygon polygon) {
// Attempt to establish a doubly-linked list of the provided shell points (should be CCW, but this will correct);
// then filter instances of intersections.0
Node outerNode = createDoublyLinkedList(XYEncodingUtils.floatArrayToDoubleArray(polygon.getPolyX()), XYEncodingUtils.floatArrayToDoubleArray(polygon.getPolyY()),
polygon.getWindingOrder(), false, 0, WindingOrder.CW);
// If an outer node hasn't been detected, the shape is malformed. (must comply with OGC SFA specification)
if(outerNode == null) {
throw new IllegalArgumentException("Malformed shape detected in Tessellator!");
}
// Determine if the specified list of points contains holes
if (polygon.numHoles() > 0) {
// Eliminate the hole triangulation.
outerNode = eliminateHoles(polygon, outerNode);
}
// If the shape crosses VERTEX_THRESHOLD, use z-order curve hashing:
final boolean mortonOptimized;
{
int threshold = VERTEX_THRESHOLD - polygon.numPoints();
for (int i = 0; threshold >= 0 && i < polygon.numHoles(); ++i) {
threshold -= polygon.getHole(i).numPoints();
}
// Link polygon nodes in Z-Order
mortonOptimized = threshold < 0;
if (mortonOptimized == true) {
sortByMorton(outerNode);
}
}
// Calculate the tessellation using the doubly LinkedList.
List result = earcutLinkedList(polygon, outerNode, new ArrayList<>(), State.INIT, mortonOptimized);
if (result.size() == 0) {
throw new IllegalArgumentException("Unable to Tessellate shape [" + polygon + "]. Possible malformed shape detected.");
}
return result;
}
/** Creates a circular doubly linked list using polygon points. The order is governed by the specified winding order */
private static final Node createDoublyLinkedList(final double[] x, final double[] y, final WindingOrder polyWindingOrder,
boolean isGeo, int startIndex, final WindingOrder windingOrder) {
Node lastNode = null;
// Link points into the circular doubly-linked list in the specified winding order
if (windingOrder == polyWindingOrder) {
for (int i = 0; i < x.length; ++i) {
lastNode = insertNode(x, y, startIndex++, i, lastNode, isGeo);
}
} else {
for (int i = x.length - 1; i >= 0; --i) {
lastNode = insertNode(x, y, startIndex++, i, lastNode, isGeo);
}
}
// if first and last node are the same then remove the end node and set lastNode to the start
if (lastNode != null && isVertexEquals(lastNode, lastNode.next)) {
removeNode(lastNode, true);
lastNode = lastNode.next;
}
// Return the last node in the Doubly-Linked List
return filterPoints(lastNode, null);
}
private static final Node eliminateHoles(final XYPolygon polygon, Node outerNode) {
// Define a list to hole a reference to each filtered hole list.
final List holeList = new ArrayList<>();
// keep a reference to the hole
final Map holeListPolygons = new HashMap<>();
// Iterate through each array of hole vertices.
XYPolygon[] holes = polygon.getHoles();
int nodeIndex = polygon.numPoints() ;
for(int i = 0; i < polygon.numHoles(); ++i) {
// create the doubly-linked hole list
Node list = createDoublyLinkedList(XYEncodingUtils.floatArrayToDoubleArray(holes[i].getPolyX()),
XYEncodingUtils.floatArrayToDoubleArray(holes[i].getPolyY()), holes[i].getWindingOrder(), false, nodeIndex, WindingOrder.CCW);
// Determine if the resulting hole polygon was successful.
if(list != null) {
// Add the leftmost vertex of the hole.
Node leftMost = fetchLeftmost(list);
holeList.add(leftMost);
holeListPolygons.put(leftMost, holes[i]);
}
nodeIndex += holes[i].numPoints();
}
return eliminateHoles(holeList, holeListPolygons, outerNode);
}
/** Links every hole into the outer loop, producing a single-ring polygon without holes. **/
private static final Node eliminateHoles(final Polygon polygon, Node outerNode) {
// Define a list to hole a reference to each filtered hole list.
final List holeList = new ArrayList<>();
// keep a reference to the hole
final Map holeListPolygons = new HashMap<>();
// Iterate through each array of hole vertices.
Polygon[] holes = polygon.getHoles();
int nodeIndex = polygon.numPoints();
for(int i = 0; i < polygon.numHoles(); ++i) {
// create the doubly-linked hole list
Node list = createDoublyLinkedList(holes[i].getPolyLons(), holes[i].getPolyLats(), holes[i].getWindingOrder(), true, nodeIndex, WindingOrder.CCW);
if (list == list.next) {
throw new IllegalArgumentException("Points are all coplanar in hole: " + holes[i]);
}
// Determine if the resulting hole polygon was successful.
if(list != null) {
// Add the leftmost vertex of the hole.
Node leftMost = fetchLeftmost(list);
holeList.add(leftMost);
holeListPolygons.put(leftMost, holes[i]);
}
nodeIndex += holes[i].numPoints();
}
return eliminateHoles(holeList, holeListPolygons, outerNode);
}
private static final Node eliminateHoles(List holeList, final Map holeListPolygons, Node outerNode) {
// Sort the hole vertices by x coordinate
holeList.sort((Node pNodeA, Node pNodeB) ->
{
double diff = pNodeA.getX() - pNodeB.getX();
if (diff == 0) {
diff = pNodeA.getY() - pNodeB.getY();
if (diff == 0) {
//same hole node
double a = Math.min(pNodeA.previous.getY(), pNodeA.next.getY());
double b = Math.min(pNodeB.previous.getY(), pNodeB.next.getY());
diff = a - b;
}
}
return diff < 0 ? -1 : diff > 0 ? 1 : 0;
});
// Process holes from left to right.
for(int i = 0; i < holeList.size(); ++i) {
// Eliminate hole triangles from the result set
final Node holeNode = holeList.get(i);
double holeMinX, holeMaxX, holeMinY, holeMaxY;
Object h = holeListPolygons.get(holeNode);
if (h instanceof Polygon) {
Polygon holePoly = (Polygon)h;
holeMinX = holePoly.minLon;
holeMaxX = holePoly.maxLon;
holeMinY = holePoly.minLat;
holeMaxY = holePoly.maxLat;
} else {
XYPolygon holePoly = (XYPolygon)h;
holeMinX = holePoly.minX;
holeMaxX = holePoly.maxX;
holeMinY = holePoly.minY;
holeMaxY = holePoly.maxY;
}
eliminateHole(holeNode, outerNode, holeMinX, holeMaxX, holeMinY, holeMaxY);
// Filter the new polygon.
outerNode = filterPoints(outerNode, outerNode.next);
}
// Return a pointer to the list.
return outerNode;
}
/** Finds a bridge between vertices that connects a hole with an outer ring, and links it */
private static final void eliminateHole(final Node holeNode, Node outerNode, double holeMinX, double holeMaxX, double holeMinY, double holeMaxY) {
// Attempt to find a common point between the HoleNode and OuterNode.
Node next = outerNode;
do {
if (Rectangle.containsPoint(next.getY(), next.getX(), holeMinY, holeMaxY, holeMinX, holeMaxX)) {
Node sharedVertex = getSharedVertex(holeNode, next);
if (sharedVertex != null) {
// Split the resulting polygon.
Node node = splitPolygon(next, sharedVertex, true);
// Filter the split nodes.
filterPoints(node, node.next);
return;
}
}
next = next.next;
} while (next != outerNode);
// Attempt to find a logical bridge between the HoleNode and OuterNode.
outerNode = fetchHoleBridge(holeNode, outerNode);
// Determine whether a hole bridge could be fetched.
if(outerNode != null) {
// compute if the bridge overlaps with a polygon edge.
boolean fromPolygon = isPointInLine(outerNode, outerNode.next, holeNode) || isPointInLine(holeNode, holeNode.next, outerNode);
// Split the resulting polygon.
Node node = splitPolygon(outerNode, holeNode, fromPolygon);
// Filter the split nodes.
filterPoints(node, node.next);
}
}
/**
* David Eberly's algorithm for finding a bridge between a hole and outer polygon
*
* see: http://www.geometrictools.com/Documentation/TriangulationByEarClipping.pdf
**/
private static final Node fetchHoleBridge(final Node holeNode, final Node outerNode) {
Node p = outerNode;
double qx = Double.NEGATIVE_INFINITY;
final double hx = holeNode.getX();
final double hy = holeNode.getY();
Node connection = null;
// 1. find a segment intersected by a ray from the hole's leftmost point to the left;
// segment's endpoint with lesser x will be potential connection point
{
do {
if (hy <= p.getY() && hy >= p.next.getY() && p.next.getY() != p.getY()) {
final double x = p.getX() + (hy - p.getY()) * (p.next.getX() - p.getX()) / (p.next.getY() - p.getY());
if (x <= hx && x > qx) {
qx = x;
if (x == hx) {
if (hy == p.getY()) return p;
if (hy == p.next.getY()) return p.next;
}
connection = p.getX() < p.next.getX() ? p : p.next;
}
}
p = p.next;
} while (p != outerNode);
}
if (connection == null) {
return null;
} else if (hx == qx) {
return connection.previous;
}
// 2. look for points inside the triangle of hole point, segment intersection, and endpoint
// its a valid connection iff there are no points found;
// otherwise choose the point of the minimum angle with the ray as the connection point
Node stop = connection;
final double mx = connection.getX();
final double my = connection.getY();
double tanMin = Double.POSITIVE_INFINITY;
double tan;
p = connection.next;
{
while (p != stop) {
if (hx >= p.getX() && p.getX() >= mx && hx != p.getX()
&& pointInEar(p.getX(), p.getY(), hy < my ? hx : qx, hy, mx, my, hy < my ? qx : hx, hy)) {
tan = Math.abs(hy - p.getY()) / (hx - p.getX()); // tangential
if (isVertexEquals(p, connection) && isLocallyInside(p, holeNode)) {
// make sure we are not crossing the polygon. This might happen when several holes have a bridge to the same polygon vertex
// and this vertex has different vertex.
boolean crosses = GeoUtils.lineCrossesLine(p.getX(), p.getY(), holeNode.getX(), holeNode.getY(),
connection.next.getX(), connection.next.getY(), connection.previous.getX(), connection.previous.getY());
if (crosses == false) {
connection = p;
tanMin = tan;
}
} else if ((tan < tanMin || (tan == tanMin && p.getX() > connection.getX())) && isLocallyInside(p, holeNode)) {
connection = p;
tanMin = tan;
}
}
p = p.next;
}
}
return connection;
}
/** Check if the provided vertex is in the polygon and return it **/
private static Node getSharedVertex(final Node polygon, final Node vertex) {
Node next = polygon;
do {
if (isVertexEquals(next, vertex)) {
// make sure we are not crossing the polygon. This might happen when several holes share the same polygon vertex.
boolean crosses = GeoUtils.lineCrossesLine(next.previous.getX(), next.previous.getY(), vertex.next.getX(), vertex.next.getY(),
next.next.getX(), next.next.getY(), vertex.previous.getX(), vertex.previous.getY());
if (crosses == false) {
return next;
}
}
next = next.next;
} while(next != polygon);
return null;
}
/** Finds the left-most hole of a polygon ring. **/
private static final Node fetchLeftmost(final Node start) {
Node node = start;
Node leftMost = start;
do {
// Determine if the current node possesses a lesser X position.
if (node.getX() < leftMost.getX() || (node.getX() == leftMost.getX() && node.getY() < leftMost.getY())) {
// Maintain a reference to this Node.
leftMost = node;
}
// Progress the search to the next node in the doubly-linked list.
node = node.next;
} while (node != start);
// Return the node with the smallest X value.
return leftMost;
}
/** Main ear slicing loop which triangulates the vertices of a polygon, provided as a doubly-linked list. **/
private static final List earcutLinkedList(Object polygon, Node currEar, final List tessellation,
State state, final boolean mortonOptimized) {
earcut : do {
if (currEar == null || currEar.previous == currEar.next) {
return tessellation;
}
Node stop = currEar;
Node prevNode;
Node nextNode;
// Iteratively slice ears
do {
prevNode = currEar.previous;
nextNode = currEar.next;
// Determine whether the current triangle must be cut off.
final boolean isReflex = area(prevNode.getX(), prevNode.getY(), currEar.getX(), currEar.getY(), nextNode.getX(), nextNode.getY()) >= 0;
if (isReflex == false && isEar(currEar, mortonOptimized) == true) {
// Compute if edges belong to the polygon
boolean abFromPolygon = prevNode.isNextEdgeFromPolygon;
boolean bcFromPolygon = currEar.isNextEdgeFromPolygon;
boolean caFromPolygon = isEdgeFromPolygon(prevNode, nextNode, mortonOptimized);
// Return the triangulated data
tessellation.add(new Triangle(prevNode, abFromPolygon, currEar, bcFromPolygon, nextNode, caFromPolygon));
// Remove the ear node.
removeNode(currEar, caFromPolygon);
// Skipping to the next node leaves fewer slither triangles.
currEar = nextNode.next;
stop = nextNode.next;
continue;
}
currEar = nextNode;
// If the whole polygon has been iterated over and no more ears can be found.
if (currEar == stop) {
switch (state) {
case INIT:
// try filtering points and slicing again
currEar = filterPoints(currEar, null);
state = State.CURE;
continue earcut;
case CURE:
// if this didn't work, try curing all small self-intersections locally
currEar = cureLocalIntersections(currEar, tessellation, mortonOptimized);
state = State.SPLIT;
continue earcut;
case SPLIT:
// as a last resort, try splitting the remaining polygon into two
if (splitEarcut(polygon, currEar, tessellation, mortonOptimized) == false) {
//we could not process all points. Tessellation failed
throw new IllegalArgumentException("Unable to Tessellate shape [" + polygon + "]. Possible malformed shape detected.");
}
break;
}
break;
}
} while (currEar.previous != currEar.next);
break;
} while (true);
// Return the calculated tessellation
return tessellation;
}
/** Determines whether a polygon node forms a valid ear with adjacent nodes. **/
private static final boolean isEar(final Node ear, final boolean mortonOptimized) {
if (mortonOptimized == true) {
return mortonIsEar(ear);
}
// make sure there aren't other points inside the potential ear
Node node = ear.next.next;
while (node != ear.previous) {
if (pointInEar(node.getX(), node.getY(), ear.previous.getX(), ear.previous.getY(), ear.getX(), ear.getY(), ear.next.getX(), ear.next.getY())
&& area(node.previous.getX(), node.previous.getY(), node.getX(), node.getY(),
node.next.getX(), node.next.getY()) >= 0) {
return false;
}
node = node.next;
}
return true;
}
/** Uses morton code for speed to determine whether or a polygon node forms a valid ear w/ adjacent nodes */
private static final boolean mortonIsEar(final Node ear) {
// triangle bbox (flip the bits so negative encoded values are < positive encoded values)
int minTX = StrictMath.min(StrictMath.min(ear.previous.x, ear.x), ear.next.x) ^ 0x80000000;
int minTY = StrictMath.min(StrictMath.min(ear.previous.y, ear.y), ear.next.y) ^ 0x80000000;
int maxTX = StrictMath.max(StrictMath.max(ear.previous.x, ear.x), ear.next.x) ^ 0x80000000;
int maxTY = StrictMath.max(StrictMath.max(ear.previous.y, ear.y), ear.next.y) ^ 0x80000000;
// z-order range for the current triangle bbox;
long minZ = BitUtil.interleave(minTX, minTY);
long maxZ = BitUtil.interleave(maxTX, maxTY);
// now make sure we don't have other points inside the potential ear;
// look for points inside the triangle in both directions
Node p = ear.previousZ;
Node n = ear.nextZ;
while (p != null && Long.compareUnsigned(p.morton, minZ) >= 0
&& n != null && Long.compareUnsigned(n.morton, maxZ) <= 0) {
if (p.idx != ear.previous.idx && p.idx != ear.next.idx &&
pointInEar(p.getX(), p.getY(), ear.previous.getX(), ear.previous.getY(), ear.getX(), ear.getY(), ear.next.getX(), ear.next.getY()) &&
area(p.previous.getX(), p.previous.getY(), p.getX(), p.getY(), p.next.getX(), p.next.getY()) >= 0) return false;
p = p.previousZ;
if (n.idx != ear.previous.idx && n.idx != ear.next.idx &&
pointInEar(n.getX(), n.getY(), ear.previous.getX(), ear.previous.getY(), ear.getX(), ear.getY(), ear.next.getX(), ear.next.getY()) &&
area(n.previous.getX(), n.previous.getY(), n.getX(), n.getY(), n.next.getX(), n.next.getY()) >= 0) return false;
n = n.nextZ;
}
// first look for points inside the triangle in decreasing z-order
while (p != null && Long.compareUnsigned(p.morton, minZ) >= 0) {
if (p.idx != ear.previous.idx && p.idx != ear.next.idx
&& pointInEar(p.getX(), p.getY(), ear.previous.getX(), ear.previous.getY(), ear.getX(), ear.getY(), ear.next.getX(), ear.next.getY())
&& area(p.previous.getX(), p.previous.getY(), p.getX(), p.getY(), p.next.getX(), p.next.getY()) >= 0) {
return false;
}
p = p.previousZ;
}
// then look for points in increasing z-order
while (n != null &&
Long.compareUnsigned(n.morton, maxZ) <= 0) {
if (n.idx != ear.previous.idx && n.idx != ear.next.idx
&& pointInEar(n.getX(), n.getY(), ear.previous.getX(), ear.previous.getY(), ear.getX(), ear.getY(), ear.next.getX(), ear.next.getY())
&& area(n.previous.getX(), n.previous.getY(), n.getX(), n.getY(), n.next.getX(), n.next.getY()) >= 0) {
return false;
}
n = n.nextZ;
}
return true;
}
/** Iterate through all polygon nodes and remove small local self-intersections **/
private static final Node cureLocalIntersections(Node startNode, final List tessellation, final boolean mortonOptimized) {
Node node = startNode;
Node nextNode;
do {
nextNode = node.next;
Node a = node.previous;
Node b = nextNode.next;
// a self-intersection where edge (v[i-1],v[i]) intersects (v[i+1],v[i+2])
if (isVertexEquals(a, b) == false
&& isIntersectingPolygon(a, a.getX(), a.getY(), b.getX(), b.getY()) == false
&& linesIntersect(a.getX(), a.getY(), node.getX(), node.getY(), nextNode.getX(), nextNode.getY(), b.getX(), b.getY())
&& isLocallyInside(a, b) && isLocallyInside(b, a)) {
// compute edges from polygon
boolean abFromPolygon = (a.next == node) ? a.isNextEdgeFromPolygon : isEdgeFromPolygon(a, node, mortonOptimized);
boolean bcFromPolygon = (node.next == b) ? node.isNextEdgeFromPolygon : isEdgeFromPolygon(node, b, mortonOptimized);
boolean caFromPolygon = (b.next == a) ? b.isNextEdgeFromPolygon : isEdgeFromPolygon(a, b, mortonOptimized);
tessellation.add(new Triangle(a, abFromPolygon, node, bcFromPolygon, b, caFromPolygon));
// Return the triangulated vertices to the tessellation
tessellation.add(new Triangle(a, abFromPolygon, node, bcFromPolygon, b, caFromPolygon));
// remove two nodes involved
removeNode(node, caFromPolygon);
removeNode(node.next, caFromPolygon);
node = startNode = b;
}
node = node.next;
} while (node != startNode);
return node;
}
/** Attempt to split a polygon and independently triangulate each side. Return true if the polygon was splitted **/
private static final boolean splitEarcut(final Object polygon, final Node start, final List tessellation, final boolean mortonOptimized) {
// Search for a valid diagonal that divides the polygon into two.
Node searchNode = start;
Node nextNode;
do {
nextNode = searchNode.next;
Node diagonal = nextNode.next;
while (diagonal != searchNode.previous) {
if(searchNode.idx != diagonal.idx && isValidDiagonal(searchNode, diagonal)) {
// Split the polygon into two at the point of the diagonal
Node splitNode = splitPolygon(searchNode, diagonal, isEdgeFromPolygon(searchNode, diagonal, mortonOptimized));
// Filter the resulting polygon.
searchNode = filterPoints(searchNode, searchNode.next);
splitNode = filterPoints(splitNode, splitNode.next);
// Attempt to earcut both of the resulting polygons
if (mortonOptimized) {
sortByMortonWithReset(searchNode);
sortByMortonWithReset(splitNode);
}
earcutLinkedList(polygon, searchNode, tessellation, State.INIT, mortonOptimized);
earcutLinkedList(polygon, splitNode, tessellation, State.INIT, mortonOptimized);
// Finish the iterative search
return true;
}
diagonal = diagonal.next;
}
searchNode = searchNode.next;
} while (searchNode != start);
return false;
}
/** Computes if edge defined by a and b overlaps with a polygon edge **/
private static boolean isEdgeFromPolygon(final Node a, final Node b, final boolean isMorton) {
if (isMorton) {
return isMortonEdgeFromPolygon(a, b);
}
Node next = a;
do {
if (isPointInLine(next, next.next, a) && isPointInLine(next, next.next, b)) {
return next.isNextEdgeFromPolygon;
}
if (isPointInLine(next, next.previous, a) && isPointInLine(next, next.previous, b)) {
return next.previous.isNextEdgeFromPolygon;
}
next = next.next;
} while(next != a);
return false;
}
/** Uses morton code for speed to determine whether or not and edge defined by a and b overlaps with a polygon edge */
private static final boolean isMortonEdgeFromPolygon(final Node a, final Node b) {
// edge bbox (flip the bits so negative encoded values are < positive encoded values)
final int minTX = StrictMath.min(a.x, b.x) ^ 0x80000000;
final int minTY = StrictMath.min(a.y, b.y) ^ 0x80000000;
final int maxTX = StrictMath.max(a.x, b.x) ^ 0x80000000;
final int maxTY = StrictMath.max(a.y, b.y) ^ 0x80000000;
// z-order range for the current edge;
final long minZ = BitUtil.interleave(minTX, minTY);
final long maxZ = BitUtil.interleave(maxTX, maxTY);
// now make sure we don't have other points inside the potential ear;
// look for points inside edge in both directions
Node p = a.previousZ;
Node n = a.nextZ;
while (p != null && Long.compareUnsigned(p.morton, minZ) >= 0
&& n != null && Long.compareUnsigned(n.morton, maxZ) <= 0) {
if (isPointInLine(p, p.next, a) && isPointInLine(p, p.next, b)) {
return p.isNextEdgeFromPolygon;
}
if (isPointInLine(p, p.previous, a) && isPointInLine(p, p.previous, b)) {
return p.previous.isNextEdgeFromPolygon;
}
p = p.previousZ;
if (isPointInLine(n, n.next, a) && isPointInLine(n, n.next, b)) {
return n.isNextEdgeFromPolygon;
}
if (isPointInLine(n, n.previous, a) && isPointInLine(n, n.previous, b)) {
return n.previous.isNextEdgeFromPolygon;
}
n = n.nextZ;
}
// first look for points inside the edge in decreasing z-order
while (p != null && Long.compareUnsigned(p.morton, minZ) >= 0) {
if (isPointInLine(p, p.next, a) && isPointInLine(p, p.next, b)) {
return p.isNextEdgeFromPolygon;
}
if (isPointInLine(p, p.previous, a) && isPointInLine(p, p.previous, b)) {
return p.previous.isNextEdgeFromPolygon;
}
p = p.previousZ;
}
// then look for points in increasing z-order
while (n != null &&
Long.compareUnsigned(n.morton, maxZ) <= 0) {
if (isPointInLine(n, n.next, a) && isPointInLine(n, n.next, b)) {
return n.isNextEdgeFromPolygon;
}
if (isPointInLine(n, n.previous, a) && isPointInLine(n, n.previous, b)) {
return n.previous.isNextEdgeFromPolygon;
}
n = n.nextZ;
}
return false;
}
private static boolean isPointInLine(final Node a, final Node b, final Node point) {
return isPointInLine(a, b, point.getX(), point.getY());
}
private static boolean isPointInLine(final Node a, final Node b, final double lon, final double lat) {
final double dxc = lon - a.getX();
final double dyc = lat - a.getY();
final double dxl = b.getX() - a.getX();
final double dyl = b.getY() - a.getY();
if (dxc * dyl - dyc * dxl == 0) {
if (Math.abs(dxl) >= Math.abs(dyl)) {
return dxl > 0 ?
a.getX() <= lon && lon <= b.getX() :
b.getX() <= lon && lon <= a.getX();
} else {
return dyl > 0 ?
a.getY() <= lat && lat <= b.getY() :
b.getY() <= lat && lat <= a.getY();
}
}
return false;
}
/** Links two polygon vertices using a bridge. **/
private static final Node splitPolygon(final Node a, final Node b, boolean edgeFromPolygon) {
final Node a2 = new Node(a);
final Node b2 = new Node(b);
final Node an = a.next;
final Node bp = b.previous;
a.next = b;
a.isNextEdgeFromPolygon = edgeFromPolygon;
a.nextZ = b;
b.previous = a;
b.previousZ = a;
a2.next = an;
a2.nextZ = an;
an.previous = a2;
an.previousZ = a2;
b2.next = a2;
b2.isNextEdgeFromPolygon = edgeFromPolygon;
b2.nextZ = a2;
a2.previous = b2;
a2.previousZ = b2;
bp.next = b2;
bp.nextZ = b2;
return b2;
}
/** Determines whether a diagonal between two polygon nodes lies within a polygon interior. (This determines the validity of the ray.) **/
private static final boolean isValidDiagonal(final Node a, final Node b) {
if (isVertexEquals(a, b)) {
//If points are equal then use it if they are valid polygons
return isCWPolygon(a, b);
}
return a.next.idx != b.idx && a.previous.idx != b.idx
&& isIntersectingPolygon(a, a.getX(), a.getY(), b.getX(), b.getY()) == false
&& isLocallyInside(a, b) && isLocallyInside(b, a)
&& middleInsert(a, a.getX(), a.getY(), b.getX(), b.getY());
}
/** Determine whether the polygon defined between node start and node end is CW */
private static boolean isCWPolygon(final Node start, final Node end) {
Node next = start;
double windingSum = 0;
do {
// compute signed area
windingSum += area(next.getX(), next.getY(), next.next.getX(), next.next.getY(), end.getX(), end.getY());
next = next.next;
} while (next.next != end);
//The polygon must be CW
return (windingSum < 0) ? true : false;
}
private static final boolean isLocallyInside(final Node a, final Node b) {
double area = area(a.previous.getX(), a.previous.getY(), a.getX(), a.getY(), a.next.getX(), a.next.getY());
if (area == 0) {
// parallel
return false;
} else if (area < 0) {
// if a is cw
return area(a.getX(), a.getY(), b.getX(), b.getY(), a.next.getX(), a.next.getY()) >= 0
&& area(a.getX(), a.getY(), a.previous.getX(), a.previous.getY(), b.getX(), b.getY()) >= 0;
} else {
// ccw
return area(a.getX(), a.getY(), b.getX(), b.getY(), a.previous.getX(), a.previous.getY()) < 0
|| area(a.getX(), a.getY(), a.next.getX(), a.next.getY(), b.getX(), b.getY()) < 0;
}
}
/** Determine whether the middle point of a polygon diagonal is contained within the polygon */
private static final boolean middleInsert(final Node start, final double x0, final double y0,
final double x1, final double y1) {
Node node = start;
Node nextNode;
boolean lIsInside = false;
final double lDx = (x0 + x1) / 2.0f;
final double lDy = (y0 + y1) / 2.0f;
do {
nextNode = node.next;
if (node.getY() > lDy != nextNode.getY() > lDy &&
lDx < (nextNode.getX() - node.getX()) * (lDy - node.getY()) / (nextNode.getY() - node.getY()) + node.getX()) {
lIsInside = !lIsInside;
}
node = node.next;
} while (node != start);
return lIsInside;
}
/** Determines if the diagonal of a polygon is intersecting with any polygon elements. **/
private static final boolean isIntersectingPolygon(final Node start, final double x0, final double y0,
final double x1, final double y1) {
Node node = start;
Node nextNode;
do {
nextNode = node.next;
if(isVertexEquals(node, x0, y0) == false && isVertexEquals(node, x1, y1) == false) {
if (linesIntersect(node.getX(), node.getY(), nextNode.getX(), nextNode.getY(), x0, y0, x1, y1)) {
return true;
}
}
node = nextNode;
} while (node != start);
return false;
}
/** Determines whether two line segments intersect. **/
public static final boolean linesIntersect(final double aX0, final double aY0, final double aX1, final double aY1,
final double bX0, final double bY0, final double bX1, final double bY1) {
return (area(aX0, aY0, aX1, aY1, bX0, bY0) > 0) != (area(aX0, aY0, aX1, aY1, bX1, bY1) > 0)
&& (area(bX0, bY0, bX1, bY1, aX0, aY0) > 0) != (area(bX0, bY0, bX1, bY1, aX1, aY1) > 0);
}
/** Interlinks polygon nodes in Z-Order. It reset the values on the z values**/
private static final void sortByMortonWithReset(Node start) {
Node next = start;
do {
next.previousZ = next.previous;
next.nextZ = next.next;
next = next.next;
} while (next != start);
sortByMorton(start);
}
/** Interlinks polygon nodes in Z-Order. **/
private static final void sortByMorton(Node start) {
start.previousZ.nextZ = null;
start.previousZ = null;
// Sort the generated ring using Z ordering.
tathamSort(start);
}
/**
* Simon Tatham's doubly-linked list O(n log n) mergesort
* see: http://www.chiark.greenend.org.uk/~sgtatham/algorithms/listsort.html
**/
private static final void tathamSort(Node list) {
Node p, q, e, tail;
int i, numMerges, pSize, qSize;
int inSize = 1;
if (list == null) {
return;
}
do {
p = list;
list = null;
tail = null;
// count number of merges in this pass
numMerges = 0;
while(p != null) {
++numMerges;
// step 'insize' places along from p
q = p;
for (i = 0, pSize = 0; i < inSize && q != null; ++i, ++pSize, q = q.nextZ);
// if q hasn't fallen off end, we have two lists to merge
qSize = inSize;
// now we have two lists; merge
while (pSize > 0 || (qSize > 0 && q != null)) {
if (pSize != 0 && (qSize == 0 || q == null || Long.compareUnsigned(p.morton, q.morton) <= 0)) {
e = p;
p = p.nextZ;
--pSize;
} else {
e = q;
q = q.nextZ;
--qSize;
}
if (tail != null) {
tail.nextZ = e;
} else {
list = e;
}
// maintain reverse pointers
e.previousZ = tail;
tail = e;
}
// now p has stepped 'insize' places along, and q has too
p = q;
}
tail.nextZ = null;
inSize *= 2;
} while (numMerges > 1);
}
/** Eliminate colinear/duplicate points from the doubly linked list */
private static final Node filterPoints(final Node start, Node end) {
if (start == null) {
return start;
}
if(end == null) {
end = start;
}
Node node = start;
Node nextNode;
Node prevNode;
boolean continueIteration;
do {
continueIteration = false;
nextNode = node.next;
prevNode = node.previous;
// we can filter points when:
if (isVertexEquals(node, nextNode) || // 1. they are the same,
isVertexEquals(prevNode, nextNode) || // 2.- each one starts and ends in each other
(prevNode.isNextEdgeFromPolygon == node.isNextEdgeFromPolygon && // 3.- they are co-linear and both edges have the same value in .isNextEdgeFromPolygon
area(prevNode.getX(), prevNode.getY(), node.getX(), node.getY(), nextNode.getX(), nextNode.getY()) == 0)) {
// Remove the node
boolean nextEdgeFromPol = prevNode.isNextEdgeFromPolygon != node.isNextEdgeFromPolygon ? true : prevNode.isNextEdgeFromPolygon;
removeNode(node, nextEdgeFromPol);
node = end = prevNode;
if (node == nextNode) {
break;
}
continueIteration = true;
} else {
node = nextNode;
}
} while (continueIteration || node != end);
return end;
}
/** Creates a node and optionally links it with a previous node in a circular doubly-linked list */
private static final Node insertNode(final double[] x, final double[] y, int index, int vertexIndex, final Node lastNode, boolean isGeo) {
final Node node = new Node(x, y, index, vertexIndex, isGeo);
if(lastNode == null) {
node.previous = node;
node.previousZ = node;
node.next = node;
node.nextZ = node;
} else {
node.next = lastNode.next;
node.nextZ = lastNode.next;
node.previous = lastNode;
node.previousZ = lastNode;
lastNode.next.previous = node;
lastNode.nextZ.previousZ = node;
lastNode.next = node;
lastNode.nextZ = node;
}
return node;
}
/** Removes a node from the doubly linked list */
private static final void removeNode(Node node, boolean edgeFromPolygon) {
node.next.previous = node.previous;
node.previous.next = node.next;
node.previous.isNextEdgeFromPolygon = edgeFromPolygon;
if (node.previousZ != null) {
node.previousZ.nextZ = node.nextZ;
}
if (node.nextZ != null) {
node.nextZ.previousZ = node.previousZ;
}
}
/** Determines if two point vertices are equal. **/
private static final boolean isVertexEquals(final Node a, final Node b) {
return isVertexEquals(a, b.getX(), b.getY());
}
/** Determines if two point vertices are equal. **/
private static final boolean isVertexEquals(final Node a, final double x, final double y) {
return a.getX() == x && a.getY() == y;
}
/** Compute signed area of triangle */
private static double area(final double aX, final double aY, final double bX, final double bY,
final double cX, final double cY) {
return (bY - aY) * (cX - bX) - (bX - aX) * (cY - bY);
}
/** Compute whether point is in a candidate ear */
private static boolean pointInEar(final double x, final double y, final double ax, final double ay,
final double bx, final double by, final double cx, final double cy) {
return (cx - x) * (ay - y) - (ax - x) * (cy - y) >= 0 &&
(ax - x) * (by - y) - (bx - x) * (ay - y) >= 0 &&
(bx - x) * (cy - y) - (cx - x) * (by - y) >= 0;
}
/** compute whether the given x, y point is in a triangle; uses the winding order method */
public static boolean pointInTriangle (double x, double y, double ax, double ay, double bx, double by, double cx, double cy) {
double minX = StrictMath.min(ax, StrictMath.min(bx, cx));
double minY = StrictMath.min(ay, StrictMath.min(by, cy));
double maxX = StrictMath.max(ax, StrictMath.max(bx, cx));
double maxY = StrictMath.max(ay, StrictMath.max(by, cy));
//check the bounding box because if the triangle is degenerated, e.g points and lines, we need to filter out
//coplanar points that are not part of the triangle.
if (x >= minX && x <= maxX && y >= minY && y <= maxY ) {
int a = orient(x, y, ax, ay, bx, by);
int b = orient(x, y, bx, by, cx, cy);
if (a == 0 || b == 0 || a < 0 == b < 0) {
int c = orient(x, y, cx, cy, ax, ay);
return c == 0 || (c < 0 == (b < 0 || a < 0));
}
return false;
} else {
return false;
}
}
/** Brute force compute if a point is in the polygon by traversing entire triangulation
* todo: speed this up using either binary tree or prefix coding (filtering by bounding box of triangle)
**/
public static final boolean pointInPolygon(final List tessellation, double lat, double lon) {
// each triangle
for (int i = 0; i < tessellation.size(); ++i) {
if (tessellation.get(i).containsPoint(lat, lon)) {
return true;
}
}
return false;
}
/** Circular Doubly-linked list used for polygon coordinates */
protected static class Node {
// node index in the linked list
private final int idx;
// vertex index in the polygon
private final int vrtxIdx;
// reference to the polygon for lat/lon values;
private final double[] polyX;
private final double[] polyY;
// encoded x value
private final int x;
// encoded y value
private final int y;
// morton code for sorting
private final long morton;
// previous node
private Node previous;
// next node
private Node next;
// previous z node
private Node previousZ;
// next z node
private Node nextZ;
// if the edge from this node to the next node is part of the polygon edges
private boolean isNextEdgeFromPolygon;
protected Node(final double[] x, final double[] y, final int index, final int vertexIndex, final boolean isGeo) {
this.idx = index;
this.vrtxIdx = vertexIndex;
this.polyX = x;
this.polyY = y;
// casting to float is safe as original values for non-geo are represented as floats
this.y = isGeo ? encodeLatitude(polyY[vrtxIdx]) : XYEncodingUtils.encode((float) polyY[vrtxIdx]);
this.x = isGeo ? encodeLongitude(polyX[vrtxIdx]) : XYEncodingUtils.encode((float) polyX[vrtxIdx]);
this.morton = BitUtil.interleave(this.x ^ 0x80000000, this.y ^ 0x80000000);
this.previous = null;
this.next = null;
this.previousZ = null;
this.nextZ = null;
this.isNextEdgeFromPolygon = true;
}
/** simple deep copy constructor */
protected Node(Node other) {
this.idx = other.idx;
this.vrtxIdx = other.vrtxIdx;
this.polyX = other.polyX;
this.polyY = other.polyY;
this.morton = other.morton;
this.x = other.x;
this.y = other.y;
this.previous = other.previous;
this.next = other.next;
this.previousZ = other.previousZ;
this.nextZ = other.nextZ;
this.isNextEdgeFromPolygon = other.isNextEdgeFromPolygon;
}
/** get the x value */
public final double getX() {
return polyX[vrtxIdx];
}
/** get the y value */
public final double getY() {
return polyY[vrtxIdx];
}
@Override
public String toString() {
StringBuilder builder = new StringBuilder();
if (this.previous == null)
builder.append("||-");
else
builder.append(this.previous.idx).append(" <- ");
builder.append(this.idx);
if (this.next == null)
builder.append(" -||");
else
builder.append(" -> ").append(this.next.idx);
return builder.toString();
}
}
/** Triangle in the tessellated mesh */
public final static class Triangle {
Node[] vertex;
boolean[] edgeFromPolygon;
protected Triangle(Node a, boolean isABfromPolygon, Node b, boolean isBCfromPolygon, Node c, boolean isCAfromPolygon) {
this.vertex = new Node[] {a, b, c};
this.edgeFromPolygon = new boolean[] {isABfromPolygon, isBCfromPolygon, isCAfromPolygon};
}
/** get quantized x value for the given vertex */
public int getEncodedX(int vertex) {
return this.vertex[vertex].x;
}
/** get quantized y value for the given vertex */
public int getEncodedY(int vertex) {
return this.vertex[vertex].y;
}
/** get y value for the given vertex */
public double getY(int vertex) {
return this.vertex[vertex].getY();
}
/** get x value for the given vertex */
public double getX(int vertex) {
return this.vertex[vertex].getX();
}
/** get if edge is shared with the polygon for the given edge */
public boolean isEdgefromPolygon(int startVertex) {
return edgeFromPolygon[startVertex];
}
/** utility method to compute whether the point is in the triangle */
protected boolean containsPoint(double lat, double lon) {
return pointInTriangle(lon, lat,
vertex[0].getX(), vertex[0].getY(),
vertex[1].getX(), vertex[1].getY(),
vertex[2].getX(), vertex[2].getY());
}
/** pretty print the triangle vertices */
public String toString() {
String result = vertex[0].x + ", " + vertex[0].y + " [" + edgeFromPolygon[0] + "] " +
vertex[1].x + ", " + vertex[1].y + " [" + edgeFromPolygon[1] + "] " +
vertex[2].x + ", " + vertex[2].y + " [" + edgeFromPolygon[2] + "]";
return result;
}
}
}