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 "io" 19 "math" 20 21 "github.com/golang/geo/r1" 22 "github.com/golang/geo/r2" 23 "github.com/golang/geo/r3" 24 "github.com/golang/geo/s1" 25 ) 26 27 // Cell is an S2 region object that represents a cell. Unlike CellIDs, 28 // it supports efficient containment and intersection tests. However, it is 29 // also a more expensive representation. 30 type Cell struct { 31 face int8 32 level int8 33 orientation int8 34 id CellID 35 uv r2.Rect 36 } 37 38 // CellFromCellID constructs a Cell corresponding to the given CellID. 39 func CellFromCellID(id CellID) Cell { 40 c := Cell{} 41 c.id = id 42 f, i, j, o := c.id.faceIJOrientation() 43 c.face = int8(f) 44 c.level = int8(c.id.Level()) 45 c.orientation = int8(o) 46 c.uv = ijLevelToBoundUV(i, j, int(c.level)) 47 return c 48 } 49 50 // CellFromPoint constructs a cell for the given Point. 51 func CellFromPoint(p Point) Cell { 52 return CellFromCellID(cellIDFromPoint(p)) 53 } 54 55 // CellFromLatLng constructs a cell for the given LatLng. 56 func CellFromLatLng(ll LatLng) Cell { 57 return CellFromCellID(CellIDFromLatLng(ll)) 58 } 59 60 // Face returns the face this cell is on. 61 func (c Cell) Face() int { 62 return int(c.face) 63 } 64 65 // oppositeFace returns the face opposite the given face. 66 func oppositeFace(face int) int { 67 return (face + 3) % 6 68 } 69 70 // Level returns the level of this cell. 71 func (c Cell) Level() int { 72 return int(c.level) 73 } 74 75 // ID returns the CellID this cell represents. 76 func (c Cell) ID() CellID { 77 return c.id 78 } 79 80 // IsLeaf returns whether this Cell is a leaf or not. 81 func (c Cell) IsLeaf() bool { 82 return c.level == MaxLevel 83 } 84 85 // SizeIJ returns the edge length of this cell in (i,j)-space. 86 func (c Cell) SizeIJ() int { 87 return sizeIJ(int(c.level)) 88 } 89 90 // SizeST returns the edge length of this cell in (s,t)-space. 91 func (c Cell) SizeST() float64 { 92 return c.id.sizeST(int(c.level)) 93 } 94 95 // Vertex returns the k-th vertex of the cell (k = 0,1,2,3) in CCW order 96 // (lower left, lower right, upper right, upper left in the UV plane). 97 func (c Cell) Vertex(k int) Point { 98 return Point{faceUVToXYZ(int(c.face), c.uv.Vertices()[k].X, c.uv.Vertices()[k].Y).Normalize()} 99 } 100 101 // Edge returns the inward-facing normal of the great circle passing through 102 // the CCW ordered edge from vertex k to vertex k+1 (mod 4) (for k = 0,1,2,3). 103 func (c Cell) Edge(k int) Point { 104 switch k { 105 case 0: 106 return Point{vNorm(int(c.face), c.uv.Y.Lo).Normalize()} // Bottom 107 case 1: 108 return Point{uNorm(int(c.face), c.uv.X.Hi).Normalize()} // Right 109 case 2: 110 return Point{vNorm(int(c.face), c.uv.Y.Hi).Mul(-1.0).Normalize()} // Top 111 default: 112 return Point{uNorm(int(c.face), c.uv.X.Lo).Mul(-1.0).Normalize()} // Left 113 } 114 } 115 116 // BoundUV returns the bounds of this cell in (u,v)-space. 117 func (c Cell) BoundUV() r2.Rect { 118 return c.uv 119 } 120 121 // Center returns the direction vector corresponding to the center in 122 // (s,t)-space of the given cell. This is the point at which the cell is 123 // divided into four subcells; it is not necessarily the centroid of the 124 // cell in (u,v)-space or (x,y,z)-space 125 func (c Cell) Center() Point { 126 return Point{c.id.rawPoint().Normalize()} 127 } 128 129 // Children returns the four direct children of this cell in traversal order 130 // and returns true. If this is a leaf cell, or the children could not be created, 131 // false is returned. 132 // The C++ method is called Subdivide. 133 func (c Cell) Children() ([4]Cell, bool) { 134 var children [4]Cell 135 136 if c.id.IsLeaf() { 137 return children, false 138 } 139 140 // Compute the cell midpoint in uv-space. 141 uvMid := c.id.centerUV() 142 143 // Create four children with the appropriate bounds. 144 cid := c.id.ChildBegin() 145 for pos := 0; pos < 4; pos++ { 146 children[pos] = Cell{ 147 face: c.face, 148 level: c.level + 1, 149 orientation: c.orientation ^ int8(posToOrientation[pos]), 150 id: cid, 151 } 152 153 // We want to split the cell in half in u and v. To decide which 154 // side to set equal to the midpoint value, we look at cell's (i,j) 155 // position within its parent. The index for i is in bit 1 of ij. 156 ij := posToIJ[c.orientation][pos] 157 i := ij >> 1 158 j := ij & 1 159 if i == 1 { 160 children[pos].uv.X.Hi = c.uv.X.Hi 161 children[pos].uv.X.Lo = uvMid.X 162 } else { 163 children[pos].uv.X.Lo = c.uv.X.Lo 164 children[pos].uv.X.Hi = uvMid.X 165 } 166 if j == 1 { 167 children[pos].uv.Y.Hi = c.uv.Y.Hi 168 children[pos].uv.Y.Lo = uvMid.Y 169 } else { 170 children[pos].uv.Y.Lo = c.uv.Y.Lo 171 children[pos].uv.Y.Hi = uvMid.Y 172 } 173 cid = cid.Next() 174 } 175 return children, true 176 } 177 178 // ExactArea returns the area of this cell as accurately as possible. 179 func (c Cell) ExactArea() float64 { 180 v0, v1, v2, v3 := c.Vertex(0), c.Vertex(1), c.Vertex(2), c.Vertex(3) 181 return PointArea(v0, v1, v2) + PointArea(v0, v2, v3) 182 } 183 184 // ApproxArea returns the approximate area of this cell. This method is accurate 185 // to within 3% percent for all cell sizes and accurate to within 0.1% for cells 186 // at level 5 or higher (i.e. squares 350km to a side or smaller on the Earth's 187 // surface). It is moderately cheap to compute. 188 func (c Cell) ApproxArea() float64 { 189 // All cells at the first two levels have the same area. 190 if c.level < 2 { 191 return c.AverageArea() 192 } 193 194 // First, compute the approximate area of the cell when projected 195 // perpendicular to its normal. The cross product of its diagonals gives 196 // the normal, and the length of the normal is twice the projected area. 197 flatArea := 0.5 * (c.Vertex(2).Sub(c.Vertex(0).Vector). 198 Cross(c.Vertex(3).Sub(c.Vertex(1).Vector)).Norm()) 199 200 // Now, compensate for the curvature of the cell surface by pretending 201 // that the cell is shaped like a spherical cap. The ratio of the 202 // area of a spherical cap to the area of its projected disc turns out 203 // to be 2 / (1 + sqrt(1 - r*r)) where r is the radius of the disc. 204 // For example, when r=0 the ratio is 1, and when r=1 the ratio is 2. 205 // Here we set Pi*r*r == flatArea to find the equivalent disc. 206 return flatArea * 2 / (1 + math.Sqrt(1-math.Min(1/math.Pi*flatArea, 1))) 207 } 208 209 // AverageArea returns the average area of cells at the level of this cell. 210 // This is accurate to within a factor of 1.7. 211 func (c Cell) AverageArea() float64 { 212 return AvgAreaMetric.Value(int(c.level)) 213 } 214 215 // IntersectsCell reports whether the intersection of this cell and the other cell is not nil. 216 func (c Cell) IntersectsCell(oc Cell) bool { 217 return c.id.Intersects(oc.id) 218 } 219 220 // ContainsCell reports whether this cell contains the other cell. 221 func (c Cell) ContainsCell(oc Cell) bool { 222 return c.id.Contains(oc.id) 223 } 224 225 // CellUnionBound computes a covering of the Cell. 226 func (c Cell) CellUnionBound() []CellID { 227 return c.CapBound().CellUnionBound() 228 } 229 230 // latitude returns the latitude of the cell vertex in radians given by (i,j), 231 // where i and j indicate the Hi (1) or Lo (0) corner. 232 func (c Cell) latitude(i, j int) float64 { 233 var u, v float64 234 switch { 235 case i == 0 && j == 0: 236 u = c.uv.X.Lo 237 v = c.uv.Y.Lo 238 case i == 0 && j == 1: 239 u = c.uv.X.Lo 240 v = c.uv.Y.Hi 241 case i == 1 && j == 0: 242 u = c.uv.X.Hi 243 v = c.uv.Y.Lo 244 case i == 1 && j == 1: 245 u = c.uv.X.Hi 246 v = c.uv.Y.Hi 247 default: 248 panic("i and/or j is out of bounds") 249 } 250 return latitude(Point{faceUVToXYZ(int(c.face), u, v)}).Radians() 251 } 252 253 // longitude returns the longitude of the cell vertex in radians given by (i,j), 254 // where i and j indicate the Hi (1) or Lo (0) corner. 255 func (c Cell) longitude(i, j int) float64 { 256 var u, v float64 257 switch { 258 case i == 0 && j == 0: 259 u = c.uv.X.Lo 260 v = c.uv.Y.Lo 261 case i == 0 && j == 1: 262 u = c.uv.X.Lo 263 v = c.uv.Y.Hi 264 case i == 1 && j == 0: 265 u = c.uv.X.Hi 266 v = c.uv.Y.Lo 267 case i == 1 && j == 1: 268 u = c.uv.X.Hi 269 v = c.uv.Y.Hi 270 default: 271 panic("i and/or j is out of bounds") 272 } 273 return longitude(Point{faceUVToXYZ(int(c.face), u, v)}).Radians() 274 } 275 276 var ( 277 poleMinLat = math.Asin(math.Sqrt(1.0/3)) - 0.5*dblEpsilon 278 ) 279 280 // RectBound returns the bounding rectangle of this cell. 281 func (c Cell) RectBound() Rect { 282 if c.level > 0 { 283 // Except for cells at level 0, the latitude and longitude extremes are 284 // attained at the vertices. Furthermore, the latitude range is 285 // determined by one pair of diagonally opposite vertices and the 286 // longitude range is determined by the other pair. 287 // 288 // We first determine which corner (i,j) of the cell has the largest 289 // absolute latitude. To maximize latitude, we want to find the point in 290 // the cell that has the largest absolute z-coordinate and the smallest 291 // absolute x- and y-coordinates. To do this we look at each coordinate 292 // (u and v), and determine whether we want to minimize or maximize that 293 // coordinate based on the axis direction and the cell's (u,v) quadrant. 294 u := c.uv.X.Lo + c.uv.X.Hi 295 v := c.uv.Y.Lo + c.uv.Y.Hi 296 var i, j int 297 if uAxis(int(c.face)).Z == 0 { 298 if u < 0 { 299 i = 1 300 } 301 } else if u > 0 { 302 i = 1 303 } 304 if vAxis(int(c.face)).Z == 0 { 305 if v < 0 { 306 j = 1 307 } 308 } else if v > 0 { 309 j = 1 310 } 311 lat := r1.IntervalFromPoint(c.latitude(i, j)).AddPoint(c.latitude(1-i, 1-j)) 312 lng := s1.EmptyInterval().AddPoint(c.longitude(i, 1-j)).AddPoint(c.longitude(1-i, j)) 313 314 // We grow the bounds slightly to make sure that the bounding rectangle 315 // contains LatLngFromPoint(P) for any point P inside the loop L defined by the 316 // four *normalized* vertices. Note that normalization of a vector can 317 // change its direction by up to 0.5 * dblEpsilon radians, and it is not 318 // enough just to add Normalize calls to the code above because the 319 // latitude/longitude ranges are not necessarily determined by diagonally 320 // opposite vertex pairs after normalization. 321 // 322 // We would like to bound the amount by which the latitude/longitude of a 323 // contained point P can exceed the bounds computed above. In the case of 324 // longitude, the normalization error can change the direction of rounding 325 // leading to a maximum difference in longitude of 2 * dblEpsilon. In 326 // the case of latitude, the normalization error can shift the latitude by 327 // up to 0.5 * dblEpsilon and the other sources of error can cause the 328 // two latitudes to differ by up to another 1.5 * dblEpsilon, which also 329 // leads to a maximum difference of 2 * dblEpsilon. 330 return Rect{lat, lng}.expanded(LatLng{s1.Angle(2 * dblEpsilon), s1.Angle(2 * dblEpsilon)}).PolarClosure() 331 } 332 333 // The 4 cells around the equator extend to +/-45 degrees latitude at the 334 // midpoints of their top and bottom edges. The two cells covering the 335 // poles extend down to +/-35.26 degrees at their vertices. The maximum 336 // error in this calculation is 0.5 * dblEpsilon. 337 var bound Rect 338 switch c.face { 339 case 0: 340 bound = Rect{r1.Interval{-math.Pi / 4, math.Pi / 4}, s1.Interval{-math.Pi / 4, math.Pi / 4}} 341 case 1: 342 bound = Rect{r1.Interval{-math.Pi / 4, math.Pi / 4}, s1.Interval{math.Pi / 4, 3 * math.Pi / 4}} 343 case 2: 344 bound = Rect{r1.Interval{poleMinLat, math.Pi / 2}, s1.FullInterval()} 345 case 3: 346 bound = Rect{r1.Interval{-math.Pi / 4, math.Pi / 4}, s1.Interval{3 * math.Pi / 4, -3 * math.Pi / 4}} 347 case 4: 348 bound = Rect{r1.Interval{-math.Pi / 4, math.Pi / 4}, s1.Interval{-3 * math.Pi / 4, -math.Pi / 4}} 349 default: 350 bound = Rect{r1.Interval{-math.Pi / 2, -poleMinLat}, s1.FullInterval()} 351 } 352 353 // Finally, we expand the bound to account for the error when a point P is 354 // converted to an LatLng to test for containment. (The bound should be 355 // large enough so that it contains the computed LatLng of any contained 356 // point, not just the infinite-precision version.) We don't need to expand 357 // longitude because longitude is calculated via a single call to math.Atan2, 358 // which is guaranteed to be semi-monotonic. 359 return bound.expanded(LatLng{s1.Angle(dblEpsilon), s1.Angle(0)}) 360 } 361 362 // CapBound returns the bounding cap of this cell. 363 func (c Cell) CapBound() Cap { 364 // We use the cell center in (u,v)-space as the cap axis. This vector is very close 365 // to GetCenter() and faster to compute. Neither one of these vectors yields the 366 // bounding cap with minimal surface area, but they are both pretty close. 367 cap := CapFromPoint(Point{faceUVToXYZ(int(c.face), c.uv.Center().X, c.uv.Center().Y).Normalize()}) 368 for k := 0; k < 4; k++ { 369 cap = cap.AddPoint(c.Vertex(k)) 370 } 371 return cap 372 } 373 374 // ContainsPoint reports whether this cell contains the given point. Note that 375 // unlike Loop/Polygon, a Cell is considered to be a closed set. This means 376 // that a point on a Cell's edge or vertex belong to the Cell and the relevant 377 // adjacent Cells too. 378 // 379 // If you want every point to be contained by exactly one Cell, 380 // you will need to convert the Cell to a Loop. 381 func (c Cell) ContainsPoint(p Point) bool { 382 // We can't just call XYZtoFaceUV, because for points that lie on the 383 // boundary between two faces (i.e. u or v is +1/-1) we need to return 384 // true for both adjacent cells. 385 // 386 // We can get away with not checking the face if the point matches the face of 387 // the cell here because, for the 4 faces adjacent to c.face, p will be 388 // projected outside the range of ([-1,1]x[-1,1]) and thus can't intersect the 389 // cell bounds (except on the face boundary which we want). 390 // 391 // For the face opposite c.face, the sign of the UV coordinates of P will be 392 // flipped so it will automatically fall outside the cell boundary as no cells 393 // cross the origin. 394 var uv r2.Point 395 var ok bool 396 if uv.X, uv.Y, ok = faceXYZToUV(int(c.face), p); !ok { 397 return false 398 } 399 400 // Expand the (u,v) bound to ensure that 401 // 402 // CellFromPoint(p).ContainsPoint(p) 403 // 404 // is always true. To do this, we need to account for the error when 405 // converting from (u,v) coordinates to (s,t) coordinates. In the 406 // normal case the total error is at most dblEpsilon. 407 return c.uv.ExpandedByMargin(dblEpsilon).ContainsPoint(uv) 408 } 409 410 // Encode encodes the Cell. 411 func (c Cell) Encode(w io.Writer) error { 412 e := &encoder{w: w} 413 c.encode(e) 414 return e.err 415 } 416 417 func (c Cell) encode(e *encoder) { 418 c.id.encode(e) 419 } 420 421 // Decode decodes the Cell. 422 func (c *Cell) Decode(r io.Reader) error { 423 d := &decoder{r: asByteReader(r)} 424 c.decode(d) 425 return d.err 426 } 427 428 func (c *Cell) decode(d *decoder) { 429 c.id.decode(d) 430 *c = CellFromCellID(c.id) 431 } 432 433 // vertexChordDist2 returns the squared chord distance from point P to the 434 // given corner vertex specified by the Hi or Lo values of each. 435 func (c Cell) vertexChordDist2(p Point, xHi, yHi bool) s1.ChordAngle { 436 x := c.uv.X.Lo 437 y := c.uv.Y.Lo 438 if xHi { 439 x = c.uv.X.Hi 440 } 441 if yHi { 442 y = c.uv.Y.Hi 443 } 444 445 return ChordAngleBetweenPoints(p, PointFromCoords(x, y, 1)) 446 } 447 448 // uEdgeIsClosest reports whether a point P is closer to the interior of the specified 449 // Cell edge (either the lower or upper edge of the Cell) or to the endpoints. 450 func (c Cell) uEdgeIsClosest(p Point, vHi bool) bool { 451 u0 := c.uv.X.Lo 452 u1 := c.uv.X.Hi 453 v := c.uv.Y.Lo 454 if vHi { 455 v = c.uv.Y.Hi 456 } 457 // These are the normals to the planes that are perpendicular to the edge 458 // and pass through one of its two endpoints. 459 dir0 := r3.Vector{v*v + 1, -u0 * v, -u0} 460 dir1 := r3.Vector{v*v + 1, -u1 * v, -u1} 461 return p.Dot(dir0) > 0 && p.Dot(dir1) < 0 462 } 463 464 // vEdgeIsClosest reports whether a point P is closer to the interior of the specified 465 // Cell edge (either the right or left edge of the Cell) or to the endpoints. 466 func (c Cell) vEdgeIsClosest(p Point, uHi bool) bool { 467 v0 := c.uv.Y.Lo 468 v1 := c.uv.Y.Hi 469 u := c.uv.X.Lo 470 if uHi { 471 u = c.uv.X.Hi 472 } 473 dir0 := r3.Vector{-u * v0, u*u + 1, -v0} 474 dir1 := r3.Vector{-u * v1, u*u + 1, -v1} 475 return p.Dot(dir0) > 0 && p.Dot(dir1) < 0 476 } 477 478 // edgeDistance reports the distance from a Point P to a given Cell edge. The point 479 // P is given by its dot product, and the uv edge by its normal in the 480 // given coordinate value. 481 func edgeDistance(ij, uv float64) s1.ChordAngle { 482 // Let P by the target point and let R be the closest point on the given 483 // edge AB. The desired distance PR can be expressed as PR^2 = PQ^2 + QR^2 484 // where Q is the point P projected onto the plane through the great circle 485 // through AB. We can compute the distance PQ^2 perpendicular to the plane 486 // from "dirIJ" (the dot product of the target point P with the edge 487 // normal) and the squared length the edge normal (1 + uv**2). 488 pq2 := (ij * ij) / (1 + uv*uv) 489 490 // We can compute the distance QR as (1 - OQ) where O is the sphere origin, 491 // and we can compute OQ^2 = 1 - PQ^2 using the Pythagorean theorem. 492 // (This calculation loses accuracy as angle POQ approaches Pi/2.) 493 qr := 1 - math.Sqrt(1-pq2) 494 return s1.ChordAngleFromSquaredLength(pq2 + qr*qr) 495 } 496 497 // distanceInternal reports the distance from the given point to the interior of 498 // the cell if toInterior is true or to the boundary of the cell otherwise. 499 func (c Cell) distanceInternal(targetXYZ Point, toInterior bool) s1.ChordAngle { 500 // All calculations are done in the (u,v,w) coordinates of this cell's face. 501 target := faceXYZtoUVW(int(c.face), targetXYZ) 502 503 // Compute dot products with all four upward or rightward-facing edge 504 // normals. dirIJ is the dot product for the edge corresponding to axis 505 // I, endpoint J. For example, dir01 is the right edge of the Cell 506 // (corresponding to the upper endpoint of the u-axis). 507 dir00 := target.X - target.Z*c.uv.X.Lo 508 dir01 := target.X - target.Z*c.uv.X.Hi 509 dir10 := target.Y - target.Z*c.uv.Y.Lo 510 dir11 := target.Y - target.Z*c.uv.Y.Hi 511 inside := true 512 if dir00 < 0 { 513 inside = false // Target is to the left of the cell 514 if c.vEdgeIsClosest(target, false) { 515 return edgeDistance(-dir00, c.uv.X.Lo) 516 } 517 } 518 if dir01 > 0 { 519 inside = false // Target is to the right of the cell 520 if c.vEdgeIsClosest(target, true) { 521 return edgeDistance(dir01, c.uv.X.Hi) 522 } 523 } 524 if dir10 < 0 { 525 inside = false // Target is below the cell 526 if c.uEdgeIsClosest(target, false) { 527 return edgeDistance(-dir10, c.uv.Y.Lo) 528 } 529 } 530 if dir11 > 0 { 531 inside = false // Target is above the cell 532 if c.uEdgeIsClosest(target, true) { 533 return edgeDistance(dir11, c.uv.Y.Hi) 534 } 535 } 536 if inside { 537 if toInterior { 538 return s1.ChordAngle(0) 539 } 540 // Although you might think of Cells as rectangles, they are actually 541 // arbitrary quadrilaterals after they are projected onto the sphere. 542 // Therefore the simplest approach is just to find the minimum distance to 543 // any of the four edges. 544 return minChordAngle(edgeDistance(-dir00, c.uv.X.Lo), 545 edgeDistance(dir01, c.uv.X.Hi), 546 edgeDistance(-dir10, c.uv.Y.Lo), 547 edgeDistance(dir11, c.uv.Y.Hi)) 548 } 549 550 // Otherwise, the closest point is one of the four cell vertices. Note that 551 // it is *not* trivial to narrow down the candidates based on the edge sign 552 // tests above, because (1) the edges don't meet at right angles and (2) 553 // there are points on the far side of the sphere that are both above *and* 554 // below the cell, etc. 555 return minChordAngle(c.vertexChordDist2(target, false, false), 556 c.vertexChordDist2(target, true, false), 557 c.vertexChordDist2(target, false, true), 558 c.vertexChordDist2(target, true, true)) 559 } 560 561 // Distance reports the distance from the cell to the given point. Returns zero if 562 // the point is inside the cell. 563 func (c Cell) Distance(target Point) s1.ChordAngle { 564 return c.distanceInternal(target, true) 565 } 566 567 // MaxDistance reports the maximum distance from the cell (including its interior) to the 568 // given point. 569 func (c Cell) MaxDistance(target Point) s1.ChordAngle { 570 // First check the 4 cell vertices. If all are within the hemisphere 571 // centered around target, the max distance will be to one of these vertices. 572 targetUVW := faceXYZtoUVW(int(c.face), target) 573 maxDist := maxChordAngle(c.vertexChordDist2(targetUVW, false, false), 574 c.vertexChordDist2(targetUVW, true, false), 575 c.vertexChordDist2(targetUVW, false, true), 576 c.vertexChordDist2(targetUVW, true, true)) 577 578 if maxDist <= s1.RightChordAngle { 579 return maxDist 580 } 581 582 // Otherwise, find the minimum distance dMin to the antipodal point and the 583 // maximum distance will be pi - dMin. 584 return s1.StraightChordAngle - c.Distance(Point{target.Mul(-1)}) 585 } 586 587 // BoundaryDistance reports the distance from the cell boundary to the given point. 588 func (c Cell) BoundaryDistance(target Point) s1.ChordAngle { 589 return c.distanceInternal(target, false) 590 } 591 592 // DistanceToEdge returns the minimum distance from the cell to the given edge AB. Returns 593 // zero if the edge intersects the cell interior. 594 func (c Cell) DistanceToEdge(a, b Point) s1.ChordAngle { 595 // Possible optimizations: 596 // - Currently the (cell vertex, edge endpoint) distances are computed 597 // twice each, and the length of AB is computed 4 times. 598 // - To fix this, refactor GetDistance(target) so that it skips calculating 599 // the distance to each cell vertex. Instead, compute the cell vertices 600 // and distances in this function, and add a low-level UpdateMinDistance 601 // that allows the XA, XB, and AB distances to be passed in. 602 // - It might also be more efficient to do all calculations in UVW-space, 603 // since this would involve transforming 2 points rather than 4. 604 605 // First, check the minimum distance to the edge endpoints A and B. 606 // (This also detects whether either endpoint is inside the cell.) 607 minDist := minChordAngle(c.Distance(a), c.Distance(b)) 608 if minDist == 0 { 609 return minDist 610 } 611 612 // Otherwise, check whether the edge crosses the cell boundary. 613 crosser := NewChainEdgeCrosser(a, b, c.Vertex(3)) 614 for i := 0; i < 4; i++ { 615 if crosser.ChainCrossingSign(c.Vertex(i)) != DoNotCross { 616 return 0 617 } 618 } 619 620 // Finally, check whether the minimum distance occurs between a cell vertex 621 // and the interior of the edge AB. (Some of this work is redundant, since 622 // it also checks the distance to the endpoints A and B again.) 623 // 624 // Note that we don't need to check the distance from the interior of AB to 625 // the interior of a cell edge, because the only way that this distance can 626 // be minimal is if the two edges cross (already checked above). 627 for i := 0; i < 4; i++ { 628 minDist, _ = UpdateMinDistance(c.Vertex(i), a, b, minDist) 629 } 630 return minDist 631 } 632 633 // MaxDistanceToEdge returns the maximum distance from the cell (including its interior) 634 // to the given edge AB. 635 func (c Cell) MaxDistanceToEdge(a, b Point) s1.ChordAngle { 636 // If the maximum distance from both endpoints to the cell is less than π/2 637 // then the maximum distance from the edge to the cell is the maximum of the 638 // two endpoint distances. 639 maxDist := maxChordAngle(c.MaxDistance(a), c.MaxDistance(b)) 640 if maxDist <= s1.RightChordAngle { 641 return maxDist 642 } 643 644 return s1.StraightChordAngle - c.DistanceToEdge(Point{a.Mul(-1)}, Point{b.Mul(-1)}) 645 } 646 647 // DistanceToCell returns the minimum distance from this cell to the given cell. 648 // It returns zero if one cell contains the other. 649 func (c Cell) DistanceToCell(target Cell) s1.ChordAngle { 650 // If the cells intersect, the distance is zero. We use the (u,v) ranges 651 // rather than CellID intersects so that cells that share a partial edge or 652 // corner are considered to intersect. 653 if c.face == target.face && c.uv.Intersects(target.uv) { 654 return 0 655 } 656 657 // Otherwise, the minimum distance always occurs between a vertex of one 658 // cell and an edge of the other cell (including the edge endpoints). This 659 // represents a total of 32 possible (vertex, edge) pairs. 660 // 661 // TODO(roberts): This could be optimized to be at least 5x faster by pruning 662 // the set of possible closest vertex/edge pairs using the faces and (u,v) 663 // ranges of both cells. 664 var va, vb [4]Point 665 for i := 0; i < 4; i++ { 666 va[i] = c.Vertex(i) 667 vb[i] = target.Vertex(i) 668 } 669 minDist := s1.InfChordAngle() 670 for i := 0; i < 4; i++ { 671 for j := 0; j < 4; j++ { 672 minDist, _ = UpdateMinDistance(va[i], vb[j], vb[(j+1)&3], minDist) 673 minDist, _ = UpdateMinDistance(vb[i], va[j], va[(j+1)&3], minDist) 674 } 675 } 676 return minDist 677 } 678 679 // MaxDistanceToCell returns the maximum distance from the cell (including its 680 // interior) to the given target cell. 681 func (c Cell) MaxDistanceToCell(target Cell) s1.ChordAngle { 682 // Need to check the antipodal target for intersection with the cell. If it 683 // intersects, the distance is the straight ChordAngle. 684 // antipodalUV is the transpose of the original UV, interpreted within the opposite face. 685 antipodalUV := r2.Rect{target.uv.Y, target.uv.X} 686 if int(c.face) == oppositeFace(int(target.face)) && c.uv.Intersects(antipodalUV) { 687 return s1.StraightChordAngle 688 } 689 690 // Otherwise, the maximum distance always occurs between a vertex of one 691 // cell and an edge of the other cell (including the edge endpoints). This 692 // represents a total of 32 possible (vertex, edge) pairs. 693 // 694 // TODO(roberts): When the maximum distance is at most π/2, the maximum is 695 // always attained between a pair of vertices, and this could be made much 696 // faster by testing each vertex pair once rather than the current 4 times. 697 var va, vb [4]Point 698 for i := 0; i < 4; i++ { 699 va[i] = c.Vertex(i) 700 vb[i] = target.Vertex(i) 701 } 702 maxDist := s1.NegativeChordAngle 703 for i := 0; i < 4; i++ { 704 for j := 0; j < 4; j++ { 705 maxDist, _ = UpdateMaxDistance(va[i], vb[j], vb[(j+1)&3], maxDist) 706 maxDist, _ = UpdateMaxDistance(vb[i], va[j], va[(j+1)&3], maxDist) 707 } 708 } 709 return maxDist 710 } 711