chromium/third_party/s2cellid/src/s2/s2coords.h

// 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: [email protected] (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 <algorithm>
#include <cmath>

#include "base/check_op.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 = …;

// The maximum index of a valid leaf cell plus one.  The range of valid leaf
// cell indices is [0..kLimitIJ-1].
int const kLimitIJ = …;  // == 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 = …;

// 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 …
#define S2_TAN_PROJECTION …
#define S2_QUADRATIC_PROJECTION …

#define S2_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) { … }

inline double UVtoST(double u) { … }

#else

#error Unknown value for S2_PROJECTION

#endif

inline double IJtoSTMin(int i) { … }

inline int STtoIJ(double s) { … }

inline double SiTitoST(unsigned int si) { … }

inline unsigned int STtoSiTi(double s) { … }

inline S2Point FaceUVtoXYZ(int face, double u, double v) { … }

inline S2Point FaceUVtoXYZ(int face, R2Point const& uv) { … }

inline void ValidFaceXYZtoUV(int face,
                             S2Point const& p,
                             double* pu,
                             double* pv) { … }

inline void ValidFaceXYZtoUV(int face, S2Point const& p, R2Point* puv) { … }

inline int GetFace(S2Point const& p) { … }

inline int XYZtoFaceUV(S2Point const& p, double* pu, double* pv) { … }

inline int XYZtoFaceUV(S2Point const& p, R2Point* puv) { … }

inline bool FaceXYZtoUV(int face, S2Point const& p, double* pu, double* pv) { … }

inline bool FaceXYZtoUV(int face, S2Point const& p, R2Point* puv) { … }

inline S2Point GetUNorm(int face, double u) { … }

inline S2Point GetVNorm(int face, double v) { … }

inline S2Point GetNorm(int face) { … }

inline S2Point GetUAxis(int face) { … }

inline S2Point GetVAxis(int face) { … }

inline S2Point GetUVWAxis(int face, int axis) { … }

inline int GetUVWFace(int face, int axis, int direction) { … }

}  // namespace S2

#endif  // S2_S2COORDS_H_