1 // Copyright 2014 Google Inc. All rights reserved. 2 // 3 // Licensed under the Apache License, Version 2.0 (the "License"); 4 // you may not use this file except in compliance with the License. 5 // You may obtain a copy of the License at 6 // 7 // http://www.apache.org/licenses/LICENSE-2.0 8 // 9 // Unless required by applicable law or agreed to in writing, software 10 // distributed under the License is distributed on an "AS IS" BASIS, 11 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 12 // See the License for the specific language governing permissions and 13 // limitations under the License. 14 15 package s2 16 17 import ( 18 "math" 19 20 "github.com/golang/geo/r3" 21 ) 22 23 // 24 // This file contains documentation of the various coordinate systems used 25 // throughout the library. Most importantly, S2 defines a framework for 26 // decomposing the unit sphere into a hierarchy of "cells". Each cell is a 27 // quadrilateral bounded by four geodesics. The top level of the hierarchy is 28 // obtained by projecting the six faces of a cube onto the unit sphere, and 29 // lower levels are obtained by subdividing each cell into four children 30 // recursively. Cells are numbered such that sequentially increasing cells 31 // follow a continuous space-filling curve over the entire sphere. The 32 // transformation is designed to make the cells at each level fairly uniform 33 // in size. 34 // 35 ////////////////////////// S2 Cell Decomposition ///////////////////////// 36 // 37 // The following methods define the cube-to-sphere projection used by 38 // the Cell decomposition. 39 // 40 // In the process of converting a latitude-longitude pair to a 64-bit cell 41 // id, the following coordinate systems are used: 42 // 43 // (id) 44 // An CellID is a 64-bit encoding of a face and a Hilbert curve position 45 // on that face. The Hilbert curve position implicitly encodes both the 46 // position of a cell and its subdivision level (see s2cellid.go). 47 // 48 // (face, i, j) 49 // Leaf-cell coordinates. "i" and "j" are integers in the range 50 // [0,(2**30)-1] that identify a particular leaf cell on the given face. 51 // The (i, j) coordinate system is right-handed on each face, and the 52 // faces are oriented such that Hilbert curves connect continuously from 53 // one face to the next. 54 // 55 // (face, s, t) 56 // Cell-space coordinates. "s" and "t" are real numbers in the range 57 // [0,1] that identify a point on the given face. For example, the point 58 // (s, t) = (0.5, 0.5) corresponds to the center of the top-level face 59 // cell. This point is also a vertex of exactly four cells at each 60 // subdivision level greater than zero. 61 // 62 // (face, si, ti) 63 // Discrete cell-space coordinates. These are obtained by multiplying 64 // "s" and "t" by 2**31 and rounding to the nearest unsigned integer. 65 // Discrete coordinates lie in the range [0,2**31]. This coordinate 66 // system can represent the edge and center positions of all cells with 67 // no loss of precision (including non-leaf cells). In binary, each 68 // coordinate of a level-k cell center ends with a 1 followed by 69 // (30 - k) 0s. The coordinates of its edges end with (at least) 70 // (31 - k) 0s. 71 // 72 // (face, u, v) 73 // Cube-space coordinates in the range [-1,1]. To make the cells at each 74 // level more uniform in size after they are projected onto the sphere, 75 // we apply a nonlinear transformation of the form u=f(s), v=f(t). 76 // The (u, v) coordinates after this transformation give the actual 77 // coordinates on the cube face (modulo some 90 degree rotations) before 78 // it is projected onto the unit sphere. 79 // 80 // (face, u, v, w) 81 // Per-face coordinate frame. This is an extension of the (face, u, v) 82 // cube-space coordinates that adds a third axis "w" in the direction of 83 // the face normal. It is always a right-handed 3D coordinate system. 84 // Cube-space coordinates can be converted to this frame by setting w=1, 85 // while (u,v,w) coordinates can be projected onto the cube face by 86 // dividing by w, i.e. (face, u/w, v/w). 87 // 88 // (x, y, z) 89 // Direction vector (Point). Direction vectors are not necessarily unit 90 // length, and are often chosen to be points on the biunit cube 91 // [-1,+1]x[-1,+1]x[-1,+1]. They can be be normalized to obtain the 92 // corresponding point on the unit sphere. 93 // 94 // (lat, lng) 95 // Latitude and longitude (LatLng). Latitudes must be between -90 and 96 // 90 degrees inclusive, and longitudes must be between -180 and 180 97 // degrees inclusive. 98 // 99 // Note that the (i, j), (s, t), (si, ti), and (u, v) coordinate systems are 100 // right-handed on all six faces. 101 // 102 // 103 // There are a number of different projections from cell-space (s,t) to 104 // cube-space (u,v): linear, quadratic, and tangent. They have the following 105 // tradeoffs: 106 // 107 // Linear - This is the fastest transformation, but also produces the least 108 // uniform cell sizes. Cell areas vary by a factor of about 5.2, with the 109 // largest cells at the center of each face and the smallest cells in 110 // the corners. 111 // 112 // Tangent - Transforming the coordinates via Atan makes the cell sizes 113 // more uniform. The areas vary by a maximum ratio of 1.4 as opposed to a 114 // maximum ratio of 5.2. However, each call to Atan is about as expensive 115 // as all of the other calculations combined when converting from points to 116 // cell ids, i.e. it reduces performance by a factor of 3. 117 // 118 // Quadratic - This is an approximation of the tangent projection that 119 // is much faster and produces cells that are almost as uniform in size. 120 // It is about 3 times faster than the tangent projection for converting 121 // cell ids to points or vice versa. Cell areas vary by a maximum ratio of 122 // about 2.1. 123 // 124 // Here is a table comparing the cell uniformity using each projection. Area 125 // Ratio is the maximum ratio over all subdivision levels of the largest cell 126 // area to the smallest cell area at that level, Edge Ratio is the maximum 127 // ratio of the longest edge of any cell to the shortest edge of any cell at 128 // the same level, and Diag Ratio is the ratio of the longest diagonal of 129 // any cell to the shortest diagonal of any cell at the same level. 130 // 131 // Area Edge Diag 132 // Ratio Ratio Ratio 133 // ----------------------------------- 134 // Linear: 5.200 2.117 2.959 135 // Tangent: 1.414 1.414 1.704 136 // Quadratic: 2.082 1.802 1.932 137 // 138 // The worst-case cell aspect ratios are about the same with all three 139 // projections. The maximum ratio of the longest edge to the shortest edge 140 // within the same cell is about 1.4 and the maximum ratio of the diagonals 141 // within the same cell is about 1.7. 142 // 143 // For Go we have chosen to use only the Quadratic approach. Other language 144 // implementations may offer other choices. 145 146 const ( 147 // maxSiTi is the maximum value of an si- or ti-coordinate. 148 // It is one shift more than MaxSize. The range of valid (si,ti) 149 // values is [0..maxSiTi]. 150 maxSiTi = MaxSize << 1 151 ) 152 153 // siTiToST converts an si- or ti-value to the corresponding s- or t-value. 154 // Value is capped at 1.0 because there is no DCHECK in Go. 155 func siTiToST(si uint32) float64 { 156 if si > maxSiTi { 157 return 1.0 158 } 159 return float64(si) / float64(maxSiTi) 160 } 161 162 // stToSiTi converts the s- or t-value to the nearest si- or ti-coordinate. 163 // The result may be outside the range of valid (si,ti)-values. Value of 164 // 0.49999999999999994 (math.NextAfter(0.5, -1)), will be incorrectly rounded up. 165 func stToSiTi(s float64) uint32 { 166 if s < 0 { 167 return uint32(s*maxSiTi - 0.5) 168 } 169 return uint32(s*maxSiTi + 0.5) 170 } 171 172 // stToUV converts an s or t value to the corresponding u or v value. 173 // This is a non-linear transformation from [-1,1] to [-1,1] that 174 // attempts to make the cell sizes more uniform. 175 // This uses what the C++ version calls 'the quadratic transform'. 176 func stToUV(s float64) float64 { 177 if s >= 0.5 { 178 return (1 / 3.) * (4*s*s - 1) 179 } 180 return (1 / 3.) * (1 - 4*(1-s)*(1-s)) 181 } 182 183 // uvToST is the inverse of the stToUV transformation. Note that it 184 // is not always true that uvToST(stToUV(x)) == x due to numerical 185 // errors. 186 func uvToST(u float64) float64 { 187 if u >= 0 { 188 return 0.5 * math.Sqrt(1+3*u) 189 } 190 return 1 - 0.5*math.Sqrt(1-3*u) 191 } 192 193 // face returns face ID from 0 to 5 containing the r. For points on the 194 // boundary between faces, the result is arbitrary but deterministic. 195 func face(r r3.Vector) int { 196 f := r.LargestComponent() 197 switch { 198 case f == r3.XAxis && r.X < 0: 199 f += 3 200 case f == r3.YAxis && r.Y < 0: 201 f += 3 202 case f == r3.ZAxis && r.Z < 0: 203 f += 3 204 } 205 return int(f) 206 } 207 208 // validFaceXYZToUV given a valid face for the given point r (meaning that 209 // dot product of r with the face normal is positive), returns 210 // the corresponding u and v values, which may lie outside the range [-1,1]. 211 func validFaceXYZToUV(face int, r r3.Vector) (float64, float64) { 212 switch face { 213 case 0: 214 return r.Y / r.X, r.Z / r.X 215 case 1: 216 return -r.X / r.Y, r.Z / r.Y 217 case 2: 218 return -r.X / r.Z, -r.Y / r.Z 219 case 3: 220 return r.Z / r.X, r.Y / r.X 221 case 4: 222 return r.Z / r.Y, -r.X / r.Y 223 } 224 return -r.Y / r.Z, -r.X / r.Z 225 } 226 227 // xyzToFaceUV converts a direction vector (not necessarily unit length) to 228 // (face, u, v) coordinates. 229 func xyzToFaceUV(r r3.Vector) (f int, u, v float64) { 230 f = face(r) 231 u, v = validFaceXYZToUV(f, r) 232 return f, u, v 233 } 234 235 // faceUVToXYZ turns face and UV coordinates into an unnormalized 3 vector. 236 func faceUVToXYZ(face int, u, v float64) r3.Vector { 237 switch face { 238 case 0: 239 return r3.Vector{1, u, v} 240 case 1: 241 return r3.Vector{-u, 1, v} 242 case 2: 243 return r3.Vector{-u, -v, 1} 244 case 3: 245 return r3.Vector{-1, -v, -u} 246 case 4: 247 return r3.Vector{v, -1, -u} 248 default: 249 return r3.Vector{v, u, -1} 250 } 251 } 252 253 // faceXYZToUV returns the u and v values (which may lie outside the range 254 // [-1, 1]) if the dot product of the point p with the given face normal is positive. 255 func faceXYZToUV(face int, p Point) (u, v float64, ok bool) { 256 switch face { 257 case 0: 258 if p.X <= 0 { 259 return 0, 0, false 260 } 261 case 1: 262 if p.Y <= 0 { 263 return 0, 0, false 264 } 265 case 2: 266 if p.Z <= 0 { 267 return 0, 0, false 268 } 269 case 3: 270 if p.X >= 0 { 271 return 0, 0, false 272 } 273 case 4: 274 if p.Y >= 0 { 275 return 0, 0, false 276 } 277 default: 278 if p.Z >= 0 { 279 return 0, 0, false 280 } 281 } 282 283 u, v = validFaceXYZToUV(face, p.Vector) 284 return u, v, true 285 } 286 287 // faceXYZtoUVW transforms the given point P to the (u,v,w) coordinate frame of the given 288 // face where the w-axis represents the face normal. 289 func faceXYZtoUVW(face int, p Point) Point { 290 // The result coordinates are simply the dot products of P with the (u,v,w) 291 // axes for the given face (see faceUVWAxes). 292 switch face { 293 case 0: 294 return Point{r3.Vector{p.Y, p.Z, p.X}} 295 case 1: 296 return Point{r3.Vector{-p.X, p.Z, p.Y}} 297 case 2: 298 return Point{r3.Vector{-p.X, -p.Y, p.Z}} 299 case 3: 300 return Point{r3.Vector{-p.Z, -p.Y, -p.X}} 301 case 4: 302 return Point{r3.Vector{-p.Z, p.X, -p.Y}} 303 default: 304 return Point{r3.Vector{p.Y, p.X, -p.Z}} 305 } 306 } 307 308 // faceSiTiToXYZ transforms the (si, ti) coordinates to a (not necessarily 309 // unit length) Point on the given face. 310 func faceSiTiToXYZ(face int, si, ti uint32) Point { 311 return Point{faceUVToXYZ(face, stToUV(siTiToST(si)), stToUV(siTiToST(ti)))} 312 } 313 314 // xyzToFaceSiTi transforms the (not necessarily unit length) Point to 315 // (face, si, ti) coordinates and the level the Point is at. 316 func xyzToFaceSiTi(p Point) (face int, si, ti uint32, level int) { 317 face, u, v := xyzToFaceUV(p.Vector) 318 si = stToSiTi(uvToST(u)) 319 ti = stToSiTi(uvToST(v)) 320 321 // If the levels corresponding to si,ti are not equal, then p is not a cell 322 // center. The si,ti values of 0 and maxSiTi need to be handled specially 323 // because they do not correspond to cell centers at any valid level; they 324 // are mapped to level -1 by the code at the end. 325 level = MaxLevel - findLSBSetNonZero64(uint64(si|maxSiTi)) 326 if level < 0 || level != MaxLevel-findLSBSetNonZero64(uint64(ti|maxSiTi)) { 327 return face, si, ti, -1 328 } 329 330 // In infinite precision, this test could be changed to ST == SiTi. However, 331 // due to rounding errors, uvToST(xyzToFaceUV(faceUVToXYZ(stToUV(...)))) is 332 // not idempotent. On the other hand, the center is computed exactly the same 333 // way p was originally computed (if it is indeed the center of a Cell); 334 // the comparison can be exact. 335 if p.Vector == faceSiTiToXYZ(face, si, ti).Normalize() { 336 return face, si, ti, level 337 } 338 339 return face, si, ti, -1 340 } 341 342 // uNorm returns the right-handed normal (not necessarily unit length) for an 343 // edge in the direction of the positive v-axis at the given u-value on 344 // the given face. (This vector is perpendicular to the plane through 345 // the sphere origin that contains the given edge.) 346 func uNorm(face int, u float64) r3.Vector { 347 switch face { 348 case 0: 349 return r3.Vector{u, -1, 0} 350 case 1: 351 return r3.Vector{1, u, 0} 352 case 2: 353 return r3.Vector{1, 0, u} 354 case 3: 355 return r3.Vector{-u, 0, 1} 356 case 4: 357 return r3.Vector{0, -u, 1} 358 default: 359 return r3.Vector{0, -1, -u} 360 } 361 } 362 363 // vNorm returns the right-handed normal (not necessarily unit length) for an 364 // edge in the direction of the positive u-axis at the given v-value on 365 // the given face. 366 func vNorm(face int, v float64) r3.Vector { 367 switch face { 368 case 0: 369 return r3.Vector{-v, 0, 1} 370 case 1: 371 return r3.Vector{0, -v, 1} 372 case 2: 373 return r3.Vector{0, -1, -v} 374 case 3: 375 return r3.Vector{v, -1, 0} 376 case 4: 377 return r3.Vector{1, v, 0} 378 default: 379 return r3.Vector{1, 0, v} 380 } 381 } 382 383 // faceUVWAxes are the U, V, and W axes for each face. 384 var faceUVWAxes = [6][3]Point{ 385 {Point{r3.Vector{0, 1, 0}}, Point{r3.Vector{0, 0, 1}}, Point{r3.Vector{1, 0, 0}}}, 386 {Point{r3.Vector{-1, 0, 0}}, Point{r3.Vector{0, 0, 1}}, Point{r3.Vector{0, 1, 0}}}, 387 {Point{r3.Vector{-1, 0, 0}}, Point{r3.Vector{0, -1, 0}}, Point{r3.Vector{0, 0, 1}}}, 388 {Point{r3.Vector{0, 0, -1}}, Point{r3.Vector{0, -1, 0}}, Point{r3.Vector{-1, 0, 0}}}, 389 {Point{r3.Vector{0, 0, -1}}, Point{r3.Vector{1, 0, 0}}, Point{r3.Vector{0, -1, 0}}}, 390 {Point{r3.Vector{0, 1, 0}}, Point{r3.Vector{1, 0, 0}}, Point{r3.Vector{0, 0, -1}}}, 391 } 392 393 // faceUVWFaces are the precomputed neighbors of each face. 394 var faceUVWFaces = [6][3][2]int{ 395 {{4, 1}, {5, 2}, {3, 0}}, 396 {{0, 3}, {5, 2}, {4, 1}}, 397 {{0, 3}, {1, 4}, {5, 2}}, 398 {{2, 5}, {1, 4}, {0, 3}}, 399 {{2, 5}, {3, 0}, {1, 4}}, 400 {{4, 1}, {3, 0}, {2, 5}}, 401 } 402 403 // uvwAxis returns the given axis of the given face. 404 func uvwAxis(face, axis int) Point { 405 return faceUVWAxes[face][axis] 406 } 407 408 // uvwFaces returns the face in the (u,v,w) coordinate system on the given axis 409 // in the given direction. 410 func uvwFace(face, axis, direction int) int { 411 return faceUVWFaces[face][axis][direction] 412 } 413 414 // uAxis returns the u-axis for the given face. 415 func uAxis(face int) Point { 416 return uvwAxis(face, 0) 417 } 418 419 // vAxis returns the v-axis for the given face. 420 func vAxis(face int) Point { 421 return uvwAxis(face, 1) 422 } 423 424 // Return the unit-length normal for the given face. 425 func unitNorm(face int) Point { 426 return uvwAxis(face, 2) 427 } 428