Mercurial > gemma
view pkg/octree/strtree.go @ 3276:75db3199f76e
client: define stretches: fixed review button for stretches with pending import
author | Markus Kottlaender <markus@intevation.de> |
---|---|
date | Wed, 15 May 2019 17:55:38 +0200 |
parents | 114979e97a6c |
children | f456ce0a6a0e |
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" "sort" ) const numEntries = 8 type STRTree struct { tin *Tin index []int32 bboxes []Box2D } func (s *STRTree) Build(t *Tin) { s.tin = t all := make([]int32, len(t.Triangles)) for i := range all { all[i] = int32(i) } s.index = append(s.index, 0) root := s.build(all) s.index[0] = root } func (s *STRTree) Clip(p *Polygon) map[int32]struct{} { removed := make(map[int32]struct{}) stack := []int32{s.index[0]} vertices := s.tin.Vertices for len(stack) > 0 { top := stack[len(stack)-1] stack = stack[:len(stack)-1] if top > 0 { // node switch p.IntersectionBox2D(s.bbox(top)) { case IntersectionInside: // all triangles are inside polygon case IntersectionOutSide: // all triangles are outside polygon s.allTriangles(top, removed) default: // mixed bag for i, n := int32(0), s.index[top+1]; i < n; i++ { stack = append(stack, s.index[top+2+i]) } } } else { // leaf top = -top - 1 for i, n := int32(0), s.index[top+1]; i < n; i++ { idx := s.index[top+2+i] ti := s.tin.Triangles[idx] t := Triangle{ vertices[ti[0]], vertices[ti[1]], vertices[ti[2]], } if p.IntersectionWithTriangle(&t) != IntersectionInside { removed[idx] = struct{}{} } } } } return removed } func (s *STRTree) allTriangles(pos int32, tris map[int32]struct{}) { stack := []int32{pos} for len(stack) > 0 { top := stack[len(stack)-1] stack = stack[:len(stack)-1] if top > 0 { // node for i, n := int32(0), s.index[top+1]; i < n; i++ { stack = append(stack, s.index[top+2+i]) } } else { // leaf top = -top - 1 for i, n := int32(0), s.index[top+1]; i < n; i++ { tris[s.index[top+2+i]] = struct{}{} } } } } func (s *STRTree) build(items []int32) int32 { sort.Slice(items, func(i, j int) bool { ti := s.tin.Triangles[items[i]] tj := s.tin.Triangles[items[j]] return s.tin.Vertices[ti[0]].X < s.tin.Vertices[tj[0]].X }) P := int(math.Ceil(float64(len(items)) / float64(numEntries))) S := int(math.Ceil(math.Sqrt(float64(P)))) slices := strSplit(items, S) leaves := strJoin( slices, S, func(i, j int32) bool { ti := s.tin.Triangles[i] tj := s.tin.Triangles[j] return s.tin.Vertices[ti[0]].Y < s.tin.Vertices[tj[0]].Y }, s.allocLeaf, ) return s.buildNodes(leaves) } func (s *STRTree) buildNodes(items []int32) int32 { if len(items) <= numEntries { return s.allocNode(items) } sort.Slice(items, func(i, j int) bool { return s.bbox(items[i]).X1 < s.bbox(items[j]).X1 }) P := int(math.Ceil(float64(len(items)) / float64(numEntries))) S := int(math.Ceil(math.Sqrt(float64(P)))) slices := strSplit(items, S) nodes := strJoin( slices, S, func(i, j int32) bool { return s.bbox(i).Y1 < s.bbox(j).Y1 }, s.allocNode, ) return s.buildNodes(nodes) } func (s *STRTree) bbox(idx int32) Box2D { if idx < 0 { // Don't care if leaf or node. idx = -idx - 1 } return s.bboxes[s.index[idx]] } func strSplit(items []int32, S int) [][]int32 { sm := S * numEntries slices := make([][]int32, S) for i := range slices { var n int if l := len(items); l < sm { n = l } else { n = sm } slices[i] = items[:n] items = items[n:] } return slices } func strJoin( slices [][]int32, S int, less func(int32, int32) bool, alloc func([]int32) int32, ) []int32 { nodes := make([]int32, 0, S*S) for _, slice := range slices { sort.Slice(slice, func(i, j int) bool { return less(slice[i], slice[j]) }) for len(slice) > 0 { var n int if l := len(slice); l >= numEntries { n = numEntries } else { n = l } nodes = append(nodes, alloc(slice[:n])) slice = slice[n:] } } return nodes } func (s *STRTree) allocNode(items []int32) int32 { pos := len(s.index) s.index = append(s.index, 0, int32(len(items))) s.index = append(s.index, items...) if len(items) > 0 { box := s.bbox(items[0]) for i := 1; i < len(items); i++ { box = box.Union(s.bbox(items[i])) } s.index[pos] = int32(s.allocBBox(box)) } return int32(pos) } func (s *STRTree) allocBBox(box Box2D) int { pos := len(s.bboxes) s.bboxes = append(s.bboxes, box) return pos } func (s *STRTree) allocLeaf(items []int32) int32 { pos := len(s.index) s.index = append(s.index, 0, int32(len(items))) s.index = append(s.index, items...) if len(items) > 0 { vertices := s.tin.Vertices ti := s.tin.Triangles[items[0]] t := Triangle{ vertices[ti[0]], vertices[ti[1]], vertices[ti[2]], } box := t.BBox() for i := 1; i < len(items); i++ { it := items[i] ti := s.tin.Triangles[it] t := Triangle{ vertices[ti[0]], vertices[ti[1]], vertices[ti[2]], } box = box.Union(t.BBox()) } s.index[pos] = int32(s.allocBBox(box)) } return int32(-(pos + 1)) }