add go-play3d logic

Author: soypat

internal/internalmath.go

@@ -1,6 +1,8 @@
 package internal
 
 const (
-	Smallfloat32 = 1e-5
-	Smallfloat64 = 1e-8
+	Smallfloat32 float32 = 1e-5
+	Smallfloat64 float64 = 1e-8
+	Largefloat32 float32 = 1e23
+	Largefloat64 float64 = 1e53
 )

ms2/line.go

@@ -0,0 +1,66 @@
+package ms2
+
+import (
+	math "github.com/chewxy/math32"
+)
+
+// Line can be used to represent a line between two points
+// or an infinite line.
+type Line [2]Vec
+
+// Interpolate takes a value between 0 and 1 to linearly
+// interpolate a point on the line.
+//
+//	Interpolate(0) returns l[0]
+//	Interpolate(1) returns l[1]
+func (ln Line) Interpolate(t float32) Vec {
+	lineDir := Sub(ln[1], ln[0])
+	return Add(ln[0], Scale(t, lineDir))
+}
+
+// DistanceInfinite returns the minimum euclidean distance of point p to the infinite line represented by l.
+func (ln Line) DistanceInfinite(point Vec) float32 {
+	// https://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html
+	p1 := ln[0]
+	p2 := ln[1]
+	num := math.Abs((p2.X-p1.X)*(p1.Y-point.Y) - (p1.X-point.X)*(p2.Y-p1.Y))
+	return num / math.Hypot(p2.X-p1.X, p2.Y-p1.Y)
+}
+
+// ClosestInfinite returns the point on the infinite line closest to the argument point.
+func (ln Line) ClosestInfinite(point Vec) Vec {
+	// https://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html
+	t := -Dot(Sub(ln[0], point), Sub(ln[1], point)) / Norm2(Sub(ln[1], ln[0]))
+	return ln.Interpolate(t)
+}
+
+// Closest returns the closest point on the line l[0]..l[1] to the argument point.
+// The integer return param is 0 or 1 when closest to vertex l[0] or l[1], respectively,
+// and is -1 for when the point is closest to the line segment.
+func (ln Line) Closest(point Vec) (closest Vec, vertexOrSegment int8) {
+	lineDir := Sub(ln[1], ln[0])
+	perpendicular := Vec{-lineDir.Y, lineDir.X}
+	perpend2 := Add(ln[1], perpendicular)
+	e2 := Line{ln[1], perpend2}.edgeEquation(point)
+	if e2 > 0 {
+		return ln[1], 0
+	}
+	perpend1 := Add(ln[0], perpendicular)
+	e1 := Line{ln[0], perpend1}.edgeEquation(point)
+	if e1 < 0 {
+		return ln[0], 1
+	}
+	e3 := ln.DistanceInfinite(point)
+	toPoint := Scale(-e3, Unit(perpendicular))
+	return Sub(point, toPoint), -1
+}
+
+// edgeEquation returns a signed distance of a point to
+// an infinite line defined by two points
+// Edge equation for a line passing through (X,Y)
+// with gradient dY/dX
+// E ( x; y ) =(x-X)*dY - (y-Y)*dX
+func (ln Line) edgeEquation(p Vec) float32 {
+	dxy := Sub(ln[1], ln[0])
+	return (p.X-ln[0].X)*dxy.Y - (p.Y-ln[0].Y)*dxy.X
+}

ms2/triangle.go

@@ -2,6 +2,7 @@ package ms2
 
 import (
 	math "github.com/chewxy/math32"
+	"github.com/soypat/geometry/internal"
 )
 
 // Triangle represents a triangle in 2D space and
@@ -17,6 +18,17 @@ func (t Triangle) Centroid() Vec {
 	return Scale(1.0/3.0, Add(Add(t[0], t[1]), t[2]))
 }
 
+// Sides returns the triangle's sides as lines:
+//
+//	[(t[0],t[1]), (t[1],t[2]), (t[2],t[0])]
+func (t Triangle) Sides() [3]Line {
+	return [3]Line{
+		{t[0], t[1]},
+		{t[1], t[2]},
+		{t[2], t[0]},
+	}
+}
+
 // Area returns the surface area of the triangle.
 func (t Triangle) Area() float32 {
 	// Heron's Formula, see https://en.wikipedia.org/wiki/Heron%27s_formula.
@@ -64,33 +76,10 @@ func (t Triangle) IsDegenerate(tol float32) bool {
 	}
 	// calculate vertex distance from longest side
 	ln := Line{t[longIdx], t[(longIdx+1)%3]}
-	dist := ln.Distance(t[(longIdx+2)%3])
+	dist := ln.DistanceInfinite(t[(longIdx+2)%3])
 	return dist <= tol
 }
 
-// Line is an infinite 3D Line
-// defined by two points on the Line.
-type Line [2]Vec
-
-// Interpolate takes a value between 0 and 1 to linearly
-// interpolate a point on the line.
-//
-//	Interpolate(0) returns l[0]
-//	Interpolate(1) returns l[1]
-func (l Line) Interpolate(t float32) Vec {
-	lineDir := Sub(l[1], l[0])
-	return Add(l[0], Scale(t, lineDir))
-}
-
-// Distance returns the minimum euclidean Distance of point p to the line.
-func (l Line) Distance(p Vec) float32 {
-	// https://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html
-	p1 := l[0]
-	p2 := l[1]
-	num := math.Abs((p2.X-p1.X)*(p1.Y-p.Y) - (p1.X-p.X)*(p2.Y-p1.Y))
-	return num / math.Hypot(p2.X-p1.X, p2.Y-p1.Y)
-}
-
 // sort performs the sort-3 algorithm and returns
 // l1, l2, l3 such that l1 ≤ l2 ≤ l3.
 func sort(a, b, c float32) (l1, l2, l3 float32) {
@@ -110,14 +99,57 @@ func sort(a, b, c float32) (l1, l2, l3 float32) {
 // orderedLengths returns the lengths of the sides of the triangle such that
 // a ≤ b ≤ c.
 func (t Triangle) orderedLengths() (a, b, c float32) {
-	s1, s2, s3 := t.sides()
+	s1, s2, s3 := t.edges()
 	l1 := Norm(s1)
 	l2 := Norm(s2)
 	l3 := Norm(s3)
 	return sort(l1, l2, l3)
 }
 
-// sides returns vectors for each of the sides of t.
-func (t Triangle) sides() (Vec, Vec, Vec) {
+// edges returns vectors for each of the edges of t.
+func (t Triangle) edges() (Vec, Vec, Vec) {
 	return Sub(t[1], t[0]), Sub(t[2], t[1]), Sub(t[0], t[2])
 }
+
+// Contains returns true if point is contained within the triangle's surface.
+func (t Triangle) Contains(point Vec) bool {
+	d1 := d2Sign(point, t[0], t[1])
+	d2 := d2Sign(point, t[1], t[2])
+	d3 := d2Sign(point, t[2], t[0])
+	has_neg := (d1 < 0) || (d2 < 0) || (d3 < 0)
+	has_pos := (d1 > 0) || (d2 > 0) || (d3 > 0)
+	return !(has_neg && has_pos)
+}
+
+func d2Sign(p1, p2, p3 Vec) float32 {
+	// TODO: replace with CopyOrientation?
+	return (p1.X-p3.X)*(p2.Y-p3.Y) - (p2.X-p3.X)*(p1.Y-p3.Y)
+}
+
+// Closest returns the point on the triangle closest to the argument point p.
+// side and vertex are non negative to flag the point is closest the non-negative side or vertex index.
+// If the point lies on the triangle it returns the same point and side=vertex=-1. side and vertex cannot be both non-negative.
+func (t Triangle) Closest(p Vec) (closest Vec, side int8, vertex int8) {
+	if t.Contains(p) {
+		return p, -1, -1
+	}
+	minDist := internal.Largefloat32
+	for j := range t {
+		nxt := (j + 1) % 3
+		edge := Line{{X: t[j].X, Y: t[j].Y}, {X: t[nxt].X, Y: t[nxt].Y}}
+		pointOnTriangle, maybeVertex := edge.Closest(p)
+		d2 := Norm2(Sub(p, pointOnTriangle))
+		if d2 < minDist {
+			if vertex < 0 {
+				vertex = -1
+				side = int8(j)
+			} else {
+				side = -1
+				vertex = (int8(j) + maybeVertex) % 3
+			}
+			minDist = d2
+			closest = pointOnTriangle
+		}
+	}
+	return closest, side, vertex
+}

ms2/vec.go

@@ -236,6 +236,7 @@ func CopyOrientation(f float32, p1, p2, p3 Vec) float32 {
 	if slope1 == slope2 {
 		return 0
 	}
+
 	return math.Copysign(f, slope2-slope1)
 }
 

ms3/line.go

@@ -0,0 +1,30 @@
+package ms3
+
+// Line can be used to represent a line between two points
+// or an infinite line.
+type Line [2]Vec
+
+// Interpolate takes a value between 0 and 1 to linearly
+// interpolate a point on the line.
+//
+//	Interpolate(0) returns l[0]
+//	Interpolate(1) returns l[1]
+func (ln Line) Interpolate(t float32) Vec {
+	lineDir := Sub(ln[1], ln[0])
+	return Add(ln[0], Scale(t, lineDir))
+}
+
+// DistanceInfinite returns the minimum euclidean Distance of point p
+// to the infinite line represented by l.
+func (ln Line) DistanceInfinite(point Vec) float32 {
+	// https://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html
+	num := Norm(Cross(Sub(point, ln[0]), Sub(point, ln[1])))
+	return num / Norm(Sub(ln[1], ln[0]))
+}
+
+// ClosestInfinite returns the point on the infinite line closest to the argument point.
+func (ln Line) ClosestInfinite(point Vec) Vec {
+	// https://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html
+	t := -Dot(Sub(ln[0], point), Sub(ln[1], point)) / Norm2(Sub(ln[1], ln[0]))
+	return ln.Interpolate(t)
+}

ms3/plane.go

@@ -0,0 +1,20 @@
+package ms3
+
+import math "github.com/chewxy/math32"
+
+// Plane represents a plane in 3D space. To safely construct one use [Triangle.Plane]
+type Plane struct {
+	// p is a point on the plane
+	p Vec
+	// n is the unit vector normal to the plane.
+	n Vec
+}
+
+func newPlane(p, n Vec) Plane {
+	return Plane{p: p, n: Unit(n)}
+}
+
+func (p Plane) distanceTo(q Vec) float32 {
+	v := Sub(q, p.p)
+	return math.Abs(Dot(v, p.n))
+}

ms3/quat.go

@@ -53,7 +53,8 @@ const (
 type Quat struct {
 	// V contains I, J and K imaginary parts.
 	I, J, K float32
-	W       float32
+	// W is the quaternion's real part.
+	W float32
 }
 
 // IJK returns I,J,K fields of q as a vector with set fields X,Y,Z, respectively.

ms3/tetrahedron.go

@@ -0,0 +1,79 @@
+package ms3
+
+import (
+	math "github.com/chewxy/math32"
+)
+
+// Tetra represents a tetrahedron in 3D space.
+type Tetra [4]Vec
+
+// Centroid returns the tetrahedron's centroid.
+func (t Tetra) Centroid() Vec {
+	return Scale(1./4., Add(t[0], Add(t[1], Add(t[2], t[3]))))
+}
+
+// Sides returns the sides as lines of the tetrahedron:
+//
+//	[(t[0],t[1]),(t[1],t[2]),(t[2],t[0]),
+//	(t[3],t[2]),(t[0],t[3]),(t[1],t[3])]
+func (t Tetra) Sides() [6]Line {
+	return [6]Line{
+		{t[0], t[1]}, {t[1], t[2]}, {t[2], t[0]},
+		{t[3], t[2]}, {t[0], t[3]}, {t[1], t[3]},
+	}
+}
+
+// Edges returns the directional vectors composing the sides of the tetrahedron as returned by [Tetra.Sides].
+func (t Tetra) Edges() [6]Vec {
+	return [6]Vec{
+		Sub(t[1], t[0]), Sub(t[2], t[1]), Sub(t[0], t[2]),
+		Sub(t[2], t[3]), Sub(t[3], t[0]), Sub(t[3], t[1]),
+	}
+}
+
+// Volume returns the volume contained in the tetrahedron.
+func (t Tetra) Volume() float32 {
+	const third = 1.0 / 3.0
+	base := Triangle{t[0], t[1], t[2]}
+	area := base.Area()
+	height := base.Plane().distanceTo(t[3])
+	return third * area * height
+}
+
+// Aspect returns the aspect ratio of the tetrahedron calculated
+// as: longestEdge/minHeight.
+func (t Tetra) Aspect() float32 {
+	e := t.longestEdge()
+	heights := t.heights()
+	hmin, _, _ := Sort(heights[0], heights[1], heights[2])
+	if heights[3] < hmin {
+		hmin = heights[3]
+	}
+	return e / hmin
+
+}
+
+func (t Tetra) longestEdge() float32 {
+	edges := t.Edges()
+	_, _, e1 := Sort(Norm2(edges[0]), Norm2(edges[1]), Norm2(edges[2]))
+	_, _, e2 := Sort(Norm2(edges[3]), Norm2(edges[4]), Norm2(edges[5]))
+	if e1 > e2 {
+		return math.Sqrt(e1)
+	}
+	return math.Sqrt(e2)
+}
+
+// heights returns the "heights" of the tetrahedron,
+// which are basically the shortest normal dropped from a vertex to the opposite face.
+func (t Tetra) heights() (alt [4]float32) {
+	for i := range t {
+		j := (i + 1) % 4
+		k := (i + 2) % 4
+		l := (i + 3) % 4
+		e1 := Sub(t[k], t[j])
+		e2 := Sub(t[l], t[j])
+		p := newPlane(t[l], Cross(e1, e2))
+		alt[i] = p.distanceTo(t[i])
+	}
+	return alt
+}