spatial/r3: add Triangle

gonum · Jun 25, 2022 · 67f3e1d · 67f3e1d
1 parent 7a3f4f5
commit 67f3e1d
Show file tree

Hide file tree

Showing 3 changed files with 424 additions and 0 deletions.
diff --git a/spatial/r3/icosahedron_example_test.go b/spatial/r3/icosahedron_example_test.go
@@ -0,0 +1,104 @@
+// Copyright ©2022 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package r3_test
+
+import (
+	"fmt"
+
+	"gonum.org/v1/gonum/spatial/r3"
+)
+
+func ExampleTriangle_icosphere() {
+	// This example generates a 3D icosphere from
+	// a starting icosahedron by subdividing surfaces.
+	// See https://schneide.blog/2016/07/15/generating-an-icosphere-in-c/.
+	const subdivisions = 5
+	// vertices is a slice of r3.Vec
+	// triangles is a slice of [3]int indices
+	// referencing the vertices.
+	vertices, triangles := icosahedron()
+	for i := 0; i < subdivisions; i++ {
+		vertices, triangles = subdivide(vertices, triangles)
+	}
+	var faces []r3.Triangle
+	for _, t := range triangles {
+		var face r3.Triangle
+		for i := 0; i < 3; i++ {
+			face[i] = vertices[t[i]]
+		}
+		faces = append(faces, face)
+	}
+	fmt.Println(faces)
+	// The 3D rendering of the icosphere is left as an exercise to the reader.
+}
+
+// edgeIdx represents an edge of the icosahedron
+type edgeIdx [2]int
+
+func subdivide(vertices []r3.Vec, triangles [][3]int) ([]r3.Vec, [][3]int) {
+	// We generate a lookup table of all newly generated vertices so as to not
+	// duplicate new vertices. edgeIdx has lower index first.
+	lookup := make(map[edgeIdx]int)
+	var result [][3]int
+	for _, triangle := range triangles {
+		var mid [3]int
+		for edge := 0; edge < 3; edge++ {
+			lookup, mid[edge], vertices = subdivideEdge(lookup, vertices, triangle[edge], triangle[(edge+1)%3])
+		}
+		newTriangles := [][3]int{
+			{triangle[0], mid[0], mid[2]},
+			{triangle[1], mid[1], mid[0]},
+			{triangle[2], mid[2], mid[1]},
+			{mid[0], mid[1], mid[2]},
+		}
+		result = append(result, newTriangles...)
+	}
+	return vertices, result
+}
+
+// subdivideEdge takes the vertices list and indices first and second which
+// refer to the edge that will be subdivided.
+// lookup is a table of all newly generated vertices from
+// previous calls to subdivideEdge so as to not duplicate vertices.
+func subdivideEdge(lookup map[edgeIdx]int, vertices []r3.Vec, first, second int) (map[edgeIdx]int, int, []r3.Vec) {
+	key := edgeIdx{first, second}
+	if first > second {
+		// Swap to ensure edgeIdx always has lower index first.
+		key[0], key[1] = key[1], key[0]
+	}
+	vertIdx, vertExists := lookup[key]
+	if !vertExists {
+		// If edge not already subdivided add
+		// new dividing vertex to lookup table.
+		edge0 := vertices[first]
+		edge1 := vertices[second]
+		point := r3.Unit(r3.Add(edge0, edge1)) // vertex at a normalized position.
+		vertices = append(vertices, point)
+		vertIdx = len(vertices) - 1
+		lookup[key] = vertIdx
+	}
+	return lookup, vertIdx, vertices
+}
+
+// icosahedron returns an icosahedron mesh.
+func icosahedron() (vertices []r3.Vec, triangles [][3]int) {
+	const (
+		radiusSqrt = 1.0 // Example designed for unit sphere generation.
+		X          = radiusSqrt * .525731112119133606
+		Z          = radiusSqrt * .850650808352039932
+		N          = 0.0
+	)
+	return []r3.Vec{
+			{-X, N, Z}, {X, N, Z}, {-X, N, -Z}, {X, N, -Z},
+			{N, Z, X}, {N, Z, -X}, {N, -Z, X}, {N, -Z, -X},
+			{Z, X, N}, {-Z, X, N}, {Z, -X, N}, {-Z, -X, N},
+		}, [][3]int{
+			{0, 1, 4}, {0, 4, 9}, {9, 4, 5}, {4, 8, 5},
+			{4, 1, 8}, {8, 1, 10}, {8, 10, 3}, {5, 8, 3},
+			{5, 3, 2}, {2, 3, 7}, {7, 3, 10}, {7, 10, 6},
+			{7, 6, 11}, {11, 6, 0}, {0, 6, 1}, {6, 10, 1},
+			{9, 11, 0}, {9, 2, 11}, {9, 5, 2}, {7, 11, 2},
+		}
+}
diff --git a/spatial/r3/triangle.go b/spatial/r3/triangle.go
@@ -0,0 +1,117 @@
+// Copyright ©2022 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package r3
+
+import "math"
+
+// Triangle represents a triangle in 3D space and
+// is composed by 3 vectors corresponding to the position
+// of each of the vertices. Ordering of these vertices
+// decides the "normal" direction.
+// Inverting ordering of two vertices inverts the resulting direction.
+type Triangle [3]Vec
+
+// Centroid returns the intersection of the three medians of the triangle
+// as a point in space.
+func (t Triangle) Centroid() Vec {
+	return Scale(1.0/3.0, Add(Add(t[0], t[1]), t[2]))
+}
+
+// Normal returns the vector with direction
+// perpendicular to the Triangle's face and magnitude
+// twice that of the Triangle's area. The ordering
+// of the triangle vertices decides the normal's resulting
+// direction. The returned vector is not normalized.
+func (t Triangle) Normal() Vec {
+	s1, s2, _ := t.sides()
+	return Cross(s1, s2)
+}
+
+// IsDegenerate returns true if all of triangle's vertices are
+// within tol distance of its longest side.
+func (t Triangle) IsDegenerate(tol float64) bool {
+	longIdx := t.longIdx()
+	// calculate vertex distance from longest side
+	ln := line{t[longIdx], t[(longIdx+1)%3]}
+	dist := ln.distance(t[(longIdx+2)%3])
+	return dist <= tol
+}
+
+// longIdx returns index of the longest side. The sides
+// of the triangles are are as follows:
+//  - Side 0 formed by vertices 0 and 1
+//  - Side 1 formed by vertices 1 and 2
+//  - Side 2 formed by vertices 0 and 2
+func (t Triangle) longIdx() int {
+	sides := [3]Vec{Sub(t[1], t[0]), Sub(t[2], t[1]), Sub(t[0], t[2])}
+	len2 := [3]float64{Norm2(sides[0]), Norm2(sides[1]), Norm2(sides[2])}
+	longLen := len2[0]
+	longIdx := 0
+	if len2[1] > longLen {
+		longLen = len2[1]
+		longIdx = 1
+	}
+	if len2[2] > longLen {
+		longIdx = 2
+	}
+	return longIdx
+}
+
+// Area returns the surface area of the triangle.
+func (t Triangle) Area() float64 {
+	// Heron's Formula, see https://en.wikipedia.org/wiki/Heron%27s_formula.
+	// Also see William M. Kahan (24 March 2000). "Miscalculating Area and Angles of a Needle-like Triangle"
+	// for more discussion. https://people.eecs.berkeley.edu/~wkahan/Triangle.pdf.
+	a, b, c := t.orderedLengths()
+	A := (c + (b + a)) * (a - (c - b))
+	A *= (a + (c - b)) * (c + (b - a))
+	return math.Sqrt(A) / 4
+}
+
+// orderedLengths returns the lengths of the sides of the triangle such that
+// a ≤ b ≤ c.
+func (t Triangle) orderedLengths() (a, b, c float64) {
+	s1, s2, s3 := t.sides()
+	l1 := Norm(s1)
+	l2 := Norm(s2)
+	l3 := Norm(s3)
+	// sort-3
+	if l2 < l1 {
+		l1, l2 = l2, l1
+	}
+	if l3 < l2 {
+		l2, l3 = l3, l2
+		if l2 < l1 {
+			l1, l2 = l2, l1
+		}
+	}
+	return l1, l2, l3
+}
+
+// sides returns vectors for each of the sides of t.
+func (t Triangle) sides() (Vec, Vec, Vec) {
+	return Sub(t[1], t[0]), Sub(t[2], t[1]), Sub(t[0], t[2])
+}
+
+// line is an infinite 3D line
+// defined by two points on the line.
+type line [2]Vec
+
+// vecOnLine takes a value between 0 and 1 to linearly
+// interpolate a point on the line.
+//  vecOnLine(0) returns l[0]
+//  vecOnLine(1) returns l[1]
+func (l line) vecOnLine(t float64) Vec {
+	lineDir := Sub(l[1], l[0])
+	return Add(l[0], Scale(t, lineDir))
+}
+
+// distance returns the minimum euclidean distance of point p
+// to the line.
+func (l line) distance(p Vec) float64 {
+	// https://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html
+	num := Norm(Cross(Sub(p, l[0]), Sub(p, l[1])))
+	return num / Norm(Sub(l[1], l[0]))
+}