Mercurial > gemma
changeset 2572:7686c7c23506
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 | eec11d3d74f9 |
| children | 4486ca003b55 |
| files | cmd/octreediff/main.go pkg/controllers/diff.go pkg/octree/slice.go pkg/octree/tree.go pkg/octree/vertex.go |
| diffstat | 5 files changed, 183 insertions(+), 139 deletions(-) [+] |
line wrap: on
line diff
--- a/cmd/octreediff/main.go Mon Mar 11 13:42:27 2019 +0100 +++ b/cmd/octreediff/main.go Mon Mar 11 14:00:49 2019 +0100 @@ -130,97 +130,6 @@ } } -type point struct { - x float64 - y float64 -} - -type pointMap map[point]float64 - -func (pm pointMap) triangulate() (*octree.Triangulation, error) { - vertices := make([]octree.Vertex, len(pm)) - var i int - for p, z := range pm { - vertices[i] = octree.Vertex{X: p.x, Y: p.y, Z: z} - i++ - } - return octree.Triangulate(vertices) -} - -func sliceWork( - vs []octree.Vertex, - dst pointMap, - fn func([]octree.Vertex, func([]octree.Vertex) []octree.Vertex), -) { - n := runtime.NumCPU() - - wg := new(sync.WaitGroup) - - slices := make(chan []octree.Vertex) - out := make(chan []octree.Vertex) - - pool := make(chan []octree.Vertex, n) - - const pageSize = 2048 - - turn := func(p []octree.Vertex) []octree.Vertex { - if p != nil { - out <- p - } - select { - case p = <-pool: - default: - p = make([]octree.Vertex, 0, pageSize) - } - return p - } - - for i := 0; i < n; i++ { - wg.Add(1) - go func() { - defer wg.Done() - for slice := range slices { - fn(slice, turn) - } - }() - } - done := make(chan struct{}) - go func() { - defer close(done) - for s := range out { - for i := range s { - v := &s[i] - key := point{x: v.X, y: v.Y} - if z, found := dst[key]; found { - dst[key] = (z + v.Z) * 0.5 - } else { - dst[key] = v.Z - } - } - select { - case pool <- s[:0:pageSize]: - default: - } - } - }() - - size := len(vs)/n + 1 - for len(vs) > 0 { - var l int - if len(vs) < size { - l = len(vs) - } else { - l = size - } - slices <- vs[:l] - vs = vs[l:] - } - close(slices) - wg.Wait() - close(out) - <-done -} - func process(bottleneck string, firstDate, secondDate time.Time) error { start := time.Now() defer func() { log.Printf("processing took %v\n", time.Since(start)) }() @@ -310,53 +219,7 @@ log.Printf("loading took %v\n", now.Sub(start)) last := now - firstVs, secondVs := first.Vertices(), second.Vertices() - - result := make(pointMap, len(firstVs)+len(secondVs)) - - sliceWork( - firstVs, - result, - func( - slice []octree.Vertex, - turn func([]octree.Vertex) []octree.Vertex, - ) { - p := turn(nil) - for i := range slice { - v := &slice[i] - if z, found := second.Value(v.X, v.Y); found { - p = append(p, octree.Vertex{v.X, v.Y, v.Z - z}) - if len(p) == cap(p) { - p = turn(p) - } - } - } - if len(p) > 0 { - turn(p) - } - }) - - sliceWork( - secondVs, - result, - func( - slice []octree.Vertex, - turn func([]octree.Vertex) []octree.Vertex, - ) { - p := turn(nil) - for i := range slice { - v := &slice[i] - if z, found := first.Value(v.X, v.Y); found { - p = append(p, octree.Vertex{v.X, v.Y, z - v.Z}) - if len(p) == cap(p) { - p = turn(p) - } - } - } - if len(p) > 0 { - turn(p) - } - }) + result := first.Diff(second) now = time.Now() log.Printf("setting in took %v\n", now.Sub(last)) @@ -384,7 +247,7 @@ log.Printf("loading clipping polygon took %v\n", now.Sub(last)) last = now - tri, err := result.triangulate() + tri, err := result.Triangulate() if err != nil { return err }
--- a/pkg/controllers/diff.go Mon Mar 11 13:42:27 2019 +0100 +++ b/pkg/controllers/diff.go Mon Mar 11 14:00:49 2019 +0100 @@ -89,6 +89,29 @@ return } + // We need a slow path implementation for this. + if minuendTree.EPSG != subtrahendTree.EPSG { + err = JSONError{ + Code: http.StatusInternalServerError, + Message: "Calculating differences between two different " + + "EPSG code octrees are not supported, yet.", + } + return + } + + start = time.Now() + points := minuendTree.Diff(subtrahendTree) + log.Printf("info: A - B took %v\n", time.Since(start)) + + start = time.Now() + var tri *octree.Triangulation + if tri, err = points.Triangulate(); err != nil { + return + } + log.Printf("triangulation took %v\n", time.Since(start)) + + _ = tri + // TODO: Implement me! return
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pkg/octree/slice.go Mon Mar 11 14:00:49 2019 +0100 @@ -0,0 +1,105 @@ +// This is Free Software under GNU Affero General Public License v >= 3.0 +// without warranty, see README.md and license for details. +// +// SPDX-License-Identifier: AGPL-3.0-or-later +// License-Filename: LICENSES/AGPL-3.0.txt +// +// Copyright (C) 2018, 2019 by via donau +// – Österreichische Wasserstraßen-Gesellschaft mbH +// Software engineering by Intevation GmbH +// +// Author(s): +// * Sascha L. Teichmann <sascha.teichmann@intevation.de> + +package octree + +import ( + "runtime" + "sync" +) + +type PointMap map[Point]float64 + +func (pm PointMap) Triangulate() (*Triangulation, error) { + vertices := make([]Vertex, len(pm)) + var i int + for p, z := range pm { + vertices[i] = Vertex{X: p.X, Y: p.Y, Z: z} + i++ + } + return Triangulate(vertices) +} + +func sliceWork( + vs []Vertex, + dst PointMap, + fn func([]Vertex, func([]Vertex) []Vertex), +) { + n := runtime.NumCPU() + + wg := new(sync.WaitGroup) + + slices := make(chan []Vertex) + out := make(chan []Vertex) + + pool := make(chan []Vertex, n) + + const pageSize = 2048 + + turn := func(p []Vertex) []Vertex { + if p != nil { + out <- p + } + select { + case p = <-pool: + default: + p = make([]Vertex, 0, pageSize) + } + return p + } + + for i := 0; i < n; i++ { + wg.Add(1) + go func() { + defer wg.Done() + for slice := range slices { + fn(slice, turn) + } + }() + } + done := make(chan struct{}) + go func() { + defer close(done) + for s := range out { + for i := range s { + v := &s[i] + key := Point{X: v.X, Y: v.Y} + if z, found := dst[key]; found { + dst[key] = (z + v.Z) * 0.5 + } else { + dst[key] = v.Z + } + } + select { + case pool <- s[:0:pageSize]: + default: + } + } + }() + + size := len(vs)/n + 1 + for len(vs) > 0 { + var l int + if len(vs) < size { + l = len(vs) + } else { + l = size + } + slices <- vs[:l] + vs = vs[l:] + } + close(slices) + wg.Wait() + close(out) + <-done +}
--- a/pkg/octree/tree.go Mon Mar 11 13:42:27 2019 +0100 +++ b/pkg/octree/tree.go Mon Mar 11 14:00:49 2019 +0100 @@ -268,3 +268,51 @@ } } } + +func (ot *Tree) Diff(other *Tree) PointMap { + + firstVs, secondVs := ot.Vertices(), other.Vertices() + + result := make(PointMap, len(firstVs)+len(secondVs)) + + sliceWork( + firstVs, + result, + func(slice []Vertex, turn func([]Vertex) []Vertex) { + p := turn(nil) + for i := range slice { + v := &slice[i] + if z, found := other.Value(v.X, v.Y); found { + p = append(p, Vertex{v.X, v.Y, v.Z - z}) + if len(p) == cap(p) { + p = turn(p) + } + } + } + if len(p) > 0 { + turn(p) + } + }) + + sliceWork( + secondVs, + result, + func( + slice []Vertex, turn func([]Vertex) []Vertex) { + p := turn(nil) + for i := range slice { + v := &slice[i] + if z, found := ot.Value(v.X, v.Y); found { + p = append(p, Vertex{v.X, v.Y, z - v.Z}) + if len(p) == cap(p) { + p = turn(p) + } + } + } + if len(p) > 0 { + turn(p) + } + }) + + return result +}