Mercurial > gemma
view pkg/octree/tree.go @ 755:1a1a8b5f2d02
Vertical traversal of octree for cross sections. WIP.
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Mon, 24 Sep 2018 18:41:29 +0200 |
parents | e9d37300ce99 |
children | 5e14000829d1 |
line wrap: on
line source
package octree import ( "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 = [][]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}, } 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 } type frame struct { pos int32 box2d } stack := []frame{{1, box2d{ot.Min.X, ot.Min.Y, ot.Max.X, ot.Max.Y}}} for len(stack) > 0 { top := stack[len(stack)-1] stack = stack[:len(stack)-1] if top.pos > 0 { // node midx, midy := top.mid() var all uint for i := 0; i < 2; i++ { for j := 0; j < 2; j++ { x, y := box.xi(i), box.yi(j) var code uint if x > midx { code |= 1 } if y > midy { code |= 2 } all |= 1 << code } } dx := top.x2 - top.x1 dy := top.y2 - top.y1 for i := 0; all != 0; i++ { if all&1 == 1 { s := scale[i] nbox := box2d{ dx*s[0] + top.x1, dy*s[1] + top.y1, dx*s[2] + top.x1, dy*s[3] + top.y1, } // TODO: Push two frames to stack. _ = nbox } all >>= 1 } } else { // leaf } } } 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 } stack := []frame{{1, ot.Min.Z, ot.Max.Z}} dupes := map[int32]struct{}{} 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) } } } } }