Mercurial > gemma

comparison pkg/octree/strtree.go @ 2512:2768f74d54ab octree-diff

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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
comparison
equal deleted inserted replaced
2505:17538d937b8c 2512:2768f74d54ab
1 // This is Free Software under GNU Affero General Public License v >= 3.0
2 // without warranty, see README.md and license for details.
3 //
4 // SPDX-License-Identifier: AGPL-3.0-or-later
5 // License-Filename: LICENSES/AGPL-3.0.txt
6 //
7 // Copyright (C) 2018 by via donau
8 // – Österreichische Wasserstraßen-Gesellschaft mbH
9 // Software engineering by Intevation GmbH
10 //
11 // Author(s):
12 // * Sascha L. Teichmann <sascha.teichmann@intevation.de>
13
14 package octree
15
16 import (
17 "math"
18 "sort"
19 )
20
21 const numEntries = 8
22
23 type STRTree struct {
24 tri *Triangulation
25 index []int32
26 triangles [][]int32
27 bboxes []Box2D
28 }
29
30 func (s *STRTree) Build(t *Triangulation) {
31
32 s.tri = t
33
34 s.triangles = t.TriangleSlices()
35
36 all := make([]int32, len(s.triangles))
37
38 for i := range all {
39 all[i] = int32(i)
40 }
41
42 s.index = append(s.index, 0)
43
44 root := s.build(all)
45
46 s.index[0] = root
47 }
48
49 func (s *STRTree) build(items []int32) int32 {
50 sort.Slice(items, func(i, j int) bool {
51 ti := s.triangles[items[i]]
52 tj := s.triangles[items[j]]
53 return s.tri.Points[ti[0]].X < s.tri.Points[tj[0]].X
54 })
55
56 P := int(math.Ceil(float64(len(items)) / float64(numEntries)))
57 S := int(math.Ceil(math.Sqrt(float64(P))))
58
59 sm := S * numEntries
60
61 slices := make([][]int32, S)
62 for i := range slices {
63 end := (i + 1) * sm
64 if end > len(items) {
65 end = len(items)
66 }
67 slices[i] = items[i*sm : end]
68 }
69
70 leaves := make([]int32, 0, S*S)
71
72 for _, slice := range slices {
73 sort.Slice(slice, func(i, j int) bool {
74 ti := s.triangles[slice[i]]
75 tj := s.triangles[slice[j]]
76 return s.tri.Points[ti[0]].X < s.tri.Points[tj[0]].X
77 })
78
79 for len(slice) > 0 {
80 n := numEntries
81 if l := len(slice); l < numEntries {
82 n = l
83 }
84 leaves = append(leaves, s.allocLeaf(slice[:n]))
85 slice = slice[n:]
86 }
87 }
88
89 return s.buildNodes(leaves)
90 }
91
92 func bboxIndex(idx int32) int32 {
93 if idx < 0 {
94 return -idx - 1
95 }
96 return idx
97 }
98
99 func (s *STRTree) buildNodes(items []int32) int32 {
100
101 if len(items) <= numEntries {
102 return s.allocNode(items)
103 }
104
105 sort.Slice(items, func(i, j int) bool {
106 bi := s.bboxes[s.index[bboxIndex(items[i])]]
107 bj := s.bboxes[s.index[bboxIndex(items[j])]]
108 return bi.X1 < bj.X1
109 })
110
111 P := int(math.Ceil(float64(len(items)) / float64(numEntries)))
112 S := int(math.Ceil(math.Sqrt(float64(P))))
113
114 sm := S * numEntries
115
116 slices := make([][]int32, S)
117 for i := range slices {
118 end := (i + 1) * sm
119 if end > len(items) {
120 end = len(items)
121 }
122 slices[i] = items[i*sm : end]
123 }
124
125 nodes := make([]int32, 0, S*S)
126
127 for _, slice := range slices {
128 sort.Slice(slice, func(i, j int) bool {
129 bi := s.bboxes[s.index[bboxIndex(slice[i])]]
130 bj := s.bboxes[s.index[bboxIndex(slice[j])]]
131 return bi.X1 < bj.X1
132 })
133
134 for len(slice) > 0 {
135 n := numEntries
136 if l := len(slice); l < numEntries {
137 n = l
138 }
139 nodes = append(nodes, s.allocNode(slice[:n]))
140 slice = slice[n:]
141 }
142 }
143
144 return s.buildNodes(nodes)
145 }
146
147 func (s *STRTree) allocNode(items []int32) int32 {
148 pos := len(s.index)
149 s.index = append(s.index, 0)
150 s.index = append(s.index, items...)
151 if len(items) > 0 {
152 box := s.bboxes[s.index[bboxIndex(items[0])]]
153 for i := 1; i < len(items); i++ {
154 it := items[i]
155 b := s.bboxes[s.index[bboxIndex(it)]]
156 box = box.Union(b)
157 }
158 s.index[pos] = int32(s.allocBBox(box))
159 }
160 return int32(pos)
161 }
162
163 func (s *STRTree) allocBBox(box Box2D) int {
164 pos := len(s.bboxes)
165 s.bboxes = append(s.bboxes, box)
166 return pos
167 }
168
169 func (s *STRTree) allocLeaf(items []int32) int32 {
170 pos := len(s.index)
171 s.index = append(s.index, 0)
172 s.index = append(s.index, items...)
173 if len(items) > 0 {
174 ti := s.triangles[items[0]]
175 t := Triangle{
176 s.tri.Points[ti[0]],
177 s.tri.Points[ti[1]],
178 s.tri.Points[ti[2]],
179 }
180 box := t.BBox()
181 for i := 1; i < len(items); i++ {
182 it := items[i]
183 ti := s.triangles[it]
184 t := Triangle{
185 s.tri.Points[ti[0]],
186 s.tri.Points[ti[1]],
187 s.tri.Points[ti[2]],
188 }
189 box = box.Union(t.BBox())
190 }
191 s.index[pos] = int32(s.allocBBox(box))
192 }
193 return -int32(pos) - 1
194 }