Mercurial > gemma

comparison pkg/octree/raster.go @ 4768:a2f16bbcc846 direct-diff

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