Mercurial > gemma
view pkg/octree/tree.go @ 3488:aa264021e2bf zpg-ldc
merged default into zpg-ldc branch.
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Mon, 27 May 2019 16:38:07 +0200 |
parents | 7686c7c23506 |
children | a6c671abbc35 |
line wrap: on
line source
// This is Free Software under GNU Affero General Public License v >= 3.0 // without warranty, see README.md and license for details. // // SPDX-License-Identifier: AGPL-3.0-or-later // License-Filename: LICENSES/AGPL-3.0.txt // // Copyright (C) 2018 by via donau // – Österreichische Wasserstraßen-Gesellschaft mbH // Software engineering by Intevation GmbH // // Author(s): // * Sascha L. Teichmann <sascha.teichmann@intevation.de> package octree import ( "math" ) // Tree is an Octree holding triangles. type Tree struct { // EPSG is the projection. EPSG uint32 vertices []Vertex triangles [][]int32 index []int32 // Min is the lower left corner of the bbox. Min Vertex // Max is the upper right corner of the bbox. Max Vertex } type boxFrame struct { pos int32 Box2D } func (ot *Tree) Vertices() []Vertex { return ot.vertices } 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}, } func (ot *Tree) Value(x, y float64) (float64, bool) { // out of bounding box if x < ot.Min.X || ot.Max.X < x || y < ot.Min.Y || ot.Max.Y < y { return 0, false } all := Box2D{ot.Min.X, ot.Min.Y, ot.Max.X, ot.Max.Y} stack := []boxFrame{{1, all}} for len(stack) > 0 { top := stack[len(stack)-1] stack = stack[:len(stack)-1] if top.pos > 0 { // node if index := ot.index[top.pos:]; len(index) > 7 { for i := 0; i < 4; i++ { a := index[i] b := index[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.Contains(x, y) { if a != 0 { stack = append(stack, boxFrame{a, nbox}) } if b != 0 { stack = append(stack, boxFrame{b, nbox}) } break } } } } else { // leaf pos := -top.pos - 1 n := ot.index[pos] indices := ot.index[pos+1 : pos+1+n] for _, idx := range indices { tri := ot.triangles[idx] t := Triangle{ ot.vertices[tri[0]], ot.vertices[tri[1]], ot.vertices[tri[2]], } if t.Contains(x, y) { return t.Plane3D().Z(x, y), true } } } } return 0, false } // Vertical does a vertical cross cut from (x1, y1) to (x2, y2). 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) 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 := []boxFrame{{1, all}} for len(stack) > 0 { top := stack[len(stack)-1] stack = stack[:len(stack)-1] if top.pos > 0 { // node if index := ot.index[top.pos:]; len(index) > 7 { for i := 0; i < 4; i++ { a := index[i] b := index[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) && nbox.IntersectsPlane(line) { if a != 0 { stack = append(stack, boxFrame{a, nbox}) } if b != 0 { stack = append(stack, boxFrame{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 { 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) if onPlane(v0) || onPlane(v1) || onPlane(v2) || sides(sides(sides(0, v0), v1), v2) == 3 { fn(&t) } dupes[idx] = struct{}{} } } } } // Horizontal does a horizontal cross cut. 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 } if pos < int32(len(ot.index)) { if index := ot.index[pos:]; len(index) > 3 { stack = append(stack, frame{index[0], min, max}, frame{index[1], min, max}, frame{index[2], min, max}, frame{index[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) } } } } } func (ot *Tree) Diff(other *Tree) PointMap { firstVs, secondVs := ot.Vertices(), other.Vertices() result := make(PointMap, len(firstVs)+len(secondVs)) sliceWork( firstVs, result, func(slice []Vertex, turn func([]Vertex) []Vertex) { p := turn(nil) for i := range slice { v := &slice[i] if z, found := other.Value(v.X, v.Y); found { p = append(p, Vertex{v.X, v.Y, v.Z - z}) if len(p) == cap(p) { p = turn(p) } } } if len(p) > 0 { turn(p) } }) sliceWork( secondVs, result, func( slice []Vertex, turn func([]Vertex) []Vertex) { p := turn(nil) for i := range slice { v := &slice[i] if z, found := ot.Value(v.X, v.Y); found { p = append(p, Vertex{v.X, v.Y, z - v.Z}) if len(p) == cap(p) { p = turn(p) } } } if len(p) > 0 { turn(p) } }) return result }