1// Copyright 2005 Google Inc. All Rights Reserved.
2
3#include "s2cell.h"
4
5#include "base/integral_types.h"
6#include "base/logging.h"
7#include "s2.h"
8#include "s2cap.h"
9#include "s2latlngrect.h"
10#include "util/math/vector2-inl.h"
11
12// Since S2Cells are copied by value, the following assertion is a reminder
13// not to add fields unnecessarily. An S2Cell currently consists of 43 data
14// bytes, one vtable pointer, plus alignment overhead. This works out to 48
15// bytes on 32 bit architectures and 56 bytes on 64 bit architectures.
16//
17// The expression below rounds up (43 + sizeof(void*)) to the nearest
18// multiple of sizeof(void*).
19COMPILE_ASSERT(sizeof(S2Cell) <= ((43+2*sizeof(void*)-1) & -sizeof(void*)),
20 S2Cell_is_getting_bloated);
21
22S2Point S2Cell::GetVertexRaw(int k) const {
23 // Vertices are returned in the order SW, SE, NE, NW.
24 return S2::FaceUVtoXYZ(face_, uv_[0][(k>>1) ^ (k&1)], uv_[1][k>>1]);
25}
26
27S2Point S2Cell::GetEdgeRaw(int k) const {
28 switch (k) {
29 case 0: return S2::GetVNorm(face_, uv_[1][0]); // South
30 case 1: return S2::GetUNorm(face_, uv_[0][1]); // East
31 case 2: return -S2::GetVNorm(face_, uv_[1][1]); // North
32 default: return -S2::GetUNorm(face_, uv_[0][0]); // West
33 }
34}
35
36void S2Cell::Init(S2CellId const& id) {
37 id_ = id;
38 int ij[2], orientation;
39 face_ = id.ToFaceIJOrientation(&ij[0], &ij[1], &orientation);
40 orientation_ = orientation; // Compress int to a byte.
41 level_ = id.level();
42 int cellSize = GetSizeIJ(); // Depends on level_.
43 for (int d = 0; d < 2; ++d) {
44 int ij_lo = ij[d] & -cellSize;
45 int ij_hi = ij_lo + cellSize;
46 uv_[d][0] = S2::STtoUV((1.0 / S2CellId::kMaxSize) * ij_lo);
47 uv_[d][1] = S2::STtoUV((1.0 / S2CellId::kMaxSize) * ij_hi);
48 }
49}
50
51bool S2Cell::Subdivide(S2Cell children[4]) const {
52 // This function is equivalent to just iterating over the child cell ids
53 // and calling the S2Cell constructor, but it is about 2.5 times faster.
54
55 if (id_.is_leaf()) return false;
56
57 // Compute the cell midpoint in uv-space.
58 Vector2_d uv_mid = id_.GetCenterUV();
59
60 // Create four children with the appropriate bounds.
61 S2CellId id = id_.child_begin();
62 for (int pos = 0; pos < 4; ++pos, id = id.next()) {
63 S2Cell *child = &children[pos];
64 child->face_ = face_;
65 child->level_ = level_ + 1;
66 child->orientation_ = orientation_ ^ S2::kPosToOrientation[pos];
67 child->id_ = id;
68 // We want to split the cell in half in "u" and "v". To decide which
69 // side to set equal to the midpoint value, we look at cell's (i,j)
70 // position within its parent. The index for "i" is in bit 1 of ij.
71 int ij = S2::kPosToIJ[orientation_][pos];
72 int i = ij >> 1;
73 int j = ij & 1;
74 child->uv_[0][i] = uv_[0][i];
75 child->uv_[0][1-i] = uv_mid[0];
76 child->uv_[1][j] = uv_[1][j];
77 child->uv_[1][1-j] = uv_mid[1];
78 }
79 return true;
80}
81
82S2Point S2Cell::GetCenterRaw() const {
83 return id_.ToPointRaw();
84}
85
86double S2Cell::AverageArea(int level) {
87 return S2::kAvgArea.GetValue(level);
88}
89
90double S2Cell::ApproxArea() const {
91 // All cells at the first two levels have the same area.
92 if (level_ < 2) return AverageArea(level_);
93
94 // First, compute the approximate area of the cell when projected
95 // perpendicular to its normal. The cross product of its diagonals gives
96 // the normal, and the length of the normal is twice the projected area.
97 double flat_area = 0.5 * (GetVertex(2) - GetVertex(0)).
98 CrossProd(GetVertex(3) - GetVertex(1)).Norm();
99
100 // Now, compensate for the curvature of the cell surface by pretending
101 // that the cell is shaped like a spherical cap. The ratio of the
102 // area of a spherical cap to the area of its projected disc turns out
103 // to be 2 / (1 + sqrt(1 - r*r)) where "r" is the radius of the disc.
104 // For example, when r=0 the ratio is 1, and when r=1 the ratio is 2.
105 // Here we set Pi*r*r == flat_area to find the equivalent disc.
106 return flat_area * 2 / (1 + sqrt(1 - min(M_1_PI * flat_area, 1.0)));
107}
108
109double S2Cell::ExactArea() const {
110 S2Point v0 = GetVertex(0);
111 S2Point v1 = GetVertex(1);
112 S2Point v2 = GetVertex(2);
113 S2Point v3 = GetVertex(3);
114 return S2::Area(v0, v1, v2) + S2::Area(v0, v2, v3);
115}
116
117S2Cell* S2Cell::Clone() const {
118 return new S2Cell(*this);
119}
120
121S2Cap S2Cell::GetCapBound() const {
122 // Use the cell center in (u,v)-space as the cap axis. This vector is
123 // very close to GetCenter() and faster to compute. Neither one of these
124 // vectors yields the bounding cap with minimal surface area, but they
125 // are both pretty close.
126 //
127 // It's possible to show that the two vertices that are furthest from
128 // the (u,v)-origin never determine the maximum cap size (this is a
129 // possible future optimization).
130
131 double u = 0.5 * (uv_[0][0] + uv_[0][1]);
132 double v = 0.5 * (uv_[1][0] + uv_[1][1]);
133 S2Cap cap = S2Cap::FromAxisHeight(S2::FaceUVtoXYZ(face_,u,v).Normalize(), 0);
134 for (int k = 0; k < 4; ++k) {
135 cap.AddPoint(GetVertex(k));
136 }
137 return cap;
138}
139
140inline double S2Cell::GetLatitude(int i, int j) const {
141 S2Point p = S2::FaceUVtoXYZ(face_, uv_[0][i], uv_[1][j]);
142 return S2LatLng::Latitude(p).radians();
143}
144
145inline double S2Cell::GetLongitude(int i, int j) const {
146 S2Point p = S2::FaceUVtoXYZ(face_, uv_[0][i], uv_[1][j]);
147 return S2LatLng::Longitude(p).radians();
148}
149
150S2LatLngRect S2Cell::GetRectBound() const {
151 if (level_ > 0) {
152 // Except for cells at level 0, the latitude and longitude extremes are
153 // attained at the vertices. Furthermore, the latitude range is
154 // determined by one pair of diagonally opposite vertices and the
155 // longitude range is determined by the other pair.
156 //
157 // We first determine which corner (i,j) of the cell has the largest
158 // absolute latitude. To maximize latitude, we want to find the point in
159 // the cell that has the largest absolute z-coordinate and the smallest
160 // absolute x- and y-coordinates. To do this we look at each coordinate
161 // (u and v), and determine whether we want to minimize or maximize that
162 // coordinate based on the axis direction and the cell's (u,v) quadrant.
163 double u = uv_[0][0] + uv_[0][1];
164 double v = uv_[1][0] + uv_[1][1];
165 int i = S2::GetUAxis(face_)[2] == 0 ? (u < 0) : (u > 0);
166 int j = S2::GetVAxis(face_)[2] == 0 ? (v < 0) : (v > 0);
167
168 // We grow the bounds slightly to make sure that the bounding rectangle
169 // also contains the normalized versions of the vertices. Note that the
170 // maximum result magnitude is Pi, with a floating-point exponent of 1.
171 // Therefore adding or subtracting 2**-51 will always change the result.
172 static double const kMaxError = 1.0 / (int64(1) << 51);
173 R1Interval lat = R1Interval::FromPointPair(GetLatitude(i, j),
174 GetLatitude(1-i, 1-j));
175 lat = lat.Expanded(kMaxError).Intersection(S2LatLngRect::FullLat());
176 if (lat.lo() == -M_PI_2 || lat.hi() == M_PI_2) {
177 return S2LatLngRect(lat, S1Interval::Full());
178 }
179 S1Interval lng = S1Interval::FromPointPair(GetLongitude(i, 1-j),
180 GetLongitude(1-i, j));
181 return S2LatLngRect(lat, lng.Expanded(kMaxError));
182 }
183
184 // The 4 cells around the equator extend to +/-45 degrees latitude at the
185 // midpoints of their top and bottom edges. The two cells covering the
186 // poles extend down to +/-35.26 degrees at their vertices.
187 static double const kPoleMinLat = asin(sqrt(1./3)); // 35.26 degrees
188
189 // The face centers are the +X, +Y, +Z, -X, -Y, -Z axes in that order.
190 DCHECK_EQ(((face_ < 3) ? 1 : -1), S2::GetNorm(face_)[face_ % 3]);
191 switch (face_) {
192 case 0: return S2LatLngRect(R1Interval(-M_PI_4, M_PI_4),
193 S1Interval(-M_PI_4, M_PI_4));
194 case 1: return S2LatLngRect(R1Interval(-M_PI_4, M_PI_4),
195 S1Interval(M_PI_4, 3*M_PI_4));
196 case 2: return S2LatLngRect(R1Interval(kPoleMinLat, M_PI_2),
197 S1Interval(-M_PI, M_PI));
198 case 3: return S2LatLngRect(R1Interval(-M_PI_4, M_PI_4),
199 S1Interval(3*M_PI_4, -3*M_PI_4));
200 case 4: return S2LatLngRect(R1Interval(-M_PI_4, M_PI_4),
201 S1Interval(-3*M_PI_4, -M_PI_4));
202 default: return S2LatLngRect(R1Interval(-M_PI_2, -kPoleMinLat),
203 S1Interval(-M_PI, M_PI));
204 }
205}
206
207bool S2Cell::MayIntersect(S2Cell const& cell) const {
208 return id_.intersects(cell.id_);
209}
210
211bool S2Cell::Contains(S2Cell const& cell) const {
212 return id_.contains(cell.id_);
213}
214
215bool S2Cell::Contains(S2Point const& p) const {
216 // We can't just call XYZtoFaceUV, because for points that lie on the
217 // boundary between two faces (i.e. u or v is +1/-1) we need to return
218 // true for both adjacent cells.
219 double u, v;
220 if (!S2::FaceXYZtoUV(face_, p, &u, &v)) return false;
221 return (u >= uv_[0][0] && u <= uv_[0][1] &&
222 v >= uv_[1][0] && v <= uv_[1][1]);
223}
224