Mercurial > gemma

comparison pkg/octree/tree.go @ 2516:1ec4c5633eb6 octree-diff

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Clip STR tree and not Octree.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Tue, 05 Mar 2019 17:08:16 +0100
parents e0a7571ac13a
children 7686c7c23506
comparison
equal deleted inserted replaced
2515:6bcaa8bf2603 2516:1ec4c5633eb6
12 // * Sascha L. Teichmann <sascha.teichmann@intevation.de> 12 // * Sascha L. Teichmann <sascha.teichmann@intevation.de>
13 13
14 package octree 14 package octree
15 15
16 import ( 16 import (
17 "fmt"
18 "log"
19 "math" 17 "math"
20 "os"
21 ) 18 )
22 19
23 // Tree is an Octree holding triangles. 20 // Tree is an Octree holding triangles.
24 type Tree struct { 21 type Tree struct {
25 // EPSG is the projection. 22 // EPSG is the projection.
47 var scale = [4][4]float64{ 44 var scale = [4][4]float64{
48 {0.0, 0.0, 0.5, 0.5}, 45 {0.0, 0.0, 0.5, 0.5},
49 {0.5, 0.0, 1.0, 0.5}, 46 {0.5, 0.0, 1.0, 0.5},
50 {0.0, 0.5, 0.5, 1.0}, 47 {0.0, 0.5, 0.5, 1.0},
51 {0.5, 0.5, 1.0, 1.0}, 48 {0.5, 0.5, 1.0, 1.0},
52 }
53
54 func (ot *Tree) Dump() error {
55 fake, err := os.Create("/tmp/dump.geojson")
56 if err != nil {
57 return err
58 }
59 defer fake.Close()
60
61 fmt.Fprintln(fake,
62 `{
63 "crs": {
64 "type": "name",
65 "properties": {
66 "name": "urn:ogc:def:crs:EPSG::32633"
67 }
68 },
69
70 "type": "FeatureCollection",
71 "features": [`)
72
73 first := true
74
75 already := make(map[int32]struct{})
76
77 stack := []int32{1}
78 for len(stack) > 0 {
79 pos := stack[len(stack)-1]
80 stack = stack[:len(stack)-1]
81
82 if pos > 0 { // node
83 index := ot.index[pos:]
84 if len(index) > 8 {
85 index = index[:8]
86 }
87 for i := range index {
88 if index[i] != 0 {
89 stack = append(stack, index[i])
90 }
91 }
92 } else {
93 pos = -pos - 1
94 n := ot.index[pos]
95 indices := ot.index[pos+1 : pos+1+n]
96 for _, index := range indices {
97 if _, found := already[index]; found {
98 continue
99 }
100 already[index] = struct{}{}
101 tri := ot.triangles[index]
102 t := Triangle{
103 ot.vertices[tri[0]],
104 ot.vertices[tri[1]],
105 ot.vertices[tri[2]],
106 }
107 if first {
108 first = false
109 } else {
110 fmt.Fprintln(fake, ",")
111 }
112 fmt.Fprintf(fake,
113 `{ "type": "Feature",
114 "properties": {
115 "id": %d
116 },
117 "geometry": {
118 "type": "Polygon",
119 "coordinates": [[
120 [ %.6f, %.6f ],
121 [ %.6f, %.6f ],
122 [ %.6f, %.6f ],
123 [ %.6f, %.6f ]]
124 ]
125 }
126 }`, index,
127 t[0].X, t[0].Y,
128 t[1].X, t[1].Y,
129 t[2].X, t[2].Y,
130 t[0].X, t[0].Y)
131 }
132 }
133 }
134
135 fmt.Fprintln(fake,
136 `]
137 }`)
138 return nil
139 }
140
141 func (ot *Tree) eliminate(removed map[int32]IntersectionType) {
142
143 var bad int
144
145 stack := []int32{1}
146 for len(stack) > 0 {
147 pos := stack[len(stack)-1]
148 stack = stack[:len(stack)-1]
149
150 if pos > 0 { // node
151 index := ot.index[pos:]
152 if len(index) > 8 {
153 index = index[:8]
154 }
155 for i := range index {
156 if index[i] != 0 {
157 stack = append(stack, index[i])
158 }
159 }
160 } else {
161 pos = -pos - 1
162 n := ot.index[pos]
163 indices := ot.index[pos+1 : pos+1+n]
164
165 for i := len(indices) - 1; i >= 0; i-- {
166 index := indices[i]
167 if t, found := removed[index]; found && t != IntersectionInside {
168 if i < len(indices)-1 {
169 copy(indices[i:], indices[i+1:])
170 }
171 indices[len(indices)-1] = 0
172 indices = indices[:len(indices)-1]
173 bad++
174 }
175 }
176 ot.index[pos] = int32(len(indices))
177 }
178 }
179
180 log.Printf("bad triangles: %d\n", bad)
181 }
182
183 func (ot *Tree) Clip(p *Polygon) {
184
185 /*
186 fake, _ := os.Create("/tmp/removed.geojson")
187 defer fake.Close()
188
189 fmt.Fprintln(fake,
190 `{
191 "crs": {
192 "type": "name",
193 "properties": {
194 "name": "urn:ogc:def:crs:EPSG::32633"
195 }
196 },
197
198 "type": "FeatureCollection",
199 "features": [`)
200
201 first := true
202 */
203
204 log.Printf("info: num triangles: %d\n", len(ot.triangles))
205
206 all := Box2D{ot.Min.X, ot.Min.Y, ot.Max.X, ot.Max.Y}
207
208 stack := []boxFrame{{1, all}}
209
210 triChecks := make(map[int32]IntersectionType)
211
212 var nodeTests int
213 var nodesClipped int
214 var trianglesClipped int
215 var nodesAllInside int
216
217 frames:
218 for len(stack) > 0 {
219 top := stack[len(stack)-1]
220 stack = stack[:len(stack)-1]
221
222 if top.pos > 0 { // node
223 nodeTests++
224 switch p.IntersectionBox2D(top.Box2D) {
225 case IntersectionInside:
226 // all inside so nothing to clip.
227 nodesAllInside++
228 continue frames
229 case IntersectionOutSide:
230 // all outside -> clip from tree.
231 index := ot.index[top.pos:]
232 if len(index) > 8 {
233 index = index[:8]
234 }
235 for i := range index {
236 if index[i] != 0 {
237 index[i] = 0
238 nodesClipped++
239 }
240 }
241 continue frames
242 default: // Overlaps
243 if index := ot.index[top.pos:]; len(index) > 7 {
244 children:
245 for i := 0; i < 4; i++ {
246 a := index[i]
247 b := index[i+4]
248 if a == 0 && b == 0 {
249 continue
250 }
251 dx := top.X2 - top.X1
252 dy := top.Y2 - top.Y1
253 nbox := Box2D{
254 dx*scale[i][0] + top.X1,
255 dy*scale[i][1] + top.Y1,
256 dx*scale[i][2] + top.X1,
257 dy*scale[i][3] + top.Y1,
258 }
259 switch p.IntersectionBox2D(nbox) {
260 case IntersectionInside:
261 // all inside so nothing to clip.
262 nodesAllInside++
263 continue children
264 case IntersectionOutSide:
265 // all are ouside -> clip from tree.
266 if a != 0 {
267 index[i] = 0
268 nodesClipped++
269 }
270 if b != 0 {
271 index[i+4] = 0
272 nodesClipped++
273 }
274 continue children
275 default: // Overlaps
276 if a != 0 {
277 stack = append(stack, boxFrame{a, nbox})
278 }
279 if b != 0 {
280 stack = append(stack, boxFrame{b, nbox})
281 }
282 }
283 }
284 }
285 }
286 } else { // leaf
287 pos := -top.pos - 1
288 n := ot.index[pos]
289 indices := ot.index[pos+1 : pos+1+n]
290 tris:
291 for i := len(indices) - 1; i >= 0; i-- {
292 triIndex := indices[i]
293 what, found := triChecks[triIndex]
294 if !found {
295 tri := ot.triangles[triIndex]
296 t := Triangle{
297 ot.vertices[tri[0]],
298 ot.vertices[tri[1]],
299 ot.vertices[tri[2]],
300 }
301 what = p.IntersectionWithTriangle(&t)
302 triChecks[triIndex] = what
303
304 /*
305
306 if what != IntersectionInside {
307 if first {
308 first = false
309 } else {
310 fmt.Fprintln(fake, ",")
311 }
312 fmt.Fprintf(fake,
313 `{ "type": "Feature",
314 "properties": {
315 "id": %d
316 },
317 "geometry": {
318 "type": "Polygon",
319 "coordinates": [[
320 [ %.6f, %.6f ],
321 [ %.6f, %.6f ],
322 [ %.6f, %.6f ],
323 [ %.6f, %.6f ]]
324 ]
325 }
326 }`, triIndex,
327 t[0].X, t[0].Y,
328 t[1].X, t[1].Y,
329 t[2].X, t[2].Y,
330 t[0].X, t[0].Y)
331 }
332 */
333 }
334 switch what {
335 case IntersectionInside:
336 // triangle inside -> stay.
337 continue tris
338 default:
339 trianglesClipped++
340 // outside or not fully covered -> remove.
341 if i < len(indices)-1 {
342 copy(indices[i:], indices[i+1:])
343 }
344 indices[len(indices)-1] = 0
345 indices = indices[:len(indices)-1]
346 }
347 }
348 ot.index[pos] = int32(len(indices))
349 }
350 }
351 /*
352 fmt.Fprintln(fake,
353 `]
354 }`)
355 */
356 log.Printf("info: node tests: %d\n", nodeTests)
357 log.Printf("info: nodes clipped: %d\n", nodesClipped)
358 log.Printf("info: nodes all inside: %d\n", nodesAllInside)
359 log.Printf("info: triangle clipped: %d\n", trianglesClipped)
360 log.Printf("info: tested triangles: %d\n", len(triChecks))
361
362 ot.eliminate(triChecks)
363 } 49 }
364 50
365 func (ot *Tree) Value(x, y float64) (float64, bool) { 51 func (ot *Tree) Value(x, y float64) (float64, bool) {
366 52
367 // out of bounding box 53 // out of bounding box