package/gonum
1
0
Fork 0
mirror of https://github.com/gonum/gonum.git synced 2025-10-17 20:51:06 +08:00
Code Issues Actions Packages Projects Releases Wiki Activity

spatial/r3: remove duplicated API

Browse Source
This commit is contained in:
Dan Kortschak
2022-02-25 18:51:02 +10:30
parent f4e711476a
commit 809af93335
7 changed files with 110 additions and 131 deletions
Download Patch File Download Diff File Expand all files Collapse all files

27
spatial/r3/mat.go
View File

@@ -4,10 +4,7 @@
package r3
import (
"gonum.org/v1/gonum/mat"
"gonum.org/v1/gonum/num/quat"
)
import "gonum.org/v1/gonum/mat"
// Mat represents a 3×3 matrix. Useful for rotation matrices and such.
// The zero value is usable as the 3×3 zero matrix.
@@ -185,28 +182,6 @@ func (m *Mat) Outer(alpha float64, x, y Vec) {
m.Set(2, 2, az*y.Z)
}
// RotationFromQuat converts the quaternion q to the corresponding orthogonal
// matrix and stores the result in the receiver. q should be normalized but this
// assumption is not checked.
//
// See https://en.wikipedia.org/wiki/Rotation_matrix#Quaternion for more
// details.
func (m *Mat) RotationFromQuat(q quat.Number) {
w, x, y, z := q.Real, q.Imag, q.Jmag, q.Kmag
x2, y2, z2 := x*x, y*y, z*z
xw, xy, xz := x*w, x*y, x*z
yw, zw, yz := y*w, z*w, y*z
m.Set(0, 0, 1-2*y2-2*z2)
m.Set(0, 1, 2*xy-2*zw)
m.Set(0, 2, 2*xz+2*yw)
m.Set(1, 0, 2*xy+2*zw)
m.Set(1, 1, 1-2*x2-2*z2)
m.Set(1, 2, 2*yz-2*xw)
m.Set(2, 0, 2*xz-2*yw)
m.Set(2, 1, 2*yz+2*xw)
m.Set(2, 2, 1-2*x2-2*y2)
}
// Det calculates the determinant of the receiver using the following formula
// ⎡a b c⎤
// m = ⎢d e f⎥

22
spatial/r3/mat_safe.go
View File

@@ -136,3 +136,25 @@ func arrayFrom(vals []float64) *array {
copy(a[:], vals)
return &a
}
// 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 (r Rotation) Mat() *Mat {
w, i, j, k := r.Real, r.Imag, r.Jmag, r.Kmag
ii := 2 * i * i
jj := 2 * j * j
kk := 2 * k * k
wi := 2 * w * i
wj := 2 * w * j
wk := 2 * w * k
ij := 2 * i * j
jk := 2 * j * k
ki := 2 * k * i
return &Mat{&array{
1 - (jj + kk), ij - wk, ki + wj,
ij + wk, 1 - (ii + kk), jk - wi,
ki - wj, jk + wi, 1 - (ii + jj),
}}
}

10
spatial/r3/mat_test.go
View File

@@ -188,7 +188,7 @@ func TestOuter(t *testing.T) {
}
}
func TestRotationFromQuat(t *testing.T) {
func TestRotationMat(t *testing.T) {
const tol = 1e-14
rnd := rand.New(rand.NewSource(1))
@@ -198,8 +198,7 @@ func TestRotationFromQuat(t *testing.T) {
q = quat.Scale(1/quat.Abs(q), q)
// Convert it to a rotation matrix R.
var r Mat
r.RotationFromQuat(q)
r := Rotation(q).Mat()
// Check that the matrix has the determinant approximately equal to 1.
diff := math.Abs(r.Det() - 1)
@@ -224,9 +223,10 @@ func TestRotationFromQuat(t *testing.T) {
func BenchmarkQuat(b *testing.B) {
rnd := rand.New(rand.NewSource(1))
m := NewMat(nil)
for i := 0; i < b.N; i++ {
q := quat.Number{Real: rnd.Float64(), Imag: rnd.Float64(), Jmag: rnd.Float64(), Kmag: rnd.Float64()}
m.RotationFromQuat(q)
if Rotation(q).Mat() == nil {
b.Fatal("nil return")
}
}
}

22
spatial/r3/mat_unsafe.go
View File

@@ -117,6 +117,28 @@ func (m *Mat) RawMatrix() blas64.General {
return blas64.General{Rows: 3, Cols: 3, Data: m.slice(), Stride: 3}
}
// 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 (r Rotation) Mat() *Mat {
w, i, j, k := r.Real, r.Imag, r.Jmag, r.Kmag
ii := 2 * i * i
jj := 2 * j * j
kk := 2 * k * k
wi := 2 * w * i
wj := 2 * w * j
wk := 2 * w * k
ij := 2 * i * j
jk := 2 * j * k
ki := 2 * k * i
return &Mat{&array{
{1 - (jj + kk), ij - wk, ki + wj},
{ij + wk, 1 - (ii + kk), jk - wi},
{ki - wj, jk + wi, 1 - (ii + jj)},
}}
}
func arrayFrom(vals []float64) *array {
// TODO(kortschak): Use array conversion when go1.16 is no longer supported.
return (*array)(unsafe.Pointer(&vals[0]))

58
spatial/r3/rotation.go Normal file
View File

@@ -0,0 +1,58 @@
// Copyright ©2021 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"
"gonum.org/v1/gonum/num/quat"
)
// TODO: possibly useful additions to the current rotation API:
// - create rotations from Euler angles (NewRotationFromEuler?)
// - create rotations from rotation matrices (NewRotationFromMatrix?)
// - return the equivalent Euler angles from a Rotation
//
// Euler angles have issues (see [1] for a discussion).
// We should think carefully before adding them in.
// [1]: http://www.euclideanspace.com/maths/geometry/rotations/conversions/quaternionToEuler/
// Rotation describes a rotation in space.
type Rotation quat.Number
// NewRotation creates a rotation by alpha, around axis.
func NewRotation(alpha float64, axis Vec) Rotation {
if alpha == 0 {
return Rotation{Real: 1}
}
q := raise(axis)
sin, cos := math.Sincos(0.5 * alpha)
q = quat.Scale(sin/quat.Abs(q), q)
q.Real += cos
if len := quat.Abs(q); len != 1 {
q = quat.Scale(1/len, q)
}
return Rotation(q)
}
// Rotate returns p rotated according to the parameters used to construct
// the receiver.
func (r Rotation) Rotate(p Vec) Vec {
if r.isIdentity() {
return p
}
qq := quat.Number(r)
pp := quat.Mul(quat.Mul(qq, raise(p)), quat.Conj(qq))
return Vec{X: pp.Imag, Y: pp.Jmag, Z: pp.Kmag}
}
func (r Rotation) isIdentity() bool {
return r == Rotation{Real: 1}
}
func raise(p Vec) quat.Number {
return quat.Number{Imag: p.X, Jmag: p.Y, Kmag: p.Z}
}

100
spatial/r3/vector.go
View File

@@ -4,13 +4,7 @@
package r3
import (
"math"
"gonum.org/v1/gonum/blas/blas64"
"gonum.org/v1/gonum/mat"
"gonum.org/v1/gonum/num/quat"
)
import "math"
// Vec is a 3D vector.
type Vec struct {
@@ -93,95 +87,3 @@ func Cos(p, q Vec) float64 {
type Box struct {
Min, Max Vec
}
// TODO: possibly useful additions to the current rotation API:
// - create rotations from Euler angles (NewRotationFromEuler?)
// - create rotations from rotation matrices (NewRotationFromMatrix?)
// - return the equivalent Euler angles from a Rotation
//
// Euler angles have issues (see [1] for a discussion).
// We should think carefully before adding them in.
// [1]: http://www.euclideanspace.com/maths/geometry/rotations/conversions/quaternionToEuler/
// Rotation describes a rotation in space.
type Rotation quat.Number
// NewRotation creates a rotation by alpha, around axis.
func NewRotation(alpha float64, axis Vec) Rotation {
if alpha == 0 {
return Rotation{Real: 1}
}
q := raise(axis)
sin, cos := math.Sincos(0.5 * alpha)
q = quat.Scale(sin/quat.Abs(q), q)
q.Real += cos
if len := quat.Abs(q); len != 1 {
q = quat.Scale(1/len, q)
}
return Rotation(q)
}
// Rotate returns p rotated according to the parameters used to construct
// the receiver.
func (r Rotation) Rotate(p Vec) Vec {
if r.isIdentity() {
return p
}
qq := quat.Number(r)
pp := quat.Mul(quat.Mul(qq, raise(p)), quat.Conj(qq))
return Vec{X: pp.Imag, Y: pp.Jmag, Z: pp.Kmag}
}
func (r Rotation) isIdentity() bool {
return r == Rotation{Real: 1}
}
func raise(p Vec) quat.Number {
return quat.Number{Imag: p.X, Jmag: p.Y, Kmag: p.Z}
}
// Matrix 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 (r Rotation) Matrix() mat.Matrix {
re, im, jm, km := r.Real, r.Imag, r.Jmag, r.Kmag
im2 := im * im
jm2 := jm * jm
km2 := km * km
rim := re * im
rjm := re * jm
rkm := re * km
ijm := im * jm
jkm := jm * km
kim := km * im
return &matrix{
1 - 2*(jm2+km2), 2 * (ijm - rkm), 2 * (kim + rjm),
2 * (ijm + rkm), 1 - 2*(im2+km2), 2 * (jkm - rim),
2 * (kim - rjm), 2 * (jkm + rim), 1 - 2*(im2+jm2),
}
}
// matrix is a 3×3 rotation matrix.
type matrix [9]float64
var (
_ mat.Matrix = (*matrix)(nil)
_ mat.RawMatrixer = (*matrix)(nil)
)
func (m *matrix) At(i, j int) float64 {
if uint(i) >= 3 {
panic(mat.ErrRowAccess)
}
if uint(j) >= 3 {
panic(mat.ErrColAccess)
}
return m[i*3+j]
}
func (m *matrix) Dims() (r, c int) { return 3, 3 }
func (m *matrix) T() mat.Matrix { return mat.Transpose{Matrix: m} }
func (m *matrix) RawMatrix() blas64.General {
return blas64.General{Rows: 3, Cols: 3, Data: m[:], Stride: 3}
}

2
spatial/r3/vector_test.go
View File

@@ -243,7 +243,7 @@ func TestRotate(t *testing.T) {
}
var gotv mat.VecDense
gotv.MulVec(NewRotation(test.alpha, test.axis).Matrix(), vecDense(test.v))
gotv.MulVec(NewRotation(test.alpha, test.axis).Mat(), vecDense(test.v))
got = vec(gotv)
if !vecApproxEqual(got, test.want, tol) {
t.Errorf(