1 // Copyright 2018 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 "github.com/golang/geo/s1" 21 ) 22 23 // PointArea returns the area of triangle ABC. This method combines two different 24 // algorithms to get accurate results for both large and small triangles. 25 // The maximum error is about 5e-15 (about 0.25 square meters on the Earth's 26 // surface), the same as GirardArea below, but unlike that method it is 27 // also accurate for small triangles. Example: when the true area is 100 28 // square meters, PointArea yields an error about 1 trillion times smaller than 29 // GirardArea. 30 // 31 // All points should be unit length, and no two points should be antipodal. 32 // The area is always positive. 33 func PointArea(a, b, c Point) float64 { 34 // This method is based on l'Huilier's theorem, 35 // 36 // tan(E/4) = sqrt(tan(s/2) tan((s-a)/2) tan((s-b)/2) tan((s-c)/2)) 37 // 38 // where E is the spherical excess of the triangle (i.e. its area), 39 // a, b, c are the side lengths, and 40 // s is the semiperimeter (a + b + c) / 2. 41 // 42 // The only significant source of error using l'Huilier's method is the 43 // cancellation error of the terms (s-a), (s-b), (s-c). This leads to a 44 // *relative* error of about 1e-16 * s / min(s-a, s-b, s-c). This compares 45 // to a relative error of about 1e-15 / E using Girard's formula, where E is 46 // the true area of the triangle. Girard's formula can be even worse than 47 // this for very small triangles, e.g. a triangle with a true area of 1e-30 48 // might evaluate to 1e-5. 49 // 50 // So, we prefer l'Huilier's formula unless dmin < s * (0.1 * E), where 51 // dmin = min(s-a, s-b, s-c). This basically includes all triangles 52 // except for extremely long and skinny ones. 53 // 54 // Since we don't know E, we would like a conservative upper bound on 55 // the triangle area in terms of s and dmin. It's possible to show that 56 // E <= k1 * s * sqrt(s * dmin), where k1 = 2*sqrt(3)/Pi (about 1). 57 // Using this, it's easy to show that we should always use l'Huilier's 58 // method if dmin >= k2 * s^5, where k2 is about 1e-2. Furthermore, 59 // if dmin < k2 * s^5, the triangle area is at most k3 * s^4, where 60 // k3 is about 0.1. Since the best case error using Girard's formula 61 // is about 1e-15, this means that we shouldn't even consider it unless 62 // s >= 3e-4 or so. 63 sa := b.stableAngle(c) 64 sb := c.stableAngle(a) 65 sc := a.stableAngle(b) 66 s := 0.5 * (sa + sb + sc) 67 if s >= 3e-4 { 68 // Consider whether Girard's formula might be more accurate. 69 dmin := s - maxAngle(sa, sb, sc) 70 if dmin < 1e-2*s*s*s*s*s { 71 // This triangle is skinny enough to use Girard's formula. 72 area := GirardArea(a, b, c) 73 if dmin < s*0.1*s1.Angle(area+5e-15) { 74 return area 75 } 76 } 77 } 78 79 // Use l'Huilier's formula. 80 return 4 * math.Atan(math.Sqrt(math.Max(0.0, 81 math.Tan(float64(0.5*s))*math.Tan(0.5*float64(s-sa))* 82 math.Tan(0.5*float64(s-sb))*math.Tan(0.5*float64(s-sc))))) 83 } 84 85 // GirardArea returns the area of the triangle computed using Girard's formula. 86 // All points should be unit length, and no two points should be antipodal. 87 // 88 // This method is about twice as fast as PointArea() but has poor relative 89 // accuracy for small triangles. The maximum error is about 5e-15 (about 90 // 0.25 square meters on the Earth's surface) and the average error is about 91 // 1e-15. These bounds apply to triangles of any size, even as the maximum 92 // edge length of the triangle approaches 180 degrees. But note that for 93 // such triangles, tiny perturbations of the input points can change the 94 // true mathematical area dramatically. 95 func GirardArea(a, b, c Point) float64 { 96 // This is equivalent to the usual Girard's formula but is slightly more 97 // accurate, faster to compute, and handles a == b == c without a special 98 // case. PointCross is necessary to get good accuracy when two of 99 // the input points are very close together. 100 ab := a.PointCross(b) 101 bc := b.PointCross(c) 102 ac := a.PointCross(c) 103 104 area := float64(ab.Angle(ac.Vector) - ab.Angle(bc.Vector) + bc.Angle(ac.Vector)) 105 if area < 0 { 106 area = 0 107 } 108 return area 109 } 110 111 // SignedArea returns a positive value for counterclockwise triangles and a negative 112 // value otherwise (similar to PointArea). 113 func SignedArea(a, b, c Point) float64 { 114 return float64(RobustSign(a, b, c)) * PointArea(a, b, c) 115 } 116 117 // Angle returns the interior angle at the vertex B in the triangle ABC. The 118 // return value is always in the range [0, pi]. All points should be 119 // normalized. Ensures that Angle(a,b,c) == Angle(c,b,a) for all a,b,c. 120 // 121 // The angle is undefined if A or C is diametrically opposite from B, and 122 // becomes numerically unstable as the length of edge AB or BC approaches 123 // 180 degrees. 124 func Angle(a, b, c Point) s1.Angle { 125 // PointCross is necessary to get good accuracy when two of the input 126 // points are very close together. 127 return a.PointCross(b).Angle(c.PointCross(b).Vector) 128 } 129 130 // TurnAngle returns the exterior angle at vertex B in the triangle ABC. The 131 // return value is positive if ABC is counterclockwise and negative otherwise. 132 // If you imagine an ant walking from A to B to C, this is the angle that the 133 // ant turns at vertex B (positive = left = CCW, negative = right = CW). 134 // This quantity is also known as the "geodesic curvature" at B. 135 // 136 // Ensures that TurnAngle(a,b,c) == -TurnAngle(c,b,a) for all distinct 137 // a,b,c. The result is undefined if (a == b || b == c), but is either 138 // -Pi or Pi if (a == c). All points should be normalized. 139 func TurnAngle(a, b, c Point) s1.Angle { 140 // We use PointCross to get good accuracy when two points are very 141 // close together, and RobustSign to ensure that the sign is correct for 142 // turns that are close to 180 degrees. 143 angle := a.PointCross(b).Angle(b.PointCross(c).Vector) 144 145 // Don't return RobustSign * angle because it is legal to have (a == c). 146 if RobustSign(a, b, c) == CounterClockwise { 147 return angle 148 } 149 return -angle 150 } 151