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 import ( 18 "math" 19 ) 20 21 // EdgeCrosser allows edges to be efficiently tested for intersection with a 22 // given fixed edge AB. It is especially efficient when testing for 23 // intersection with an edge chain connecting vertices v0, v1, v2, ... 24 // 25 // Example usage: 26 // 27 // func CountIntersections(a, b Point, edges []Edge) int { 28 // count := 0 29 // crosser := NewEdgeCrosser(a, b) 30 // for _, edge := range edges { 31 // if crosser.CrossingSign(&edge.First, &edge.Second) != DoNotCross { 32 // count++ 33 // } 34 // } 35 // return count 36 // } 37 type EdgeCrosser struct { 38 a Point 39 b Point 40 aXb Point 41 42 // To reduce the number of calls to expensiveSign, we compute an 43 // outward-facing tangent at A and B if necessary. If the plane 44 // perpendicular to one of these tangents separates AB from CD (i.e., one 45 // edge on each side) then there is no intersection. 46 aTangent Point // Outward-facing tangent at A. 47 bTangent Point // Outward-facing tangent at B. 48 49 // The fields below are updated for each vertex in the chain. 50 c Point // Previous vertex in the vertex chain. 51 acb Direction // The orientation of triangle ACB. 52 } 53 54 // NewEdgeCrosser returns an EdgeCrosser with the fixed edge AB. 55 func NewEdgeCrosser(a, b Point) *EdgeCrosser { 56 norm := a.PointCross(b) 57 return &EdgeCrosser{ 58 a: a, 59 b: b, 60 aXb: Point{a.Cross(b.Vector)}, 61 aTangent: Point{a.Cross(norm.Vector)}, 62 bTangent: Point{norm.Cross(b.Vector)}, 63 } 64 } 65 66 // CrossingSign reports whether the edge AB intersects the edge CD. If any two 67 // vertices from different edges are the same, returns MaybeCross. If either edge 68 // is degenerate (A == B or C == D), returns either DoNotCross or MaybeCross. 69 // 70 // Properties of CrossingSign: 71 // 72 // (1) CrossingSign(b,a,c,d) == CrossingSign(a,b,c,d) 73 // (2) CrossingSign(c,d,a,b) == CrossingSign(a,b,c,d) 74 // (3) CrossingSign(a,b,c,d) == MaybeCross if a==c, a==d, b==c, b==d 75 // (3) CrossingSign(a,b,c,d) == DoNotCross or MaybeCross if a==b or c==d 76 // 77 // Note that if you want to check an edge against a chain of other edges, 78 // it is slightly more efficient to use the single-argument version 79 // ChainCrossingSign below. 80 func (e *EdgeCrosser) CrossingSign(c, d Point) Crossing { 81 if c != e.c { 82 e.RestartAt(c) 83 } 84 return e.ChainCrossingSign(d) 85 } 86 87 // EdgeOrVertexCrossing reports whether if CrossingSign(c, d) > 0, or AB and 88 // CD share a vertex and VertexCrossing(a, b, c, d) is true. 89 // 90 // This method extends the concept of a "crossing" to the case where AB 91 // and CD have a vertex in common. The two edges may or may not cross, 92 // according to the rules defined in VertexCrossing above. The rules 93 // are designed so that point containment tests can be implemented simply 94 // by counting edge crossings. Similarly, determining whether one edge 95 // chain crosses another edge chain can be implemented by counting. 96 func (e *EdgeCrosser) EdgeOrVertexCrossing(c, d Point) bool { 97 if c != e.c { 98 e.RestartAt(c) 99 } 100 return e.EdgeOrVertexChainCrossing(d) 101 } 102 103 // NewChainEdgeCrosser is a convenience constructor that uses AB as the fixed edge, 104 // and C as the first vertex of the vertex chain (equivalent to calling RestartAt(c)). 105 // 106 // You don't need to use this or any of the chain functions unless you're trying to 107 // squeeze out every last drop of performance. Essentially all you are saving is a test 108 // whether the first vertex of the current edge is the same as the second vertex of the 109 // previous edge. 110 func NewChainEdgeCrosser(a, b, c Point) *EdgeCrosser { 111 e := NewEdgeCrosser(a, b) 112 e.RestartAt(c) 113 return e 114 } 115 116 // RestartAt sets the current point of the edge crosser to be c. 117 // Call this method when your chain 'jumps' to a new place. 118 // The argument must point to a value that persists until the next call. 119 func (e *EdgeCrosser) RestartAt(c Point) { 120 e.c = c 121 e.acb = -triageSign(e.a, e.b, e.c) 122 } 123 124 // ChainCrossingSign is like CrossingSign, but uses the last vertex passed to one of 125 // the crossing methods (or RestartAt) as the first vertex of the current edge. 126 func (e *EdgeCrosser) ChainCrossingSign(d Point) Crossing { 127 // For there to be an edge crossing, the triangles ACB, CBD, BDA, DAC must 128 // all be oriented the same way (CW or CCW). We keep the orientation of ACB 129 // as part of our state. When each new point D arrives, we compute the 130 // orientation of BDA and check whether it matches ACB. This checks whether 131 // the points C and D are on opposite sides of the great circle through AB. 132 133 // Recall that triageSign is invariant with respect to rotating its 134 // arguments, i.e. ABD has the same orientation as BDA. 135 bda := triageSign(e.a, e.b, d) 136 if e.acb == -bda && bda != Indeterminate { 137 // The most common case -- triangles have opposite orientations. Save the 138 // current vertex D as the next vertex C, and also save the orientation of 139 // the new triangle ACB (which is opposite to the current triangle BDA). 140 e.c = d 141 e.acb = -bda 142 return DoNotCross 143 } 144 return e.crossingSign(d, bda) 145 } 146 147 // EdgeOrVertexChainCrossing is like EdgeOrVertexCrossing, but uses the last vertex 148 // passed to one of the crossing methods (or RestartAt) as the first vertex of the current edge. 149 func (e *EdgeCrosser) EdgeOrVertexChainCrossing(d Point) bool { 150 // We need to copy e.c since it is clobbered by ChainCrossingSign. 151 c := e.c 152 switch e.ChainCrossingSign(d) { 153 case DoNotCross: 154 return false 155 case Cross: 156 return true 157 } 158 return VertexCrossing(e.a, e.b, c, d) 159 } 160 161 // crossingSign handle the slow path of CrossingSign. 162 func (e *EdgeCrosser) crossingSign(d Point, bda Direction) Crossing { 163 // Compute the actual result, and then save the current vertex D as the next 164 // vertex C, and save the orientation of the next triangle ACB (which is 165 // opposite to the current triangle BDA). 166 defer func() { 167 e.c = d 168 e.acb = -bda 169 }() 170 171 // At this point, a very common situation is that A,B,C,D are four points on 172 // a line such that AB does not overlap CD. (For example, this happens when 173 // a line or curve is sampled finely, or when geometry is constructed by 174 // computing the union of S2CellIds.) Most of the time, we can determine 175 // that AB and CD do not intersect using the two outward-facing 176 // tangents at A and B (parallel to AB) and testing whether AB and CD are on 177 // opposite sides of the plane perpendicular to one of these tangents. This 178 // is moderately expensive but still much cheaper than expensiveSign. 179 180 // The error in RobustCrossProd is insignificant. The maximum error in 181 // the call to CrossProd (i.e., the maximum norm of the error vector) is 182 // (0.5 + 1/sqrt(3)) * dblEpsilon. The maximum error in each call to 183 // DotProd below is dblEpsilon. (There is also a small relative error 184 // term that is insignificant because we are comparing the result against a 185 // constant that is very close to zero.) 186 maxError := (1.5 + 1/math.Sqrt(3)) * dblEpsilon 187 if (e.c.Dot(e.aTangent.Vector) > maxError && d.Dot(e.aTangent.Vector) > maxError) || (e.c.Dot(e.bTangent.Vector) > maxError && d.Dot(e.bTangent.Vector) > maxError) { 188 return DoNotCross 189 } 190 191 // Otherwise, eliminate the cases where two vertices from different edges are 192 // equal. (These cases could be handled in the code below, but we would rather 193 // avoid calling ExpensiveSign if possible.) 194 if e.a == e.c || e.a == d || e.b == e.c || e.b == d { 195 return MaybeCross 196 } 197 198 // Eliminate the cases where an input edge is degenerate. (Note that in 199 // most cases, if CD is degenerate then this method is not even called 200 // because acb and bda have different signs.) 201 if e.a == e.b || e.c == d { 202 return DoNotCross 203 } 204 205 // Otherwise it's time to break out the big guns. 206 if e.acb == Indeterminate { 207 e.acb = -expensiveSign(e.a, e.b, e.c) 208 } 209 if bda == Indeterminate { 210 bda = expensiveSign(e.a, e.b, d) 211 } 212 213 if bda != e.acb { 214 return DoNotCross 215 } 216 217 cbd := -RobustSign(e.c, d, e.b) 218 if cbd != e.acb { 219 return DoNotCross 220 } 221 dac := RobustSign(e.c, d, e.a) 222 if dac != e.acb { 223 return DoNotCross 224 } 225 return Cross 226 } 227