1// Copyright 2005 Google Inc. All Rights Reserved.
2
3#include "base/integral_types.h"
4#include "base/logging.h"
5#include "s2.h"
6#include "s2cap.h"
7#include "s2cell.h"
8#include "s2latlngrect.h"
9
10namespace {
11
12// Multiply a positive number by this constant to ensure that the result
13// of a floating point operation is at least as large as the true
14// infinite-precision result.
15double const kRoundUp = 1.0 + 1.0 / (uint64(1) << 52);
16
17// Return the cap height corresponding to the given non-negative cap angle in
18// radians. Cap angles of Pi radians or larger yield a full cap.
19double GetHeightForAngle(double radians) {
20 DCHECK_GE(radians, 0);
21
22 // Caps of Pi radians or more are full.
23 if (radians >= M_PI) return 2;
24
25 // The height of the cap can be computed as 1 - cos(radians), but this isn't
26 // very accurate for angles close to zero (where cos(radians) is almost 1).
27 // Computing it as 2 * (sin(radians / 2) ** 2) gives much better precision.
28 double d = sin(0.5 * radians);
29 return 2 * d * d;
30}
31
32} // namespace
33
34S2Cap S2Cap::FromAxisAngle(S2Point const& axis, S1Angle const& angle) {
35 DCHECK(S2::IsUnitLength(axis));
36 DCHECK_GE(angle.radians(), 0);
37 return S2Cap(axis, GetHeightForAngle(angle.radians()));
38}
39
40S1Angle S2Cap::angle() const {
41 // This could also be computed as acos(1 - height_), but the following
42 // formula is much more accurate when the cap height is small. It
43 // follows from the relationship h = 1 - cos(theta) = 2 sin^2(theta/2).
44 if (is_empty()) return S1Angle::Radians(-1);
45 return S1Angle::Radians(2 * asin(sqrt(0.5 * height_)));
46}
47
48S2Cap S2Cap::Complement() const {
49 // The complement of a full cap is an empty cap, not a singleton.
50 // Also make sure that the complement of an empty cap has height 2.
51 double height = is_full() ? -1 : 2 - max(height_, 0.0);
52 return S2Cap::FromAxisHeight(-axis_, height);
53}
54
55bool S2Cap::Contains(S2Cap const& other) const {
56 if (is_full() || other.is_empty()) return true;
57 return angle().radians() >= axis_.Angle(other.axis_) +
58 other.angle().radians();
59}
60
61bool S2Cap::Intersects(S2Cap const& other) const {
62 if (is_empty() || other.is_empty()) return false;
63
64 return (angle().radians() + other.angle().radians() >=
65 axis_.Angle(other.axis_));
66}
67
68bool S2Cap::InteriorIntersects(S2Cap const& other) const {
69 // Make sure this cap has an interior and the other cap is non-empty.
70 if (height_ <= 0 || other.is_empty()) return false;
71
72 return (angle().radians() + other.angle().radians() >
73 axis_.Angle(other.axis_));
74}
75
76void S2Cap::AddPoint(S2Point const& p) {
77 // Compute the squared chord length, then convert it into a height.
78 DCHECK(S2::IsUnitLength(p));
79 if (is_empty()) {
80 axis_ = p;
81 height_ = 0;
82 } else {
83 // To make sure that the resulting cap actually includes this point,
84 // we need to round up the distance calculation. That is, after
85 // calling cap.AddPoint(p), cap.Contains(p) should be true.
86 double dist2 = (axis_ - p).Norm2();
87 height_ = max(height_, kRoundUp * 0.5 * dist2);
88 }
89}
90
91void S2Cap::AddCap(S2Cap const& other) {
92 if (is_empty()) {
93 *this = other;
94 } else {
95 // See comments for AddPoint(). This could be optimized by doing the
96 // calculation in terms of cap heights rather than cap opening angles.
97 double radians = axis_.Angle(other.axis_) + other.angle().radians();
98 height_ = max(height_, kRoundUp * GetHeightForAngle(radians));
99 }
100}
101
102S2Cap S2Cap::Expanded(S1Angle const& distance) const {
103 DCHECK_GE(distance.radians(), 0);
104 if (is_empty()) return Empty();
105 return FromAxisAngle(axis_, angle() + distance);
106}
107
108S2Cap* S2Cap::Clone() const {
109 return new S2Cap(*this);
110}
111
112S2Cap S2Cap::GetCapBound() const {
113 return *this;
114}
115
116S2LatLngRect S2Cap::GetRectBound() const {
117 if (is_empty()) return S2LatLngRect::Empty();
118
119 // Convert the axis to a (lat,lng) pair, and compute the cap angle.
120 S2LatLng axis_ll(axis_);
121 double cap_angle = angle().radians();
122
123 bool all_longitudes = false;
124 double lat[2], lng[2];
125 lng[0] = -M_PI;
126 lng[1] = M_PI;
127
128 // Check whether cap includes the south pole.
129 lat[0] = axis_ll.lat().radians() - cap_angle;
130 if (lat[0] <= -M_PI_2) {
131 lat[0] = -M_PI_2;
132 all_longitudes = true;
133 }
134 // Check whether cap includes the north pole.
135 lat[1] = axis_ll.lat().radians() + cap_angle;
136 if (lat[1] >= M_PI_2) {
137 lat[1] = M_PI_2;
138 all_longitudes = true;
139 }
140 if (!all_longitudes) {
141 // Compute the range of longitudes covered by the cap. We use the law
142 // of sines for spherical triangles. Consider the triangle ABC where
143 // A is the north pole, B is the center of the cap, and C is the point
144 // of tangency between the cap boundary and a line of longitude. Then
145 // C is a right angle, and letting a,b,c denote the sides opposite A,B,C,
146 // we have sin(a)/sin(A) = sin(c)/sin(C), or sin(A) = sin(a)/sin(c).
147 // Here "a" is the cap angle, and "c" is the colatitude (90 degrees
148 // minus the latitude). This formula also works for negative latitudes.
149 //
150 // The formula for sin(a) follows from the relationship h = 1 - cos(a).
151
152 double sin_a = sqrt(height_ * (2 - height_));
153 double sin_c = cos(axis_ll.lat().radians());
154 if (sin_a <= sin_c) {
155 double angle_A = asin(sin_a / sin_c);
156 lng[0] = drem(axis_ll.lng().radians() - angle_A, 2 * M_PI);
157 lng[1] = drem(axis_ll.lng().radians() + angle_A, 2 * M_PI);
158 }
159 }
160 return S2LatLngRect(R1Interval(lat[0], lat[1]),
161 S1Interval(lng[0], lng[1]));
162}
163
164bool S2Cap::Intersects(S2Cell const& cell, S2Point const* vertices) const {
165 // Return true if this cap intersects any point of 'cell' excluding its
166 // vertices (which are assumed to already have been checked).
167
168 // If the cap is a hemisphere or larger, the cell and the complement of the
169 // cap are both convex. Therefore since no vertex of the cell is contained,
170 // no other interior point of the cell is contained either.
171 if (height_ >= 1) return false;
172
173 // We need to check for empty caps due to the axis check just below.
174 if (is_empty()) return false;
175
176 // Optimization: return true if the cell contains the cap axis. (This
177 // allows half of the edge checks below to be skipped.)
178 if (cell.Contains(axis_)) return true;
179
180 // At this point we know that the cell does not contain the cap axis,
181 // and the cap does not contain any cell vertex. The only way that they
182 // can intersect is if the cap intersects the interior of some edge.
183
184 double sin2_angle = height_ * (2 - height_); // sin^2(cap_angle)
185 for (int k = 0; k < 4; ++k) {
186 S2Point edge = cell.GetEdgeRaw(k);
187 double dot = axis_.DotProd(edge);
188 if (dot > 0) {
189 // The axis is in the interior half-space defined by the edge. We don't
190 // need to consider these edges, since if the cap intersects this edge
191 // then it also intersects the edge on the opposite side of the cell
192 // (because we know the axis is not contained with the cell).
193 continue;
194 }
195 // The Norm2() factor is necessary because "edge" is not normalized.
196 if (dot * dot > sin2_angle * edge.Norm2()) {
197 return false; // Entire cap is on the exterior side of this edge.
198 }
199 // Otherwise, the great circle containing this edge intersects
200 // the interior of the cap. We just need to check whether the point
201 // of closest approach occurs between the two edge endpoints.
202 S2Point dir = edge.CrossProd(axis_);
203 if (dir.DotProd(vertices[k]) < 0 && dir.DotProd(vertices[(k+1)&3]) > 0)
204 return true;
205 }
206 return false;
207}
208
209bool S2Cap::Contains(S2Cell const& cell) const {
210 // If the cap does not contain all cell vertices, return false.
211 // We check the vertices before taking the Complement() because we can't
212 // accurately represent the complement of a very small cap (a height
213 // of 2-epsilon is rounded off to 2).
214 S2Point vertices[4];
215 for (int k = 0; k < 4; ++k) {
216 vertices[k] = cell.GetVertex(k);
217 if (!Contains(vertices[k])) return false;
218 }
219 // Otherwise, return true if the complement of the cap does not intersect
220 // the cell. (This test is slightly conservative, because technically we
221 // want Complement().InteriorIntersects() here.)
222 return !Complement().Intersects(cell, vertices);
223}
224
225bool S2Cap::MayIntersect(S2Cell const& cell) const {
226 // If the cap contains any cell vertex, return true.
227 S2Point vertices[4];
228 for (int k = 0; k < 4; ++k) {
229 vertices[k] = cell.GetVertex(k);
230 if (Contains(vertices[k])) return true;
231 }
232 return Intersects(cell, vertices);
233}
234
235bool S2Cap::Contains(S2Point const& p) const {
236 DCHECK(S2::IsUnitLength(p));
237 return (axis_ - p).Norm2() <= 2 * height_;
238}
239
240bool S2Cap::InteriorContains(S2Point const& p) const {
241 DCHECK(S2::IsUnitLength(p));
242 return is_full() || (axis_ - p).Norm2() < 2 * height_;
243}
244
245bool S2Cap::operator==(S2Cap const& other) const {
246 return (axis_ == other.axis_ && height_ == other.height_) ||
247 (is_empty() && other.is_empty()) ||
248 (is_full() && other.is_full());
249}
250
251bool S2Cap::ApproxEquals(S2Cap const& other, double max_error) {
252 return (S2::ApproxEquals(axis_, other.axis_, max_error) &&
253 fabs(height_ - other.height_) <= max_error) ||
254 (is_empty() && other.height_ <= max_error) ||
255 (other.is_empty() && height_ <= max_error) ||
256 (is_full() && other.height_ >= 2 - max_error) ||
257 (other.is_full() && height_ >= 2 - max_error);
258}
259
260ostream& operator<<(ostream& os, S2Cap const& cap) {
261 return os << "[Axis=" << cap.axis() << ", Angle=" << cap.angle() << "]";
262}
263