blob: 5b0961bbbdd7875f0fc253f7024084b66bd1afab [file] [log] [blame]
 // Copyright 2005 Google Inc. All Rights Reserved. // // 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. // // Author: ericv@google.com (Eric Veach) // // This file contains documentation of the various coordinate systems used // throughout the library. Most importantly, S2 defines a framework for // decomposing the unit sphere into a hierarchy of "cells". Each cell is a // quadrilateral bounded by four geodesics. The top level of the hierarchy is // obtained by projecting the six faces of a cube onto the unit sphere, and // lower levels are obtained by subdividing each cell into four children // recursively. Cells are numbered such that sequentially increasing cells // follow a continuous space-filling curve over the entire sphere. The // transformation is designed to make the cells at each level fairly uniform // in size. // // ////////////////////////// S2Cell Decomposition ///////////////////////// // // The following methods define the cube-to-sphere projection used by // the S2Cell decomposition. // // In the process of converting a latitude-longitude pair to a 64-bit cell // id, the following coordinate systems are used: // // (id) // An S2CellId is a 64-bit encoding of a face and a Hilbert curve position // on that face. The Hilbert curve position implicitly encodes both the // position of a cell and its subdivision level (see s2cellid.h). // // (face, i, j) // Leaf-cell coordinates. "i" and "j" are integers in the range // [0,(2**30)-1] that identify a particular leaf cell on the given face. // The (i, j) coordinate system is right-handed on each face, and the // faces are oriented such that Hilbert curves connect continuously from // one face to the next. // // (face, s, t) // Cell-space coordinates. "s" and "t" are real numbers in the range // [0,1] that identify a point on the given face. For example, the point // (s, t) = (0.5, 0.5) corresponds to the center of the top-level face // cell. This point is also a vertex of exactly four cells at each // subdivision level greater than zero. // // (face, si, ti) // Discrete cell-space coordinates. These are obtained by multiplying // "s" and "t" by 2**31 and rounding to the nearest unsigned integer. // Discrete coordinates lie in the range [0,2**31]. This coordinate // system can represent the edge and center positions of all cells with // no loss of precision (including non-leaf cells). In binary, each // coordinate of a level-k cell center ends with a 1 followed by // (30 - k) 0s. The coordinates of its edges end with (at least) // (31 - k) 0s. // // (face, u, v) // Cube-space coordinates in the range [-1,1]. To make the cells at each // level more uniform in size after they are projected onto the sphere, // we apply a nonlinear transformation of the form u=f(s), v=f(t). // The (u, v) coordinates after this transformation give the actual // coordinates on the cube face (modulo some 90 degree rotations) before // it is projected onto the unit sphere. // // (face, u, v, w) // Per-face coordinate frame. This is an extension of the (face, u, v) // cube-space coordinates that adds a third axis "w" in the direction of // the face normal. It is always a right-handed 3D coordinate system. // Cube-space coordinates can be converted to this frame by setting w=1, // while (u,v,w) coordinates can be projected onto the cube face by // dividing by w, i.e. (face, u/w, v/w). // // (x, y, z) // Direction vector (S2Point). Direction vectors are not necessarily unit // length, and are often chosen to be points on the biunit cube // [-1,+1]x[-1,+1]x[-1,+1]. They can be be normalized to obtain the // corresponding point on the unit sphere. // // (lat, lng) // Latitude and longitude (S2LatLng). Latitudes must be between -90 and // 90 degrees inclusive, and longitudes must be between -180 and 180 // degrees inclusive. // // Note that the (i, j), (s, t), (si, ti), and (u, v) coordinate systems are // right-handed on all six faces. #ifndef S2_S2COORDS_H_ #define S2_S2COORDS_H_ #include #include #include "base/logging.h" #include "s2/r2.h" #include "s2/s2coords-internal.h" #include "s2/s2point.h" #include "s2/util/math/mathutil.h" // S2 is a namespace for constants and simple utility functions that are used // throughout the S2 library. The name "S2" is derived from the mathematical // symbol for the two-dimensional unit sphere (note that the "2" refers to the // dimension of the surface, not the space it is embedded in). namespace S2 { // This is the number of levels needed to specify a leaf cell. This // constant is defined here so that the S2::Metric class and the conversion // functions below can be implemented without including s2cellid.h. Please // see s2cellid.h for other useful constants and conversion functions. int const kMaxCellLevel = 30; // The maximum index of a valid leaf cell plus one. The range of valid leaf // cell indices is [0..kLimitIJ-1]. int const kLimitIJ = 1 << kMaxCellLevel; // == S2CellId::kMaxSize // The maximum value of an si- or ti-coordinate. The range of valid (si,ti) // values is [0..kMaxSiTi]. unsigned int const kMaxSiTi = 1U << (kMaxCellLevel + 1); // Convert an s- or t-value to the corresponding u- or v-value. This is // a non-linear transformation from [-1,1] to [-1,1] that attempts to // make the cell sizes more uniform. double STtoUV(double s); // The inverse of the STtoUV transformation. Note that it is not always // true that UVtoST(STtoUV(x)) == x due to numerical errors. double UVtoST(double u); // Convert the i- or j-index of a leaf cell to the minimum corresponding s- // or t-value contained by that cell. The argument must be in the range // [0..2**30], i.e. up to one position beyond the normal range of valid leaf // cell indices. double IJtoSTMin(int i); // Return the i- or j-index of the leaf cell containing the given // s- or t-value. If the argument is outside the range spanned by valid // leaf cell indices, return the index of the closest valid leaf cell (i.e., // return values are clamped to the range of valid leaf cell indices). int STtoIJ(double s); // Convert an si- or ti-value to the corresponding s- or t-value. double SiTitoST(unsigned int si); // Return the si- or ti-coordinate that is nearest to the given s- or // t-value. The result may be outside the range of valid (si,ti)-values. unsigned int STtoSiTi(double s); // Convert (face, u, v) coordinates to a direction vector (not // necessarily unit length). S2Point FaceUVtoXYZ(int face, double u, double v); S2Point FaceUVtoXYZ(int face, R2Point const& uv); // If the dot product of p with the given face normal is positive, // set the corresponding u and v values (which may lie outside the range // [-1,1]) and return true. Otherwise return false. bool FaceXYZtoUV(int face, S2Point const& p, double* pu, double* pv); bool FaceXYZtoUV(int face, S2Point const& p, R2Point* puv); // Given a *valid* face for the given point p (meaning that dot product // of p with the face normal is positive), return the corresponding // u and v values (which may lie outside the range [-1,1]). void ValidFaceXYZtoUV(int face, S2Point const& p, double* pu, double* pv); void ValidFaceXYZtoUV(int face, S2Point const& p, R2Point* puv); // Transform the given point P to the (u,v,w) coordinate frame of the given // face (where the w-axis represents the face normal). S2Point FaceXYZtoUVW(int face, S2Point const& p); // Return the face containing the given direction vector. (For points on // the boundary between faces, the result is arbitrary but repeatable.) int GetFace(S2Point const& p); // Convert a direction vector (not necessarily unit length) to // (face, u, v) coordinates. int XYZtoFaceUV(S2Point const& p, double* pu, double* pv); int XYZtoFaceUV(S2Point const& p, R2Point* puv); // Convert a direction vector (not necessarily unit length) to // (face, si, ti) coordinates and, if p is exactly equal to the center of a // cell, return the level of this cell (-1 otherwise). int XYZtoFaceSiTi(S2Point const& p, int* face, unsigned int* si, unsigned int* ti); // Convert (face, si, ti) coordinates to a direction vector (not necessarily // unit length). S2Point FaceSiTitoXYZ(int face, unsigned int si, unsigned int ti); // Return the right-handed normal (not necessarily unit length) for an // edge in the direction of the positive v-axis at the given u-value on // the given face. (This vector is perpendicular to the plane through // the sphere origin that contains the given edge.) S2Point GetUNorm(int face, double u); // Return the right-handed normal (not necessarily unit length) for an // edge in the direction of the positive u-axis at the given v-value on // the given face. S2Point GetVNorm(int face, double v); // Return the unit-length normal, u-axis, or v-axis for the given face. S2Point GetNorm(int face); S2Point GetUAxis(int face); S2Point GetVAxis(int face); // Return the given axis of the given face (u=0, v=1, w=2). S2Point GetUVWAxis(int face, int axis); // With respect to the (u,v,w) coordinate system of a given face, return the // face that lies in the given direction (negative=0, positive=1) of the // given axis (u=0, v=1, w=2). For example, GetUVWFace(4, 0, 1) returns the // face that is adjacent to face 4 in the positive u-axis direction. int GetUVWFace(int face, int axis, int direction); ////////////////// Implementation details follow //////////////////// // We have implemented three different projections from cell-space (s,t) to // cube-space (u,v): linear, quadratic, and tangent. They have the following // tradeoffs: // // Linear - This is the fastest transformation, but also produces the least // uniform cell sizes. Cell areas vary by a factor of about 5.2, with the // largest cells at the center of each face and the smallest cells in // the corners. // // Tangent - Transforming the coordinates via atan() makes the cell sizes // more uniform. The areas vary by a maximum ratio of 1.4 as opposed to a // maximum ratio of 5.2. However, each call to atan() is about as expensive // as all of the other calculations combined when converting from points to // cell ids, i.e. it reduces performance by a factor of 3. // // Quadratic - This is an approximation of the tangent projection that // is much faster and produces cells that are almost as uniform in size. // It is about 3 times faster than the tangent projection for converting // cell ids to points or vice versa. Cell areas vary by a maximum ratio of // about 2.1. // // Here is a table comparing the cell uniformity using each projection. "Area // ratio" is the maximum ratio over all subdivision levels of the largest cell // area to the smallest cell area at that level, "edge ratio" is the maximum // ratio of the longest edge of any cell to the shortest edge of any cell at // the same level, and "diag ratio" is the ratio of the longest diagonal of // any cell to the shortest diagonal of any cell at the same level. "ToPoint" // and "FromPoint" are the times in microseconds required to convert cell ids // to and from points (unit vectors) respectively. "ToPointRaw" is the time // to convert to a non-unit-length vector, which is all that is needed for // some purposes. // // Area Edge Diag ToPointRaw ToPoint FromPoint // Ratio Ratio Ratio (microseconds) // ------------------------------------------------------------------- // Linear: 5.200 2.117 2.959 0.020 0.087 0.085 // Tangent: 1.414 1.414 1.704 0.237 0.299 0.258 // Quadratic: 2.082 1.802 1.932 0.033 0.096 0.108 // // The worst-case cell aspect ratios are about the same with all three // projections. The maximum ratio of the longest edge to the shortest edge // within the same cell is about 1.4 and the maximum ratio of the diagonals // within the same cell is about 1.7. // // This data was produced using s2cell_test and s2cellid_test. #define S2_LINEAR_PROJECTION 0 #define S2_TAN_PROJECTION 1 #define S2_QUADRATIC_PROJECTION 2 #define S2_PROJECTION S2_QUADRATIC_PROJECTION #if S2_PROJECTION == S2_LINEAR_PROJECTION inline double STtoUV(double s) { return 2 * s - 1; } inline double UVtoST(double u) { return 0.5 * (u + 1); } #elif S2_PROJECTION == S2_TAN_PROJECTION inline double STtoUV(double s) { // Unfortunately, tan(M_PI_4) is slightly less than 1.0. This isn't due to // a flaw in the implementation of tan(), it's because the derivative of // tan(x) at x=pi/4 is 2, and it happens that the two adjacent floating // point numbers on either side of the infinite-precision value of pi/4 have // tangents that are slightly below and slightly above 1.0 when rounded to // the nearest double-precision result. s = std::tan(M_PI_2 * s - M_PI_4); return s + (1.0 / (GG_LONGLONG(1) << 53)) * s; } inline double UVtoST(double u) { volatile double a = std::atan(u); return (2 * M_1_PI) * (a + M_PI_4); } #elif S2_PROJECTION == S2_QUADRATIC_PROJECTION inline double STtoUV(double s) { if (s >= 0.5) return (1 / 3.) * (4 * s * s - 1); else return (1 / 3.) * (1 - 4 * (1 - s) * (1 - s)); } inline double UVtoST(double u) { if (u >= 0) return 0.5 * std::sqrt(1 + 3 * u); else return 1 - 0.5 * std::sqrt(1 - 3 * u); } #else #error Unknown value for S2_PROJECTION #endif inline double IJtoSTMin(int i) { DCHECK(i >= 0 && i <= kLimitIJ); return (1.0 / kLimitIJ) * i; } inline int STtoIJ(double s) { return std::max( 0, std::min(kLimitIJ - 1, MathUtil::FastIntRound(kLimitIJ * s - 0.5))); } inline double SiTitoST(unsigned int si) { DCHECK(si >= 0 && si <= kMaxSiTi); return (1.0 / kMaxSiTi) * si; } inline unsigned int STtoSiTi(double s) { // kMaxSiTi == 2^31, so the result doesn't fit in an int32_t when s == 1. return static_cast(MathUtil::FastInt64Round(s * kMaxSiTi)); } inline S2Point FaceUVtoXYZ(int face, double u, double v) { switch (face) { case 0: return S2Point(1, u, v); case 1: return S2Point(-u, 1, v); case 2: return S2Point(-u, -v, 1); case 3: return S2Point(-1, -v, -u); case 4: return S2Point(v, -1, -u); default: return S2Point(v, u, -1); } } inline S2Point FaceUVtoXYZ(int face, R2Point const& uv) { return FaceUVtoXYZ(face, uv, uv); } inline void ValidFaceXYZtoUV(int face, S2Point const& p, double* pu, double* pv) { DCHECK_GT(p.DotProd(GetNorm(face)), 0); switch (face) { case 0: *pu = p / p; *pv = p / p; break; case 1: *pu = -p / p; *pv = p / p; break; case 2: *pu = -p / p; *pv = -p / p; break; case 3: *pu = p / p; *pv = p / p; break; case 4: *pu = p / p; *pv = -p / p; break; default: *pu = -p / p; *pv = -p / p; break; } } inline void ValidFaceXYZtoUV(int face, S2Point const& p, R2Point* puv) { ValidFaceXYZtoUV(face, p, &(*puv), &(*puv)); } inline int GetFace(S2Point const& p) { int face = p.LargestAbsComponent(); if (p[face] < 0) face += 3; return face; } inline int XYZtoFaceUV(S2Point const& p, double* pu, double* pv) { int face = GetFace(p); ValidFaceXYZtoUV(face, p, pu, pv); return face; } inline int XYZtoFaceUV(S2Point const& p, R2Point* puv) { return XYZtoFaceUV(p, &(*puv), &(*puv)); } inline bool FaceXYZtoUV(int face, S2Point const& p, double* pu, double* pv) { if (face < 3) { if (p[face] <= 0) return false; } else { if (p[face - 3] >= 0) return false; } ValidFaceXYZtoUV(face, p, pu, pv); return true; } inline bool FaceXYZtoUV(int face, S2Point const& p, R2Point* puv) { return FaceXYZtoUV(face, p, &(*puv), &(*puv)); } inline S2Point GetUNorm(int face, double u) { switch (face) { case 0: return S2Point(u, -1, 0); case 1: return S2Point(1, u, 0); case 2: return S2Point(1, 0, u); case 3: return S2Point(-u, 0, 1); case 4: return S2Point(0, -u, 1); default: return S2Point(0, -1, -u); } } inline S2Point GetVNorm(int face, double v) { switch (face) { case 0: return S2Point(-v, 0, 1); case 1: return S2Point(0, -v, 1); case 2: return S2Point(0, -1, -v); case 3: return S2Point(v, -1, 0); case 4: return S2Point(1, v, 0); default: return S2Point(1, 0, v); } } inline S2Point GetNorm(int face) { return GetUVWAxis(face, 2); } inline S2Point GetUAxis(int face) { return GetUVWAxis(face, 0); } inline S2Point GetVAxis(int face) { return GetUVWAxis(face, 1); } inline S2Point GetUVWAxis(int face, int axis) { double const* p = internal::kFaceUVWAxes[face][axis]; return S2Point(p, p, p); } inline int GetUVWFace(int face, int axis, int direction) { DCHECK(face >= 0 && face <= 5); DCHECK(axis >= 0 && axis <= 2); DCHECK(direction >= 0 && direction <= 1); return internal::kFaceUVWFaces[face][axis][direction]; } } // namespace S2 #endif // S2_S2COORDS_H_