Michael Bolin | 6fb4657 | 2011-09-22 15:24:39 -0400 | [diff] [blame] | 1 | /* |
| 2 | * Copyright 2005 Google Inc. |
| 3 | * |
| 4 | * Licensed under the Apache License, Version 2.0 (the "License"); |
| 5 | * you may not use this file except in compliance with the License. |
| 6 | * You may obtain a copy of the License at |
| 7 | * |
| 8 | * http://www.apache.org/licenses/LICENSE-2.0 |
| 9 | * |
| 10 | * Unless required by applicable law or agreed to in writing, software |
| 11 | * distributed under the License is distributed on an "AS IS" BASIS, |
| 12 | * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| 13 | * See the License for the specific language governing permissions and |
| 14 | * limitations under the License. |
| 15 | */ |
| 16 | package com.google.common.geometry; |
| 17 | |
| 18 | |
| 19 | /** |
| 20 | * An S2Cell is an S2Region object that represents a cell. Unlike S2CellIds, it |
| 21 | * supports efficient containment and intersection tests. However, it is also a |
| 22 | * more expensive representation. |
| 23 | * |
| 24 | */ |
| 25 | |
David Beaumont | 0095e3d | 2011-10-31 18:43:31 -0400 | [diff] [blame] | 26 | public final strictfp class S2Cell implements S2Region { |
Michael Bolin | 6fb4657 | 2011-09-22 15:24:39 -0400 | [diff] [blame] | 27 | |
| 28 | private static final int MAX_CELL_SIZE = 1 << S2CellId.MAX_LEVEL; |
| 29 | |
Michael Bolin | 6fb4657 | 2011-09-22 15:24:39 -0400 | [diff] [blame] | 30 | byte face; |
| 31 | byte level; |
| 32 | byte orientation; |
| 33 | S2CellId cellId; |
| 34 | double[][] uv = new double[2][2]; |
| 35 | |
| 36 | /** |
| 37 | * Default constructor used only internally. |
| 38 | */ |
| 39 | S2Cell() { |
| 40 | } |
| 41 | |
| 42 | /** |
| 43 | * An S2Cell always corresponds to a particular S2CellId. The other |
| 44 | * constructors are just convenience methods. |
| 45 | */ |
| 46 | public S2Cell(S2CellId id) { |
| 47 | init(id); |
| 48 | } |
| 49 | |
| 50 | // This is a static method in order to provide named parameters. |
| 51 | public static S2Cell fromFacePosLevel(int face, byte pos, int level) { |
| 52 | return new S2Cell(S2CellId.fromFacePosLevel(face, pos, level)); |
| 53 | } |
| 54 | |
| 55 | // Convenience methods. |
| 56 | public S2Cell(S2Point p) { |
| 57 | init(S2CellId.fromPoint(p)); |
| 58 | } |
| 59 | |
| 60 | public S2Cell(S2LatLng ll) { |
| 61 | init(S2CellId.fromLatLng(ll)); |
| 62 | } |
| 63 | |
| 64 | |
| 65 | public S2CellId id() { |
| 66 | return cellId; |
| 67 | } |
| 68 | |
| 69 | public int face() { |
| 70 | return face; |
| 71 | } |
| 72 | |
| 73 | public byte level() { |
| 74 | return level; |
| 75 | } |
| 76 | |
| 77 | public byte orientation() { |
| 78 | return orientation; |
| 79 | } |
| 80 | |
| 81 | public boolean isLeaf() { |
| 82 | return level == S2CellId.MAX_LEVEL; |
| 83 | } |
| 84 | |
| 85 | public S2Point getVertex(int k) { |
| 86 | return S2Point.normalize(getVertexRaw(k)); |
| 87 | } |
| 88 | |
| 89 | /** |
| 90 | * Return the k-th vertex of the cell (k = 0,1,2,3). Vertices are returned in |
| 91 | * CCW order. The points returned by GetVertexRaw are not necessarily unit |
| 92 | * length. |
| 93 | */ |
| 94 | public S2Point getVertexRaw(int k) { |
| 95 | // Vertices are returned in the order SW, SE, NE, NW. |
| 96 | return S2Projections.faceUvToXyz(face, uv[0][(k >> 1) ^ (k & 1)], uv[1][k >> 1]); |
| 97 | } |
| 98 | |
| 99 | public S2Point getEdge(int k) { |
| 100 | return S2Point.normalize(getEdgeRaw(k)); |
| 101 | } |
| 102 | |
| 103 | public S2Point getEdgeRaw(int k) { |
| 104 | switch (k) { |
| 105 | case 0: |
| 106 | return S2Projections.getVNorm(face, uv[1][0]); // South |
| 107 | case 1: |
| 108 | return S2Projections.getUNorm(face, uv[0][1]); // East |
| 109 | case 2: |
| 110 | return S2Point.neg(S2Projections.getVNorm(face, uv[1][1])); // North |
| 111 | default: |
| 112 | return S2Point.neg(S2Projections.getUNorm(face, uv[0][0])); // West |
| 113 | } |
| 114 | } |
| 115 | |
| 116 | /** |
| 117 | * Return the inward-facing normal of the great circle passing through the |
| 118 | * edge from vertex k to vertex k+1 (mod 4). The normals returned by |
| 119 | * GetEdgeRaw are not necessarily unit length. |
| 120 | * |
| 121 | * If this is not a leaf cell, set children[0..3] to the four children of |
| 122 | * this cell (in traversal order) and return true. Otherwise returns false. |
| 123 | * This method is equivalent to the following: |
| 124 | * |
| 125 | * for (pos=0, id=child_begin(); id != child_end(); id = id.next(), ++pos) |
| 126 | * children[i] = S2Cell(id); |
| 127 | * |
| 128 | * except that it is more than two times faster. |
| 129 | */ |
| 130 | public boolean subdivide(S2Cell children[]) { |
| 131 | // This function is equivalent to just iterating over the child cell ids |
| 132 | // and calling the S2Cell constructor, but it is about 2.5 times faster. |
| 133 | |
| 134 | if (cellId.isLeaf()) { |
| 135 | return false; |
| 136 | } |
| 137 | |
| 138 | // Compute the cell midpoint in uv-space. |
| 139 | R2Vector uvMid = getCenterUV(); |
| 140 | |
| 141 | // Create four children with the appropriate bounds. |
| 142 | S2CellId id = cellId.childBegin(); |
| 143 | for (int pos = 0; pos < 4; ++pos, id = id.next()) { |
| 144 | S2Cell child = children[pos]; |
| 145 | child.face = face; |
| 146 | child.level = (byte) (level + 1); |
David Beaumont | 89b8770 | 2011-10-31 18:44:08 -0400 | [diff] [blame] | 147 | child.orientation = (byte) (orientation ^ S2.posToOrientation(pos)); |
Michael Bolin | 6fb4657 | 2011-09-22 15:24:39 -0400 | [diff] [blame] | 148 | child.cellId = id; |
David Beaumont | 89b8770 | 2011-10-31 18:44:08 -0400 | [diff] [blame] | 149 | int ij = S2.posToIJ(orientation, pos); |
Michael Bolin | 6fb4657 | 2011-09-22 15:24:39 -0400 | [diff] [blame] | 150 | for (int d = 0; d < 2; ++d) { |
| 151 | // The dimension 0 index (i/u) is in bit 1 of ij. |
| 152 | int m = 1 - ((ij >> (1 - d)) & 1); |
| 153 | child.uv[d][m] = uvMid.get(d); |
| 154 | child.uv[d][1 - m] = uv[d][1 - m]; |
| 155 | } |
| 156 | } |
| 157 | return true; |
| 158 | } |
| 159 | |
| 160 | /** |
| 161 | * Return the direction vector corresponding to the center in (s,t)-space of |
| 162 | * the given cell. This is the point at which the cell is divided into four |
| 163 | * subcells; it is not necessarily the centroid of the cell in (u,v)-space or |
| 164 | * (x,y,z)-space. The point returned by GetCenterRaw is not necessarily unit |
| 165 | * length. |
| 166 | */ |
| 167 | public S2Point getCenter() { |
| 168 | return S2Point.normalize(getCenterRaw()); |
| 169 | } |
| 170 | |
| 171 | public S2Point getCenterRaw() { |
| 172 | return cellId.toPointRaw(); |
| 173 | } |
| 174 | |
| 175 | /** |
| 176 | * Return the center of the cell in (u,v) coordinates (see {@code |
| 177 | * S2Projections}). Note that the center of the cell is defined as the point |
| 178 | * at which it is recursively subdivided into four children; in general, it is |
| 179 | * not at the midpoint of the (u,v) rectangle covered by the cell |
| 180 | */ |
| 181 | public R2Vector getCenterUV() { |
| 182 | MutableInteger i = new MutableInteger(0); |
| 183 | MutableInteger j = new MutableInteger(0); |
Michael Bolin | 6fb4657 | 2011-09-22 15:24:39 -0400 | [diff] [blame] | 184 | cellId.toFaceIJOrientation(i, j, null); |
| 185 | int cellSize = 1 << (S2CellId.MAX_LEVEL - level); |
| 186 | |
David Beaumont | 0095e3d | 2011-10-31 18:43:31 -0400 | [diff] [blame] | 187 | // TODO(dbeaumont): Figure out a better naming of the variables here (and elsewhere). |
| 188 | int si = (i.intValue() & -cellSize) * 2 + cellSize - MAX_CELL_SIZE; |
| 189 | double x = S2Projections.stToUV((1.0 / MAX_CELL_SIZE) * si); |
Michael Bolin | 6fb4657 | 2011-09-22 15:24:39 -0400 | [diff] [blame] | 190 | |
David Beaumont | 0095e3d | 2011-10-31 18:43:31 -0400 | [diff] [blame] | 191 | int sj = (j.intValue() & -cellSize) * 2 + cellSize - MAX_CELL_SIZE; |
| 192 | double y = S2Projections.stToUV((1.0 / MAX_CELL_SIZE) * sj); |
Michael Bolin | 6fb4657 | 2011-09-22 15:24:39 -0400 | [diff] [blame] | 193 | |
David Beaumont | 0095e3d | 2011-10-31 18:43:31 -0400 | [diff] [blame] | 194 | return new R2Vector(x, y); |
Michael Bolin | 6fb4657 | 2011-09-22 15:24:39 -0400 | [diff] [blame] | 195 | } |
| 196 | |
| 197 | /** |
| 198 | * Return the average area for cells at the given level. |
| 199 | */ |
| 200 | public static double averageArea(int level) { |
| 201 | return S2Projections.AVG_AREA.getValue(level); |
| 202 | } |
| 203 | |
| 204 | /** |
| 205 | * Return the average area of cells at this level. This is accurate to within |
| 206 | * a factor of 1.7 (for S2_QUADRATIC_PROJECTION) and is extremely cheap to |
| 207 | * compute. |
| 208 | */ |
| 209 | public double averageArea() { |
| 210 | return averageArea(level); |
| 211 | } |
| 212 | |
| 213 | /** |
| 214 | * Return the approximate area of this cell. This method is accurate to within |
| 215 | * 3% percent for all cell sizes and accurate to within 0.1% for cells at |
| 216 | * level 5 or higher (i.e. 300km square or smaller). It is moderately cheap to |
| 217 | * compute. |
| 218 | */ |
| 219 | public double approxArea() { |
| 220 | |
| 221 | // All cells at the first two levels have the same area. |
| 222 | if (level < 2) { |
| 223 | return averageArea(level); |
| 224 | } |
| 225 | |
| 226 | // First, compute the approximate area of the cell when projected |
| 227 | // perpendicular to its normal. The cross product of its diagonals gives |
| 228 | // the normal, and the length of the normal is twice the projected area. |
| 229 | double flatArea = 0.5 * S2Point.crossProd( |
| 230 | S2Point.sub(getVertex(2), getVertex(0)), S2Point.sub(getVertex(3), getVertex(1))).norm(); |
| 231 | |
| 232 | // Now, compensate for the curvature of the cell surface by pretending |
| 233 | // that the cell is shaped like a spherical cap. The ratio of the |
| 234 | // area of a spherical cap to the area of its projected disc turns out |
| 235 | // to be 2 / (1 + sqrt(1 - r*r)) where "r" is the radius of the disc. |
| 236 | // For example, when r=0 the ratio is 1, and when r=1 the ratio is 2. |
| 237 | // Here we set Pi*r*r == flat_area to find the equivalent disc. |
| 238 | return flatArea * 2 / (1 + Math.sqrt(1 - Math.min(S2.M_1_PI * flatArea, 1.0))); |
| 239 | } |
| 240 | |
| 241 | /** |
| 242 | * Return the area of this cell as accurately as possible. This method is more |
| 243 | * expensive but it is accurate to 6 digits of precision even for leaf cells |
| 244 | * (whose area is approximately 1e-18). |
| 245 | */ |
| 246 | public double exactArea() { |
| 247 | S2Point v0 = getVertex(0); |
| 248 | S2Point v1 = getVertex(1); |
| 249 | S2Point v2 = getVertex(2); |
| 250 | S2Point v3 = getVertex(3); |
| 251 | return S2.area(v0, v1, v2) + S2.area(v0, v2, v3); |
| 252 | } |
| 253 | |
| 254 | // ////////////////////////////////////////////////////////////////////// |
| 255 | // S2Region interface (see {@code S2Region} for details): |
| 256 | |
| 257 | @Override |
| 258 | public S2Region clone() { |
| 259 | S2Cell clone = new S2Cell(); |
| 260 | clone.face = this.face; |
| 261 | clone.level = this.level; |
| 262 | clone.orientation = this.orientation; |
| 263 | clone.uv = this.uv.clone(); |
| 264 | |
| 265 | return clone; |
| 266 | } |
| 267 | |
| 268 | @Override |
| 269 | public S2Cap getCapBound() { |
| 270 | // Use the cell center in (u,v)-space as the cap axis. This vector is |
| 271 | // very close to GetCenter() and faster to compute. Neither one of these |
| 272 | // vectors yields the bounding cap with minimal surface area, but they |
| 273 | // are both pretty close. |
| 274 | // |
| 275 | // It's possible to show that the two vertices that are furthest from |
| 276 | // the (u,v)-origin never determine the maximum cap size (this is a |
| 277 | // possible future optimization). |
| 278 | |
| 279 | double u = 0.5 * (uv[0][0] + uv[0][1]); |
| 280 | double v = 0.5 * (uv[1][0] + uv[1][1]); |
| 281 | S2Cap cap = S2Cap.fromAxisHeight(S2Point.normalize(S2Projections.faceUvToXyz(face, u, v)), 0); |
| 282 | for (int k = 0; k < 4; ++k) { |
| 283 | cap = cap.addPoint(getVertex(k)); |
| 284 | } |
| 285 | return cap; |
| 286 | } |
| 287 | |
| 288 | // We grow the bounds slightly to make sure that the bounding rectangle |
| 289 | // also contains the normalized versions of the vertices. Note that the |
| 290 | // maximum result magnitude is Pi, with a floating-point exponent of 1. |
| 291 | // Therefore adding or subtracting 2**-51 will always change the result. |
| 292 | private static final double MAX_ERROR = 1.0 / (1L << 51); |
| 293 | |
| 294 | // The 4 cells around the equator extend to +/-45 degrees latitude at the |
| 295 | // midpoints of their top and bottom edges. The two cells covering the |
| 296 | // poles extend down to +/-35.26 degrees at their vertices. |
| 297 | // adding kMaxError (as opposed to the C version) because of asin and atan2 |
| 298 | // roundoff errors |
| 299 | private static final double POLE_MIN_LAT = Math.asin(Math.sqrt(1.0 / 3.0)) - MAX_ERROR; |
| 300 | // 35.26 degrees |
| 301 | |
| 302 | |
| 303 | @Override |
| 304 | public S2LatLngRect getRectBound() { |
| 305 | if (level > 0) { |
| 306 | // Except for cells at level 0, the latitude and longitude extremes are |
| 307 | // attained at the vertices. Furthermore, the latitude range is |
| 308 | // determined by one pair of diagonally opposite vertices and the |
| 309 | // longitude range is determined by the other pair. |
| 310 | // |
| 311 | // We first determine which corner (i,j) of the cell has the largest |
| 312 | // absolute latitude. To maximize latitude, we want to find the point in |
| 313 | // the cell that has the largest absolute z-coordinate and the smallest |
| 314 | // absolute x- and y-coordinates. To do this we look at each coordinate |
| 315 | // (u and v), and determine whether we want to minimize or maximize that |
| 316 | // coordinate based on the axis direction and the cell's (u,v) quadrant. |
| 317 | double u = uv[0][0] + uv[0][1]; |
| 318 | double v = uv[1][0] + uv[1][1]; |
| 319 | int i = S2Projections.getUAxis(face).z == 0 ? (u < 0 ? 1 : 0) : (u > 0 ? 1 : 0); |
| 320 | int j = S2Projections.getVAxis(face).z == 0 ? (v < 0 ? 1 : 0) : (v > 0 ? 1 : 0); |
| 321 | |
| 322 | |
| 323 | R1Interval lat = R1Interval.fromPointPair(getLatitude(i, j), getLatitude(1 - i, 1 - j)); |
| 324 | lat = lat.expanded(MAX_ERROR).intersection(S2LatLngRect.fullLat()); |
| 325 | if (lat.lo() == -S2.M_PI_2 || lat.hi() == S2.M_PI_2) { |
| 326 | return new S2LatLngRect(lat, S1Interval.full()); |
| 327 | } |
| 328 | S1Interval lng = S1Interval.fromPointPair(getLongitude(i, 1 - j), getLongitude(1 - i, j)); |
| 329 | return new S2LatLngRect(lat, lng.expanded(MAX_ERROR)); |
| 330 | } |
| 331 | |
| 332 | |
| 333 | // The face centers are the +X, +Y, +Z, -X, -Y, -Z axes in that order. |
| 334 | // assert (S2Projections.getNorm(face).get(face % 3) == ((face < 3) ? 1 : -1)); |
| 335 | switch (face) { |
| 336 | case 0: |
| 337 | return new S2LatLngRect( |
| 338 | new R1Interval(-S2.M_PI_4, S2.M_PI_4), new S1Interval(-S2.M_PI_4, S2.M_PI_4)); |
| 339 | case 1: |
| 340 | return new S2LatLngRect( |
| 341 | new R1Interval(-S2.M_PI_4, S2.M_PI_4), new S1Interval(S2.M_PI_4, 3 * S2.M_PI_4)); |
| 342 | case 2: |
| 343 | return new S2LatLngRect( |
| 344 | new R1Interval(POLE_MIN_LAT, S2.M_PI_2), new S1Interval(-S2.M_PI, S2.M_PI)); |
| 345 | case 3: |
| 346 | return new S2LatLngRect( |
| 347 | new R1Interval(-S2.M_PI_4, S2.M_PI_4), new S1Interval(3 * S2.M_PI_4, -3 * S2.M_PI_4)); |
| 348 | case 4: |
| 349 | return new S2LatLngRect( |
| 350 | new R1Interval(-S2.M_PI_4, S2.M_PI_4), new S1Interval(-3 * S2.M_PI_4, -S2.M_PI_4)); |
| 351 | default: |
| 352 | return new S2LatLngRect( |
| 353 | new R1Interval(-S2.M_PI_2, -POLE_MIN_LAT), new S1Interval(-S2.M_PI, S2.M_PI)); |
| 354 | } |
| 355 | |
| 356 | } |
| 357 | |
| 358 | @Override |
| 359 | public boolean mayIntersect(S2Cell cell) { |
| 360 | return cellId.intersects(cell.cellId); |
| 361 | } |
| 362 | |
| 363 | public boolean contains(S2Point p) { |
| 364 | // We can't just call XYZtoFaceUV, because for points that lie on the |
| 365 | // boundary between two faces (i.e. u or v is +1/-1) we need to return |
| 366 | // true for both adjacent cells. |
David Beaumont | 0095e3d | 2011-10-31 18:43:31 -0400 | [diff] [blame] | 367 | R2Vector uvPoint = S2Projections.faceXyzToUv(face, p); |
| 368 | if (uvPoint == null) { |
Michael Bolin | 6fb4657 | 2011-09-22 15:24:39 -0400 | [diff] [blame] | 369 | return false; |
| 370 | } |
David Beaumont | c04b68b | 2011-11-08 11:05:51 -0500 | [diff] [blame] | 371 | return (uvPoint.x() >= uv[0][0] && uvPoint.x() <= uv[0][1] |
| 372 | && uvPoint.y() >= uv[1][0] && uvPoint.y() <= uv[1][1]); |
Michael Bolin | 6fb4657 | 2011-09-22 15:24:39 -0400 | [diff] [blame] | 373 | } |
| 374 | |
| 375 | // The point 'p' does not need to be normalized. |
| 376 | @Override |
| 377 | public boolean contains(S2Cell cell) { |
| 378 | return cellId.contains(cell.cellId); |
| 379 | } |
| 380 | |
| 381 | private void init(S2CellId id) { |
| 382 | cellId = id; |
| 383 | MutableInteger ij[] = new MutableInteger[2]; |
| 384 | MutableInteger mOrientation = new MutableInteger(0); |
| 385 | |
| 386 | for (int d = 0; d < 2; ++d) { |
| 387 | ij[d] = new MutableInteger(0); |
| 388 | } |
| 389 | |
| 390 | face = (byte) id.toFaceIJOrientation(ij[0], ij[1], mOrientation); |
| 391 | orientation = (byte) mOrientation.intValue(); // Compress int to a byte. |
| 392 | level = (byte) id.level(); |
| 393 | int cellSize = 1 << (S2CellId.MAX_LEVEL - level); |
| 394 | for (int d = 0; d < 2; ++d) { |
| 395 | // Compute the cell bounds in scaled (i,j) coordinates. |
| 396 | int sijLo = (ij[d].intValue() & -cellSize) * 2 - MAX_CELL_SIZE; |
| 397 | int sijHi = sijLo + cellSize * 2; |
| 398 | uv[d][0] = S2Projections.stToUV((1.0 / MAX_CELL_SIZE) * sijLo); |
| 399 | uv[d][1] = S2Projections.stToUV((1.0 / MAX_CELL_SIZE) * sijHi); |
| 400 | } |
| 401 | } |
| 402 | |
| 403 | |
| 404 | // Internal method that does the actual work in the constructors. |
| 405 | |
| 406 | private double getLatitude(int i, int j) { |
| 407 | S2Point p = S2Projections.faceUvToXyz(face, uv[0][i], uv[1][j]); |
| 408 | return Math.atan2(p.z, Math.sqrt(p.x * p.x + p.y * p.y)); |
| 409 | } |
| 410 | |
| 411 | private double getLongitude(int i, int j) { |
| 412 | S2Point p = S2Projections.faceUvToXyz(face, uv[0][i], uv[1][j]); |
| 413 | return Math.atan2(p.y, p.x); |
| 414 | } |
| 415 | |
| 416 | // Return the latitude or longitude of the cell vertex given by (i,j), |
| 417 | // where "i" and "j" are either 0 or 1. |
| 418 | |
| 419 | @Override |
| 420 | public String toString() { |
| 421 | return "[" + face + ", " + level + ", " + orientation + ", " + cellId + "]"; |
| 422 | } |
| 423 | |
| 424 | @Override |
| 425 | public int hashCode() { |
| 426 | int value = 17; |
| 427 | value = 37 * (37 * (37 * value + face) + orientation) + level; |
| 428 | return 37 * value + id().hashCode(); |
| 429 | } |
| 430 | |
| 431 | @Override |
| 432 | public boolean equals(Object that) { |
| 433 | if (that instanceof S2Cell) { |
| 434 | S2Cell thatCell = (S2Cell) that; |
| 435 | return this.face == thatCell.face && this.level == thatCell.level |
| 436 | && this.orientation == thatCell.orientation && this.cellId.equals(thatCell.cellId); |
| 437 | } |
| 438 | return false; |
| 439 | } |
| 440 | |
| 441 | } |