<commit_before>/*package geom contains routines for computing geometric quantities.

Contains implementations of algorithms described in Platis & Theoharis, 2015
as well as Schneider & Eberly.

The calling signatures might be more convoluted than they have to be because
I was too worried about memory consumption when I started it. I should include
examples.
*/
package geom

import (
	"math"
)

// Vec is a three dimensional vector. (Duh!)
type Vec [3]float32

// PluckerVec represents a ray. If P and L are the position of the ray's 
// origin and the unit vector representing its direction, respectively, then
// U = L and V = L cross P.
type PluckerVec struct {
	U, V Vec
}

// AnchoredPluckerVec is a Plucker vector which also saves the position of
// the ray's origin.
type AnchoredPluckerVec struct {
	PluckerVec
	P Vec
}

// Init initializes a Plucker vector given a ray origin, P, and a unit
// direction vector, L.
func (p *PluckerVec) Init(P, L *Vec) {
	p.U = *L
	
	p.V[0] = -P[1]*L[2] + P[2]*L[1]
    p.V[1] = -P[2]*L[0] + P[0]*L[2]
    p.V[2] = -P[0]*L[1] + P[1]*L[0]
}

// InitFromSegment initialized a Plucker vector which corresponds to a ray
// pointing from the position vector P1 to the position vector P2.
func (p *PluckerVec) InitFromSegment(P1, P2 *Vec) {
	var sum float32
	for i := 0; i < 3; i++ {
		p.U[i] = P2[i] - P1[i]
		sum += p.U[i]*p.U[i]
	}
	sum = float32(math.Sqrt(float64(sum)))
	//for i := 0; i < 3; i++ { p.U[i] /= sum }

	p.V[0] = -P1[1]*p.U[2] + P1[2]*p.U[1]
    p.V[1] = -P1[2]*p.U[0] + P1[0]*p.U[2]
    p.V[2] = -P1[0]*p.U[1] + P1[1]*p.U[0]
}

// Translate translates a Plucker vector along the given vector.
func (p *PluckerVec) Translate(dx *Vec) {
	p.V[0] += -dx[1]*p.U[2] + dx[2]*p.U[1]
    p.V[1] += -dx[2]*p.U[0] + dx[0]*p.U[2]
    p.V[2] += -dx[0]*p.U[1] + dx[1]*p.U[0]
}

// Dot computes the permuted inner product of p1 and p2, i.e.
// p1.U*p2.V + p1.V*p2.U.
func (p1 *PluckerVec) Dot(p2 *PluckerVec, flip bool) float32 {
	var sum float32
	for i := 0; i < 3; i++ {
		sum += p1.U[i]*p2.V[i] + p1.V[i]*p2.U[i]
	}
	if flip {
		return sum
	} else {
		return -sum
	}
}

// Dot computes the permuted inner product of p1 and p2, i.e.
// p1.U*p2.V + p1.V*p2.U and also returns a sign flag of -1, 0, or +1 if
// that product is negative, zero, or positive, respectively.
func (p1 *PluckerVec) SignDot(p2 *PluckerVec, flip bool) (float32, int) {
	dot := p1.Dot(p2, flip)
	if dot == 0 {
		return dot, 0
	} else if dot > 0 {
		return dot, +1
	} else {
		return dot, -1
	}
}

// Init initializes an anchored Plucker vector given a ray origin, P, and a
// unit direction vector, L.
func (ap *AnchoredPluckerVec) Init(P, L *Vec) {
	ap.PluckerVec.Init(P, L)
	ap.P = *P
}

// InitFromSegment initialized a Plucker vector which corresponds to a ray
// pointing from the position vector P1 to the position vector P2.
func (ap *AnchoredPluckerVec) InitFromSegment(P1, P2 *Vec) {
	ap.PluckerVec.InitFromSegment(P1, P2)
	ap.P = *P1
}

// Translate translates a Plucker vector along the given vector.
func (ap *AnchoredPluckerVec) Translate(dx *Vec) {
	ap.PluckerVec.Translate(dx)
	for i := 0; i < 3; i++ { ap.P[i] += dx[i] }
}

// Tetra is a tetrahedron. (Duh!)
//
// Face ordering is:
// F0(V3, V2, V1)
// F1(V2, V3, V0)
// F2(V1, V0, V3)
// F3(V0, V1, V2)
type Tetra [4]Vec

var tetraIdxs = [4][3]int {
	[3]int{ 3, 2, 1 },
	[3]int{ 2, 3, 0 },
	[3]int{ 1, 0, 3 },
	[3]int{ 0, 1, 2 },
}

// VertexIdx returns the index into the given tetrahedron corresponding to
// the specified face and vertex.
func (_ *Tetra) VertexIdx(face, vertex int) int {
	return tetraIdxs[face][vertex]
}

// Orient arranges tetrahedron points so that all faces point outward for
// dir = +1 and inward for dir = -1.
func (t *Tetra) Orient(dir int) {
	v, w, n := Vec{}, Vec{}, Vec{}
	for i := 0; i < 3; i++ {
		v[i] = t[1][i] - t[0][i]
		w[i] = t[2][i] - t[0][i]
	}
	n[0] = -v[1]*w[2] + v[2]*w[1]
    n[1] = -v[2]*w[0] + v[0]*w[2]
    n[2] = -v[0]*w[1] + v[1]*w[0]

	var dot float32
	for i := 0; i < 3; i++ {
		dot += n[i] * (t[3][i] - t[0][i])
	}

	if (dot < 0 && dir == -1) || (dot > 0 && dir == +1) {
		t[0], t[1] = t[1], t[0]
	}
}

// Translate translates a tetrahedron by the given vector.
func (t *Tetra) Translate(dx *Vec) {
	for i := 0; i < 4; i++ {
		for j := 0; j < 3; j++ {
			t[i][j] += dx[j]
		}
	}
}

// Sphere is exactly what you think it is.
type Sphere struct {
	X, Y, Z, R float32
}

// Intersect returns true if the two spheres intersect and false otherwise.
func (s1 *Sphere) Intersect(s2 *Sphere) bool {
	dx, dy, dz, dr := s1.X - s2.X, s1.Y - s2.Y, s1.Z - s2.Z, s1.R - s2.R
	return dr*dr > dx*dx + dy*dy + dz*dz
}

// BoundingSphere draws a bounding sphere aorund the given tetrahedron.
func (t *Tetra) BoundingSphere(sph *Sphere) {
	bx := (t[0][0] + t[1][0] + t[2][0] + t[3][0]) / 4
	by := (t[0][1] + t[1][1] + t[2][1] + t[3][1]) / 4
	bz := (t[0][2] + t[1][2] + t[2][2] + t[3][2]) / 4

	dx, dy, dz := bx-t[0][0], by-t[0][1], bz-t[0][2]
	maxRSqr := dx*dx + dy*dy + dz*dz
	for i := 1; i < 4; i++ {
		dx, dy, dz = bx-t[i][0], by-t[i][1], bz-t[i][2]
		rSqr := dx*dx + dy*dy + dz*dz
		if rSqr > maxRSqr { maxRSqr = rSqr }
	}

	sph.X, sph.Y, sph.Z = bx, by, bz
	sph.R = float32(math.Sqrt(float64(maxRSqr)))
}

// TetraFaceBary contains information specifying the barycentric coordinates<commit_msg>Fixed bug in how sphere intersection was detected.<commit_after>/*package geom contains routines for computing geometric quantities.

Contains implementations of algorithms described in Platis & Theoharis, 2015
as well as Schneider & Eberly.

The calling signatures might be more convoluted than they have to be because
I was too worried about memory consumption when I started it. I should include
examples.
*/
package geom

import (
	"math"
)

// Vec is a three dimensional vector. (Duh!)
type Vec [3]float32

// PluckerVec represents a ray. If P and L are the position of the ray's 
// origin and the unit vector representing its direction, respectively, then
// U = L and V = L cross P.
type PluckerVec struct {
	U, V Vec
}

// AnchoredPluckerVec is a Plucker vector which also saves the position of
// the ray's origin.
type AnchoredPluckerVec struct {
	PluckerVec
	P Vec
}

// Init initializes a Plucker vector given a ray origin, P, and a unit
// direction vector, L.
func (p *PluckerVec) Init(P, L *Vec) {
	p.U = *L
	
	p.V[0] = -P[1]*L[2] + P[2]*L[1]
    p.V[1] = -P[2]*L[0] + P[0]*L[2]
    p.V[2] = -P[0]*L[1] + P[1]*L[0]
}

// InitFromSegment initialized a Plucker vector which corresponds to a ray
// pointing from the position vector P1 to the position vector P2.
func (p *PluckerVec) InitFromSegment(P1, P2 *Vec) {
	var sum float32
	for i := 0; i < 3; i++ {
		p.U[i] = P2[i] - P1[i]
		sum += p.U[i]*p.U[i]
	}
	sum = float32(math.Sqrt(float64(sum)))
	//for i := 0; i < 3; i++ { p.U[i] /= sum }

	p.V[0] = -P1[1]*p.U[2] + P1[2]*p.U[1]
    p.V[1] = -P1[2]*p.U[0] + P1[0]*p.U[2]
    p.V[2] = -P1[0]*p.U[1] + P1[1]*p.U[0]
}

// Translate translates a Plucker vector along the given vector.
func (p *PluckerVec) Translate(dx *Vec) {
	p.V[0] += -dx[1]*p.U[2] + dx[2]*p.U[1]
    p.V[1] += -dx[2]*p.U[0] + dx[0]*p.U[2]
    p.V[2] += -dx[0]*p.U[1] + dx[1]*p.U[0]
}

// Dot computes the permuted inner product of p1 and p2, i.e.
// p1.U*p2.V + p1.V*p2.U.
func (p1 *PluckerVec) Dot(p2 *PluckerVec, flip bool) float32 {
	var sum float32
	for i := 0; i < 3; i++ {
		sum += p1.U[i]*p2.V[i] + p1.V[i]*p2.U[i]
	}
	if flip {
		return sum
	} else {
		return -sum
	}
}

// Dot computes the permuted inner product of p1 and p2, i.e.
// p1.U*p2.V + p1.V*p2.U and also returns a sign flag of -1, 0, or +1 if
// that product is negative, zero, or positive, respectively.
func (p1 *PluckerVec) SignDot(p2 *PluckerVec, flip bool) (float32, int) {
	dot := p1.Dot(p2, flip)
	if dot == 0 {
		return dot, 0
	} else if dot > 0 {
		return dot, +1
	} else {
		return dot, -1
	}
}

// Init initializes an anchored Plucker vector given a ray origin, P, and a
// unit direction vector, L.
func (ap *AnchoredPluckerVec) Init(P, L *Vec) {
	ap.PluckerVec.Init(P, L)
	ap.P = *P
}

// InitFromSegment initialized a Plucker vector which corresponds to a ray
// pointing from the position vector P1 to the position vector P2.
func (ap *AnchoredPluckerVec) InitFromSegment(P1, P2 *Vec) {
	ap.PluckerVec.InitFromSegment(P1, P2)
	ap.P = *P1
}

// Translate translates a Plucker vector along the given vector.
func (ap *AnchoredPluckerVec) Translate(dx *Vec) {
	ap.PluckerVec.Translate(dx)
	for i := 0; i < 3; i++ { ap.P[i] += dx[i] }
}

// Tetra is a tetrahedron. (Duh!)
//
// Face ordering is:
// F0(V3, V2, V1)
// F1(V2, V3, V0)
// F2(V1, V0, V3)
// F3(V0, V1, V2)
type Tetra [4]Vec

var tetraIdxs = [4][3]int {
	[3]int{ 3, 2, 1 },
	[3]int{ 2, 3, 0 },
	[3]int{ 1, 0, 3 },
	[3]int{ 0, 1, 2 },
}

// VertexIdx returns the index into the given tetrahedron corresponding to
// the specified face and vertex.
func (_ *Tetra) VertexIdx(face, vertex int) int {
	return tetraIdxs[face][vertex]
}

// Orient arranges tetrahedron points so that all faces point outward for
// dir = +1 and inward for dir = -1.
func (t *Tetra) Orient(dir int) {
	v, w, n := Vec{}, Vec{}, Vec{}
	for i := 0; i < 3; i++ {
		v[i] = t[1][i] - t[0][i]
		w[i] = t[2][i] - t[0][i]
	}
	n[0] = -v[1]*w[2] + v[2]*w[1]
    n[1] = -v[2]*w[0] + v[0]*w[2]
    n[2] = -v[0]*w[1] + v[1]*w[0]

	var dot float32
	for i := 0; i < 3; i++ {
		dot += n[i] * (t[3][i] - t[0][i])
	}

	if (dot < 0 && dir == -1) || (dot > 0 && dir == +1) {
		t[0], t[1] = t[1], t[0]
	}
}

// Translate translates a tetrahedron by the given vector.
func (t *Tetra) Translate(dx *Vec) {
	for i := 0; i < 4; i++ {
		for j := 0; j < 3; j++ {
			t[i][j] += dx[j]
		}
	}
}

// Sphere is exactly what you think it is.
type Sphere struct {
	X, Y, Z, R float32
}

// Intersect returns true if the two spheres intersect and false otherwise.
func (s1 *Sphere) Intersect(s2 *Sphere) bool {
	dx, dy, dz, dr := s1.X - s2.X, s1.Y - s2.Y, s1.Z - s2.Z, s1.R + s2.R
	return dr*dr > dx*dx + dy*dy + dz*dz
}

// BoundingSphere draws a bounding sphere aorund the given tetrahedron.
func (t *Tetra) BoundingSphere(sph *Sphere) {
	bx := (t[0][0] + t[1][0] + t[2][0] + t[3][0]) / 4
	by := (t[0][1] + t[1][1] + t[2][1] + t[3][1]) / 4
	bz := (t[0][2] + t[1][2] + t[2][2] + t[3][2]) / 4

	dx, dy, dz := bx-t[0][0], by-t[0][1], bz-t[0][2]
	maxRSqr := dx*dx + dy*dy + dz*dz
	for i := 1; i < 4; i++ {
		dx, dy, dz = bx-t[i][0], by-t[i][1], bz-t[i][2]
		rSqr := dx*dx + dy*dy + dz*dz
		if rSqr > maxRSqr { maxRSqr = rSqr }
	}

	sph.X, sph.Y, sph.Z = bx, by, bz
	sph.R = float32(math.Sqrt(float64(maxRSqr)))
}

// TetraFaceBary contains information specifying the barycentric coordinates
{'fraction_non_alphanumeric': 0.1871427232686721, 'fraction_numerical': 0.047043388623371755, 'mean_word_length': 3.0826319816373373, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2094144', 'n_tokens_mistral': 4998, 'n_tokens_neox': 4693, 'n_words': 1438}