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 26 public final strictfp class S2Cell implements S2Region { 27 28 private static final int MAX_CELL_SIZE = 1 << S2CellId.MAX_LEVEL; 29 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 */ S2Cell()39 S2Cell() { 40 } 41 42 /** 43 * An S2Cell always corresponds to a particular S2CellId. The other 44 * constructors are just convenience methods. 45 */ S2Cell(S2CellId id)46 public S2Cell(S2CellId id) { 47 init(id); 48 } 49 50 // This is a static method in order to provide named parameters. fromFacePosLevel(int face, byte pos, int level)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. S2Cell(S2Point p)56 public S2Cell(S2Point p) { 57 init(S2CellId.fromPoint(p)); 58 } 59 S2Cell(S2LatLng ll)60 public S2Cell(S2LatLng ll) { 61 init(S2CellId.fromLatLng(ll)); 62 } 63 64 id()65 public S2CellId id() { 66 return cellId; 67 } 68 face()69 public int face() { 70 return face; 71 } 72 level()73 public byte level() { 74 return level; 75 } 76 orientation()77 public byte orientation() { 78 return orientation; 79 } 80 isLeaf()81 public boolean isLeaf() { 82 return level == S2CellId.MAX_LEVEL; 83 } 84 getVertex(int k)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 */ getVertexRaw(int k)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 getEdge(int k)99 public S2Point getEdge(int k) { 100 return S2Point.normalize(getEdgeRaw(k)); 101 } 102 getEdgeRaw(int k)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 */ subdivide(S2Cell children[])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); 147 child.orientation = (byte) (orientation ^ S2.posToOrientation(pos)); 148 child.cellId = id; 149 int ij = S2.posToIJ(orientation, pos); 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 */ getCenter()167 public S2Point getCenter() { 168 return S2Point.normalize(getCenterRaw()); 169 } 170 getCenterRaw()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 */ getCenterUV()181 public R2Vector getCenterUV() { 182 MutableInteger i = new MutableInteger(0); 183 MutableInteger j = new MutableInteger(0); 184 cellId.toFaceIJOrientation(i, j, null); 185 int cellSize = 1 << (S2CellId.MAX_LEVEL - level); 186 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); 190 191 int sj = (j.intValue() & -cellSize) * 2 + cellSize - MAX_CELL_SIZE; 192 double y = S2Projections.stToUV((1.0 / MAX_CELL_SIZE) * sj); 193 194 return new R2Vector(x, y); 195 } 196 197 /** 198 * Return the average area for cells at the given level. 199 */ averageArea(int level)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 */ averageArea()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 */ approxArea()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 */ exactArea()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 clone()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 getCapBound()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 getRectBound()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 mayIntersect(S2Cell cell)359 public boolean mayIntersect(S2Cell cell) { 360 return cellId.intersects(cell.cellId); 361 } 362 contains(S2Point p)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. 367 R2Vector uvPoint = S2Projections.faceXyzToUv(face, p); 368 if (uvPoint == null) { 369 return false; 370 } 371 return (uvPoint.x() >= uv[0][0] && uvPoint.x() <= uv[0][1] 372 && uvPoint.y() >= uv[1][0] && uvPoint.y() <= uv[1][1]); 373 } 374 375 // The point 'p' does not need to be normalized. 376 @Override contains(S2Cell cell)377 public boolean contains(S2Cell cell) { 378 return cellId.contains(cell.cellId); 379 } 380 init(S2CellId id)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 getLatitude(int i, int j)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 getLongitude(int i, int j)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 toString()420 public String toString() { 421 return "[" + face + ", " + level + ", " + orientation + ", " + cellId + "]"; 422 } 423 424 @Override hashCode()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 equals(Object that)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 } 442