package r3

import "gonum.org/v1/gonum/spatial/r3"

Package r3 provides 3D vectors and boxes and operations on them.

Example (Slerp)

Spherically interpolate between two quaternions to obtain a rotation.

Code:play  

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 %+.2f\n", t, rotated)
	}

}

Output:

0.00 {+1.00 -0.00 +1.41}
0.10 {+1.00 -0.33 +1.38}
0.20 {+1.00 -0.64 +1.26}
0.30 {+1.00 -0.92 +1.08}
0.40 {+1.00 -1.14 +0.83}
0.50 {+1.00 -1.31 +0.54}
0.60 {+1.00 -1.40 +0.22}
0.70 {+1.00 -1.41 -0.11}
0.80 {+1.00 -1.34 -0.44}
0.90 {+1.00 -1.21 -0.74}
1.00 {+1.00 -1.00 -1.00}

Index

Examples

Functions

func Cos 

func Cos(p, q Vec) float64

Cos returns the cosine of the opening angle between p and q.

func Divergence 

func Divergence(p, step Vec, field func(Vec) Vec) float64

Divergence returns the divergence of the vector field at the point p, approximated using finite differences with the given step sizes.

func Dot 

func Dot(p, q Vec) float64

Dot returns the dot product p·q.

func Norm 

func Norm(p Vec) float64

Norm returns the Euclidean norm of p 

|p| = sqrt(p_x^2 + p_y^2 + p_z^2).

func Norm2 

func Norm2(p Vec) float64

Norm2 returns the Euclidean squared norm of p 

|p|^2 = p_x^2 + p_y^2 + p_z^2.

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 

func NewBox(x0, y0, z0, x1, y1, z1 float64) Box

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 

func (a Box) Add(v Vec) Box

Add adds v to the bounding box components. It is the equivalent of translating the Box by v.

func (Box) Canon 

func (a Box) Canon() Box

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) Center 

func (a Box) Center() Vec

Center returns the center of the Box.

func (Box) Contains 

func (a Box) Contains(v Vec) bool

Contains returns true if v is contained within the bounds of the Box.

func (Box) Empty 

func (a Box) Empty() bool

Empty returns true if a Box's volume is zero or if a Min component is greater than its Max component.

func (Box) Scale 

func (a Box) Scale(scale Vec) Box

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) Size 

func (a Box) Size() Vec

Size returns the size of the Box.

func (Box) Union 

func (a Box) Union(b Box) Box

Union returns a box enclosing both the receiver and argument Boxes.

func (Box) Vertices 

func (a Box) Vertices() []Vec

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 

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 Eye 

func Eye() *Mat

Eye returns the 3×3 Identity matrix

func NewMat 

func NewMat(val []float64) *Mat

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 Skew 

func Skew(v Vec) (M *Mat)

Skew returns the 3×3 skew symmetric matrix (right hand system) of v. 

                ⎡ 0 -z  y⎤
Skew({x,y,z}) = ⎢ z  0 -x⎥
                ⎣-y  x  0⎦

Deprecated: use Mat.Skew()

func (*Mat) Add 

func (m *Mat) Add(a, b mat.Matrix)

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 

func (m *Mat) At(i, j int) float64

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 

func (m *Mat) CloneFrom(a mat.Matrix)

CloneFrom makes a copy of a into the receiver m. Mat expects a 3×3 input matrix.

func (*Mat) Det 

func (m *Mat) Det() float64

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 

func (m *Mat) Dims() (r, c int)

Dims returns the number of rows and columns of this matrix. This method will always return 3×3 for a Mat.

func (*Mat) Hessian 

func (m *Mat) Hessian(p, step Vec, field func(Vec) float64)

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 

func (m *Mat) Jacobian(p, step Vec, field func(Vec) Vec)

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 

func (m *Mat) Mul(a, b mat.Matrix)

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) MulVec 

func (m *Mat) MulVec(v Vec) Vec

MulVec returns the matrix-vector product M⋅v.

func (*Mat) MulVecTrans 

func (m *Mat) MulVecTrans(v Vec) Vec

MulVecTrans returns the matrix-vector product Mᵀ⋅v.

func (*Mat) Outer 

func (m *Mat) Outer(alpha float64, x, y Vec)

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 

func (m *Mat) RawMatrix() blas64.General

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 

func (m *Mat) Scale(f float64, a mat.Matrix)

Scale multiplies the elements of a by f, placing the result in the receiver.

See the mat.Scaler interface for more information.

func (*Mat) Set 

func (m *Mat) Set(i, j int, v float64)

Set sets the element at row i, column j to the value v.

func (*Mat) Skew 

func (m *Mat) Skew(v Vec)

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 

func (m *Mat) Sub(a, b mat.Matrix)

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 

func (m *Mat) T() mat.Matrix

T returns the transpose of Mat. Changes in the receiver will be reflected in the returned matrix.

func (*Mat) VecCol 

func (m *Mat) VecCol(j int) Vec

VecCol returns the elements in the jth column of the receiver.

func (*Mat) VecRow 

func (m *Mat) VecRow(i int) Vec

VecRow returns the elements in the ith row of the receiver.

type Rotation 

type Rotation quat.Number

Rotation describes a rotation in space.

Example (EulerAngles)

Code:play  

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 

func NewRotation(alpha float64, axis Vec) Rotation

NewRotation creates a rotation by alpha, around axis.

func (Rotation) Mat 

func (r Rotation) Mat() *Mat

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.

func (Rotation) Rotate 

func (r Rotation) Rotate(p Vec) Vec

Rotate returns p rotated according to the parameters used to construct the receiver.

type Triangle 

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)

Code:play  

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},
		}
}

func (Triangle) Area 

func (t Triangle) Area() float64

Area returns the surface area of the triangle.

func (Triangle) Centroid 

func (t Triangle) Centroid() Vec

Centroid returns the intersection of the three medians of the triangle as a point in space.

func (Triangle) IsDegenerate 

func (t Triangle) IsDegenerate(tol float64) bool

IsDegenerate returns true if all of triangle's vertices are within tol distance of its longest side.

func (Triangle) Normal 

func (t Triangle) Normal() Vec

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.

type Vec 

type Vec struct {
	X, Y, Z float64
}

Vec is a 3D vector.

func Add 

func Add(p, q Vec) Vec

Add returns the vector sum of p and q.

func Cross 

func Cross(p, q Vec) Vec

Cross returns the cross product p×q.

func Gradient 

func Gradient(p, step Vec, field func(Vec) float64) Vec

Gradient returns the gradient of the scalar field at the point p, approximated using finite differences with the given step sizes.

func Rotate 

func Rotate(p Vec, alpha float64, axis Vec) Vec

Rotate returns a new vector, rotated by alpha around the provided axis.

func Scale 

func Scale(f float64, p Vec) Vec

Scale returns the vector p scaled by f.

func Sub 

func Sub(p, q Vec) Vec

Sub returns the vector sum of p and -q.

func Unit 

func Unit(p Vec) Vec

Unit returns the unit vector colinear to p. Unit returns {NaN,NaN,NaN} for the zero vector.

Source Files

box.go doc.go mat.go mat_unsafe.go rotation.go triangle.go vector.go

Version
v0.15.1 (latest)
Published
Aug 16, 2024
Platform
linux/amd64
Imports
5 packages
Last checked
12 hours ago –

Tools for package owners.