fix bug for semi-arc polygon generation not completing full arc

Author: soypat

md2/polygon.go

@@ -225,7 +225,7 @@ func arcCenterFrom2points(p1, p2 Vec, r float64) (Vec, float64, error) {
 	} else if maxChordLen-chordLen <= -semiArcTol {
 		return Vec{}, 0, errBadArc
 	} else if ms1.EqualWithinAbs(sinTheta, 1, semiArcTol) {
-		return Scale(0.5, Add(p1, p2)), math.Copysign(math.Pi/2, r), nil
+		return Scale(0.5, Add(p1, p2)), math.Copysign(math.Pi, r), nil
 	}
 	chordThetaDiv2 := math.Asin(sinTheta)
 	diffTo90 := chordThetaDiv2 - math.Pi/2

md2/polygon_test.go

@@ -5,9 +5,12 @@
 package md2
 
 import (
+	"strconv"
 	"testing"
 
 	math "math"
+	"github.com/soypat/geometry/internal"
+	ms1 "github.com/soypat/geometry/md1"
 )
 
 var testoffsets = []Vec{{-1, -2}, {-2, 1}, {2, -1}, {}, {1, 0}, {0, 1}, {1, 1}}
@@ -187,3 +190,44 @@ func TestArc_invalidArc(t *testing.T) {
 		}
 	}
 }
+
+func TestArcUniformity(t *testing.T) {
+	const tol = internal.Smallfloat64
+	var cases = []struct {
+		start, end Vec
+		radius     float64
+	}{
+		{start: Vec{0, 0}, end: Vec{0, 4}, radius: 2},
+	}
+	var poly PolygonBuilder
+	var buf []Vec
+	var err error
+	for itest, test := range cases {
+		for facets := 3; facets <= 3; facets++ {
+			t.Run("case="+strconv.Itoa(itest)+" facets="+strconv.Itoa(facets), func(t *testing.T) {
+				poly.Reset()
+				poly.Add(test.start)
+				poly.Add(test.end).Arc(test.radius, facets)
+				buf, err = poly.AppendVecs(buf[:0])
+				if err != nil {
+					t.Fatal(err)
+				}
+				center, angle, err := arcCenterFrom2points(test.start, test.end, test.radius)
+				if err != nil {
+					t.Errorf("getting center: %s", err)
+				}
+				anglePerFacet := angle / float64(facets)
+				dir0 := Sub(test.start, center)
+				for i, v := range buf {
+					dir := Sub(v, center)
+					cosangle := Cos(dir, dir0)
+					gotAngle := math.Acos(cosangle)
+					wantAngle := float64(i) * anglePerFacet
+					if !ms1.EqualWithinAbs(gotAngle, wantAngle, tol) {
+						t.Errorf("%d/%d facet-dir want angle swept %f, got %f", i+1, len(buf), wantAngle, gotAngle)
+					}
+				}
+			})
+		}
+	}
+}

md3/md3_test.go

@@ -6,10 +6,11 @@ package md3
 
 import (
 	"fmt"
-	"math"
 	"math/rand"
 	"testing"
 
+	math "math"
+
 	"github.com/soypat/geometry/internal"
 )
 
@@ -36,14 +37,16 @@ func TestRotationConversions(t *testing.T) {
 		Rotation(2*math.Pi, Vec{X: 1}),
 		Rotation(0, Vec{X: 1}),
 	}
+	maxAngle := 5 * math.Pi / 180
+	angleTol := math.Cos(maxAngle / 2)
 	rng := newRNG(1)
 	for _, rot := range rotations {
 		angle, axis := rot.Rotation()
 		gotrot := Rotation(angle, axis)
 		if !EqualQuat(rot, gotrot, tol) {
 			t.Errorf("want %v, got %v", rot, gotrot)
 		}
-		if !rot.EqualOrientation(gotrot, tol) {
+		if !rot.EqualOrientation(gotrot, angleTol) {
 			t.Errorf("not equal orientations %v, got %v", rot, gotrot)
 		}
 		m3 := rot.RotationMat3()

md3/rotation.go

@@ -124,13 +124,14 @@ func RotationBetweenVecs(start, dest Vec) Quat {
 }
 
 // EqualOrientation returns whether the quaternions represents the same orientation with a given tolerence
-func (q1 Quat) EqualOrientation(q2 Quat, tolsq float64) bool {
-	n1sq, n2sq := q1.Dot(q1), q2.Dot(q2)
+func (q1 Quat) EqualOrientation(q2 Quat, tol float64) bool {
+	n1sq := q1.Dot(q1)
+	n2sq := q2.Dot(q2)
 	if n1sq == 0 || n2sq == 0 {
 		return false // Degenerate quaternion.
 	}
-	dot := math.Abs(q1.Dot(q2) / (n1sq * n2sq))
-	return dot > 1-tolsq
+	d := q1.Dot(q2)
+	return d*d >= tol*tol*(n1sq*n2sq)
 }
 
 /*

ms2/polygon.go

@@ -221,7 +221,7 @@ func arcCenterFrom2points(p1, p2 Vec, r float32) (Vec, float32, error) {
 	} else if maxChordLen-chordLen <= -semiArcTol {
 		return Vec{}, 0, errBadArc
 	} else if ms1.EqualWithinAbs(sinTheta, 1, semiArcTol) {
-		return Scale(0.5, Add(p1, p2)), math.Copysign(math.Pi/2, r), nil
+		return Scale(0.5, Add(p1, p2)), math.Copysign(math.Pi, r), nil
 	}
 	chordThetaDiv2 := math.Asin(sinTheta)
 	diffTo90 := chordThetaDiv2 - math.Pi/2

ms2/polygon_test.go

@@ -1,9 +1,12 @@
 package ms2
 
 import (
+	"strconv"
 	"testing"
 
 	math "github.com/chewxy/math32"
+	"github.com/soypat/geometry/internal"
+	"github.com/soypat/geometry/ms1"
 )
 
 var testoffsets = []Vec{{-1, -2}, {-2, 1}, {2, -1}, {}, {1, 0}, {0, 1}, {1, 1}}
@@ -183,3 +186,44 @@ func TestArc_invalidArc(t *testing.T) {
 		}
 	}
 }
+
+func TestArcUniformity(t *testing.T) {
+	const tol = internal.Smallfloat32
+	var cases = []struct {
+		start, end Vec
+		radius     float32
+	}{
+		{start: Vec{0, 0}, end: Vec{0, 4}, radius: 2},
+	}
+	var poly PolygonBuilder
+	var buf []Vec
+	var err error
+	for itest, test := range cases {
+		for _, facets := range []int{2, 3, 5, 9} {
+			t.Run("case="+strconv.Itoa(itest)+" facets="+strconv.Itoa(facets), func(t *testing.T) {
+				poly.Reset()
+				poly.Add(test.start)
+				poly.Add(test.end).Arc(test.radius, facets)
+				buf, err = poly.AppendVecs(buf[:0])
+				if err != nil {
+					t.Fatal(err)
+				}
+				center, angle, err := arcCenterFrom2points(test.start, test.end, test.radius)
+				if err != nil {
+					t.Errorf("getting center: %s", err)
+				}
+				anglePerFacet := angle / float32(facets)
+				dir0 := Sub(test.start, center)
+				for i, v := range buf {
+					dir := Sub(v, center)
+					cosangle := Cos(dir, dir0)
+					gotAngle := math.Acos(cosangle)
+					wantAngle := float32(i) * anglePerFacet
+					if !ms1.EqualWithinAbs(gotAngle, wantAngle, tol) {
+						t.Errorf("%d/%d facet-dir want angle swept %f, got %f", i+1, len(buf), wantAngle, gotAngle)
+					}
+				}
+			})
+		}
+	}
+}