Mercurial > gemma

comparison cmd/octreediff/main.go @ 2572:7686c7c23506

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Morphological differences: Moved some code into octree package.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Mon, 11 Mar 2019 14:00:49 +0100
parents 1ec4c5633eb6
children 2833ff156cb2
comparison
equal deleted inserted replaced
2571:eec11d3d74f9 2572:7686c7c23506
128 if err != nil { 128 if err != nil {
129 log.Fatalf("error: %v\n", err) 129 log.Fatalf("error: %v\n", err)
130 } 130 }
131 } 131 }
132 132
133 type point struct {
134 x float64
135 y float64
136 }
137
138 type pointMap map[point]float64
139
140 func (pm pointMap) triangulate() (*octree.Triangulation, error) {
141 vertices := make([]octree.Vertex, len(pm))
142 var i int
143 for p, z := range pm {
144 vertices[i] = octree.Vertex{X: p.x, Y: p.y, Z: z}
145 i++
146 }
147 return octree.Triangulate(vertices)
148 }
149
150 func sliceWork(
151 vs []octree.Vertex,
152 dst pointMap,
153 fn func([]octree.Vertex, func([]octree.Vertex) []octree.Vertex),
154 ) {
155 n := runtime.NumCPU()
156
157 wg := new(sync.WaitGroup)
158
159 slices := make(chan []octree.Vertex)
160 out := make(chan []octree.Vertex)
161
162 pool := make(chan []octree.Vertex, n)
163
164 const pageSize = 2048
165
166 turn := func(p []octree.Vertex) []octree.Vertex {
167 if p != nil {
168 out <- p
169 }
170 select {
171 case p = <-pool:
172 default:
173 p = make([]octree.Vertex, 0, pageSize)
174 }
175 return p
176 }
177
178 for i := 0; i < n; i++ {
179 wg.Add(1)
180 go func() {
181 defer wg.Done()
182 for slice := range slices {
183 fn(slice, turn)
184 }
185 }()
186 }
187 done := make(chan struct{})
188 go func() {
189 defer close(done)
190 for s := range out {
191 for i := range s {
192 v := &s[i]
193 key := point{x: v.X, y: v.Y}
194 if z, found := dst[key]; found {
195 dst[key] = (z + v.Z) * 0.5
196 } else {
197 dst[key] = v.Z
198 }
199 }
200 select {
201 case pool <- s[:0:pageSize]:
202 default:
203 }
204 }
205 }()
206
207 size := len(vs)/n + 1
208 for len(vs) > 0 {
209 var l int
210 if len(vs) < size {
211 l = len(vs)
212 } else {
213 l = size
214 }
215 slices <- vs[:l]
216 vs = vs[l:]
217 }
218 close(slices)
219 wg.Wait()
220 close(out)
221 <-done
222 }
223
224 func process(bottleneck string, firstDate, secondDate time.Time) error { 133 func process(bottleneck string, firstDate, secondDate time.Time) error {
225 start := time.Now() 134 start := time.Now()
226 defer func() { log.Printf("processing took %v\n", time.Since(start)) }() 135 defer func() { log.Printf("processing took %v\n", time.Since(start)) }()
227 136
228 ctx := context.Background() 137 ctx := context.Background()
308 217
309 now := time.Now() 218 now := time.Now()
310 log.Printf("loading took %v\n", now.Sub(start)) 219 log.Printf("loading took %v\n", now.Sub(start))
311 last := now 220 last := now
312 221
313 firstVs, secondVs := first.Vertices(), second.Vertices() 222 result := first.Diff(second)
314
315 result := make(pointMap, len(firstVs)+len(secondVs))
316
317 sliceWork(
318 firstVs,
319 result,
320 func(
321 slice []octree.Vertex,
322 turn func([]octree.Vertex) []octree.Vertex,
323 ) {
324 p := turn(nil)
325 for i := range slice {
326 v := &slice[i]
327 if z, found := second.Value(v.X, v.Y); found {
328 p = append(p, octree.Vertex{v.X, v.Y, v.Z - z})
329 if len(p) == cap(p) {
330 p = turn(p)
331 }
332 }
333 }
334 if len(p) > 0 {
335 turn(p)
336 }
337 })
338
339 sliceWork(
340 secondVs,
341 result,
342 func(
343 slice []octree.Vertex,
344 turn func([]octree.Vertex) []octree.Vertex,
345 ) {
346 p := turn(nil)
347 for i := range slice {
348 v := &slice[i]
349 if z, found := first.Value(v.X, v.Y); found {
350 p = append(p, octree.Vertex{v.X, v.Y, z - v.Z})
351 if len(p) == cap(p) {
352 p = turn(p)
353 }
354 }
355 }
356 if len(p) > 0 {
357 turn(p)
358 }
359 })
360 223
361 now = time.Now() 224 now = time.Now()
362 log.Printf("setting in took %v\n", now.Sub(last)) 225 log.Printf("setting in took %v\n", now.Sub(last))
363 last = now 226 last = now
364 log.Printf("num points: %d\n", len(result)) 227 log.Printf("num points: %d\n", len(result))
382 245
383 now = time.Now() 246 now = time.Now()
384 log.Printf("loading clipping polygon took %v\n", now.Sub(last)) 247 log.Printf("loading clipping polygon took %v\n", now.Sub(last))
385 last = now 248 last = now
386 249
387 tri, err := result.triangulate() 250 tri, err := result.Triangulate()
388 if err != nil { 251 if err != nil {
389 return err 252 return err
390 } 253 }
391 254
392 now = time.Now() 255 now = time.Now()