Mercurial > gemma
view pkg/octree/strtree.go @ 2512:2768f74d54ab octree-diff
Added an STRTree for the triangulation. Should be better suited for filtering out the clipped triangles.
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Tue, 05 Mar 2019 13:05:50 +0100 |
parents | |
children | 1534df518201 |
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 { tri *Triangulation index []int32 triangles [][]int32 bboxes []Box2D } func (s *STRTree) Build(t *Triangulation) { s.tri = t s.triangles = t.TriangleSlices() all := make([]int32, len(s.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) build(items []int32) int32 { sort.Slice(items, func(i, j int) bool { ti := s.triangles[items[i]] tj := s.triangles[items[j]] return s.tri.Points[ti[0]].X < s.tri.Points[tj[0]].X }) P := int(math.Ceil(float64(len(items)) / float64(numEntries))) S := int(math.Ceil(math.Sqrt(float64(P)))) sm := S * numEntries slices := make([][]int32, S) for i := range slices { end := (i + 1) * sm if end > len(items) { end = len(items) } slices[i] = items[i*sm : end] } leaves := make([]int32, 0, S*S) for _, slice := range slices { sort.Slice(slice, func(i, j int) bool { ti := s.triangles[slice[i]] tj := s.triangles[slice[j]] return s.tri.Points[ti[0]].X < s.tri.Points[tj[0]].X }) for len(slice) > 0 { n := numEntries if l := len(slice); l < numEntries { n = l } leaves = append(leaves, s.allocLeaf(slice[:n])) slice = slice[n:] } } return s.buildNodes(leaves) } func bboxIndex(idx int32) int32 { if idx < 0 { return -idx - 1 } return idx } func (s *STRTree) buildNodes(items []int32) int32 { if len(items) <= numEntries { return s.allocNode(items) } sort.Slice(items, func(i, j int) bool { bi := s.bboxes[s.index[bboxIndex(items[i])]] bj := s.bboxes[s.index[bboxIndex(items[j])]] return bi.X1 < bj.X1 }) P := int(math.Ceil(float64(len(items)) / float64(numEntries))) S := int(math.Ceil(math.Sqrt(float64(P)))) sm := S * numEntries slices := make([][]int32, S) for i := range slices { end := (i + 1) * sm if end > len(items) { end = len(items) } slices[i] = items[i*sm : end] } nodes := make([]int32, 0, S*S) for _, slice := range slices { sort.Slice(slice, func(i, j int) bool { bi := s.bboxes[s.index[bboxIndex(slice[i])]] bj := s.bboxes[s.index[bboxIndex(slice[j])]] return bi.X1 < bj.X1 }) for len(slice) > 0 { n := numEntries if l := len(slice); l < numEntries { n = l } nodes = append(nodes, s.allocNode(slice[:n])) slice = slice[n:] } } return s.buildNodes(nodes) } func (s *STRTree) allocNode(items []int32) int32 { pos := len(s.index) s.index = append(s.index, 0) s.index = append(s.index, items...) if len(items) > 0 { box := s.bboxes[s.index[bboxIndex(items[0])]] for i := 1; i < len(items); i++ { it := items[i] b := s.bboxes[s.index[bboxIndex(it)]] box = box.Union(b) } 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) s.index = append(s.index, items...) if len(items) > 0 { ti := s.triangles[items[0]] t := Triangle{ s.tri.Points[ti[0]], s.tri.Points[ti[1]], s.tri.Points[ti[2]], } box := t.BBox() for i := 1; i < len(items); i++ { it := items[i] ti := s.triangles[it] t := Triangle{ s.tri.Points[ti[0]], s.tri.Points[ti[1]], s.tri.Points[ti[2]], } box = box.Union(t.BBox()) } s.index[pos] = int32(s.allocBBox(box)) } return -int32(pos) - 1 }