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 "sort" 22 23 "github.com/golang/geo/r3" 24 "github.com/golang/geo/s1" 25 ) 26 27 // Point represents a point on the unit sphere as a normalized 3D vector. 28 // Fields should be treated as read-only. Use one of the factory methods for creation. 29 type Point struct { 30 r3.Vector 31 } 32 33 // sortPoints sorts the slice of Points in place. 34 func sortPoints(e []Point) { 35 sort.Sort(points(e)) 36 } 37 38 // points implements the Sort interface for slices of Point. 39 type points []Point 40 41 func (p points) Len() int { return len(p) } 42 func (p points) Swap(i, j int) { p[i], p[j] = p[j], p[i] } 43 func (p points) Less(i, j int) bool { return p[i].Cmp(p[j].Vector) == -1 } 44 45 // PointFromCoords creates a new normalized point from coordinates. 46 // 47 // This always returns a valid point. If the given coordinates can not be normalized 48 // the origin point will be returned. 49 // 50 // This behavior is different from the C++ construction of a S2Point from coordinates 51 // (i.e. S2Point(x, y, z)) in that in C++ they do not Normalize. 52 func PointFromCoords(x, y, z float64) Point { 53 if x == 0 && y == 0 && z == 0 { 54 return OriginPoint() 55 } 56 return Point{r3.Vector{x, y, z}.Normalize()} 57 } 58 59 // OriginPoint returns a unique "origin" on the sphere for operations that need a fixed 60 // reference point. In particular, this is the "point at infinity" used for 61 // point-in-polygon testing (by counting the number of edge crossings). 62 // 63 // It should *not* be a point that is commonly used in edge tests in order 64 // to avoid triggering code to handle degenerate cases (this rules out the 65 // north and south poles). It should also not be on the boundary of any 66 // low-level S2Cell for the same reason. 67 func OriginPoint() Point { 68 return Point{r3.Vector{-0.0099994664350250197, 0.0025924542609324121, 0.99994664350250195}} 69 } 70 71 // PointCross returns a Point that is orthogonal to both p and op. This is similar to 72 // p.Cross(op) (the true cross product) except that it does a better job of 73 // ensuring orthogonality when the Point is nearly parallel to op, it returns 74 // a non-zero result even when p == op or p == -op and the result is a Point. 75 // 76 // It satisfies the following properties (f == PointCross): 77 // 78 // (1) f(p, op) != 0 for all p, op 79 // (2) f(op,p) == -f(p,op) unless p == op or p == -op 80 // (3) f(-p,op) == -f(p,op) unless p == op or p == -op 81 // (4) f(p,-op) == -f(p,op) unless p == op or p == -op 82 func (p Point) PointCross(op Point) Point { 83 // NOTE(dnadasi): In the C++ API the equivalent method here was known as "RobustCrossProd", 84 // but PointCross more accurately describes how this method is used. 85 x := p.Add(op.Vector).Cross(op.Sub(p.Vector)) 86 87 // Compare exactly to the 0 vector. 88 if x == (r3.Vector{}) { 89 // The only result that makes sense mathematically is to return zero, but 90 // we find it more convenient to return an arbitrary orthogonal vector. 91 return Point{p.Ortho()} 92 } 93 94 return Point{x} 95 } 96 97 // OrderedCCW returns true if the edges OA, OB, and OC are encountered in that 98 // order while sweeping CCW around the point O. 99 // 100 // You can think of this as testing whether A <= B <= C with respect to the 101 // CCW ordering around O that starts at A, or equivalently, whether B is 102 // contained in the range of angles (inclusive) that starts at A and extends 103 // CCW to C. Properties: 104 // 105 // (1) If OrderedCCW(a,b,c,o) && OrderedCCW(b,a,c,o), then a == b 106 // (2) If OrderedCCW(a,b,c,o) && OrderedCCW(a,c,b,o), then b == c 107 // (3) If OrderedCCW(a,b,c,o) && OrderedCCW(c,b,a,o), then a == b == c 108 // (4) If a == b or b == c, then OrderedCCW(a,b,c,o) is true 109 // (5) Otherwise if a == c, then OrderedCCW(a,b,c,o) is false 110 func OrderedCCW(a, b, c, o Point) bool { 111 sum := 0 112 if RobustSign(b, o, a) != Clockwise { 113 sum++ 114 } 115 if RobustSign(c, o, b) != Clockwise { 116 sum++ 117 } 118 if RobustSign(a, o, c) == CounterClockwise { 119 sum++ 120 } 121 return sum >= 2 122 } 123 124 // Distance returns the angle between two points. 125 func (p Point) Distance(b Point) s1.Angle { 126 return p.Vector.Angle(b.Vector) 127 } 128 129 // ApproxEqual reports whether the two points are similar enough to be equal. 130 func (p Point) ApproxEqual(other Point) bool { 131 return p.approxEqual(other, s1.Angle(epsilon)) 132 } 133 134 // approxEqual reports whether the two points are within the given epsilon. 135 func (p Point) approxEqual(other Point, eps s1.Angle) bool { 136 return p.Vector.Angle(other.Vector) <= eps 137 } 138 139 // ChordAngleBetweenPoints constructs a ChordAngle corresponding to the distance 140 // between the two given points. The points must be unit length. 141 func ChordAngleBetweenPoints(x, y Point) s1.ChordAngle { 142 return s1.ChordAngle(math.Min(4.0, x.Sub(y.Vector).Norm2())) 143 } 144 145 // regularPoints generates a slice of points shaped as a regular polygon with 146 // the numVertices vertices, all located on a circle of the specified angular radius 147 // around the center. The radius is the actual distance from center to each vertex. 148 func regularPoints(center Point, radius s1.Angle, numVertices int) []Point { 149 return regularPointsForFrame(getFrame(center), radius, numVertices) 150 } 151 152 // regularPointsForFrame generates a slice of points shaped as a regular polygon 153 // with numVertices vertices, all on a circle of the specified angular radius around 154 // the center. The radius is the actual distance from the center to each vertex. 155 func regularPointsForFrame(frame matrix3x3, radius s1.Angle, numVertices int) []Point { 156 // We construct the loop in the given frame coordinates, with the center at 157 // (0, 0, 1). For a loop of radius r, the loop vertices have the form 158 // (x, y, z) where x^2 + y^2 = sin(r) and z = cos(r). The distance on the 159 // sphere (arc length) from each vertex to the center is acos(cos(r)) = r. 160 z := math.Cos(radius.Radians()) 161 r := math.Sin(radius.Radians()) 162 radianStep := 2 * math.Pi / float64(numVertices) 163 var vertices []Point 164 165 for i := 0; i < numVertices; i++ { 166 angle := float64(i) * radianStep 167 p := Point{r3.Vector{r * math.Cos(angle), r * math.Sin(angle), z}} 168 vertices = append(vertices, Point{fromFrame(frame, p).Normalize()}) 169 } 170 171 return vertices 172 } 173 174 // CapBound returns a bounding cap for this point. 175 func (p Point) CapBound() Cap { 176 return CapFromPoint(p) 177 } 178 179 // RectBound returns a bounding latitude-longitude rectangle from this point. 180 func (p Point) RectBound() Rect { 181 return RectFromLatLng(LatLngFromPoint(p)) 182 } 183 184 // ContainsCell returns false as Points do not contain any other S2 types. 185 func (p Point) ContainsCell(c Cell) bool { return false } 186 187 // IntersectsCell reports whether this Point intersects the given cell. 188 func (p Point) IntersectsCell(c Cell) bool { 189 return c.ContainsPoint(p) 190 } 191 192 // ContainsPoint reports if this Point contains the other Point. 193 // (This method is named to satisfy the Region interface.) 194 func (p Point) ContainsPoint(other Point) bool { 195 return p.Contains(other) 196 } 197 198 // CellUnionBound computes a covering of the Point. 199 func (p Point) CellUnionBound() []CellID { 200 return p.CapBound().CellUnionBound() 201 } 202 203 // Contains reports if this Point contains the other Point. 204 // (This method matches all other s2 types where the reflexive Contains 205 // method does not contain the type's name.) 206 func (p Point) Contains(other Point) bool { return p == other } 207 208 // Encode encodes the Point. 209 func (p Point) Encode(w io.Writer) error { 210 e := &encoder{w: w} 211 p.encode(e) 212 return e.err 213 } 214 215 func (p Point) encode(e *encoder) { 216 e.writeInt8(encodingVersion) 217 e.writeFloat64(p.X) 218 e.writeFloat64(p.Y) 219 e.writeFloat64(p.Z) 220 } 221 222 // Decode decodes the Point. 223 func (p *Point) Decode(r io.Reader) error { 224 d := &decoder{r: asByteReader(r)} 225 p.decode(d) 226 return d.err 227 } 228 229 func (p *Point) decode(d *decoder) { 230 version := d.readInt8() 231 if d.err != nil { 232 return 233 } 234 if version != encodingVersion { 235 d.err = fmt.Errorf("only version %d is supported", encodingVersion) 236 return 237 } 238 p.X = d.readFloat64() 239 p.Y = d.readFloat64() 240 p.Z = d.readFloat64() 241 } 242 243 // Ortho returns a unit-length vector that is orthogonal to "a". Satisfies 244 // Ortho(-a) = -Ortho(a) for all a. 245 // 246 // Note that Vector3 also defines an "Ortho" method, but this one is 247 // preferred for use in S2 code because it explicitly tries to avoid result 248 // coordinates that are zero. (This is a performance optimization that 249 // reduces the amount of time spent in functions that handle degeneracies.) 250 func Ortho(a Point) Point { 251 temp := r3.Vector{0.012, 0.0053, 0.00457} 252 switch a.LargestComponent() { 253 case r3.XAxis: 254 temp.Z = 1 255 case r3.YAxis: 256 temp.X = 1 257 default: 258 temp.Y = 1 259 } 260 return Point{a.Cross(temp).Normalize()} 261 } 262 263 // referenceDir returns a unit-length vector to use as the reference direction for 264 // deciding whether a polygon with semi-open boundaries contains the given vertex "a" 265 // (see ContainsVertexQuery). The result is unit length and is guaranteed 266 // to be different from the given point "a". 267 func (p Point) referenceDir() Point { 268 return Ortho(p) 269 } 270 271 // Rotate the given point about the given axis by the given angle. p and 272 // axis must be unit length; angle has no restrictions (e.g., it can be 273 // positive, negative, greater than 360 degrees, etc). 274 func Rotate(p, axis Point, angle s1.Angle) Point { 275 // Let M be the plane through P that is perpendicular to axis, and let 276 // center be the point where M intersects axis. We construct a 277 // right-handed orthogonal frame (dx, dy, center) such that dx is the 278 // vector from center to P, and dy has the same length as dx. The 279 // result can then be expressed as (cos(angle)*dx + sin(angle)*dy + center). 280 center := axis.Mul(p.Dot(axis.Vector)) 281 dx := p.Sub(center) 282 dy := axis.Cross(p.Vector) 283 // Mathematically the result is unit length, but normalization is necessary 284 // to ensure that numerical errors don't accumulate. 285 return Point{dx.Mul(math.Cos(angle.Radians())).Add(dy.Mul(math.Sin(angle.Radians()))).Add(center).Normalize()} 286 } 287 288 // stableAngle reports the angle between two vectors with better precision when 289 // the two are nearly parallel. 290 // 291 // The .Angle() member function uses atan(|AxB|, A.B) to compute the angle 292 // between A and B, which can lose about half its precision when A and B are 293 // nearly (anti-)parallel. 294 // 295 // Kahan provides a much more stable form: 296 // 297 // 2*atan2(| A*|B| - |A|*B |, | A*|B| + |A|*B |) 298 // 299 // Since Points are unit magnitude by construction we can simplify further: 300 // 301 // 2*atan2(|A-B|,|A+B|) 302 // 303 // This likely can't replace Vectors Angle since it requires four magnitude 304 // calculations, each of which takes 5 operations + a square root, plus 6 305 // operations to find the sum and difference of the vectors, for a total of 26 + 306 // 4 square roots. Vectors Angle requires 19 + 1 square root. 307 // 308 // Since we always have unit vectors, we can elide two of those magnitude 309 // calculations for a total of 16 + 2 square roots which is competitive with 310 // Vectors Angle performance. 311 // 312 // Reference: Kahan, W. (2006, Jan 11). "How Futile are Mindless Assessments of 313 // Roundoff in Floating-Point Computation?" (p. 47). 314 // https://people.eecs.berkeley.edu/~wkahan/Mindless.pdf 315 // 316 // The 2 points must be normalized. 317 func (p Point) stableAngle(o Point) s1.Angle { 318 return s1.Angle(2 * math.Atan2(p.Sub(o.Vector).Norm(), p.Add(o.Vector).Norm())) 319 } 320