Some new additions, some leaner implementations and one breaking chan…

…ge(!) (#35)

* Breaking change(!) cross vec2 (float32 float64) now return a scalar value, not a new vector. Cross-product in 2D is not well-defined but it is changed to the 2D version stated at https://mathworld.wolfram.com/CrossProduct.html .

Add Sinus and Cosine for angle between vectors (vec2 vec3 float64 float32)
Changed implementation for Angle (vec2 vec3 float32 float64) to utilize Cosine implementation (with math.Acos).
Changed implementation for left and right winding (vec2 float32 float64)

Duplicated test file for quaternion (float32 float64)

* doc fix

* rename fmath.Sqrtf to Sqrt and fix build for non amd64

* fix quaternion tests. They were failing due to precision errors in float values.

---------

Co-authored-by: Erik Unger <[email protected]>
Co-authored-by: Erik Unger <[email protected]>

float64/quaternion/quaternion_test.go

@@ -0,0 +1,92 @@
+package quaternion
+
+import (
+	"fmt"
+	"math"
+	"testing"
+
+	"github.com/ungerik/go3d/float64/vec3"
+)
+
+// RotateVec3 rotates v by the rotation represented by the quaternion.
+func rotateAndNormalizeVec3(quat *T, v *vec3.T) {
+	qv := T{v[0], v[1], v[2], 0}
+	inv := quat.Inverted()
+	q := Mul3(quat, &qv, &inv)
+	v[0] = q[0]
+	v[1] = q[1]
+	v[2] = q[2]
+}
+
+func TestQuaternionRotateVec3(t *testing.T) {
+	eularAngles := []vec3.T{
+		{90, 20, 21},
+		{-90, 0, 0},
+		{28, 1043, -38},
+	}
+	vecs := []vec3.T{
+		{2, 3, 4},
+		{1, 3, -2},
+		{-6, 2, 9},
+	}
+	for _, vec := range vecs {
+		for _, eularAngle := range eularAngles {
+			func() {
+				q := FromEulerAngles(eularAngle[1]*math.Pi/180.0, eularAngle[0]*math.Pi/180.0, eularAngle[2]*math.Pi/180.0)
+				vec_r1 := vec
+				vec_r2 := vec
+				magSqr := vec_r1.LengthSqr()
+				rotateAndNormalizeVec3(&q, &vec_r2)
+				q.RotateVec3(&vec_r1)
+				vecd := q.RotatedVec3(&vec)
+				magSqr2 := vec_r1.LengthSqr()
+
+				if !vecd.PracticallyEquals(&vec_r1, 0.000000000000001) {
+					t.Logf("test case %v rotates %v failed - vector rotation: %+v, %+v\n", eularAngle, vec, vecd, vec_r1)
+					t.Fail()
+				}
+
+				angle := vec3.Angle(&vec_r1, &vec_r2)
+				length := math.Abs(magSqr - magSqr2)
+
+				if angle > 0.0000001 {
+					t.Logf("test case %v rotates %v failed - angle difference to large\n", eularAngle, vec)
+					t.Logf("vectors: %+v, %+v\n", vec_r1, vec_r2)
+					t.Logf("angle: %v\n", angle)
+					t.Fail()
+				}
+
+				if length > 0.000000000001 {
+					t.Logf("test case %v rotates %v failed - squared length difference to large\n", eularAngle, vec)
+					t.Logf("vectors: %+v %+v\n", vec_r1, vec_r2)
+					t.Logf("squared lengths: %v, %v\n", magSqr, magSqr2)
+					t.Fail()
+				}
+			}()
+		}
+	}
+}
+
+func TestToEulerAngles(t *testing.T) {
+	specialValues := []float64{-5, -math.Pi, -2, -math.Pi / 2, 0, math.Pi / 2, 2.4, math.Pi, 3.9}
+	for _, x := range specialValues {
+		for _, y := range specialValues {
+			for _, z := range specialValues {
+				quat1 := FromEulerAngles(y, x, z)
+				ry, rx, rz := quat1.ToEulerAngles()
+				quat2 := FromEulerAngles(ry, rx, rz)
+				// quat must be equivalent
+				const e64 = 1e-14
+				cond1 := math.Abs(quat1[0]-quat2[0]) < e64 && math.Abs(quat1[1]-quat2[1]) < e64 && math.Abs(quat1[2]-quat2[2]) < e64 && math.Abs(quat1[3]-quat2[3]) < e64
+				cond2 := math.Abs(quat1[0]+quat2[0]) < e64 && math.Abs(quat1[1]+quat2[1]) < e64 && math.Abs(quat1[2]+quat2[2]) < e64 && math.Abs(quat1[3]+quat2[3]) < e64
+				if !cond1 && !cond2 {
+					fmt.Printf("test case %v, %v, %v failed\n", x, y, z)
+					fmt.Printf("result is %v, %v, %v\n", rx, ry, rz)
+					fmt.Printf("quat1 is %v\n", quat1)
+					fmt.Printf("quat2 is %v\n", quat2)
+					t.Fail()
+				}
+			}
+		}
+	}
+}

float64/vec2/vec2.go

@@ -104,6 +104,11 @@ func (vec *T) PracticallyEquals(compareVector *T, allowedDelta float64) bool {
 		(math.Abs(vec[1]-compareVector[1]) <= allowedDelta)
 }
 
+// PracticallyEquals compares two values if they are equal with each other within a delta tolerance.
+func PracticallyEquals(v1, v2, allowedDelta float64) bool {
+	return math.Abs(v1-v2) <= allowedDelta
+}
+
 // Invert inverts the vector.
 func (vec *T) Invert() *T {
 	vec[0] = -vec[0]
@@ -134,7 +139,7 @@ func (vec *T) Normalize() *T {
 	if sl == 0 || sl == 1 {
 		return vec
 	}
-	return vec.Scale(1 / math.Sqrt(sl))
+	return vec.Scale(1.0 / math.Sqrt(sl))
 }
 
 // Normalized returns a unit length normalized copy of the vector.
@@ -144,16 +149,18 @@ func (vec *T) Normalized() T {
 	return v
 }
 
-// Normal returns an orthogonal vector.
+// Normal returns a new normalized orthogonal vector.
 // The normal is orthogonal clockwise to the vector.
+// See also function Rotate90DegRight.
 func (vec *T) Normal() T {
 	n := *vec
 	n[0], n[1] = n[1], -n[0]
 	return *n.Normalize()
 }
 
-// Normal returns an orthogonal vector.
-// The normal is orthogonal counter clockwise to the vector.
+// NormalCCW returns a new normalized orthogonal vector.
+// The normal is orthogonal counterclockwise to the vector.
+// See also function Rotate90DegLeft.
 func (vec *T) NormalCCW() T {
 	n := *vec
 	n[0], n[1] = -n[1], n[0]
@@ -218,21 +225,46 @@ func (vec *T) RotateAroundPoint(point *T, angle float64) *T {
 }
 
 // Rotate90DegLeft rotates the vector 90 degrees left (counter-clockwise).
+// See also function NormalCCW.
 func (vec *T) Rotate90DegLeft() *T {
-	temp := vec[0]
-	vec[0] = -vec[1]
-	vec[1] = temp
+	vec[0], vec[1] = -vec[1], vec[0]
 	return vec
 }
 
 // Rotate90DegRight rotates the vector 90 degrees right (clockwise).
+// See also function Normal.
 func (vec *T) Rotate90DegRight() *T {
-	temp := vec[0]
-	vec[0] = vec[1]
-	vec[1] = -temp
+	vec[0], vec[1] = vec[1], -vec[0]
 	return vec
 }
 
+// Sinus returns the sinus value of the (shortest/smallest) angle between the two vectors a and b.
+// The returned sine value is in the range -1.0 ≤ value ≤ 1.0.
+// The angle is always considered to be in the range 0 to Pi radians and thus the sine value returned is always positive.
+func Sinus(a, b *T) float64 {
+	v := Cross(a, b) / math.Sqrt(a.LengthSqr()*b.LengthSqr())
+
+	if v > 1.0 {
+		return 1.0
+	} else if v < -1.0 {
+		return -1.0
+	}
+	return v
+}
+
+// Cosine returns the cosine value of the angle between the two vectors.
+// The returned cosine value is in the range -1.0 ≤ value ≤ 1.0.
+func Cosine(a, b *T) float64 {
+	v := Dot(a, b) / math.Sqrt(a.LengthSqr()*b.LengthSqr())
+
+	if v > 1.0 {
+		return 1.0
+	} else if v < -1.0 {
+		return -1.0
+	}
+	return v
+}
+
 // Angle returns the counter-clockwise angle of the vector from the x axis.
 func (vec *T) Angle() float64 {
 	return math.Atan2(vec[1], vec[0])
@@ -258,36 +290,30 @@ func Dot(a, b *T) float64 {
 	return a[0]*b[0] + a[1]*b[1]
 }
 
-// Cross returns the cross product of two vectors.
-func Cross(a, b *T) T {
-	return T{
-		a[1]*b[0] - a[0]*b[1],
-		a[0]*b[1] - a[1]*b[0],
-	}
+// Cross returns the "cross product" of two vectors.
+// In 2D space it is a scalar value.
+// It is the same as the determinant value of the 2D matrix constructed by the two vectors.
+// Cross product in 2D is not well-defined but this is the implementation stated at https://mathworld.wolfram.com/CrossProduct.html .
+func Cross(a, b *T) float64 {
+	return a[0]*b[1] - a[1]*b[0]
 }
 
-// Angle returns the angle between two vectors.
+// Angle returns the angle value of the (shortest/smallest) angle between the two vectors a and b.
+// The returned value is in the range 0 ≤ angle ≤ Pi radians.
 func Angle(a, b *T) float64 {
-	v := Dot(a, b) / (a.Length() * b.Length())
-	// prevent NaN
-	if v > 1. {
-		v = v - 2
-	} else if v < -1. {
-		v = v + 2
-	}
-	return math.Acos(v)
+	return math.Acos(Cosine(a, b))
 }
 
 // IsLeftWinding returns if the angle from a to b is left winding.
+// Two parallell or anti parallell vectors will give a false result.
 func IsLeftWinding(a, b *T) bool {
-	ab := b.Rotated(-a.Angle())
-	return ab.Angle() > 0
+	return Cross(a, b) > 0 // It's really the sign changing part of the Sinus(a, b) function
 }
 
 // IsRightWinding returns if the angle from a to b is right winding.
+// Two parallell or anti parallell vectors will give a false result.
 func IsRightWinding(a, b *T) bool {
-	ab := b.Rotated(-a.Angle())
-	return ab.Angle() < 0
+	return Cross(a, b) < 0 // It's really the sign changing part of the Sinus(a, b) function
 }
 
 // Min returns the component wise minimum of two vectors.
@@ -316,7 +342,7 @@ func Max(a, b *T) T {
 
 // Interpolate interpolates between a and b at t (0,1).
 func Interpolate(a, b *T, t float64) T {
-	t1 := 1 - t
+	t1 := 1.0 - t
 	return T{
 		a[0]*t1 + b[0]*t,
 		a[1]*t1 + b[1]*t,

float64/vec2/vec2_test.go

@@ -2,6 +2,7 @@ package vec2
 
 import (
 	"math"
+	"strconv"
 	"testing"
 )
 
@@ -214,3 +215,126 @@ func TestMuled(t *testing.T) {
 		t.Fail()
 	}
 }
+
+func TestAngle(t *testing.T) {
+	radFor45deg := math.Pi / 4.0
+	testSetups := []struct {
+		a, b          T
+		expectedAngle float64
+		name          string
+	}{
+		{a: T{1, 0}, b: T{1, 0}, expectedAngle: 0 * radFor45deg, name: "0/360 degree angle, equal/parallell vectors"},
+		{a: T{1, 0}, b: T{1, 1}, expectedAngle: 1 * radFor45deg, name: "45 degree angle"},
+		{a: T{1, 0}, b: T{0, 1}, expectedAngle: 2 * radFor45deg, name: "90 degree angle, orthogonal vectors"},
+		{a: T{1, 0}, b: T{-1, 1}, expectedAngle: 3 * radFor45deg, name: "135 degree angle"},
+		{a: T{1, 0}, b: T{-1, 0}, expectedAngle: 4 * radFor45deg, name: "180 degree angle, inverted/anti parallell vectors"},
+		{a: T{1, 0}, b: T{-1, -1}, expectedAngle: (8 - 5) * radFor45deg, name: "225 degree angle"},
+		{a: T{1, 0}, b: T{0, -1}, expectedAngle: (8 - 6) * radFor45deg, name: "270 degree angle, orthogonal vectors"},
+		{a: T{1, 0}, b: T{1, -1}, expectedAngle: (8 - 7) * radFor45deg, name: "315 degree angle"},
+	}
+
+	for _, testSetup := range testSetups {
+		t.Run(testSetup.name, func(t *testing.T) {
+			angle := Angle(&testSetup.a, &testSetup.b)
+
+			if !PracticallyEquals(angle, testSetup.expectedAngle, 0.00000001) {
+				t.Errorf("Angle expected to be %f but was %f for test \"%s\".", testSetup.expectedAngle, angle, testSetup.name)
+			}
+		})
+	}
+}
+
+func TestCosine(t *testing.T) {
+	radFor45deg := math.Pi / 4.0
+	testSetups := []struct {
+		a, b           T
+		expectedCosine float64
+		name           string
+	}{
+		{a: T{1, 0}, b: T{1, 0}, expectedCosine: math.Cos(0 * radFor45deg), name: "0/360 degree angle, equal/parallell vectors"},
+		{a: T{1, 0}, b: T{1, 1}, expectedCosine: math.Cos(1 * radFor45deg), name: "45 degree angle"},
+		{a: T{1, 0}, b: T{0, 1}, expectedCosine: math.Cos(2 * radFor45deg), name: "90 degree angle, orthogonal vectors"},
+		{a: T{1, 0}, b: T{-1, 1}, expectedCosine: math.Cos(3 * radFor45deg), name: "135 degree angle"},
+		{a: T{1, 0}, b: T{-1, 0}, expectedCosine: math.Cos(4 * radFor45deg), name: "180 degree angle, inverted/anti parallell vectors"},
+		{a: T{1, 0}, b: T{-1, -1}, expectedCosine: math.Cos(5 * radFor45deg), name: "225 degree angle"},
+		{a: T{1, 0}, b: T{0, -1}, expectedCosine: math.Cos(6 * radFor45deg), name: "270 degree angle, orthogonal vectors"},
+		{a: T{1, 0}, b: T{1, -1}, expectedCosine: math.Cos(7 * radFor45deg), name: "315 degree angle"},
+	}
+
+	for _, testSetup := range testSetups {
+		t.Run(testSetup.name, func(t *testing.T) {
+			cos := Cosine(&testSetup.a, &testSetup.b)
+
+			if !PracticallyEquals(cos, testSetup.expectedCosine, 0.00000001) {
+				t.Errorf("Cosine expected to be %f but was %f for test \"%s\".", testSetup.expectedCosine, cos, testSetup.name)
+			}
+		})
+	}
+}
+
+func TestSinus(t *testing.T) {
+	radFor45deg := math.Pi / 4.0
+	testSetups := []struct {
+		a, b         T
+		expectedSine float64
+		name         string
+	}{
+		{a: T{1, 0}, b: T{1, 0}, expectedSine: math.Sin(0 * radFor45deg), name: "0/360 degree angle, equal/parallell vectors"},
+		{a: T{1, 0}, b: T{1, 1}, expectedSine: math.Sin(1 * radFor45deg), name: "45 degree angle"},
+		{a: T{1, 0}, b: T{0, 1}, expectedSine: math.Sin(2 * radFor45deg), name: "90 degree angle, orthogonal vectors"},
+		{a: T{1, 0}, b: T{-1, 1}, expectedSine: math.Sin(3 * radFor45deg), name: "135 degree angle"},
+		{a: T{1, 0}, b: T{-1, 0}, expectedSine: math.Sin(4 * radFor45deg), name: "180 degree angle, inverted/anti parallell vectors"},
+		{a: T{1, 0}, b: T{-1, -1}, expectedSine: math.Sin(5 * radFor45deg), name: "225 degree angle"},
+		{a: T{1, 0}, b: T{0, -1}, expectedSine: math.Sin(6 * radFor45deg), name: "270 degree angle, orthogonal vectors"},
+		{a: T{1, 0}, b: T{1, -1}, expectedSine: math.Sin(7 * radFor45deg), name: "315 degree angle"},
+	}
+
+	for _, testSetup := range testSetups {
+		t.Run(testSetup.name, func(t *testing.T) {
+			sin := Sinus(&testSetup.a, &testSetup.b)
+
+			if !PracticallyEquals(sin, testSetup.expectedSine, 0.00000001) {
+				t.Errorf("Sine expected to be %f but was %f for test \"%s\".", testSetup.expectedSine, sin, testSetup.name)
+			}
+		})
+	}
+}
+
+func TestLeftRightWinding(t *testing.T) {
+	a := T{1.0, 0.0}
+
+	for angle := 0; angle <= 360; angle += 15 {
+		rad := (math.Pi / 180.0) * float64(angle)
+
+		bx := clampDecimals(math.Cos(rad), 4)
+		by := clampDecimals(math.Sin(rad), 4)
+		b := T{bx, by}
+
+		t.Run("left winding angle "+strconv.Itoa(angle), func(t *testing.T) {
+			lw := IsLeftWinding(&a, &b)
+			rw := IsRightWinding(&a, &b)
+
+			if angle%180 == 0 {
+				// No winding at 0, 180 and 360 degrees
+				if lw || rw {
+					t.Errorf("Neither left or right winding should be true on angle %d. Left winding=%t, right winding=%t", angle, lw, rw)
+				}
+			} else if angle < 180 {
+				// Left winding at 0 < angle < 180
+				if !lw || rw {
+					t.Errorf("Left winding should be true (not right winding) on angle %d. Left winding=%t, right winding=%t", angle, lw, rw)
+				}
+			} else if angle > 180 {
+				// Right winding at 180 < angle < 360
+				if lw || !rw {
+					t.Errorf("Right winding should be true (not left winding) on angle %d. Left winding=%t, right winding=%t", angle, lw, rw)
+				}
+			}
+		})
+	}
+}
+
+func clampDecimals(decimalValue float64, amountDecimals float64) float64 {
+	factor := math.Pow(10, amountDecimals)
+	return math.Round(decimalValue*factor) / factor
+}