com.google.common.geometry.S2Cell Maven / Gradle / Ivy
The newest version!
/*
* Copyright 2005 Google Inc.
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package com.google.common.geometry;
/**
* An S2Cell is an S2Region object that represents a cell. Unlike S2CellIds, it
* supports efficient containment and intersection tests. However, it is also a
* more expensive representation.
*
*/
public final strictfp class S2Cell implements S2Region {
private static final int MAX_CELL_SIZE = 1 << S2CellId.MAX_LEVEL;
byte face;
byte level;
byte orientation;
S2CellId cellId;
double[][] uv = new double[2][2];
/**
* Default constructor used only internally.
*/
S2Cell() {
}
/**
* An S2Cell always corresponds to a particular S2CellId. The other
* constructors are just convenience methods.
*/
public S2Cell(S2CellId id) {
init(id);
}
// This is a static method in order to provide named parameters.
public static S2Cell fromFacePosLevel(int face, byte pos, int level) {
return new S2Cell(S2CellId.fromFacePosLevel(face, pos, level));
}
// Convenience methods.
public S2Cell(S2Point p) {
init(S2CellId.fromPoint(p));
}
public S2Cell(S2LatLng ll) {
init(S2CellId.fromLatLng(ll));
}
public S2CellId id() {
return cellId;
}
public int face() {
return face;
}
public byte level() {
return level;
}
public byte orientation() {
return orientation;
}
public boolean isLeaf() {
return level == S2CellId.MAX_LEVEL;
}
public S2Point getVertex(int k) {
return S2Point.normalize(getVertexRaw(k));
}
/**
* Return the k-th vertex of the cell (k = 0,1,2,3). Vertices are returned in
* CCW order. The points returned by GetVertexRaw are not necessarily unit
* length.
*/
public S2Point getVertexRaw(int k) {
// Vertices are returned in the order SW, SE, NE, NW.
return S2Projections.faceUvToXyz(face, uv[0][(k >> 1) ^ (k & 1)], uv[1][k >> 1]);
}
public S2Point getEdge(int k) {
return S2Point.normalize(getEdgeRaw(k));
}
public S2Point getEdgeRaw(int k) {
switch (k) {
case 0:
return S2Projections.getVNorm(face, uv[1][0]); // South
case 1:
return S2Projections.getUNorm(face, uv[0][1]); // East
case 2:
return S2Point.neg(S2Projections.getVNorm(face, uv[1][1])); // North
default:
return S2Point.neg(S2Projections.getUNorm(face, uv[0][0])); // West
}
}
/**
* Return the inward-facing normal of the great circle passing through the
* edge from vertex k to vertex k+1 (mod 4). The normals returned by
* GetEdgeRaw are not necessarily unit length.
*
* If this is not a leaf cell, set children[0..3] to the four children of
* this cell (in traversal order) and return true. Otherwise returns false.
* This method is equivalent to the following:
*
* for (pos=0, id=child_begin(); id != child_end(); id = id.next(), ++pos)
* children[i] = S2Cell(id);
*
* except that it is more than two times faster.
*/
public boolean subdivide(S2Cell children[]) {
// This function is equivalent to just iterating over the child cell ids
// and calling the S2Cell constructor, but it is about 2.5 times faster.
if (cellId.isLeaf()) {
return false;
}
// Compute the cell midpoint in uv-space.
R2Vector uvMid = getCenterUV();
// Create four children with the appropriate bounds.
S2CellId id = cellId.childBegin();
for (int pos = 0; pos < 4; ++pos, id = id.next()) {
S2Cell child = children[pos];
child.face = face;
child.level = (byte) (level + 1);
child.orientation = (byte) (orientation ^ S2.posToOrientation(pos));
child.cellId = id;
int ij = S2.posToIJ(orientation, pos);
for (int d = 0; d < 2; ++d) {
// The dimension 0 index (i/u) is in bit 1 of ij.
int m = 1 - ((ij >> (1 - d)) & 1);
child.uv[d][m] = uvMid.get(d);
child.uv[d][1 - m] = uv[d][1 - m];
}
}
return true;
}
/**
* Return the direction vector corresponding to the center in (s,t)-space of
* the given cell. This is the point at which the cell is divided into four
* subcells; it is not necessarily the centroid of the cell in (u,v)-space or
* (x,y,z)-space. The point returned by GetCenterRaw is not necessarily unit
* length.
*/
public S2Point getCenter() {
return S2Point.normalize(getCenterRaw());
}
public S2Point getCenterRaw() {
return cellId.toPointRaw();
}
/**
* Return the center of the cell in (u,v) coordinates (see {@code
* S2Projections}). Note that the center of the cell is defined as the point
* at which it is recursively subdivided into four children; in general, it is
* not at the midpoint of the (u,v) rectangle covered by the cell
*/
public R2Vector getCenterUV() {
MutableInteger i = new MutableInteger(0);
MutableInteger j = new MutableInteger(0);
cellId.toFaceIJOrientation(i, j, null);
int cellSize = 1 << (S2CellId.MAX_LEVEL - level);
// TODO(dbeaumont): Figure out a better naming of the variables here (and elsewhere).
int si = (i.intValue() & -cellSize) * 2 + cellSize - MAX_CELL_SIZE;
double x = S2Projections.stToUV((1.0 / MAX_CELL_SIZE) * si);
int sj = (j.intValue() & -cellSize) * 2 + cellSize - MAX_CELL_SIZE;
double y = S2Projections.stToUV((1.0 / MAX_CELL_SIZE) * sj);
return new R2Vector(x, y);
}
/**
* Return the average area for cells at the given level.
*/
public static double averageArea(int level) {
return S2Projections.AVG_AREA.getValue(level);
}
/**
* Return the average area of cells at this level. This is accurate to within
* a factor of 1.7 (for S2_QUADRATIC_PROJECTION) and is extremely cheap to
* compute.
*/
public double averageArea() {
return averageArea(level);
}
/**
* Return the approximate area of this cell. This method is accurate to within
* 3% percent for all cell sizes and accurate to within 0.1% for cells at
* level 5 or higher (i.e. 300km square or smaller). It is moderately cheap to
* compute.
*/
public double approxArea() {
// All cells at the first two levels have the same area.
if (level < 2) {
return averageArea(level);
}
// First, compute the approximate area of the cell when projected
// perpendicular to its normal. The cross product of its diagonals gives
// the normal, and the length of the normal is twice the projected area.
double flatArea = 0.5 * S2Point.crossProd(
S2Point.sub(getVertex(2), getVertex(0)), S2Point.sub(getVertex(3), getVertex(1))).norm();
// Now, compensate for the curvature of the cell surface by pretending
// that the cell is shaped like a spherical cap. The ratio of the
// area of a spherical cap to the area of its projected disc turns out
// to be 2 / (1 + sqrt(1 - r*r)) where "r" is the radius of the disc.
// For example, when r=0 the ratio is 1, and when r=1 the ratio is 2.
// Here we set Pi*r*r == flat_area to find the equivalent disc.
return flatArea * 2 / (1 + Math.sqrt(1 - Math.min(S2.M_1_PI * flatArea, 1.0)));
}
/**
* Return the area of this cell as accurately as possible. This method is more
* expensive but it is accurate to 6 digits of precision even for leaf cells
* (whose area is approximately 1e-18).
*/
public double exactArea() {
S2Point v0 = getVertex(0);
S2Point v1 = getVertex(1);
S2Point v2 = getVertex(2);
S2Point v3 = getVertex(3);
return S2.area(v0, v1, v2) + S2.area(v0, v2, v3);
}
// //////////////////////////////////////////////////////////////////////
// S2Region interface (see {@code S2Region} for details):
@Override
public S2Region clone() {
S2Cell clone = new S2Cell();
clone.face = this.face;
clone.level = this.level;
clone.orientation = this.orientation;
clone.uv = this.uv.clone();
return clone;
}
@Override
public S2Cap getCapBound() {
// Use the cell center in (u,v)-space as the cap axis. This vector is
// very close to GetCenter() and faster to compute. Neither one of these
// vectors yields the bounding cap with minimal surface area, but they
// are both pretty close.
//
// It's possible to show that the two vertices that are furthest from
// the (u,v)-origin never determine the maximum cap size (this is a
// possible future optimization).
double u = 0.5 * (uv[0][0] + uv[0][1]);
double v = 0.5 * (uv[1][0] + uv[1][1]);
S2Cap cap = S2Cap.fromAxisHeight(S2Point.normalize(S2Projections.faceUvToXyz(face, u, v)), 0);
for (int k = 0; k < 4; ++k) {
cap = cap.addPoint(getVertex(k));
}
return cap;
}
// We grow the bounds slightly to make sure that the bounding rectangle
// also contains the normalized versions of the vertices. Note that the
// maximum result magnitude is Pi, with a floating-point exponent of 1.
// Therefore adding or subtracting 2**-51 will always change the result.
private static final double MAX_ERROR = 1.0 / (1L << 51);
// The 4 cells around the equator extend to +/-45 degrees latitude at the
// midpoints of their top and bottom edges. The two cells covering the
// poles extend down to +/-35.26 degrees at their vertices.
// adding kMaxError (as opposed to the C version) because of asin and atan2
// roundoff errors
private static final double POLE_MIN_LAT = Math.asin(Math.sqrt(1.0 / 3.0)) - MAX_ERROR;
// 35.26 degrees
@Override
public S2LatLngRect getRectBound() {
if (level > 0) {
// Except for cells at level 0, the latitude and longitude extremes are
// attained at the vertices. Furthermore, the latitude range is
// determined by one pair of diagonally opposite vertices and the
// longitude range is determined by the other pair.
//
// We first determine which corner (i,j) of the cell has the largest
// absolute latitude. To maximize latitude, we want to find the point in
// the cell that has the largest absolute z-coordinate and the smallest
// absolute x- and y-coordinates. To do this we look at each coordinate
// (u and v), and determine whether we want to minimize or maximize that
// coordinate based on the axis direction and the cell's (u,v) quadrant.
double u = uv[0][0] + uv[0][1];
double v = uv[1][0] + uv[1][1];
int i = S2Projections.getUAxis(face).z == 0 ? (u < 0 ? 1 : 0) : (u > 0 ? 1 : 0);
int j = S2Projections.getVAxis(face).z == 0 ? (v < 0 ? 1 : 0) : (v > 0 ? 1 : 0);
R1Interval lat = R1Interval.fromPointPair(getLatitude(i, j), getLatitude(1 - i, 1 - j));
lat = lat.expanded(MAX_ERROR).intersection(S2LatLngRect.fullLat());
if (lat.lo() == -S2.M_PI_2 || lat.hi() == S2.M_PI_2) {
return new S2LatLngRect(lat, S1Interval.full());
}
S1Interval lng = S1Interval.fromPointPair(getLongitude(i, 1 - j), getLongitude(1 - i, j));
return new S2LatLngRect(lat, lng.expanded(MAX_ERROR));
}
// The face centers are the +X, +Y, +Z, -X, -Y, -Z axes in that order.
// assert (S2Projections.getNorm(face).get(face % 3) == ((face < 3) ? 1 : -1));
switch (face) {
case 0:
return new S2LatLngRect(
new R1Interval(-S2.M_PI_4, S2.M_PI_4), new S1Interval(-S2.M_PI_4, S2.M_PI_4));
case 1:
return new S2LatLngRect(
new R1Interval(-S2.M_PI_4, S2.M_PI_4), new S1Interval(S2.M_PI_4, 3 * S2.M_PI_4));
case 2:
return new S2LatLngRect(
new R1Interval(POLE_MIN_LAT, S2.M_PI_2), new S1Interval(-S2.M_PI, S2.M_PI));
case 3:
return new S2LatLngRect(
new R1Interval(-S2.M_PI_4, S2.M_PI_4), new S1Interval(3 * S2.M_PI_4, -3 * S2.M_PI_4));
case 4:
return new S2LatLngRect(
new R1Interval(-S2.M_PI_4, S2.M_PI_4), new S1Interval(-3 * S2.M_PI_4, -S2.M_PI_4));
default:
return new S2LatLngRect(
new R1Interval(-S2.M_PI_2, -POLE_MIN_LAT), new S1Interval(-S2.M_PI, S2.M_PI));
}
}
@Override
public boolean mayIntersect(S2Cell cell) {
return cellId.intersects(cell.cellId);
}
public boolean contains(S2Point p) {
// We can't just call XYZtoFaceUV, because for points that lie on the
// boundary between two faces (i.e. u or v is +1/-1) we need to return
// true for both adjacent cells.
R2Vector uvPoint = S2Projections.faceXyzToUv(face, p);
if (uvPoint == null) {
return false;
}
return (uvPoint.x() >= uv[0][0] && uvPoint.x() <= uv[0][1]
&& uvPoint.y() >= uv[1][0] && uvPoint.y() <= uv[1][1]);
}
// The point 'p' does not need to be normalized.
@Override
public boolean contains(S2Cell cell) {
return cellId.contains(cell.cellId);
}
private void init(S2CellId id) {
cellId = id;
MutableInteger ij[] = new MutableInteger[2];
MutableInteger mOrientation = new MutableInteger(0);
for (int d = 0; d < 2; ++d) {
ij[d] = new MutableInteger(0);
}
face = (byte) id.toFaceIJOrientation(ij[0], ij[1], mOrientation);
orientation = (byte) mOrientation.intValue(); // Compress int to a byte.
level = (byte) id.level();
int cellSize = 1 << (S2CellId.MAX_LEVEL - level);
for (int d = 0; d < 2; ++d) {
// Compute the cell bounds in scaled (i,j) coordinates.
int sijLo = (ij[d].intValue() & -cellSize) * 2 - MAX_CELL_SIZE;
int sijHi = sijLo + cellSize * 2;
uv[d][0] = S2Projections.stToUV((1.0 / MAX_CELL_SIZE) * sijLo);
uv[d][1] = S2Projections.stToUV((1.0 / MAX_CELL_SIZE) * sijHi);
}
}
// Internal method that does the actual work in the constructors.
private double getLatitude(int i, int j) {
S2Point p = S2Projections.faceUvToXyz(face, uv[0][i], uv[1][j]);
return Math.atan2(p.z, Math.sqrt(p.x * p.x + p.y * p.y));
}
private double getLongitude(int i, int j) {
S2Point p = S2Projections.faceUvToXyz(face, uv[0][i], uv[1][j]);
return Math.atan2(p.y, p.x);
}
// Return the latitude or longitude of the cell vertex given by (i,j),
// where "i" and "j" are either 0 or 1.
@Override
public String toString() {
return "[" + face + ", " + level + ", " + orientation + ", " + cellId + "]";
}
@Override
public int hashCode() {
int value = 17;
value = 37 * (37 * (37 * value + face) + orientation) + level;
return 37 * value + id().hashCode();
}
@Override
public boolean equals(Object that) {
if (that instanceof S2Cell) {
S2Cell thatCell = (S2Cell) that;
return this.face == thatCell.face && this.level == thatCell.level
&& this.orientation == thatCell.orientation && this.cellId.equals(thatCell.cellId);
}
return false;
}
}
© 2015 - 2025 Weber Informatics LLC | Privacy Policy