1 // Copyright 2014 Google Inc. All rights reserved. 2 // 3 // Licensed under the Apache License, Version 2.0 (the "License"); 4 // you may not use this file except in compliance with the License. 5 // You may obtain a copy of the License at 6 // 7 // http://www.apache.org/licenses/LICENSE-2.0 8 // 9 // Unless required by applicable law or agreed to in writing, software 10 // distributed under the License is distributed on an "AS IS" BASIS, 11 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 12 // See the License for the specific language governing permissions and 13 // limitations under the License. 14 15 package s2 16 17 import ( 18 "fmt" 19 "io" 20 "math" 21 22 "github.com/golang/geo/r1" 23 "github.com/golang/geo/s1" 24 ) 25 26 var ( 27 // centerPoint is the default center for Caps 28 centerPoint = PointFromCoords(1.0, 0, 0) 29 ) 30 31 // Cap represents a disc-shaped region defined by a center and radius. 32 // Technically this shape is called a "spherical cap" (rather than disc) 33 // because it is not planar; the cap represents a portion of the sphere that 34 // has been cut off by a plane. The boundary of the cap is the circle defined 35 // by the intersection of the sphere and the plane. For containment purposes, 36 // the cap is a closed set, i.e. it contains its boundary. 37 // 38 // For the most part, you can use a spherical cap wherever you would use a 39 // disc in planar geometry. The radius of the cap is measured along the 40 // surface of the sphere (rather than the straight-line distance through the 41 // interior). Thus a cap of radius π/2 is a hemisphere, and a cap of radius 42 // π covers the entire sphere. 43 // 44 // The center is a point on the surface of the unit sphere. (Hence the need for 45 // it to be of unit length.) 46 // 47 // A cap can also be defined by its center point and height. The height is the 48 // distance from the center point to the cutoff plane. There is also support for 49 // "empty" and "full" caps, which contain no points and all points respectively. 50 // 51 // Here are some useful relationships between the cap height (h), the cap 52 // radius (r), the maximum chord length from the cap's center (d), and the 53 // radius of cap's base (a). 54 // 55 // h = 1 - cos(r) 56 // = 2 * sin^2(r/2) 57 // d^2 = 2 * h 58 // = a^2 + h^2 59 // 60 // The zero value of Cap is an invalid cap. Use EmptyCap to get a valid empty cap. 61 type Cap struct { 62 center Point 63 radius s1.ChordAngle 64 } 65 66 // CapFromPoint constructs a cap containing a single point. 67 func CapFromPoint(p Point) Cap { 68 return CapFromCenterChordAngle(p, 0) 69 } 70 71 // CapFromCenterAngle constructs a cap with the given center and angle. 72 func CapFromCenterAngle(center Point, angle s1.Angle) Cap { 73 return CapFromCenterChordAngle(center, s1.ChordAngleFromAngle(angle)) 74 } 75 76 // CapFromCenterChordAngle constructs a cap where the angle is expressed as an 77 // s1.ChordAngle. This constructor is more efficient than using an s1.Angle. 78 func CapFromCenterChordAngle(center Point, radius s1.ChordAngle) Cap { 79 return Cap{ 80 center: center, 81 radius: radius, 82 } 83 } 84 85 // CapFromCenterHeight constructs a cap with the given center and height. A 86 // negative height yields an empty cap; a height of 2 or more yields a full cap. 87 // The center should be unit length. 88 func CapFromCenterHeight(center Point, height float64) Cap { 89 return CapFromCenterChordAngle(center, s1.ChordAngleFromSquaredLength(2*height)) 90 } 91 92 // CapFromCenterArea constructs a cap with the given center and surface area. 93 // Note that the area can also be interpreted as the solid angle subtended by the 94 // cap (because the sphere has unit radius). A negative area yields an empty cap; 95 // an area of 4*π or more yields a full cap. 96 func CapFromCenterArea(center Point, area float64) Cap { 97 return CapFromCenterChordAngle(center, s1.ChordAngleFromSquaredLength(area/math.Pi)) 98 } 99 100 // EmptyCap returns a cap that contains no points. 101 func EmptyCap() Cap { 102 return CapFromCenterChordAngle(centerPoint, s1.NegativeChordAngle) 103 } 104 105 // FullCap returns a cap that contains all points. 106 func FullCap() Cap { 107 return CapFromCenterChordAngle(centerPoint, s1.StraightChordAngle) 108 } 109 110 // IsValid reports whether the Cap is considered valid. 111 func (c Cap) IsValid() bool { 112 return c.center.Vector.IsUnit() && c.radius <= s1.StraightChordAngle 113 } 114 115 // IsEmpty reports whether the cap is empty, i.e. it contains no points. 116 func (c Cap) IsEmpty() bool { 117 return c.radius < 0 118 } 119 120 // IsFull reports whether the cap is full, i.e. it contains all points. 121 func (c Cap) IsFull() bool { 122 return c.radius == s1.StraightChordAngle 123 } 124 125 // Center returns the cap's center point. 126 func (c Cap) Center() Point { 127 return c.center 128 } 129 130 // Height returns the height of the cap. This is the distance from the center 131 // point to the cutoff plane. 132 func (c Cap) Height() float64 { 133 return float64(0.5 * c.radius) 134 } 135 136 // Radius returns the cap radius as an s1.Angle. (Note that the cap angle 137 // is stored internally as a ChordAngle, so this method requires a trigonometric 138 // operation and may yield a slightly different result than the value passed 139 // to CapFromCenterAngle). 140 func (c Cap) Radius() s1.Angle { 141 return c.radius.Angle() 142 } 143 144 // Area returns the surface area of the Cap on the unit sphere. 145 func (c Cap) Area() float64 { 146 return 2.0 * math.Pi * math.Max(0, c.Height()) 147 } 148 149 // Contains reports whether this cap contains the other. 150 func (c Cap) Contains(other Cap) bool { 151 // In a set containment sense, every cap contains the empty cap. 152 if c.IsFull() || other.IsEmpty() { 153 return true 154 } 155 return c.radius >= ChordAngleBetweenPoints(c.center, other.center).Add(other.radius) 156 } 157 158 // Intersects reports whether this cap intersects the other cap. 159 // i.e. whether they have any points in common. 160 func (c Cap) Intersects(other Cap) bool { 161 if c.IsEmpty() || other.IsEmpty() { 162 return false 163 } 164 165 return c.radius.Add(other.radius) >= ChordAngleBetweenPoints(c.center, other.center) 166 } 167 168 // InteriorIntersects reports whether this caps interior intersects the other cap. 169 func (c Cap) InteriorIntersects(other Cap) bool { 170 // Make sure this cap has an interior and the other cap is non-empty. 171 if c.radius <= 0 || other.IsEmpty() { 172 return false 173 } 174 175 return c.radius.Add(other.radius) > ChordAngleBetweenPoints(c.center, other.center) 176 } 177 178 // ContainsPoint reports whether this cap contains the point. 179 func (c Cap) ContainsPoint(p Point) bool { 180 return ChordAngleBetweenPoints(c.center, p) <= c.radius 181 } 182 183 // InteriorContainsPoint reports whether the point is within the interior of this cap. 184 func (c Cap) InteriorContainsPoint(p Point) bool { 185 return c.IsFull() || ChordAngleBetweenPoints(c.center, p) < c.radius 186 } 187 188 // Complement returns the complement of the interior of the cap. A cap and its 189 // complement have the same boundary but do not share any interior points. 190 // The complement operator is not a bijection because the complement of a 191 // singleton cap (containing a single point) is the same as the complement 192 // of an empty cap. 193 func (c Cap) Complement() Cap { 194 if c.IsFull() { 195 return EmptyCap() 196 } 197 if c.IsEmpty() { 198 return FullCap() 199 } 200 201 return CapFromCenterChordAngle(Point{c.center.Mul(-1)}, s1.StraightChordAngle.Sub(c.radius)) 202 } 203 204 // CapBound returns a bounding spherical cap. This is not guaranteed to be exact. 205 func (c Cap) CapBound() Cap { 206 return c 207 } 208 209 // RectBound returns a bounding latitude-longitude rectangle. 210 // The bounds are not guaranteed to be tight. 211 func (c Cap) RectBound() Rect { 212 if c.IsEmpty() { 213 return EmptyRect() 214 } 215 216 capAngle := c.Radius().Radians() 217 allLongitudes := false 218 lat := r1.Interval{ 219 Lo: latitude(c.center).Radians() - capAngle, 220 Hi: latitude(c.center).Radians() + capAngle, 221 } 222 lng := s1.FullInterval() 223 224 // Check whether cap includes the south pole. 225 if lat.Lo <= -math.Pi/2 { 226 lat.Lo = -math.Pi / 2 227 allLongitudes = true 228 } 229 230 // Check whether cap includes the north pole. 231 if lat.Hi >= math.Pi/2 { 232 lat.Hi = math.Pi / 2 233 allLongitudes = true 234 } 235 236 if !allLongitudes { 237 // Compute the range of longitudes covered by the cap. We use the law 238 // of sines for spherical triangles. Consider the triangle ABC where 239 // A is the north pole, B is the center of the cap, and C is the point 240 // of tangency between the cap boundary and a line of longitude. Then 241 // C is a right angle, and letting a,b,c denote the sides opposite A,B,C, 242 // we have sin(a)/sin(A) = sin(c)/sin(C), or sin(A) = sin(a)/sin(c). 243 // Here "a" is the cap angle, and "c" is the colatitude (90 degrees 244 // minus the latitude). This formula also works for negative latitudes. 245 // 246 // The formula for sin(a) follows from the relationship h = 1 - cos(a). 247 sinA := c.radius.Sin() 248 sinC := math.Cos(latitude(c.center).Radians()) 249 if sinA <= sinC { 250 angleA := math.Asin(sinA / sinC) 251 lng.Lo = math.Remainder(longitude(c.center).Radians()-angleA, math.Pi*2) 252 lng.Hi = math.Remainder(longitude(c.center).Radians()+angleA, math.Pi*2) 253 } 254 } 255 return Rect{lat, lng} 256 } 257 258 // Equal reports whether this cap is equal to the other cap. 259 func (c Cap) Equal(other Cap) bool { 260 return (c.radius == other.radius && c.center == other.center) || 261 (c.IsEmpty() && other.IsEmpty()) || 262 (c.IsFull() && other.IsFull()) 263 } 264 265 // ApproxEqual reports whether this cap is equal to the other cap within the given tolerance. 266 func (c Cap) ApproxEqual(other Cap) bool { 267 const epsilon = 1e-14 268 r2 := float64(c.radius) 269 otherR2 := float64(other.radius) 270 return c.center.ApproxEqual(other.center) && 271 math.Abs(r2-otherR2) <= epsilon || 272 c.IsEmpty() && otherR2 <= epsilon || 273 other.IsEmpty() && r2 <= epsilon || 274 c.IsFull() && otherR2 >= 2-epsilon || 275 other.IsFull() && r2 >= 2-epsilon 276 } 277 278 // AddPoint increases the cap if necessary to include the given point. If this cap is empty, 279 // then the center is set to the point with a zero height. p must be unit-length. 280 func (c Cap) AddPoint(p Point) Cap { 281 if c.IsEmpty() { 282 c.center = p 283 c.radius = 0 284 return c 285 } 286 287 // After calling cap.AddPoint(p), cap.Contains(p) must be true. However 288 // we don't need to do anything special to achieve this because Contains() 289 // does exactly the same distance calculation that we do here. 290 if newRad := ChordAngleBetweenPoints(c.center, p); newRad > c.radius { 291 c.radius = newRad 292 } 293 return c 294 } 295 296 // AddCap increases the cap height if necessary to include the other cap. If this cap is empty, 297 // it is set to the other cap. 298 func (c Cap) AddCap(other Cap) Cap { 299 if c.IsEmpty() { 300 return other 301 } 302 if other.IsEmpty() { 303 return c 304 } 305 306 // We round up the distance to ensure that the cap is actually contained. 307 // TODO(roberts): Do some error analysis in order to guarantee this. 308 dist := ChordAngleBetweenPoints(c.center, other.center).Add(other.radius) 309 if newRad := dist.Expanded(dblEpsilon * float64(dist)); newRad > c.radius { 310 c.radius = newRad 311 } 312 return c 313 } 314 315 // Expanded returns a new cap expanded by the given angle. If the cap is empty, 316 // it returns an empty cap. 317 func (c Cap) Expanded(distance s1.Angle) Cap { 318 if c.IsEmpty() { 319 return EmptyCap() 320 } 321 return CapFromCenterChordAngle(c.center, c.radius.Add(s1.ChordAngleFromAngle(distance))) 322 } 323 324 func (c Cap) String() string { 325 return fmt.Sprintf("[Center=%v, Radius=%f]", c.center.Vector, c.Radius().Degrees()) 326 } 327 328 // radiusToHeight converts an s1.Angle into the height of the cap. 329 func radiusToHeight(r s1.Angle) float64 { 330 if r.Radians() < 0 { 331 return float64(s1.NegativeChordAngle) 332 } 333 if r.Radians() >= math.Pi { 334 return float64(s1.RightChordAngle) 335 } 336 return float64(0.5 * s1.ChordAngleFromAngle(r)) 337 338 } 339 340 // ContainsCell reports whether the cap contains the given cell. 341 func (c Cap) ContainsCell(cell Cell) bool { 342 // If the cap does not contain all cell vertices, return false. 343 var vertices [4]Point 344 for k := 0; k < 4; k++ { 345 vertices[k] = cell.Vertex(k) 346 if !c.ContainsPoint(vertices[k]) { 347 return false 348 } 349 } 350 // Otherwise, return true if the complement of the cap does not intersect the cell. 351 return !c.Complement().intersects(cell, vertices) 352 } 353 354 // IntersectsCell reports whether the cap intersects the cell. 355 func (c Cap) IntersectsCell(cell Cell) bool { 356 // If the cap contains any cell vertex, return true. 357 var vertices [4]Point 358 for k := 0; k < 4; k++ { 359 vertices[k] = cell.Vertex(k) 360 if c.ContainsPoint(vertices[k]) { 361 return true 362 } 363 } 364 return c.intersects(cell, vertices) 365 } 366 367 // intersects reports whether the cap intersects any point of the cell excluding 368 // its vertices (which are assumed to already have been checked). 369 func (c Cap) intersects(cell Cell, vertices [4]Point) bool { 370 // If the cap is a hemisphere or larger, the cell and the complement of the cap 371 // are both convex. Therefore since no vertex of the cell is contained, no other 372 // interior point of the cell is contained either. 373 if c.radius >= s1.RightChordAngle { 374 return false 375 } 376 377 // We need to check for empty caps due to the center check just below. 378 if c.IsEmpty() { 379 return false 380 } 381 382 // Optimization: return true if the cell contains the cap center. This allows half 383 // of the edge checks below to be skipped. 384 if cell.ContainsPoint(c.center) { 385 return true 386 } 387 388 // At this point we know that the cell does not contain the cap center, and the cap 389 // does not contain any cell vertex. The only way that they can intersect is if the 390 // cap intersects the interior of some edge. 391 sin2Angle := c.radius.Sin2() 392 for k := 0; k < 4; k++ { 393 edge := cell.Edge(k).Vector 394 dot := c.center.Vector.Dot(edge) 395 if dot > 0 { 396 // The center is in the interior half-space defined by the edge. We do not need 397 // to consider these edges, since if the cap intersects this edge then it also 398 // intersects the edge on the opposite side of the cell, because the center is 399 // not contained with the cell. 400 continue 401 } 402 403 // The Norm2() factor is necessary because "edge" is not normalized. 404 if dot*dot > sin2Angle*edge.Norm2() { 405 return false 406 } 407 408 // Otherwise, the great circle containing this edge intersects the interior of the cap. We just 409 // need to check whether the point of closest approach occurs between the two edge endpoints. 410 dir := edge.Cross(c.center.Vector) 411 if dir.Dot(vertices[k].Vector) < 0 && dir.Dot(vertices[(k+1)&3].Vector) > 0 { 412 return true 413 } 414 } 415 return false 416 } 417 418 // CellUnionBound computes a covering of the Cap. In general the covering 419 // consists of at most 4 cells except for very large caps, which may need 420 // up to 6 cells. The output is not sorted. 421 func (c Cap) CellUnionBound() []CellID { 422 // TODO(roberts): The covering could be made quite a bit tighter by mapping 423 // the cap to a rectangle in (i,j)-space and finding a covering for that. 424 425 // Find the maximum level such that the cap contains at most one cell vertex 426 // and such that CellID.AppendVertexNeighbors() can be called. 427 level := MinWidthMetric.MaxLevel(c.Radius().Radians()) - 1 428 429 // If level < 0, more than three face cells are required. 430 if level < 0 { 431 cellIDs := make([]CellID, 6) 432 for face := 0; face < 6; face++ { 433 cellIDs[face] = CellIDFromFace(face) 434 } 435 return cellIDs 436 } 437 // The covering consists of the 4 cells at the given level that share the 438 // cell vertex that is closest to the cap center. 439 return cellIDFromPoint(c.center).VertexNeighbors(level) 440 } 441 442 // Centroid returns the true centroid of the cap multiplied by its surface area 443 // The result lies on the ray from the origin through the cap's center, but it 444 // is not unit length. Note that if you just want the "surface centroid", i.e. 445 // the normalized result, then it is simpler to call Center. 446 // 447 // The reason for multiplying the result by the cap area is to make it 448 // easier to compute the centroid of more complicated shapes. The centroid 449 // of a union of disjoint regions can be computed simply by adding their 450 // Centroid() results. Caveat: for caps that contain a single point 451 // (i.e., zero radius), this method always returns the origin (0, 0, 0). 452 // This is because shapes with no area don't affect the centroid of a 453 // union whose total area is positive. 454 func (c Cap) Centroid() Point { 455 // From symmetry, the centroid of the cap must be somewhere on the line 456 // from the origin to the center of the cap on the surface of the sphere. 457 // When a sphere is divided into slices of constant thickness by a set of 458 // parallel planes, all slices have the same surface area. This implies 459 // that the radial component of the centroid is simply the midpoint of the 460 // range of radial distances spanned by the cap. That is easily computed 461 // from the cap height. 462 if c.IsEmpty() { 463 return Point{} 464 } 465 r := 1 - 0.5*c.Height() 466 return Point{c.center.Mul(r * c.Area())} 467 } 468 469 // Union returns the smallest cap which encloses this cap and other. 470 func (c Cap) Union(other Cap) Cap { 471 // If the other cap is larger, swap c and other for the rest of the computations. 472 if c.radius < other.radius { 473 c, other = other, c 474 } 475 476 if c.IsFull() || other.IsEmpty() { 477 return c 478 } 479 480 // TODO: This calculation would be more efficient using s1.ChordAngles. 481 cRadius := c.Radius() 482 otherRadius := other.Radius() 483 distance := c.center.Distance(other.center) 484 if cRadius >= distance+otherRadius { 485 return c 486 } 487 488 resRadius := 0.5 * (distance + cRadius + otherRadius) 489 resCenter := InterpolateAtDistance(0.5*(distance-cRadius+otherRadius), c.center, other.center) 490 return CapFromCenterAngle(resCenter, resRadius) 491 } 492 493 // Encode encodes the Cap. 494 func (c Cap) Encode(w io.Writer) error { 495 e := &encoder{w: w} 496 c.encode(e) 497 return e.err 498 } 499 500 func (c Cap) encode(e *encoder) { 501 e.writeFloat64(c.center.X) 502 e.writeFloat64(c.center.Y) 503 e.writeFloat64(c.center.Z) 504 e.writeFloat64(float64(c.radius)) 505 } 506 507 // Decode decodes the Cap. 508 func (c *Cap) Decode(r io.Reader) error { 509 d := &decoder{r: asByteReader(r)} 510 c.decode(d) 511 return d.err 512 } 513 514 func (c *Cap) decode(d *decoder) { 515 c.center.X = d.readFloat64() 516 c.center.Y = d.readFloat64() 517 c.center.Z = d.readFloat64() 518 c.radius = s1.ChordAngle(d.readFloat64()) 519 } 520