Mercurial > gemma

comparison cmd/isoareas/algo.go @ 4551:b5b23b6d76f1 iso-areas

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Move own algorith to separate file.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Tue, 01 Oct 2019 11:07:33 +0200
parents
children
comparison
equal deleted inserted replaced
4550:aa2d0006e742 4551:b5b23b6d76f1
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 main
15
16 import (
17 "bufio"
18 "container/list"
19 "log"
20 "math"
21 "math/bits"
22 "os"
23 "runtime"
24 "sync"
25
26 "gemma.intevation.de/gemma/pkg/octree"
27 )
28
29 type indexedArc struct {
30 arc int32
31 index int32
32 }
33
34 func glue(a, b octree.LineStringZ) octree.LineStringZ {
35
36 a1, a2 := a[0], a[len(a)-1]
37 b1, b2 := b[0], b[len(b)-1]
38
39 if a1.EpsEquals2D(b2) {
40 c := make(octree.LineStringZ, len(a)-1+len(b))
41 copy(c, b)
42 copy(c[len(b):], a[1:])
43 return c
44 }
45
46 if a2.EpsEquals2D(b1) {
47 c := make(octree.LineStringZ, len(a)-1+len(b))
48 copy(c, a)
49 copy(c[len(a):], b[1:])
50 return c
51 }
52
53 if a1.EpsEquals2D(b1) {
54 c := make(octree.LineStringZ, len(a)-1+len(b))
55 copy(c, b)
56 c[:len(b)].Reverse()
57 copy(c[len(b):], a[1:])
58 return c
59 }
60
61 if a2.EpsEquals2D(b2) {
62 c := make(octree.LineStringZ, len(a)-1+len(b))
63 copy(c, a)
64 copy(c[len(a):], b[:len(b)-1])
65 c[len(a):].Reverse()
66 return c
67 }
68
69 return nil
70 }
71
72 func connectArcs(tri *octree.Triangulation, cuts []indexedArc, arcs *[]octree.LineStringZ) {
73
74 unique := map[int32]struct{}{}
75 for i := range cuts {
76 unique[cuts[i].arc] = struct{}{}
77 }
78 before := len(unique)
79
80 origLen := int32(len(*arcs))
81
82 for i := range cuts {
83 if cuts[i].arc >= origLen {
84 // already has a connected arc assigned.
85 continue
86 }
87
88 line := (*arcs)[cuts[i].arc]
89
90 var modified []int
91 visited := map[int32]struct{}{}
92
93 var recursive func(int32)
94 recursive = func(idx int32) {
95 visited[idx] = struct{}{}
96
97 ns := neighbors(tri, idx)
98 for _, n := range ns {
99 n /= 3
100 if _, already := visited[n]; already {
101 continue
102 }
103
104 arcIdx := findArc(n, cuts)
105 if arcIdx == -1 {
106 continue
107 }
108 arc := cuts[arcIdx].arc
109 if arc >= origLen {
110 // already a new arc.
111 continue
112 }
113 oline := (*arcs)[arc]
114
115 nline := glue(line, oline)
116 if nline == nil {
117 // not joinable
118 continue
119 }
120 line = nline
121 modified = append(modified, arcIdx)
122 recursive(cuts[arcIdx].index)
123 }
124 }
125 recursive(cuts[i].index)
126 if len(modified) > 0 {
127 // alloc a new line an fix the references.
128 nidx := int32(len(*arcs))
129 *arcs = append(*arcs, line)
130 cuts[i].arc = nidx
131 for _, j := range modified {
132 cuts[j].arc = nidx
133 }
134 }
135 }
136
137 unique = map[int32]struct{}{}
138 for i := range cuts {
139 unique[cuts[i].arc] = struct{}{}
140 }
141 log.Printf("unique arcs: before: %d after %d\n",
142 before, len(unique))
143 }
144
145 func intersectTriangles(
146 tri *octree.Triangulation,
147 heights []float64,
148 min, max octree.Vertex,
149 ) [][]octree.LineStringZ {
150
151 cuts := make([][]indexedArc, len(heights))
152
153 classes := make([][]int32, len(heights)+1)
154
155 var arcs []octree.LineStringZ
156
157 all:
158 for i := int32(0); i < int32(len(tri.Triangles)); i += 3 {
159 ti := tri.Triangles[i : i+3]
160 v0 := tri.Points[ti[0]]
161 v1 := tri.Points[ti[1]]
162 v2 := tri.Points[ti[2]]
163 min := math.Min(v0.Z, math.Min(v1.Z, v2.Z))
164
165 if min > heights[len(heights)-1] {
166 classes[len(classes)-1] = append(classes[len(classes)-1], i/3)
167 continue
168 }
169 max := math.Max(v0.Z, math.Max(v1.Z, v2.Z))
170
171 for j, h := range heights {
172
173 var l float64
174 if j > 0 {
175 l = heights[j-1]
176 } else {
177 l = -math.MaxFloat64
178 }
179
180 if l < min && max < h {
181 classes[j] = append(classes[j], i/3)
182 continue all
183 }
184 if min > h || max < h {
185 continue
186 }
187 t := octree.Triangle{v0, v1, v2}
188 cut := t.IntersectHorizontal(h)
189 if len(cut) >= 2 {
190 arc := int32(len(arcs))
191 arcs = append(arcs, cut)
192 cuts[j] = append(cuts[j], indexedArc{arc: arc, index: i / 3})
193 }
194 }
195 }
196
197 // connect the arcs in a cut list to longer arcs.
198
199 for _, c := range cuts {
200 connectArcs(tri, c, &arcs)
201 }
202
203 func() {
204 out, err := os.Create("arcs.svg")
205 if err != nil {
206 log.Printf("err: %v\n", err)
207 return
208 }
209 defer func() {
210 out.Close()
211 log.Println("writing arcs done")
212 }()
213
214 buf := bufio.NewWriter(out)
215
216 dumpArcsSVG(
217 buf,
218 min, max,
219 cuts,
220 arcs)
221 buf.Flush()
222 }()
223
224 result := make([][]octree.LineStringZ, len(classes))
225
226 jobs := make(chan int)
227
228 var wg sync.WaitGroup
229
230 worker := func() {
231 defer wg.Done()
232
233 pb := polygonBuilder{open: list.New()}
234
235 for i := range jobs {
236 usedArcs := map[int32]struct{}{}
237 var dupes int
238 var isolated, inside, found int
239 c := classes[i]
240
241 allInClass:
242 for _, idx := range c {
243 ns := neighbors(tri, idx)
244 mask := where(ns, c)
245 switch bits.OnesCount8(mask) {
246 case 3:
247 // Totally insides do not contribute to the geometries.
248 inside++
249 continue allInClass
250 case 0:
251 // Isolated are areas w/o connections to the rest.
252 isolated++
253 ti := tri.Triangles[idx*3 : idx*3+3]
254 pb.addTriangle(
255 tri.Points[ti[0]],
256 tri.Points[ti[1]],
257 tri.Points[ti[2]])
258 continue allInClass
259 default:
260 ti := tri.Triangles[idx*3 : idx*3+3]
261 for j := 0; j < 3; j++ {
262 if (mask & (1 << j)) == 0 {
263
264 var curr octree.LineStringZ
265
266 for l := i - 1; l <= i; l++ {
267 if l < 0 || l >= len(cuts) {
268 continue
269 }
270 arcIdx := findArc(ns[j]/3, cuts[l])
271 if arcIdx == -1 {
272 continue
273 }
274 found++
275 aIdx := cuts[l][arcIdx].arc
276 if _, already := usedArcs[aIdx]; already {
277 dupes++
278 continue
279 }
280 usedArcs[aIdx] = struct{}{}
281
282 curr = arcs[aIdx]
283 break
284 }
285
286 if curr == nil {
287 k := (j + 1) % 3
288 front := tri.Points[ti[j]]
289 back := tri.Points[ti[k]]
290 curr = octree.LineStringZ{front, back}
291 }
292
293 segment:
294 for e := pb.open.Front(); e != nil; e = e.Next() {
295 line := e.Value.(octree.LineStringZ)
296 if nline := glue(curr, line); nline != nil {
297 curr = nline
298 pb.open.Remove(e)
299 goto segment
300 }
301 } // all open
302
303 // check if closed
304 if curr[0].EpsEquals2D(curr[len(curr)-1]) {
305 if !curr.CCW() {
306 curr.Reverse()
307 }
308 pb.polygons = append(pb.polygons, curr)
309 } else {
310 pb.open.PushFront(curr)
311 }
312 } // for all border parts
313 }
314 }
315 }
316
317 log.Printf("\t%d: inside: %d / isolated: %d open: %d closed: %d dupes: %d found: %d\n",
318 i, inside, isolated, pb.open.Len(), len(pb.polygons), dupes, found)
319 /*
320 for e := pb.open.Front(); e != nil; e = e.Next() {
321 line := e.Value.(octree.LineStringZ)
322 if !line.CCW() {
323 line.Reverse()
324 }
325 pb.polygons = append(pb.polygons, line)
326 }
327 */
328
329 result[i] = pb.polygons
330 pb.reset()
331 }
332 }
333
334 for n := runtime.NumCPU(); n >= 0; n-- {
335 wg.Add(1)
336 go worker()
337 }
338
339 log.Println("inside classes:")
340 for i := range classes {
341 jobs <- i
342 }
343 close(jobs)
344
345 wg.Wait()
346
347 log.Println("cuts:")
348 for i, c := range cuts {
349 log.Printf("\t%.3f: %d\n", heights[i], len(c))
350 }
351
352 return result
353
354 // TODO: sew segments together.
355
356 }
357
358 type polygonBuilder struct {
359 polygons []octree.LineStringZ
360
361 open *list.List
362 }
363
364 func (pb *polygonBuilder) addTriangle(v0, v1, v2 octree.Vertex) {
365 polygon := octree.LineStringZ{v0, v1, v2, v0}
366 pb.polygons = append(pb.polygons, polygon)
367 }
368
369 func (pb *polygonBuilder) reset() {
370 pb.polygons = pb.polygons[:0]
371 pb.open.Init()
372 }
373
374 func neighbors(t *octree.Triangulation, idx int32) []int32 {
375 idx *= 3
376 return t.Halfedges[idx : idx+3]
377 }
378
379 func findArc(needle int32, haystack []indexedArc) int {
380 lo, hi := 0, len(haystack)-1
381 for lo <= hi {
382 mid := (hi-lo)/2 + lo
383 switch v := haystack[mid].index; {
384 case v < needle:
385 lo = mid + 1
386 case v > needle:
387 hi = mid - 1
388 default:
389 return mid
390 }
391 }
392 return -1
393 }
394
395 func contains(needle int32, haystack []int32) bool {
396 lo, hi := 0, len(haystack)-1
397 for lo <= hi {
398 mid := (hi-lo)/2 + lo
399 switch v := haystack[mid]; {
400 case v < needle:
401 lo = mid + 1
402 case v > needle:
403 hi = mid - 1
404 default:
405 return true
406 }
407 }
408 return false
409 }
410
411 func where(neighbors, indices []int32) byte {
412 var mask byte
413 for i, n := range neighbors {
414 if n >= 0 && contains(n/3, indices) {
415 mask |= 1 << i
416 }
417 }
418 return mask
419 }