Merge branch 'GITHUB-11883' into main

Pulling in changes to address ticket 11883.
This commit is contained in:
Karl David Wright 2022-11-25 16:32:02 -05:00
commit 5c4896321d
4 changed files with 177 additions and 43 deletions

View File

@ -42,14 +42,16 @@ class GeoConvexPolygon extends GeoBasePolygon {
/** A list of edges */
protected SidedPlane[] edges = null;
/** A list of edge starting bounding planes */
protected SidedPlane[] startBounds = null;
/** A list of edge ending bounding planes */
protected SidedPlane[] endBounds = null;
/** The set of notable points for each edge */
protected GeoPoint[][] notableEdgePoints = null;
/** A point which is on the boundary of the polygon */
protected GeoPoint[] edgePoints = null;
/** Set to true when the polygon is complete */
protected boolean isDone = false;
/** A bounds object for each sided plane */
protected Map<SidedPlane, Membership> eitherBounds = null;
/** Map from edge to its previous non-coplanar brother */
protected Map<SidedPlane, SidedPlane> prevBrotherMap = null;
/** Map from edge to its next non-coplanar brother */
@ -213,6 +215,8 @@ class GeoConvexPolygon extends GeoBasePolygon {
// Time to construct the planes. If the polygon is truly convex, then any adjacent point
// to a segment can provide an interior measurement.
edges = new SidedPlane[points.size()];
startBounds = new SidedPlane[points.size()];
endBounds = new SidedPlane[points.size()];
notableEdgePoints = new GeoPoint[points.size()][];
for (int i = 0; i < points.size(); i++) {
@ -235,11 +239,12 @@ class GeoConvexPolygon extends GeoBasePolygon {
final GeoPoint check = points.get(endPointIndex);
final SidedPlane sp = new SidedPlane(check, start, end);
edges[i] = sp;
startBounds[i] = SidedPlane.constructSidedPlaneFromOnePoint(end, sp, start);
endBounds[i] = SidedPlane.constructSidedPlaneFromOnePoint(start, sp, end);
notableEdgePoints[i] = new GeoPoint[] {start, end};
}
// For each edge, create a bounds object.
eitherBounds = CollectionUtil.newHashMap(edges.length);
prevBrotherMap = CollectionUtil.newHashMap(edges.length);
nextBrotherMap = CollectionUtil.newHashMap(edges.length);
for (int edgeIndex = 0; edgeIndex < edges.length; edgeIndex++) {
@ -273,7 +278,6 @@ class GeoConvexPolygon extends GeoBasePolygon {
"Convex polygon has a side that is more than 180 degrees");
}
}
eitherBounds.put(edge, new EitherBound(edges[bound1Index], edges[bound2Index]));
// When we are done with this cycle, we'll need to build the intersection bound for each edge
// and its brother. For now, keep track of the relationships.
nextBrotherMap.put(edge, edges[bound1Index]);
@ -411,7 +415,13 @@ class GeoConvexPolygon extends GeoBasePolygon {
// System.err.println("Checking convex edge " + edge
// + " for intersection against plane " + p);
if (edge.intersects(
planetModel, p, notablePoints, points, bounds, eitherBounds.get(edge))) {
planetModel,
p,
notablePoints,
points,
bounds,
startBounds[edgeIndex],
endBounds[edgeIndex])) {
// System.err.println(" intersects!");
return true;
}
@ -436,7 +446,7 @@ class GeoConvexPolygon extends GeoBasePolygon {
final SidedPlane edge = edges[edgeIndex];
final GeoPoint[] points = this.notableEdgePoints[edgeIndex];
if (!isInternalEdges.get(edgeIndex)) {
if (shape.intersects(edge, points, eitherBounds.get(edge))) {
if (shape.intersects(edge, points, startBounds[edgeIndex], endBounds[edgeIndex])) {
return true;
}
}
@ -451,39 +461,6 @@ class GeoConvexPolygon extends GeoBasePolygon {
return false;
}
/** A membership implementation representing polygon edges that must apply. */
protected static class EitherBound implements Membership {
protected final SidedPlane sideBound1;
protected final SidedPlane sideBound2;
/**
* Constructor.
*
* @param sideBound1 is the first side bound.
* @param sideBound2 is the second side bound.
*/
public EitherBound(final SidedPlane sideBound1, final SidedPlane sideBound2) {
this.sideBound1 = sideBound1;
this.sideBound2 = sideBound2;
}
@Override
public boolean isWithin(final Vector v) {
return sideBound1.isWithin(v) && sideBound2.isWithin(v);
}
@Override
public boolean isWithin(final double x, final double y, final double z) {
return sideBound1.isWithin(x, y, z) && sideBound2.isWithin(x, y, z);
}
@Override
public String toString() {
return "(" + sideBound1 + "," + sideBound2 + ")";
}
}
@Override
public void getBounds(Bounds bounds) {
// Because of holes, we don't want to use superclass method
@ -512,8 +489,9 @@ class GeoConvexPolygon extends GeoBasePolygon {
}
// Add planes with membership.
for (final SidedPlane edge : edges) {
bounds.addPlane(planetModel, edge, eitherBounds.get(edge));
for (int edgeIndex = 0; edgeIndex < edges.length; edgeIndex++) {
final SidedPlane edge = edges[edgeIndex];
bounds.addPlane(planetModel, edge, startBounds[edgeIndex], endBounds[edgeIndex]);
final SidedPlane nextEdge = nextBrotherMap.get(edge);
bounds.addIntersection(
planetModel, edge, nextEdge, prevBrotherMap.get(edge), nextBrotherMap.get(nextEdge));
@ -530,10 +508,11 @@ class GeoConvexPolygon extends GeoBasePolygon {
minimumDistance = newDist;
}
}
for (final SidedPlane edgePlane : edges) {
for (int edgeIndex = 0; edgeIndex < edges.length; edgeIndex++) {
final SidedPlane edgePlane = edges[edgeIndex];
final double newDist =
distanceStyle.computeDistance(
planetModel, edgePlane, x, y, z, eitherBounds.get(edgePlane));
planetModel, edgePlane, x, y, z, startBounds[edgeIndex], endBounds[edgeIndex]);
if (newDist < minimumDistance) {
minimumDistance = newDist;
}

View File

@ -126,6 +126,123 @@ public class Plane extends Vector {
: Math.nextDown(basePlane.D - MINIMUM_RESOLUTION));
}
/**
* Given a plane and one point that is on that plane, find a perpendicular plane that goes through
* both the origin and the point.
*
* @param plane is the original plane
* @param M is the point on that plane
* @return a plane that goes through M, the origin, and is perpendicular to the original plane
*/
public static Plane constructPerpendicularCenterPlaneOnePoint(final Plane plane, final Vector M) {
// Original plane:
// A0x + B0y + C0z = 0 (D0 = 0)
assert plane.evaluateIsZero(M);
final double A0 = plane.x;
final double B0 = plane.y;
final double C0 = plane.z;
// Second plane equations:
// A1Mx + B1My + C1Mz + D1 = 0
// A0*A1 + B0*B1 + C0*C1 = 0
// A1^2 + B1^2 + C1^2 = 1
// D1 = 0.0 because it goes through origin.
// Basic strategy: Pick a variable and set it to 1.
final double a1Denom = C0 * M.y - B0 * M.z;
final double b1Denom = C0 * M.x - A0 * M.z;
final double c1Denom = B0 * M.x - A0 * M.y;
final double A1;
final double B1;
final double C1;
if (Math.abs(a1Denom) >= Math.abs(b1Denom) && Math.abs(a1Denom) >= Math.abs(c1Denom)) {
A1 = 1.0;
if (Math.abs(M.y) >= Math.abs(M.z)) {
// e.g. A1 = 1.
// Then:
//
// Mx + B1 My + C1 Mz = 0
// B1 = (-Mx - C1 Mz) / My
// A0 + B0 * (-Mx - C1Mz)/My + C0*C1 = 0
// A0 * My - B0 * Mx - B0 * C1 * Mz + C0 * C1 * My = 0
// C1 (C0 * My - B0 * Mz) = B0 * Mx - A0 * My
// C1 = (B0 * Mx - A0 * My) / (C0 * My - B0 * Mz)
C1 = (B0 * M.x - A0 * M.y) / a1Denom;
B1 = (-M.x - C1 * M.z) / M.y;
} else {
// Alternative substitution sequence:
// C1 = (-Mx - B1 My) / Mz
// A0 + B0 * B1 + C0 * (-Mx - B1 My) / Mz = 0
// A0 * Mz + B0 * B1 * Mz - C0 * Mx - C0 * B1 * My = 0
// B1 (B0 * Mz - C0 * My) = C0 * Mx - A0 * Mz
// B1 = (C0 * Mx - A0 * Mz) / (B0 * Mz - C0 * My)
B1 = (A0 * M.z - C0 * M.x) / a1Denom;
C1 = (-M.x - B1 * M.y) / M.z;
}
} else if (Math.abs(b1Denom) >= Math.abs(a1Denom) && Math.abs(b1Denom) >= Math.abs(c1Denom)) {
B1 = 1.0;
if (Math.abs(M.x) >= Math.abs(M.z)) {
// B1 = 1
// Then:
//
// A1 * Mx + My + C1 * Mz = 0
// A1 = (-My - C1 * Mz) / Mx
// A0 * (-My - C1 * Mz) / Mx + B0 + C0 * C1 = 0
// -A0 * My - C1 * Mz + B0 * Mx + C0 * C1 * Mx = 0
// C1 (C0 * Mx - A0 * Mz) = A0 * My - B0 * Mx
// C1 = (A0 * My - B0 * Mx) / (C0 * Mx - A0 * Mz)
C1 = (A0 * M.y - B0 * M.x) / b1Denom;
A1 = (-M.y - C1 * M.z) / M.x;
} else {
// Alternative:
// C1 = (-My - A1 * Mx) / Mz
// A0 * A1 + B0 + C0 * (-My - A1 * Mx) / Mz = 0
// A0 * A1 * Mz + B0 * Mz - C0 * My - C0 * A1 * Mx = 0
// A1 (A0 * Mz - C0 * Mx) = C0 * My - B0 * Mz
// A1 = (C0 * My - B0 * Mz) / (A0 * Mz - C0 * Mx)
A1 = (B0 * M.z - C0 * M.y) / b1Denom;
C1 = (-M.y - A1 * M.x) / M.z;
}
} else if (Math.abs(c1Denom) >= Math.abs(a1Denom) && Math.abs(c1Denom) >= Math.abs(b1Denom)) {
C1 = 1.0;
if (Math.abs(M.x) >= Math.abs(M.y)) {
// C1 = 1
// Then:
//
// A1 * Mx + B1 * My + Mz = 0
// A1 = (-Mz - B1 * My)/Mx
// A0 * (-Mz - B1 * My)/Mx + B0 * B1 + C0 = 0
// - A0 * Mz - A0 * B1 * My + B0 * B1 * Mx + C0 * Mx = 0
// B1 (B0 * Mx - A0 * My) = A0 * Mz - C0 * Mx
// B1 = (A0 * Mz - C0 * Mx) / (B0 * Mx - A0 * My)
B1 = (A0 * M.z - C0 * M.x) / c1Denom;
A1 = (-M.z - B1 * M.y) / M.x;
} else {
// Alternative:
// B1 = (-Mz - A1 * Mx) / My
// A0 * A1 + B0 * (-Mz - A1 * Mx) / My + C0 = 0
// A0 * A1 * My - B0 * Mz - B0 * A1 * Mx + C0 * My = 0
// A1 (A0 * My - B0 * Mx) = B0 * Mz - C0 * My
// A1 = (B0 * Mz - C0 * My) / (A0 * My - B0 * Mx)
A1 = (C0 * M.y - B0 * M.z) / c1Denom;
B1 = (-M.z - A1 * M.x) / M.y;
}
} else {
throw new IllegalArgumentException("Cannot find perpendicular plane as requested");
}
// Normalize the vector
final double normFactor = 1.0 / Math.sqrt(A1 * A1 + B1 * B1 + C1 * C1);
final Vector v = new Vector(A1 * normFactor, B1 * normFactor, C1 * normFactor);
final Plane rval = new Plane(v, -(v.x * M.x + v.y * M.y + v.z * M.z));
return rval;
}
/**
* Given two points, construct a plane that goes through them and the origin. Then, find a plane
* that is perpendicular to that which also goes through the original two points. This is useful

View File

@ -237,6 +237,17 @@ public class SidedPlane extends Plane implements Membership {
return new SidedPlane(insidePoint, plane.x, plane.y, plane.z, plane.D);
}
/**
* Construct sided plane from a plane and one point. This finds a plane perpendicular to the
* passed-in plane, and goes through both the origin and the point.
*/
public static SidedPlane constructSidedPlaneFromOnePoint(
final Vector insidePoint, final Plane plane, final Vector intersectionPoint) {
final Plane newPlane =
Plane.constructPerpendicularCenterPlaneOnePoint(plane, intersectionPoint);
return new SidedPlane(insidePoint, newPlane.x, newPlane.y, newPlane.z, newPlane.D);
}
/** Construct a sided plane from three points. */
public static SidedPlane constructNormalizedThreePointSidedPlane(
final Vector insidePoint, final Vector point1, final Vector point2, final Vector point3) {

View File

@ -25,6 +25,33 @@ import org.junit.Test;
public class TestGeoPolygon extends LuceneTestCase {
@Test
public void testH3CellsWrongIntersection() {
final List<GeoPoint> points1 = new ArrayList<>();
addToList(points1, PlanetModel.SPHERE, -64.2102198418716, -39.14233318389477);
addToList(points1, PlanetModel.SPHERE, -64.21016450005413, -39.142267144439614);
addToList(points1, PlanetModel.SPHERE, -64.21021077465937, -39.1421844504783);
addToList(points1, PlanetModel.SPHERE, -64.21031239098929, -39.142167795785085);
addToList(points1, PlanetModel.SPHERE, -64.2103677330169, -39.14223383513698);
addToList(points1, PlanetModel.SPHERE, -64.21032145850448, -39.14231652928534);
final List<GeoPoint> points2 = new ArrayList<>();
addToList(points2, PlanetModel.SPHERE, -64.20991499254879, -39.14238314705201);
addToList(points2, PlanetModel.SPHERE, -64.20985965132967, -39.14231710755475);
addToList(points2, PlanetModel.SPHERE, -64.20990592624541, -39.142234413886875);
addToList(points2, PlanetModel.SPHERE, -64.21000754228744, -39.142217759529174);
addToList(points2, PlanetModel.SPHERE, -64.21006288371667, -39.14228379892317);
addToList(points2, PlanetModel.SPHERE, -64.21001660889377, -39.14236649277811);
final GeoPolygon polygon1 = GeoPolygonFactory.makeGeoPolygon(PlanetModel.SPHERE, points1);
final GeoPolygon polygon2 = GeoPolygonFactory.makeGeoPolygon(PlanetModel.SPHERE, points2);
assertFalse(polygon1.intersects(polygon2));
}
private static void addToList(
List<GeoPoint> points, PlanetModel planetModel, double lon, double lat) {
points.add(new GeoPoint(planetModel, Geo3DUtil.fromDegrees(lat), Geo3DUtil.fromDegrees(lon)));
}
@Test
public void testPolygonPointFiltering() {
final GeoPoint point1 = new GeoPoint(PlanetModel.WGS84, 1.0, 2.0);