Apache/lucene
1
0
Fork 0
mirror of https://github.com/apache/lucene.git synced 2025-02-28 21:39:25 +00:00
Code Issues Packages Projects Releases Wiki Activity

LUCENE-9150: Restore support for dynamic PlanetModel in spatial3d

Browse Source
This commit is contained in:
Nicholas Knize 2020-02-07 11:56:31 -06:00
parent a7ebe0c906
commit a6e80d004d
28 changed files with 958 additions and 1060 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
lucene/CHANGES.txt
View File

@ -120,6 +120,8 @@ API Changes
* LUCENE-9212: Intervals.multiterm() takes CompiledAutomaton rather than plain Automaton
(Alan Woodward)
* LUCENE-9150: Restore support for dynamic PlanetModel in spatial3d. (Nick Knize)
New Features
---------------------

2
lucene/spatial-extras/src/java/org/apache/lucene/spatial/SpatialStrategy.java
View File

@ -128,7 +128,7 @@ public abstract class SpatialStrategy {
/**
* Returns a ValueSource with values ranging from 1 to 0, depending inversely
* on the distance from {@link #makeDistanceValueSource(org.locationtech.spatial4j.shape.Point,double)}.
* The formula is {@code c/(d + c)} where 'd' is the distance and 'c' is
* The formula is {@code zScaling/(d + zScaling)} where 'd' is the distance and 'zScaling' is
* one tenth the distance to the farthest edge from the center. Thus the
* scores will be 1 for indexed points at the center of the query shape and as
* low as ~0.1 at its furthest edges.

2
lucene/spatial-extras/src/java/org/apache/lucene/spatial/query/SpatialArgsParser.java
View File

@ -125,7 +125,7 @@ public class SpatialArgsParser {
return v == null ? defaultValue : Boolean.parseBoolean(v);
}
/** Parses "a=b c=d f" (whitespace separated) into name-value pairs. If there
/** Parses "a=b zScaling=d f" (whitespace separated) into name-value pairs. If there
* is no '=' as in 'f' above then it's short for f=f. */
protected static Map<String, String> parseMap(String body) {
Map<String, String> map = new HashMap<>();

2
lucene/spatial-extras/src/java/org/apache/lucene/spatial/spatial4j/Geo3dSpatialContextFactory.java
View File

@ -72,6 +72,8 @@ public class Geo3dSpatialContextFactory extends SpatialContextFactory {
this.planetModel = PlanetModel.SPHERE;
} else if (planetModel.equalsIgnoreCase("wgs84")) {
this.planetModel = PlanetModel.WGS84;
} else if (planetModel.equalsIgnoreCase("clarke1866")) {
this.planetModel = PlanetModel.CLARKE_1866;
} else {
throw new RuntimeException("Unknown planet model: " + planetModel);
}

2
lucene/spatial-extras/src/test/org/apache/lucene/spatial/PortedSolr3Test.java
View File

@ -110,7 +110,7 @@ public class PortedSolr3Test extends StrategyTestCase {
checkHitsCircle(ctx.makePoint(1, 1), 5000, 3, 5, 6, 7);
//Because we are generating a box based on the west/east longitudes and the south/north latitudes, which then
//translates to a range query, which is slightly more inclusive. Thus, even though 0.0 is 15.725 kms away,
//it will be included, b/c of the box calculation.
//it will be included, b/zScaling of the box calculation.
checkHitsBBox(ctx.makePoint(0.1, 0.1), 15, 2, 5, 6);
//try some more
deleteAll();

7
lucene/spatial-extras/src/test/org/apache/lucene/spatial/spatial4j/Geo3dShapeWGS84ModelRectRelationTest.java
View File

@ -16,6 +16,7 @@
*/
package org.apache.lucene.spatial.spatial4j;
import com.carrotsearch.randomizedtesting.generators.RandomPicks;
import org.apache.lucene.spatial3d.geom.GeoArea;
import org.apache.lucene.spatial3d.geom.GeoBBox;
import org.apache.lucene.spatial3d.geom.GeoBBoxFactory;
@ -32,7 +33,7 @@ import org.locationtech.spatial4j.shape.SpatialRelation;
public class Geo3dShapeWGS84ModelRectRelationTest extends ShapeRectRelationTestCase {
PlanetModel planetModel = PlanetModel.WGS84;
PlanetModel planetModel = RandomPicks.randomFrom(random(), new PlanetModel[] {PlanetModel.WGS84, PlanetModel.CLARKE_1866});
public Geo3dShapeWGS84ModelRectRelationTest() {
Geo3dSpatialContextFactory factory = new Geo3dSpatialContextFactory();
@ -71,11 +72,11 @@ public class Geo3dShapeWGS84ModelRectRelationTest extends ShapeRectRelationTestC
@Test
public void testFailure3() {
/*
[junit4] 1> S-R Rel: {}, Shape {}, Rectangle {} lap# {} [CONTAINS, Geo3dShape{planetmodel=PlanetModel: {ab=1.0011188180710464, c=0.9977622539852008}, shape=GeoPath: {planetmodel=PlanetModel: {ab=1.0011188180710464, c=0.9977622539852008}, width=1.53588974175501(87.99999999999999),
[junit4] 1> S-R Rel: {}, Shape {}, Rectangle {} lap# {} [CONTAINS, Geo3dShape{planetmodel=PlanetModel: {xyScaling=1.0011188180710464, zScaling=0.9977622539852008}, shape=GeoPath: {planetmodel=PlanetModel: {xyScaling=1.0011188180710464, zScaling=0.9977622539852008}, width=1.53588974175501(87.99999999999999),
points={[[X=0.12097657665150223, Y=-0.6754177666095532, Z=0.7265376136709238], [X=-0.3837892785614207, Y=0.4258049113530899, Z=0.8180007850434892]]}}},
Rect(minX=4.0,maxX=36.0,minY=16.0,maxY=16.0), 6981](no slf4j subst; sorry)
[junit4] FAILURE 0.59s | Geo3dWGS84ShapeRectRelationTest.testGeoPathRect <<<
[junit4] > Throwable #1: java.lang.AssertionError: Geo3dShape{planetmodel=PlanetModel: {ab=1.0011188180710464, c=0.9977622539852008}, shape=GeoPath: {planetmodel=PlanetModel: {ab=1.0011188180710464, c=0.9977622539852008}, width=1.53588974175501(87.99999999999999),
[junit4] > Throwable #1: java.lang.AssertionError: Geo3dShape{planetmodel=PlanetModel: {xyScaling=1.0011188180710464, zScaling=0.9977622539852008}, shape=GeoPath: {planetmodel=PlanetModel: {xyScaling=1.0011188180710464, zScaling=0.9977622539852008}, width=1.53588974175501(87.99999999999999),
points={[[X=0.12097657665150223, Y=-0.6754177666095532, Z=0.7265376136709238], [X=-0.3837892785614207, Y=0.4258049113530899, Z=0.8180007850434892]]}}} intersect Pt(x=23.81626064835212,y=16.0)
[junit4] > at __randomizedtesting.SeedInfo.seed([2595268DA3F13FEA:6CC30D8C83453E5D]:0)
[junit4] > at org.apache.lucene.spatial.spatial4j.RandomizedShapeTestCase._assertIntersect(RandomizedShapeTestCase.java:168)

252
lucene/spatial3d/src/java/org/apache/lucene/spatial3d/Geo3DDocValuesField.java
View File

@ -48,43 +48,14 @@ import org.apache.lucene.spatial3d.geom.GeoOutsideDistance;
* @see Geo3DPoint
*/
public class Geo3DDocValuesField extends Field {
// These are the multiplicative constants we need to use to arrive at values that fit in 21 bits.
// The formula we use to go from double to encoded value is: Math.floor((value - minimum) * factor + 0.5)
// If we plug in maximum for value, we should get 0x1FFFFF.
// So, 0x1FFFFF = Math.floor((maximum - minimum) * factor + 0.5)
// We factor out the 0.5 and Math.floor by stating instead:
// 0x1FFFFF = (maximum - minimum) * factor
// So, factor = 0x1FFFFF / (maximum - minimum)
private final static double inverseMaximumValue = 1.0 / (double)(0x1FFFFF);
private final static double inverseXFactor = (PlanetModel.WGS84.getMaximumXValue() - PlanetModel.WGS84.getMinimumXValue()) * inverseMaximumValue;
private final static double inverseYFactor = (PlanetModel.WGS84.getMaximumYValue() - PlanetModel.WGS84.getMinimumYValue()) * inverseMaximumValue;
private final static double inverseZFactor = (PlanetModel.WGS84.getMaximumZValue() - PlanetModel.WGS84.getMinimumZValue()) * inverseMaximumValue;
private final static double xFactor = 1.0 / inverseXFactor;
private final static double yFactor = 1.0 / inverseYFactor;
private final static double zFactor = 1.0 / inverseZFactor;
// Fudge factor for step adjustments. This is here solely to handle inaccuracies in bounding boxes
// that occur because of quantization. For unknown reasons, the fudge factor needs to be
// 10.0 rather than 1.0. See LUCENE-7430.
private final static double STEP_FUDGE = 10.0;
// These values are the delta between a value and the next value in each specific dimension
private final static double xStep = inverseXFactor * STEP_FUDGE;
private final static double yStep = inverseYFactor * STEP_FUDGE;
private final static double zStep = inverseZFactor * STEP_FUDGE;
private final PlanetModel planetModel;
/**
* Type for a Geo3DDocValuesField
* <p>
* Each value stores a 64-bit long where the three values (x, y, and z) are given
* 21 bits each. Each 21-bit value represents the maximum extent in that dimension
* for the WGS84 planet model.
* for the defined planet model.
*/
public static final FieldType TYPE = new FieldType();
static {
@ -98,8 +69,8 @@ public class Geo3DDocValuesField extends Field {
* @param point is the point.
* @throws IllegalArgumentException if the field name is null or the point is out of bounds
*/
public Geo3DDocValuesField(final String name, final GeoPoint point) {
super(name, TYPE);
public Geo3DDocValuesField(final String name, final GeoPoint point, final PlanetModel planetModel) {
this(name, TYPE, planetModel);
setLocationValue(point);
}
@ -111,18 +82,23 @@ public class Geo3DDocValuesField extends Field {
* @param z is the z value for the point.
* @throws IllegalArgumentException if the field name is null or x, y, or z are out of bounds
*/
public Geo3DDocValuesField(final String name, final double x, final double y, final double z) {
super(name, TYPE);
public Geo3DDocValuesField(final String name, final double x, final double y, final double z, final PlanetModel planetModel) {
this(name, TYPE, planetModel);
setLocationValue(x, y, z);
}
private Geo3DDocValuesField(final String name, final FieldType type, final PlanetModel planetModel) {
super(name, TYPE);
this.planetModel = planetModel;
}
/**
* Change the values of this field
* @param point is the point.
* @throws IllegalArgumentException if the point is out of bounds
*/
public void setLocationValue(final GeoPoint point) {
fieldsData = Long.valueOf(encodePoint(point));
fieldsData = Long.valueOf(planetModel.getDocValueEncoder().encodePoint(point));
}
/**
@ -133,158 +109,7 @@ public class Geo3DDocValuesField extends Field {
* @throws IllegalArgumentException if x, y, or z are out of bounds
*/
public void setLocationValue(final double x, final double y, final double z) {
fieldsData = Long.valueOf(encodePoint(x, y, z));
}
/** Encode a point.
* @param point is the point
* @return the encoded long
*/
public static long encodePoint(final GeoPoint point) {
return encodePoint(point.x, point.y, point.z);
}
/** Encode a point.
* @param x is the x value
* @param y is the y value
* @param z is the z value
* @return the encoded long
*/
public static long encodePoint(final double x, final double y, final double z) {
int XEncoded = encodeX(x);
int YEncoded = encodeY(y);
int ZEncoded = encodeZ(z);
return
(((long)(XEncoded & 0x1FFFFF)) << 42) |
(((long)(YEncoded & 0x1FFFFF)) << 21) |
((long)(ZEncoded & 0x1FFFFF));
}
/** Decode GeoPoint value from long docvalues value.
* @param docValue is the doc values value.
* @return the GeoPoint.
*/
public static GeoPoint decodePoint(final long docValue) {
return new GeoPoint(decodeX(((int)(docValue >> 42)) & 0x1FFFFF),
decodeY(((int)(docValue >> 21)) & 0x1FFFFF),
decodeZ(((int)(docValue)) & 0x1FFFFF));
}
/** Decode X value from long docvalues value.
* @param docValue is the doc values value.
* @return the x value.
*/
public static double decodeXValue(final long docValue) {
return decodeX(((int)(docValue >> 42)) & 0x1FFFFF);
}
/** Decode Y value from long docvalues value.
* @param docValue is the doc values value.
* @return the y value.
*/
public static double decodeYValue(final long docValue) {
return decodeY(((int)(docValue >> 21)) & 0x1FFFFF);
}
/** Decode Z value from long docvalues value.
* @param docValue is the doc values value.
* @return the z value.
*/
public static double decodeZValue(final long docValue) {
return decodeZ(((int)(docValue)) & 0x1FFFFF);
}
/** Round the provided X value down, by encoding it, decrementing it, and unencoding it.
* @param startValue is the starting value.
* @return the rounded value.
*/
public static double roundDownX(final double startValue) {
return startValue - xStep;
}
/** Round the provided X value up, by encoding it, incrementing it, and unencoding it.
* @param startValue is the starting value.
* @return the rounded value.
*/
public static double roundUpX(final double startValue) {
return startValue + xStep;
}
/** Round the provided Y value down, by encoding it, decrementing it, and unencoding it.
* @param startValue is the starting value.
* @return the rounded value.
*/
public static double roundDownY(final double startValue) {
return startValue - yStep;
}
/** Round the provided Y value up, by encoding it, incrementing it, and unencoding it.
* @param startValue is the starting value.
* @return the rounded value.
*/
public static double roundUpY(final double startValue) {
return startValue + yStep;
}
/** Round the provided Z value down, by encoding it, decrementing it, and unencoding it.
* @param startValue is the starting value.
* @return the rounded value.
*/
public static double roundDownZ(final double startValue) {
return startValue - zStep;
}
/** Round the provided Z value up, by encoding it, incrementing it, and unencoding it.
* @param startValue is the starting value.
* @return the rounded value.
*/
public static double roundUpZ(final double startValue) {
return startValue + zStep;
}
// For encoding/decoding, we generally want the following behavior:
// (1) If you encode the maximum value or the minimum value, the resulting int fits in 21 bits.
// (2) If you decode an encoded value, you get back the original value for both the minimum and maximum planet model values.
// (3) Rounding occurs such that a small delta from the minimum and maximum planet model values still returns the same
// values -- that is, these are in the center of the range of input values that should return the minimum or maximum when decoded
private static int encodeX(final double x) {
if (x > PlanetModel.WGS84.getMaximumXValue()) {
throw new IllegalArgumentException("x value exceeds WGS84 maximum");
} else if (x < PlanetModel.WGS84.getMinimumXValue()) {
throw new IllegalArgumentException("x value less than WGS84 minimum");
}
return (int)Math.floor((x - PlanetModel.WGS84.getMinimumXValue()) * xFactor + 0.5);
}
private static double decodeX(final int x) {
return x * inverseXFactor + PlanetModel.WGS84.getMinimumXValue();
}
private static int encodeY(final double y) {
if (y > PlanetModel.WGS84.getMaximumYValue()) {
throw new IllegalArgumentException("y value exceeds WGS84 maximum");
} else if (y < PlanetModel.WGS84.getMinimumYValue()) {
throw new IllegalArgumentException("y value less than WGS84 minimum");
}
return (int)Math.floor((y - PlanetModel.WGS84.getMinimumYValue()) * yFactor + 0.5);
}
private static double decodeY(final int y) {
return y * inverseYFactor + PlanetModel.WGS84.getMinimumYValue();
}
private static int encodeZ(final double z) {
if (z > PlanetModel.WGS84.getMaximumZValue()) {
throw new IllegalArgumentException("z value exceeds WGS84 maximum");
} else if (z < PlanetModel.WGS84.getMinimumZValue()) {
throw new IllegalArgumentException("z value less than WGS84 minimum");
}
return (int)Math.floor((z - PlanetModel.WGS84.getMinimumZValue()) * zFactor + 0.5);
}
private static double decodeZ(final int z) {
return z * inverseZFactor + PlanetModel.WGS84.getMinimumZValue();
fieldsData = Long.valueOf(planetModel.getDocValueEncoder().encodePoint(x, y, z));
}
/** helper: checks a fieldinfo and throws exception if its definitely not a Geo3DDocValuesField */
@ -307,11 +132,11 @@ public class Geo3DDocValuesField extends Field {
long currentValue = (Long)fieldsData;
result.append(decodeXValue(currentValue));
result.append(planetModel.getDocValueEncoder().decodeXValue(currentValue));
result.append(',');
result.append(decodeYValue(currentValue));
result.append(planetModel.getDocValueEncoder().decodeYValue(currentValue));
result.append(',');
result.append(decodeZValue(currentValue));
result.append(planetModel.getDocValueEncoder().decodeZValue(currentValue));
result.append('>');
return result.toString();
@ -335,9 +160,9 @@ public class Geo3DDocValuesField extends Field {
* @return SortField ordering documents by distance
* @throws IllegalArgumentException if {@code field} is null or circle has invalid coordinates.
*/
public static SortField newDistanceSort(final String field, final double latitude, final double longitude, final double maxRadiusMeters) {
final GeoDistanceShape shape = Geo3DUtil.fromDistance(latitude, longitude, maxRadiusMeters);
return new Geo3DPointSortField(field, shape);
public static SortField newDistanceSort(final String field, final double latitude, final double longitude, final double maxRadiusMeters, final PlanetModel planetModel) {
final GeoDistanceShape shape = Geo3DUtil.fromDistance(planetModel, latitude, longitude, maxRadiusMeters);
return new Geo3DPointSortField(field, planetModel, shape);
}
/**
@ -358,9 +183,9 @@ public class Geo3DDocValuesField extends Field {
* @return SortField ordering documents by distance
* @throws IllegalArgumentException if {@code field} is null or path has invalid coordinates.
*/
public static SortField newPathSort(final String field, final double[] pathLatitudes, final double[] pathLongitudes, final double pathWidthMeters) {
final GeoDistanceShape shape = Geo3DUtil.fromPath(pathLatitudes, pathLongitudes, pathWidthMeters);
return new Geo3DPointSortField(field, shape);
public static SortField newPathSort(final String field, final double[] pathLatitudes, final double[] pathLongitudes, final double pathWidthMeters, final PlanetModel planetModel) {
final GeoDistanceShape shape = Geo3DUtil.fromPath(planetModel, pathLatitudes, pathLongitudes, pathWidthMeters);
return new Geo3DPointSortField(field, planetModel, shape);
}
// Outside distances
@ -384,9 +209,9 @@ public class Geo3DDocValuesField extends Field {
* @return SortField ordering documents by distance
* @throws IllegalArgumentException if {@code field} is null or location has invalid coordinates.
*/
public static SortField newOutsideDistanceSort(final String field, final double latitude, final double longitude, final double maxRadiusMeters) {
final GeoOutsideDistance shape = Geo3DUtil.fromDistance(latitude, longitude, maxRadiusMeters);
return new Geo3DPointOutsideSortField(field, shape);
public static SortField newOutsideDistanceSort(final String field, final double latitude, final double longitude, final double maxRadiusMeters, final PlanetModel planetModel) {
final GeoOutsideDistance shape = Geo3DUtil.fromDistance(planetModel, latitude, longitude, maxRadiusMeters);
return new Geo3DPointOutsideSortField(field, planetModel, shape);
}
/**
@ -409,9 +234,9 @@ public class Geo3DDocValuesField extends Field {
* @return SortField ordering documents by distance
* @throws IllegalArgumentException if {@code field} is null or box has invalid coordinates.
*/
public static SortField newOutsideBoxSort(final String field, final double minLatitude, final double maxLatitude, final double minLongitude, final double maxLongitude) {
final GeoOutsideDistance shape = Geo3DUtil.fromBox(minLatitude, maxLatitude, minLongitude, maxLongitude);
return new Geo3DPointOutsideSortField(field, shape);
public static SortField newOutsideBoxSort(final String field, final double minLatitude, final double maxLatitude, final double minLongitude, final double maxLongitude, final PlanetModel planetModel) {
final GeoOutsideDistance shape = Geo3DUtil.fromBox(planetModel, minLatitude, maxLatitude, minLongitude, maxLongitude);
return new Geo3DPointOutsideSortField(field, planetModel, shape);
}
/**
@ -431,9 +256,9 @@ public class Geo3DDocValuesField extends Field {
* @return SortField ordering documents by distance
* @throws IllegalArgumentException if {@code field} is null or polygon has invalid coordinates.
*/
public static SortField newOutsidePolygonSort(final String field, final Polygon... polygons) {
final GeoOutsideDistance shape = Geo3DUtil.fromPolygon(polygons);
return new Geo3DPointOutsideSortField(field, shape);
public static SortField newOutsidePolygonSort(final String field, final PlanetModel planetModel, final Polygon... polygons) {
final GeoOutsideDistance shape = Geo3DUtil.fromPolygon(planetModel, polygons);
return new Geo3DPointOutsideSortField(field, planetModel, shape);
}
/**
@ -454,9 +279,9 @@ public class Geo3DDocValuesField extends Field {
* @return SortField ordering documents by distance
* @throws IllegalArgumentException if {@code field} is null or polygon has invalid coordinates.
*/
public static SortField newOutsideLargePolygonSort(final String field, final Polygon... polygons) {
final GeoOutsideDistance shape = Geo3DUtil.fromLargePolygon(polygons);
return new Geo3DPointOutsideSortField(field, shape);
public static SortField newOutsideLargePolygonSort(final String field, final PlanetModel planetModel, final Polygon... polygons) {
final GeoOutsideDistance shape = Geo3DUtil.fromLargePolygon(planetModel, polygons);
return new Geo3DPointOutsideSortField(field, planetModel, shape);
}
/**
@ -478,9 +303,8 @@ public class Geo3DDocValuesField extends Field {
* @return SortField ordering documents by distance
* @throws IllegalArgumentException if {@code field} is null or path has invalid coordinates.
*/
public static SortField newOutsidePathSort(final String field, final double[] pathLatitudes, final double[] pathLongitudes, final double pathWidthMeters) {
final GeoOutsideDistance shape = Geo3DUtil.fromPath(pathLatitudes, pathLongitudes, pathWidthMeters);
return new Geo3DPointOutsideSortField(field, shape);
public static SortField newOutsidePathSort(final String field, final double[] pathLatitudes, final double[] pathLongitudes, final double pathWidthMeters, final PlanetModel planetModel) {
final GeoOutsideDistance shape = Geo3DUtil.fromPath(planetModel, pathLatitudes, pathLongitudes, pathWidthMeters);
return new Geo3DPointOutsideSortField(field, planetModel, shape);
}
}

99
lucene/spatial3d/src/java/org/apache/lucene/spatial3d/Geo3DPoint.java
View File

@ -41,6 +41,9 @@ import org.apache.lucene.util.NumericUtils;
* @lucene.experimental */
public final class Geo3DPoint extends Field {
/** Planet Model for this Geo3DPoint */
protected final PlanetModel planetModel;
/** Indexing {@link FieldType}. */
public static final FieldType TYPE = new FieldType();
static {
@ -48,24 +51,54 @@ public final class Geo3DPoint extends Field {
TYPE.freeze();
}
/**
* Creates a new Geo3DPoint field with the specified latitude, longitude (in degrees).
/**
* Creates a new Geo3DPoint field with the specified latitude, longitude (in degrees), with default WGS84 PlanetModel.
*
* @throws IllegalArgumentException if the field name is null or latitude or longitude are out of bounds
*/
public Geo3DPoint(String name, double latitude, double longitude) {
public Geo3DPoint(String name, double lat, double lon) {
this(name, PlanetModel.WGS84, lat, lon);
}
/**
* Creates a new Geo3DPoint field with the specified x,y,z, using default WGS84 planet model.
*
* @throws IllegalArgumentException if the field name is null or latitude or longitude are out of bounds
*/
public Geo3DPoint(String name, double x, double y, double z) {
this(name, PlanetModel.WGS84, x, y, z);
}
/**
* Creates a new Geo3DPoint field with the specified x,y,z, and given planet model.
*
* @throws IllegalArgumentException if the field name is null or latitude or longitude are out of bounds
*/
public Geo3DPoint(String name, PlanetModel planetModel, double x, double y, double z) {
super(name, TYPE);
this.planetModel = planetModel;
fillFieldsData(planetModel, x, y, z);
}
/**
* Creates a new Geo3DPoint field with the specified latitude, longitude (in degrees), given a planet model.
*
* @throws IllegalArgumentException if the field name is null or latitude or longitude are out of bounds
*/
public Geo3DPoint(String name, PlanetModel planetModel, double latitude, double longitude) {
super(name, TYPE);
GeoUtils.checkLatitude(latitude);
GeoUtils.checkLongitude(longitude);
this.planetModel = planetModel;
// Translate latitude/longitude to x,y,z:
final GeoPoint point = new GeoPoint(PlanetModel.WGS84, Geo3DUtil.fromDegrees(latitude), Geo3DUtil.fromDegrees(longitude));
fillFieldsData(point.x, point.y, point.z);
final GeoPoint point = new GeoPoint(planetModel, Geo3DUtil.fromDegrees(latitude), Geo3DUtil.fromDegrees(longitude));
fillFieldsData(planetModel, point.x, point.y, point.z);
}
/**
* Create a query for matching points within the specified distance of the supplied location.
* @param field field name. must not be null. Note that because
* {@link PlanetModel#WGS84} is used, this query is approximate and may have up
* @param field field name. must not be null. Note that if
* {@link PlanetModel#WGS84} is used, the query is approximate and may have up
* to 0.5% error.
*
* @param latitude latitude at the center: must be within standard +/-90 coordinate bounds.
@ -74,8 +107,8 @@ public final class Geo3DPoint extends Field {
* @return query matching points within this distance
* @throws IllegalArgumentException if {@code field} is null, location has invalid coordinates, or radius is invalid.
*/
public static Query newDistanceQuery(final String field, final double latitude, final double longitude, final double radiusMeters) {
final GeoShape shape = Geo3DUtil.fromDistance(latitude, longitude, radiusMeters);
public static Query newDistanceQuery(final String field, final PlanetModel planetModel, final double latitude, final double longitude, final double radiusMeters) {
final GeoShape shape = Geo3DUtil.fromDistance(planetModel, latitude, longitude, radiusMeters);
return newShapeQuery(field, shape);
}
@ -91,8 +124,8 @@ public final class Geo3DPoint extends Field {
* @return query matching points within this box
* @throws IllegalArgumentException if {@code field} is null, or the box has invalid coordinates.
*/
public static Query newBoxQuery(final String field, final double minLatitude, final double maxLatitude, final double minLongitude, final double maxLongitude) {
final GeoShape shape = Geo3DUtil.fromBox(minLatitude, maxLatitude, minLongitude, maxLongitude);
public static Query newBoxQuery(final String field, final PlanetModel planetModel, final double minLatitude, final double maxLatitude, final double minLongitude, final double maxLongitude) {
final GeoShape shape = Geo3DUtil.fromBox(planetModel, minLatitude, maxLatitude, minLongitude, maxLongitude);
return newShapeQuery(field, shape);
}
@ -105,8 +138,8 @@ public final class Geo3DPoint extends Field {
* @param polygons is the list of polygons to use to construct the query; must be at least one.
* @return query matching points within this polygon
*/
public static Query newPolygonQuery(final String field, final Polygon... polygons) {
final GeoShape shape = Geo3DUtil.fromPolygon(polygons);
public static Query newPolygonQuery(final String field, final PlanetModel planetModel, final Polygon... polygons) {
final GeoShape shape = Geo3DUtil.fromPolygon(planetModel, polygons);
return newShapeQuery(field, shape);
}
@ -119,8 +152,8 @@ public final class Geo3DPoint extends Field {
* @param polygons is the list of polygons to use to construct the query; must be at least one.
* @return query matching points within this polygon
*/
public static Query newLargePolygonQuery(final String field, final Polygon... polygons) {
final GeoShape shape = Geo3DUtil.fromLargePolygon(polygons);
public static Query newLargePolygonQuery(final String field, PlanetModel planetModel, final Polygon... polygons) {
final GeoShape shape = Geo3DUtil.fromLargePolygon(planetModel, polygons);
return newShapeQuery(field, shape);
}
@ -133,39 +166,29 @@ public final class Geo3DPoint extends Field {
* @param pathWidthMeters width of the path in meters.
* @return query matching points within this polygon
*/
public static Query newPathQuery(final String field, final double[] pathLatitudes, final double[] pathLongitudes, final double pathWidthMeters) {
final GeoShape shape = Geo3DUtil.fromPath(pathLatitudes, pathLongitudes, pathWidthMeters);
public static Query newPathQuery(final String field, final double[] pathLatitudes, final double[] pathLongitudes, final double pathWidthMeters, PlanetModel planetModel) {
final GeoShape shape = Geo3DUtil.fromPath(planetModel, pathLatitudes, pathLongitudes, pathWidthMeters);
return newShapeQuery(field, shape);
}
/**
* Creates a new Geo3DPoint field with the specified x,y,z.
*
* @throws IllegalArgumentException if the field name is null or latitude or longitude are out of bounds
*/
public Geo3DPoint(String name, double x, double y, double z) {
super(name, TYPE);
fillFieldsData(x, y, z);
}
private void fillFieldsData(double x, double y, double z) {
private void fillFieldsData(PlanetModel planetModel, double x, double y, double z) {
byte[] bytes = new byte[12];
encodeDimension(x, bytes, 0);
encodeDimension(y, bytes, Integer.BYTES);
encodeDimension(z, bytes, 2*Integer.BYTES);
encodeDimension(x, bytes, 0, planetModel);
encodeDimension(y, bytes, Integer.BYTES, planetModel);
encodeDimension(z, bytes, 2*Integer.BYTES, planetModel);
fieldsData = new BytesRef(bytes);
}
// public helper methods (e.g. for queries)
/** Encode single dimension */
public static void encodeDimension(double value, byte bytes[], int offset) {
NumericUtils.intToSortableBytes(Geo3DUtil.encodeValue(value), bytes, offset);
public static void encodeDimension(double value, byte bytes[], int offset, PlanetModel planetModel) {
NumericUtils.intToSortableBytes(planetModel.encodeValue(value), bytes, offset);
}
/** Decode single dimension */
public static double decodeDimension(byte value[], int offset) {
return Geo3DUtil.decodeValue(NumericUtils.sortableBytesToInt(value, offset));
public static double decodeDimension(byte value[], int offset, PlanetModel planetModel) {
return planetModel.decodeValue(NumericUtils.sortableBytesToInt(value, offset));
}
/** Returns a query matching all points inside the provided shape.
@ -186,9 +209,9 @@ public final class Geo3DPoint extends Field {
result.append(':');
BytesRef bytes = (BytesRef) fieldsData;
result.append(" x=").append(decodeDimension(bytes.bytes, bytes.offset));
result.append(" y=").append(decodeDimension(bytes.bytes, bytes.offset + Integer.BYTES));
result.append(" z=").append(decodeDimension(bytes.bytes, bytes.offset + 2 * Integer.BYTES));
result.append(" x=").append(decodeDimension(bytes.bytes, bytes.offset, this.planetModel));
result.append(" y=").append(decodeDimension(bytes.bytes, bytes.offset + Integer.BYTES, this.planetModel));
result.append(" z=").append(decodeDimension(bytes.bytes, bytes.offset + 2 * Integer.BYTES, this.planetModel));
result.append('>');
return result.toString();
}

18
lucene/spatial3d/src/java/org/apache/lucene/spatial3d/Geo3DPointDistanceComparator.java
View File

@ -42,6 +42,7 @@ class Geo3DPointDistanceComparator extends FieldComparator<Double> implements Le
final String field;
final GeoDistanceShape distanceShape;
final private PlanetModel planetModel;
final double[] values;
double bottomDistance;
@ -53,9 +54,10 @@ class Geo3DPointDistanceComparator extends FieldComparator<Double> implements Le
// the number of times setBottom has been called (adversary protection)
int setBottomCounter = 0;
public Geo3DPointDistanceComparator(String field, final GeoDistanceShape distanceShape, int numHits) {
public Geo3DPointDistanceComparator(String field, final PlanetModel planetModel, final GeoDistanceShape distanceShape, int numHits) {
this.field = field;
this.distanceShape = distanceShape;
this.planetModel = planetModel;
this.values = new double[numHits];
}
@ -105,9 +107,9 @@ class Geo3DPointDistanceComparator extends FieldComparator<Double> implements Le
// Test against bounds.
// First we need to decode...
final double x = Geo3DDocValuesField.decodeXValue(encoded);
final double y = Geo3DDocValuesField.decodeYValue(encoded);
final double z = Geo3DDocValuesField.decodeZValue(encoded);
final double x = planetModel.getDocValueEncoder().decodeXValue(encoded);
final double y = planetModel.getDocValueEncoder().decodeYValue(encoded);
final double z = planetModel.getDocValueEncoder().decodeZValue(encoded);
if (x > priorityQueueBounds.getMaximumX() ||
x < priorityQueueBounds.getMinimumX() ||
@ -142,7 +144,7 @@ class Geo3DPointDistanceComparator extends FieldComparator<Double> implements Le
@Override
public Double value(int slot) {
// Return the arc distance
return Double.valueOf(values[slot] * PlanetModel.WGS84_MEAN);
return Double.valueOf(values[slot] * planetModel.getMeanRadiusMeters());
}
@Override
@ -160,9 +162,9 @@ class Geo3DPointDistanceComparator extends FieldComparator<Double> implements Le
for (int i = 0; i < numValues; i++) {
final long encoded = currentDocs.nextValue();
final double distance = distanceShape.computeDistance(DistanceStyle.ARC,
Geo3DDocValuesField.decodeXValue(encoded),
Geo3DDocValuesField.decodeYValue(encoded),
Geo3DDocValuesField.decodeZValue(encoded));
planetModel.getDocValueEncoder().decodeXValue(encoded),
planetModel.getDocValueEncoder().decodeYValue(encoded),
planetModel.getDocValueEncoder().decodeZValue(encoded));
minValue = Math.min(minValue, distance);
}
}

18
lucene/spatial3d/src/java/org/apache/lucene/spatial3d/Geo3DPointOutsideDistanceComparator.java
View File

@ -37,14 +37,16 @@ class Geo3DPointOutsideDistanceComparator extends FieldComparator<Double> implem
final String field;
final GeoOutsideDistance distanceShape;
final private PlanetModel planetModel;
final double[] values;
double bottomDistance;
double topValue;
SortedNumericDocValues currentDocs;
public Geo3DPointOutsideDistanceComparator(String field, final GeoOutsideDistance distanceShape, int numHits) {
public Geo3DPointOutsideDistanceComparator(String field, final PlanetModel planetModel, final GeoOutsideDistance distanceShape, int numHits) {
this.field = field;
this.planetModel = planetModel;
this.distanceShape = distanceShape;
this.values = new double[numHits];
}
@ -85,9 +87,9 @@ class Geo3DPointOutsideDistanceComparator extends FieldComparator<Double> implem
// Test against bounds.
// First we need to decode...
final double x = Geo3DDocValuesField.decodeXValue(encoded);
final double y = Geo3DDocValuesField.decodeYValue(encoded);
final double z = Geo3DDocValuesField.decodeZValue(encoded);
final double x = planetModel.getDocValueEncoder().decodeXValue(encoded);
final double y = planetModel.getDocValueEncoder().decodeYValue(encoded);
final double z = planetModel.getDocValueEncoder().decodeZValue(encoded);
cmp = Math.max(cmp, Double.compare(bottomDistance, distanceShape.computeOutsideDistance(DistanceStyle.ARC, x, y, z)));
}
@ -113,7 +115,7 @@ class Geo3DPointOutsideDistanceComparator extends FieldComparator<Double> implem
@Override
public Double value(int slot) {
// Return the arc distance
return Double.valueOf(values[slot] * PlanetModel.WGS84_MEAN);
return Double.valueOf(values[slot] * planetModel.getMeanRadiusMeters());
}
@Override
@ -131,9 +133,9 @@ class Geo3DPointOutsideDistanceComparator extends FieldComparator<Double> implem
for (int i = 0; i < numValues; i++) {
final long encoded = currentDocs.nextValue();
final double distance = distanceShape.computeOutsideDistance(DistanceStyle.ARC,
Geo3DDocValuesField.decodeXValue(encoded),
Geo3DDocValuesField.decodeYValue(encoded),
Geo3DDocValuesField.decodeZValue(encoded));
planetModel.getDocValueEncoder().decodeXValue(encoded),
planetModel.getDocValueEncoder().decodeYValue(encoded),
planetModel.getDocValueEncoder().decodeZValue(encoded));
minValue = Math.min(minValue, distance);
}
}

7
lucene/spatial3d/src/java/org/apache/lucene/spatial3d/Geo3DPointOutsideSortField.java
View File

@ -19,14 +19,16 @@ package org.apache.lucene.spatial3d;
import org.apache.lucene.search.FieldComparator;
import org.apache.lucene.search.SortField;
import org.apache.lucene.spatial3d.geom.GeoOutsideDistance;
import org.apache.lucene.spatial3d.geom.PlanetModel;
/**
* Sorts by outside distance from an origin location.
*/
final class Geo3DPointOutsideSortField extends SortField {
final GeoOutsideDistance distanceShape;
final private PlanetModel planetModel;
Geo3DPointOutsideSortField(final String field, final GeoOutsideDistance distanceShape) {
Geo3DPointOutsideSortField(final String field, final PlanetModel planetModel, final GeoOutsideDistance distanceShape) {
super(field, SortField.Type.CUSTOM);
if (field == null) {
throw new IllegalArgumentException("field must not be null");
@ -34,13 +36,14 @@ final class Geo3DPointOutsideSortField extends SortField {
if (distanceShape == null) {
throw new IllegalArgumentException("distanceShape must not be null");
}
this.planetModel = planetModel;
this.distanceShape = distanceShape;
setMissingValue(Double.POSITIVE_INFINITY);
}
@Override
public FieldComparator<?> getComparator(int numHits, int sortPos) {
return new Geo3DPointOutsideDistanceComparator(getField(), distanceShape, numHits);
return new Geo3DPointOutsideDistanceComparator(getField(), planetModel, distanceShape, numHits);
}
@Override

7
lucene/spatial3d/src/java/org/apache/lucene/spatial3d/Geo3DPointSortField.java
View File

@ -19,14 +19,16 @@ package org.apache.lucene.spatial3d;
import org.apache.lucene.search.FieldComparator;
import org.apache.lucene.search.SortField;
import org.apache.lucene.spatial3d.geom.GeoDistanceShape;
import org.apache.lucene.spatial3d.geom.PlanetModel;
/**
* Sorts by distance from an origin location.
*/
final class Geo3DPointSortField extends SortField {
final GeoDistanceShape distanceShape;
final PlanetModel planetModel;
Geo3DPointSortField(final String field, final GeoDistanceShape distanceShape) {
Geo3DPointSortField(final String field, final PlanetModel planetModel, final GeoDistanceShape distanceShape) {
super(field, SortField.Type.CUSTOM);
if (field == null) {
throw new IllegalArgumentException("field must not be null");
@ -35,12 +37,13 @@ final class Geo3DPointSortField extends SortField {
throw new IllegalArgumentException("distanceShape must not be null");
}
this.distanceShape = distanceShape;
this.planetModel = planetModel;
setMissingValue(Double.POSITIVE_INFINITY);
}
@Override
public FieldComparator<?> getComparator(int numHits, int sortPos) {
return new Geo3DPointDistanceComparator(getField(), distanceShape, numHits);
return new Geo3DPointDistanceComparator(getField(), planetModel, distanceShape, numHits);
}
@Override

126
lucene/spatial3d/src/java/org/apache/lucene/spatial3d/Geo3DUtil.java
View File

@ -36,119 +36,57 @@ import java.util.ArrayList;
class Geo3DUtil {
/** How many radians are in one earth surface meter */
final static double RADIANS_PER_METER = 1.0 / PlanetModel.WGS84_MEAN;
/** How many radians are in one degree */
final static double RADIANS_PER_DEGREE = Math.PI / 180.0;
private static final double MAX_VALUE = PlanetModel.WGS84.getMaximumMagnitude();
private static final int BITS = 32;
private static final double MUL = (0x1L<<BITS)/(2*MAX_VALUE);
static final double DECODE = getNextSafeDouble(1/MUL);
static final int MIN_ENCODED_VALUE = encodeValue(-MAX_VALUE);
static final int MAX_ENCODED_VALUE = encodeValue(MAX_VALUE);
public static int encodeValue(double x) {
if (x > MAX_VALUE) {
throw new IllegalArgumentException("value=" + x + " is out-of-bounds (greater than WGS84's planetMax=" + MAX_VALUE + ")");
}
if (x < -MAX_VALUE) {
throw new IllegalArgumentException("value=" + x + " is out-of-bounds (less than than WGS84's -planetMax=" + -MAX_VALUE + ")");
}
long result = (long) Math.floor(x / DECODE);
assert result >= Integer.MIN_VALUE;
assert result <= Integer.MAX_VALUE;
return (int) result;
}
public static double decodeValue(int x) {
double result;
if (x == MIN_ENCODED_VALUE) {
// We must special case this, because -MAX_VALUE is not guaranteed to land precisely at a floor value, and we don't ever want to
// return a value outside of the planet's range (I think?). The max value is "safe" because we floor during encode:
result = -MAX_VALUE;
} else if (x == MAX_ENCODED_VALUE) {
result = MAX_VALUE;
} else {
// We decode to the center value; this keeps the encoding stable
result = (x+0.5) * DECODE;
}
assert result >= -MAX_VALUE && result <= MAX_VALUE;
return result;
}
/** Returns smallest double that would encode to int x. */
// NOTE: keep this package private!!
static double decodeValueFloor(int x) {
assert x <= MAX_ENCODED_VALUE && x >= MIN_ENCODED_VALUE;
if (x == MIN_ENCODED_VALUE) {
return -MAX_VALUE;
// NOTE: remains in this class to keep method package private!!
static double decodeValueFloor(int x, PlanetModel planetModel) {
assert x <= planetModel.MAX_ENCODED_VALUE && x >= planetModel.MIN_ENCODED_VALUE;
if (x == planetModel.MIN_ENCODED_VALUE) {
return -planetModel.MAX_VALUE;
}
return x * DECODE;
}
/** Returns a double value >= x such that if you multiply that value by an int, and then
* divide it by that int again, you get precisely the same value back */
private static double getNextSafeDouble(double x) {
// Move to double space:
long bits = Double.doubleToLongBits(x);
// Make sure we are beyond the actual maximum value:
bits += Integer.MAX_VALUE;
// Clear the bottom 32 bits:
bits &= ~((long) Integer.MAX_VALUE);
// Convert back to double:
double result = Double.longBitsToDouble(bits);
assert result > x;
return result;
return x * planetModel.DECODE;
}
/** Returns largest double that would encode to int x. */
// NOTE: keep this package private!!
static double decodeValueCeil(int x) {
assert x <= MAX_ENCODED_VALUE && x >= MIN_ENCODED_VALUE;
if (x == MAX_ENCODED_VALUE) {
return MAX_VALUE;
static double decodeValueCeil(int x, PlanetModel planetModel) {
assert x <= planetModel.MAX_ENCODED_VALUE && x >= planetModel.MIN_ENCODED_VALUE;
if (x == planetModel.MAX_ENCODED_VALUE) {
return planetModel.MAX_VALUE;
}
return Math.nextDown((x+1) * DECODE);
return Math.nextDown((x+1) * planetModel.DECODE);
}
/** Converts degress to radians */
static double fromDegrees(final double degrees) {
return degrees * RADIANS_PER_DEGREE;
}
/** Converts earth-surface meters to radians */
static double fromMeters(final double meters) {
return meters * RADIANS_PER_METER;
}
/**
* Convert a set of Polygon objects into a GeoPolygon.
* @param polygons are the Polygon objects.
* @return the GeoPolygon.
*/
static GeoPolygon fromPolygon(final Polygon... polygons) {
static GeoPolygon fromPolygon(final PlanetModel planetModel, final Polygon... polygons) {
//System.err.println("Creating polygon...");
if (polygons.length < 1) {
throw new IllegalArgumentException("need at least one polygon");
}
final GeoPolygon shape;
if (polygons.length == 1) {
final GeoPolygon component = fromPolygon(polygons[0]);
final GeoPolygon component = fromPolygon(planetModel, polygons[0]);
if (component == null) {
// Polygon is degenerate
shape = new GeoCompositePolygon(PlanetModel.WGS84);
shape = new GeoCompositePolygon(planetModel);
} else {
shape = component;
}
} else {
final GeoCompositePolygon poly = new GeoCompositePolygon(PlanetModel.WGS84);
final GeoCompositePolygon poly = new GeoCompositePolygon(planetModel);
for (final Polygon p : polygons) {
final GeoPolygon component = fromPolygon(p);
final GeoPolygon component = fromPolygon(planetModel, p);
if (component != null) {
poly.addShape(component);
}
@ -165,11 +103,11 @@ class Geo3DUtil {
* @param polygons is the list of polygons to convert.
* @return the large GeoPolygon.
*/
static GeoPolygon fromLargePolygon(final Polygon... polygons) {
static GeoPolygon fromLargePolygon(final PlanetModel planetModel, final Polygon... polygons) {
if (polygons.length < 1) {
throw new IllegalArgumentException("need at least one polygon");
}
return GeoPolygonFactory.makeLargeGeoPolygon(PlanetModel.WGS84, convertToDescription(polygons));
return GeoPolygonFactory.makeLargeGeoPolygon(planetModel, convertToDescription(planetModel, polygons));
}
/**
@ -179,7 +117,7 @@ class Geo3DUtil {
* @param pathWidthMeters width of the path in meters.
* @return the path.
*/
static GeoPath fromPath(final double[] pathLatitudes, final double[] pathLongitudes, final double pathWidthMeters) {
static GeoPath fromPath(final PlanetModel planetModel, final double[] pathLatitudes, final double[] pathLongitudes, final double pathWidthMeters) {
if (pathLatitudes.length != pathLongitudes.length) {
throw new IllegalArgumentException("same number of latitudes and longitudes required");
}
@ -189,9 +127,9 @@ class Geo3DUtil {
final double longitude = pathLongitudes[i];
GeoUtils.checkLatitude(latitude);
GeoUtils.checkLongitude(longitude);
points[i] = new GeoPoint(PlanetModel.WGS84, fromDegrees(latitude), fromDegrees(longitude));
points[i] = new GeoPoint(planetModel, fromDegrees(latitude), fromDegrees(longitude));
}
return GeoPathFactory.makeGeoPath(PlanetModel.WGS84, fromMeters(pathWidthMeters), points);
return GeoPathFactory.makeGeoPath(planetModel, planetModel.fromMeters(pathWidthMeters), points);
}
/**
@ -201,10 +139,10 @@ class Geo3DUtil {
* @param radiusMeters maximum distance from the center in meters: must be non-negative and finite.
* @return the circle.
*/
static GeoCircle fromDistance(final double latitude, final double longitude, final double radiusMeters) {
static GeoCircle fromDistance(final PlanetModel planetModel, final double latitude, final double longitude, final double radiusMeters) {
GeoUtils.checkLatitude(latitude);
GeoUtils.checkLongitude(longitude);
return GeoCircleFactory.makeGeoCircle(PlanetModel.WGS84, fromDegrees(latitude), fromDegrees(longitude), fromMeters(radiusMeters));
return GeoCircleFactory.makeGeoCircle(planetModel, fromDegrees(latitude), fromDegrees(longitude), planetModel.fromMeters(radiusMeters));
}
/**
@ -215,12 +153,12 @@ class Geo3DUtil {
* @param maxLongitude longitude upper bound: must be within standard +/-180 coordinate bounds.
* @return the box.
*/
static GeoBBox fromBox(final double minLatitude, final double maxLatitude, final double minLongitude, final double maxLongitude) {
static GeoBBox fromBox(final PlanetModel planetModel, final double minLatitude, final double maxLatitude, final double minLongitude, final double maxLongitude) {
GeoUtils.checkLatitude(minLatitude);
GeoUtils.checkLongitude(minLongitude);
GeoUtils.checkLatitude(maxLatitude);
GeoUtils.checkLongitude(maxLongitude);
return GeoBBoxFactory.makeGeoBBox(PlanetModel.WGS84,
return GeoBBoxFactory.makeGeoBBox(planetModel,
Geo3DUtil.fromDegrees(maxLatitude), Geo3DUtil.fromDegrees(minLatitude), Geo3DUtil.fromDegrees(minLongitude), Geo3DUtil.fromDegrees(maxLongitude));
}
@ -230,14 +168,14 @@ class Geo3DUtil {
* @param polygon is the Polygon object.
* @return the GeoPolygon.
*/
private static GeoPolygon fromPolygon(final Polygon polygon) {
private static GeoPolygon fromPolygon(final PlanetModel planetModel, final Polygon polygon) {
// First, assemble the "holes". The geo3d convention is to use the same polygon sense on the inner ring as the
// outer ring, so we process these recursively with reverseMe flipped.
final Polygon[] theHoles = polygon.getHoles();
final List<GeoPolygon> holeList = new ArrayList<>(theHoles.length);
for (final Polygon hole : theHoles) {
//System.out.println("Hole: "+hole);
final GeoPolygon component = fromPolygon(hole);
final GeoPolygon component = fromPolygon(planetModel, hole);
if (component != null) {
holeList.add(component);
}
@ -252,10 +190,10 @@ class Geo3DUtil {
// We skip the last point anyway because the API requires it to be repeated, and geo3d doesn't repeat it.
for (int i = 0; i < polyLats.length - 1; i++) {
final int index = polyLats.length - 2 - i;
points.add(new GeoPoint(PlanetModel.WGS84, fromDegrees(polyLats[index]), fromDegrees(polyLons[index])));
points.add(new GeoPoint(planetModel, fromDegrees(polyLats[index]), fromDegrees(polyLons[index])));
}
//System.err.println(" building polygon with "+points.size()+" points...");
final GeoPolygon rval = GeoPolygonFactory.makeGeoPolygon(PlanetModel.WGS84, points, holeList);
final GeoPolygon rval = GeoPolygonFactory.makeGeoPolygon(planetModel, points, holeList);
//System.err.println(" ...done");
return rval;
}
@ -265,11 +203,11 @@ class Geo3DUtil {
* @param polygons is the list of polygons to convert.
* @return the list of polygon descriptions.
*/
private static List<GeoPolygonFactory.PolygonDescription> convertToDescription(final Polygon... polygons) {
private static List<GeoPolygonFactory.PolygonDescription> convertToDescription(final PlanetModel planetModel, final Polygon... polygons) {
final List<GeoPolygonFactory.PolygonDescription> descriptions = new ArrayList<>(polygons.length);
for (final Polygon polygon : polygons) {
final Polygon[] theHoles = polygon.getHoles();
final List<GeoPolygonFactory.PolygonDescription> holes = convertToDescription(theHoles);
final List<GeoPolygonFactory.PolygonDescription> holes = convertToDescription(planetModel, theHoles);
// Now do the polygon itself
final double[] polyLats = polygon.getPolyLats();
@ -280,7 +218,7 @@ class Geo3DUtil {
// We skip the last point anyway because the API requires it to be repeated, and geo3d doesn't repeat it.
for (int i = 0; i < polyLats.length - 1; i++) {
final int index = polyLats.length - 2 - i;
points.add(new GeoPoint(PlanetModel.WGS84, fromDegrees(polyLats[index]), fromDegrees(polyLons[index])));
points.add(new GeoPoint(planetModel, fromDegrees(polyLats[index]), fromDegrees(polyLons[index])));
}
descriptions.add(new GeoPolygonFactory.PolygonDescription(points, holes));

9
lucene/spatial3d/src/java/org/apache/lucene/spatial3d/PointInGeo3DShapeQuery.java
View File

@ -29,9 +29,7 @@ import org.apache.lucene.search.QueryVisitor;
import org.apache.lucene.search.ScoreMode;
import org.apache.lucene.search.Scorer;
import org.apache.lucene.search.Weight;
import org.apache.lucene.spatial3d.geom.BasePlanetObject;
import org.apache.lucene.spatial3d.geom.GeoShape;
import org.apache.lucene.spatial3d.geom.PlanetModel;
import org.apache.lucene.spatial3d.geom.XYZBounds;
import org.apache.lucene.util.Accountable;
import org.apache.lucene.util.DocIdSetBuilder;
@ -56,13 +54,6 @@ final class PointInGeo3DShapeQuery extends Query implements Accountable {
this.shape = shape;
this.shapeBounds = new XYZBounds();
shape.getBounds(shapeBounds);
if (shape instanceof BasePlanetObject) {
BasePlanetObject planetObject = (BasePlanetObject) shape;
if (planetObject.getPlanetModel().equals(PlanetModel.WGS84) == false) {
throw new IllegalArgumentException("this qurey requires PlanetModel.WGS84, but got: " + planetObject.getPlanetModel());
}
}
}
@Override

35
lucene/spatial3d/src/java/org/apache/lucene/spatial3d/PointInShapeIntersectVisitor.java
View File

@ -22,7 +22,7 @@ import org.apache.lucene.index.PointValues.Relation;
import org.apache.lucene.spatial3d.geom.GeoArea;
import org.apache.lucene.spatial3d.geom.GeoAreaFactory;
import org.apache.lucene.spatial3d.geom.GeoShape;
import org.apache.lucene.spatial3d.geom.PlanetModel;
import org.apache.lucene.spatial3d.geom.PlanetModel.DocValueEncoder;
import org.apache.lucene.spatial3d.geom.XYZBounds;
import org.apache.lucene.util.DocIdSetBuilder;
import org.apache.lucene.util.NumericUtils;
@ -43,12 +43,13 @@ class PointInShapeIntersectVisitor implements IntersectVisitor {
XYZBounds bounds) {
this.hits = hits;
this.shape = shape;
this.minimumX = Geo3DDocValuesField.roundDownX(bounds.getMinimumX());
this.maximumX = Geo3DDocValuesField.roundUpX(bounds.getMaximumX());
this.minimumY = Geo3DDocValuesField.roundDownY(bounds.getMinimumY());
this.maximumY = Geo3DDocValuesField.roundUpY(bounds.getMaximumY());
this.minimumZ = Geo3DDocValuesField.roundDownZ(bounds.getMinimumZ());
this.maximumZ = Geo3DDocValuesField.roundUpZ(bounds.getMaximumZ());
DocValueEncoder docValueEncoder = shape.getPlanetModel().getDocValueEncoder();
this.minimumX = docValueEncoder.roundDownX(bounds.getMinimumX());
this.maximumX = docValueEncoder.roundUpX(bounds.getMaximumX());
this.minimumY = docValueEncoder.roundDownY(bounds.getMinimumY());
this.maximumY = docValueEncoder.roundUpY(bounds.getMaximumY());
this.minimumZ = docValueEncoder.roundDownZ(bounds.getMinimumZ());
this.maximumZ = docValueEncoder.roundUpZ(bounds.getMaximumZ());
}
@Override
@ -64,9 +65,9 @@ class PointInShapeIntersectVisitor implements IntersectVisitor {
@Override
public void visit(int docID, byte[] packedValue) {
assert packedValue.length == 12;
double x = Geo3DPoint.decodeDimension(packedValue, 0);
double y = Geo3DPoint.decodeDimension(packedValue, Integer.BYTES);
double z = Geo3DPoint.decodeDimension(packedValue, 2 * Integer.BYTES);
double x = Geo3DPoint.decodeDimension(packedValue, 0, shape.getPlanetModel());
double y = Geo3DPoint.decodeDimension(packedValue, Integer.BYTES, shape.getPlanetModel());
double z = Geo3DPoint.decodeDimension(packedValue, 2 * Integer.BYTES, shape.getPlanetModel());
if (x >= minimumX && x <= maximumX &&
y >= minimumY && y <= maximumY &&
z >= minimumZ && z <= maximumZ) {
@ -82,12 +83,12 @@ class PointInShapeIntersectVisitor implements IntersectVisitor {
// here are inclusive, we need to extend the bounds to the largest un-quantized values that
// could quantize into these bounds. The encoding (Geo3DUtil.encodeValue) does
// a Math.round from double to long, so e.g. 1.4 -> 1, and -1.4 -> -1:
double xMin = Geo3DUtil.decodeValueFloor(NumericUtils.sortableBytesToInt(minPackedValue, 0));
double xMax = Geo3DUtil.decodeValueCeil(NumericUtils.sortableBytesToInt(maxPackedValue, 0));
double yMin = Geo3DUtil.decodeValueFloor(NumericUtils.sortableBytesToInt(minPackedValue, 1 * Integer.BYTES));
double yMax = Geo3DUtil.decodeValueCeil(NumericUtils.sortableBytesToInt(maxPackedValue, 1 * Integer.BYTES));
double zMin = Geo3DUtil.decodeValueFloor(NumericUtils.sortableBytesToInt(minPackedValue, 2 * Integer.BYTES));
double zMax = Geo3DUtil.decodeValueCeil(NumericUtils.sortableBytesToInt(maxPackedValue, 2 * Integer.BYTES));
double xMin = Geo3DUtil.decodeValueFloor(NumericUtils.sortableBytesToInt(minPackedValue, 0), shape.getPlanetModel());
double xMax = Geo3DUtil.decodeValueCeil(NumericUtils.sortableBytesToInt(maxPackedValue, 0), shape.getPlanetModel());
double yMin = Geo3DUtil.decodeValueFloor(NumericUtils.sortableBytesToInt(minPackedValue, 1 * Integer.BYTES), shape.getPlanetModel());
double yMax = Geo3DUtil.decodeValueCeil(NumericUtils.sortableBytesToInt(maxPackedValue, 1 * Integer.BYTES), shape.getPlanetModel());
double zMin = Geo3DUtil.decodeValueFloor(NumericUtils.sortableBytesToInt(minPackedValue, 2 * Integer.BYTES), shape.getPlanetModel());
double zMax = Geo3DUtil.decodeValueCeil(NumericUtils.sortableBytesToInt(maxPackedValue, 2 * Integer.BYTES), shape.getPlanetModel());
//System.out.println(" compare: x=" + cellXMin + "-" + cellXMax + " y=" + cellYMin + "-" + cellYMax + " z=" + cellZMin + "-" + cellZMax);
assert xMin <= xMax;
@ -102,7 +103,7 @@ class PointInShapeIntersectVisitor implements IntersectVisitor {
}
// Quick test failed so do slower one...
GeoArea xyzSolid = GeoAreaFactory.makeGeoArea(PlanetModel.WGS84, xMin, xMax, yMin, yMax, zMin, zMax);
GeoArea xyzSolid = GeoAreaFactory.makeGeoArea(shape.getPlanetModel(), xMin, xMax, yMin, yMax, zMin, zMax);
switch(xyzSolid.getRelationship(shape)) {
case GeoArea.CONTAINS:

26
lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/GeoComplexPolygon.java
View File

@ -300,8 +300,8 @@ class GeoComplexPolygon extends GeoBasePolygon {
if (testPointFixedYAbovePlane != null && testPointFixedYBelowPlane != null && fixedXAbovePlane != null && fixedXBelowPlane != null) {
//check if planes intersects inside world
final double checkAbove = 4.0 * (fixedXAbovePlane.D * fixedXAbovePlane.D * planetModel.inverseAbSquared + testPointFixedYAbovePlane.D * testPointFixedYAbovePlane.D * planetModel.inverseAbSquared - 1.0);
final double checkBelow = 4.0 * (fixedXBelowPlane.D * fixedXBelowPlane.D * planetModel.inverseAbSquared + testPointFixedYBelowPlane.D * testPointFixedYBelowPlane.D * planetModel.inverseAbSquared - 1.0);
final double checkAbove = 4.0 * (fixedXAbovePlane.D * fixedXAbovePlane.D * planetModel.inverseXYScalingSquared + testPointFixedYAbovePlane.D * testPointFixedYAbovePlane.D * planetModel.inverseXYScalingSquared - 1.0);
final double checkBelow = 4.0 * (fixedXBelowPlane.D * fixedXBelowPlane.D * planetModel.inverseXYScalingSquared + testPointFixedYBelowPlane.D * testPointFixedYBelowPlane.D * planetModel.inverseXYScalingSquared - 1.0);
if (checkAbove < Vector.MINIMUM_RESOLUTION_SQUARED && checkBelow < Vector.MINIMUM_RESOLUTION_SQUARED) {
//System.out.println(" Looking for intersections between travel and test point planes...");
final GeoPoint[] XIntersectionsY = travelPlaneFixedX.findIntersections(planetModel, testPointFixedYPlane);
@ -325,8 +325,8 @@ class GeoComplexPolygon extends GeoBasePolygon {
}
if (testPointFixedZAbovePlane != null && testPointFixedZBelowPlane != null && fixedXAbovePlane != null && fixedXBelowPlane != null) {
//check if planes intersects inside world
final double checkAbove = 4.0 * (fixedXAbovePlane.D * fixedXAbovePlane.D * planetModel.inverseAbSquared + testPointFixedZAbovePlane.D * testPointFixedZAbovePlane.D * planetModel.inverseCSquared - 1.0);
final double checkBelow = 4.0 * (fixedXBelowPlane.D * fixedXBelowPlane.D * planetModel.inverseAbSquared + testPointFixedZBelowPlane.D * testPointFixedZBelowPlane.D * planetModel.inverseCSquared - 1.0);
final double checkAbove = 4.0 * (fixedXAbovePlane.D * fixedXAbovePlane.D * planetModel.inverseXYScalingSquared + testPointFixedZAbovePlane.D * testPointFixedZAbovePlane.D * planetModel.inverseZScalingSquared - 1.0);
final double checkBelow = 4.0 * (fixedXBelowPlane.D * fixedXBelowPlane.D * planetModel.inverseXYScalingSquared + testPointFixedZBelowPlane.D * testPointFixedZBelowPlane.D * planetModel.inverseZScalingSquared - 1.0);
if (checkAbove < Vector.MINIMUM_RESOLUTION_SQUARED && checkBelow < Vector.MINIMUM_RESOLUTION_SQUARED) {
//System.out.println(" Looking for intersections between travel and test point planes...");
final GeoPoint[] XIntersectionsZ = travelPlaneFixedX.findIntersections(planetModel, testPointFixedZPlane);
@ -349,8 +349,8 @@ class GeoComplexPolygon extends GeoBasePolygon {
}
if (testPointFixedXAbovePlane != null && testPointFixedXBelowPlane != null && fixedYAbovePlane != null && fixedYBelowPlane != null) {
//check if planes intersects inside world
final double checkAbove = 4.0 * (testPointFixedXAbovePlane.D * testPointFixedXAbovePlane.D * planetModel.inverseAbSquared + fixedYAbovePlane.D * fixedYAbovePlane.D * planetModel.inverseAbSquared - 1.0);
final double checkBelow = 4.0 * (testPointFixedXBelowPlane.D * testPointFixedXBelowPlane.D * planetModel.inverseAbSquared + fixedYBelowPlane.D * fixedYBelowPlane.D * planetModel.inverseAbSquared - 1.0);
final double checkAbove = 4.0 * (testPointFixedXAbovePlane.D * testPointFixedXAbovePlane.D * planetModel.inverseXYScalingSquared + fixedYAbovePlane.D * fixedYAbovePlane.D * planetModel.inverseXYScalingSquared - 1.0);
final double checkBelow = 4.0 * (testPointFixedXBelowPlane.D * testPointFixedXBelowPlane.D * planetModel.inverseXYScalingSquared + fixedYBelowPlane.D * fixedYBelowPlane.D * planetModel.inverseXYScalingSquared - 1.0);
if (checkAbove < Vector.MINIMUM_RESOLUTION_SQUARED && checkBelow < Vector.MINIMUM_RESOLUTION_SQUARED) {
//System.out.println(" Looking for intersections between travel and test point planes...");
final GeoPoint[] YIntersectionsX = travelPlaneFixedY.findIntersections(planetModel, testPointFixedXPlane);
@ -373,8 +373,8 @@ class GeoComplexPolygon extends GeoBasePolygon {
}
if (testPointFixedZAbovePlane != null && testPointFixedZBelowPlane != null && fixedYAbovePlane != null && fixedYBelowPlane != null) {
//check if planes intersects inside world
final double checkAbove = 4.0 * (testPointFixedZAbovePlane.D * testPointFixedZAbovePlane.D * planetModel.inverseCSquared + fixedYAbovePlane.D * fixedYAbovePlane.D * planetModel.inverseAbSquared - 1.0);
final double checkBelow = 4.0 * (testPointFixedZBelowPlane.D * testPointFixedZBelowPlane.D * planetModel.inverseCSquared + fixedYBelowPlane.D * fixedYBelowPlane.D * planetModel.inverseAbSquared - 1.0);
final double checkAbove = 4.0 * (testPointFixedZAbovePlane.D * testPointFixedZAbovePlane.D * planetModel.inverseZScalingSquared + fixedYAbovePlane.D * fixedYAbovePlane.D * planetModel.inverseXYScalingSquared - 1.0);
final double checkBelow = 4.0 * (testPointFixedZBelowPlane.D * testPointFixedZBelowPlane.D * planetModel.inverseZScalingSquared + fixedYBelowPlane.D * fixedYBelowPlane.D * planetModel.inverseXYScalingSquared - 1.0);
if (checkAbove < Vector.MINIMUM_RESOLUTION_SQUARED && checkBelow < Vector.MINIMUM_RESOLUTION_SQUARED) {
//System.out.println(" Looking for intersections between travel and test point planes...");
final GeoPoint[] YIntersectionsZ = travelPlaneFixedY.findIntersections(planetModel, testPointFixedZPlane);
@ -397,8 +397,8 @@ class GeoComplexPolygon extends GeoBasePolygon {
}
if (testPointFixedXAbovePlane != null && testPointFixedXBelowPlane != null && fixedZAbovePlane != null && fixedZBelowPlane != null) {
//check if planes intersects inside world
final double checkAbove = 4.0 * (testPointFixedXAbovePlane.D * testPointFixedXAbovePlane.D * planetModel.inverseAbSquared + fixedZAbovePlane.D * fixedZAbovePlane.D * planetModel.inverseCSquared - 1.0);
final double checkBelow = 4.0 * (testPointFixedXBelowPlane.D * testPointFixedXBelowPlane.D * planetModel.inverseAbSquared + fixedZBelowPlane.D * fixedZBelowPlane.D * planetModel.inverseCSquared - 1.0);
final double checkAbove = 4.0 * (testPointFixedXAbovePlane.D * testPointFixedXAbovePlane.D * planetModel.inverseXYScalingSquared + fixedZAbovePlane.D * fixedZAbovePlane.D * planetModel.inverseZScalingSquared - 1.0);
final double checkBelow = 4.0 * (testPointFixedXBelowPlane.D * testPointFixedXBelowPlane.D * planetModel.inverseXYScalingSquared + fixedZBelowPlane.D * fixedZBelowPlane.D * planetModel.inverseZScalingSquared - 1.0);
if (checkAbove < Vector.MINIMUM_RESOLUTION_SQUARED && checkBelow < Vector.MINIMUM_RESOLUTION_SQUARED) {
//System.out.println(" Looking for intersections between travel and test point planes...");
final GeoPoint[] ZIntersectionsX = travelPlaneFixedZ.findIntersections(planetModel, testPointFixedXPlane);
@ -421,8 +421,8 @@ class GeoComplexPolygon extends GeoBasePolygon {
}
if (testPointFixedYAbovePlane != null && testPointFixedYBelowPlane != null && fixedZAbovePlane != null && fixedZBelowPlane != null) {
//check if planes intersects inside world
final double checkAbove = 4.0 * (testPointFixedYAbovePlane.D * testPointFixedYAbovePlane.D * planetModel.inverseAbSquared + fixedZAbovePlane.D * fixedZAbovePlane.D * planetModel.inverseCSquared - 1.0);
final double checkBelow = 4.0 * (testPointFixedYBelowPlane.D * testPointFixedYBelowPlane.D * planetModel.inverseAbSquared + fixedZBelowPlane.D * fixedZBelowPlane.D * planetModel.inverseCSquared - 1.0);
final double checkAbove = 4.0 * (testPointFixedYAbovePlane.D * testPointFixedYAbovePlane.D * planetModel.inverseXYScalingSquared + fixedZAbovePlane.D * fixedZAbovePlane.D * planetModel.inverseZScalingSquared - 1.0);
final double checkBelow = 4.0 * (testPointFixedYBelowPlane.D * testPointFixedYBelowPlane.D * planetModel.inverseXYScalingSquared + fixedZBelowPlane.D * fixedZBelowPlane.D * planetModel.inverseZScalingSquared - 1.0);
if (checkAbove < Vector.MINIMUM_RESOLUTION_SQUARED && checkBelow < Vector.MINIMUM_RESOLUTION_SQUARED) {
//System.out.println(" Looking for intersections between travel and test point planes...");
final GeoPoint[] ZIntersectionsY = travelPlaneFixedZ.findIntersections(planetModel, testPointFixedYPlane);
@ -1466,7 +1466,7 @@ class GeoComplexPolygon extends GeoBasePolygon {
final GeoPoint[] belowAbove = travelBelowPlane.findIntersections(planetModel, testPointAbovePlane, intersectionBound1, intersectionBound2);
assert belowAbove != null : "Below + above should not be coplanar";
assert ((aboveAbove.length > 0)?1:0) + ((aboveBelow.length > 0)?1:0) + ((belowBelow.length > 0)?1:0) + ((belowAbove.length > 0)?1:0) == 1 : "Can be exactly one inside point, instead was: aa="+aboveAbove.length+" ab=" + aboveBelow.length+" bb="+ belowBelow.length+" ba=" + belowAbove.length;
assert ((aboveAbove.length > 0)?1:0) + ((aboveBelow.length > 0)?1:0) + ((belowBelow.length > 0)?1:0) + ((belowAbove.length > 0)?1:0) == 1 : "Can be exactly one inside point, instead was: aa="+aboveAbove.length+" xyScaling=" + aboveBelow.length+" bb="+ belowBelow.length+" ba=" + belowAbove.length;
final GeoPoint[] insideInsidePoints;
if (aboveAbove.length > 0) {

2
lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/GeoExactCircle.java
View File

@ -81,7 +81,7 @@ class GeoExactCircle extends GeoBaseCircle {
final GeoPoint westPoint = planetModel.surfacePointOnBearing(center, radius, Math.PI * 1.5);
final GeoPoint edgePoint;
if (planetModel.c > planetModel.ab) {
if (planetModel.zScaling > planetModel.xyScaling) {
// z can be greater than x or y, so ellipse is longer in height than width
slices.add(new ApproximationSlice(center, eastPoint, Math.PI * 0.5, westPoint, Math.PI * -0.5, northPoint, 0.0, true));
slices.add(new ApproximationSlice(center, westPoint, Math.PI * 1.5, eastPoint, Math.PI * 0.5, southPoint, Math.PI, true));

2
lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/GeoWorld.java
View File

@ -36,7 +36,7 @@ class GeoWorld extends GeoBaseBBox {
*/
public GeoWorld(final PlanetModel planetModel) {
super(planetModel);
originPoint = new GeoPoint(planetModel.ab, 1.0, 0.0, 0.0);
originPoint = new GeoPoint(planetModel.xyScaling, 1.0, 0.0, 0.0);
}
/** Constructor.

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

@ -722,7 +722,7 @@ public class Plane extends Vector {
// ... can be solved by Cramer's rule:
// x = det(S0 R0 / S1 R1) / det(Q0 R0 / Q1 R1)
// y = det(Q0 S0 / Q1 S1) / det(Q0 R0 / Q1 R1)
// ... where det( a b / c d ) = ad - bc, so:
// ... where det( a b / zScaling d ) = ad - bc, so:
// x = (S0 * R1 - R0 * S1) / (Q0 * R1 - R0 * Q1)
// y = (Q0 * S1 - S0 * Q1) / (Q0 * R1 - R0 * Q1)
double x0;
@ -766,16 +766,16 @@ public class Plane extends Vector {
// Once an intersecting line is determined, the next step is to intersect that line with the ellipsoid, which
// will yield zero, one, or two points.
// The ellipsoid equation: 1,0 = x^2/a^2 + y^2/b^2 + z^2/c^2
// 1.0 = (At+A0)^2/a^2 + (Bt+B0)^2/b^2 + (Ct+C0)^2/c^2
// A^2 t^2 / a^2 + 2AA0t / a^2 + A0^2 / a^2 + B^2 t^2 / b^2 + 2BB0t / b^2 + B0^2 / b^2 + C^2 t^2 / c^2 + 2CC0t / c^2 + C0^2 / c^2 - 1,0 = 0.0
// [A^2 / a^2 + B^2 / b^2 + C^2 / c^2] t^2 + [2AA0 / a^2 + 2BB0 / b^2 + 2CC0 / c^2] t + [A0^2 / a^2 + B0^2 / b^2 + C0^2 / c^2 - 1,0] = 0.0
// The ellipsoid equation: 1,0 = x^2/a^2 + y^2/b^2 + z^2/zScaling^2
// 1.0 = (At+A0)^2/a^2 + (Bt+B0)^2/b^2 + (Ct+C0)^2/zScaling^2
// A^2 t^2 / a^2 + 2AA0t / a^2 + A0^2 / a^2 + B^2 t^2 / b^2 + 2BB0t / b^2 + B0^2 / b^2 + C^2 t^2 / zScaling^2 + 2CC0t / zScaling^2 + C0^2 / zScaling^2 - 1,0 = 0.0
// [A^2 / a^2 + B^2 / b^2 + C^2 / zScaling^2] t^2 + [2AA0 / a^2 + 2BB0 / b^2 + 2CC0 / zScaling^2] t + [A0^2 / a^2 + B0^2 / b^2 + C0^2 / zScaling^2 - 1,0] = 0.0
// Use the quadratic formula to determine t values and candidate point(s)
final double A = lineVectorX * lineVectorX * planetModel.inverseAbSquared +
lineVectorY * lineVectorY * planetModel.inverseAbSquared +
lineVectorZ * lineVectorZ * planetModel.inverseCSquared;
final double B = 2.0 * (lineVectorX * x0 * planetModel.inverseAbSquared + lineVectorY * y0 * planetModel.inverseAbSquared + lineVectorZ * z0 * planetModel.inverseCSquared);
final double C = x0 * x0 * planetModel.inverseAbSquared + y0 * y0 * planetModel.inverseAbSquared + z0 * z0 * planetModel.inverseCSquared - 1.0;
final double A = lineVectorX * lineVectorX * planetModel.inverseXYScalingSquared +
lineVectorY * lineVectorY * planetModel.inverseXYScalingSquared +
lineVectorZ * lineVectorZ * planetModel.inverseZScalingSquared;
final double B = 2.0 * (lineVectorX * x0 * planetModel.inverseXYScalingSquared + lineVectorY * y0 * planetModel.inverseXYScalingSquared + lineVectorZ * z0 * planetModel.inverseZScalingSquared);
final double C = x0 * x0 * planetModel.inverseXYScalingSquared + y0 * y0 * planetModel.inverseXYScalingSquared + z0 * z0 * planetModel.inverseZScalingSquared - 1.0;
final double BsquaredMinus = B * B - 4.0 * A * C;
if (Math.abs(BsquaredMinus) < MINIMUM_RESOLUTION_SQUARED) {
@ -898,7 +898,7 @@ public class Plane extends Vector {
// ... can be solved by Cramer's rule:
// x = det(S0 R0 / S1 R1) / det(Q0 R0 / Q1 R1)
// y = det(Q0 S0 / Q1 S1) / det(Q0 R0 / Q1 R1)
// ... where det( a b / c d ) = ad - bc, so:
// ... where det( a b / zScaling d ) = ad - bc, so:
// x = (S0 * R1 - R0 * S1) / (Q0 * R1 - R0 * Q1)
// y = (Q0 * S1 - S0 * Q1) / (Q0 * R1 - R0 * Q1)
double x0;
@ -939,16 +939,16 @@ public class Plane extends Vector {
// Once an intersecting line is determined, the next step is to intersect that line with the ellipsoid, which
// will yield zero, one, or two points.
// The ellipsoid equation: 1,0 = x^2/a^2 + y^2/b^2 + z^2/c^2
// 1.0 = (At+A0)^2/a^2 + (Bt+B0)^2/b^2 + (Ct+C0)^2/c^2
// A^2 t^2 / a^2 + 2AA0t / a^2 + A0^2 / a^2 + B^2 t^2 / b^2 + 2BB0t / b^2 + B0^2 / b^2 + C^2 t^2 / c^2 + 2CC0t / c^2 + C0^2 / c^2 - 1,0 = 0.0
// [A^2 / a^2 + B^2 / b^2 + C^2 / c^2] t^2 + [2AA0 / a^2 + 2BB0 / b^2 + 2CC0 / c^2] t + [A0^2 / a^2 + B0^2 / b^2 + C0^2 / c^2 - 1,0] = 0.0
// The ellipsoid equation: 1,0 = x^2/a^2 + y^2/b^2 + z^2/zScaling^2
// 1.0 = (At+A0)^2/a^2 + (Bt+B0)^2/b^2 + (Ct+C0)^2/zScaling^2
// A^2 t^2 / a^2 + 2AA0t / a^2 + A0^2 / a^2 + B^2 t^2 / b^2 + 2BB0t / b^2 + B0^2 / b^2 + C^2 t^2 / zScaling^2 + 2CC0t / zScaling^2 + C0^2 / zScaling^2 - 1,0 = 0.0
// [A^2 / a^2 + B^2 / b^2 + C^2 / zScaling^2] t^2 + [2AA0 / a^2 + 2BB0 / b^2 + 2CC0 / zScaling^2] t + [A0^2 / a^2 + B0^2 / b^2 + C0^2 / zScaling^2 - 1,0] = 0.0
// Use the quadratic formula to determine t values and candidate point(s)
final double A = lineVectorX * lineVectorX * planetModel.inverseAbSquared +
lineVectorY * lineVectorY * planetModel.inverseAbSquared +
lineVectorZ * lineVectorZ * planetModel.inverseCSquared;
final double B = 2.0 * (lineVectorX * x0 * planetModel.inverseAbSquared + lineVectorY * y0 * planetModel.inverseAbSquared + lineVectorZ * z0 * planetModel.inverseCSquared);
final double C = x0 * x0 * planetModel.inverseAbSquared + y0 * y0 * planetModel.inverseAbSquared + z0 * z0 * planetModel.inverseCSquared - 1.0;
final double A = lineVectorX * lineVectorX * planetModel.inverseXYScalingSquared +
lineVectorY * lineVectorY * planetModel.inverseXYScalingSquared +
lineVectorZ * lineVectorZ * planetModel.inverseZScalingSquared;
final double B = 2.0 * (lineVectorX * x0 * planetModel.inverseXYScalingSquared + lineVectorY * y0 * planetModel.inverseXYScalingSquared + lineVectorZ * z0 * planetModel.inverseZScalingSquared);
final double C = x0 * x0 * planetModel.inverseXYScalingSquared + y0 * y0 * planetModel.inverseXYScalingSquared + z0 * z0 * planetModel.inverseZScalingSquared - 1.0;
final double BsquaredMinus = B * B - 4.0 * A * C;
if (Math.abs(BsquaredMinus) < MINIMUM_RESOLUTION_SQUARED) {
@ -1048,7 +1048,7 @@ public class Plane extends Vector {
// ... can be solved by Cramer's rule:
// x = det(S0 R0 / S1 R1) / det(Q0 R0 / Q1 R1)
// y = det(Q0 S0 / Q1 S1) / det(Q0 R0 / Q1 R1)
// ... where det( a b / c d ) = ad - bc, so:
// ... where det( a b / zScaling d ) = ad - bc, so:
// x = (S0 * R1 - R0 * S1) / (Q0 * R1 - R0 * Q1)
// y = (Q0 * S1 - S0 * Q1) / (Q0 * R1 - R0 * Q1)
// We try to maximize the determinant in the denominator
@ -1138,16 +1138,16 @@ public class Plane extends Vector {
final Membership... bounds) {
// Once an intersecting line is determined, the next step is to intersect that line with the ellipsoid, which
// will yield zero, one, or two points.
// The ellipsoid equation: 1,0 = x^2/a^2 + y^2/b^2 + z^2/c^2
// 1.0 = (At+A0)^2/a^2 + (Bt+B0)^2/b^2 + (Ct+C0)^2/c^2
// A^2 t^2 / a^2 + 2AA0t / a^2 + A0^2 / a^2 + B^2 t^2 / b^2 + 2BB0t / b^2 + B0^2 / b^2 + C^2 t^2 / c^2 + 2CC0t / c^2 + C0^2 / c^2 - 1,0 = 0.0
// [A^2 / a^2 + B^2 / b^2 + C^2 / c^2] t^2 + [2AA0 / a^2 + 2BB0 / b^2 + 2CC0 / c^2] t + [A0^2 / a^2 + B0^2 / b^2 + C0^2 / c^2 - 1,0] = 0.0
// The ellipsoid equation: 1,0 = x^2/a^2 + y^2/b^2 + z^2/zScaling^2
// 1.0 = (At+A0)^2/a^2 + (Bt+B0)^2/b^2 + (Ct+C0)^2/zScaling^2
// A^2 t^2 / a^2 + 2AA0t / a^2 + A0^2 / a^2 + B^2 t^2 / b^2 + 2BB0t / b^2 + B0^2 / b^2 + C^2 t^2 / zScaling^2 + 2CC0t / zScaling^2 + C0^2 / zScaling^2 - 1,0 = 0.0
// [A^2 / a^2 + B^2 / b^2 + C^2 / zScaling^2] t^2 + [2AA0 / a^2 + 2BB0 / b^2 + 2CC0 / zScaling^2] t + [A0^2 / a^2 + B0^2 / b^2 + C0^2 / zScaling^2 - 1,0] = 0.0
// Use the quadratic formula to determine t values and candidate point(s)
final double A = lineVectorX * lineVectorX * planetModel.inverseAbSquared +
lineVectorY * lineVectorY * planetModel.inverseAbSquared +
lineVectorZ * lineVectorZ * planetModel.inverseCSquared;
final double B = 2.0 * (lineVectorX * x0 * planetModel.inverseAbSquared + lineVectorY * y0 * planetModel.inverseAbSquared + lineVectorZ * z0 * planetModel.inverseCSquared);
final double C = x0 * x0 * planetModel.inverseAbSquared + y0 * y0 * planetModel.inverseAbSquared + z0 * z0 * planetModel.inverseCSquared - 1.0;
final double A = lineVectorX * lineVectorX * planetModel.inverseXYScalingSquared +
lineVectorY * lineVectorY * planetModel.inverseXYScalingSquared +
lineVectorZ * lineVectorZ * planetModel.inverseZScalingSquared;
final double B = 2.0 * (lineVectorX * x0 * planetModel.inverseXYScalingSquared + lineVectorY * y0 * planetModel.inverseXYScalingSquared + lineVectorZ * z0 * planetModel.inverseZScalingSquared);
final double C = x0 * x0 * planetModel.inverseXYScalingSquared + y0 * y0 * planetModel.inverseXYScalingSquared + z0 * z0 * planetModel.inverseZScalingSquared - 1.0;
final double BsquaredMinus = B * B - 4.0 * A * C;
if (Math.abs(BsquaredMinus) < MINIMUM_RESOLUTION_SQUARED) {
@ -1213,7 +1213,7 @@ public class Plane extends Vector {
throw new RuntimeException("Intersection point not on original plane; point="+point+", plane="+this);
if (!q.evaluateIsZero(point))
throw new RuntimeException("Intersection point not on intersected plane; point="+point+", plane="+q);
if (Math.abs(point.x * point.x * planetModel.inverseASquared + point.y * point.y * planetModel.inverseBSquared + point.z * point.z * planetModel.inverseCSquared - 1.0) >= MINIMUM_RESOLUTION)
if (Math.abs(point.x * point.x * planetModel.inverseASquared + point.y * point.y * planetModel.inverseBSquared + point.z * point.z * planetModel.inverseZScalingSquared - 1.0) >= MINIMUM_RESOLUTION)
throw new RuntimeException("Intersection point not on ellipsoid; point="+point);
}
*/
@ -1283,9 +1283,9 @@ public class Plane extends Vector {
}
// First, compute common subexpressions
final double k = 1.0 / ((x*x + y*y)*planetModel.ab*planetModel.ab + z*z*planetModel.c*planetModel.c);
final double abSquared = planetModel.ab * planetModel.ab;
final double cSquared = planetModel.c * planetModel.c;
final double k = 1.0 / ((x*x + y*y)*planetModel.xyScaling *planetModel.xyScaling + z*z*planetModel.zScaling *planetModel.zScaling);
final double abSquared = planetModel.xyScaling * planetModel.xyScaling;
final double cSquared = planetModel.zScaling * planetModel.zScaling;
final double ASquared = A * A;
final double BSquared = B * B;
final double CSquared = C * C;
@ -1298,11 +1298,11 @@ public class Plane extends Vector {
//
// For this, we need grad(F(x,y,z)) = (dF/dx, dF/dy, dF/dz).
//
// Minimize and maximize f(x,y,z) = x, with respect to g(x,y,z) = Ax + By + Cz - D and h(x,y,z) = x^2/ab^2 + y^2/ab^2 + z^2/c^2 - 1
// Minimize and maximize f(x,y,z) = x, with respect to g(x,y,z) = Ax + By + Cz - D and h(x,y,z) = x^2/xyScaling^2 + y^2/xyScaling^2 + z^2/zScaling^2 - 1
//
// grad(f(x,y,z)) = (1,0,0)
// grad(g(x,y,z)) = (A,B,C)
// grad(h(x,y,z)) = (2x/ab^2,2y/ab^2,2z/c^2)
// grad(h(x,y,z)) = (2x/xyScaling^2,2y/xyScaling^2,2z/zScaling^2)
//
// Equations we need to simultaneously solve:
//
@ -1311,49 +1311,49 @@ public class Plane extends Vector {
// h(x,y,z) = 0
//
// Equations:
// 1 = l*A + m*2x/ab^2
// 0 = l*B + m*2y/ab^2
// 0 = l*C + m*2z/c^2
// 1 = l*A + m*2x/xyScaling^2
// 0 = l*B + m*2y/xyScaling^2
// 0 = l*C + m*2z/zScaling^2
// Ax + By + Cz + D = 0
// x^2/ab^2 + y^2/ab^2 + z^2/c^2 - 1 = 0
// x^2/xyScaling^2 + y^2/xyScaling^2 + z^2/zScaling^2 - 1 = 0
//
// Solve for x,y,z in terms of (l, m):
//
// x = ((1 - l*A) * ab^2 ) / (2 * m)
// y = (-l*B * ab^2) / ( 2 * m)
// z = (-l*C * c^2)/ (2 * m)
// x = ((1 - l*A) * xyScaling^2 ) / (2 * m)
// y = (-l*B * xyScaling^2) / ( 2 * m)
// z = (-l*C * zScaling^2)/ (2 * m)
//
// Two equations, two unknowns:
//
// A * (((1 - l*A) * ab^2 ) / (2 * m)) + B * ((-l*B * ab^2) / ( 2 * m)) + C * ((-l*C * c^2)/ (2 * m)) + D = 0
// A * (((1 - l*A) * xyScaling^2 ) / (2 * m)) + B * ((-l*B * xyScaling^2) / ( 2 * m)) + C * ((-l*C * zScaling^2)/ (2 * m)) + D = 0
//
// and
//
// (((1 - l*A) * ab^2 ) / (2 * m))^2/ab^2 + ((-l*B * ab^2) / ( 2 * m))^2/ab^2 + ((-l*C * c^2)/ (2 * m))^2/c^2 - 1 = 0
// (((1 - l*A) * xyScaling^2 ) / (2 * m))^2/xyScaling^2 + ((-l*B * xyScaling^2) / ( 2 * m))^2/xyScaling^2 + ((-l*C * zScaling^2)/ (2 * m))^2/zScaling^2 - 1 = 0
//
// Simple: solve for l and m, then find x from it.
//
// (a) Use first equation to find l in terms of m.
//
// A * (((1 - l*A) * ab^2 ) / (2 * m)) + B * ((-l*B * ab^2) / ( 2 * m)) + C * ((-l*C * c^2)/ (2 * m)) + D = 0
// A * ((1 - l*A) * ab^2 ) + B * (-l*B * ab^2) + C * (-l*C * c^2) + D * 2 * m = 0
// A * ab^2 - l*A^2* ab^2 - B^2 * l * ab^2 - C^2 * l * c^2 + D * 2 * m = 0
// - l *(A^2* ab^2 + B^2 * ab^2 + C^2 * c^2) + (A * ab^2 + D * 2 * m) = 0
// l = (A * ab^2 + D * 2 * m) / (A^2* ab^2 + B^2 * ab^2 + C^2 * c^2)
// l = A * ab^2 / (A^2* ab^2 + B^2 * ab^2 + C^2 * c^2) + m * 2 * D / (A^2* ab^2 + B^2 * ab^2 + C^2 * c^2)
// A * (((1 - l*A) * xyScaling^2 ) / (2 * m)) + B * ((-l*B * xyScaling^2) / ( 2 * m)) + C * ((-l*C * zScaling^2)/ (2 * m)) + D = 0
// A * ((1 - l*A) * xyScaling^2 ) + B * (-l*B * xyScaling^2) + C * (-l*C * zScaling^2) + D * 2 * m = 0
// A * xyScaling^2 - l*A^2* xyScaling^2 - B^2 * l * xyScaling^2 - C^2 * l * zScaling^2 + D * 2 * m = 0
// - l *(A^2* xyScaling^2 + B^2 * xyScaling^2 + C^2 * zScaling^2) + (A * xyScaling^2 + D * 2 * m) = 0
// l = (A * xyScaling^2 + D * 2 * m) / (A^2* xyScaling^2 + B^2 * xyScaling^2 + C^2 * zScaling^2)
// l = A * xyScaling^2 / (A^2* xyScaling^2 + B^2 * xyScaling^2 + C^2 * zScaling^2) + m * 2 * D / (A^2* xyScaling^2 + B^2 * xyScaling^2 + C^2 * zScaling^2)
//
// For convenience:
//
// k = 1.0 / (A^2* ab^2 + B^2 * ab^2 + C^2 * c^2)
// k = 1.0 / (A^2* xyScaling^2 + B^2 * xyScaling^2 + C^2 * zScaling^2)
//
// Then:
//
// l = A * ab^2 * k + m * 2 * D * k
// l = k * (A*ab^2 + m*2*D)
// l = A * xyScaling^2 * k + m * 2 * D * k
// l = k * (A*xyScaling^2 + m*2*D)
//
// For further convenience:
//
// q = A*ab^2*k
// q = A*xyScaling^2*k
// r = 2*D*k
//
// l = (r*m + q)
@ -1361,23 +1361,23 @@ public class Plane extends Vector {
//
// (b) Simplify the second equation before substitution
//
// (((1 - l*A) * ab^2 ) / (2 * m))^2/ab^2 + ((-l*B * ab^2) / ( 2 * m))^2/ab^2 + ((-l*C * c^2)/ (2 * m))^2/c^2 - 1 = 0
// ((1 - l*A) * ab^2 )^2/ab^2 + (-l*B * ab^2)^2/ab^2 + (-l*C * c^2)^2/c^2 = 4 * m^2
// (1 - l*A)^2 * ab^2 + (-l*B)^2 * ab^2 + (-l*C)^2 * c^2 = 4 * m^2
// (1 - 2*l*A + l^2*A^2) * ab^2 + l^2*B^2 * ab^2 + l^2*C^2 * c^2 = 4 * m^2
// ab^2 - 2*A*ab^2*l + A^2*ab^2*l^2 + B^2*ab^2*l^2 + C^2*c^2*l^2 - 4*m^2 = 0
// (((1 - l*A) * xyScaling^2 ) / (2 * m))^2/xyScaling^2 + ((-l*B * xyScaling^2) / ( 2 * m))^2/xyScaling^2 + ((-l*C * zScaling^2)/ (2 * m))^2/zScaling^2 - 1 = 0
// ((1 - l*A) * xyScaling^2 )^2/xyScaling^2 + (-l*B * xyScaling^2)^2/xyScaling^2 + (-l*C * zScaling^2)^2/zScaling^2 = 4 * m^2
// (1 - l*A)^2 * xyScaling^2 + (-l*B)^2 * xyScaling^2 + (-l*C)^2 * zScaling^2 = 4 * m^2
// (1 - 2*l*A + l^2*A^2) * xyScaling^2 + l^2*B^2 * xyScaling^2 + l^2*C^2 * zScaling^2 = 4 * m^2
// xyScaling^2 - 2*A*xyScaling^2*l + A^2*xyScaling^2*l^2 + B^2*xyScaling^2*l^2 + C^2*zScaling^2*l^2 - 4*m^2 = 0
//
// (c) Substitute for l, l^2
// (zScaling) Substitute for l, l^2
//
// ab^2 - 2*A*ab^2*(r*m + q) + A^2*ab^2*(r^2 * m^2 + 2*r*m*q + q^2) + B^2*ab^2*(r^2 * m^2 + 2*r*m*q + q^2) + C^2*c^2*(r^2 * m^2 + 2*r*m*q + q^2) - 4*m^2 = 0
// ab^2 - 2*A*ab^2*r*m - 2*A*ab^2*q + A^2*ab^2*r^2*m^2 + 2*A^2*ab^2*r*q*m +
// A^2*ab^2*q^2 + B^2*ab^2*r^2*m^2 + 2*B^2*ab^2*r*q*m + B^2*ab^2*q^2 + C^2*c^2*r^2*m^2 + 2*C^2*c^2*r*q*m + C^2*c^2*q^2 - 4*m^2 = 0
// xyScaling^2 - 2*A*xyScaling^2*(r*m + q) + A^2*xyScaling^2*(r^2 * m^2 + 2*r*m*q + q^2) + B^2*xyScaling^2*(r^2 * m^2 + 2*r*m*q + q^2) + C^2*zScaling^2*(r^2 * m^2 + 2*r*m*q + q^2) - 4*m^2 = 0
// xyScaling^2 - 2*A*xyScaling^2*r*m - 2*A*xyScaling^2*q + A^2*xyScaling^2*r^2*m^2 + 2*A^2*xyScaling^2*r*q*m +
// A^2*xyScaling^2*q^2 + B^2*xyScaling^2*r^2*m^2 + 2*B^2*xyScaling^2*r*q*m + B^2*xyScaling^2*q^2 + C^2*zScaling^2*r^2*m^2 + 2*C^2*zScaling^2*r*q*m + C^2*zScaling^2*q^2 - 4*m^2 = 0
//
// (d) Group
//
// m^2 * [A^2*ab^2*r^2 + B^2*ab^2*r^2 + C^2*c^2*r^2 - 4] +
// m * [- 2*A*ab^2*r + 2*A^2*ab^2*r*q + 2*B^2*ab^2*r*q + 2*C^2*c^2*r*q] +
// [ab^2 - 2*A*ab^2*q + A^2*ab^2*q^2 + B^2*ab^2*q^2 + C^2*c^2*q^2] = 0
// m^2 * [A^2*xyScaling^2*r^2 + B^2*xyScaling^2*r^2 + C^2*zScaling^2*r^2 - 4] +
// m * [- 2*A*xyScaling^2*r + 2*A^2*xyScaling^2*r*q + 2*B^2*xyScaling^2*r*q + 2*C^2*zScaling^2*r*q] +
// [xyScaling^2 - 2*A*xyScaling^2*q + A^2*xyScaling^2*q^2 + B^2*xyScaling^2*q^2 + C^2*zScaling^2*q^2] = 0
// Useful subexpressions for this bound
final double q = A*abSquared*k;
@ -1396,14 +1396,14 @@ public class Plane extends Vector {
// Valid?
if (Math.abs(m) >= MINIMUM_RESOLUTION) {
final double l = r * m + q;
// x = ((1 - l*A) * ab^2 ) / (2 * m)
// y = (-l*B * ab^2) / ( 2 * m)
// z = (-l*C * c^2)/ (2 * m)
// x = ((1 - l*A) * xyScaling^2 ) / (2 * m)
// y = (-l*B * xyScaling^2) / ( 2 * m)
// z = (-l*C * zScaling^2)/ (2 * m)
final double denom0 = 0.5 / m;
final GeoPoint thePoint = new GeoPoint((1.0-l*A) * abSquared * denom0, -l*B * abSquared * denom0, -l*C * cSquared * denom0);
//Math is not quite accurate enough for this
//assert planetModel.pointOnSurface(thePoint): "Point: "+thePoint+"; Planetmodel="+planetModel+"; A="+A+" B="+B+" C="+C+" D="+D+" planetfcn="+
// (thePoint.x*thePoint.x*planetModel.inverseAb*planetModel.inverseAb + thePoint.y*thePoint.y*planetModel.inverseAb*planetModel.inverseAb + thePoint.z*thePoint.z*planetModel.inverseC*planetModel.inverseC);
// (thePoint.x*thePoint.x*planetModel.inverseXYScaling*planetModel.inverseXYScaling + thePoint.y*thePoint.y*planetModel.inverseXYScaling*planetModel.inverseXYScaling + thePoint.z*thePoint.z*planetModel.inverseZScaling*planetModel.inverseZScaling);
//assert evaluateIsZero(thePoint): "Evaluation of point: "+evaluate(thePoint);
addPoint(boundsInfo, bounds, thePoint);
} else {
@ -1421,18 +1421,18 @@ public class Plane extends Vector {
if (Math.abs(m1) >= MINIMUM_RESOLUTION || Math.abs(m2) >= MINIMUM_RESOLUTION) {
final double l1 = r * m1 + q;
final double l2 = r * m2 + q;
// x = ((1 - l*A) * ab^2 ) / (2 * m)
// y = (-l*B * ab^2) / ( 2 * m)
// z = (-l*C * c^2)/ (2 * m)
// x = ((1 - l*A) * xyScaling^2 ) / (2 * m)
// y = (-l*B * xyScaling^2) / ( 2 * m)
// z = (-l*C * zScaling^2)/ (2 * m)
final double denom1 = 0.5 / m1;
final double denom2 = 0.5 / m2;
final GeoPoint thePoint1 = new GeoPoint((1.0-l1*A) * abSquared * denom1, -l1*B * abSquared * denom1, -l1*C * cSquared * denom1);
final GeoPoint thePoint2 = new GeoPoint((1.0-l2*A) * abSquared * denom2, -l2*B * abSquared * denom2, -l2*C * cSquared * denom2);
//Math is not quite accurate enough for this
//assert planetModel.pointOnSurface(thePoint1): "Point1: "+thePoint1+"; Planetmodel="+planetModel+"; A="+A+" B="+B+" C="+C+" D="+D+" planetfcn="+
// (thePoint1.x*thePoint1.x*planetModel.inverseAb*planetModel.inverseAb + thePoint1.y*thePoint1.y*planetModel.inverseAb*planetModel.inverseAb + thePoint1.z*thePoint1.z*planetModel.inverseC*planetModel.inverseC);
// (thePoint1.x*thePoint1.x*planetModel.inverseXYScaling*planetModel.inverseXYScaling + thePoint1.y*thePoint1.y*planetModel.inverseXYScaling*planetModel.inverseXYScaling + thePoint1.z*thePoint1.z*planetModel.inverseZScaling*planetModel.inverseZScaling);
//assert planetModel.pointOnSurface(thePoint2): "Point1: "+thePoint2+"; Planetmodel="+planetModel+"; A="+A+" B="+B+" C="+C+" D="+D+" planetfcn="+
// (thePoint2.x*thePoint2.x*planetModel.inverseAb*planetModel.inverseAb + thePoint2.y*thePoint2.y*planetModel.inverseAb*planetModel.inverseAb + thePoint2.z*thePoint2.z*planetModel.inverseC*planetModel.inverseC);
// (thePoint2.x*thePoint2.x*planetModel.inverseXYScaling*planetModel.inverseXYScaling + thePoint2.y*thePoint2.y*planetModel.inverseXYScaling*planetModel.inverseXYScaling + thePoint2.z*thePoint2.z*planetModel.inverseZScaling*planetModel.inverseZScaling);
//assert evaluateIsZero(thePoint1): "Evaluation of point1: "+evaluate(thePoint1);
//assert evaluateIsZero(thePoint2): "Evaluation of point2: "+evaluate(thePoint2);
addPoint(boundsInfo, bounds, thePoint1);
@ -1445,17 +1445,17 @@ public class Plane extends Vector {
// No solutions
}
} else if (Math.abs(b) > MINIMUM_RESOLUTION_SQUARED) {
// a = 0, so m = - c / b
// a = 0, so m = - zScaling / b
final double m = -c / b;
final double l = r * m + q;
// x = ((1 - l*A) * ab^2 ) / (2 * m)
// y = (-l*B * ab^2) / ( 2 * m)
// z = (-l*C * c^2)/ (2 * m)
// x = ((1 - l*A) * xyScaling^2 ) / (2 * m)
// y = (-l*B * xyScaling^2) / ( 2 * m)
// z = (-l*C * zScaling^2)/ (2 * m)
final double denom0 = 0.5 / m;
final GeoPoint thePoint = new GeoPoint((1.0-l*A) * abSquared * denom0, -l*B * abSquared * denom0, -l*C * cSquared * denom0);
//Math is not quite accurate enough for this
//assert planetModel.pointOnSurface(thePoint): "Point: "+thePoint+"; Planetmodel="+planetModel+"; A="+A+" B="+B+" C="+C+" D="+D+" planetfcn="+
// (thePoint.x*thePoint.x*planetModel.inverseAb*planetModel.inverseAb + thePoint.y*thePoint.y*planetModel.inverseAb*planetModel.inverseAb + thePoint.z*thePoint.z*planetModel.inverseC*planetModel.inverseC);
// (thePoint.x*thePoint.x*planetModel.inverseXYScaling*planetModel.inverseXYScaling + thePoint.y*thePoint.y*planetModel.inverseXYScaling*planetModel.inverseXYScaling + thePoint.z*thePoint.z*planetModel.inverseZScaling*planetModel.inverseZScaling);
//assert evaluateIsZero(thePoint): "Evaluation of point: "+evaluate(thePoint);
addPoint(boundsInfo, bounds, thePoint);
} else {
@ -1469,11 +1469,11 @@ public class Plane extends Vector {
//
// For this, we need grad(F(x,y,z)) = (dF/dx, dF/dy, dF/dz).
//
// Minimize and maximize f(x,y,z) = y, with respect to g(x,y,z) = Ax + By + Cz - D and h(x,y,z) = x^2/ab^2 + y^2/ab^2 + z^2/c^2 - 1
// Minimize and maximize f(x,y,z) = y, with respect to g(x,y,z) = Ax + By + Cz - D and h(x,y,z) = x^2/xyScaling^2 + y^2/xyScaling^2 + z^2/zScaling^2 - 1
//
// grad(f(x,y,z)) = (0,1,0)
// grad(g(x,y,z)) = (A,B,C)
// grad(h(x,y,z)) = (2x/ab^2,2y/ab^2,2z/c^2)
// grad(h(x,y,z)) = (2x/xyScaling^2,2y/xyScaling^2,2z/zScaling^2)
//
// Equations we need to simultaneously solve:
//
@ -1482,49 +1482,49 @@ public class Plane extends Vector {
// h(x,y,z) = 0
//
// Equations:
// 0 = l*A + m*2x/ab^2
// 1 = l*B + m*2y/ab^2
// 0 = l*C + m*2z/c^2
// 0 = l*A + m*2x/xyScaling^2
// 1 = l*B + m*2y/xyScaling^2
// 0 = l*C + m*2z/zScaling^2
// Ax + By + Cz + D = 0
// x^2/ab^2 + y^2/ab^2 + z^2/c^2 - 1 = 0
// x^2/xyScaling^2 + y^2/xyScaling^2 + z^2/zScaling^2 - 1 = 0
//
// Solve for x,y,z in terms of (l, m):
//
// x = (-l*A * ab^2 ) / (2 * m)
// y = ((1 - l*B) * ab^2) / ( 2 * m)
// z = (-l*C * c^2)/ (2 * m)
// x = (-l*A * xyScaling^2 ) / (2 * m)
// y = ((1 - l*B) * xyScaling^2) / ( 2 * m)
// z = (-l*C * zScaling^2)/ (2 * m)
//
// Two equations, two unknowns:
//
// A * ((-l*A * ab^2 ) / (2 * m)) + B * (((1 - l*B) * ab^2) / ( 2 * m)) + C * ((-l*C * c^2)/ (2 * m)) + D = 0
// A * ((-l*A * xyScaling^2 ) / (2 * m)) + B * (((1 - l*B) * xyScaling^2) / ( 2 * m)) + C * ((-l*C * zScaling^2)/ (2 * m)) + D = 0
//
// and
//
// ((-l*A * ab^2 ) / (2 * m))^2/ab^2 + (((1 - l*B) * ab^2) / ( 2 * m))^2/ab^2 + ((-l*C * c^2)/ (2 * m))^2/c^2 - 1 = 0
// ((-l*A * xyScaling^2 ) / (2 * m))^2/xyScaling^2 + (((1 - l*B) * xyScaling^2) / ( 2 * m))^2/xyScaling^2 + ((-l*C * zScaling^2)/ (2 * m))^2/zScaling^2 - 1 = 0
//
// Simple: solve for l and m, then find y from it.
//
// (a) Use first equation to find l in terms of m.
//
// A * ((-l*A * ab^2 ) / (2 * m)) + B * (((1 - l*B) * ab^2) / ( 2 * m)) + C * ((-l*C * c^2)/ (2 * m)) + D = 0
// A * (-l*A * ab^2 ) + B * ((1-l*B) * ab^2) + C * (-l*C * c^2) + D * 2 * m = 0
// -A^2*l*ab^2 + B*ab^2 - l*B^2*ab^2 - C^2*l*c^2 + D*2*m = 0
// - l *(A^2* ab^2 + B^2 * ab^2 + C^2 * c^2) + (B * ab^2 + D * 2 * m) = 0
// l = (B * ab^2 + D * 2 * m) / (A^2* ab^2 + B^2 * ab^2 + C^2 * c^2)
// l = B * ab^2 / (A^2* ab^2 + B^2 * ab^2 + C^2 * c^2) + m * 2 * D / (A^2* ab^2 + B^2 * ab^2 + C^2 * c^2)
// A * ((-l*A * xyScaling^2 ) / (2 * m)) + B * (((1 - l*B) * xyScaling^2) / ( 2 * m)) + C * ((-l*C * zScaling^2)/ (2 * m)) + D = 0
// A * (-l*A * xyScaling^2 ) + B * ((1-l*B) * xyScaling^2) + C * (-l*C * zScaling^2) + D * 2 * m = 0
// -A^2*l*xyScaling^2 + B*xyScaling^2 - l*B^2*xyScaling^2 - C^2*l*zScaling^2 + D*2*m = 0
// - l *(A^2* xyScaling^2 + B^2 * xyScaling^2 + C^2 * zScaling^2) + (B * xyScaling^2 + D * 2 * m) = 0
// l = (B * xyScaling^2 + D * 2 * m) / (A^2* xyScaling^2 + B^2 * xyScaling^2 + C^2 * zScaling^2)
// l = B * xyScaling^2 / (A^2* xyScaling^2 + B^2 * xyScaling^2 + C^2 * zScaling^2) + m * 2 * D / (A^2* xyScaling^2 + B^2 * xyScaling^2 + C^2 * zScaling^2)
//
// For convenience:
//
// k = 1.0 / (A^2* ab^2 + B^2 * ab^2 + C^2 * c^2)
// k = 1.0 / (A^2* xyScaling^2 + B^2 * xyScaling^2 + C^2 * zScaling^2)
//
// Then:
//
// l = B * ab^2 * k + m * 2 * D * k
// l = k * (B*ab^2 + m*2*D)
// l = B * xyScaling^2 * k + m * 2 * D * k
// l = k * (B*xyScaling^2 + m*2*D)
//
// For further convenience:
//
// q = B*ab^2*k
// q = B*xyScaling^2*k
// r = 2*D*k
//
// l = (r*m + q)
@ -1532,23 +1532,23 @@ public class Plane extends Vector {
//
// (b) Simplify the second equation before substitution
//
// ((-l*A * ab^2 ) / (2 * m))^2/ab^2 + (((1 - l*B) * ab^2) / ( 2 * m))^2/ab^2 + ((-l*C * c^2)/ (2 * m))^2/c^2 - 1 = 0
// (-l*A * ab^2 )^2/ab^2 + ((1 - l*B) * ab^2)^2/ab^2 + (-l*C * c^2)^2/c^2 = 4 * m^2
// (-l*A)^2 * ab^2 + (1 - l*B)^2 * ab^2 + (-l*C)^2 * c^2 = 4 * m^2
// l^2*A^2 * ab^2 + (1 - 2*l*B + l^2*B^2) * ab^2 + l^2*C^2 * c^2 = 4 * m^2
// A^2*ab^2*l^2 + ab^2 - 2*B*ab^2*l + B^2*ab^2*l^2 + C^2*c^2*l^2 - 4*m^2 = 0
// ((-l*A * xyScaling^2 ) / (2 * m))^2/xyScaling^2 + (((1 - l*B) * xyScaling^2) / ( 2 * m))^2/xyScaling^2 + ((-l*C * zScaling^2)/ (2 * m))^2/zScaling^2 - 1 = 0
// (-l*A * xyScaling^2 )^2/xyScaling^2 + ((1 - l*B) * xyScaling^2)^2/xyScaling^2 + (-l*C * zScaling^2)^2/zScaling^2 = 4 * m^2
// (-l*A)^2 * xyScaling^2 + (1 - l*B)^2 * xyScaling^2 + (-l*C)^2 * zScaling^2 = 4 * m^2
// l^2*A^2 * xyScaling^2 + (1 - 2*l*B + l^2*B^2) * xyScaling^2 + l^2*C^2 * zScaling^2 = 4 * m^2
// A^2*xyScaling^2*l^2 + xyScaling^2 - 2*B*xyScaling^2*l + B^2*xyScaling^2*l^2 + C^2*zScaling^2*l^2 - 4*m^2 = 0
//
// (c) Substitute for l, l^2
// (zScaling) Substitute for l, l^2
//
// A^2*ab^2*(r^2 * m^2 + 2*r*m*q + q^2) + ab^2 - 2*B*ab^2*(r*m + q) + B^2*ab^2*(r^2 * m^2 + 2*r*m*q + q^2) + C^2*c^2*(r^2 * m^2 + 2*r*m*q + q^2) - 4*m^2 = 0
// A^2*ab^2*r^2*m^2 + 2*A^2*ab^2*r*q*m + A^2*ab^2*q^2 + ab^2 - 2*B*ab^2*r*m - 2*B*ab^2*q + B^2*ab^2*r^2*m^2 +
// 2*B^2*ab^2*r*q*m + B^2*ab^2*q^2 + C^2*c^2*r^2*m^2 + 2*C^2*c^2*r*q*m + C^2*c^2*q^2 - 4*m^2 = 0
// A^2*xyScaling^2*(r^2 * m^2 + 2*r*m*q + q^2) + xyScaling^2 - 2*B*xyScaling^2*(r*m + q) + B^2*xyScaling^2*(r^2 * m^2 + 2*r*m*q + q^2) + C^2*zScaling^2*(r^2 * m^2 + 2*r*m*q + q^2) - 4*m^2 = 0
// A^2*xyScaling^2*r^2*m^2 + 2*A^2*xyScaling^2*r*q*m + A^2*xyScaling^2*q^2 + xyScaling^2 - 2*B*xyScaling^2*r*m - 2*B*xyScaling^2*q + B^2*xyScaling^2*r^2*m^2 +
// 2*B^2*xyScaling^2*r*q*m + B^2*xyScaling^2*q^2 + C^2*zScaling^2*r^2*m^2 + 2*C^2*zScaling^2*r*q*m + C^2*zScaling^2*q^2 - 4*m^2 = 0
//
// (d) Group
//
// m^2 * [A^2*ab^2*r^2 + B^2*ab^2*r^2 + C^2*c^2*r^2 - 4] +
// m * [2*A^2*ab^2*r*q - 2*B*ab^2*r + 2*B^2*ab^2*r*q + 2*C^2*c^2*r*q] +
// [A^2*ab^2*q^2 + ab^2 - 2*B*ab^2*q + B^2*ab^2*q^2 + C^2*c^2*q^2] = 0
// m^2 * [A^2*xyScaling^2*r^2 + B^2*xyScaling^2*r^2 + C^2*zScaling^2*r^2 - 4] +
// m * [2*A^2*xyScaling^2*r*q - 2*B*xyScaling^2*r + 2*B^2*xyScaling^2*r*q + 2*C^2*zScaling^2*r*q] +
// [A^2*xyScaling^2*q^2 + xyScaling^2 - 2*B*xyScaling^2*q + B^2*xyScaling^2*q^2 + C^2*zScaling^2*q^2] = 0
//System.err.println(" computing Y bound");
@ -1569,14 +1569,14 @@ public class Plane extends Vector {
// Valid?
if (Math.abs(m) >= MINIMUM_RESOLUTION) {
final double l = r * m + q;
// x = (-l*A * ab^2 ) / (2 * m)
// y = ((1.0-l*B) * ab^2) / ( 2 * m)
// z = (-l*C * c^2)/ (2 * m)
// x = (-l*A * xyScaling^2 ) / (2 * m)
// y = ((1.0-l*B) * xyScaling^2) / ( 2 * m)
// z = (-l*C * zScaling^2)/ (2 * m)
final double denom0 = 0.5 / m;
final GeoPoint thePoint = new GeoPoint(-l*A * abSquared * denom0, (1.0-l*B) * abSquared * denom0, -l*C * cSquared * denom0);
//Math is not quite accurate enough for this
//assert planetModel.pointOnSurface(thePoint): "Point: "+thePoint+"; Planetmodel="+planetModel+"; A="+A+" B="+B+" C="+C+" D="+D+" planetfcn="+
// (thePoint1.x*thePoint.x*planetModel.inverseAb*planetModel.inverseAb + thePoint.y*thePoint.y*planetModel.inverseAb*planetModel.inverseAb + thePoint.z*thePoint.z*planetModel.inverseC*planetModel.inverseC);
// (thePoint1.x*thePoint.x*planetModel.inverseXYScaling*planetModel.inverseXYScaling + thePoint.y*thePoint.y*planetModel.inverseXYScaling*planetModel.inverseXYScaling + thePoint.z*thePoint.z*planetModel.inverseZScaling*planetModel.inverseZScaling);
//assert evaluateIsZero(thePoint): "Evaluation of point: "+evaluate(thePoint);
addPoint(boundsInfo, bounds, thePoint);
} else {
@ -1594,18 +1594,18 @@ public class Plane extends Vector {
if (Math.abs(m1) >= MINIMUM_RESOLUTION || Math.abs(m2) >= MINIMUM_RESOLUTION) {
final double l1 = r * m1 + q;
final double l2 = r * m2 + q;
// x = (-l*A * ab^2 ) / (2 * m)
// y = ((1.0-l*B) * ab^2) / ( 2 * m)
// z = (-l*C * c^2)/ (2 * m)
// x = (-l*A * xyScaling^2 ) / (2 * m)
// y = ((1.0-l*B) * xyScaling^2) / ( 2 * m)
// z = (-l*C * zScaling^2)/ (2 * m)
final double denom1 = 0.5 / m1;
final double denom2 = 0.5 / m2;
final GeoPoint thePoint1 = new GeoPoint(-l1*A * abSquared * denom1, (1.0-l1*B) * abSquared * denom1, -l1*C * cSquared * denom1);
final GeoPoint thePoint2 = new GeoPoint(-l2*A * abSquared * denom2, (1.0-l2*B) * abSquared * denom2, -l2*C * cSquared * denom2);
//Math is not quite accurate enough for this
//assert planetModel.pointOnSurface(thePoint1): "Point1: "+thePoint1+"; Planetmodel="+planetModel+"; A="+A+" B="+B+" C="+C+" D="+D+" planetfcn="+
// (thePoint1.x*thePoint1.x*planetModel.inverseAb*planetModel.inverseAb + thePoint1.y*thePoint1.y*planetModel.inverseAb*planetModel.inverseAb + thePoint1.z*thePoint1.z*planetModel.inverseC*planetModel.inverseC);
// (thePoint1.x*thePoint1.x*planetModel.inverseXYScaling*planetModel.inverseXYScaling + thePoint1.y*thePoint1.y*planetModel.inverseXYScaling*planetModel.inverseXYScaling + thePoint1.z*thePoint1.z*planetModel.inverseZScaling*planetModel.inverseZScaling);
//assert planetModel.pointOnSurface(thePoint2): "Point2: "+thePoint2+"; Planetmodel="+planetModel+"; A="+A+" B="+B+" C="+C+" D="+D+" planetfcn="+
// (thePoint2.x*thePoint2.x*planetModel.inverseAb*planetModel.inverseAb + thePoint2.y*thePoint2.y*planetModel.inverseAb*planetModel.inverseAb + thePoint2.z*thePoint2.z*planetModel.inverseC*planetModel.inverseC);
// (thePoint2.x*thePoint2.x*planetModel.inverseXYScaling*planetModel.inverseXYScaling + thePoint2.y*thePoint2.y*planetModel.inverseXYScaling*planetModel.inverseXYScaling + thePoint2.z*thePoint2.z*planetModel.inverseZScaling*planetModel.inverseZScaling);
//assert evaluateIsZero(thePoint1): "Evaluation of point1: "+evaluate(thePoint1);
//assert evaluateIsZero(thePoint2): "Evaluation of point2: "+evaluate(thePoint2);
addPoint(boundsInfo, bounds, thePoint1);
@ -1618,17 +1618,17 @@ public class Plane extends Vector {
// No solutions
}
} else if (Math.abs(b) > MINIMUM_RESOLUTION_SQUARED) {
// a = 0, so m = - c / b
// a = 0, so m = - zScaling / b
final double m = -c / b;
final double l = r * m + q;
// x = ( -l*A * ab^2 ) / (2 * m)
// y = ((1-l*B) * ab^2) / ( 2 * m)
// z = (-l*C * c^2)/ (2 * m)
// x = ( -l*A * xyScaling^2 ) / (2 * m)
// y = ((1-l*B) * xyScaling^2) / ( 2 * m)
// z = (-l*C * zScaling^2)/ (2 * m)
final double denom0 = 0.5 / m;
final GeoPoint thePoint = new GeoPoint(-l*A * abSquared * denom0, (1.0-l*B) * abSquared * denom0, -l*C * cSquared * denom0);
//Math is not quite accurate enough for this
//assert planetModel.pointOnSurface(thePoint): "Point: "+thePoint+"; Planetmodel="+planetModel+"; A="+A+" B="+B+" C="+C+" D="+D+" planetfcn="+
// (thePoint.x*thePoint.x*planetModel.inverseAb*planetModel.inverseAb + thePoint.y*thePoint.y*planetModel.inverseAb*planetModel.inverseAb + thePoint.z*thePoint.z*planetModel.inverseC*planetModel.inverseC);
// (thePoint.x*thePoint.x*planetModel.inverseXYScaling*planetModel.inverseXYScaling + thePoint.y*thePoint.y*planetModel.inverseXYScaling*planetModel.inverseXYScaling + thePoint.z*thePoint.z*planetModel.inverseZScaling*planetModel.inverseZScaling);
//assert evaluateIsZero(thePoint): "Evaluation of point: "+evaluate(thePoint);
addPoint(boundsInfo, bounds, thePoint);
} else {
@ -1734,9 +1734,9 @@ public class Plane extends Vector {
// Group:
// y^2 * [B^2/a^2 + A^2/b^2] + y [2BD/a^2] + [D^2/a^2-A^2] = 0
a = B * B * planetModel.inverseAbSquared + A * A * planetModel.inverseAbSquared;
b = 2.0 * B * D * planetModel.inverseAbSquared;
c = D * D * planetModel.inverseAbSquared - A * A;
a = B * B * planetModel.inverseXYScalingSquared + A * A * planetModel.inverseXYScalingSquared;
b = 2.0 * B * D * planetModel.inverseXYScalingSquared;
c = D * D * planetModel.inverseXYScalingSquared - A * A;
double sqrtClause = b * b - 4.0 * a * c;
@ -1767,9 +1767,9 @@ public class Plane extends Vector {
// Use equation suitable for B != 0
// Since I != 0, we rewrite:
// y = (-Ax - D)/B
a = B * B * planetModel.inverseAbSquared + A * A * planetModel.inverseAbSquared;
b = 2.0 * A * D * planetModel.inverseAbSquared;
c = D * D * planetModel.inverseAbSquared - B * B;
a = B * B * planetModel.inverseXYScalingSquared + A * A * planetModel.inverseXYScalingSquared;
b = 2.0 * A * D * planetModel.inverseXYScalingSquared;
c = D * D * planetModel.inverseXYScalingSquared - B * B;
double sqrtClause = b * b - 4.0 * a * c;
@ -1807,21 +1807,21 @@ public class Plane extends Vector {
// From plane:
// z = (-Ax - By - D) / C
// From ellipsoid:
// x^2/a^2 + y^2/b^2 + [(-Ax - By - D) / C]^2/c^2 = 1
// x^2/a^2 + y^2/b^2 + [(-Ax - By - D) / C]^2/zScaling^2 = 1
// Simplify/expand:
// C^2*x^2/a^2 + C^2*y^2/b^2 + (-Ax - By - D)^2/c^2 = C^2
// C^2*x^2/a^2 + C^2*y^2/b^2 + (-Ax - By - D)^2/zScaling^2 = C^2
//
// x^2 * C^2/a^2 + y^2 * C^2/b^2 + x^2 * A^2/c^2 + ABxy/c^2 + ADx/c^2 + ABxy/c^2 + y^2 * B^2/c^2 + BDy/c^2 + ADx/c^2 + BDy/c^2 + D^2/c^2 = C^2
// x^2 * C^2/a^2 + y^2 * C^2/b^2 + x^2 * A^2/zScaling^2 + ABxy/zScaling^2 + ADx/zScaling^2 + ABxy/zScaling^2 + y^2 * B^2/zScaling^2 + BDy/zScaling^2 + ADx/zScaling^2 + BDy/zScaling^2 + D^2/zScaling^2 = C^2
// Group:
// [A^2/c^2 + C^2/a^2] x^2 + [B^2/c^2 + C^2/b^2] y^2 + [2AB/c^2]xy + [2AD/c^2]x + [2BD/c^2]y + [D^2/c^2-C^2] = 0
// [A^2/zScaling^2 + C^2/a^2] x^2 + [B^2/zScaling^2 + C^2/b^2] y^2 + [2AB/zScaling^2]xy + [2AD/zScaling^2]x + [2BD/zScaling^2]y + [D^2/zScaling^2-C^2] = 0
// For convenience, introduce post-projection coefficient variables to make life easier.
// E x^2 + F y^2 + G xy + H x + I y + J = 0
double E = A * A * planetModel.inverseCSquared + C * C * planetModel.inverseAbSquared;
double F = B * B * planetModel.inverseCSquared + C * C * planetModel.inverseAbSquared;
double G = 2.0 * A * B * planetModel.inverseCSquared;
double H = 2.0 * A * D * planetModel.inverseCSquared;
double I = 2.0 * B * D * planetModel.inverseCSquared;
double J = D * D * planetModel.inverseCSquared - C * C;
double E = A * A * planetModel.inverseZScalingSquared + C * C * planetModel.inverseXYScalingSquared;
double F = B * B * planetModel.inverseZScalingSquared + C * C * planetModel.inverseXYScalingSquared;
double G = 2.0 * A * B * planetModel.inverseZScalingSquared;
double H = 2.0 * A * D * planetModel.inverseZScalingSquared;
double I = 2.0 * B * D * planetModel.inverseZScalingSquared;
double J = D * D * planetModel.inverseZScalingSquared - C * C;
//System.err.println("E = " + E + " F = " + F + " G = " + G + " H = "+ H + " I = " + I + " J = " + J);
@ -1875,7 +1875,7 @@ public class Plane extends Vector {
b = 4.0 * E * I * J - 2.0 * G * H * J;
c = 4.0 * E * J * J - J * H * H;
//System.out.println("a="+a+" b="+b+" c="+c);
//System.out.println("a="+a+" b="+b+" zScaling="+zScaling);
double sqrtClause = b * b - 4.0 * a * c;
//System.out.println("sqrtClause="+sqrtClause);
@ -1927,7 +1927,7 @@ public class Plane extends Vector {
b = 4.0 * F * H * J - 2.0 * G * I * J;
c = 4.0 * F * J * J - J * I * I;
//System.out.println("a="+a+" b="+b+" c="+c);
//System.out.println("a="+a+" b="+b+" zScaling="+zScaling);
double sqrtClause = b * b - 4.0 * a * c;
//System.out.println("sqrtClause="+sqrtClause);
if (Math.abs(sqrtClause) < MINIMUM_RESOLUTION_CUBED) {
@ -2044,7 +2044,7 @@ public class Plane extends Vector {
// ... can be solved by Cramer's rule:
// x = det(S0 R0 / S1 R1) / det(Q0 R0 / Q1 R1)
// y = det(Q0 S0 / Q1 S1) / det(Q0 R0 / Q1 R1)
// ... where det( a b / c d ) = ad - bc, so:
// ... where det( a b / zScaling d ) = ad - bc, so:
// x = (S0 * R1 - R0 * S1) / (Q0 * R1 - R0 * Q1)
// y = (Q0 * S1 - S0 * Q1) / (Q0 * R1 - R0 * Q1)
double x0;
@ -2088,16 +2088,16 @@ public class Plane extends Vector {
// Once an intersecting line is determined, the next step is to intersect that line with the ellipsoid, which
// will yield zero, one, or two points.
// The ellipsoid equation: 1,0 = x^2/a^2 + y^2/b^2 + z^2/c^2
// 1.0 = (At+A0)^2/a^2 + (Bt+B0)^2/b^2 + (Ct+C0)^2/c^2
// A^2 t^2 / a^2 + 2AA0t / a^2 + A0^2 / a^2 + B^2 t^2 / b^2 + 2BB0t / b^2 + B0^2 / b^2 + C^2 t^2 / c^2 + 2CC0t / c^2 + C0^2 / c^2 - 1,0 = 0.0
// [A^2 / a^2 + B^2 / b^2 + C^2 / c^2] t^2 + [2AA0 / a^2 + 2BB0 / b^2 + 2CC0 / c^2] t + [A0^2 / a^2 + B0^2 / b^2 + C0^2 / c^2 - 1,0] = 0.0
// The ellipsoid equation: 1,0 = x^2/a^2 + y^2/b^2 + z^2/zScaling^2
// 1.0 = (At+A0)^2/a^2 + (Bt+B0)^2/b^2 + (Ct+C0)^2/zScaling^2
// A^2 t^2 / a^2 + 2AA0t / a^2 + A0^2 / a^2 + B^2 t^2 / b^2 + 2BB0t / b^2 + B0^2 / b^2 + C^2 t^2 / zScaling^2 + 2CC0t / zScaling^2 + C0^2 / zScaling^2 - 1,0 = 0.0
// [A^2 / a^2 + B^2 / b^2 + C^2 / zScaling^2] t^2 + [2AA0 / a^2 + 2BB0 / b^2 + 2CC0 / zScaling^2] t + [A0^2 / a^2 + B0^2 / b^2 + C0^2 / zScaling^2 - 1,0] = 0.0
// Use the quadratic formula to determine t values and candidate point(s)
final double A = lineVectorX * lineVectorX * planetModel.inverseAbSquared +
lineVectorY * lineVectorY * planetModel.inverseAbSquared +
lineVectorZ * lineVectorZ * planetModel.inverseCSquared;
final double B = 2.0 * (lineVectorX * x0 * planetModel.inverseAbSquared + lineVectorY * y0 * planetModel.inverseAbSquared + lineVectorZ * z0 * planetModel.inverseCSquared);
final double C = x0 * x0 * planetModel.inverseAbSquared + y0 * y0 * planetModel.inverseAbSquared + z0 * z0 * planetModel.inverseCSquared - 1.0;
final double A = lineVectorX * lineVectorX * planetModel.inverseXYScalingSquared +
lineVectorY * lineVectorY * planetModel.inverseXYScalingSquared +
lineVectorZ * lineVectorZ * planetModel.inverseZScalingSquared;
final double B = 2.0 * (lineVectorX * x0 * planetModel.inverseXYScalingSquared + lineVectorY * y0 * planetModel.inverseXYScalingSquared + lineVectorZ * z0 * planetModel.inverseZScalingSquared);
final double C = x0 * x0 * planetModel.inverseXYScalingSquared + y0 * y0 * planetModel.inverseXYScalingSquared + z0 * z0 * planetModel.inverseZScalingSquared - 1.0;
final double BsquaredMinus = B * B - 4.0 * A * C;
if (Math.abs(BsquaredMinus) < MINIMUM_RESOLUTION_SQUARED) {
@ -2236,7 +2236,7 @@ public class Plane extends Vector {
// ... can be solved by Cramer's rule:
// x = det(S0 R0 / S1 R1) / det(Q0 R0 / Q1 R1)
// y = det(Q0 S0 / Q1 S1) / det(Q0 R0 / Q1 R1)
// ... where det( a b / c d ) = ad - bc, so:
// ... where det( a b / zScaling d ) = ad - bc, so:
// x = (S0 * R1 - R0 * S1) / (Q0 * R1 - R0 * Q1)
// y = (Q0 * S1 - S0 * Q1) / (Q0 * R1 - R0 * Q1)
double x0;
@ -2280,16 +2280,16 @@ public class Plane extends Vector {
// Once an intersecting line is determined, the next step is to intersect that line with the ellipsoid, which
// will yield zero, one, or two points.
// The ellipsoid equation: 1,0 = x^2/a^2 + y^2/b^2 + z^2/c^2
// 1.0 = (At+A0)^2/a^2 + (Bt+B0)^2/b^2 + (Ct+C0)^2/c^2
// A^2 t^2 / a^2 + 2AA0t / a^2 + A0^2 / a^2 + B^2 t^2 / b^2 + 2BB0t / b^2 + B0^2 / b^2 + C^2 t^2 / c^2 + 2CC0t / c^2 + C0^2 / c^2 - 1,0 = 0.0
// [A^2 / a^2 + B^2 / b^2 + C^2 / c^2] t^2 + [2AA0 / a^2 + 2BB0 / b^2 + 2CC0 / c^2] t + [A0^2 / a^2 + B0^2 / b^2 + C0^2 / c^2 - 1,0] = 0.0
// The ellipsoid equation: 1,0 = x^2/a^2 + y^2/b^2 + z^2/zScaling^2
// 1.0 = (At+A0)^2/a^2 + (Bt+B0)^2/b^2 + (Ct+C0)^2/zScaling^2
// A^2 t^2 / a^2 + 2AA0t / a^2 + A0^2 / a^2 + B^2 t^2 / b^2 + 2BB0t / b^2 + B0^2 / b^2 + C^2 t^2 / zScaling^2 + 2CC0t / zScaling^2 + C0^2 / zScaling^2 - 1,0 = 0.0
// [A^2 / a^2 + B^2 / b^2 + C^2 / zScaling^2] t^2 + [2AA0 / a^2 + 2BB0 / b^2 + 2CC0 / zScaling^2] t + [A0^2 / a^2 + B0^2 / b^2 + C0^2 / zScaling^2 - 1,0] = 0.0
// Use the quadratic formula to determine t values and candidate point(s)
final double A = lineVectorX * lineVectorX * planetModel.inverseAbSquared +
lineVectorY * lineVectorY * planetModel.inverseAbSquared +
lineVectorZ * lineVectorZ * planetModel.inverseCSquared;
final double B = 2.0 * (lineVectorX * x0 * planetModel.inverseAbSquared + lineVectorY * y0 * planetModel.inverseAbSquared + lineVectorZ * z0 * planetModel.inverseCSquared);
final double C = x0 * x0 * planetModel.inverseAbSquared + y0 * y0 * planetModel.inverseAbSquared + z0 * z0 * planetModel.inverseCSquared - 1.0;
final double A = lineVectorX * lineVectorX * planetModel.inverseXYScalingSquared +
lineVectorY * lineVectorY * planetModel.inverseXYScalingSquared +
lineVectorZ * lineVectorZ * planetModel.inverseZScalingSquared;
final double B = 2.0 * (lineVectorX * x0 * planetModel.inverseXYScalingSquared + lineVectorY * y0 * planetModel.inverseXYScalingSquared + lineVectorZ * z0 * planetModel.inverseZScalingSquared);
final double C = x0 * x0 * planetModel.inverseXYScalingSquared + y0 * y0 * planetModel.inverseXYScalingSquared + z0 * z0 * planetModel.inverseZScalingSquared - 1.0;
final double BsquaredMinus = B * B - 4.0 * A * C;
if (Math.abs(BsquaredMinus) < MINIMUM_RESOLUTION_SQUARED) {

493
lucene/spatial3d/src/java/org/apache/lucene/spatial3d/geom/PlanetModel.java
View File

@ -29,41 +29,47 @@ public class PlanetModel implements SerializableObject {
/** Planet model corresponding to sphere. */
public static final PlanetModel SPHERE = new PlanetModel(1.0,1.0);
/** Mean radius */
/** Planet model corresponding to WGS84 ellipsoid*/
// see http://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf
public static final double WGS84_MEAN = 6371008.7714;
/** Polar radius */
public static final double WGS84_POLAR = 6356752.314245;
/** Equatorial radius */
public static final double WGS84_EQUATORIAL = 6378137.0;
/** Planet model corresponding to WGS84 */
public static final PlanetModel WGS84 = new PlanetModel(WGS84_EQUATORIAL/WGS84_MEAN,
WGS84_POLAR/WGS84_MEAN);
public static final PlanetModel WGS84 = new PlanetModel(6378137.0d, 6356752.314245d, 298.257223563d);
/** Planet model corresponding to Clarke 1866 ellipsoid*/
// see https://georepository.com/ellipsoid_7008/Clarke-1866.html
public static final PlanetModel CLARKE_1866 = new PlanetModel(6378206.4d, 6356583.8d, 294.9786982d);
// Surface of the planet:
// x^2/a^2 + y^2/b^2 + z^2/c^2 = 1.0
// Scaling factors are a,b,c. geo3d can only support models where a==b, so use ab instead.
// x^2/a^2 + y^2/b^2 + z^2/zScaling^2 = 1.0
// Scaling factors are a,b,zScaling. geo3d can only support models where a==b, so use xyScaling instead.
/** Semi-major axis */
public double a;
/** Semi-minor axis */
public double b;
/** Flattening factor */
public double f;
/** The x/y scaling factor */
public final double ab;
public final double xyScaling;
/** The z scaling factor */
public final double c;
/** The inverse of ab */
public final double inverseAb;
/** The inverse of c */
public final double inverseC;
/** The square of the inverse of ab */
public final double inverseAbSquared;
/** The square of the inverse of c */
public final double inverseCSquared;
/** The flattening value */
public final double flattening;
public final double zScaling;
/** The inverse of xyScaling */
public final double inverseXYScaling;
/** The inverse of zScaling */
public final double inverseZScaling;
/** The square of the inverse of xyScaling */
public final double inverseXYScalingSquared;
/** The square of the inverse of zScaling */
public final double inverseZScalingSquared;
/** The scaled flattening value */
public final double scaledFlattening;
/** The square ratio */
public final double squareRatio;
/** The scale of the planet */
public final double scale;
/** The inverse of scale */
public final double inverseScale;
/** The mean radius of the planet model */
// Computed as (2a + b) / 3 from: "Geodetic Reference System 1980" by H. Moritz
// ftp://athena.fsv.cvut.cz/ZFG/grs80-Moritz.pdf
public double r1;
// We do NOT include radius, because all computations in geo3d are in radians, not meters.
@ -84,29 +90,79 @@ public class PlanetModel implements SerializableObject {
public final GeoPoint MAX_Y_POLE;
/** Minimum surface distance between poles */
public final double minimumPoleDistance;
/** Constructor.
* @param ab is the x/y scaling factor.
* @param c is the z scaling factor.
// ENCODING / DECODING CONSTANTS
/** bit space for integer encoding */
private static final int BITS = 32;
/** maximum magnitude value for *this* planet model */
public final double MAX_VALUE;
/** numeric space (buckets) for mapping double values into integer range */
private final double MUL;
/** scalar value used to decode from integer back into double space*/
public final double DECODE;
/** Max encoded value */
public final int MAX_ENCODED_VALUE;
/** Min encoded value */
public final int MIN_ENCODED_VALUE;
/** computed radians per meter for *this* planet model */
public final double RADIANS_PER_METER;
/** utility class used to encode/decode from lat/lon (decimal degrees) into doc value integers */
public final DocValueEncoder docValueEncoder;
/**
* Construct a Planet Model from the semi major axis, semi minor axis, and radius of curvature.
*
* NOTE: These parameters are only used to compute the scaling factors used in all computations. If
* scaling factors are already known, use {@code PlanetModel(final double xyScaling, final double zScaling)} instead
*
* @param semiMajorAxis is the semi major axis (in meters) defined as 'a' in projection formulae.
* @param semiMinorAxis is the semi minor axis (in meters) defined as 'b' in projection formulae.
* @param flattening is the flattening factor defined as a - b / a.
*/
public PlanetModel(final double ab, final double c) {
this.ab = ab;
this.c = c;
this.inverseAb = 1.0 / ab;
this.inverseC = 1.0 / c;
this.flattening = (ab - c) * inverseAb;
this.squareRatio = (ab * ab - c * c) / (c * c);
this.inverseAbSquared = inverseAb * inverseAb;
this.inverseCSquared = inverseC * inverseC;
this.NORTH_POLE = new GeoPoint(c, 0.0, 0.0, 1.0, Math.PI * 0.5, 0.0);
this.SOUTH_POLE = new GeoPoint(c, 0.0, 0.0, -1.0, -Math.PI * 0.5, 0.0);
this.MIN_X_POLE = new GeoPoint(ab, -1.0, 0.0, 0.0, 0.0, -Math.PI);
this.MAX_X_POLE = new GeoPoint(ab, 1.0, 0.0, 0.0, 0.0, 0.0);
this.MIN_Y_POLE = new GeoPoint(ab, 0.0, -1.0, 0.0, 0.0, -Math.PI * 0.5);
this.MAX_Y_POLE = new GeoPoint(ab, 0.0, 1.0, 0.0, 0.0, Math.PI * 0.5);
this.scale = (2.0 * ab + c)/3.0;
public PlanetModel(final double semiMajorAxis, final double semiMinorAxis, final double flattening) { // final double zScaling) {
this(semiMajorAxis, semiMinorAxis, flattening, ((2d * semiMajorAxis) + semiMinorAxis) / 3d);
}
private PlanetModel(final double a, final double b, final double f, final double r1) {
this(a / r1, b / r1);
this.a = a;
this.b = b;
this.f = f;
this.r1 = r1;
}
/** Constructor.
* @param xyScaling is the x/y scaling factor.
* @param zScaling is the z scaling factor.
*/
public PlanetModel(final double xyScaling, final double zScaling) {
this.xyScaling = xyScaling;
this.zScaling = zScaling;
this.inverseXYScaling = 1.0 / xyScaling;
this.inverseZScaling = 1.0 / zScaling;
this.scaledFlattening = (xyScaling - zScaling) * inverseXYScaling;
this.squareRatio = (xyScaling * xyScaling - zScaling * zScaling) / (zScaling * zScaling);
this.inverseXYScalingSquared = inverseXYScaling * inverseXYScaling;
this.inverseZScalingSquared = inverseZScaling * inverseZScaling;
this.NORTH_POLE = new GeoPoint(zScaling, 0.0, 0.0, 1.0, Math.PI * 0.5, 0.0);
this.SOUTH_POLE = new GeoPoint(zScaling, 0.0, 0.0, -1.0, -Math.PI * 0.5, 0.0);
this.MIN_X_POLE = new GeoPoint(xyScaling, -1.0, 0.0, 0.0, 0.0, -Math.PI);
this.MAX_X_POLE = new GeoPoint(xyScaling, 1.0, 0.0, 0.0, 0.0, 0.0);
this.MIN_Y_POLE = new GeoPoint(xyScaling, 0.0, -1.0, 0.0, 0.0, -Math.PI * 0.5);
this.MAX_Y_POLE = new GeoPoint(xyScaling, 0.0, 1.0, 0.0, 0.0, Math.PI * 0.5);
this.scale = (2.0 * xyScaling + zScaling)/3.0;
this.inverseScale = 1.0 / scale;
this.minimumPoleDistance = Math.min(surfaceDistance(NORTH_POLE, SOUTH_POLE), surfaceDistance(MIN_X_POLE, MAX_X_POLE));
this.MAX_VALUE = getMaximumMagnitude();
this.MUL = (0x1L << BITS) / (2 * this.MAX_VALUE);
this.DECODE = getNextSafeDouble(1/MUL);
this.MIN_ENCODED_VALUE = encodeValue(-MAX_VALUE);
this.MAX_ENCODED_VALUE = encodeValue(MAX_VALUE);
this.RADIANS_PER_METER = 1.0 / xyScaling;
this.docValueEncoder = new DocValueEncoder(this);
}
/** Deserialization constructor.
@ -118,71 +174,141 @@ public class PlanetModel implements SerializableObject {
@Override
public void write(final OutputStream outputStream) throws IOException {
SerializableObject.writeDouble(outputStream, ab);
SerializableObject.writeDouble(outputStream, c);
SerializableObject.writeDouble(outputStream, xyScaling);
SerializableObject.writeDouble(outputStream, zScaling);
}
/** Does this planet model describe a sphere?
*@return true if so.
*/
public boolean isSphere() {
return this.ab == this.c;
return this.xyScaling == this.zScaling;
}
/** Find the minimum magnitude of all points on the ellipsoid.
* @return the minimum magnitude for the planet.
*/
public double getMinimumMagnitude() {
return Math.min(this.ab, this.c);
return Math.min(this.xyScaling, this.zScaling);
}
/** Find the maximum magnitude of all points on the ellipsoid.
* @return the maximum magnitude for the planet.
*/
public double getMaximumMagnitude() {
return Math.max(this.ab, this.c);
return Math.max(this.xyScaling, this.zScaling);
}
/** Find the minimum x value.
*@return the minimum X value.
*/
public double getMinimumXValue() {
return -this.ab;
return -this.xyScaling;
}
/** Find the maximum x value.
*@return the maximum X value.
*/
public double getMaximumXValue() {
return this.ab;
return this.xyScaling;
}
/** Find the minimum y value.
*@return the minimum Y value.
*/
public double getMinimumYValue() {
return -this.ab;
return -this.xyScaling;
}
/** Find the maximum y value.
*@return the maximum Y value.
*/
public double getMaximumYValue() {
return this.ab;
return this.xyScaling;
}
/** Find the minimum z value.
*@return the minimum Z value.
*/
public double getMinimumZValue() {
return -this.c;
return -this.zScaling;
}
/** Find the maximum z value.
*@return the maximum Z value.
*/
public double getMaximumZValue() {
return this.c;
return this.zScaling;
}
/** return the calculated mean radius (in meters) */
public double getMeanRadiusMeters() {
return this.r1;
}
/** encode the provided value from double to integer space */
public int encodeValue(double x) {
if (x > getMaximumMagnitude()) {
throw new IllegalArgumentException("value=" + x + " is out-of-bounds (greater than planetMax=" + getMaximumMagnitude() + ")");
}
if (x == getMaximumMagnitude()) {
x = Math.nextDown(x);
}
if (x < -getMaximumMagnitude()) {
throw new IllegalArgumentException("value=" + x + " is out-of-bounds (less than than -planetMax=" + -getMaximumMagnitude() + ")");
}
long result = (long) Math.floor(x / DECODE);
assert result >= Integer.MIN_VALUE;
assert result <= Integer.MAX_VALUE;
return (int) result;
}
/**
* Decodes a given integer back into the radian value according to the defined planet model
*/
public double decodeValue(int x) {
double result;
if (x == MIN_ENCODED_VALUE) {
// We must special case this, because -MAX_VALUE is not guaranteed to land precisely at a floor value, and we don't ever want to
// return a value outside of the planet's range (I think?). The max value is "safe" because we floor during encode:
result = -MAX_VALUE;
} else if (x == MAX_ENCODED_VALUE) {
result = MAX_VALUE;
} else {
// We decode to the center value; this keeps the encoding stable
result = (x+0.5) * DECODE;
}
assert result >= -MAX_VALUE && result <= MAX_VALUE;
return result;
}
/** return reference to the DocValueEncoder used to encode/decode Geo3DDocValues */
public DocValueEncoder getDocValueEncoder() {
return this.docValueEncoder;
}
/** Returns a double value >= x such that if you multiply that value by an int, and then
* divide it by that int again, you get precisely the same value back */
private static double getNextSafeDouble(double x) {
// Move to double space:
long bits = Double.doubleToLongBits(x);
// Make sure we are beyond the actual maximum value:
bits += Integer.MAX_VALUE;
// Clear the bottom 32 bits:
bits &= ~((long) Integer.MAX_VALUE);
// Convert back to double:
double result = Double.longBitsToDouble(bits);
assert result >= x;
return result;
}
/** Converts earth-surface meters to radians */
public double fromMeters(final double meters) {
return meters * RADIANS_PER_METER;
}
/** Check if point is on surface.
@ -200,8 +326,8 @@ public class PlanetModel implements SerializableObject {
*/
public boolean pointOnSurface(final double x, final double y, final double z) {
// Equation of planet surface is:
// x^2 / a^2 + y^2 / b^2 + z^2 / c^2 - 1 = 0
return Math.abs(x * x * inverseAb * inverseAb + y * y * inverseAb * inverseAb + z * z * inverseC * inverseC - 1.0) < Vector.MINIMUM_RESOLUTION;
// x^2 / a^2 + y^2 / b^2 + z^2 / zScaling^2 - 1 = 0
return Math.abs(x * x * inverseXYScaling * inverseXYScaling + y * y * inverseXYScaling * inverseXYScaling + z * z * inverseZScaling * inverseZScaling - 1.0) < Vector.MINIMUM_RESOLUTION;
}
/** Check if point is outside surface.
@ -219,8 +345,8 @@ public class PlanetModel implements SerializableObject {
*/
public boolean pointOutside(final double x, final double y, final double z) {
// Equation of planet surface is:
// x^2 / a^2 + y^2 / b^2 + z^2 / c^2 - 1 = 0
return (x * x + y * y) * inverseAb * inverseAb + z * z * inverseC * inverseC - 1.0 > Vector.MINIMUM_RESOLUTION;
// x^2 / a^2 + y^2 / b^2 + z^2 / zScaling^2 - 1 = 0
return (x * x + y * y) * inverseXYScaling * inverseXYScaling + z * z * inverseZScaling * inverseZScaling - 1.0 > Vector.MINIMUM_RESOLUTION;
}
/** Compute a GeoPoint that's scaled to actually be on the planet surface.
@ -239,12 +365,12 @@ public class PlanetModel implements SerializableObject {
*/
public GeoPoint createSurfacePoint(final double x, final double y, final double z) {
// The equation of the surface is:
// (x^2 / a^2 + y^2 / b^2 + z^2 / c^2) = 1
// (x^2 / a^2 + y^2 / b^2 + z^2 / zScaling^2) = 1
// We will need to scale the passed-in x, y, z values:
// ((tx)^2 / a^2 + (ty)^2 / b^2 + (tz)^2 / c^2) = 1
// t^2 * (x^2 / a^2 + y^2 / b^2 + z^2 / c^2) = 1
// t = sqrt ( 1 / (x^2 / a^2 + y^2 / b^2 + z^2 / c^2))
final double t = Math.sqrt(1.0 / (x*x*inverseAbSquared + y*y*inverseAbSquared + z*z*inverseCSquared));
// ((tx)^2 / a^2 + (ty)^2 / b^2 + (tz)^2 / zScaling^2) = 1
// t^2 * (x^2 / a^2 + y^2 / b^2 + z^2 / zScaling^2) = 1
// t = sqrt ( 1 / (x^2 / a^2 + y^2 / b^2 + z^2 / zScaling^2))
final double t = Math.sqrt(1.0 / (x*x* inverseXYScalingSquared + y*y* inverseXYScalingSquared + z*z* inverseZScalingSquared));
return new GeoPoint(t*x, t*y, t*z);
}
@ -258,9 +384,9 @@ public class PlanetModel implements SerializableObject {
final double B0 = (pt1.y + pt2.y) * 0.5;
final double C0 = (pt1.z + pt2.z) * 0.5;
final double denom = inverseAbSquared * A0 * A0 +
inverseAbSquared * B0 * B0 +
inverseCSquared * C0 * C0;
final double denom = inverseXYScalingSquared * A0 * A0 +
inverseXYScalingSquared * B0 * B0 +
inverseZScalingSquared * C0 * C0;
if(denom < Vector.MINIMUM_RESOLUTION) {
// Bisection is undefined
@ -280,8 +406,8 @@ public class PlanetModel implements SerializableObject {
*/
public double surfaceDistance(final GeoPoint pt1, final GeoPoint pt2) {
final double L = pt2.getLongitude() - pt1.getLongitude();
final double U1 = Math.atan((1.0-flattening) * Math.tan(pt1.getLatitude()));
final double U2 = Math.atan((1.0-flattening) * Math.tan(pt2.getLatitude()));
final double U1 = Math.atan((1.0- scaledFlattening) * Math.tan(pt1.getLatitude()));
final double U2 = Math.atan((1.0- scaledFlattening) * Math.tan(pt2.getLatitude()));
final double sinU1 = Math.sin(U1);
final double cosU1 = Math.cos(U1);
@ -324,9 +450,9 @@ public class PlanetModel implements SerializableObject {
if (Double.isNaN(cos2SigmaM))
cos2SigmaM = 0.0; // equatorial line: cosSqAlpha=0
C = flattening / 16.0 * cosSqAlpha * (4.0 + flattening * (4.0 - 3.0 * cosSqAlpha));
C = scaledFlattening / 16.0 * cosSqAlpha * (4.0 + scaledFlattening * (4.0 - 3.0 * cosSqAlpha));
lambdaP = lambda;
lambda = L + (1.0 - C) * flattening * sinAlpha *
lambda = L + (1.0 - C) * scaledFlattening * sinAlpha *
(sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1.0 + 2.0 * cos2SigmaM *cos2SigmaM)));
} while (Math.abs(lambda-lambdaP) >= Vector.MINIMUM_RESOLUTION && ++iterLimit < 100);
final double uSq = cosSqAlpha * this.squareRatio;
@ -335,7 +461,7 @@ public class PlanetModel implements SerializableObject {
final double deltaSigma = B * sinSigma * (cos2SigmaM + B / 4.0 * (cosSigma * (-1.0 + 2.0 * cos2SigmaM * cos2SigmaM)-
B / 6.0 * cos2SigmaM * (-3.0 + 4.0 * sinSigma * sinSigma) * (-3.0 + 4.0 * cos2SigmaM * cos2SigmaM)));
return c * inverseScale * A * (sigma - deltaSigma);
return zScaling * inverseScale * A * (sigma - deltaSigma);
}
/** Compute new point given original point, a bearing direction, and an adjusted angle (as would be computed by
@ -357,7 +483,7 @@ public class PlanetModel implements SerializableObject {
double sinα1 = Math.sin(bearing);
double cosα1 = Math.cos(bearing);
double tanU1 = (1.0 - flattening) * Math.tan(lat);
double tanU1 = (1.0 - scaledFlattening) * Math.tan(lat);
double cosU1 = 1.0 / Math.sqrt((1.0 + tanU1 * tanU1));
double sinU1 = tanU1 * cosU1;
@ -373,7 +499,7 @@ public class PlanetModel implements SerializableObject {
double cosσ;
double Δσ;
double σ = dist / (c * inverseScale * A);
double σ = dist / (zScaling * inverseScale * A);
double σʹ;
double iterations = 0;
do {
@ -383,30 +509,233 @@ public class PlanetModel implements SerializableObject {
Δσ = B * sinσ * (cos2σM + B / 4.0 * (cosσ * (-1.0 + 2.0 * cos2σM * cos2σM) -
B / 6.0 * cos2σM * (-3.0 + 4.0 * sinσ * sinσ) * (-3.0 + 4.0 * cos2σM * cos2σM)));
σʹ = σ;
σ = dist / (c * inverseScale * A) + Δσ;
σ = dist / (zScaling * inverseScale * A) + Δσ;
} while (Math.abs(σ - σʹ) >= Vector.MINIMUM_RESOLUTION && ++iterations < 100);
double x = sinU1 * sinσ - cosU1 * cosσ * cosα1;
double φ2 = Math.atan2(sinU1 * cosσ + cosU1 * sinσ * cosα1, (1.0 - flattening) * Math.sqrt(sinα * sinα + x * x));
double φ2 = Math.atan2(sinU1 * cosσ + cosU1 * sinσ * cosα1, (1.0 - scaledFlattening) * Math.sqrt(sinα * sinα + x * x));
double λ = Math.atan2(sinσ * sinα1, cosU1 * cosσ - sinU1 * sinσ * cosα1);
double C = flattening / 16.0 * cosSqα * (4.0 + flattening * (4.0 - 3.0 * cosSqα));
double L = λ - (1.0 - C) * flattening * sinα *
double C = scaledFlattening / 16.0 * cosSqα * (4.0 + scaledFlattening * (4.0 - 3.0 * cosSqα));
double L = λ - (1.0 - C) * scaledFlattening * sinα *
(σ + C * sinσ * (cos2σM + C * cosσ * (-1.0 + 2.0 * cos2σM * cos2σM)));
double λ2 = (lon + L + 3.0 * Math.PI) % (2.0 * Math.PI) - Math.PI; // normalise to -180..+180
return new GeoPoint(this, φ2, λ2);
}
/** Utility class for encoding / decoding from lat/lon (decimal degrees) into sortable doc value numerics (integers) */
public static class DocValueEncoder {
private final PlanetModel planetModel;
// These are the multiplicative constants we need to use to arrive at values that fit in 21 bits.
// The formula we use to go from double to encoded value is: Math.floor((value - minimum) * factor + 0.5)
// If we plug in maximum for value, we should get 0x1FFFFF.
// So, 0x1FFFFF = Math.floor((maximum - minimum) * factor + 0.5)
// We factor out the 0.5 and Math.floor by stating instead:
// 0x1FFFFF = (maximum - minimum) * factor
// So, factor = 0x1FFFFF / (maximum - minimum)
private final static double inverseMaximumValue = 1.0 / (double)(0x1FFFFF);
private final double inverseXFactor;
private final double inverseYFactor;
private final double inverseZFactor;
private final double xFactor;
private final double yFactor;
private final double zFactor;
// Fudge factor for step adjustments. This is here solely to handle inaccuracies in bounding boxes
// that occur because of quantization. For unknown reasons, the fudge factor needs to be
// 10.0 rather than 1.0. See LUCENE-7430.
private final static double STEP_FUDGE = 10.0;
// These values are the delta between a value and the next value in each specific dimension
private final double xStep;
private final double yStep;
private final double zStep;
/** construct an encoder/decoder instance from the provided PlanetModel definition */
private DocValueEncoder(final PlanetModel planetModel) {
this.planetModel = planetModel;
this.inverseXFactor = (planetModel.getMaximumXValue() - planetModel.getMinimumXValue()) * inverseMaximumValue;
this.inverseYFactor = (planetModel.getMaximumYValue() - planetModel.getMinimumYValue()) * inverseMaximumValue;
this.inverseZFactor = (planetModel.getMaximumZValue() - planetModel.getMinimumZValue()) * inverseMaximumValue;
this.xFactor = 1.0 / inverseXFactor;
this.yFactor = 1.0 / inverseYFactor;
this.zFactor = 1.0 / inverseZFactor;
this.xStep = inverseXFactor * STEP_FUDGE;
this.yStep = inverseYFactor * STEP_FUDGE;
this.zStep = inverseZFactor * STEP_FUDGE;
}
/** Encode a point.
* @param point is the point
* @return the encoded long
*/
public long encodePoint(final GeoPoint point) {
return encodePoint(point.x, point.y, point.z);
}
/** Encode a point.
* @param x is the x value
* @param y is the y value
* @param z is the z value
* @return the encoded long
*/
public long encodePoint(final double x, final double y, final double z) {
int XEncoded = encodeX(x);
int YEncoded = encodeY(y);
int ZEncoded = encodeZ(z);
return
(((long)(XEncoded & 0x1FFFFF)) << 42) |
(((long)(YEncoded & 0x1FFFFF)) << 21) |
((long)(ZEncoded & 0x1FFFFF));
}
/** Decode GeoPoint value from long docvalues value.
* @param docValue is the doc values value.
* @return the GeoPoint.
*/
public GeoPoint decodePoint(final long docValue) {
return new GeoPoint(decodeX(((int)(docValue >> 42)) & 0x1FFFFF),
decodeY(((int)(docValue >> 21)) & 0x1FFFFF),
decodeZ(((int)(docValue)) & 0x1FFFFF));
}
/** Decode X value from long docvalues value.
* @param docValue is the doc values value.
* @return the x value.
*/
public double decodeXValue(final long docValue) {
return decodeX(((int)(docValue >> 42)) & 0x1FFFFF);
}
/** Decode Y value from long docvalues value.
* @param docValue is the doc values value.
* @return the y value.
*/
public double decodeYValue(final long docValue) {
return decodeY(((int)(docValue >> 21)) & 0x1FFFFF);
}
/** Decode Z value from long docvalues value.
* @param docValue is the doc values value.
* @return the z value.
*/
public double decodeZValue(final long docValue) {
return decodeZ(((int)(docValue)) & 0x1FFFFF);
}
/** Round the provided X value down, by encoding it, decrementing it, and unencoding it.
* @param startValue is the starting value.
* @return the rounded value.
*/
public double roundDownX(final double startValue) {
return startValue - xStep;
}
/** Round the provided X value up, by encoding it, incrementing it, and unencoding it.
* @param startValue is the starting value.
* @return the rounded value.
*/
public double roundUpX(final double startValue) {
return startValue + xStep;
}
/** Round the provided Y value down, by encoding it, decrementing it, and unencoding it.
* @param startValue is the starting value.
* @return the rounded value.
*/
public double roundDownY(final double startValue) {
return startValue - yStep;
}
/** Round the provided Y value up, by encoding it, incrementing it, and unencoding it.
* @param startValue is the starting value.
* @return the rounded value.
*/
public double roundUpY(final double startValue) {
return startValue + yStep;
}
/** Round the provided Z value down, by encoding it, decrementing it, and unencoding it.
* @param startValue is the starting value.
* @return the rounded value.
*/
public double roundDownZ(final double startValue) {
return startValue - zStep;
}
/** Round the provided Z value up, by encoding it, incrementing it, and unencoding it.
* @param startValue is the starting value.
* @return the rounded value.
*/
public double roundUpZ(final double startValue) {
return startValue + zStep;
}
// For encoding/decoding, we generally want the following behavior:
// (1) If you encode the maximum value or the minimum value, the resulting int fits in 21 bits.
// (2) If you decode an encoded value, you get back the original value for both the minimum and maximum planet model values.
// (3) Rounding occurs such that a small delta from the minimum and maximum planet model values still returns the same
// values -- that is, these are in the center of the range of input values that should return the minimum or maximum when decoded
private int encodeX(final double x) {
if (x > planetModel.getMaximumXValue()) {
throw new IllegalArgumentException("x value exceeds planet model maximum");
} else if (x < planetModel.getMinimumXValue()) {
throw new IllegalArgumentException("x value less than planet model minimum");
}
return (int)Math.floor((x - planetModel.getMinimumXValue()) * xFactor + 0.5);
}
private double decodeX(final int x) {
return x * inverseXFactor + planetModel.getMinimumXValue();
}
private int encodeY(final double y) {
if (y > planetModel.getMaximumYValue()) {
throw new IllegalArgumentException("y value exceeds planet model maximum");
} else if (y < planetModel.getMinimumYValue()) {
throw new IllegalArgumentException("y value less than planet model minimum");
}
return (int)Math.floor((y - planetModel.getMinimumYValue()) * yFactor + 0.5);
}
private double decodeY(final int y) {
return y * inverseYFactor + planetModel.getMinimumYValue();
}
private int encodeZ(final double z) {
if (z > planetModel.getMaximumZValue()) {
throw new IllegalArgumentException("z value exceeds planet model maximum");
} else if (z < planetModel.getMinimumZValue()) {
throw new IllegalArgumentException("z value less than planet model minimum");
}
return (int)Math.floor((z - planetModel.getMinimumZValue()) * zFactor + 0.5);
}
private double decodeZ(final int z) {
return z * inverseZFactor + planetModel.getMinimumZValue();
}
}
@Override
public boolean equals(final Object o) {
if (!(o instanceof PlanetModel))
return false;
final PlanetModel other = (PlanetModel)o;
return ab == other.ab && c == other.c;
return xyScaling == other.xyScaling && zScaling == other.zScaling;
}
@Override
public int hashCode() {
return Double.hashCode(ab) + Double.hashCode(c);
return Double.hashCode(xyScaling) + Double.hashCode(zScaling);
}
@Override
@ -415,8 +744,10 @@ public class PlanetModel implements SerializableObject {
return "PlanetModel.SPHERE";
} else if (this.equals(WGS84)) {
return "PlanetModel.WGS84";
} else if (this.equals(CLARKE_1866)) {
return "PlanetModel.CLARKE_1866";
} else {
return "PlanetModel(ab="+ab+" c="+c+")";
return "PlanetModel(xyScaling="+ xyScaling +" zScaling="+ zScaling +")";
}
}
}

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

@ -562,7 +562,7 @@ public class Vector {
* @return a magnitude value for that (x,y,z) that projects the vector onto the specified ellipsoid.
*/
static double computeDesiredEllipsoidMagnitude(final PlanetModel planetModel, final double x, final double y, final double z) {
return 1.0 / Math.sqrt(x*x*planetModel.inverseAbSquared + y*y*planetModel.inverseAbSquared + z*z*planetModel.inverseCSquared);
return 1.0 / Math.sqrt(x*x*planetModel.inverseXYScalingSquared + y*y*planetModel.inverseXYScalingSquared + z*z*planetModel.inverseZScalingSquared);
}
/** Compute the desired magnitude of a unit vector projected to a given
@ -572,7 +572,7 @@ public class Vector {
* @return a magnitude value for that z value that projects the vector onto the specified ellipsoid.
*/
static double computeDesiredEllipsoidMagnitude(final PlanetModel planetModel, final double z) {
return 1.0 / Math.sqrt((1.0-z*z)*planetModel.inverseAbSquared + z*z*planetModel.inverseCSquared);
return 1.0 / Math.sqrt((1.0-z*z)*planetModel.inverseXYScalingSquared + z*z*planetModel.inverseZScalingSquared);
}
@Override

20
lucene/spatial3d/src/test/org/apache/lucene/spatial3d/TestGeo3DDocValues.java
View File

@ -16,6 +16,7 @@
*/
package org.apache.lucene.spatial3d;
import com.carrotsearch.randomizedtesting.generators.RandomPicks;
import org.apache.lucene.spatial3d.geom.GeoPoint;
import org.apache.lucene.spatial3d.geom.PlanetModel;
@ -35,17 +36,18 @@ public class TestGeo3DDocValues extends LuceneTestCase {
}
void checkPointEncoding(final double latitude, final double longitude) {
final GeoPoint point = new GeoPoint(PlanetModel.WGS84, Geo3DUtil.fromDegrees(latitude), Geo3DUtil.fromDegrees(longitude));
long pointValue = Geo3DDocValuesField.encodePoint(point);
final double x = Geo3DDocValuesField.decodeXValue(pointValue);
final double y = Geo3DDocValuesField.decodeYValue(pointValue);
final double z = Geo3DDocValuesField.decodeZValue(pointValue);
PlanetModel planetModel = RandomPicks.randomFrom(random(), new PlanetModel[] {PlanetModel.WGS84, PlanetModel.CLARKE_1866});
final GeoPoint point = new GeoPoint(planetModel, Geo3DUtil.fromDegrees(latitude), Geo3DUtil.fromDegrees(longitude));
long pointValue = planetModel.getDocValueEncoder().encodePoint(point);
final double x = planetModel.getDocValueEncoder().decodeXValue(pointValue);
final double y = planetModel.getDocValueEncoder().decodeYValue(pointValue);
final double z = planetModel.getDocValueEncoder().decodeZValue(pointValue);
final GeoPoint pointR = new GeoPoint(x,y,z);
// Check whether stable
pointValue = Geo3DDocValuesField.encodePoint(x, y, z);
assertEquals(x, Geo3DDocValuesField.decodeXValue(pointValue), 0.0);
assertEquals(y, Geo3DDocValuesField.decodeYValue(pointValue), 0.0);
assertEquals(z, Geo3DDocValuesField.decodeZValue(pointValue), 0.0);
pointValue = planetModel.getDocValueEncoder().encodePoint(x, y, z);
assertEquals(x, planetModel.getDocValueEncoder().decodeXValue(pointValue), 0.0);
assertEquals(y, planetModel.getDocValueEncoder().decodeYValue(pointValue), 0.0);
assertEquals(z, planetModel.getDocValueEncoder().decodeZValue(pointValue), 0.0);
// Check whether has some relationship with original point
assertEquals(0.0, point.arcDistance(pointR), 0.02);
}

270
lucene/spatial3d/src/test/org/apache/lucene/spatial3d/TestGeo3DPoint.java
View File

@ -26,6 +26,7 @@ import java.util.List;
import java.util.Random;
import java.util.Set;
import com.carrotsearch.randomizedtesting.generators.RandomPicks;
import org.apache.lucene.codecs.Codec;
import org.apache.lucene.codecs.FilterCodec;
import org.apache.lucene.codecs.PointsFormat;
@ -84,6 +85,10 @@ import com.carrotsearch.randomizedtesting.generators.RandomNumbers;
public class TestGeo3DPoint extends LuceneTestCase {
protected PlanetModel randomPlanetModel() {
return RandomPicks.randomFrom(random(), new PlanetModel[] {PlanetModel.WGS84, PlanetModel.CLARKE_1866, PlanetModel.SPHERE});
}
private static Codec getCodec() {
if (Codec.getDefault().getName().equals("Lucene84")) {
int maxPointsInLeafNode = TestUtil.nextInt(random(), 16, 2048);
@ -119,13 +124,14 @@ public class TestGeo3DPoint extends LuceneTestCase {
iwc.setCodec(getCodec());
IndexWriter w = new IndexWriter(dir, iwc);
Document doc = new Document();
doc.add(new Geo3DPoint("field", 50.7345267, -97.5303555));
PlanetModel planetModel = randomPlanetModel();
doc.add(new Geo3DPoint("field", planetModel, 50.7345267, -97.5303555));
w.addDocument(doc);
IndexReader r = DirectoryReader.open(w);
// We can't wrap with "exotic" readers because the query must see the BKD3DDVFormat:
IndexSearcher s = newSearcher(r, false);
assertEquals(1, s.search(Geo3DPoint.newShapeQuery("field",
GeoCircleFactory.makeGeoCircle(PlanetModel.WGS84, toRadians(50), toRadians(-97), Math.PI/180.)), 1).totalHits.value);
GeoCircleFactory.makeGeoCircle(planetModel, toRadians(50), toRadians(-97), Math.PI/180.)), 1).totalHits.value);
w.close();
r.close();
dir.close();
@ -144,11 +150,13 @@ public class TestGeo3DPoint extends LuceneTestCase {
final int yMinEnc, yMaxEnc;
final int zMinEnc, zMaxEnc;
final int splitCount;
final PlanetModel planetModel;
public Cell(Cell parent,
int xMinEnc, int xMaxEnc,
int yMinEnc, int yMaxEnc,
int zMinEnc, int zMaxEnc,
final PlanetModel planetModel,
int splitCount) {
this.parent = parent;
this.xMinEnc = xMinEnc;
@ -159,13 +167,14 @@ public class TestGeo3DPoint extends LuceneTestCase {
this.zMaxEnc = zMaxEnc;
this.cellID = nextCellID++;
this.splitCount = splitCount;
this.planetModel = planetModel;
}
/** Returns true if the quantized point lies within this cell, inclusive on all bounds. */
public boolean contains(GeoPoint point) {
int docX = Geo3DUtil.encodeValue(point.x);
int docY = Geo3DUtil.encodeValue(point.y);
int docZ = Geo3DUtil.encodeValue(point.z);
int docX = planetModel.encodeValue(point.x);
int docY = planetModel.encodeValue(point.y);
int docZ = planetModel.encodeValue(point.z);
return docX >= xMinEnc && docX <= xMaxEnc &&
docY >= yMinEnc && docY <= yMaxEnc &&
@ -178,12 +187,12 @@ public class TestGeo3DPoint extends LuceneTestCase {
}
}
private static double quantize(double xyzValue) {
return Geo3DUtil.decodeValue(Geo3DUtil.encodeValue(xyzValue));
private static double quantize(double xyzValue, final PlanetModel planetModel) {
return planetModel.decodeValue(planetModel.encodeValue(xyzValue));
}
private static GeoPoint quantize(GeoPoint point) {
return new GeoPoint(quantize(point.x), quantize(point.y), quantize(point.z));
private static GeoPoint quantize(GeoPoint point, final PlanetModel planetModel) {
return new GeoPoint(quantize(point.x, planetModel), quantize(point.y, planetModel), quantize(point.z, planetModel));
}
/** Tests consistency of GeoArea.getRelationship vs GeoShape.isWithin */
@ -196,9 +205,10 @@ public class TestGeo3DPoint extends LuceneTestCase {
GeoPoint[] docs = new GeoPoint[numDocs];
GeoPoint[] unquantizedDocs = new GeoPoint[numDocs];
PlanetModel planetModel = PlanetModel.CLARKE_1866;
for(int docID=0;docID<numDocs;docID++) {
unquantizedDocs[docID] = new GeoPoint(PlanetModel.WGS84, toRadians(GeoTestUtil.nextLatitude()), toRadians(GeoTestUtil.nextLongitude()));
docs[docID] = quantize(unquantizedDocs[docID]);
unquantizedDocs[docID] = new GeoPoint(planetModel, toRadians(GeoTestUtil.nextLatitude()), toRadians(GeoTestUtil.nextLongitude()));
docs[docID] = quantize(unquantizedDocs[docID], planetModel);
if (VERBOSE) {
System.out.println(" doc=" + docID + ": " + docs[docID] + "; unquantized: "+unquantizedDocs[docID]);
}
@ -209,7 +219,7 @@ public class TestGeo3DPoint extends LuceneTestCase {
int recurseDepth = RandomNumbers.randomIntBetween(random(), 5, 15);
for(int iter=0;iter<iters;iter++) {
GeoShape shape = randomShape();
GeoShape shape = randomShape(planetModel);
StringWriter sw = new StringWriter();
PrintWriter log = new PrintWriter(sw, true);
@ -223,12 +233,13 @@ public class TestGeo3DPoint extends LuceneTestCase {
// Start with the root cell that fully contains the shape:
Cell root = new Cell(null,
encodeValueLenient(bounds.getMinimumX()),
encodeValueLenient(bounds.getMaximumX()),
encodeValueLenient(bounds.getMinimumY()),
encodeValueLenient(bounds.getMaximumY()),
encodeValueLenient(bounds.getMinimumZ()),
encodeValueLenient(bounds.getMaximumZ()),
encodeValueLenient(bounds.getMinimumX(), planetModel),
encodeValueLenient(bounds.getMaximumX(), planetModel),
encodeValueLenient(bounds.getMinimumY(), planetModel),
encodeValueLenient(bounds.getMaximumY(), planetModel),
encodeValueLenient(bounds.getMinimumZ(), planetModel),
encodeValueLenient(bounds.getMaximumZ(), planetModel),
planetModel,
0);
if (VERBOSE) {
@ -293,16 +304,15 @@ public class TestGeo3DPoint extends LuceneTestCase {
}
}
} else {
GeoArea xyzSolid = GeoAreaFactory.makeGeoArea(PlanetModel.WGS84,
Geo3DUtil.decodeValueFloor(cell.xMinEnc), Geo3DUtil.decodeValueCeil(cell.xMaxEnc),
Geo3DUtil.decodeValueFloor(cell.yMinEnc), Geo3DUtil.decodeValueCeil(cell.yMaxEnc),
Geo3DUtil.decodeValueFloor(cell.zMinEnc), Geo3DUtil.decodeValueCeil(cell.zMaxEnc));
GeoArea xyzSolid = GeoAreaFactory.makeGeoArea(planetModel,
Geo3DUtil.decodeValueFloor(cell.xMinEnc, planetModel), Geo3DUtil.decodeValueCeil(cell.xMaxEnc, planetModel),
Geo3DUtil.decodeValueFloor(cell.yMinEnc, planetModel), Geo3DUtil.decodeValueCeil(cell.yMaxEnc, planetModel),
Geo3DUtil.decodeValueFloor(cell.zMinEnc, planetModel), Geo3DUtil.decodeValueCeil(cell.zMaxEnc, planetModel));
if (VERBOSE) {
log.println(" minx="+Geo3DUtil.decodeValueFloor(cell.xMinEnc)+" maxx="+Geo3DUtil.decodeValueCeil(cell.xMaxEnc)+
" miny="+Geo3DUtil.decodeValueFloor(cell.yMinEnc)+" maxy="+Geo3DUtil.decodeValueCeil(cell.yMaxEnc)+
" minz="+Geo3DUtil.decodeValueFloor(cell.zMinEnc)+" maxz="+Geo3DUtil.decodeValueCeil(cell.zMaxEnc));
log.println(" minx="+Geo3DUtil.decodeValueFloor(cell.xMinEnc, planetModel)+" maxx="+Geo3DUtil.decodeValueCeil(cell.xMaxEnc, planetModel)+
" miny="+Geo3DUtil.decodeValueFloor(cell.yMinEnc, planetModel)+" maxy="+Geo3DUtil.decodeValueCeil(cell.yMaxEnc, planetModel)+
" minz="+Geo3DUtil.decodeValueFloor(cell.zMinEnc, planetModel)+" maxz="+Geo3DUtil.decodeValueCeil(cell.zMaxEnc, planetModel));
}
switch (xyzSolid.getRelationship(shape)) {
@ -364,11 +374,13 @@ public class TestGeo3DPoint extends LuceneTestCase {
cell.xMinEnc, splitValue,
cell.yMinEnc, cell.yMaxEnc,
cell.zMinEnc, cell.zMaxEnc,
planetModel,
cell.splitCount+1);
Cell cell2 = new Cell(cell,
splitValue, cell.xMaxEnc,
cell.yMinEnc, cell.yMaxEnc,
cell.zMinEnc, cell.zMaxEnc,
planetModel,
cell.splitCount+1);
if (VERBOSE) {
log.println(" split cell1: " + cell1);
@ -390,11 +402,13 @@ public class TestGeo3DPoint extends LuceneTestCase {
cell.xMinEnc, cell.xMaxEnc,
cell.yMinEnc, splitValue,
cell.zMinEnc, cell.zMaxEnc,
planetModel,
cell.splitCount+1);
Cell cell2 = new Cell(cell,
cell.xMinEnc, cell.xMaxEnc,
splitValue, cell.yMaxEnc,
cell.zMinEnc, cell.zMaxEnc,
planetModel,
cell.splitCount+1);
if (VERBOSE) {
log.println(" split cell1: " + cell1);
@ -416,11 +430,13 @@ public class TestGeo3DPoint extends LuceneTestCase {
cell.xMinEnc, cell.xMaxEnc,
cell.yMinEnc, cell.yMaxEnc,
cell.zMinEnc, splitValue,
planetModel,
cell.splitCount+1);
Cell cell2 = new Cell(cell,
cell.xMinEnc, cell.xMaxEnc,
cell.yMinEnc, cell.yMaxEnc,
splitValue, cell.zMaxEnc,
planetModel,
cell.splitCount+1);
if (VERBOSE) {
log.println(" split cell1: " + cell1);
@ -543,7 +559,7 @@ public class TestGeo3DPoint extends LuceneTestCase {
}
}
verify(lats, lons);
verify(lats, lons, randomPlanetModel());
}
public void testPolygonOrdering() {
@ -551,14 +567,12 @@ public class TestGeo3DPoint extends LuceneTestCase {
51.204382859999996, 50.89947531437482, 50.8093624806861,50.8093624806861, 50.89947531437482, 51.204382859999996, 51.51015366140113, 51.59953838204167, 51.59953838204167, 51.51015366140113, 51.204382859999996};
final double[] lons = new double[] {
0.8747711978759765, 0.6509219832137298, 0.35960265165247807, 0.10290284834752167, -0.18841648321373008, -0.41226569787597667, -0.18960465285650027, 0.10285893781346236, 0.35964656218653757, 0.6521101528565002, 0.8747711978759765};
final Query q = Geo3DPoint.newPolygonQuery("point", new Polygon(lats, lons));
final Query q = Geo3DPoint.newPolygonQuery("point", randomPlanetModel(), new Polygon(lats, lons));
//System.out.println(q);
assertTrue(!q.toString().contains("GeoConcavePolygon"));
}
private static final double MEAN_EARTH_RADIUS_METERS = PlanetModel.WGS84_MEAN;
private static Query random3DQuery(final String field) {
private static Query random3DQuery(final String field, final PlanetModel planetModel) {
while (true) {
final int shapeType = random().nextInt(5);
switch (shapeType) {
@ -566,8 +580,8 @@ public class TestGeo3DPoint extends LuceneTestCase {
// Large polygons
final boolean isClockwise = random().nextDouble() < 0.5;
try {
final Query q = Geo3DPoint.newLargePolygonQuery(field, makePoly(PlanetModel.WGS84,
new GeoPoint(PlanetModel.WGS84, toRadians(GeoTestUtil.nextLatitude()), toRadians(GeoTestUtil.nextLongitude())),
final Query q = Geo3DPoint.newLargePolygonQuery(field, planetModel, makePoly(planetModel,
new GeoPoint(planetModel, toRadians(GeoTestUtil.nextLatitude()), toRadians(GeoTestUtil.nextLongitude())),
isClockwise,
true));
//System.err.println("Generated: "+q);
@ -582,8 +596,8 @@ public class TestGeo3DPoint extends LuceneTestCase {
// Polygons
final boolean isClockwise = random().nextDouble() < 0.5;
try {
final Query q = Geo3DPoint.newPolygonQuery(field, makePoly(PlanetModel.WGS84,
new GeoPoint(PlanetModel.WGS84, toRadians(GeoTestUtil.nextLatitude()), toRadians(GeoTestUtil.nextLongitude())),
final Query q = Geo3DPoint.newPolygonQuery(field, planetModel, makePoly(planetModel,
new GeoPoint(planetModel, toRadians(GeoTestUtil.nextLatitude()), toRadians(GeoTestUtil.nextLongitude())),
isClockwise,
true));
//System.err.println("Generated: "+q);
@ -596,9 +610,9 @@ public class TestGeo3DPoint extends LuceneTestCase {
case 1: {
// Circles
final double widthMeters = random().nextDouble() * Math.PI * MEAN_EARTH_RADIUS_METERS;
final double widthMeters = random().nextDouble() * Math.PI * planetModel.getMeanRadiusMeters() /*MEAN_EARTH_RADIUS_METERS*/;
try {
return Geo3DPoint.newDistanceQuery(field, GeoTestUtil.nextLatitude(), GeoTestUtil.nextLongitude(), widthMeters);
return Geo3DPoint.newDistanceQuery(field, planetModel, GeoTestUtil.nextLatitude(), GeoTestUtil.nextLongitude(), widthMeters);
} catch (IllegalArgumentException e) {
continue;
}
@ -608,7 +622,7 @@ public class TestGeo3DPoint extends LuceneTestCase {
// Rectangles
final Rectangle r = GeoTestUtil.nextBox();
try {
return Geo3DPoint.newBoxQuery(field, r.minLat, r.maxLat, r.minLon, r.maxLon);
return Geo3DPoint.newBoxQuery(field, planetModel, r.minLat, r.maxLat, r.minLon, r.maxLon);
} catch (IllegalArgumentException e) {
continue;
}
@ -618,7 +632,7 @@ public class TestGeo3DPoint extends LuceneTestCase {
// Paths
// TBD: Need to rework generation to be realistic
final int pointCount = random().nextInt(5) + 1;
final double width = random().nextDouble() * Math.PI * 0.5 * MEAN_EARTH_RADIUS_METERS;
final double width = random().nextDouble() * Math.PI * 0.5 * planetModel.getMeanRadiusMeters(); /* MEAN_EARTH_RADIUS_METERS;*/
final double[] latitudes = new double[pointCount];
final double[] longitudes = new double[pointCount];
for (int i = 0; i < pointCount; i++) {
@ -626,7 +640,7 @@ public class TestGeo3DPoint extends LuceneTestCase {
longitudes[i] = GeoTestUtil.nextLongitude();
}
try {
return Geo3DPoint.newPathQuery(field, latitudes, longitudes, width);
return Geo3DPoint.newPathQuery(field, latitudes, longitudes, width, planetModel);
} catch (IllegalArgumentException e) {
// This is what happens when we create a shape that is invalid. Although it is conceivable that there are cases where
// the exception is thrown incorrectly, we aren't going to be able to do that in this random test.
@ -642,7 +656,7 @@ public class TestGeo3DPoint extends LuceneTestCase {
}
// Poached from Geo3dRptTest.randomShape:
private static GeoShape randomShape() {
private static GeoShape randomShape(final PlanetModel planetModel) {
while (true) {
final int shapeType = random().nextInt(4);
switch (shapeType) {
@ -651,11 +665,11 @@ public class TestGeo3DPoint extends LuceneTestCase {
final int vertexCount = random().nextInt(3) + 3;
final List<GeoPoint> geoPoints = new ArrayList<>();
while (geoPoints.size() < vertexCount) {
final GeoPoint gPt = new GeoPoint(PlanetModel.WGS84, toRadians(GeoTestUtil.nextLatitude()), toRadians(GeoTestUtil.nextLongitude()));
final GeoPoint gPt = new GeoPoint(planetModel, toRadians(GeoTestUtil.nextLatitude()), toRadians(GeoTestUtil.nextLongitude()));
geoPoints.add(gPt);
}
try {
final GeoShape rval = GeoPolygonFactory.makeGeoPolygon(PlanetModel.WGS84, geoPoints);
final GeoShape rval = GeoPolygonFactory.makeGeoPolygon(planetModel, geoPoints);
if (rval == null) {
// Degenerate polygon
continue;
@ -677,7 +691,7 @@ public class TestGeo3DPoint extends LuceneTestCase {
double angle = random().nextDouble() * Math.PI/2.0;
try {
return GeoCircleFactory.makeGeoCircle(PlanetModel.WGS84, lat, lon, angle);
return GeoCircleFactory.makeGeoCircle(planetModel, lat, lon, angle);
} catch (IllegalArgumentException iae) {
// angle is too small; try again:
continue;
@ -701,7 +715,7 @@ public class TestGeo3DPoint extends LuceneTestCase {
lon1 = x;
}
return GeoBBoxFactory.makeGeoBBox(PlanetModel.WGS84, lat1, lat0, lon0, lon1);
return GeoBBoxFactory.makeGeoBBox(planetModel, lat1, lat0, lon0, lon1);
}
case 3: {
@ -710,10 +724,10 @@ public class TestGeo3DPoint extends LuceneTestCase {
final double width = toRadians(random().nextInt(89)+1);
final GeoPoint[] points = new GeoPoint[pointCount];
for (int i = 0; i < pointCount; i++) {
points[i] = new GeoPoint(PlanetModel.WGS84, toRadians(GeoTestUtil.nextLatitude()), toRadians(GeoTestUtil.nextLongitude()));
points[i] = new GeoPoint(planetModel, toRadians(GeoTestUtil.nextLatitude()), toRadians(GeoTestUtil.nextLongitude()));
}
try {
return GeoPathFactory.makeGeoPath(PlanetModel.WGS84, width, points);
return GeoPathFactory.makeGeoPath(planetModel, width, points);
} catch (IllegalArgumentException e) {
// This is what happens when we create a shape that is invalid. Although it is conceivable that there are cases where
// the exception is thrown incorrectly, we aren't going to be able to do that in this random test.
@ -727,7 +741,7 @@ public class TestGeo3DPoint extends LuceneTestCase {
}
}
private static void verify(double[] lats, double[] lons) throws Exception {
private static void verify(double[] lats, double[] lons, final PlanetModel planetModel) throws Exception {
IndexWriterConfig iwc = newIndexWriterConfig();
GeoPoint[] points = new GeoPoint[lats.length];
@ -737,8 +751,8 @@ public class TestGeo3DPoint extends LuceneTestCase {
for(int i=0;i<lats.length;i++) {
if (Double.isNaN(lats[i]) == false) {
//System.out.println("lats[" + i + "] = " + lats[i]);
unquantizedPoints[i] = new GeoPoint(PlanetModel.WGS84, toRadians(lats[i]), toRadians(lons[i]));
points[i] = quantize(unquantizedPoints[i]);
unquantizedPoints[i] = new GeoPoint(planetModel, toRadians(lats[i]), toRadians(lons[i]));
points[i] = quantize(unquantizedPoints[i], planetModel);
}
}
@ -763,7 +777,7 @@ public class TestGeo3DPoint extends LuceneTestCase {
doc.add(new NumericDocValuesField("id", id));
GeoPoint point = points[id];
if (point != null) {
doc.add(new Geo3DPoint("point", point.x, point.y, point.z));
doc.add(new Geo3DPoint("point", planetModel, point.x, point.y, point.z));
}
w.addDocument(doc);
if (id > 0 && random().nextInt(100) == 42) {
@ -799,7 +813,7 @@ public class TestGeo3DPoint extends LuceneTestCase {
}
*/
Query query = random3DQuery("point"); // Geo3DPoint.newShapeQuery("point", shape);
Query query = random3DQuery("point", planetModel); // Geo3DPoint.newShapeQuery("point", shape);
if (VERBOSE) {
System.err.println(" using query: " + query);
@ -842,7 +856,7 @@ public class TestGeo3DPoint extends LuceneTestCase {
GeoShape shape = ((PointInGeo3DShapeQuery)query).getShape();
XYZBounds bounds = new XYZBounds();
shape.getBounds(bounds);
XYZSolid solid = XYZSolidFactory.makeXYZSolid(PlanetModel.WGS84, bounds.getMinimumX(), bounds.getMaximumX(), bounds.getMinimumY(), bounds.getMaximumY(), bounds.getMinimumZ(), bounds.getMaximumZ());
XYZSolid solid = XYZSolidFactory.makeXYZSolid(planetModel, bounds.getMinimumX(), bounds.getMaximumX(), bounds.getMinimumY(), bounds.getMaximumY(), bounds.getMinimumZ(), bounds.getMaximumZ());
boolean expected = ((deleted.contains(id) == false) && shape.isWithin(point));
if (hits.get(docID) != expected) {
@ -855,9 +869,9 @@ public class TestGeo3DPoint extends LuceneTestCase {
b.append(" shape=" + shape + "\n");
b.append(" bounds=" + bounds + "\n");
b.append(" world bounds=(" +
" minX=" + PlanetModel.WGS84.getMinimumXValue() + " maxX=" + PlanetModel.WGS84.getMaximumXValue() +
" minY=" + PlanetModel.WGS84.getMinimumYValue() + " maxY=" + PlanetModel.WGS84.getMaximumYValue() +
" minZ=" + PlanetModel.WGS84.getMinimumZValue() + " maxZ=" + PlanetModel.WGS84.getMaximumZValue() + "\n");
" minX=" + planetModel.getMinimumXValue() + " maxX=" + planetModel.getMaximumXValue() +
" minY=" + planetModel.getMinimumYValue() + " maxY=" + planetModel.getMaximumYValue() +
" minZ=" + planetModel.getMinimumZValue() + " maxZ=" + planetModel.getMaximumZValue() + "\n");
b.append(" quantized point=" + point + " within shape? "+shape.isWithin(point)+" within bounds? "+solid.isWithin(point)+"\n");
b.append(" unquantized point=" + unquantizedPoint + " within shape? "+shape.isWithin(unquantizedPoint)+" within bounds? "+solid.isWithin(unquantizedPoint)+"\n");
b.append(" docID=" + docID + " deleted?=" + deleted.contains(id) + "\n");
@ -876,7 +890,7 @@ public class TestGeo3DPoint extends LuceneTestCase {
public void testToString() {
// Don't compare entire strings because Java 9 and Java 8 have slightly different values
Geo3DPoint point = new Geo3DPoint("point", 44.244272, 7.769736);
Geo3DPoint point = new Geo3DPoint("point", randomPlanetModel(), 44.244272, 7.769736);
final String stringToCompare = "Geo3DPoint <point: x=";
assertEquals(stringToCompare, point.toString().substring(0,stringToCompare.length()));
}
@ -893,7 +907,8 @@ public class TestGeo3DPoint extends LuceneTestCase {
}
public void testEquals() {
GeoShape shape = randomShape();
PlanetModel planetModel = randomPlanetModel();
GeoShape shape = randomShape(planetModel);
Query q = Geo3DPoint.newShapeQuery("point", shape);
assertEquals(q, Geo3DPoint.newShapeQuery("point", shape));
assertFalse(q.equals(Geo3DPoint.newShapeQuery("point2", shape)));
@ -901,14 +916,14 @@ public class TestGeo3DPoint extends LuceneTestCase {
// make a different random shape:
GeoShape shape2;
do {
shape2 = randomShape();
shape2 = randomShape(planetModel);
} while (shape.equals(shape2));
assertFalse(q.equals(Geo3DPoint.newShapeQuery("point", shape2)));
}
public void testComplexPolygons() {
final PlanetModel pm = PlanetModel.WGS84;
final PlanetModel pm = randomPlanetModel();
// Pick a random pole
final GeoPoint randomPole = new GeoPoint(pm, toRadians(GeoTestUtil.nextLatitude()), toRadians(GeoTestUtil.nextLongitude()));
int iters = atLeast(100);
@ -1137,33 +1152,33 @@ public class TestGeo3DPoint extends LuceneTestCase {
}
public void testEncodeDecodeCeil() throws Exception {
PlanetModel planetModel = randomPlanetModel();
// just for testing quantization error
final double ENCODING_TOLERANCE = Geo3DUtil.DECODE;
final double ENCODING_TOLERANCE = planetModel.DECODE;
int iters = atLeast(10000);
for(int iter=0;iter<iters;iter++) {
GeoPoint point = new GeoPoint(PlanetModel.WGS84, toRadians(GeoTestUtil.nextLatitude()), toRadians(GeoTestUtil.nextLongitude()));
double xEnc = Geo3DUtil.decodeValue(Geo3DUtil.encodeValue(point.x));
GeoPoint point = new GeoPoint(planetModel, toRadians(GeoTestUtil.nextLatitude()), toRadians(GeoTestUtil.nextLongitude()));
double xEnc = planetModel.decodeValue(planetModel.encodeValue(point.x));
assertEquals("x=" + point.x + " xEnc=" + xEnc + " diff=" + (point.x - xEnc), point.x, xEnc, ENCODING_TOLERANCE);
double yEnc = Geo3DUtil.decodeValue(Geo3DUtil.encodeValue(point.y));
double yEnc = planetModel.decodeValue(planetModel.encodeValue(point.y));
assertEquals("y=" + point.y + " yEnc=" + yEnc + " diff=" + (point.y - yEnc), point.y, yEnc, ENCODING_TOLERANCE);
double zEnc = Geo3DUtil.decodeValue(Geo3DUtil.encodeValue(point.z));
double zEnc = planetModel.decodeValue(planetModel.encodeValue(point.z));
assertEquals("z=" + point.z + " zEnc=" + zEnc + " diff=" + (point.z - zEnc), point.z, zEnc, ENCODING_TOLERANCE);
}
// check edge/interesting cases explicitly
double planetMax = PlanetModel.WGS84.getMaximumMagnitude();
double planetMax = planetModel.getMaximumMagnitude();
for (double value : new double[] {0.0, -planetMax, planetMax}) {
assertEquals(value, Geo3DUtil.decodeValue(Geo3DUtil.encodeValue(value)), ENCODING_TOLERANCE);
assertEquals(value, planetModel.decodeValue(planetModel.encodeValue(value)), ENCODING_TOLERANCE);
}
}
/** make sure values always go down: this is important for edge case consistency */
public void testEncodeDecodeRoundsDown() throws Exception {
PlanetModel planetModel = randomPlanetModel();
int iters = atLeast(1000);
for(int iter=0;iter<iters;iter++) {
final double latBase = GeoTestUtil.nextLatitude();
@ -1175,10 +1190,10 @@ public class TestGeo3DPoint extends LuceneTestCase {
for (int i = 0; i < 1000; i++) {
lat = Math.min(90, Math.nextUp(lat));
lon = Math.min(180, Math.nextUp(lon));
GeoPoint point = new GeoPoint(PlanetModel.WGS84, toRadians(lat), toRadians(lon));
GeoPoint pointEnc = new GeoPoint(Geo3DUtil.decodeValueFloor(Geo3DUtil.encodeValue(point.x)),
Geo3DUtil.decodeValueFloor(Geo3DUtil.encodeValue(point.y)),
Geo3DUtil.decodeValueFloor(Geo3DUtil.encodeValue(point.z)));
GeoPoint point = new GeoPoint(planetModel, toRadians(lat), toRadians(lon));
GeoPoint pointEnc = new GeoPoint(Geo3DUtil.decodeValueFloor(planetModel.encodeValue(point.x), planetModel),
Geo3DUtil.decodeValueFloor(planetModel.encodeValue(point.y), planetModel),
Geo3DUtil.decodeValueFloor(planetModel.encodeValue(point.z), planetModel));
assertTrue(pointEnc.x <= point.x);
assertTrue(pointEnc.y <= point.y);
assertTrue(pointEnc.z <= point.z);
@ -1190,10 +1205,10 @@ public class TestGeo3DPoint extends LuceneTestCase {
for (int i = 0; i < 1000; i++) {
lat = Math.max(-90, Math.nextDown(lat));
lon = Math.max(-180, Math.nextDown(lon));
GeoPoint point = new GeoPoint(PlanetModel.WGS84, toRadians(lat), toRadians(lon));
GeoPoint pointEnc = new GeoPoint(Geo3DUtil.decodeValueFloor(Geo3DUtil.encodeValue(point.x)),
Geo3DUtil.decodeValueFloor(Geo3DUtil.encodeValue(point.y)),
Geo3DUtil.decodeValueFloor(Geo3DUtil.encodeValue(point.z)));
GeoPoint point = new GeoPoint(planetModel, toRadians(lat), toRadians(lon));
GeoPoint pointEnc = new GeoPoint(Geo3DUtil.decodeValueFloor(planetModel.encodeValue(point.x), planetModel),
Geo3DUtil.decodeValueFloor(planetModel.encodeValue(point.y), planetModel),
Geo3DUtil.decodeValueFloor(planetModel.encodeValue(point.z), planetModel));
assertTrue(pointEnc.x <= point.x);
assertTrue(pointEnc.y <= point.y);
assertTrue(pointEnc.z <= point.z);
@ -1202,37 +1217,39 @@ public class TestGeo3DPoint extends LuceneTestCase {
}
public void testMinValueQuantization(){
int encoded = Geo3DUtil.MIN_ENCODED_VALUE;
double minValue= -PlanetModel.WGS84.getMaximumMagnitude();
PlanetModel planetModel = randomPlanetModel();
int encoded = planetModel.MIN_ENCODED_VALUE;
double minValue= -planetModel.getMaximumMagnitude();
//Normal encoding
double decoded = Geo3DUtil.decodeValue(encoded);
double decoded = planetModel.decodeValue(encoded);
assertEquals(minValue, decoded, 0d);
assertEquals(encoded, Geo3DUtil.encodeValue(decoded));
assertEquals(encoded, planetModel.encodeValue(decoded));
//Encoding floor
double decodedFloor = Geo3DUtil.decodeValueFloor(encoded);
double decodedFloor = Geo3DUtil.decodeValueFloor(encoded, planetModel);
assertEquals(minValue, decodedFloor, 0d);
assertEquals(encoded, Geo3DUtil.encodeValue(decodedFloor));
assertEquals(encoded, planetModel.encodeValue(decodedFloor));
//Encoding ceiling
double decodedCeiling = Geo3DUtil.decodeValueCeil(encoded);
double decodedCeiling = Geo3DUtil.decodeValueCeil(encoded, planetModel);
assertTrue(decodedCeiling > minValue);
assertEquals(encoded, Geo3DUtil.encodeValue(decodedCeiling));
assertEquals(encoded, planetModel.encodeValue(decodedCeiling));
}
public void testMaxValueQuantization(){
int encoded = Geo3DUtil.MAX_ENCODED_VALUE;
double maxValue= PlanetModel.WGS84.getMaximumMagnitude();
PlanetModel planetModel = randomPlanetModel();
int encoded = planetModel.MAX_ENCODED_VALUE;
double maxValue= planetModel.getMaximumMagnitude();
//Normal encoding
double decoded = Geo3DUtil.decodeValue(encoded);
double decoded = planetModel.decodeValue(encoded);
assertEquals(maxValue, decoded, 0d);
assertEquals(encoded, Geo3DUtil.encodeValue(decoded));
assertEquals(encoded, planetModel.encodeValue(decoded));
//Encoding floor
double decodedFloor = Geo3DUtil.decodeValueFloor(encoded);
double decodedFloor = Geo3DUtil.decodeValueFloor(encoded, planetModel);
assertTrue(decodedFloor < maxValue);
assertEquals(encoded, Geo3DUtil.encodeValue(decodedFloor));
assertEquals(encoded, planetModel.encodeValue(decodedFloor));
//Encoding ceiling
double decodedCeiling = Geo3DUtil.decodeValueCeil(encoded);
double decodedCeiling = Geo3DUtil.decodeValueCeil(encoded, planetModel);
assertEquals(maxValue, decodedCeiling, 0d);
assertEquals(encoded, Geo3DUtil.encodeValue(decodedCeiling));
assertEquals(encoded, planetModel.encodeValue(decodedCeiling));
}
// poached from TestGeoEncodingUtils.testLatitudeQuantization:
@ -1244,33 +1261,34 @@ public class TestGeo3DPoint extends LuceneTestCase {
*/
public void testQuantization() throws Exception {
Random random = random();
PlanetModel planetModel = randomPlanetModel();
for (int i = 0; i < 10000; i++) {
int encoded = random.nextInt();
if (encoded <= Geo3DUtil.MIN_ENCODED_VALUE) {
if (encoded <= planetModel.MIN_ENCODED_VALUE) {
continue;
}
if (encoded >= Geo3DUtil.MAX_ENCODED_VALUE) {
if (encoded >= planetModel.MAX_ENCODED_VALUE) {
continue;
}
double min = encoded * Geo3DUtil.DECODE;
double decoded = Geo3DUtil.decodeValueFloor(encoded);
double min = encoded * planetModel.DECODE;
double decoded = Geo3DUtil.decodeValueFloor(encoded, planetModel);
// should exactly equal expected value
assertEquals(min, decoded, 0.0D);
// should round-trip
assertEquals(encoded, Geo3DUtil.encodeValue(decoded));
assertEquals(encoded, planetModel.encodeValue(decoded));
// test within the range
if (encoded != Integer.MAX_VALUE) {
// this is the next representable value
// all double values between [min .. max) should encode to the current integer
double max = min + Geo3DUtil.DECODE;
assertEquals(max, Geo3DUtil.decodeValueFloor(encoded+1), 0.0D);
assertEquals(encoded+1, Geo3DUtil.encodeValue(max));
double max = min + planetModel.DECODE;
assertEquals(max, Geo3DUtil.decodeValueFloor(encoded+1, planetModel), 0.0D);
assertEquals(encoded+1, planetModel.encodeValue(max));
// first and last doubles in range that will be quantized
double minEdge = Math.nextUp(min);
double maxEdge = Math.nextDown(max);
assertEquals(encoded, Geo3DUtil.encodeValue(minEdge));
assertEquals(encoded, Geo3DUtil.encodeValue(maxEdge));
assertEquals(encoded, planetModel.encodeValue(minEdge));
assertEquals(encoded, planetModel.encodeValue(maxEdge));
// check random values within the double range
long minBits = NumericUtils.doubleToSortableLong(minEdge);
@ -1278,30 +1296,30 @@ public class TestGeo3DPoint extends LuceneTestCase {
for (int j = 0; j < 100; j++) {
double value = NumericUtils.sortableLongToDouble(TestUtil.nextLong(random, minBits, maxBits));
// round down
assertEquals(encoded, Geo3DUtil.encodeValue(value));
assertEquals(encoded, planetModel.encodeValue(value));
}
}
}
}
public void testEncodeDecodeIsStable() throws Exception {
PlanetModel planetModel = randomPlanetModel();
int iters = atLeast(1000);
for(int iter=0;iter<iters;iter++) {
double lat = GeoTestUtil.nextLatitude();
double lon = GeoTestUtil.nextLongitude();
GeoPoint point = new GeoPoint(PlanetModel.WGS84, toRadians(lat), toRadians(lon));
GeoPoint point = new GeoPoint(planetModel, toRadians(lat), toRadians(lon));
// encode point
GeoPoint pointEnc = new GeoPoint(Geo3DUtil.decodeValue(Geo3DUtil.encodeValue(point.x)),
Geo3DUtil.decodeValue(Geo3DUtil.encodeValue(point.y)),
Geo3DUtil.decodeValue(Geo3DUtil.encodeValue(point.z)));
GeoPoint pointEnc = new GeoPoint(planetModel.decodeValue(planetModel.encodeValue(point.x)),
planetModel.decodeValue(planetModel.encodeValue(point.y)),
planetModel.decodeValue(planetModel.encodeValue(point.z)));
// encode it again (double encode)
GeoPoint pointEnc2 = new GeoPoint(Geo3DUtil.decodeValue(Geo3DUtil.encodeValue(pointEnc.x)),
Geo3DUtil.decodeValue(Geo3DUtil.encodeValue(pointEnc.y)),
Geo3DUtil.decodeValue(Geo3DUtil.encodeValue(pointEnc.z)));
GeoPoint pointEnc2 = new GeoPoint(planetModel.decodeValue(planetModel.encodeValue(pointEnc.x)),
planetModel.decodeValue(planetModel.encodeValue(pointEnc.y)),
planetModel.decodeValue(planetModel.encodeValue(pointEnc.z)));
//System.out.println("TEST " + iter + ":\n point =" + point + "\n pointEnc =" + pointEnc + "\n pointEnc2=" + pointEnc2);
assertEquals(pointEnc.x, pointEnc2.x, 0.0);
@ -1311,14 +1329,14 @@ public class TestGeo3DPoint extends LuceneTestCase {
}
/** Clips the incoming value to the allowed min/max range before encoding, instead of throwing an exception. */
private static int encodeValueLenient(double x) {
double planetMax = PlanetModel.WGS84.getMaximumMagnitude();
private static int encodeValueLenient(double x, PlanetModel planetModel) {
double planetMax = planetModel.getMaximumMagnitude();
if (x > planetMax) {
x = planetMax;
} else if (x < -planetMax) {
x = -planetMax;
}
return Geo3DUtil.encodeValue(x);
return planetModel.encodeValue(x);
}
private static class ExplainingVisitor implements IntersectVisitor {
@ -1381,9 +1399,9 @@ public class TestGeo3DPoint extends LuceneTestCase {
}
} else {
if (docID == targetDocID) {
double x = Geo3DPoint.decodeDimension(packedValue, 0);
double y = Geo3DPoint.decodeDimension(packedValue, Integer.BYTES);
double z = Geo3DPoint.decodeDimension(packedValue, 2 * Integer.BYTES);
double x = Geo3DPoint.decodeDimension(packedValue, 0, shape.getPlanetModel());
double y = Geo3DPoint.decodeDimension(packedValue, Integer.BYTES, shape.getPlanetModel());
double z = Geo3DPoint.decodeDimension(packedValue, 2 * Integer.BYTES, shape.getPlanetModel());
b.append("leaf visit docID=" + docID + " x=" + x + " y=" + y + " z=" + z + "\n");
in.visit(docID, packedValue);
}
@ -1448,13 +1466,13 @@ public class TestGeo3DPoint extends LuceneTestCase {
@Override
public String toString() {
double xMin = Geo3DUtil.decodeValueFloor(NumericUtils.sortableBytesToInt(minPackedValue, 0));
double xMax = Geo3DUtil.decodeValueCeil(NumericUtils.sortableBytesToInt(maxPackedValue, 0));
double yMin = Geo3DUtil.decodeValueFloor(NumericUtils.sortableBytesToInt(minPackedValue, 1 * Integer.BYTES));
double yMax = Geo3DUtil.decodeValueCeil(NumericUtils.sortableBytesToInt(maxPackedValue, 1 * Integer.BYTES));
double zMin = Geo3DUtil.decodeValueFloor(NumericUtils.sortableBytesToInt(minPackedValue, 2 * Integer.BYTES));
double zMax = Geo3DUtil.decodeValueCeil(NumericUtils.sortableBytesToInt(maxPackedValue, 2 * Integer.BYTES));
final XYZSolid xyzSolid = XYZSolidFactory.makeXYZSolid(PlanetModel.WGS84, xMin, xMax, yMin, yMax, zMin, zMax);
double xMin = Geo3DUtil.decodeValueFloor(NumericUtils.sortableBytesToInt(minPackedValue, 0), shape.getPlanetModel());
double xMax = Geo3DUtil.decodeValueCeil(NumericUtils.sortableBytesToInt(maxPackedValue, 0), shape.getPlanetModel());
double yMin = Geo3DUtil.decodeValueFloor(NumericUtils.sortableBytesToInt(minPackedValue, 1 * Integer.BYTES), shape.getPlanetModel());
double yMax = Geo3DUtil.decodeValueCeil(NumericUtils.sortableBytesToInt(maxPackedValue, 1 * Integer.BYTES), shape.getPlanetModel());
double zMin = Geo3DUtil.decodeValueFloor(NumericUtils.sortableBytesToInt(minPackedValue, 2 * Integer.BYTES), shape.getPlanetModel());
double zMax = Geo3DUtil.decodeValueCeil(NumericUtils.sortableBytesToInt(maxPackedValue, 2 * Integer.BYTES), shape.getPlanetModel());
final XYZSolid xyzSolid = XYZSolidFactory.makeXYZSolid(shape.getPlanetModel(), xMin, xMax, yMin, yMax, zMin, zMax);
final int relationship = xyzSolid.getRelationship(shape);
final boolean pointWithinCell = xyzSolid.isWithin(targetDocPoint);
final boolean scaledWithinCell = xyzSolid.isWithin(scaledDocPoint);

245
lucene/spatial3d/src/test/org/apache/lucene/spatial3d/geom/Geo3DUtil.java
View File

@ -16,89 +16,13 @@
*/
package org.apache.lucene.spatial3d.geom;
import org.apache.lucene.geo.Polygon;
import org.apache.lucene.geo.GeoUtils;
import java.util.List;
import java.util.ArrayList;
class Geo3DUtil {
/** How many radians are in one earth surface meter */
final static double RADIANS_PER_METER = 1.0 / PlanetModel.WGS84_MEAN;
/** How many radians are in one degree */
final static double RADIANS_PER_DEGREE = Math.PI / 180.0;
/** How many degrees in a radian */
final static double DEGREES_PER_RADIAN = 180.0 / Math.PI;
private static final double MAX_VALUE = PlanetModel.WGS84.getMaximumMagnitude();
private static final int BITS = 32;
private static final double MUL = (0x1L<<BITS)/(2*MAX_VALUE);
static final double DECODE = getNextSafeDouble(1/MUL);
static final int MIN_ENCODED_VALUE = encodeValue(-MAX_VALUE);
static final int MAX_ENCODED_VALUE = encodeValue(MAX_VALUE);
public static int encodeValue(double x) {
if (x > MAX_VALUE) {
throw new IllegalArgumentException("value=" + x + " is out-of-bounds (greater than WGS84's planetMax=" + MAX_VALUE + ")");
}
if (x < -MAX_VALUE) {
throw new IllegalArgumentException("value=" + x + " is out-of-bounds (less than than WGS84's -planetMax=" + -MAX_VALUE + ")");
}
long result = (long) Math.floor(x / DECODE);
assert result >= Integer.MIN_VALUE;
assert result <= Integer.MAX_VALUE;
return (int) result;
}
public static double decodeValue(int x) {
double result;
if (x == MIN_ENCODED_VALUE) {
// We must special case this, because -MAX_VALUE is not guaranteed to land precisely at a floor value, and we don't ever want to
// return a value outside of the planet's range (I think?). The max value is "safe" because we floor during encode:
result = -MAX_VALUE;
} else if (x == MAX_ENCODED_VALUE) {
result = MAX_VALUE;
} else {
// We decode to the center value; this keeps the encoding stable
result = (x+0.5) * DECODE;
}
assert result >= -MAX_VALUE && result <= MAX_VALUE;
return result;
}
/** Returns smallest double that would encode to int x. */
// NOTE: keep this package private!!
static double decodeValueFloor(int x) {
return x * DECODE;
}
/** Returns a double value >= x such that if you multiply that value by an int, and then
* divide it by that int again, you get precisely the same value back */
private static double getNextSafeDouble(double x) {
// Move to double space:
long bits = Double.doubleToLongBits(x);
// Make sure we are beyond the actual maximum value:
bits += Integer.MAX_VALUE;
// Clear the bottom 32 bits:
bits &= ~((long) Integer.MAX_VALUE);
// Convert back to double:
double result = Double.longBitsToDouble(bits);
assert result > x;
return result;
}
/** Returns largest double that would encode to int x. */
// NOTE: keep this package private!!
static double decodeValueCeil(int x) {
assert x < Integer.MAX_VALUE;
return Math.nextDown((x+1) * DECODE);
}
/** Converts degress to radians */
static double fromDegrees(final double degrees) {
return degrees * RADIANS_PER_DEGREE;
@ -108,173 +32,4 @@ class Geo3DUtil {
static double toDegrees(final double radians) {
return radians * DEGREES_PER_RADIAN;
}
/** Converts earth-surface meters to radians */
static double fromMeters(final double meters) {
return meters * RADIANS_PER_METER;
}
/**
* Convert a set of Polygon objects into a GeoPolygon.
* @param polygons are the Polygon objects.
* @return the GeoPolygon.
*/
static GeoPolygon fromPolygon(final Polygon... polygons) {
//System.err.println("Creating polygon...");
if (polygons.length < 1) {
throw new IllegalArgumentException("need at least one polygon");
}
final GeoPolygon shape;
if (polygons.length == 1) {
final GeoPolygon component = fromPolygon(polygons[0]);
if (component == null) {
// Polygon is degenerate
shape = new GeoCompositePolygon(PlanetModel.WGS84);
} else {
shape = component;
}
} else {
final GeoCompositePolygon poly = new GeoCompositePolygon(PlanetModel.WGS84);
for (final Polygon p : polygons) {
final GeoPolygon component = fromPolygon(p);
if (component != null) {
poly.addShape(component);
}
}
shape = poly;
}
return shape;
//System.err.println("...done");
}
/**
* Convert a Polygon object to a large GeoPolygon.
* @param polygons is the list of polygons to convert.
* @return the large GeoPolygon.
*/
static GeoPolygon fromLargePolygon(final Polygon... polygons) {
if (polygons.length < 1) {
throw new IllegalArgumentException("need at least one polygon");
}
return GeoPolygonFactory.makeLargeGeoPolygon(PlanetModel.WGS84, convertToDescription(polygons));
}
/**
* Convert input parameters to a path.
* @param pathLatitudes latitude values for points of the path: must be within standard +/-90 coordinate bounds.
* @param pathLongitudes longitude values for points of the path: must be within standard +/-180 coordinate bounds.
* @param pathWidthMeters width of the path in meters.
* @return the path.
*/
static GeoPath fromPath(final double[] pathLatitudes, final double[] pathLongitudes, final double pathWidthMeters) {
if (pathLatitudes.length != pathLongitudes.length) {
throw new IllegalArgumentException("same number of latitudes and longitudes required");
}
final GeoPoint[] points = new GeoPoint[pathLatitudes.length];
for (int i = 0; i < pathLatitudes.length; i++) {
final double latitude = pathLatitudes[i];
final double longitude = pathLongitudes[i];
GeoUtils.checkLatitude(latitude);
GeoUtils.checkLongitude(longitude);
points[i] = new GeoPoint(PlanetModel.WGS84, fromDegrees(latitude), fromDegrees(longitude));
}
return GeoPathFactory.makeGeoPath(PlanetModel.WGS84, fromMeters(pathWidthMeters), points);
}
/**
* Convert input parameters to a circle.
* @param latitude latitude at the center: must be within standard +/-90 coordinate bounds.
* @param longitude longitude at the center: must be within standard +/-180 coordinate bounds.
* @param radiusMeters maximum distance from the center in meters: must be non-negative and finite.
* @return the circle.
*/
static GeoCircle fromDistance(final double latitude, final double longitude, final double radiusMeters) {
GeoUtils.checkLatitude(latitude);
GeoUtils.checkLongitude(longitude);
return GeoCircleFactory.makeGeoCircle(PlanetModel.WGS84, fromDegrees(latitude), fromDegrees(longitude), fromMeters(radiusMeters));
}
/**
* Convert input parameters to a box.
* @param minLatitude latitude lower bound: must be within standard +/-90 coordinate bounds.
* @param maxLatitude latitude upper bound: must be within standard +/-90 coordinate bounds.
* @param minLongitude longitude lower bound: must be within standard +/-180 coordinate bounds.
* @param maxLongitude longitude upper bound: must be within standard +/-180 coordinate bounds.
* @return the box.
*/
static GeoBBox fromBox(final double minLatitude, final double maxLatitude, final double minLongitude, final double maxLongitude) {
GeoUtils.checkLatitude(minLatitude);
GeoUtils.checkLongitude(minLongitude);
GeoUtils.checkLatitude(maxLatitude);
GeoUtils.checkLongitude(maxLongitude);
return GeoBBoxFactory.makeGeoBBox(PlanetModel.WGS84,
Geo3DUtil.fromDegrees(maxLatitude), Geo3DUtil.fromDegrees(minLatitude), Geo3DUtil.fromDegrees(minLongitude), Geo3DUtil.fromDegrees(maxLongitude));
}
/**
* Convert a Polygon object into a GeoPolygon.
* This method uses
* @param polygon is the Polygon object.
* @return the GeoPolygon.
*/
private static GeoPolygon fromPolygon(final Polygon polygon) {
// First, assemble the "holes". The geo3d convention is to use the same polygon sense on the inner ring as the
// outer ring, so we process these recursively with reverseMe flipped.
final Polygon[] theHoles = polygon.getHoles();
final List<GeoPolygon> holeList = new ArrayList<>(theHoles.length);
for (final Polygon hole : theHoles) {
//System.out.println("Hole: "+hole);
final GeoPolygon component = fromPolygon(hole);
if (component != null) {
holeList.add(component);
}
}
// Now do the polygon itself
final double[] polyLats = polygon.getPolyLats();
final double[] polyLons = polygon.getPolyLons();
// I presume the arguments have already been checked
final List<GeoPoint> points = new ArrayList<>(polyLats.length-1);
// We skip the last point anyway because the API requires it to be repeated, and geo3d doesn't repeat it.
for (int i = 0; i < polyLats.length - 1; i++) {
final int index = polyLats.length - 2 - i;
points.add(new GeoPoint(PlanetModel.WGS84, fromDegrees(polyLats[index]), fromDegrees(polyLons[index])));
}
//System.err.println(" building polygon with "+points.size()+" points...");
final GeoPolygon rval = GeoPolygonFactory.makeGeoPolygon(PlanetModel.WGS84, points, holeList);
//System.err.println(" ...done");
return rval;
}
/**
* Convert a list of polygons to a list of polygon descriptions.
* @param polygons is the list of polygons to convert.
* @return the list of polygon descriptions.
*/
private static List<GeoPolygonFactory.PolygonDescription> convertToDescription(final Polygon... polygons) {
final List<GeoPolygonFactory.PolygonDescription> descriptions = new ArrayList<>(polygons.length);
for (final Polygon polygon : polygons) {
final Polygon[] theHoles = polygon.getHoles();
final List<GeoPolygonFactory.PolygonDescription> holes = convertToDescription(theHoles);
// Now do the polygon itself
final double[] polyLats = polygon.getPolyLats();
final double[] polyLons = polygon.getPolyLons();
// I presume the arguments have already been checked
final List<GeoPoint> points = new ArrayList<>(polyLats.length-1);
// We skip the last point anyway because the API requires it to be repeated, and geo3d doesn't repeat it.
for (int i = 0; i < polyLats.length - 1; i++) {
final int index = polyLats.length - 2 - i;
points.add(new GeoPoint(PlanetModel.WGS84, fromDegrees(polyLats[index]), fromDegrees(polyLons[index])));
}
descriptions.add(new GeoPolygonFactory.PolygonDescription(points, holes));
}
return descriptions;
}
}

6
lucene/spatial3d/src/test/org/apache/lucene/spatial3d/geom/GeoCircleTest.java
View File

@ -172,7 +172,7 @@ public class GeoCircleTest extends LuceneTestCase {
// Sixth BKD discovered failure
c = GeoCircleFactory.makeGeoCircle(PlanetModel.WGS84,-0.006450320645814321,0.004660694205115142,0.00489710732634323);
//xyzb = new XYZBounds();
//c.getBounds(xyzb);
//zScaling.getBounds(xyzb);
//System.err.println("xmin="+xyzb.getMinimumX()+", xmax="+xyzb.getMaximumX()+",ymin="+xyzb.getMinimumY()+", ymax="+xyzb.getMaximumY()+",zmin="+xyzb.getMinimumZ()+", zmax="+xyzb.getMaximumZ());
//xmin=1.0010356621420726, xmax=1.0011141249179447,ymin=-2.5326643901354566E-4, ymax=0.009584741915757169,zmin=-0.011359874956269283, zmax=-0.0015549504447452225
area = GeoAreaFactory.makeGeoArea(PlanetModel.WGS84,1.0010822580620098,1.0010945779732867,0.007079167343247293,0.007541006774427837,-0.0021855011220022575,-0.001896122718181518);
@ -183,9 +183,9 @@ public class GeoCircleTest extends LuceneTestCase {
//p2 = new GeoPoint(PlanetModel.WGS84,-0.002164069780096702, 0.007505617500830066);
p2 = new GeoPoint(PlanetModel.WGS84,p1.getLatitude(),p1.getLongitude());
assertTrue(PlanetModel.WGS84.pointOnSurface(p2));
assertTrue(!c.isWithin(p2));
assertTrue(!zScaling.isWithin(p2));
assertTrue(!area.isWithin(p2));
assertTrue(!c.isWithin(p1));
assertTrue(!zScaling.isWithin(p1));
assertTrue(PlanetModel.WGS84.pointOnSurface(p1)); // This fails
assertTrue(!area.isWithin(p1)); // This fails
*/

6
lucene/spatial3d/src/test/org/apache/lucene/spatial3d/geom/GeoExactCircleTest.java
View File

@ -55,7 +55,7 @@ public class GeoExactCircleTest extends RandomGeo3dShapeGenerator{
@Test
public void testSurfacePointOnBearingScale(){
PlanetModel p1 = PlanetModel.WGS84;
PlanetModel p2 = new PlanetModel(0.5 * PlanetModel.WGS84.ab, 0.5 * PlanetModel.WGS84.c );
PlanetModel p2 = new PlanetModel(0.5 * PlanetModel.WGS84.xyScaling, 0.5 * PlanetModel.WGS84.zScaling);
GeoPoint point1P1 = new GeoPoint(p1, 0, 0);
GeoPoint point2P1 = new GeoPoint(p1, 1, 1);
GeoPoint point1P2 = new GeoPoint(p2, point1P1.getLatitude(), point1P1.getLongitude());
@ -92,7 +92,7 @@ public class GeoExactCircleTest extends RandomGeo3dShapeGenerator{
@Repeat(iterations = 100)
public void RandomPointBearingCardinalTest(){
//surface distance calculations methods start not converging when
//planet flattening > 0.4
//planet scaledFlattening > 0.4
PlanetModel planetModel;
do {
double ab = random().nextDouble() * 2;
@ -102,7 +102,7 @@ public class GeoExactCircleTest extends RandomGeo3dShapeGenerator{
} else {
planetModel = new PlanetModel(c, ab);
}
} while (Math.abs(planetModel.flattening) > 0.4);
} while (Math.abs(planetModel.scaledFlattening) > 0.4);
GeoPoint center = randomGeoPoint(planetModel);
double radius = random().nextDouble() * 0.9 * planetModel.minimumPoleDistance;
checkBearingPoint(planetModel, center, radius, 0);

2
lucene/spatial3d/src/test/org/apache/lucene/spatial3d/geom/GeoPathTest.java
View File

@ -204,7 +204,7 @@ public class GeoPathTest extends LuceneTestCase {
assertTrue(relationship == GeoArea.WITHIN || relationship == GeoArea.OVERLAPS);
assertTrue(area.isWithin(point));
// No longer true due to fixed GeoStandardPath waypoints.
//assertTrue(c.isWithin(point));
//assertTrue(zScaling.isWithin(point));
c = new GeoStandardPath(PlanetModel.WGS84, 0.6894050545377601);
c.addPoint(-0.0788176065762948, 0.9431251741731624);

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

@ -122,7 +122,7 @@ public class GeoPolygonTest extends LuceneTestCase {
GeoPolygonFactory.PolygonDescription pd = new GeoPolygonFactory.PolygonDescription(points);
c = GeoPolygonFactory.makeGeoPolygon(PlanetModel.SPHERE, pd);
//System.out.println(c);
//System.out.println(zScaling);
// Middle point should NOT be within!!
gp = new GeoPoint(PlanetModel.SPHERE, 0.0, -0.5);
@ -143,7 +143,7 @@ public class GeoPolygonTest extends LuceneTestCase {
pd = new GeoPolygonFactory.PolygonDescription(points);
c = GeoPolygonFactory.makeGeoPolygon(PlanetModel.SPHERE, pd);
//System.out.println(c);
//System.out.println(zScaling);
// Middle point should be within!!
gp = new GeoPoint(PlanetModel.SPHERE, 0.0, -0.5);
@ -251,7 +251,7 @@ public class GeoPolygonTest extends LuceneTestCase {
shapes.add(pd);
c = GeoPolygonFactory.makeLargeGeoPolygon(PlanetModel.SPHERE, shapes);
//System.out.println("Large polygon = "+c);
//System.out.println("Large polygon = "+zScaling);
// Sample some points within
gp = new GeoPoint(PlanetModel.SPHERE, 0.0, -0.45);
@ -452,10 +452,10 @@ public class GeoPolygonTest extends LuceneTestCase {
@Test
public void testGeoPolygonBoundsCase2() {
// [junit4] 1> TEST: iter=23 shape=GeoCompositeMembershipShape: {[GeoConvexPolygon: {planetmodel=PlanetModel(ab=0.7563871189161702 c=1.2436128810838298), points=
// [junit4] 1> TEST: iter=23 shape=GeoCompositeMembershipShape: {[GeoConvexPolygon: {planetmodel=PlanetModel(xyScaling=0.7563871189161702 zScaling=1.2436128810838298), points=
// [[lat=0.014071770744627236, lon=0.011030818292803128],
// [lat=0.006772117088906782, lon=-0.0012531892445234592],
// [lat=0.0022201615609504792, lon=0.005941293187389326]]}, GeoConcavePolygon: {planetmodel=PlanetModel(ab=0.7563871189161702 c=1.2436128810838298), points=
// [lat=0.0022201615609504792, lon=0.005941293187389326]]}, GeoConcavePolygon: {planetmodel=PlanetModel(xyScaling=0.7563871189161702 zScaling=1.2436128810838298), points=
// [[lat=-0.005507100238396111, lon=-0.008487706131259667],
// [lat=0.014071770744627236, lon=0.011030818292803128],
// [lat=0.0022201615609504792, lon=0.005941293187389326]]}]}
@ -476,7 +476,7 @@ public class GeoPolygonTest extends LuceneTestCase {
BitSet p2bits = new BitSet();
p2bits.set(1, true);
c.addShape(new GeoConcavePolygon(pm, points2, p2bits, false));
//System.out.println(c);
//System.out.println(zScaling);
// [junit4] 1> point=[lat=0.003540694517552105, lon=-9.99517927901697E-4]
// [junit4] 1> quantized=[X=0.7563849869428783, Y=-7.560204674780763E-4, Z=0.0026781405884151086]
@ -511,10 +511,10 @@ doc=906 added here:
shape:
[junit4] 1> TEST: iter=18 shape=GeoCompositeMembershipShape: {[GeoConvexPolygon: {
planetmodel=PlanetModel(ab=0.8568069516722363 c=1.1431930483277637), points=
planetmodel=PlanetModel(xyScaling=0.8568069516722363 zScaling=1.1431930483277637), points=
[[lat=1.1577814487635816, lon=1.6283601832010004],
[lat=0.6664570999069251, lon=2.0855825542851574],
[lat=-0.23953537010974632, lon=1.8498724094352876]]}, GeoConcavePolygon: {planetmodel=PlanetModel(ab=0.8568069516722363 c=1.1431930483277637), points=
[lat=-0.23953537010974632, lon=1.8498724094352876]]}, GeoConcavePolygon: {planetmodel=PlanetModel(xyScaling=0.8568069516722363 zScaling=1.1431930483277637), points=
[[lat=1.1577814487635816, lon=1.6283601832010004],
[lat=-0.23953537010974632, lon=1.8498724094352876],
[lat=-1.1766904875978805, lon=-2.1346828411344436]]}]}
@ -535,14 +535,14 @@ shape:
BitSet p2bits = new BitSet();
p2bits.set(1, true);
c.addShape(new GeoConcavePolygon(pm, points2, p2bits, false));
//System.out.println(c);
//System.out.println(zScaling);
GeoPoint point = new GeoPoint(pm, -0.9825762558001477, 2.4832136904725273);
GeoPoint quantizedPoint = new GeoPoint(-0.4505446160475436, 0.34850109186970535, -0.8539966368663765);
GeoArea xyzSolid = GeoAreaFactory.makeGeoArea(pm,
-0.6107484000858642, -0.39518364125756916, -0.8568069517709872, 0.8568069517709872, -1.1431930485939341, 1.1431930485939341);
//System.out.println("relationship = "+xyzSolid.getRelationship(c));
//System.out.println("relationship = "+xyzSolid.getRelationship(zScaling));
assertTrue(xyzSolid.getRelationship(c) == GeoArea.OVERLAPS);
}
@ -1272,7 +1272,7 @@ shape:
final GeoPoint negativeX = new GeoPoint(PlanetModel.WGS84, 0.0, Math.PI);
final GeoPoint negativeY = new GeoPoint(PlanetModel.WGS84, 0.0, -Math.PI * 0.5);
final GeoPoint positiveY = new GeoPoint(PlanetModel.WGS84, 0.0, Math.PI * 0.5);
final GeoPoint testPoint = new GeoPoint(-0.07416172733314662, 0.5686488061136892, 0.8178445379402641);
final GeoPoint testPoint = new GeoPoint(-0.074161727332972, 0.5686488061123504, 0.8178445379383386);
// Construct a standard polygon first to see what that does. This winds up being a large polygon under the covers.
GeoPolygon standard = GeoPolygonFactory.makeGeoPolygon(PlanetModel.WGS84, pd);