Apache
/
lucene
mirror of https://github.com/apache/lucene.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

LUCENE-8157: Do a better job of handling coplanar points in polygon construction. Thanks to Ignacio Vera for his help on this code.

This commit is contained in:
Karl Wright 2018-02-10 07:53:36 -05:00
parent 11a23a9029
commit 1fe34f1e89
4 changed files with 174 additions and 45 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

60
lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/GeoPolygonFactory.java
View File

@ -494,8 +494,7 @@ public class GeoPolygonFactory {
//all coplanar
return null;
}
Plane nextPlane = new Plane(startPoint, nextPoint);
if (Math.abs(nextPlane.evaluate(endPoint)) > Vector.MINIMUM_RESOLUTION + leniencyValue) {
if (!Plane.arePointsCoplanar(startPoint, endPoint, nextPoint)) {
//no coplanar.
break;
}
@ -521,7 +520,6 @@ public class GeoPolygonFactory {
return safePath;
}
/** Pick a random pole that has a good chance of being inside the polygon described by the points.
* @param generator is the random number generator to use.
* @param planetModel is the planet model to use.
@ -1079,12 +1077,19 @@ public class GeoPolygonFactory {
break;
}
final Edge newLastEdge = edgeBuffer.getNext(lastEdge);
if (Plane.arePointsCoplanar(lastEdge.startPoint, lastEdge.endPoint, newLastEdge.endPoint)) {
break;
}
if (isWithin(newLastEdge.endPoint, includedEdges)) {
//System.out.println(" maybe can extend to next edge");
// Found a candidate for extension. But do some other checks first. Basically, we need to know if we construct a polygon
// here will overlap with other remaining points?
final SidedPlane returnBoundary;
if (firstEdge.startPoint != newLastEdge.endPoint) {
if (Plane.arePointsCoplanar(firstEdge.endPoint, firstEdge.startPoint, newLastEdge.endPoint) ||
Plane.arePointsCoplanar(firstEdge.startPoint, newLastEdge.endPoint, newLastEdge.startPoint)) {
break;
}
returnBoundary = new SidedPlane(firstEdge.endPoint, firstEdge.startPoint, newLastEdge.endPoint);
} else {
returnBoundary = null;
@ -1135,12 +1140,19 @@ public class GeoPolygonFactory {
break;
}
final Edge newFirstEdge = edgeBuffer.getPrevious(firstEdge);
if (Plane.arePointsCoplanar(newFirstEdge.startPoint, newFirstEdge.endPoint, firstEdge.endPoint)) {
break;
}
if (isWithin(newFirstEdge.startPoint, includedEdges)) {
//System.out.println(" maybe can extend to previous edge");
// Found a candidate for extension. But do some other checks first. Basically, we need to know if we construct a polygon
// here will overlap with other remaining points?
final SidedPlane returnBoundary;
if (newFirstEdge.startPoint != lastEdge.endPoint) {
if(Plane.arePointsCoplanar(lastEdge.startPoint, lastEdge.endPoint, newFirstEdge.startPoint) ||
Plane.arePointsCoplanar(lastEdge.endPoint, newFirstEdge.startPoint, newFirstEdge.endPoint)) {
break;
}
returnBoundary = new SidedPlane(lastEdge.startPoint, lastEdge.endPoint, newFirstEdge.startPoint);
} else {
returnBoundary = null;
@ -1225,27 +1237,6 @@ public class GeoPolygonFactory {
edge = edgeBuffer.getNext(edge);
}
returnIsInternal = lastEdge.isInternal;
// Look for coplanarity; abort if so
for (int i = 0; i < points.size(); i++) {
final GeoPoint start = points.get(i);
final GeoPoint end = points.get(getLegalIndex(i + 1, points.size()));
// We have to find the next point that is not on the plane between start and end.
// If there is no such point, it's an error.
final Plane planeToFind = new Plane(start, end);
int endPointIndex = -1;
for (int j = 0; j < points.size(); j++) {
final int index = getLegalIndex(j + i + 2, points.size());
if (!planeToFind.evaluateIsZero(points.get(index))) {
endPointIndex = index;
break;
}
}
if (endPointIndex == -1) {
return false;
}
}
edgeBuffer.clear();
} else {
// Build the return edge (internal, of course)
@ -1270,27 +1261,6 @@ public class GeoPolygonFactory {
}
edge = edgeBuffer.getNext(edge);
}
// Look for coplanarity; abort if so
for (int i = 0; i < points.size(); i++) {
final GeoPoint start = points.get(i);
final GeoPoint end = points.get(getLegalIndex(i + 1, points.size()));
// We have to find the next point that is not on the plane between start and end.
// If there is no such point, it's an error.
final Plane planeToFind = new Plane(start, end);
int endPointIndex = -1;
for (int j = 0; j < points.size(); j++) {
final int index = getLegalIndex(j + i + 2, points.size());
if (!planeToFind.evaluateIsZero(points.get(index))) {
endPointIndex = index;
break;
}
}
if (endPointIndex == -1) {
return false;
}
}
// Modify the edge buffer
edgeBuffer.replace(edges, returnEdge);
}

15
lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/Plane.java
View File

@ -666,6 +666,21 @@ public class Plane extends Vector {
}
return findCrossings(planetModel, q, bounds, NO_BOUNDS);
}
/**
* Checks if three points are coplanar in any of the three planes they can describe.
* The planes are all assumed to go through the origin.
*
* @param A The first point.
* @param B The second point.
* @param C The third point
* @return true if provided points are coplanar in any of the three planes they can describe.
*/
public static boolean arePointsCoplanar(final GeoPoint A, final GeoPoint B, final GeoPoint C) {
return Vector.crossProductEvaluateIsZero(A, B, C) ||
Vector.crossProductEvaluateIsZero(A, C, B) ||
Vector.crossProductEvaluateIsZero(B, C, A);
}
/**
* Find the intersection points between two planes, given a set of bounds.

73
lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/Vector.java
View File

@ -180,6 +180,79 @@ public class Vector {
return new Vector(x * normFactor, y * normFactor, z * normFactor);
}
/**
* Evaluate the cross product of two vectors against a point.
* If the dot product of the resultant vector resolves to "zero", then
* return true.
* @param A is the first vector to use for the cross product.
* @param B is the second vector to use for the cross product.
* @param point is the point to evaluate.
* @return true if we get a zero dot product.
*/
public static boolean crossProductEvaluateIsZero(final Vector A, final Vector B, final Vector point) {
// Include Gram-Schmidt in-line so we avoid creating objects unnecessarily
// Compute the naive perpendicular
final double thisX = A.y * B.z - A.z * B.y;
final double thisY = A.z * B.x - A.x * B.z;
final double thisZ = A.x * B.y - A.y * B.x;
final double magnitude = magnitude(thisX, thisY, thisZ);
if (magnitude < MINIMUM_RESOLUTION) {
return true;
}
final double inverseMagnitude = 1.0/magnitude;
double normalizeX = thisX * inverseMagnitude;
double normalizeY = thisY * inverseMagnitude;
double normalizeZ = thisZ * inverseMagnitude;
// For a plane to work, the dot product between the normal vector
// and the points needs to be less than the minimum resolution.
// This is sometimes not true for points that are very close. Therefore
// we need to adjust
int i = 0;
while (true) {
final double currentDotProdA = A.x * normalizeX + A.y * normalizeY + A.z * normalizeZ;
final double currentDotProdB = B.x * normalizeX + B.y * normalizeY + B.z * normalizeZ;
if (Math.abs(currentDotProdA) < MINIMUM_RESOLUTION && Math.abs(currentDotProdB) < MINIMUM_RESOLUTION) {
break;
}
// Converge on the one that has largest dot product
final double currentVectorX;
final double currentVectorY;
final double currentVectorZ;
final double currentDotProd;
if (Math.abs(currentDotProdA) > Math.abs(currentDotProdB)) {
currentVectorX = A.x;
currentVectorY = A.y;
currentVectorZ = A.z;
currentDotProd = currentDotProdA;
} else {
currentVectorX = B.x;
currentVectorY = B.y;
currentVectorZ = B.z;
currentDotProd = currentDotProdB;
}
// Adjust
normalizeX = normalizeX - currentDotProd * currentVectorX;
normalizeY = normalizeY - currentDotProd * currentVectorY;
normalizeZ = normalizeZ - currentDotProd * currentVectorZ;
// Normalize
final double correctedMagnitude = magnitude(normalizeX, normalizeY, normalizeZ);
final double inverseCorrectedMagnitude = 1.0 / correctedMagnitude;
normalizeX = normalizeX * inverseCorrectedMagnitude;
normalizeY = normalizeY * inverseCorrectedMagnitude;
normalizeZ = normalizeZ * inverseCorrectedMagnitude;
//This is probably not needed as the method seems to converge
//quite quickly. But it is safer to have a way out.
if (i++ > 10) {
throw new IllegalArgumentException("Plane could not be constructed! Could not find a normal vector.");
}
}
return Math.abs(normalizeX * point.x + normalizeY * point.y + normalizeZ * point.z) < MINIMUM_RESOLUTION;
}
/**
* Do a dot product.
*

71
lucene/spatial3d/src/test/org/apache/lucene/spatial3d/geom/RandomGeoPolygonTest.java Normal file
View File

@ -0,0 +1,71 @@
package org.apache.lucene.spatial3d.geom;
import java.util.ArrayList;
import java.util.List;
import com.carrotsearch.randomizedtesting.annotations.Repeat;
import org.junit.Test;
/**
* Random test for polygons.
*/
public class RandomGeoPolygonTest extends RandomGeo3dShapeGenerator {
@Test
@Repeat(iterations = 10)
public void testRandomLUCENE8157() {
final PlanetModel planetModel = randomPlanetModel();
final GeoPoint startPoint = randomGeoPoint(planetModel);
double d = random().nextDouble();
final double distanceSmall = d * 1e-9 + Vector.MINIMUM_ANGULAR_RESOLUTION;
final double distanceBig = d * 1e-7 + Vector.MINIMUM_ANGULAR_RESOLUTION ;
final double bearing = random().nextDouble() * Math.PI;
GeoPoint point1 = planetModel.surfacePointOnBearing(startPoint, distanceSmall, bearing*1.001);
GeoPoint point2 = planetModel.surfacePointOnBearing(startPoint, distanceBig, bearing);
GeoPoint point3 = planetModel.surfacePointOnBearing(startPoint, distanceBig, bearing - 0.5 * Math.PI);
List<GeoPoint> points = new ArrayList<>();
points.add(startPoint);
points.add(point1);
points.add(point2);
points.add(point3);
try {
GeoPolygon polygon = GeoPolygonFactory.makeGeoPolygon(planetModel, points);
assertTrue(polygon != null);
}
catch(Exception e) {
fail(points.toString());
}
}
public void testLUCENE8157() {
GeoPoint point1 = new GeoPoint(PlanetModel.SPHERE, 0.281855362988772, -0.7816673189809037);
GeoPoint point2 = new GeoPoint(PlanetModel.SPHERE, 0.28185536309057774, -0.7816673188511931);
GeoPoint point3 = new GeoPoint(PlanetModel.SPHERE, 0.28186535556824205, -0.7816546103463846);
GeoPoint point4 = new GeoPoint(PlanetModel.SPHERE, 0.28186757010406716, -0.7816777221140381);
List<GeoPoint> points = new ArrayList<>();
points.add(point1);
points.add(point2);
points.add(point3);
points.add(point4);
try {
GeoPolygon polygon = GeoPolygonFactory.makeGeoPolygon(PlanetModel.SPHERE, points);
}
catch(Exception e) {
fail(points.toString());
}
}
@Test
public void testCoplanarityTilePolygon() {
//POLYGON((-90.55764 -0.34907,-90.55751 -0.34868,-90.55777 -0.34842,-90.55815 -0.34766,-90.55943 -0.34842, -90.55918 -0.34842,-90.55764 -0.34907))
List<GeoPoint> points = new ArrayList<>();
points.add(new GeoPoint(PlanetModel.SPHERE, Geo3DUtil.fromDegrees(-0.34907), Geo3DUtil.fromDegrees(-90.55764)));
points.add(new GeoPoint(PlanetModel.SPHERE, Geo3DUtil.fromDegrees(-0.34868), Geo3DUtil.fromDegrees(-90.55751)));
points.add(new GeoPoint(PlanetModel.SPHERE, Geo3DUtil.fromDegrees(-0.34842), Geo3DUtil.fromDegrees(-90.55777)));
points.add(new GeoPoint(PlanetModel.SPHERE, Geo3DUtil.fromDegrees(-0.34766), Geo3DUtil.fromDegrees(-90.55815)));
points.add(new GeoPoint(PlanetModel.SPHERE, Geo3DUtil.fromDegrees(-0.34842), Geo3DUtil.fromDegrees(-90.55943)));
points.add(new GeoPoint(PlanetModel.SPHERE, Geo3DUtil.fromDegrees(-0.34842), Geo3DUtil.fromDegrees(-90.55918)));
GeoCompositePolygon polygon = (GeoCompositePolygon)GeoPolygonFactory.makeGeoPolygon(PlanetModel.SPHERE, points);
assertTrue(polygon.size() == 3);
}
}