Mercurial > gemma
comparison pkg/mesh/raster.go @ 4827:f4abfd0ee8ad remove-octree-debris
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 } |