Skip to content

Commit e412daa

Browse files
foglemanfogleman
committed
initial commit
0 parents  commit e412daa

File tree

8 files changed

+434
-0
lines changed

8 files changed

+434
-0
lines changed

‎cmd/simplify/main.go

+17
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,17 @@
1+
package main
2+
3+
import (
4+
"fmt"
5+
"os"
6+
7+
"github.com/fogleman/simplify"
8+
)
9+
10+
func main() {
11+
path := os.Args[1]
12+
mesh, err := simplify.LoadBinarySTL(path)
13+
if err != nil {
14+
panic(err)
15+
}
16+
fmt.Println(len(mesh.Triangles))
17+
}

‎edge.go

+5
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
package simplify
2+
3+
type Edge struct {
4+
A, B Vector
5+
}

‎matrix.go

+122
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,122 @@
1+
package simplify
2+
3+
type Matrix struct {
4+
x00, x01, x02, x03 float64
5+
x10, x11, x12, x13 float64
6+
x20, x21, x22, x23 float64
7+
x30, x31, x32, x33 float64
8+
}
9+
10+
func Identity() Matrix {
11+
return Matrix{
12+
1, 0, 0, 0,
13+
0, 1, 0, 0,
14+
0, 0, 1, 0,
15+
0, 0, 0, 1}
16+
}
17+
18+
func (m Matrix) QuadricError(v Vector) float64 {
19+
return (v.X*m.x00*v.X + v.Y*m.x10*v.X + v.Z*m.x20*v.X + m.x30*v.X +
20+
v.X*m.x01*v.Y + v.Y*m.x11*v.Y + v.Z*m.x21*v.Y + m.x31*v.Y +
21+
v.X*m.x02*v.Z + v.Y*m.x12*v.Z + v.Z*m.x22*v.Z + m.x32*v.Z +
22+
v.X*m.x03 + v.Y*m.x13 + v.Z*m.x23 + m.x33)
23+
}
24+
25+
func (m Matrix) QuadricVector() Vector {
26+
n := Matrix{
27+
m.x00, m.x01, m.x02, m.x03,
28+
m.x10, m.x11, m.x12, m.x13,
29+
m.x20, m.x21, m.x22, m.x23,
30+
0, 0, 0, 1,
31+
}
32+
return n.Inverse().MulPosition(Vector{})
33+
}
34+
35+
func (a Matrix) Add(b Matrix) Matrix {
36+
return Matrix{
37+
a.x00 + b.x00, a.x10 + b.x10, a.x20 + b.x20, a.x30 + b.x30,
38+
a.x01 + b.x01, a.x11 + b.x11, a.x21 + b.x21, a.x31 + b.x31,
39+
a.x02 + b.x02, a.x12 + b.x12, a.x22 + b.x22, a.x32 + b.x32,
40+
a.x03 + b.x03, a.x13 + b.x13, a.x23 + b.x23, a.x33 + b.x33,
41+
}
42+
}
43+
44+
func (a Matrix) Mul(b Matrix) Matrix {
45+
m := Matrix{}
46+
m.x00 = a.x00*b.x00 + a.x01*b.x10 + a.x02*b.x20 + a.x03*b.x30
47+
m.x10 = a.x10*b.x00 + a.x11*b.x10 + a.x12*b.x20 + a.x13*b.x30
48+
m.x20 = a.x20*b.x00 + a.x21*b.x10 + a.x22*b.x20 + a.x23*b.x30
49+
m.x30 = a.x30*b.x00 + a.x31*b.x10 + a.x32*b.x20 + a.x33*b.x30
50+
m.x01 = a.x00*b.x01 + a.x01*b.x11 + a.x02*b.x21 + a.x03*b.x31
51+
m.x11 = a.x10*b.x01 + a.x11*b.x11 + a.x12*b.x21 + a.x13*b.x31
52+
m.x21 = a.x20*b.x01 + a.x21*b.x11 + a.x22*b.x21 + a.x23*b.x31
53+
m.x31 = a.x30*b.x01 + a.x31*b.x11 + a.x32*b.x21 + a.x33*b.x31
54+
m.x02 = a.x00*b.x02 + a.x01*b.x12 + a.x02*b.x22 + a.x03*b.x32
55+
m.x12 = a.x10*b.x02 + a.x11*b.x12 + a.x12*b.x22 + a.x13*b.x32
56+
m.x22 = a.x20*b.x02 + a.x21*b.x12 + a.x22*b.x22 + a.x23*b.x32
57+
m.x32 = a.x30*b.x02 + a.x31*b.x12 + a.x32*b.x22 + a.x33*b.x32
58+
m.x03 = a.x00*b.x03 + a.x01*b.x13 + a.x02*b.x23 + a.x03*b.x33
59+
m.x13 = a.x10*b.x03 + a.x11*b.x13 + a.x12*b.x23 + a.x13*b.x33
60+
m.x23 = a.x20*b.x03 + a.x21*b.x13 + a.x22*b.x23 + a.x23*b.x33
61+
m.x33 = a.x30*b.x03 + a.x31*b.x13 + a.x32*b.x23 + a.x33*b.x33
62+
return m
63+
}
64+
65+
func (a Matrix) MulPosition(b Vector) Vector {
66+
x := a.x00*b.X + a.x01*b.Y + a.x02*b.Z + a.x03
67+
y := a.x10*b.X + a.x11*b.Y + a.x12*b.Z + a.x13
68+
z := a.x20*b.X + a.x21*b.Y + a.x22*b.Z + a.x23
69+
return Vector{x, y, z}
70+
}
71+
72+
func (a Matrix) MulDirection(b Vector) Vector {
73+
x := a.x00*b.X + a.x01*b.Y + a.x02*b.Z
74+
y := a.x10*b.X + a.x11*b.Y + a.x12*b.Z
75+
z := a.x20*b.X + a.x21*b.Y + a.x22*b.Z
76+
return Vector{x, y, z}.Normalize()
77+
}
78+
79+
func (a Matrix) Transpose() Matrix {
80+
return Matrix{
81+
a.x00, a.x10, a.x20, a.x30,
82+
a.x01, a.x11, a.x21, a.x31,
83+
a.x02, a.x12, a.x22, a.x32,
84+
a.x03, a.x13, a.x23, a.x33}
85+
}
86+
87+
func (a Matrix) Determinant() float64 {
88+
return (a.x00*a.x11*a.x22*a.x33 - a.x00*a.x11*a.x23*a.x32 +
89+
a.x00*a.x12*a.x23*a.x31 - a.x00*a.x12*a.x21*a.x33 +
90+
a.x00*a.x13*a.x21*a.x32 - a.x00*a.x13*a.x22*a.x31 -
91+
a.x01*a.x12*a.x23*a.x30 + a.x01*a.x12*a.x20*a.x33 -
92+
a.x01*a.x13*a.x20*a.x32 + a.x01*a.x13*a.x22*a.x30 -
93+
a.x01*a.x10*a.x22*a.x33 + a.x01*a.x10*a.x23*a.x32 +
94+
a.x02*a.x13*a.x20*a.x31 - a.x02*a.x13*a.x21*a.x30 +
95+
a.x02*a.x10*a.x21*a.x33 - a.x02*a.x10*a.x23*a.x31 +
96+
a.x02*a.x11*a.x23*a.x30 - a.x02*a.x11*a.x20*a.x33 -
97+
a.x03*a.x10*a.x21*a.x32 + a.x03*a.x10*a.x22*a.x31 -
98+
a.x03*a.x11*a.x22*a.x30 + a.x03*a.x11*a.x20*a.x32 -
99+
a.x03*a.x12*a.x20*a.x31 + a.x03*a.x12*a.x21*a.x30)
100+
}
101+
102+
func (a Matrix) Inverse() Matrix {
103+
m := Matrix{}
104+
d := a.Determinant()
105+
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
106+
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
107+
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
108+
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
109+
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
110+
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
111+
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
112+
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
113+
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
114+
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
115+
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
116+
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
117+
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
118+
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
119+
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
120+
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
121+
return m
122+
}

‎mesh.go

+54
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,54 @@
1+
package simplify
2+
3+
type Mesh struct {
4+
Triangles []*Triangle
5+
}
6+
7+
func NewMesh(triangles []*Triangle) *Mesh {
8+
mesh := &Mesh{triangles}
9+
mesh.Initialize()
10+
return mesh
11+
}
12+
13+
func (mesh *Mesh) Initialize() {
14+
quadrics := make(map[Vector]Matrix)
15+
triangles := make(map[Vector][]*Triangle)
16+
edges := make(map[Edge]bool)
17+
for _, t := range mesh.Triangles {
18+
m := t.Quadric()
19+
quadrics[t.V1] = quadrics[t.V1].Add(m)
20+
quadrics[t.V2] = quadrics[t.V2].Add(m)
21+
quadrics[t.V3] = quadrics[t.V3].Add(m)
22+
triangles[t.V1] = append(triangles[t.V1], t)
23+
triangles[t.V2] = append(triangles[t.V2], t)
24+
triangles[t.V3] = append(triangles[t.V3], t)
25+
edges[Edge{t.V1, t.V2}] = true
26+
edges[Edge{t.V2, t.V3}] = true
27+
edges[Edge{t.V3, t.V1}] = true
28+
edges[Edge{t.V2, t.V1}] = true
29+
edges[Edge{t.V3, t.V2}] = true
30+
edges[Edge{t.V1, t.V3}] = true
31+
}
32+
// for i := 0; i < 100; i++ {
33+
// for edge := range edges {
34+
// // fmt.Println(edge)
35+
// a, b := edge.A, edge.B
36+
// // c := a.Midpoint(b)
37+
// m := Matrix{}
38+
// for _, t := range triangles[a] {
39+
// m = m.Add(t.Quadric())
40+
// }
41+
// for _, t := range triangles[b] {
42+
// m = m.Add(t.Quadric())
43+
// }
44+
// // d := m.QuadricVector()
45+
// // fmt.Println(d)
46+
// // fmt.Println(m.QuadricError(a), m.QuadricError(b), m.QuadricError(c), m.QuadricError(d))
47+
// // fmt.Println()
48+
// }
49+
// }
50+
}
51+
52+
func (m *Mesh) SaveBinarySTL(path string) error {
53+
return SaveBinarySTL(path, m)
54+
}

‎stl.go

+99
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,99 @@
1+
package simplify
2+
3+
import (
4+
"bufio"
5+
"encoding/binary"
6+
"os"
7+
"strings"
8+
)
9+
10+
type STLHeader struct {
11+
_ [80]uint8
12+
Count uint32
13+
}
14+
15+
type STLTriangle struct {
16+
_, V1, V2, V3 [3]float32
17+
_ uint16
18+
}
19+
20+
func LoadBinarySTL(path string) (*Mesh, error) {
21+
file, err := os.Open(path)
22+
if err != nil {
23+
return nil, err
24+
}
25+
defer file.Close()
26+
header := STLHeader{}
27+
if err := binary.Read(file, binary.LittleEndian, &header); err != nil {
28+
return nil, err
29+
}
30+
count := int(header.Count)
31+
triangles := make([]*Triangle, count)
32+
for i := 0; i < count; i++ {
33+
d := STLTriangle{}
34+
if err := binary.Read(file, binary.LittleEndian, &d); err != nil {
35+
return nil, err
36+
}
37+
v1 := Vector{float64(d.V1[0]), float64(d.V1[1]), float64(d.V1[2])}
38+
v2 := Vector{float64(d.V2[0]), float64(d.V2[1]), float64(d.V2[2])}
39+
v3 := Vector{float64(d.V3[0]), float64(d.V3[1]), float64(d.V3[2])}
40+
triangles[i] = NewTriangle(v1, v2, v3)
41+
}
42+
return NewMesh(triangles), nil
43+
}
44+
45+
func SaveBinarySTL(path string, mesh *Mesh) error {
46+
file, err := os.Create(path)
47+
if err != nil {
48+
return err
49+
}
50+
defer file.Close()
51+
header := STLHeader{}
52+
header.Count = uint32(len(mesh.Triangles))
53+
if err := binary.Write(file, binary.LittleEndian, &header); err != nil {
54+
return err
55+
}
56+
for _, triangle := range mesh.Triangles {
57+
d := STLTriangle{}
58+
d.V1[0] = float32(triangle.V1.X)
59+
d.V1[1] = float32(triangle.V1.Y)
60+
d.V1[2] = float32(triangle.V1.Z)
61+
d.V2[0] = float32(triangle.V2.X)
62+
d.V2[1] = float32(triangle.V2.Y)
63+
d.V2[2] = float32(triangle.V2.Z)
64+
d.V3[0] = float32(triangle.V3.X)
65+
d.V3[1] = float32(triangle.V3.Y)
66+
d.V3[2] = float32(triangle.V3.Z)
67+
if err := binary.Write(file, binary.LittleEndian, &d); err != nil {
68+
return err
69+
}
70+
}
71+
return nil
72+
}
73+
74+
func LoadSTL(path string) (*Mesh, error) {
75+
file, err := os.Open(path)
76+
if err != nil {
77+
return nil, err
78+
}
79+
defer file.Close()
80+
var vertexes []Vector
81+
scanner := bufio.NewScanner(file)
82+
for scanner.Scan() {
83+
line := scanner.Text()
84+
fields := strings.Fields(line)
85+
if len(fields) == 4 && fields[0] == "vertex" {
86+
f := parseFloats(fields[1:])
87+
v := Vector{f[0], f[1], f[2]}
88+
vertexes = append(vertexes, v)
89+
}
90+
}
91+
var triangles []*Triangle
92+
for i := 0; i < len(vertexes); i += 3 {
93+
v1 := vertexes[i+0]
94+
v2 := vertexes[i+1]
95+
v3 := vertexes[i+2]
96+
triangles = append(triangles, NewTriangle(v1, v2, v3))
97+
}
98+
return NewMesh(triangles), scanner.Err()
99+
}

‎triangle.go

+28
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,28 @@
1+
package simplify
2+
3+
type Triangle struct {
4+
V1, V2, V3 Vector
5+
}
6+
7+
func NewTriangle(v1, v2, v3 Vector) *Triangle {
8+
t := Triangle{}
9+
t.V1 = v1
10+
t.V2 = v2
11+
t.V3 = v3
12+
return &t
13+
}
14+
15+
func (t *Triangle) Quadric() Matrix {
16+
e1 := t.V2.Sub(t.V1)
17+
e2 := t.V3.Sub(t.V1)
18+
n := e1.Cross(e2).Normalize()
19+
x, y, z := t.V1.X, t.V1.Y, t.V1.Z
20+
a, b, c := n.X, n.Y, n.Z
21+
d := -a*x - b*y - c*z
22+
return Matrix{
23+
a * a, a * b, a * c, a * d,
24+
a * b, b * b, b * c, b * d,
25+
a * c, b * c, c * c, c * d,
26+
a * d, b * d, c * d, d * d,
27+
}
28+
}

‎util.go

+12
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,12 @@
1+
package simplify
2+
3+
import "strconv"
4+
5+
func parseFloats(items []string) []float64 {
6+
result := make([]float64, len(items))
7+
for i, item := range items {
8+
f, _ := strconv.ParseFloat(item, 64)
9+
result[i] = f
10+
}
11+
return result
12+
}

0 commit comments

Comments
 (0)