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