Mercurial > gemma
view pkg/octree/areas.go @ 4727:1d48be66ce69 stack-polygons
Fixed swapped coordinates ... doh!
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Thu, 17 Oct 2019 21:31:03 +0200 |
parents | c91e759007da |
children | 76dbeab4a0d6 |
line wrap: on
line source
// 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 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 ( "log" "math" "runtime" "sync" "time" "github.com/fogleman/contourmap" "gemma.intevation.de/gemma/pkg/common" "gemma.intevation.de/gemma/pkg/wkb" ) func GenerateRandomVertices( n int, min, max Vertex, eval func(float64, float64) (float64, bool), callback func([]Vertex), ) { var wg sync.WaitGroup jobs := make(chan int) out := make(chan []Vertex) done := make(chan struct{}) cpus := runtime.NumCPU() free := make(chan []Vertex, cpus) for i := 0; i < cpus; i++ { wg.Add(1) go func() { defer wg.Done() xRange := common.Random(min.X, max.X) yRange := common.Random(min.Y, max.Y) for size := range jobs { var vertices []Vertex select { case vertices = <-free: default: vertices = make([]Vertex, 0, 1000) } for len(vertices) < size { x, y := xRange(), yRange() if z, ok := eval(x, y); ok { vertices = append(vertices, Vertex{X: x, Y: y, Z: z}) } } out <- vertices } }() } go func() { defer close(done) for vertices := range out { callback(vertices) select { case free <- vertices[:0]: default: } } }() for remain := n; remain > 0; { if remain > 1000 { jobs <- 1000 remain -= 1000 } else { jobs <- remain remain = 0 } } close(jobs) wg.Wait() close(out) <-done } func TraceAreas( heights []float64, cellSize float64, min, max Vertex, eval func(float64, float64) (float64, bool), ) []wkb.MultiPolygonGeom { width := max.X - min.X height := max.Y - min.Y log.Printf("info: Width/Height: %.2f / %.2f\n", width, height) xcells := int(math.Ceil(width / cellSize)) ycells := int(math.Ceil(height / cellSize)) log.Printf("info: Raster size: (%d, %d)\n", xcells, ycells) start := time.Now() // Add border for closing raster := make([]float64, (xcells+2)*(ycells+2)) // prefill for no data const nodata = -math.MaxFloat64 for i := range raster { raster[i] = nodata } // rasterize the height model var wg sync.WaitGroup rows := make(chan int) rasterRow := func() { defer wg.Done() quat := 0.25 * cellSize for i := range rows { pos := (i+1)*(xcells+2) + 1 row := raster[pos : pos+xcells] py := min.Y + float64(i)*cellSize + cellSize/2 px := min.X + cellSize/2 y1 := py - quat y2 := py + quat for j := range row { var n int var sum float64 if v, ok := eval(px-quat, y1); ok { sum = v n = 1 } if v, ok := eval(px-quat, y2); ok { sum += v n++ } if v, ok := eval(px+quat, y1); ok { sum += v n++ } if v, ok := eval(px+quat, y2); ok { sum += v n++ } if n > 0 { row[j] = sum / float64(n) } px += cellSize } } } for n := runtime.NumCPU(); n >= 1; n-- { wg.Add(1) go rasterRow() } for i := 0; i < ycells; i++ { rows <- i } close(rows) wg.Wait() log.Printf("info: Rastering took %v\n", time.Since(start)) start = time.Now() tracer := contourmap.FromFloat64s(xcells+2, ycells+2, raster) areas := make([]wkb.MultiPolygonGeom, len(heights)) // TODO: Check if this correct! reprojX := common.Linear(0.5, min.X, 1.5, min.X+cellSize) reprojY := common.Linear(0.5, min.Y, 1.5, min.Y+cellSize) reproj := func(x, y float64) (float64, float64) { return reprojX(x), reprojY(y) } cnts := make(chan int) doContours := func() { defer wg.Done() for hIdx := range cnts { c := tracer.Contours(heights[hIdx]) if len(c) == 0 { continue } areas[hIdx] = buildMultipolygon(c, reproj) } } for n := runtime.NumCPU(); n >= 1; n-- { wg.Add(1) go doContours() } for i := range heights { cnts <- i } close(cnts) wg.Wait() log.Printf("info: Tracing areas took %v\n", time.Since(start)) return areas } func contourBBox(cnt contourmap.Contour) Box2D { minX, minY := math.MaxFloat64, math.MaxFloat64 maxX, maxY := -math.MaxFloat64, -math.MaxFloat64 for _, p := range cnt { if p.X < minX { minX = p.X } if p.X > maxX { maxX = p.X } if p.Y < minY { minY = p.Y } if p.Y > maxY { maxY = p.Y } } return Box2D{X1: minX, X2: maxX, Y1: minY, Y2: maxY} } type bboxNode struct { box Box2D cnt contourmap.Contour children []*bboxNode } type bboxForest struct { roots []*bboxNode } func (bn *bboxNode) insert(cnt contourmap.Contour, box Box2D) { // check if roots are inside new var nr *bboxNode for i, r := range bn.children { if r.box.Inside(box) { if nr == nil { nr = &bboxNode{box: box, cnt: cnt} } nr.children = append(nr.children, r) bn.children[i] = nil } } // we have a new child if nr != nil { // compact the list for i := len(bn.children) - 1; i >= 0; i-- { if bn.children[i] == nil { if i < len(bn.children)-1 { copy(bn.children[i:], bn.children[i+1:]) } bn.children[len(bn.children)-1] = nil bn.children = bn.children[:len(bn.children)-1] } } bn.children = append(bn.children, nr) return } var inserted bool // check if new is inside an old for _, r := range bn.children { if box.Inside(r.box) { r.insert(cnt, box) inserted = true } } // its a new root node. if !inserted { nr = &bboxNode{box: box, cnt: cnt} bn.children = append(bn.children, nr) } } func (bf *bboxForest) insert(cnt contourmap.Contour) { box := contourBBox(cnt) // check if roots are inside new var nr *bboxNode for i, r := range bf.roots { if r.box.Inside(box) { if nr == nil { nr = &bboxNode{box: box, cnt: cnt} } nr.children = append(nr.children, r) bf.roots[i] = nil } } // we have a new root if nr != nil { // compact the list for i := len(bf.roots) - 1; i >= 0; i-- { if bf.roots[i] == nil { if i < len(bf.roots)-1 { copy(bf.roots[i:], bf.roots[i+1:]) } bf.roots[len(bf.roots)-1] = nil bf.roots = bf.roots[:len(bf.roots)-1] } } bf.roots = append(bf.roots, nr) return } var inserted bool // check if new is inside an old for _, r := range bf.roots { if box.Inside(r.box) { r.insert(cnt, box) inserted = true } } // its a new root node. if !inserted { nr = &bboxNode{box: box, cnt: cnt} bf.roots = append(bf.roots, nr) } } type bboxOutFunc func(contourmap.Contour, []contourmap.Contour) func (bn *bboxNode) generate(shell bool, out bboxOutFunc) { if shell && len(bn.children) == 0 { out(bn.cnt, nil) return } } func (bf *bboxForest) generate(out bboxOutFunc) { for _, r := range bf.roots { r.generate(true, out) } } func buildMultipolygon( cnts []contourmap.Contour, reproj func(float64, float64) (float64, float64), ) wkb.MultiPolygonGeom { var bf bboxForest for _, cnt := range cnts { bf.insert(cnt) } log.Printf("cnts: %d roots: %d\n", len(cnts), len(bf.roots)) var mp wkb.MultiPolygonGeom out := func(sh contourmap.Contour, hls []contourmap.Contour) { polygon := make(wkb.PolygonGeom, 1+len(hls)) // Handle shell shell := make(wkb.LinearRingGeom, len(sh)) for j, pt := range sh { x, y := reproj(pt.X, pt.Y) shell[j] = wkb.PointGeom{X: x, Y: y} } if !shell.CCW() { shell.Reverse() } polygon[0] = shell // handle holes for i, hl := range hls { hole := make(wkb.LinearRingGeom, len(hl)) for j, pt := range hl { x, y := reproj(pt.X, pt.Y) hole[j] = wkb.PointGeom{X: x, Y: y} } if hole.CCW() { hole.Reverse() } polygon[1+i] = hole } mp = append(mp, polygon) } bf.generate(out) log.Printf("generated: %d\n", len(mp)) return mp }