Mercurial > gemma

comparison pkg/octree/areas.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 5c80a33edd44
children
comparison
equal deleted inserted replaced
4767:f66e5b2fa894 4768:a2f16bbcc846
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 "log"
18 "math"
19 "runtime" 17 "runtime"
20 "sync" 18 "sync"
21 "time"
22
23 "github.com/fogleman/contourmap"
24 19
25 "gemma.intevation.de/gemma/pkg/common" 20 "gemma.intevation.de/gemma/pkg/common"
26 "gemma.intevation.de/gemma/pkg/wkb"
27 ) 21 )
28 22
29 func GenerateRandomVertices( 23 func GenerateRandomVertices(
30 n int, 24 n int,
31 min, max Vertex, 25 min, max Vertex,
90 close(jobs) 84 close(jobs)
91 wg.Wait() 85 wg.Wait()
92 close(out) 86 close(out)
93 <-done 87 <-done
94 } 88 }
95
96 func TraceAreas(
97 heights []float64,
98 cellSize float64,
99 min, max Vertex,
100 eval func(float64, float64) (float64, bool),
101 ) []wkb.MultiPolygonGeom {
102
103 width := max.X - min.X
104 height := max.Y - min.Y
105
106 log.Printf("info: Width/Height: %.2f / %.2f\n", width, height)
107
108 xcells := int(math.Ceil(width / cellSize))
109 ycells := int(math.Ceil(height / cellSize))
110
111 log.Printf("info: Raster size: (%d, %d)\n", xcells, ycells)
112
113 start := time.Now()
114
115 // Add border for closing
116 raster := make([]float64, (xcells+2)*(ycells+2))
117
118 // prefill for no data
119 const nodata = -math.MaxFloat64
120 for i := range raster {
121 raster[i] = nodata
122 }
123
124 // rasterize the height model
125
126 var wg sync.WaitGroup
127
128 rows := make(chan int)
129
130 rasterRow := func() {
131 defer wg.Done()
132 quat := 0.25 * cellSize
133 for i := range rows {
134 pos := (i+1)*(xcells+2) + 1
135 row := raster[pos : pos+xcells]
136 py := min.Y + float64(i)*cellSize + cellSize/2
137 px := min.X + cellSize/2
138 y1 := py - quat
139 y2 := py + quat
140 for j := range row {
141 var n int
142 var sum float64
143
144 if v, ok := eval(px-quat, y1); ok {
145 sum = v
146 n = 1
147 }
148 if v, ok := eval(px-quat, y2); ok {
149 sum += v
150 n++
151 }
152 if v, ok := eval(px+quat, y1); ok {
153 sum += v
154 n++
155 }
156 if v, ok := eval(px+quat, y2); ok {
157 sum += v
158 n++
159 }
160
161 if n > 0 {
162 row[j] = sum / float64(n)
163 }
164 px += cellSize
165 }
166 }
167 }
168
169 for n := runtime.NumCPU(); n >= 1; n-- {
170 wg.Add(1)
171 go rasterRow()
172 }
173
174 for i := 0; i < ycells; i++ {
175 rows <- i
176 }
177 close(rows)
178
179 wg.Wait()
180 log.Printf("info: Rastering took %v\n", time.Since(start))
181
182 start = time.Now()
183
184 tracer := contourmap.FromFloat64s(xcells+2, ycells+2, raster)
185
186 areas := make([]wkb.MultiPolygonGeom, len(heights))
187
188 // TODO: Check if this correct!
189 reprojX := common.Linear(0.5, min.X, 1.5, min.X+cellSize)
190 reprojY := common.Linear(0.5, min.Y, 1.5, min.Y+cellSize)
191
192 cnts := make(chan int)
193
194 doContours := func() {
195 defer wg.Done()
196 for hIdx := range cnts {
197 c := tracer.Contours(heights[hIdx])
198
199 if len(c) == 0 {
200 continue
201 }
202
203 areas[hIdx] = buildMultipolygon(c, reprojX, reprojY)
204 }
205 }
206
207 for n := runtime.NumCPU(); n >= 1; n-- {
208 wg.Add(1)
209 go doContours()
210 }
211
212 for i := range heights {
213 cnts <- i
214 }
215 close(cnts)
216
217 wg.Wait()
218 log.Printf("info: Tracing areas took %v\n", time.Since(start))
219
220 return areas
221 }
222
223 type contour []contourmap.Point
224
225 type bboxNode struct {
226 box Box2D
227 cnt contour
228 children []*bboxNode
229 }
230
231 func (cnt contour) contains(o contour) bool {
232 return contains(cnt, o[0].X, o[0].Y) ||
233 contains(cnt, o[len(o)/2].X, o[len(o)/2].Y)
234 }
235
236 func (cnt contour) length() int {
237 return len(cnt)
238 }
239
240 func (cnt contour) point(i int) (float64, float64) {
241 return cnt[i].X, cnt[i].Y
242 }
243
244 func (cnt contour) bbox() Box2D {
245 minX, minY := math.MaxFloat64, math.MaxFloat64
246 maxX, maxY := -math.MaxFloat64, -math.MaxFloat64
247
248 for _, p := range cnt {
249 if p.X < minX {
250 minX = p.X
251 }
252 if p.X > maxX {
253 maxX = p.X
254 }
255 if p.Y < minY {
256 minY = p.Y
257 }
258 if p.Y > maxY {
259 maxY = p.Y
260 }
261 }
262 return Box2D{X1: minX, X2: maxX, Y1: minY, Y2: maxY}
263 }
264
265 func (bn *bboxNode) insert(cnt contour, box Box2D) {
266 // check if children are inside new
267 var nr *bboxNode
268
269 for i, r := range bn.children {
270 if r.box.Inside(box) && cnt.contains(r.cnt) {
271 if nr == nil {
272 nr = &bboxNode{box: box, cnt: cnt}
273 }
274 nr.children = append(nr.children, r)
275 bn.children[i] = nil
276 }
277 }
278
279 // we have a new child
280 if nr != nil {
281 // compact the list
282 for i := len(bn.children) - 1; i >= 0; i-- {
283 if bn.children[i] == nil {
284 if i < len(bn.children)-1 {
285 copy(bn.children[i:], bn.children[i+1:])
286 }
287 bn.children[len(bn.children)-1] = nil
288 bn.children = bn.children[:len(bn.children)-1]
289 }
290 }
291 bn.children = append(bn.children, nr)
292 return
293 }
294
295 // check if new is inside an old
296 for _, r := range bn.children {
297 if box.Inside(r.box) && r.cnt.contains(cnt) {
298 r.insert(cnt, box)
299 return
300 }
301 }
302
303 // its a new child node.
304 nr = &bboxNode{box: box, cnt: cnt}
305 bn.children = append(bn.children, nr)
306 }
307
308 func (bn *bboxNode) insertRoot(cnt contour) {
309 bn.insert(cnt, cnt.bbox())
310 }
311
312 type bboxOutFunc func(contour, []contour)
313
314 func (bn *bboxNode) generate(out bboxOutFunc) {
315
316 var grands []*bboxNode
317
318 holes := make([]contour, len(bn.children))
319
320 for i, ch := range bn.children {
321 holes[i] = ch.cnt
322 grands = append(grands, ch.children...)
323 }
324 out(bn.cnt, holes)
325
326 // the grand children are new polygons.
327 for _, grand := range grands {
328 grand.generate(out)
329 }
330 }
331
332 func (bn *bboxNode) generateRoot(out bboxOutFunc) {
333 for _, r := range bn.children {
334 r.generate(out)
335 }
336 }
337
338 func buildMultipolygon(
339 cnts []contourmap.Contour,
340 reprojX, reprojY func(float64) float64,
341 ) wkb.MultiPolygonGeom {
342
343 var forest bboxNode
344
345 for _, cnt := range cnts {
346 forest.insertRoot(contour(cnt))
347 }
348
349 //log.Printf("cnts: %d roots: %d\n", len(cnts), len(bf.roots))
350
351 var mp wkb.MultiPolygonGeom
352
353 out := func(sh contour, hls []contour) {
354
355 polygon := make(wkb.PolygonGeom, 1+len(hls))
356
357 // Handle shell
358 shell := make(wkb.LinearRingGeom, len(sh))
359 for j, pt := range sh {
360 shell[j] = wkb.PointGeom{
361 X: reprojX(pt.X),
362 Y: reprojY(pt.Y),
363 }
364 }
365 if shell.CCW() {
366 shell.Reverse()
367 }
368 polygon[0] = shell
369
370 // handle holes
371 for i, hl := range hls {
372 hole := make(wkb.LinearRingGeom, len(hl))
373 for j, pt := range hl {
374 hole[j] = wkb.PointGeom{
375 X: reprojX(pt.X),
376 Y: reprojY(pt.Y),
377 }
378 }
379 if !hole.CCW() {
380 hole.Reverse()
381 }
382 polygon[1+i] = hole
383 }
384
385 mp = append(mp, polygon)
386 }
387
388 forest.generateRoot(out)
389
390 return mp
391 }