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 import com.google.common.geometry.S2.Metric; 19 20 /** 21 * This class specifies the details of how the cube faces are projected onto the 22 * unit sphere. This includes getting the face ordering and orientation correct 23 * so that sequentially increasing cell ids follow a continuous space-filling 24 * curve over the entire sphere, and defining the transformation from cell-space 25 * to cube-space (see s2.h) in order to make the cells more uniform in size. 26 * 27 * 28 * We have implemented three different projections from cell-space (s,t) to 29 * cube-space (u,v): linear, quadratic, and tangent. They have the following 30 * tradeoffs: 31 * 32 * Linear - This is the fastest transformation, but also produces the least 33 * uniform cell sizes. Cell areas vary by a factor of about 5.2, with the 34 * largest cells at the center of each face and the smallest cells in the 35 * corners. 36 * 37 * Tangent - Transforming the coordinates via atan() makes the cell sizes more 38 * uniform. The areas vary by a maximum ratio of 1.4 as opposed to a maximum 39 * ratio of 5.2. However, each call to atan() is about as expensive as all of 40 * the other calculations combined when converting from points to cell ids, i.e. 41 * it reduces performance by a factor of 3. 42 * 43 * Quadratic - This is an approximation of the tangent projection that is much 44 * faster and produces cells that are almost as uniform in size. It is about 3 45 * times faster than the tangent projection for converting cell ids to points, 46 * and 2 times faster for converting points to cell ids. Cell areas vary by a 47 * maximum ratio of about 2.1. 48 * 49 * Here is a table comparing the cell uniformity using each projection. "Area 50 * ratio" is the maximum ratio over all subdivision levels of the largest cell 51 * area to the smallest cell area at that level, "edge ratio" is the maximum 52 * ratio of the longest edge of any cell to the shortest edge of any cell at the 53 * same level, and "diag ratio" is the ratio of the longest diagonal of any cell 54 * to the shortest diagonal of any cell at the same level. "ToPoint" and 55 * "FromPoint" are the times in microseconds required to convert cell ids to and 56 * from points (unit vectors) respectively. 57 * 58 * Area Edge Diag ToPoint FromPoint Ratio Ratio Ratio (microseconds) 59 * ------------------------------------------------------- Linear: 5.200 2.117 60 * 2.959 0.103 0.123 Tangent: 1.414 1.414 1.704 0.290 0.306 Quadratic: 2.082 61 * 1.802 1.932 0.116 0.161 62 * 63 * The worst-case cell aspect ratios are about the same with all three 64 * projections. The maximum ratio of the longest edge to the shortest edge 65 * within the same cell is about 1.4 and the maximum ratio of the diagonals 66 * within the same cell is about 1.7. 67 * 68 * This data was produced using s2cell_unittest and s2cellid_unittest. 69 * 70 */ 71 72 public final strictfp class S2Projections { 73 public enum Projections { 74 S2_LINEAR_PROJECTION, S2_TAN_PROJECTION, S2_QUADRATIC_PROJECTION 75 } 76 77 private static final Projections S2_PROJECTION = Projections.S2_QUADRATIC_PROJECTION; 78 79 // All of the values below were obtained by a combination of hand analysis and 80 // Mathematica. In general, S2_TAN_PROJECTION produces the most uniform 81 // shapes and sizes of cells, S2_LINEAR_PROJECTION is considerably worse, and 82 // S2_QUADRATIC_PROJECTION is somewhere in between (but generally closer to 83 // the tangent projection than the linear one). 84 85 86 // The minimum area of any cell at level k is at least MIN_AREA.GetValue(k), 87 // and the maximum is at most MAX_AREA.GetValue(k). The average area of all 88 // cells at level k is exactly AVG_AREA.GetValue(k). 89 public static final Metric MIN_AREA = new Metric(2, 90 S2_PROJECTION == Projections.S2_LINEAR_PROJECTION ? 1 / (3 * Math.sqrt(3)) : // 0.192 91 S2_PROJECTION == Projections.S2_TAN_PROJECTION ? (S2.M_PI * S2.M_PI) 92 / (16 * S2.M_SQRT2) : // 0.436 93 S2_PROJECTION == Projections.S2_QUADRATIC_PROJECTION 94 ? 2 * S2.M_SQRT2 / 9 : // 0.314 95 0); 96 public static final Metric MAX_AREA = new Metric(2, 97 S2_PROJECTION == Projections.S2_LINEAR_PROJECTION ? 1 : // 1.000 98 S2_PROJECTION == Projections.S2_TAN_PROJECTION ? S2.M_PI * S2.M_PI / 16 : // 0.617 99 S2_PROJECTION == Projections.S2_QUADRATIC_PROJECTION 100 ? 0.65894981424079037 : // 0.659 101 0); 102 public static final Metric AVG_AREA = new Metric(2, S2.M_PI / 6); // 0.524) 103 104 105 // Each cell is bounded by four planes passing through its four edges and 106 // the center of the sphere. These metrics relate to the angle between each 107 // pair of opposite bounding planes, or equivalently, between the planes 108 // corresponding to two different s-values or two different t-values. For 109 // example, the maximum angle between opposite bounding planes for a cell at 110 // level k is MAX_ANGLE_SPAN.GetValue(k), and the average angle span for all 111 // cells at level k is approximately AVG_ANGLE_SPAN.GetValue(k). 112 public static final Metric MIN_ANGLE_SPAN = new Metric(1, 113 S2_PROJECTION == Projections.S2_LINEAR_PROJECTION ? 0.5 : // 0.500 114 S2_PROJECTION == Projections.S2_TAN_PROJECTION ? S2.M_PI / 4 : // 0.785 115 S2_PROJECTION == Projections.S2_QUADRATIC_PROJECTION ? 2. / 3 : // 0.667 116 0); 117 public static final Metric MAX_ANGLE_SPAN = new Metric(1, 118 S2_PROJECTION == Projections.S2_LINEAR_PROJECTION ? 1 : // 1.000 119 S2_PROJECTION == Projections.S2_TAN_PROJECTION ? S2.M_PI / 4 : // 0.785 120 S2_PROJECTION == Projections.S2_QUADRATIC_PROJECTION 121 ? 0.85244858959960922 : // 0.852 122 0); 123 public static final Metric AVG_ANGLE_SPAN = new Metric(1, S2.M_PI / 4); // 0.785 124 125 126 // The width of geometric figure is defined as the distance between two 127 // parallel bounding lines in a given direction. For cells, the minimum 128 // width is always attained between two opposite edges, and the maximum 129 // width is attained between two opposite vertices. However, for our 130 // purposes we redefine the width of a cell as the perpendicular distance 131 // between a pair of opposite edges. A cell therefore has two widths, one 132 // in each direction. The minimum width according to this definition agrees 133 // with the classic geometric one, but the maximum width is different. (The 134 // maximum geometric width corresponds to MAX_DIAG defined below.) 135 // 136 // For a cell at level k, the distance between opposite edges is at least 137 // MIN_WIDTH.GetValue(k) and at most MAX_WIDTH.GetValue(k). The average 138 // width in both directions for all cells at level k is approximately 139 // AVG_WIDTH.GetValue(k). 140 // 141 // The width is useful for bounding the minimum or maximum distance from a 142 // point on one edge of a cell to the closest point on the opposite edge. 143 // For example, this is useful when "growing" regions by a fixed distance. 144 public static final Metric MIN_WIDTH = new Metric(1, 145 (S2Projections.S2_PROJECTION == Projections.S2_LINEAR_PROJECTION ? 1 / Math.sqrt(6) : // 0.408 146 S2_PROJECTION == Projections.S2_TAN_PROJECTION ? S2.M_PI / (4 * S2.M_SQRT2) : // 0.555 147 S2_PROJECTION == Projections.S2_QUADRATIC_PROJECTION ? S2.M_SQRT2 / 3 : // 0.471 148 0)); 149 150 public static final Metric MAX_WIDTH = new Metric(1, MAX_ANGLE_SPAN.deriv()); 151 public static final Metric AVG_WIDTH = new Metric(1, 152 S2_PROJECTION == Projections.S2_LINEAR_PROJECTION ? 0.70572967292222848 : // 0.706 153 S2_PROJECTION == Projections.S2_TAN_PROJECTION ? 0.71865931946258044 : // 0.719 154 S2_PROJECTION == Projections.S2_QUADRATIC_PROJECTION 155 ? 0.71726183644304969 : // 0.717 156 0); 157 158 // The minimum edge length of any cell at level k is at least 159 // MIN_EDGE.GetValue(k), and the maximum is at most MAX_EDGE.GetValue(k). 160 // The average edge length is approximately AVG_EDGE.GetValue(k). 161 // 162 // The edge length metrics can also be used to bound the minimum, maximum, 163 // or average distance from the center of one cell to the center of one of 164 // its edge neighbors. In particular, it can be used to bound the distance 165 // between adjacent cell centers along the space-filling Hilbert curve for 166 // cells at any given level. 167 public static final Metric MIN_EDGE = new Metric(1, 168 S2_PROJECTION == Projections.S2_LINEAR_PROJECTION ? S2.M_SQRT2 / 3 : // 0.471 169 S2_PROJECTION == Projections.S2_TAN_PROJECTION ? S2.M_PI / (4 * S2.M_SQRT2) : // 0.555 170 S2_PROJECTION == Projections.S2_QUADRATIC_PROJECTION ? S2.M_SQRT2 / 3 : // 0.471 171 0); 172 public static final Metric MAX_EDGE = new Metric(1, MAX_ANGLE_SPAN.deriv()); 173 public static final Metric AVG_EDGE = new Metric(1, 174 S2_PROJECTION == Projections.S2_LINEAR_PROJECTION ? 0.72001709647780182 : // 0.720 175 S2_PROJECTION == Projections.S2_TAN_PROJECTION ? 0.73083351627336963 : // 0.731 176 S2_PROJECTION == Projections.S2_QUADRATIC_PROJECTION 177 ? 0.72960687319305303 : // 0.730 178 0); 179 180 181 // The minimum diagonal length of any cell at level k is at least 182 // MIN_DIAG.GetValue(k), and the maximum is at most MAX_DIAG.GetValue(k). 183 // The average diagonal length is approximately AVG_DIAG.GetValue(k). 184 // 185 // The maximum diagonal also happens to be the maximum diameter of any cell, 186 // and also the maximum geometric width (see the discussion above). So for 187 // example, the distance from an arbitrary point to the closest cell center 188 // at a given level is at most half the maximum diagonal length. 189 public static final Metric MIN_DIAG = new Metric(1, 190 S2_PROJECTION == Projections.S2_LINEAR_PROJECTION ? S2.M_SQRT2 / 3 : // 0.471 191 S2_PROJECTION == Projections.S2_TAN_PROJECTION ? S2.M_PI / (3 * S2.M_SQRT2) : // 0.740 192 S2_PROJECTION == Projections.S2_QUADRATIC_PROJECTION 193 ? 4 * S2.M_SQRT2 / 9 : // 0.629 194 0); 195 public static final Metric MAX_DIAG = new Metric(1, 196 S2_PROJECTION == Projections.S2_LINEAR_PROJECTION ? S2.M_SQRT2 : // 1.414 197 S2_PROJECTION == Projections.S2_TAN_PROJECTION ? S2.M_PI / Math.sqrt(6) : // 1.283 198 S2_PROJECTION == Projections.S2_QUADRATIC_PROJECTION 199 ? 1.2193272972170106 : // 1.219 200 0); 201 public static final Metric AVG_DIAG = new Metric(1, 202 S2_PROJECTION == Projections.S2_LINEAR_PROJECTION ? 1.0159089332094063 : // 1.016 203 S2_PROJECTION == Projections.S2_TAN_PROJECTION ? 1.0318115985978178 : // 1.032 204 S2_PROJECTION == Projections.S2_QUADRATIC_PROJECTION 205 ? 1.03021136949923584 : // 1.030 206 0); 207 208 // This is the maximum edge aspect ratio over all cells at any level, where 209 // the edge aspect ratio of a cell is defined as the ratio of its longest 210 // edge length to its shortest edge length. 211 public static final double MAX_EDGE_ASPECT = 212 S2_PROJECTION == Projections.S2_LINEAR_PROJECTION ? S2.M_SQRT2 : // 1.414 213 S2_PROJECTION == Projections.S2_TAN_PROJECTION ? S2.M_SQRT2 : // 1.414 214 S2_PROJECTION == Projections.S2_QUADRATIC_PROJECTION ? 1.44261527445268292 : // 1.443 215 0; 216 217 // This is the maximum diagonal aspect ratio over all cells at any level, 218 // where the diagonal aspect ratio of a cell is defined as the ratio of its 219 // longest diagonal length to its shortest diagonal length. 220 public static final double MAX_DIAG_ASPECT = Math.sqrt(3); // 1.732 221 stToUV(double s)222 public static double stToUV(double s) { 223 switch (S2_PROJECTION) { 224 case S2_LINEAR_PROJECTION: 225 return s; 226 case S2_TAN_PROJECTION: 227 // Unfortunately, tan(M_PI_4) is slightly less than 1.0. This isn't due 228 // to 229 // a flaw in the implementation of tan(), it's because the derivative of 230 // tan(x) at x=pi/4 is 2, and it happens that the two adjacent floating 231 // point numbers on either side of the infinite-precision value of pi/4 232 // have 233 // tangents that are slightly below and slightly above 1.0 when rounded 234 // to 235 // the nearest double-precision result. 236 s = Math.tan(S2.M_PI_4 * s); 237 return s + (1.0 / (1L << 53)) * s; 238 case S2_QUADRATIC_PROJECTION: 239 if (s >= 0) { 240 return (1 / 3.) * ((1 + s) * (1 + s) - 1); 241 } else { 242 return (1 / 3.) * (1 - (1 - s) * (1 - s)); 243 } 244 default: 245 throw new IllegalStateException("Invalid value for S2_PROJECTION"); 246 } 247 } 248 uvToST(double u)249 public static double uvToST(double u) { 250 switch (S2_PROJECTION) { 251 case S2_LINEAR_PROJECTION: 252 return u; 253 case S2_TAN_PROJECTION: 254 return (4 * S2.M_1_PI) * Math.atan(u); 255 case S2_QUADRATIC_PROJECTION: 256 if (u >= 0) { 257 return Math.sqrt(1 + 3 * u) - 1; 258 } else { 259 return 1 - Math.sqrt(1 - 3 * u); 260 } 261 default: 262 throw new IllegalStateException("Invalid value for S2_PROJECTION"); 263 } 264 } 265 266 267 /** 268 * Convert (face, u, v) coordinates to a direction vector (not necessarily 269 * unit length). 270 */ faceUvToXyz(int face, double u, double v)271 public static S2Point faceUvToXyz(int face, double u, double v) { 272 switch (face) { 273 case 0: 274 return new S2Point(1, u, v); 275 case 1: 276 return new S2Point(-u, 1, v); 277 case 2: 278 return new S2Point(-u, -v, 1); 279 case 3: 280 return new S2Point(-1, -v, -u); 281 case 4: 282 return new S2Point(v, -1, -u); 283 default: 284 return new S2Point(v, u, -1); 285 } 286 } 287 validFaceXyzToUv(int face, S2Point p)288 public static R2Vector validFaceXyzToUv(int face, S2Point p) { 289 // assert (p.dotProd(faceUvToXyz(face, 0, 0)) > 0); 290 double pu; 291 double pv; 292 switch (face) { 293 case 0: 294 pu = p.y / p.x; 295 pv = p.z / p.x; 296 break; 297 case 1: 298 pu = -p.x / p.y; 299 pv = p.z / p.y; 300 break; 301 case 2: 302 pu = -p.x / p.z; 303 pv = -p.y / p.z; 304 break; 305 case 3: 306 pu = p.z / p.x; 307 pv = p.y / p.x; 308 break; 309 case 4: 310 pu = p.z / p.y; 311 pv = -p.x / p.y; 312 break; 313 default: 314 pu = -p.y / p.z; 315 pv = -p.x / p.z; 316 break; 317 } 318 return new R2Vector(pu, pv); 319 } 320 xyzToFace(S2Point p)321 public static int xyzToFace(S2Point p) { 322 int face = p.largestAbsComponent(); 323 if (p.get(face) < 0) { 324 face += 3; 325 } 326 return face; 327 } 328 faceXyzToUv(int face, S2Point p)329 public static R2Vector faceXyzToUv(int face, S2Point p) { 330 if (face < 3) { 331 if (p.get(face) <= 0) { 332 return null; 333 } 334 } else { 335 if (p.get(face - 3) >= 0) { 336 return null; 337 } 338 } 339 return validFaceXyzToUv(face, p); 340 } 341 getUNorm(int face, double u)342 public static S2Point getUNorm(int face, double u) { 343 switch (face) { 344 case 0: 345 return new S2Point(u, -1, 0); 346 case 1: 347 return new S2Point(1, u, 0); 348 case 2: 349 return new S2Point(1, 0, u); 350 case 3: 351 return new S2Point(-u, 0, 1); 352 case 4: 353 return new S2Point(0, -u, 1); 354 default: 355 return new S2Point(0, -1, -u); 356 } 357 } 358 getVNorm(int face, double v)359 public static S2Point getVNorm(int face, double v) { 360 switch (face) { 361 case 0: 362 return new S2Point(-v, 0, 1); 363 case 1: 364 return new S2Point(0, -v, 1); 365 case 2: 366 return new S2Point(0, -1, -v); 367 case 3: 368 return new S2Point(v, -1, 0); 369 case 4: 370 return new S2Point(1, v, 0); 371 default: 372 return new S2Point(1, 0, v); 373 } 374 } 375 getNorm(int face)376 public static S2Point getNorm(int face) { 377 return faceUvToXyz(face, 0, 0); 378 } 379 getUAxis(int face)380 public static S2Point getUAxis(int face) { 381 switch (face) { 382 case 0: 383 return new S2Point(0, 1, 0); 384 case 1: 385 return new S2Point(-1, 0, 0); 386 case 2: 387 return new S2Point(-1, 0, 0); 388 case 3: 389 return new S2Point(0, 0, -1); 390 case 4: 391 return new S2Point(0, 0, -1); 392 default: 393 return new S2Point(0, 1, 0); 394 } 395 } 396 getVAxis(int face)397 public static S2Point getVAxis(int face) { 398 switch (face) { 399 case 0: 400 return new S2Point(0, 0, 1); 401 case 1: 402 return new S2Point(0, 0, 1); 403 case 2: 404 return new S2Point(0, -1, 0); 405 case 3: 406 return new S2Point(0, -1, 0); 407 case 4: 408 return new S2Point(1, 0, 0); 409 default: 410 return new S2Point(1, 0, 0); 411 } 412 } 413 414 // Don't instantiate S2Projections()415 private S2Projections() { 416 } 417 } 418