1 // Copyright 2017 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 // This file contains a collection of methods for: 18 // 19 // (1) Robustly clipping geodesic edges to the faces of the S2 biunit cube 20 // (see s2stuv), and 21 // 22 // (2) Robustly clipping 2D edges against 2D rectangles. 23 // 24 // These functions can be used to efficiently find the set of CellIDs that 25 // are intersected by a geodesic edge (e.g., see CrossingEdgeQuery). 26 27 import ( 28 "math" 29 30 "github.com/golang/geo/r1" 31 "github.com/golang/geo/r2" 32 "github.com/golang/geo/r3" 33 ) 34 35 const ( 36 // edgeClipErrorUVCoord is the maximum error in a u- or v-coordinate 37 // compared to the exact result, assuming that the points A and B are in 38 // the rectangle [-1,1]x[1,1] or slightly outside it (by 1e-10 or less). 39 edgeClipErrorUVCoord = 2.25 * dblEpsilon 40 41 // edgeClipErrorUVDist is the maximum distance from a clipped point to 42 // the corresponding exact result. It is equal to the error in a single 43 // coordinate because at most one coordinate is subject to error. 44 edgeClipErrorUVDist = 2.25 * dblEpsilon 45 46 // faceClipErrorRadians is the maximum angle between a returned vertex 47 // and the nearest point on the exact edge AB. It is equal to the 48 // maximum directional error in PointCross, plus the error when 49 // projecting points onto a cube face. 50 faceClipErrorRadians = 3 * dblEpsilon 51 52 // faceClipErrorDist is the same angle expressed as a maximum distance 53 // in (u,v)-space. In other words, a returned vertex is at most this far 54 // from the exact edge AB projected into (u,v)-space. 55 faceClipErrorUVDist = 9 * dblEpsilon 56 57 // faceClipErrorUVCoord is the maximum angle between a returned vertex 58 // and the nearest point on the exact edge AB expressed as the maximum error 59 // in an individual u- or v-coordinate. In other words, for each 60 // returned vertex there is a point on the exact edge AB whose u- and 61 // v-coordinates differ from the vertex by at most this amount. 62 faceClipErrorUVCoord = 9.0 * (1.0 / math.Sqrt2) * dblEpsilon 63 64 // intersectsRectErrorUVDist is the maximum error when computing if a point 65 // intersects with a given Rect. If some point of AB is inside the 66 // rectangle by at least this distance, the result is guaranteed to be true; 67 // if all points of AB are outside the rectangle by at least this distance, 68 // the result is guaranteed to be false. This bound assumes that rect is 69 // a subset of the rectangle [-1,1]x[-1,1] or extends slightly outside it 70 // (e.g., by 1e-10 or less). 71 intersectsRectErrorUVDist = 3 * math.Sqrt2 * dblEpsilon 72 ) 73 74 // ClipToFace returns the (u,v) coordinates for the portion of the edge AB that 75 // intersects the given face, or false if the edge AB does not intersect. 76 // This method guarantees that the clipped vertices lie within the [-1,1]x[-1,1] 77 // cube face rectangle and are within faceClipErrorUVDist of the line AB, but 78 // the results may differ from those produced by FaceSegments. 79 func ClipToFace(a, b Point, face int) (aUV, bUV r2.Point, intersects bool) { 80 return ClipToPaddedFace(a, b, face, 0.0) 81 } 82 83 // ClipToPaddedFace returns the (u,v) coordinates for the portion of the edge AB that 84 // intersects the given face, but rather than clipping to the square [-1,1]x[-1,1] 85 // in (u,v) space, this method clips to [-R,R]x[-R,R] where R=(1+padding). 86 // Padding must be non-negative. 87 func ClipToPaddedFace(a, b Point, f int, padding float64) (aUV, bUV r2.Point, intersects bool) { 88 // Fast path: both endpoints are on the given face. 89 if face(a.Vector) == f && face(b.Vector) == f { 90 au, av := validFaceXYZToUV(f, a.Vector) 91 bu, bv := validFaceXYZToUV(f, b.Vector) 92 return r2.Point{au, av}, r2.Point{bu, bv}, true 93 } 94 95 // Convert everything into the (u,v,w) coordinates of the given face. Note 96 // that the cross product *must* be computed in the original (x,y,z) 97 // coordinate system because PointCross (unlike the mathematical cross 98 // product) can produce different results in different coordinate systems 99 // when one argument is a linear multiple of the other, due to the use of 100 // symbolic perturbations. 101 normUVW := pointUVW(faceXYZtoUVW(f, a.PointCross(b))) 102 aUVW := pointUVW(faceXYZtoUVW(f, a)) 103 bUVW := pointUVW(faceXYZtoUVW(f, b)) 104 105 // Padding is handled by scaling the u- and v-components of the normal. 106 // Letting R=1+padding, this means that when we compute the dot product of 107 // the normal with a cube face vertex (such as (-1,-1,1)), we will actually 108 // compute the dot product with the scaled vertex (-R,-R,1). This allows 109 // methods such as intersectsFace, exitAxis, etc, to handle padding 110 // with no further modifications. 111 scaleUV := 1 + padding 112 scaledN := pointUVW{r3.Vector{X: scaleUV * normUVW.X, Y: scaleUV * normUVW.Y, Z: normUVW.Z}} 113 if !scaledN.intersectsFace() { 114 return aUV, bUV, false 115 } 116 117 // TODO(roberts): This is a workaround for extremely small vectors where some 118 // loss of precision can occur in Normalize causing underflow. When PointCross 119 // is updated to work around this, this can be removed. 120 if math.Max(math.Abs(normUVW.X), math.Max(math.Abs(normUVW.Y), math.Abs(normUVW.Z))) < math.Ldexp(1, -511) { 121 normUVW = pointUVW{normUVW.Mul(math.Ldexp(1, 563))} 122 } 123 124 normUVW = pointUVW{normUVW.Normalize()} 125 126 aTan := pointUVW{normUVW.Cross(aUVW.Vector)} 127 bTan := pointUVW{bUVW.Cross(normUVW.Vector)} 128 129 // As described in clipDestination, if the sum of the scores from clipping the two 130 // endpoints is 3 or more, then the segment does not intersect this face. 131 aUV, aScore := clipDestination(bUVW, aUVW, pointUVW{scaledN.Mul(-1)}, bTan, aTan, scaleUV) 132 bUV, bScore := clipDestination(aUVW, bUVW, scaledN, aTan, bTan, scaleUV) 133 134 return aUV, bUV, aScore+bScore < 3 135 } 136 137 // ClipEdge returns the portion of the edge defined by AB that is contained by the 138 // given rectangle. If there is no intersection, false is returned and aClip and bClip 139 // are undefined. 140 func ClipEdge(a, b r2.Point, clip r2.Rect) (aClip, bClip r2.Point, intersects bool) { 141 // Compute the bounding rectangle of AB, clip it, and then extract the new 142 // endpoints from the clipped bound. 143 bound := r2.RectFromPoints(a, b) 144 if bound, intersects = clipEdgeBound(a, b, clip, bound); !intersects { 145 return aClip, bClip, false 146 } 147 ai := 0 148 if a.X > b.X { 149 ai = 1 150 } 151 aj := 0 152 if a.Y > b.Y { 153 aj = 1 154 } 155 156 return bound.VertexIJ(ai, aj), bound.VertexIJ(1-ai, 1-aj), true 157 } 158 159 // The three functions below (sumEqual, intersectsFace, intersectsOppositeEdges) 160 // all compare a sum (u + v) to a third value w. They are implemented in such a 161 // way that they produce an exact result even though all calculations are done 162 // with ordinary floating-point operations. Here are the principles on which these 163 // functions are based: 164 // 165 // A. If u + v < w in floating-point, then u + v < w in exact arithmetic. 166 // 167 // B. If u + v < w in exact arithmetic, then at least one of the following 168 // expressions is true in floating-point: 169 // u + v < w 170 // u < w - v 171 // v < w - u 172 // 173 // Proof: By rearranging terms and substituting ">" for "<", we can assume 174 // that all values are non-negative. Now clearly "w" is not the smallest 175 // value, so assume WLOG that "u" is the smallest. We want to show that 176 // u < w - v in floating-point. If v >= w/2, the calculation of w - v is 177 // exact since the result is smaller in magnitude than either input value, 178 // so the result holds. Otherwise we have u <= v < w/2 and w - v >= w/2 179 // (even in floating point), so the result also holds. 180 181 // sumEqual reports whether u + v == w exactly. 182 func sumEqual(u, v, w float64) bool { 183 return (u+v == w) && (u == w-v) && (v == w-u) 184 } 185 186 // pointUVW represents a Point in (u,v,w) coordinate space of a cube face. 187 type pointUVW Point 188 189 // intersectsFace reports whether a given directed line L intersects the cube face F. 190 // The line L is defined by its normal N in the (u,v,w) coordinates of F. 191 func (p pointUVW) intersectsFace() bool { 192 // L intersects the [-1,1]x[-1,1] square in (u,v) if and only if the dot 193 // products of N with the four corner vertices (-1,-1,1), (1,-1,1), (1,1,1), 194 // and (-1,1,1) do not all have the same sign. This is true exactly when 195 // |Nu| + |Nv| >= |Nw|. The code below evaluates this expression exactly. 196 u := math.Abs(p.X) 197 v := math.Abs(p.Y) 198 w := math.Abs(p.Z) 199 200 // We only need to consider the cases where u or v is the smallest value, 201 // since if w is the smallest then both expressions below will have a 202 // positive LHS and a negative RHS. 203 return (v >= w-u) && (u >= w-v) 204 } 205 206 // intersectsOppositeEdges reports whether a directed line L intersects two 207 // opposite edges of a cube face F. This includes the case where L passes 208 // exactly through a corner vertex of F. The directed line L is defined 209 // by its normal N in the (u,v,w) coordinates of F. 210 func (p pointUVW) intersectsOppositeEdges() bool { 211 // The line L intersects opposite edges of the [-1,1]x[-1,1] (u,v) square if 212 // and only exactly two of the corner vertices lie on each side of L. This 213 // is true exactly when ||Nu| - |Nv|| >= |Nw|. The code below evaluates this 214 // expression exactly. 215 u := math.Abs(p.X) 216 v := math.Abs(p.Y) 217 w := math.Abs(p.Z) 218 219 // If w is the smallest, the following line returns an exact result. 220 if math.Abs(u-v) != w { 221 return math.Abs(u-v) >= w 222 } 223 224 // Otherwise u - v = w exactly, or w is not the smallest value. In either 225 // case the following returns the correct result. 226 if u >= v { 227 return u-w >= v 228 } 229 return v-w >= u 230 } 231 232 // axis represents the possible results of exitAxis. 233 type axis int 234 235 const ( 236 axisU axis = iota 237 axisV 238 ) 239 240 // exitAxis reports which axis the directed line L exits the cube face F on. 241 // The directed line L is represented by its CCW normal N in the (u,v,w) coordinates 242 // of F. It returns axisU if L exits through the u=-1 or u=+1 edge, and axisV if L exits 243 // through the v=-1 or v=+1 edge. Either result is acceptable if L exits exactly 244 // through a corner vertex of the cube face. 245 func (p pointUVW) exitAxis() axis { 246 if p.intersectsOppositeEdges() { 247 // The line passes through through opposite edges of the face. 248 // It exits through the v=+1 or v=-1 edge if the u-component of N has a 249 // larger absolute magnitude than the v-component. 250 if math.Abs(p.X) >= math.Abs(p.Y) { 251 return axisV 252 } 253 return axisU 254 } 255 256 // The line passes through through two adjacent edges of the face. 257 // It exits the v=+1 or v=-1 edge if an even number of the components of N 258 // are negative. We test this using signbit() rather than multiplication 259 // to avoid the possibility of underflow. 260 var x, y, z int 261 if math.Signbit(p.X) { 262 x = 1 263 } 264 if math.Signbit(p.Y) { 265 y = 1 266 } 267 if math.Signbit(p.Z) { 268 z = 1 269 } 270 271 if x^y^z == 0 { 272 return axisV 273 } 274 return axisU 275 } 276 277 // exitPoint returns the UV coordinates of the point where a directed line L (represented 278 // by the CCW normal of this point), exits the cube face this point is derived from along 279 // the given axis. 280 func (p pointUVW) exitPoint(a axis) r2.Point { 281 if a == axisU { 282 u := -1.0 283 if p.Y > 0 { 284 u = 1.0 285 } 286 return r2.Point{u, (-u*p.X - p.Z) / p.Y} 287 } 288 289 v := -1.0 290 if p.X < 0 { 291 v = 1.0 292 } 293 return r2.Point{(-v*p.Y - p.Z) / p.X, v} 294 } 295 296 // clipDestination returns a score which is used to indicate if the clipped edge AB 297 // on the given face intersects the face at all. This function returns the score for 298 // the given endpoint, which is an integer ranging from 0 to 3. If the sum of the scores 299 // from both of the endpoints is 3 or more, then edge AB does not intersect this face. 300 // 301 // First, it clips the line segment AB to find the clipped destination B' on a given 302 // face. (The face is specified implicitly by expressing *all arguments* in the (u,v,w) 303 // coordinates of that face.) Second, it partially computes whether the segment AB 304 // intersects this face at all. The actual condition is fairly complicated, but it 305 // turns out that it can be expressed as a "score" that can be computed independently 306 // when clipping the two endpoints A and B. 307 func clipDestination(a, b, scaledN, aTan, bTan pointUVW, scaleUV float64) (r2.Point, int) { 308 var uv r2.Point 309 310 // Optimization: if B is within the safe region of the face, use it. 311 maxSafeUVCoord := 1 - faceClipErrorUVCoord 312 if b.Z > 0 { 313 uv = r2.Point{b.X / b.Z, b.Y / b.Z} 314 if math.Max(math.Abs(uv.X), math.Abs(uv.Y)) <= maxSafeUVCoord { 315 return uv, 0 316 } 317 } 318 319 // Otherwise find the point B' where the line AB exits the face. 320 uv = scaledN.exitPoint(scaledN.exitAxis()).Mul(scaleUV) 321 322 p := pointUVW(Point{r3.Vector{uv.X, uv.Y, 1.0}}) 323 324 // Determine if the exit point B' is contained within the segment. We do this 325 // by computing the dot products with two inward-facing tangent vectors at A 326 // and B. If either dot product is negative, we say that B' is on the "wrong 327 // side" of that point. As the point B' moves around the great circle AB past 328 // the segment endpoint B, it is initially on the wrong side of B only; as it 329 // moves further it is on the wrong side of both endpoints; and then it is on 330 // the wrong side of A only. If the exit point B' is on the wrong side of 331 // either endpoint, we can't use it; instead the segment is clipped at the 332 // original endpoint B. 333 // 334 // We reject the segment if the sum of the scores of the two endpoints is 3 335 // or more. Here is what that rule encodes: 336 // - If B' is on the wrong side of A, then the other clipped endpoint A' 337 // must be in the interior of AB (otherwise AB' would go the wrong way 338 // around the circle). There is a similar rule for A'. 339 // - If B' is on the wrong side of either endpoint (and therefore we must 340 // use the original endpoint B instead), then it must be possible to 341 // project B onto this face (i.e., its w-coordinate must be positive). 342 // This rule is only necessary to handle certain zero-length edges (A=B). 343 score := 0 344 if p.Sub(a.Vector).Dot(aTan.Vector) < 0 { 345 score = 2 // B' is on wrong side of A. 346 } else if p.Sub(b.Vector).Dot(bTan.Vector) < 0 { 347 score = 1 // B' is on wrong side of B. 348 } 349 350 if score > 0 { // B' is not in the interior of AB. 351 if b.Z <= 0 { 352 score = 3 // B cannot be projected onto this face. 353 } else { 354 uv = r2.Point{b.X / b.Z, b.Y / b.Z} 355 } 356 } 357 358 return uv, score 359 } 360 361 // updateEndpoint returns the interval with the specified endpoint updated to 362 // the given value. If the value lies beyond the opposite endpoint, nothing is 363 // changed and false is returned. 364 func updateEndpoint(bound r1.Interval, highEndpoint bool, value float64) (r1.Interval, bool) { 365 if !highEndpoint { 366 if bound.Hi < value { 367 return bound, false 368 } 369 if bound.Lo < value { 370 bound.Lo = value 371 } 372 return bound, true 373 } 374 375 if bound.Lo > value { 376 return bound, false 377 } 378 if bound.Hi > value { 379 bound.Hi = value 380 } 381 return bound, true 382 } 383 384 // clipBoundAxis returns the clipped versions of the bounding intervals for the given 385 // axes for the line segment from (a0,a1) to (b0,b1) so that neither extends beyond the 386 // given clip interval. negSlope is a precomputed helper variable that indicates which 387 // diagonal of the bounding box is spanned by AB; it is false if AB has positive slope, 388 // and true if AB has negative slope. If the clipping interval doesn't overlap the bounds, 389 // false is returned. 390 func clipBoundAxis(a0, b0 float64, bound0 r1.Interval, a1, b1 float64, bound1 r1.Interval, 391 negSlope bool, clip r1.Interval) (bound0c, bound1c r1.Interval, updated bool) { 392 393 if bound0.Lo < clip.Lo { 394 // If the upper bound is below the clips lower bound, there is nothing to do. 395 if bound0.Hi < clip.Lo { 396 return bound0, bound1, false 397 } 398 // narrow the intervals lower bound to the clip bound. 399 bound0.Lo = clip.Lo 400 if bound1, updated = updateEndpoint(bound1, negSlope, interpolateFloat64(clip.Lo, a0, b0, a1, b1)); !updated { 401 return bound0, bound1, false 402 } 403 } 404 405 if bound0.Hi > clip.Hi { 406 // If the lower bound is above the clips upper bound, there is nothing to do. 407 if bound0.Lo > clip.Hi { 408 return bound0, bound1, false 409 } 410 // narrow the intervals upper bound to the clip bound. 411 bound0.Hi = clip.Hi 412 if bound1, updated = updateEndpoint(bound1, !negSlope, interpolateFloat64(clip.Hi, a0, b0, a1, b1)); !updated { 413 return bound0, bound1, false 414 } 415 } 416 return bound0, bound1, true 417 } 418 419 // edgeIntersectsRect reports whether the edge defined by AB intersects the 420 // given closed rectangle to within the error bound. 421 func edgeIntersectsRect(a, b r2.Point, r r2.Rect) bool { 422 // First check whether the bounds of a Rect around AB intersects the given rect. 423 if !r.Intersects(r2.RectFromPoints(a, b)) { 424 return false 425 } 426 427 // Otherwise AB intersects the rect if and only if all four vertices of rect 428 // do not lie on the same side of the extended line AB. We test this by finding 429 // the two vertices of rect with minimum and maximum projections onto the normal 430 // of AB, and computing their dot products with the edge normal. 431 n := b.Sub(a).Ortho() 432 433 i := 0 434 if n.X >= 0 { 435 i = 1 436 } 437 j := 0 438 if n.Y >= 0 { 439 j = 1 440 } 441 442 max := n.Dot(r.VertexIJ(i, j).Sub(a)) 443 min := n.Dot(r.VertexIJ(1-i, 1-j).Sub(a)) 444 445 return (max >= 0) && (min <= 0) 446 } 447 448 // clippedEdgeBound returns the bounding rectangle of the portion of the edge defined 449 // by AB intersected by clip. The resulting bound may be empty. This is a convenience 450 // function built on top of clipEdgeBound. 451 func clippedEdgeBound(a, b r2.Point, clip r2.Rect) r2.Rect { 452 bound := r2.RectFromPoints(a, b) 453 if b1, intersects := clipEdgeBound(a, b, clip, bound); intersects { 454 return b1 455 } 456 return r2.EmptyRect() 457 } 458 459 // clipEdgeBound clips an edge AB to sequence of rectangles efficiently. 460 // It represents the clipped edges by their bounding boxes rather than as a pair of 461 // endpoints. Specifically, let A'B' be some portion of an edge AB, and let bound be 462 // a tight bound of A'B'. This function returns the bound that is a tight bound 463 // of A'B' intersected with a given rectangle. If A'B' does not intersect clip, 464 // it returns false and the original bound. 465 func clipEdgeBound(a, b r2.Point, clip, bound r2.Rect) (r2.Rect, bool) { 466 // negSlope indicates which diagonal of the bounding box is spanned by AB: it 467 // is false if AB has positive slope, and true if AB has negative slope. This is 468 // used to determine which interval endpoints need to be updated each time 469 // the edge is clipped. 470 negSlope := (a.X > b.X) != (a.Y > b.Y) 471 472 b0x, b0y, up1 := clipBoundAxis(a.X, b.X, bound.X, a.Y, b.Y, bound.Y, negSlope, clip.X) 473 if !up1 { 474 return bound, false 475 } 476 b1y, b1x, up2 := clipBoundAxis(a.Y, b.Y, b0y, a.X, b.X, b0x, negSlope, clip.Y) 477 if !up2 { 478 return r2.Rect{b0x, b0y}, false 479 } 480 return r2.Rect{X: b1x, Y: b1y}, true 481 } 482 483 // interpolateFloat64 returns a value with the same combination of a1 and b1 as the 484 // given value x is of a and b. This function makes the following guarantees: 485 // - If x == a, then x1 = a1 (exactly). 486 // - If x == b, then x1 = b1 (exactly). 487 // - If a <= x <= b, then a1 <= x1 <= b1 (even if a1 == b1). 488 // 489 // This requires a != b. 490 func interpolateFloat64(x, a, b, a1, b1 float64) float64 { 491 // To get results that are accurate near both A and B, we interpolate 492 // starting from the closer of the two points. 493 if math.Abs(a-x) <= math.Abs(b-x) { 494 return a1 + (b1-a1)*(x-a)/(b-a) 495 } 496 return b1 + (a1-b1)*(x-b)/(a-b) 497 } 498 499 // FaceSegment represents an edge AB clipped to an S2 cube face. It is 500 // represented by a face index and a pair of (u,v) coordinates. 501 type FaceSegment struct { 502 face int 503 a, b r2.Point 504 } 505 506 // FaceSegments subdivides the given edge AB at every point where it crosses the 507 // boundary between two S2 cube faces and returns the corresponding FaceSegments. 508 // The segments are returned in order from A toward B. The input points must be 509 // unit length. 510 // 511 // This function guarantees that the returned segments form a continuous path 512 // from A to B, and that all vertices are within faceClipErrorUVDist of the 513 // line AB. All vertices lie within the [-1,1]x[-1,1] cube face rectangles. 514 // The results are consistent with Sign, i.e. the edge is well-defined even its 515 // endpoints are antipodal. 516 // TODO(roberts): Extend the implementation of PointCross so that this is true. 517 func FaceSegments(a, b Point) []FaceSegment { 518 var segment FaceSegment 519 520 // Fast path: both endpoints are on the same face. 521 var aFace, bFace int 522 aFace, segment.a.X, segment.a.Y = xyzToFaceUV(a.Vector) 523 bFace, segment.b.X, segment.b.Y = xyzToFaceUV(b.Vector) 524 if aFace == bFace { 525 segment.face = aFace 526 return []FaceSegment{segment} 527 } 528 529 // Starting at A, we follow AB from face to face until we reach the face 530 // containing B. The following code is designed to ensure that we always 531 // reach B, even in the presence of numerical errors. 532 // 533 // First we compute the normal to the plane containing A and B. This normal 534 // becomes the ultimate definition of the line AB; it is used to resolve all 535 // questions regarding where exactly the line goes. Unfortunately due to 536 // numerical errors, the line may not quite intersect the faces containing 537 // the original endpoints. We handle this by moving A and/or B slightly if 538 // necessary so that they are on faces intersected by the line AB. 539 ab := a.PointCross(b) 540 541 aFace, segment.a = moveOriginToValidFace(aFace, a, ab, segment.a) 542 bFace, segment.b = moveOriginToValidFace(bFace, b, Point{ab.Mul(-1)}, segment.b) 543 544 // Now we simply follow AB from face to face until we reach B. 545 var segments []FaceSegment 546 segment.face = aFace 547 bSaved := segment.b 548 549 for face := aFace; face != bFace; { 550 // Complete the current segment by finding the point where AB 551 // exits the current face. 552 z := faceXYZtoUVW(face, ab) 553 n := pointUVW{z.Vector} 554 555 exitAxis := n.exitAxis() 556 segment.b = n.exitPoint(exitAxis) 557 segments = append(segments, segment) 558 559 // Compute the next face intersected by AB, and translate the exit 560 // point of the current segment into the (u,v) coordinates of the 561 // next face. This becomes the first point of the next segment. 562 exitXyz := faceUVToXYZ(face, segment.b.X, segment.b.Y) 563 face = nextFace(face, segment.b, exitAxis, n, bFace) 564 exitUvw := faceXYZtoUVW(face, Point{exitXyz}) 565 segment.face = face 566 segment.a = r2.Point{exitUvw.X, exitUvw.Y} 567 } 568 // Finish the last segment. 569 segment.b = bSaved 570 return append(segments, segment) 571 } 572 573 // moveOriginToValidFace updates the origin point to a valid face if necessary. 574 // Given a line segment AB whose origin A has been projected onto a given cube 575 // face, determine whether it is necessary to project A onto a different face 576 // instead. This can happen because the normal of the line AB is not computed 577 // exactly, so that the line AB (defined as the set of points perpendicular to 578 // the normal) may not intersect the cube face containing A. Even if it does 579 // intersect the face, the exit point of the line from that face may be on 580 // the wrong side of A (i.e., in the direction away from B). If this happens, 581 // we reproject A onto the adjacent face where the line AB approaches A most 582 // closely. This moves the origin by a small amount, but never more than the 583 // error tolerances. 584 func moveOriginToValidFace(face int, a, ab Point, aUV r2.Point) (int, r2.Point) { 585 // Fast path: if the origin is sufficiently far inside the face, it is 586 // always safe to use it. 587 const maxSafeUVCoord = 1 - faceClipErrorUVCoord 588 if math.Max(math.Abs((aUV).X), math.Abs((aUV).Y)) <= maxSafeUVCoord { 589 return face, aUV 590 } 591 592 // Otherwise check whether the normal AB even intersects this face. 593 z := faceXYZtoUVW(face, ab) 594 n := pointUVW{z.Vector} 595 if n.intersectsFace() { 596 // Check whether the point where the line AB exits this face is on the 597 // wrong side of A (by more than the acceptable error tolerance). 598 uv := n.exitPoint(n.exitAxis()) 599 exit := faceUVToXYZ(face, uv.X, uv.Y) 600 aTangent := ab.Normalize().Cross(a.Vector) 601 602 // We can use the given face. 603 if exit.Sub(a.Vector).Dot(aTangent) >= -faceClipErrorRadians { 604 return face, aUV 605 } 606 } 607 608 // Otherwise we reproject A to the nearest adjacent face. (If line AB does 609 // not pass through a given face, it must pass through all adjacent faces.) 610 var dir int 611 if math.Abs((aUV).X) >= math.Abs((aUV).Y) { 612 // U-axis 613 if aUV.X > 0 { 614 dir = 1 615 } 616 face = uvwFace(face, 0, dir) 617 } else { 618 // V-axis 619 if aUV.Y > 0 { 620 dir = 1 621 } 622 face = uvwFace(face, 1, dir) 623 } 624 625 aUV.X, aUV.Y = validFaceXYZToUV(face, a.Vector) 626 aUV.X = math.Max(-1.0, math.Min(1.0, aUV.X)) 627 aUV.Y = math.Max(-1.0, math.Min(1.0, aUV.Y)) 628 629 return face, aUV 630 } 631 632 // nextFace returns the next face that should be visited by FaceSegments, given that 633 // we have just visited face and we are following the line AB (represented 634 // by its normal N in the (u,v,w) coordinates of that face). The other 635 // arguments include the point where AB exits face, the corresponding 636 // exit axis, and the target face containing the destination point B. 637 func nextFace(face int, exit r2.Point, axis axis, n pointUVW, targetFace int) int { 638 // this bit is to work around C++ cleverly casting bools to ints for you. 639 exitA := exit.X 640 exit1MinusA := exit.Y 641 642 if axis == axisV { 643 exitA = exit.Y 644 exit1MinusA = exit.X 645 } 646 exitAPos := 0 647 if exitA > 0 { 648 exitAPos = 1 649 } 650 exit1MinusAPos := 0 651 if exit1MinusA > 0 { 652 exit1MinusAPos = 1 653 } 654 655 // We return the face that is adjacent to the exit point along the given 656 // axis. If line AB exits *exactly* through a corner of the face, there are 657 // two possible next faces. If one is the target face containing B, then 658 // we guarantee that we advance to that face directly. 659 // 660 // The three conditions below check that (1) AB exits approximately through 661 // a corner, (2) the adjacent face along the non-exit axis is the target 662 // face, and (3) AB exits *exactly* through the corner. (The sumEqual 663 // code checks whether the dot product of (u,v,1) and n is exactly zero.) 664 if math.Abs(exit1MinusA) == 1 && 665 uvwFace(face, int(1-axis), exit1MinusAPos) == targetFace && 666 sumEqual(exit.X*n.X, exit.Y*n.Y, -n.Z) { 667 return targetFace 668 } 669 670 // Otherwise return the face that is adjacent to the exit point in the 671 // direction of the exit axis. 672 return uvwFace(face, int(axis), exitAPos) 673 } 674