Documentation
¶
Overview ¶
Package r3 provides 3D vectors and boxes and operations on them.
Example (Slerp) ¶
Spherically interpolate between two quaternions to obtain a rotation.
package main
import (
"fmt"
"math"
"gonum.org/v1/gonum/num/quat"
"gonum.org/v1/gonum/spatial/r3"
)
// slerp returns the spherical interpolation between q0 and q1
// for t in [0,1]; 0 corresponds to q0 and 1 corresponds to q1.
func slerp(r0, r1 r3.Rotation, t float64) r3.Rotation {
q0 := quat.Number(r0)
q1 := quat.Number(r1)
// Based on Simo Särkkä "Notes on Quaternions" Eq. 35
// p(t) = (q1 ∗ q0^−1) ^ t ∗ q0
// https://users.aalto.fi/~ssarkka/pub/quat.pdf
q1 = quat.Mul(q1, quat.Inv(q0))
q1 = quat.PowReal(q1, t)
return r3.Rotation(quat.Mul(q1, q0))
}
// Spherically interpolate between two quaternions to obtain a rotation.
func main() {
const steps = 10
// An initial rotation of pi/4 around the x-axis (45 degrees).
initialRot := r3.NewRotation(math.Pi/4, r3.Vec{X: 1})
// Final rotation is pi around the x-axis (180 degrees).
finalRot := r3.NewRotation(math.Pi, r3.Vec{X: 1})
// The vector we are rotating is (1, 1, 1).
// The result should then be (1, -1, -1) when t=1 (finalRot) since we invert the y and z axes.
v := r3.Vec{X: 1, Y: 1, Z: 1}
for i := 0.0; i <= steps; i++ {
t := i / steps
rotated := slerp(initialRot, finalRot, t).Rotate(v)
fmt.Printf("%.2f %.4g\n", t, rotated)
}
}
Output: 0.00 {1 -1.11e-16 1.414} 0.10 {1 -0.3301 1.375} 0.20 {1 -0.642 1.26} 0.30 {1 -0.9185 1.075} 0.40 {1 -1.144 0.8313} 0.50 {1 -1.307 0.5412} 0.60 {1 -1.397 0.2212} 0.70 {1 -1.41 -0.111} 0.80 {1 -1.345 -0.437} 0.90 {1 -1.206 -0.7389} 1.00 {1 -1 -1}
Index ¶
- func Cos(p, q Vec) float64
- func Divergence(p, step Vec, field func(Vec) Vec) float64
- func Dot(p, q Vec) float64
- func Norm(p Vec) float64
- func Norm2(p Vec) float64
- type Box
- type Mat
- func (m *Mat) Add(a, b mat.Matrix)
- func (m *Mat) At(i, j int) float64
- func (m *Mat) CloneFrom(a mat.Matrix)
- func (m *Mat) Det() float64
- func (m *Mat) Dims() (r, c int)
- func (m *Mat) Hessian(p, step Vec, field func(Vec) float64)
- func (m *Mat) Jacobian(p, step Vec, field func(Vec) Vec)
- func (m *Mat) Mul(a, b mat.Matrix)
- func (m *Mat) MulVec(v Vec) Vec
- func (m *Mat) MulVecTrans(v Vec) Vec
- func (m *Mat) Outer(alpha float64, x, y Vec)
- func (m *Mat) RawMatrix() blas64.General
- func (m *Mat) Scale(f float64, a mat.Matrix)
- func (m *Mat) Set(i, j int, v float64)
- func (m *Mat) Skew(v Vec)
- func (m *Mat) Sub(a, b mat.Matrix)
- func (m *Mat) T() mat.Matrix
- func (m *Mat) VecCol(j int) Vec
- func (m *Mat) VecRow(i int) Vec
- type Rotation
- type Triangle
- type Vec
Examples ¶
Constants ¶
This section is empty.
Variables ¶
This section is empty.
Functions ¶
func Divergence ¶ added in v0.12.0
Divergence returns the divergence of the vector field at the point p, approximated using finite differences with the given step sizes.
Types ¶
type Box ¶
type Box struct {
Min, Max Vec
}
Box is a 3D bounding box. Well formed Boxes Min components are smaller than Max components.
func NewBox ¶ added in v0.12.0
NewBox is shorthand for Box{Min:Vec{x0,y0,z0}, Max:Vec{x1,y1,z1}}. The sides are swapped so that the resulting Box is well formed.
func (Box) Add ¶ added in v0.12.0
Add adds v to the bounding box components. It is the equivalent of translating the Box by v.
func (Box) Canon ¶ added in v0.12.0
Canon returns the canonical version of a. The returned Box has minimum and maximum coordinates swapped if necessary so that it is well-formed.
func (Box) Contains ¶ added in v0.12.0
Contains returns true if v is contained within the bounds of the Box.
func (Box) Empty ¶ added in v0.12.0
IsEmpty returns true if a Box's volume is zero or if a Min component is greater than its Max component.
func (Box) Scale ¶ added in v0.12.0
Scale returns a new Box scaled by a size vector around its center. The scaling is done element wise which is to say the Box's X dimension is scaled by scale.X. Negative elements of scale are interpreted as zero.
func (Box) Union ¶ added in v0.12.0
Union returns a box enclosing both the receiver and argument Boxes.
func (Box) Vertices ¶ added in v0.12.0
Vertices returns a slice of the 8 vertices corresponding to each of the Box's corners.
Ordering of vertices 0-3 is CCW in the XY plane starting at box minimum. Ordering of vertices 4-7 is CCW in the XY plane starting at box minimum for X and Y values and maximum Z value.
Edges for the box can be constructed with the following indices:
edges := [12][2]int{
{0, 1}, {1, 2}, {2, 3}, {3, 0},
{4, 5}, {5, 6}, {6, 7}, {7, 4},
{0, 4}, {1, 5}, {2, 6}, {3, 7},
}
type Mat ¶ added in v0.11.0
type Mat struct {
// contains filtered or unexported fields
}
Mat represents a 3×3 matrix. Useful for rotation matrices and such. The zero value is usable as the 3×3 zero matrix.
func NewMat ¶ added in v0.11.0
NewMat returns a new 3×3 matrix Mat type and populates its elements with values passed as argument in row-major form. If val argument is nil then NewMat returns a matrix filled with zeros.
func (*Mat) Add ¶ added in v0.11.0
Add adds a and b element-wise, placing the result in the receiver. Add will panic if the two matrices do not have the same shape.
func (*Mat) At ¶ added in v0.11.0
At returns the value of a matrix element at row i, column j. At expects indices in the range [0,2]. It will panic if i or j are out of bounds for the matrix.
func (*Mat) CloneFrom ¶ added in v0.11.0
CloneFrom makes a copy of a into the receiver m. Mat expects a 3×3 input matrix.
func (*Mat) Det ¶ added in v0.11.0
Det calculates the determinant of the receiver using the following formula
⎡a b c⎤
m = ⎢d e f⎥
⎣g h i⎦
det(m) = a(ei − fh) − b(di − fg) + c(dh − eg)
func (*Mat) Dims ¶ added in v0.11.0
Dims returns the number of rows and columns of this matrix. This method will always return 3×3 for a Mat.
func (*Mat) Hessian ¶ added in v0.12.0
Hessian sets the receiver to the Hessian matrix of the scalar field at the point p, approximated using finite differences with the given step sizes. The field is evaluated at points in the area surrounding p by adding at most 2 components of step to p. Hessian expects the field's second partial derivatives are all continuous for correct results.
func (*Mat) Jacobian ¶ added in v0.12.0
Jacobian sets the receiver to the Jacobian matrix of the vector field at the point p, approximated using finite differences with the given step sizes.
The Jacobian matrix J is the matrix of all first-order partial derivatives of f. If f maps an 3-dimensional vector x to an 3-dimensional vector y = f(x), J is a 3×3 matrix whose elements are given as
J_{i,j} = ∂f_i/∂x_j,
or expanded out
[ ∂f_1/∂x_1 ∂f_1/∂x_2 ∂f_1/∂x_3 ]
J = [ ∂f_2/∂x_1 ∂f_2/∂x_2 ∂f_2/∂x_3 ]
[ ∂f_3/∂x_1 ∂f_3/∂x_2 ∂f_3/∂x_3 ]
Jacobian expects the field's first order partial derivatives are all continuous for correct results.
func (*Mat) Mul ¶ added in v0.11.0
Mul takes the matrix product of a and b, placing the result in the receiver. If the number of columns in a does not equal 3, Mul will panic.
func (*Mat) MulVecTrans ¶ added in v0.11.0
MulVecTrans returns the matrix-vector product Mᵀ⋅v.
func (*Mat) Outer ¶ added in v0.11.0
Outer calculates the outer product of the vectors x and y, where x and y are treated as column vectors, and stores the result in the receiver.
m = alpha * x * yᵀ
func (*Mat) RawMatrix ¶ added in v0.11.0
RawMatrix returns the blas representation of the matrix with the backing data of this matrix. Changes to the returned matrix will be reflected in the receiver.
func (*Mat) Scale ¶ added in v0.11.0
Scale multiplies the elements of a by f, placing the result in the receiver.
See the mat.Scaler interface for more information.
func (*Mat) Skew ¶ added in v0.12.0
Skew sets the receiver to the 3×3 skew symmetric matrix (right hand system) of v.
⎡ 0 -z y⎤
Skew({x,y,z}) = ⎢ z 0 -x⎥
⎣-y x 0⎦
func (*Mat) Sub ¶ added in v0.11.0
Sub subtracts the matrix b from a, placing the result in the receiver. Sub will panic if the two matrices do not have the same shape.
func (*Mat) T ¶ added in v0.11.0
T returns the transpose of Mat. Changes in the receiver will be reflected in the returned matrix.
type Rotation ¶ added in v0.9.0
Rotation describes a rotation in space.
Example (EulerAngles) ¶
package main
import (
"fmt"
"math"
"gonum.org/v1/gonum/num/quat"
"gonum.org/v1/gonum/spatial/r3"
)
// euler returns an r3.Rotation that corresponds to the Euler
// angles alpha, beta and gamma which are rotations around the x,
// y and z axes respectively. The order of rotations is x, y, z;
// there are many conventions for this ordering.
func euler(alpha, beta, gamma float64) r3.Rotation {
// Note that this function can be algebraically simplified
// to reduce floating point operations, but is left in this
// form for clarity.
var rot1, rot2, rot3 quat.Number
rot1.Imag, rot1.Real = math.Sincos(alpha / 2) // x-axis rotation
rot2.Jmag, rot2.Real = math.Sincos(beta / 2) // y-axis rotation
rot3.Kmag, rot3.Real = math.Sincos(gamma / 2) // z-axis rotation
return r3.Rotation(quat.Mul(rot3, quat.Mul(rot2, rot1))) // order of rotations
}
func main() {
// It is possible to interconvert between the quaternion representation
// of a rotation and Euler angles, but this leads to problems.
//
// The first of these is that there are a variety of conventions for
// application of the rotations.
//
// The more serious consequence of using Euler angles is that it is
// possible to put the rotation system into a singularity which results
// in loss of degrees of freedom and so causes gimbal lock. This happens
// when the second axis to be rotated around is rotated to 𝝿/2.
//
// See https://en.wikipedia.org/wiki/Euler_angles for more details.
pt := r3.Vec{1, 0, 0}
// For the Euler conversion function in this example, the second rotation
// is around the y-axis.
const singularY = math.Pi / 2
arb := math.Pi / 4
fmt.Printf("rotate around x-axis: %.2f\n", euler(arb, 0, 0).Rotate(pt))
fmt.Printf("rotate around y-axis: %.2f\n", euler(0, arb, 0).Rotate(pt))
fmt.Printf("rotate around z-axis: %.2f\n", euler(0, 0, arb).Rotate(pt))
fmt.Printf("rotate around x+y-axes: %.2f\n", euler(arb, arb, 0).Rotate(pt))
fmt.Printf("rotate around x+z-axes: %.2f\n", euler(arb, 0, arb).Rotate(pt))
fmt.Printf("rotate around y+z-axes: %.2f\n", euler(0, arb, arb).Rotate(pt))
fmt.Printf("rotate around y-axis to singularity: %.2f\n", euler(0, singularY, 0).Rotate(pt))
fmt.Printf("rotate around x+y-axes with singularity → gimbal lock: %.2f\n", euler(arb, singularY, 0).Rotate(pt))
fmt.Printf("rotate around z+y-axes with singularity → gimbal lock: %.2f\n", euler(0, singularY, arb).Rotate(pt))
fmt.Printf("rotate around all-axes with singularity → gimbal lock: %.2f\n", euler(arb, singularY, arb).Rotate(pt))
}
Output: rotate around x-axis: {1.00 0.00 0.00} rotate around y-axis: {0.71 0.00 -0.71} rotate around z-axis: {0.71 0.71 0.00} rotate around x+y-axes: {0.71 0.00 -0.71} rotate around x+z-axes: {0.71 0.71 0.00} rotate around y+z-axes: {0.50 0.50 -0.71} rotate around y-axis to singularity: {0.00 0.00 -1.00} rotate around x+y-axes with singularity → gimbal lock: {0.00 0.00 -1.00} rotate around z+y-axes with singularity → gimbal lock: {0.00 0.00 -1.00} rotate around all-axes with singularity → gimbal lock: {0.00 0.00 -1.00}
func NewRotation ¶ added in v0.9.0
NewRotation creates a rotation by alpha, around axis.
func (Rotation) Mat ¶ added in v0.11.0
Mat returns a 3×3 rotation matrix corresponding to the receiver. It may be used to perform rotations on a 3-vector or to apply the rotation to a 3×n matrix of column vectors. If the receiver is not a unit quaternion, the returned matrix will not be a pure rotation.
type Triangle ¶ added in v0.12.0
type Triangle [3]Vec
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.
Example (Icosphere) ¶
package main
import (
"fmt"
"gonum.org/v1/gonum/spatial/r3"
)
func main() {
// 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},
}
}
Output:
func (Triangle) Centroid ¶ added in v0.12.0
Centroid returns the intersection of the three medians of the triangle as a point in space.
func (Triangle) IsDegenerate ¶ added in v0.12.0
IsDegenerate returns true if all of triangle's vertices are within tol distance of its longest side.
type Vec ¶
type Vec struct {
X, Y, Z float64
}
Vec is a 3D vector.
func Gradient ¶ added in v0.12.0
Gradient returns the gradient of the scalar field at the point p, approximated using finite differences with the given step sizes.
func Rotate ¶ added in v0.9.0
Rotate returns a new vector, rotated by alpha around the provided axis.
Source Files