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

LUCENE-8136: Adopt iterative convergence for construction of planes from two vectors. Thanks, Ignacio Vera!

This commit is contained in:
Karl Wright 2018-01-25 07:03:45 -05:00
parent 653935bbdf
commit ad131bde32
4 changed files with 160 additions and 31 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/GeoPoint.java
View File

@ -189,6 +189,14 @@ public class GeoPoint extends Vector implements SerializableObject {
return mag;
}
/** Compute whether point matches another.
*@param p is the other point.
*@return true if the same.
*/
public boolean isIdentical(final GeoPoint p) {
return isIdentical(p.x, p.y, p.z);
}
/** Compute whether point matches another.
*@param x is the x value
*@param y is the y value

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

@ -27,7 +27,7 @@ public class Vector {
* Values that are all considered to be essentially zero have a magnitude
* less than this.
*/
public static final double MINIMUM_RESOLUTION = 1.5e-12;
public static final double MINIMUM_RESOLUTION = 1.0e-12;
/**
* Angular version of minimum resolution.
*/
@ -72,20 +72,78 @@ public class Vector {
* @param BZ is the Z value of the second
*/
public Vector(final Vector A, final double BX, final double BY, final double BZ) {
// We're really looking at two vectors and computing a perpendicular one from that.
// We can think of this as having three points -- the origin, and two points that aren't the origin.
// Normally, we can compute the perpendicular vector this way:
// x = u2v3 - u3v2
// y = u3v1 - u1v3
// z = u1v2 - u2v1
// Sometimes that produces a plane that does not contain the original three points, however, due to
// numerical precision issues. Then we continue making the answer more precise using the
// Gram-Schmidt process: https://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process
// Compute the naive perpendicular
final double thisX = A.y * BZ - A.z * BY;
final double thisY = A.z * BX - A.x * BZ;
final double thisZ = A.x * BY - A.y * BX;
final double magnitude = magnitude(thisX, thisY, thisZ);
if (Math.abs(magnitude) < MINIMUM_RESOLUTION) {
if (magnitude < MINIMUM_RESOLUTION) {
throw new IllegalArgumentException("Degenerate/parallel vector constructed");
}
final double inverseMagnitude = 1.0 / magnitude;
this.x = thisX * inverseMagnitude;
this.y = thisY * inverseMagnitude;
this.z = thisZ * inverseMagnitude;
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 = BX * normalizeX + BY * normalizeY + BZ * 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 = BX;
currentVectorY = BY;
currentVectorZ = BZ;
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.");
}
}
this.x = normalizeX;
this.y = normalizeY;
this.z = normalizeZ;
}
/**
@ -98,20 +156,7 @@ public class Vector {
* @param B is the second
*/
public Vector(final Vector A, final Vector B) {
// x = u2v3 - u3v2
// y = u3v1 - u1v3
// z = u1v2 - u2v1
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 (Math.abs(magnitude) < MINIMUM_RESOLUTION) {
throw new IllegalArgumentException("Degenerate/parallel vector constructed");
}
final double inverseMagnitude = 1.0 / magnitude;
this.x = thisX * inverseMagnitude;
this.y = thisY * inverseMagnitude;
this.z = thisZ * inverseMagnitude;
this(A, B.x, B.y, B.z);
}
/** Compute a magnitude of an x,y,z value.

24
lucene/spatial3d/src/test/org/apache/lucene/spatial3d/geom/GeoPolygonTest.java
View File

@ -1035,17 +1035,19 @@ shape:
GeoPoint point1 = new GeoPoint(PlanetModel.SPHERE, Geo3DUtil.fromDegrees(-23.434456), Geo3DUtil.fromDegrees(14.459204));
GeoPoint point2 = new GeoPoint(PlanetModel.SPHERE, Geo3DUtil.fromDegrees(-23.43394), Geo3DUtil.fromDegrees(14.459206));
GeoPoint check = new GeoPoint(PlanetModel.SPHERE, Geo3DUtil.fromDegrees(-23.434067), Geo3DUtil.fromDegrees(14.458927));
SidedPlane plane = new SidedPlane(check, point1, point2);
assertTrue(plane.isWithin(check));
assertTrue(plane.isWithin(point1));
assertTrue(plane.isWithin(point2));
//POLYGON((14.459204 -23.434456, 14.459206 -23.43394,14.458647 -23.434196, 14.458646 -23.434452,14.459204 -23.434456))
List<GeoPoint> points = new ArrayList<>();
points.add(new GeoPoint(PlanetModel.SPHERE, Geo3DUtil.fromDegrees(-23.434456), Geo3DUtil.fromDegrees(14.459204)));
points.add(new GeoPoint(PlanetModel.SPHERE, Geo3DUtil.fromDegrees( -23.43394), Geo3DUtil.fromDegrees(14.459206)));
points.add(new GeoPoint(PlanetModel.SPHERE, Geo3DUtil.fromDegrees(-23.434196), Geo3DUtil.fromDegrees(14.458647)));
points.add(new GeoPoint(PlanetModel.SPHERE, Geo3DUtil.fromDegrees(-23.434452), Geo3DUtil.fromDegrees(14.458646)));
GeoPolygonFactory.makeGeoPolygon(PlanetModel.SPHERE, points);
if (!point1.isIdentical(point2) && !check.isIdentical(point1) && !check.isIdentical(point2)) {
SidedPlane plane = new SidedPlane(check, point1, point2);
assertTrue(plane.isWithin(check));
assertTrue(plane.isWithin(point1));
assertTrue(plane.isWithin(point2));
//POLYGON((14.459204 -23.434456, 14.459206 -23.43394,14.458647 -23.434196, 14.458646 -23.434452,14.459204 -23.434456))
List<GeoPoint> points = new ArrayList<>();
points.add(new GeoPoint(PlanetModel.SPHERE, Geo3DUtil.fromDegrees(-23.434456), Geo3DUtil.fromDegrees(14.459204)));
points.add(new GeoPoint(PlanetModel.SPHERE, Geo3DUtil.fromDegrees( -23.43394), Geo3DUtil.fromDegrees(14.459206)));
points.add(new GeoPoint(PlanetModel.SPHERE, Geo3DUtil.fromDegrees(-23.434196), Geo3DUtil.fromDegrees(14.458647)));
points.add(new GeoPoint(PlanetModel.SPHERE, Geo3DUtil.fromDegrees(-23.434452), Geo3DUtil.fromDegrees(14.458646)));
GeoPolygonFactory.makeGeoPolygon(PlanetModel.SPHERE, points);
}
}

74
lucene/spatial3d/src/test/org/apache/lucene/spatial3d/geom/RandomPlaneTest.java Normal file
View File

@ -0,0 +1,74 @@
/*
* 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.spatial3d.geom;
import java.util.ArrayList;
import java.util.List;
import com.carrotsearch.randomizedtesting.annotations.Repeat;
import org.junit.Test;
/**
* Random test for planes.
*/
public class RandomPlaneTest extends RandomGeo3dShapeGenerator {
@Test
@Repeat(iterations = 10)
public void testPlaneAccuracy() {
PlanetModel planetModel = randomPlanetModel();
GeoPoint point1 = randomGeoPoint(planetModel);
for (int i= 0; i < 1000; i++) {
double dist = random().nextDouble() * Vector.MINIMUM_ANGULAR_RESOLUTION + Vector.MINIMUM_ANGULAR_RESOLUTION;
double bearing = random().nextDouble() * 2 * Math.PI;
GeoPoint point2 = planetModel.surfacePointOnBearing(point1, dist, bearing );
GeoPoint check = randomGeoPoint(planetModel);
if (!point1.isNumericallyIdentical(point2)) {
SidedPlane plane = new SidedPlane(check, point1, point2);
String msg = dist + " point 1: " + point1 + ", point 2: " + point2 + " , check: " + check;
assertTrue(msg, plane.isWithin(check));
assertTrue(msg, plane.isWithin(point2));
assertTrue(msg, plane.isWithin(point1));
}
else {
assertFalse("numerically identical", true);
}
}
}
/*
@Test
@Repeat(iterations = 10)
public void testPolygonAccuracy() {
PlanetModel planetModel = randomPlanetModel();
GeoPoint point1 = randomGeoPoint(planetModel);
for (int i= 0; i < 1000; i++) {
double dist = random().nextDouble() * 1e-6 + Vector.MINIMUM_ANGULAR_RESOLUTION;
GeoPoint point2 = planetModel.surfacePointOnBearing(point1, dist, 0);
GeoPoint point3 = planetModel.surfacePointOnBearing(point1, dist, 0.5 * Math.PI);
List<GeoPoint> points = new ArrayList<>();
points.add(point1);
points.add(point2);
points.add(point3);
GeoPolygonFactory.makeGeoPolygon(planetModel, points);
}
}
*/
}