gtsocial-umbx

centroids.go (5964B)

---

 // Copyright 2018 Google Inc. All rights reserved.
 //
 // Licensed under the Apache License, Version 2.0 (the "License");
 // you may not use this file except in compliance with the License.
 // You may obtain a copy of the License at
 //
 //     http://www.apache.org/licenses/LICENSE-2.0
 //
 // Unless required by applicable law or agreed to in writing, software
 // distributed under the License is distributed on an "AS IS" BASIS,
 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 // See the License for the specific language governing permissions and
 // limitations under the License.
 
 package s2
 
 import (
 	"math"
 
 	"github.com/golang/geo/r3"
 )
 
 // There are several notions of the "centroid" of a triangle. First, there
 // is the planar centroid, which is simply the centroid of the ordinary
 // (non-spherical) triangle defined by the three vertices. Second, there is
 // the surface centroid, which is defined as the intersection of the three
 // medians of the spherical triangle. It is possible to show that this
 // point is simply the planar centroid projected to the surface of the
 // sphere. Finally, there is the true centroid (mass centroid), which is
 // defined as the surface integral over the spherical triangle of (x,y,z)
 // divided by the triangle area. This is the point that the triangle would
 // rotate around if it was spinning in empty space.
 //
 // The best centroid for most purposes is the true centroid. Unlike the
 // planar and surface centroids, the true centroid behaves linearly as
 // regions are added or subtracted. That is, if you split a triangle into
 // pieces and compute the average of their centroids (weighted by triangle
 // area), the result equals the centroid of the original triangle. This is
 // not true of the other centroids.
 //
 // Also note that the surface centroid may be nowhere near the intuitive
 // "center" of a spherical triangle. For example, consider the triangle
 // with vertices A=(1,eps,0), B=(0,0,1), C=(-1,eps,0) (a quarter-sphere).
 // The surface centroid of this triangle is at S=(0, 2*eps, 1), which is
 // within a distance of 2*eps of the vertex B. Note that the median from A
 // (the segment connecting A to the midpoint of BC) passes through S, since
 // this is the shortest path connecting the two endpoints. On the other
 // hand, the true centroid is at M=(0, 0.5, 0.5), which when projected onto
 // the surface is a much more reasonable interpretation of the "center" of
 // this triangle.
 //
 
 // TrueCentroid returns the true centroid of the spherical triangle ABC
 // multiplied by the signed area of spherical triangle ABC. The reasons for
 // multiplying by the signed area are (1) this is the quantity that needs to be
 // summed to compute the centroid of a union or difference of triangles, and
 // (2) it's actually easier to calculate this way. All points must have unit length.
 //
 // Note that the result of this function is defined to be Point(0, 0, 0) if
 // the triangle is degenerate.
 func TrueCentroid(a, b, c Point) Point {
 	// Use Distance to get accurate results for small triangles.
 	ra := float64(1)
 	if sa := float64(b.Distance(c)); sa != 0 {
 		ra = sa / math.Sin(sa)
 	}
 	rb := float64(1)
 	if sb := float64(c.Distance(a)); sb != 0 {
 		rb = sb / math.Sin(sb)
 	}
 	rc := float64(1)
 	if sc := float64(a.Distance(b)); sc != 0 {
 		rc = sc / math.Sin(sc)
 	}
 
 	// Now compute a point M such that:
 	//
 	//  [Ax Ay Az] [Mx]                       [ra]
 	//  [Bx By Bz] [My]  = 0.5 * det(A,B,C) * [rb]
 	//  [Cx Cy Cz] [Mz]                       [rc]
 	//
 	// To improve the numerical stability we subtract the first row (A) from the
 	// other two rows; this reduces the cancellation error when A, B, and C are
 	// very close together. Then we solve it using Cramer's rule.
 	//
 	// The result is the true centroid of the triangle multiplied by the
 	// triangle's area.
 	//
 	// This code still isn't as numerically stable as it could be.
 	// The biggest potential improvement is to compute B-A and C-A more
 	// accurately so that (B-A)x(C-A) is always inside triangle ABC.
 	x := r3.Vector{a.X, b.X - a.X, c.X - a.X}
 	y := r3.Vector{a.Y, b.Y - a.Y, c.Y - a.Y}
 	z := r3.Vector{a.Z, b.Z - a.Z, c.Z - a.Z}
 	r := r3.Vector{ra, rb - ra, rc - ra}
 
 	return Point{r3.Vector{y.Cross(z).Dot(r), z.Cross(x).Dot(r), x.Cross(y).Dot(r)}.Mul(0.5)}
 }
 
 // EdgeTrueCentroid returns the true centroid of the spherical geodesic edge AB
 // multiplied by the length of the edge AB. As with triangles, the true centroid
 // of a collection of line segments may be computed simply by summing the result
 // of this method for each segment.
 //
 // Note that the planar centroid of a line segment is simply 0.5 * (a + b),
 // while the surface centroid is (a + b).Normalize(). However neither of
 // these values is appropriate for computing the centroid of a collection of
 // edges (such as a polyline).
 //
 // Also note that the result of this function is defined to be Point(0, 0, 0)
 // if the edge is degenerate.
 func EdgeTrueCentroid(a, b Point) Point {
 	// The centroid (multiplied by length) is a vector toward the midpoint
 	// of the edge, whose length is twice the sine of half the angle between
 	// the two vertices. Defining theta to be this angle, we have:
 	vDiff := a.Sub(b.Vector) // Length == 2*sin(theta)
 	vSum := a.Add(b.Vector)  // Length == 2*cos(theta)
 	sin2 := vDiff.Norm2()
 	cos2 := vSum.Norm2()
 	if cos2 == 0 {
 		return Point{} // Ignore antipodal edges.
 	}
 	return Point{vSum.Mul(math.Sqrt(sin2 / cos2))} // Length == 2*sin(theta)
 }
 
 // PlanarCentroid returns the centroid of the planar triangle ABC. This can be
 // normalized to unit length to obtain the "surface centroid" of the corresponding
 // spherical triangle, i.e. the intersection of the three medians. However, note
 // that for large spherical triangles the surface centroid may be nowhere near
 // the intuitive "center".
 func PlanarCentroid(a, b, c Point) Point {
 	return Point{a.Add(b.Vector).Add(c.Vector).Mul(1. / 3)}
 }