Mercurial > gemma
view cmd/octree2contour/vertex.go @ 687:be90ab542aa7 octree
octree: contouring: Do the math to calculate the intersection points of the triangles and the planes.
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Thu, 20 Sep 2018 11:32:03 +0200 |
parents | 120a82bd9953 |
children | 614135d69823 |
line wrap: on
line source
package main type vertex struct { x float64 y float64 z float64 } func (v *vertex) minimize(w vertex) { if w.x < v.x { v.x = w.x } if w.y < v.y { v.y = w.y } if w.z < v.z { v.z = w.z } } func (v *vertex) maximize(w vertex) { if w.x > v.x { v.x = w.x } if w.y > v.y { v.y = w.y } if w.z > v.z { v.z = w.z } } func (v vertex) sub(w vertex) vertex { return vertex{ v.x - w.x, v.y - w.y, v.z - w.z, } } func (v vertex) add(w vertex) vertex { return vertex{ v.x + w.x, v.y + w.y, v.z + w.z, } } func (v vertex) scale(s float64) vertex { return vertex{ s * v.x, s * v.y, s * v.z, } } func interpolate(v1, v2 vertex) func(vertex) vertex { v2 = v2.sub(v1) return func(s vertex) vertex { return vertex{ v2.x*s.x + v1.x, v2.y*s.y + v1.y, v2.z*s.z + v1.z, } } } func (a vertex) less(b vertex) bool { return a.x < b.x || a.y < b.y || a.z < b.z } type line [2]vertex func newLine(p1, p2 vertex) line { return line{ p2.sub(p1), p1, } } func (l line) eval(t float64) vertex { return l[0].scale(t).add(l[1]) } func (l line) intersectH(h float64) vertex { t := (h - l[1].z) / l[0].z return l.eval(t) } type triangle [3]vertex func side(z, h float64) int { switch { case z < h: return -1 case z > h: return +1 } return 0 } func (t *triangle) intersectH(h float64, points []vertex) []vertex { sides := [3]int{ side(t[0].z, h), side(t[1].z, h), side(t[2].z, h), } for i := 0; i < 3; i++ { j := (i + 1) % 3 si := sides[i] sj := sides[j] switch { case si == 0 && sj == 0: // both on plane points = append(points, t[i], t[j]) case si == 0 && sj != 0: // first on plane points = append(points, t[i]) case si != 0 && sj == 0: // second on plane points = append(points, t[j]) case si == sj: // both on same side default: // real intersection v := newLine(t[i], t[j]).intersectH(h) points = append(points, v) } } return points }