initial commit

cmd/simplify/main.go

@@ -0,0 +1,17 @@
+package main
+
+import (
+	"fmt"
+	"os"
+
+	"github.com/fogleman/simplify"
+)
+
+func main() {
+	path := os.Args[1]
+	mesh, err := simplify.LoadBinarySTL(path)
+	if err != nil {
+		panic(err)
+	}
+	fmt.Println(len(mesh.Triangles))
+}

edge.go

@@ -0,0 +1,5 @@
+package simplify
+
+type Edge struct {
+	A, B Vector
+}

matrix.go

@@ -0,0 +1,122 @@
+package simplify
+
+type Matrix struct {
+	x00, x01, x02, x03 float64
+	x10, x11, x12, x13 float64
+	x20, x21, x22, x23 float64
+	x30, x31, x32, x33 float64
+}
+
+func Identity() Matrix {
+	return Matrix{
+		1, 0, 0, 0,
+		0, 1, 0, 0,
+		0, 0, 1, 0,
+		0, 0, 0, 1}
+}
+
+func (m Matrix) QuadricError(v Vector) float64 {
+	return (v.X*m.x00*v.X + v.Y*m.x10*v.X + v.Z*m.x20*v.X + m.x30*v.X +
+		v.X*m.x01*v.Y + v.Y*m.x11*v.Y + v.Z*m.x21*v.Y + m.x31*v.Y +
+		v.X*m.x02*v.Z + v.Y*m.x12*v.Z + v.Z*m.x22*v.Z + m.x32*v.Z +
+		v.X*m.x03 + v.Y*m.x13 + v.Z*m.x23 + m.x33)
+}
+
+func (m Matrix) QuadricVector() Vector {
+	n := Matrix{
+		m.x00, m.x01, m.x02, m.x03,
+		m.x10, m.x11, m.x12, m.x13,
+		m.x20, m.x21, m.x22, m.x23,
+		0, 0, 0, 1,
+	}
+	return n.Inverse().MulPosition(Vector{})
+}
+
+func (a Matrix) Add(b Matrix) Matrix {
+	return Matrix{
+		a.x00 + b.x00, a.x10 + b.x10, a.x20 + b.x20, a.x30 + b.x30,
+		a.x01 + b.x01, a.x11 + b.x11, a.x21 + b.x21, a.x31 + b.x31,
+		a.x02 + b.x02, a.x12 + b.x12, a.x22 + b.x22, a.x32 + b.x32,
+		a.x03 + b.x03, a.x13 + b.x13, a.x23 + b.x23, a.x33 + b.x33,
+	}
+}
+
+func (a Matrix) Mul(b Matrix) Matrix {
+	m := Matrix{}
+	m.x00 = a.x00*b.x00 + a.x01*b.x10 + a.x02*b.x20 + a.x03*b.x30
+	m.x10 = a.x10*b.x00 + a.x11*b.x10 + a.x12*b.x20 + a.x13*b.x30
+	m.x20 = a.x20*b.x00 + a.x21*b.x10 + a.x22*b.x20 + a.x23*b.x30
+	m.x30 = a.x30*b.x00 + a.x31*b.x10 + a.x32*b.x20 + a.x33*b.x30
+	m.x01 = a.x00*b.x01 + a.x01*b.x11 + a.x02*b.x21 + a.x03*b.x31
+	m.x11 = a.x10*b.x01 + a.x11*b.x11 + a.x12*b.x21 + a.x13*b.x31
+	m.x21 = a.x20*b.x01 + a.x21*b.x11 + a.x22*b.x21 + a.x23*b.x31
+	m.x31 = a.x30*b.x01 + a.x31*b.x11 + a.x32*b.x21 + a.x33*b.x31
+	m.x02 = a.x00*b.x02 + a.x01*b.x12 + a.x02*b.x22 + a.x03*b.x32
+	m.x12 = a.x10*b.x02 + a.x11*b.x12 + a.x12*b.x22 + a.x13*b.x32
+	m.x22 = a.x20*b.x02 + a.x21*b.x12 + a.x22*b.x22 + a.x23*b.x32
+	m.x32 = a.x30*b.x02 + a.x31*b.x12 + a.x32*b.x22 + a.x33*b.x32
+	m.x03 = a.x00*b.x03 + a.x01*b.x13 + a.x02*b.x23 + a.x03*b.x33
+	m.x13 = a.x10*b.x03 + a.x11*b.x13 + a.x12*b.x23 + a.x13*b.x33
+	m.x23 = a.x20*b.x03 + a.x21*b.x13 + a.x22*b.x23 + a.x23*b.x33
+	m.x33 = a.x30*b.x03 + a.x31*b.x13 + a.x32*b.x23 + a.x33*b.x33
+	return m
+}
+
+func (a Matrix) MulPosition(b Vector) Vector {
+	x := a.x00*b.X + a.x01*b.Y + a.x02*b.Z + a.x03
+	y := a.x10*b.X + a.x11*b.Y + a.x12*b.Z + a.x13
+	z := a.x20*b.X + a.x21*b.Y + a.x22*b.Z + a.x23
+	return Vector{x, y, z}
+}
+
+func (a Matrix) MulDirection(b Vector) Vector {
+	x := a.x00*b.X + a.x01*b.Y + a.x02*b.Z
+	y := a.x10*b.X + a.x11*b.Y + a.x12*b.Z
+	z := a.x20*b.X + a.x21*b.Y + a.x22*b.Z
+	return Vector{x, y, z}.Normalize()
+}
+
+func (a Matrix) Transpose() Matrix {
+	return Matrix{
+		a.x00, a.x10, a.x20, a.x30,
+		a.x01, a.x11, a.x21, a.x31,
+		a.x02, a.x12, a.x22, a.x32,
+		a.x03, a.x13, a.x23, a.x33}
+}
+
+func (a Matrix) Determinant() float64 {
+	return (a.x00*a.x11*a.x22*a.x33 - a.x00*a.x11*a.x23*a.x32 +
+		a.x00*a.x12*a.x23*a.x31 - a.x00*a.x12*a.x21*a.x33 +
+		a.x00*a.x13*a.x21*a.x32 - a.x00*a.x13*a.x22*a.x31 -
+		a.x01*a.x12*a.x23*a.x30 + a.x01*a.x12*a.x20*a.x33 -
+		a.x01*a.x13*a.x20*a.x32 + a.x01*a.x13*a.x22*a.x30 -
+		a.x01*a.x10*a.x22*a.x33 + a.x01*a.x10*a.x23*a.x32 +
+		a.x02*a.x13*a.x20*a.x31 - a.x02*a.x13*a.x21*a.x30 +
+		a.x02*a.x10*a.x21*a.x33 - a.x02*a.x10*a.x23*a.x31 +
+		a.x02*a.x11*a.x23*a.x30 - a.x02*a.x11*a.x20*a.x33 -
+		a.x03*a.x10*a.x21*a.x32 + a.x03*a.x10*a.x22*a.x31 -
+		a.x03*a.x11*a.x22*a.x30 + a.x03*a.x11*a.x20*a.x32 -
+		a.x03*a.x12*a.x20*a.x31 + a.x03*a.x12*a.x21*a.x30)
+}
+
+func (a Matrix) Inverse() Matrix {
+	m := Matrix{}
+	d := a.Determinant()
+	m.x00 = (a.x12*a.x23*a.x31 - a.x13*a.x22*a.x31 + a.x13*a.x21*a.x32 - a.x11*a.x23*a.x32 - a.x12*a.x21*a.x33 + a.x11*a.x22*a.x33) / d
+	m.x01 = (a.x03*a.x22*a.x31 - a.x02*a.x23*a.x31 - a.x03*a.x21*a.x32 + a.x01*a.x23*a.x32 + a.x02*a.x21*a.x33 - a.x01*a.x22*a.x33) / d
+	m.x02 = (a.x02*a.x13*a.x31 - a.x03*a.x12*a.x31 + a.x03*a.x11*a.x32 - a.x01*a.x13*a.x32 - a.x02*a.x11*a.x33 + a.x01*a.x12*a.x33) / d
+	m.x03 = (a.x03*a.x12*a.x21 - a.x02*a.x13*a.x21 - a.x03*a.x11*a.x22 + a.x01*a.x13*a.x22 + a.x02*a.x11*a.x23 - a.x01*a.x12*a.x23) / d
+	m.x10 = (a.x13*a.x22*a.x30 - a.x12*a.x23*a.x30 - a.x13*a.x20*a.x32 + a.x10*a.x23*a.x32 + a.x12*a.x20*a.x33 - a.x10*a.x22*a.x33) / d
+	m.x11 = (a.x02*a.x23*a.x30 - a.x03*a.x22*a.x30 + a.x03*a.x20*a.x32 - a.x00*a.x23*a.x32 - a.x02*a.x20*a.x33 + a.x00*a.x22*a.x33) / d
+	m.x12 = (a.x03*a.x12*a.x30 - a.x02*a.x13*a.x30 - a.x03*a.x10*a.x32 + a.x00*a.x13*a.x32 + a.x02*a.x10*a.x33 - a.x00*a.x12*a.x33) / d
+	m.x13 = (a.x02*a.x13*a.x20 - a.x03*a.x12*a.x20 + a.x03*a.x10*a.x22 - a.x00*a.x13*a.x22 - a.x02*a.x10*a.x23 + a.x00*a.x12*a.x23) / d
+	m.x20 = (a.x11*a.x23*a.x30 - a.x13*a.x21*a.x30 + a.x13*a.x20*a.x31 - a.x10*a.x23*a.x31 - a.x11*a.x20*a.x33 + a.x10*a.x21*a.x33) / d
+	m.x21 = (a.x03*a.x21*a.x30 - a.x01*a.x23*a.x30 - a.x03*a.x20*a.x31 + a.x00*a.x23*a.x31 + a.x01*a.x20*a.x33 - a.x00*a.x21*a.x33) / d
+	m.x22 = (a.x01*a.x13*a.x30 - a.x03*a.x11*a.x30 + a.x03*a.x10*a.x31 - a.x00*a.x13*a.x31 - a.x01*a.x10*a.x33 + a.x00*a.x11*a.x33) / d
+	m.x23 = (a.x03*a.x11*a.x20 - a.x01*a.x13*a.x20 - a.x03*a.x10*a.x21 + a.x00*a.x13*a.x21 + a.x01*a.x10*a.x23 - a.x00*a.x11*a.x23) / d
+	m.x30 = (a.x12*a.x21*a.x30 - a.x11*a.x22*a.x30 - a.x12*a.x20*a.x31 + a.x10*a.x22*a.x31 + a.x11*a.x20*a.x32 - a.x10*a.x21*a.x32) / d
+	m.x31 = (a.x01*a.x22*a.x30 - a.x02*a.x21*a.x30 + a.x02*a.x20*a.x31 - a.x00*a.x22*a.x31 - a.x01*a.x20*a.x32 + a.x00*a.x21*a.x32) / d
+	m.x32 = (a.x02*a.x11*a.x30 - a.x01*a.x12*a.x30 - a.x02*a.x10*a.x31 + a.x00*a.x12*a.x31 + a.x01*a.x10*a.x32 - a.x00*a.x11*a.x32) / d
+	m.x33 = (a.x01*a.x12*a.x20 - a.x02*a.x11*a.x20 + a.x02*a.x10*a.x21 - a.x00*a.x12*a.x21 - a.x01*a.x10*a.x22 + a.x00*a.x11*a.x22) / d
+	return m
+}

mesh.go

@@ -0,0 +1,54 @@
+package simplify
+
+type Mesh struct {
+	Triangles []*Triangle
+}
+
+func NewMesh(triangles []*Triangle) *Mesh {
+	mesh := &Mesh{triangles}
+	mesh.Initialize()
+	return mesh
+}
+
+func (mesh *Mesh) Initialize() {
+	quadrics := make(map[Vector]Matrix)
+	triangles := make(map[Vector][]*Triangle)
+	edges := make(map[Edge]bool)
+	for _, t := range mesh.Triangles {
+		m := t.Quadric()
+		quadrics[t.V1] = quadrics[t.V1].Add(m)
+		quadrics[t.V2] = quadrics[t.V2].Add(m)
+		quadrics[t.V3] = quadrics[t.V3].Add(m)
+		triangles[t.V1] = append(triangles[t.V1], t)
+		triangles[t.V2] = append(triangles[t.V2], t)
+		triangles[t.V3] = append(triangles[t.V3], t)
+		edges[Edge{t.V1, t.V2}] = true
+		edges[Edge{t.V2, t.V3}] = true
+		edges[Edge{t.V3, t.V1}] = true
+		edges[Edge{t.V2, t.V1}] = true
+		edges[Edge{t.V3, t.V2}] = true
+		edges[Edge{t.V1, t.V3}] = true
+	}
+	// for i := 0; i < 100; i++ {
+	// 	for edge := range edges {
+	// 		// fmt.Println(edge)
+	// 		a, b := edge.A, edge.B
+	// 		// c := a.Midpoint(b)
+	// 		m := Matrix{}
+	// 		for _, t := range triangles[a] {
+	// 			m = m.Add(t.Quadric())
+	// 		}
+	// 		for _, t := range triangles[b] {
+	// 			m = m.Add(t.Quadric())
+	// 		}
+	// 		// d := m.QuadricVector()
+	// 		// fmt.Println(d)
+	// 		// fmt.Println(m.QuadricError(a), m.QuadricError(b), m.QuadricError(c), m.QuadricError(d))
+	// 		// fmt.Println()
+	// 	}
+	// }
+}
+
+func (m *Mesh) SaveBinarySTL(path string) error {
+	return SaveBinarySTL(path, m)
+}

stl.go

@@ -0,0 +1,99 @@
+package simplify
+
+import (
+	"bufio"
+	"encoding/binary"
+	"os"
+	"strings"
+)
+
+type STLHeader struct {
+	_     [80]uint8
+	Count uint32
+}
+
+type STLTriangle struct {
+	_, V1, V2, V3 [3]float32
+	_             uint16
+}
+
+func LoadBinarySTL(path string) (*Mesh, error) {
+	file, err := os.Open(path)
+	if err != nil {
+		return nil, err
+	}
+	defer file.Close()
+	header := STLHeader{}
+	if err := binary.Read(file, binary.LittleEndian, &header); err != nil {
+		return nil, err
+	}
+	count := int(header.Count)
+	triangles := make([]*Triangle, count)
+	for i := 0; i < count; i++ {
+		d := STLTriangle{}
+		if err := binary.Read(file, binary.LittleEndian, &d); err != nil {
+			return nil, err
+		}
+		v1 := Vector{float64(d.V1[0]), float64(d.V1[1]), float64(d.V1[2])}
+		v2 := Vector{float64(d.V2[0]), float64(d.V2[1]), float64(d.V2[2])}
+		v3 := Vector{float64(d.V3[0]), float64(d.V3[1]), float64(d.V3[2])}
+		triangles[i] = NewTriangle(v1, v2, v3)
+	}
+	return NewMesh(triangles), nil
+}
+
+func SaveBinarySTL(path string, mesh *Mesh) error {
+	file, err := os.Create(path)
+	if err != nil {
+		return err
+	}
+	defer file.Close()
+	header := STLHeader{}
+	header.Count = uint32(len(mesh.Triangles))
+	if err := binary.Write(file, binary.LittleEndian, &header); err != nil {
+		return err
+	}
+	for _, triangle := range mesh.Triangles {
+		d := STLTriangle{}
+		d.V1[0] = float32(triangle.V1.X)
+		d.V1[1] = float32(triangle.V1.Y)
+		d.V1[2] = float32(triangle.V1.Z)
+		d.V2[0] = float32(triangle.V2.X)
+		d.V2[1] = float32(triangle.V2.Y)
+		d.V2[2] = float32(triangle.V2.Z)
+		d.V3[0] = float32(triangle.V3.X)
+		d.V3[1] = float32(triangle.V3.Y)
+		d.V3[2] = float32(triangle.V3.Z)
+		if err := binary.Write(file, binary.LittleEndian, &d); err != nil {
+			return err
+		}
+	}
+	return nil
+}
+
+func LoadSTL(path string) (*Mesh, error) {
+	file, err := os.Open(path)
+	if err != nil {
+		return nil, err
+	}
+	defer file.Close()
+	var vertexes []Vector
+	scanner := bufio.NewScanner(file)
+	for scanner.Scan() {
+		line := scanner.Text()
+		fields := strings.Fields(line)
+		if len(fields) == 4 && fields[0] == "vertex" {
+			f := parseFloats(fields[1:])
+			v := Vector{f[0], f[1], f[2]}
+			vertexes = append(vertexes, v)
+		}
+	}
+	var triangles []*Triangle
+	for i := 0; i < len(vertexes); i += 3 {
+		v1 := vertexes[i+0]
+		v2 := vertexes[i+1]
+		v3 := vertexes[i+2]
+		triangles = append(triangles, NewTriangle(v1, v2, v3))
+	}
+	return NewMesh(triangles), scanner.Err()
+}

triangle.go

@@ -0,0 +1,28 @@
+package simplify
+
+type Triangle struct {
+	V1, V2, V3 Vector
+}
+
+func NewTriangle(v1, v2, v3 Vector) *Triangle {
+	t := Triangle{}
+	t.V1 = v1
+	t.V2 = v2
+	t.V3 = v3
+	return &t
+}
+
+func (t *Triangle) Quadric() Matrix {
+	e1 := t.V2.Sub(t.V1)
+	e2 := t.V3.Sub(t.V1)
+	n := e1.Cross(e2).Normalize()
+	x, y, z := t.V1.X, t.V1.Y, t.V1.Z
+	a, b, c := n.X, n.Y, n.Z
+	d := -a*x - b*y - c*z
+	return Matrix{
+		a * a, a * b, a * c, a * d,
+		a * b, b * b, b * c, b * d,
+		a * c, b * c, c * c, c * d,
+		a * d, b * d, c * d, d * d,
+	}
+}

util.go

@@ -0,0 +1,12 @@
+package simplify
+
+import "strconv"
+
+func parseFloats(items []string) []float64 {
+	result := make([]float64, len(items))
+	for i, item := range items {
+		f, _ := strconv.ParseFloat(item, 64)
+		result[i] = f
+	}
+	return result
+}