| /* |
| * 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; |
| } |
| |
| } |