Mercurial > gemma

comparison pkg/mesh/raster.go @ 4827:f4abfd0ee8ad remove-octree-debris

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Renamed octree package to mesh as there is no octree any more.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Tue, 05 Nov 2019 14:30:22 +0100
parents pkg/octree/raster.go@a2f16bbcc846
children 34bc6041e61e
comparison
equal deleted inserted replaced
4826:ec5afada70ec 4827:f4abfd0ee8ad
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) 2019 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 mesh
15
16 import (
17 "log"
18 "math"
19 "runtime"
20 "sync"
21 "time"
22
23 "github.com/fogleman/contourmap"
24
25 "gemma.intevation.de/gemma/pkg/common"
26 "gemma.intevation.de/gemma/pkg/wkb"
27 )
28
29 type Raster struct {
30 BBox Box2D
31 CellSize float64
32 XCells int
33 YCells int
34 Cells []float64
35 }
36
37 const noData = -math.MaxFloat64
38
39 func NewRaster(bbox Box2D, cellSize float64) *Raster {
40
41 width, height := bbox.Size()
42
43 log.Printf("info raster extent: %.2f / %.2f", width, height)
44
45 xCells := int(math.Ceil(width / cellSize))
46 yCells := int(math.Ceil(height / cellSize))
47
48 log.Printf("info raster size: %d / %d\n", xCells, yCells)
49
50 size := (xCells + 2) * (yCells + 2)
51 cells := make([]float64, size)
52 for i := range cells {
53 cells[i] = noData
54 }
55 return &Raster{
56 BBox: bbox,
57 CellSize: cellSize,
58 XCells: xCells,
59 YCells: yCells,
60 Cells: cells,
61 }
62 }
63
64 func (r *Raster) Rasterize(eval func(float64, float64) (float64, bool)) {
65 var wg sync.WaitGroup
66
67 rows := make(chan int)
68
69 rasterRow := func() {
70 defer wg.Done()
71 quat := 0.25 * r.CellSize
72 for i := range rows {
73 pos := (i+1)*(r.XCells+2) + 1
74 row := r.Cells[pos : pos+r.XCells]
75 py := r.BBox.Y1 + float64(i)*r.CellSize + r.CellSize/2
76 px := r.BBox.X1 + r.CellSize/2
77 y1 := py - quat
78 y2 := py + quat
79 for j := range row {
80 var n int
81 var sum float64
82
83 if v, ok := eval(px-quat, y1); ok {
84 sum = v
85 n = 1
86 }
87 if v, ok := eval(px-quat, y2); ok {
88 sum += v
89 n++
90 }
91 if v, ok := eval(px+quat, y1); ok {
92 sum += v
93 n++
94 }
95 if v, ok := eval(px+quat, y2); ok {
96 sum += v
97 n++
98 }
99
100 if n > 0 {
101 row[j] = sum / float64(n)
102 }
103 px += r.CellSize
104 }
105 }
106 }
107
108 for n := runtime.NumCPU(); n >= 1; n-- {
109 wg.Add(1)
110 go rasterRow()
111 }
112
113 for i := 0; i < r.YCells; i++ {
114 rows <- i
115 }
116 close(rows)
117 }
118
119 func (r *Raster) Diff(eval func(float64, float64) (float64, bool)) {
120 var wg sync.WaitGroup
121
122 rows := make(chan int)
123
124 rasterRow := func() {
125 defer wg.Done()
126 quat := 0.25 * r.CellSize
127 for i := range rows {
128 pos := (i+1)*(r.XCells+2) + 1
129 row := r.Cells[pos : pos+r.XCells]
130 py := r.BBox.Y1 + float64(i)*r.CellSize + r.CellSize/2
131 px := r.BBox.X1 + r.CellSize/2
132 y1 := py - quat
133 y2 := py + quat
134 for j, old := range row {
135 // only diff where values
136 if old == noData {
137 px += r.CellSize
138 continue
139 }
140 var n int
141 var sum float64
142
143 if v, ok := eval(px-quat, y1); ok {
144 sum = v
145 n = 1
146 }
147 if v, ok := eval(px-quat, y2); ok {
148 sum += v
149 n++
150 }
151 if v, ok := eval(px+quat, y1); ok {
152 sum += v
153 n++
154 }
155 if v, ok := eval(px+quat, y2); ok {
156 sum += v
157 n++
158 }
159
160 if n > 0 {
161 row[j] -= sum / float64(n)
162 } else {
163 row[j] = noData
164 }
165
166 px += r.CellSize
167 }
168 }
169 }
170
171 for n := runtime.NumCPU(); n >= 1; n-- {
172 wg.Add(1)
173 go rasterRow()
174 }
175
176 for i := 0; i < r.YCells; i++ {
177 rows <- i
178 }
179 close(rows)
180 }
181
182 func (r *Raster) ZExtent() (float64, float64, bool) {
183 min, max := math.MaxFloat64, -math.MaxFloat64
184 for _, v := range r.Cells {
185 if v == noData {
186 continue
187 }
188 if v < min {
189 min = v
190 }
191 if v > max {
192 max = v
193 }
194 }
195 return min, max, min != math.MaxFloat64
196 }
197
198 func (r *Raster) Trace(heights []float64) []wkb.MultiPolygonGeom {
199 start := time.Now()
200
201 tracer := contourmap.FromFloat64s(r.XCells+2, r.YCells+2, r.Cells)
202
203 areas := make([]wkb.MultiPolygonGeom, len(heights))
204
205 reprojX := common.Linear(0.5, r.BBox.X1, 1.5, r.BBox.X1+r.CellSize)
206 reprojY := common.Linear(0.5, r.BBox.Y1, 1.5, r.BBox.Y1+r.CellSize)
207
208 var wg sync.WaitGroup
209
210 cnts := make(chan int)
211
212 doContours := func() {
213 defer wg.Done()
214 for hIdx := range cnts {
215 if c := tracer.Contours(heights[hIdx]); len(c) > 0 {
216 areas[hIdx] = buildMultipolygon(c, reprojX, reprojY)
217 }
218 }
219 }
220
221 for n := runtime.NumCPU(); n >= 1; n-- {
222 wg.Add(1)
223 go doContours()
224 }
225
226 for i := range heights {
227 cnts <- i
228 }
229 close(cnts)
230
231 wg.Wait()
232 log.Printf("info: Tracing areas took %v\n", time.Since(start))
233
234 return areas
235 }
236
237 type contour []contourmap.Point
238
239 type bboxNode struct {
240 box Box2D
241 cnt contour
242 children []*bboxNode
243 }
244
245 func (cnt contour) contains(o contour) bool {
246 return contains(cnt, o[0].X, o[0].Y) ||
247 contains(cnt, o[len(o)/2].X, o[len(o)/2].Y)
248 }
249
250 func (cnt contour) length() int {
251 return len(cnt)
252 }
253
254 func (cnt contour) point(i int) (float64, float64) {
255 return cnt[i].X, cnt[i].Y
256 }
257
258 func (cnt contour) bbox() Box2D {
259 minX, minY := math.MaxFloat64, math.MaxFloat64
260 maxX, maxY := -math.MaxFloat64, -math.MaxFloat64
261
262 for _, p := range cnt {
263 if p.X < minX {
264 minX = p.X
265 }
266 if p.X > maxX {
267 maxX = p.X
268 }
269 if p.Y < minY {
270 minY = p.Y
271 }
272 if p.Y > maxY {
273 maxY = p.Y
274 }
275 }
276 return Box2D{X1: minX, X2: maxX, Y1: minY, Y2: maxY}
277 }
278
279 func (bn *bboxNode) insert(cnt contour, box Box2D) {
280 // check if children are inside new
281 var nr *bboxNode
282
283 for i, r := range bn.children {
284 if r.box.Inside(box) && cnt.contains(r.cnt) {
285 if nr == nil {
286 nr = &bboxNode{box: box, cnt: cnt}
287 }
288 nr.children = append(nr.children, r)
289 bn.children[i] = nil
290 }
291 }
292
293 // we have a new child
294 if nr != nil {
295 // compact the list
296 for i := len(bn.children) - 1; i >= 0; i-- {
297 if bn.children[i] == nil {
298 if i < len(bn.children)-1 {
299 copy(bn.children[i:], bn.children[i+1:])
300 }
301 bn.children[len(bn.children)-1] = nil
302 bn.children = bn.children[:len(bn.children)-1]
303 }
304 }
305 bn.children = append(bn.children, nr)
306 return
307 }
308
309 // check if new is inside an old
310 for _, r := range bn.children {
311 if box.Inside(r.box) && r.cnt.contains(cnt) {
312 r.insert(cnt, box)
313 return
314 }
315 }
316
317 // its a new child node.
318 nr = &bboxNode{box: box, cnt: cnt}
319 bn.children = append(bn.children, nr)
320 }
321
322 func (bn *bboxNode) insertRoot(cnt contour) {
323 bn.insert(cnt, cnt.bbox())
324 }
325
326 type bboxOutFunc func(contour, []contour)
327
328 func (bn *bboxNode) generate(out bboxOutFunc) {
329
330 var grands []*bboxNode
331
332 holes := make([]contour, len(bn.children))
333
334 for i, ch := range bn.children {
335 holes[i] = ch.cnt
336 grands = append(grands, ch.children...)
337 }
338 out(bn.cnt, holes)
339
340 // the grand children are new polygons.
341 for _, grand := range grands {
342 grand.generate(out)
343 }
344 }
345
346 func (bn *bboxNode) generateRoot(out bboxOutFunc) {
347 for _, r := range bn.children {
348 r.generate(out)
349 }
350 }
351
352 func buildMultipolygon(
353 cnts []contourmap.Contour,
354 reprojX, reprojY func(float64) float64,
355 ) wkb.MultiPolygonGeom {
356
357 var forest bboxNode
358
359 for _, cnt := range cnts {
360 forest.insertRoot(contour(cnt))
361 }
362
363 //log.Printf("cnts: %d roots: %d\n", len(cnts), len(bf.roots))
364
365 var mp wkb.MultiPolygonGeom
366
367 out := func(sh contour, hls []contour) {
368
369 polygon := make(wkb.PolygonGeom, 1+len(hls))
370
371 // Handle shell
372 shell := make(wkb.LinearRingGeom, len(sh))
373 for j, pt := range sh {
374 shell[j] = wkb.PointGeom{
375 X: reprojX(pt.X),
376 Y: reprojY(pt.Y),
377 }
378 }
379 if shell.CCW() {
380 shell.Reverse()
381 }
382 polygon[0] = shell
383
384 // handle holes
385 for i, hl := range hls {
386 hole := make(wkb.LinearRingGeom, len(hl))
387 for j, pt := range hl {
388 hole[j] = wkb.PointGeom{
389 X: reprojX(pt.X),
390 Y: reprojY(pt.Y),
391 }
392 }
393 if !hole.CCW() {
394 hole.Reverse()
395 }
396 polygon[1+i] = hole
397 }
398
399 mp = append(mp, polygon)
400 }
401
402 forest.generateRoot(out)
403
404 return mp
405 }