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/r3" 21 ) 22 23 // There are several notions of the "centroid" of a triangle. First, there 24 // is the planar centroid, which is simply the centroid of the ordinary 25 // (non-spherical) triangle defined by the three vertices. Second, there is 26 // the surface centroid, which is defined as the intersection of the three 27 // medians of the spherical triangle. It is possible to show that this 28 // point is simply the planar centroid projected to the surface of the 29 // sphere. Finally, there is the true centroid (mass centroid), which is 30 // defined as the surface integral over the spherical triangle of (x,y,z) 31 // divided by the triangle area. This is the point that the triangle would 32 // rotate around if it was spinning in empty space. 33 // 34 // The best centroid for most purposes is the true centroid. Unlike the 35 // planar and surface centroids, the true centroid behaves linearly as 36 // regions are added or subtracted. That is, if you split a triangle into 37 // pieces and compute the average of their centroids (weighted by triangle 38 // area), the result equals the centroid of the original triangle. This is 39 // not true of the other centroids. 40 // 41 // Also note that the surface centroid may be nowhere near the intuitive 42 // "center" of a spherical triangle. For example, consider the triangle 43 // with vertices A=(1,eps,0), B=(0,0,1), C=(-1,eps,0) (a quarter-sphere). 44 // The surface centroid of this triangle is at S=(0, 2*eps, 1), which is 45 // within a distance of 2*eps of the vertex B. Note that the median from A 46 // (the segment connecting A to the midpoint of BC) passes through S, since 47 // this is the shortest path connecting the two endpoints. On the other 48 // hand, the true centroid is at M=(0, 0.5, 0.5), which when projected onto 49 // the surface is a much more reasonable interpretation of the "center" of 50 // this triangle. 51 // 52 53 // TrueCentroid returns the true centroid of the spherical triangle ABC 54 // multiplied by the signed area of spherical triangle ABC. The reasons for 55 // multiplying by the signed area are (1) this is the quantity that needs to be 56 // summed to compute the centroid of a union or difference of triangles, and 57 // (2) it's actually easier to calculate this way. All points must have unit length. 58 // 59 // Note that the result of this function is defined to be Point(0, 0, 0) if 60 // the triangle is degenerate. 61 func TrueCentroid(a, b, c Point) Point { 62 // Use Distance to get accurate results for small triangles. 63 ra := float64(1) 64 if sa := float64(b.Distance(c)); sa != 0 { 65 ra = sa / math.Sin(sa) 66 } 67 rb := float64(1) 68 if sb := float64(c.Distance(a)); sb != 0 { 69 rb = sb / math.Sin(sb) 70 } 71 rc := float64(1) 72 if sc := float64(a.Distance(b)); sc != 0 { 73 rc = sc / math.Sin(sc) 74 } 75 76 // Now compute a point M such that: 77 // 78 // [Ax Ay Az] [Mx] [ra] 79 // [Bx By Bz] [My] = 0.5 * det(A,B,C) * [rb] 80 // [Cx Cy Cz] [Mz] [rc] 81 // 82 // To improve the numerical stability we subtract the first row (A) from the 83 // other two rows; this reduces the cancellation error when A, B, and C are 84 // very close together. Then we solve it using Cramer's rule. 85 // 86 // The result is the true centroid of the triangle multiplied by the 87 // triangle's area. 88 // 89 // This code still isn't as numerically stable as it could be. 90 // The biggest potential improvement is to compute B-A and C-A more 91 // accurately so that (B-A)x(C-A) is always inside triangle ABC. 92 x := r3.Vector{a.X, b.X - a.X, c.X - a.X} 93 y := r3.Vector{a.Y, b.Y - a.Y, c.Y - a.Y} 94 z := r3.Vector{a.Z, b.Z - a.Z, c.Z - a.Z} 95 r := r3.Vector{ra, rb - ra, rc - ra} 96 97 return Point{r3.Vector{y.Cross(z).Dot(r), z.Cross(x).Dot(r), x.Cross(y).Dot(r)}.Mul(0.5)} 98 } 99 100 // EdgeTrueCentroid returns the true centroid of the spherical geodesic edge AB 101 // multiplied by the length of the edge AB. As with triangles, the true centroid 102 // of a collection of line segments may be computed simply by summing the result 103 // of this method for each segment. 104 // 105 // Note that the planar centroid of a line segment is simply 0.5 * (a + b), 106 // while the surface centroid is (a + b).Normalize(). However neither of 107 // these values is appropriate for computing the centroid of a collection of 108 // edges (such as a polyline). 109 // 110 // Also note that the result of this function is defined to be Point(0, 0, 0) 111 // if the edge is degenerate. 112 func EdgeTrueCentroid(a, b Point) Point { 113 // The centroid (multiplied by length) is a vector toward the midpoint 114 // of the edge, whose length is twice the sine of half the angle between 115 // the two vertices. Defining theta to be this angle, we have: 116 vDiff := a.Sub(b.Vector) // Length == 2*sin(theta) 117 vSum := a.Add(b.Vector) // Length == 2*cos(theta) 118 sin2 := vDiff.Norm2() 119 cos2 := vSum.Norm2() 120 if cos2 == 0 { 121 return Point{} // Ignore antipodal edges. 122 } 123 return Point{vSum.Mul(math.Sqrt(sin2 / cos2))} // Length == 2*sin(theta) 124 } 125 126 // PlanarCentroid returns the centroid of the planar triangle ABC. This can be 127 // normalized to unit length to obtain the "surface centroid" of the corresponding 128 // spherical triangle, i.e. the intersection of the three medians. However, note 129 // that for large spherical triangles the surface centroid may be nowhere near 130 // the intuitive "center". 131 func PlanarCentroid(a, b, c Point) Point { 132 return Point{a.Add(b.Vector).Add(c.Vector).Mul(1. / 3)} 133 } 134