gotosocial/vendor/github.com/golang/geo/s2/centroids.go
Tobi Smethurst 98263a7de6
Grand test fixup (#138)
* start fixing up tests

* fix up tests + automate with drone

* fiddle with linting

* messing about with drone.yml

* some more fiddling

* hmmm

* add cache

* add vendor directory

* verbose

* ci updates

* update some little things

* update sig
2021-08-12 21:03:24 +02:00

134 lines
5.8 KiB
Go

// 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)}
}