Mercurial > gemma
view pkg/octree/tree.go @ 758:0f3ba8bfa641
Simplified vertical traversal of octree.
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Tue, 25 Sep 2018 04:38:40 +0200 |
parents | ef3dfe7037b3 |
children | 46fe2ae761e8 |
line wrap: on
line source
package octree import ( "log" "math" ) type Tree struct { EPSG uint32 vertices []Vertex triangles [][]int32 index []int32 Min Vertex Max Vertex } type box2d struct { x1 float64 y1 float64 x2 float64 y2 float64 } func (a box2d) intersects(b box2d) bool { return !(a.x2 < a.x1 || a.x2 < b.x1 || a.y2 < a.y1 || a.y2 < b.y1) } func (a box2d) mid() (float64, float64) { return (a.x2-a.x1)*0.5 + a.x1, (a.y2-a.y1)*0.5 + a.y1 } func (a box2d) xi(i int) float64 { if i == 0 { return a.x1 } return a.x2 } func (a box2d) yi(i int) float64 { if i == 0 { return a.y1 } return a.y2 } var scale = [4][4]float64{ {0.0, 0.0, 0.5, 0.5}, {0.5, 0.0, 1.0, 0.5}, {0.0, 0.5, 0.5, 1.0}, {0.5, 0.5, 1.0, 1.0}, } type plane2d struct { a float64 b float64 c float64 } func newPlane2d(x1, y1, x2, y2 float64) plane2d { a := x2 - x1 b := y2 - y1 l := math.Sqrt(a*a + b*b) a /= l b /= l // a*x1 + b*y1 + c = 0 c := -(a*x1 + b*y1) return plane2d{a, b, c} } func (p plane2d) eval(x, y float64) float64 { return p.a*x + p.b*y + p.c } func (ot *Tree) Vertical(x1, y1, x2, y2 float64, fn func(*Triangle)) { box := box2d{ x1: math.Min(x1, x2), y1: math.Min(y1, y2), x2: math.Max(x1, x2), y2: math.Max(y1, y2), } // out of bounding box if box.x2 < ot.Min.X || ot.Max.X < box.x1 || box.y2 < ot.Min.Y || ot.Max.Y < box.y1 { return } line := newPlane2d(x1, y1, x2, y2) type frame struct { pos int32 box2d } dupes := map[int32]struct{}{} all := box2d{ot.Min.X, ot.Min.Y, ot.Max.X, ot.Max.Y} log.Printf("area: %f\n", (box.x2-box.x1)*(box.y2-box.y1)) log.Printf("all: %f\n", (all.x2-all.x1)*(all.y2-all.y1)) stack := []frame{{1, all}} var lineRejected int var dupeRejected int for len(stack) > 0 { top := stack[len(stack)-1] stack = stack[:len(stack)-1] if top.pos > 0 { // node for i := int32(0); i < 4; i++ { a := ot.index[top.pos+i] b := ot.index[top.pos+i+4] if a == 0 && b == 0 { continue } dx := top.x2 - top.x1 dy := top.y2 - top.y1 nbox := box2d{ dx*scale[i][0] + top.x1, dy*scale[i][1] + top.y1, dx*scale[i][2] + top.x1, dy*scale[i][3] + top.y1, } if nbox.intersects(box) { if a != 0 { stack = append(stack, frame{a, nbox}) } if b != 0 { stack = append(stack, frame{b, nbox}) } } } } else { // leaf pos := -top.pos - 1 n := ot.index[pos] indices := ot.index[pos+1 : pos+1+n] for _, idx := range indices { if _, found := dupes[idx]; found { dupeRejected++ continue } tri := ot.triangles[idx] t := Triangle{ ot.vertices[tri[0]], ot.vertices[tri[1]], ot.vertices[tri[2]], } v0 := line.eval(t[0].X, t[0].Y) v1 := line.eval(t[1].X, t[1].Y) v2 := line.eval(t[2].X, t[2].Y) sides := func(s int, b bool) int { if b { return s | 2 } return s | 1 } const eps = 1e-05 if math.Abs(v0) < eps || math.Abs(v1) < eps || math.Abs(v2) < eps || sides(sides(sides(0, math.Signbit(v0)), math.Signbit(v1)), math.Signbit(v2)) == 3 { fn(&t) } else { lineRejected++ } dupes[idx] = struct{}{} } } } // XXX: This value seems to high! log.Printf("rejected by line test: %d %d %d\n", lineRejected, len(dupes), dupeRejected) } func (ot *Tree) Horizontal(h float64, fn func(*Triangle)) { if h < ot.Min.Z || ot.Max.Z < h { return } type frame struct { pos int32 min float64 max float64 } dupes := map[int32]struct{}{} stack := []frame{{1, ot.Min.Z, ot.Max.Z}} for len(stack) > 0 { top := stack[len(stack)-1] stack = stack[:len(stack)-1] pos := top.pos if pos == 0 { continue } min, max := top.min, top.max if pos > 0 { // node if mid := (max-min)*0.5 + min; h >= mid { pos += 4 // nodes with z-bit set min = mid } else { max = mid } stack = append(stack, frame{ot.index[pos+0], min, max}, frame{ot.index[pos+1], min, max}, frame{ot.index[pos+2], min, max}, frame{ot.index[pos+3], min, max}) } else { // leaf pos = -pos - 1 n := ot.index[pos] //log.Printf("%d %d %d\n", pos, n, len(ot.index)) indices := ot.index[pos+1 : pos+1+n] for _, idx := range indices { if _, found := dupes[idx]; found { continue } tri := ot.triangles[idx] t := Triangle{ ot.vertices[tri[0]], ot.vertices[tri[1]], ot.vertices[tri[2]], } if !(math.Min(t[0].Z, math.Min(t[1].Z, t[2].Z)) > h || math.Max(t[0].Z, math.Max(t[1].Z, t[2].Z)) < h) { dupes[idx] = struct{}{} fn(&t) } } } } }