Mercurial > gemma
view pkg/mesh/raster.go @ 5520:05db984d3db1
Improve performance of bottleneck area calculation
Avoid buffer calculations by replacing them with simple distance comparisons
and calculate the boundary of the result geometry only once per iteration.
In some edge cases with very large numbers of iterations, this reduced
the runtime of a bottleneck import by a factor of more than twenty.
author | Tom Gottfried <tom@intevation.de> |
---|---|
date | Thu, 21 Oct 2021 19:50:39 +0200 |
parents | 5f47eeea988d |
children | 1222b777f51f |
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) 2019 by via donau // – Österreichische Wasserstraßen-Gesellschaft mbH // Software engineering by Intevation GmbH // // Author(s): // * Sascha L. Teichmann <sascha.teichmann@intevation.de> package mesh import ( "math" "runtime" "sync" "time" "github.com/fogleman/contourmap" "gemma.intevation.de/gemma/pkg/common" "gemma.intevation.de/gemma/pkg/log" "gemma.intevation.de/gemma/pkg/wkb" ) type Raster struct { BBox Box2D CellSize float64 XCells int YCells int Cells []float64 } const noData = -math.MaxFloat64 func NewRaster(bbox Box2D, cellSize float64) *Raster { width, height := bbox.Size() log.Infof("raster extent: %.2f / %.2f", width, height) xCells := int(math.Ceil(width / cellSize)) yCells := int(math.Ceil(height / cellSize)) log.Infof("raster size: %d / %d\n", xCells, yCells) size := (xCells + 2) * (yCells + 2) cells := make([]float64, size) for i := range cells { cells[i] = noData } return &Raster{ BBox: bbox, CellSize: cellSize, XCells: xCells, YCells: yCells, Cells: cells, } } func (r *Raster) Rasterize(eval func(float64, float64) (float64, bool)) { var wg sync.WaitGroup rows := make(chan int) rasterRow := func() { defer wg.Done() quat := 0.25 * r.CellSize for i := range rows { pos := (i+1)*(r.XCells+2) + 1 row := r.Cells[pos : pos+r.XCells] py := r.BBox.Y1 + float64(i)*r.CellSize + r.CellSize/2 px := r.BBox.X1 + r.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 += r.CellSize } } } for n := runtime.NumCPU(); n >= 1; n-- { wg.Add(1) go rasterRow() } for i := 0; i < r.YCells; i++ { rows <- i } close(rows) wg.Wait() } func (r *Raster) Diff(eval func(float64, float64) (float64, bool)) { var wg sync.WaitGroup rows := make(chan int) rasterRow := func() { defer wg.Done() quat := 0.25 * r.CellSize for i := range rows { pos := (i+1)*(r.XCells+2) + 1 row := r.Cells[pos : pos+r.XCells] py := r.BBox.Y1 + float64(i)*r.CellSize + r.CellSize/2 px := r.BBox.X1 + r.CellSize/2 y1 := py - quat y2 := py + quat for j, old := range row { // only diff where values if old == noData { px += r.CellSize continue } 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) } else { row[j] = noData } px += r.CellSize } } } for n := runtime.NumCPU(); n >= 1; n-- { wg.Add(1) go rasterRow() } for i := 0; i < r.YCells; i++ { rows <- i } close(rows) wg.Wait() } func (r *Raster) ZExtent() (float64, float64, bool) { min, max := math.MaxFloat64, -math.MaxFloat64 for _, v := range r.Cells { if v == noData { continue } if v < min { min = v } if v > max { max = v } } return min, max, min != math.MaxFloat64 } func (r *Raster) Trace(heights ClassBreaks) []wkb.MultiPolygonGeom { start := time.Now() tracer := contourmap.FromFloat64s(r.XCells+2, r.YCells+2, r.Cells) areas := make([]wkb.MultiPolygonGeom, len(heights)) reprojX := common.Linear(0.5, r.BBox.X1, 1.5, r.BBox.X1+r.CellSize) reprojY := common.Linear(0.5, r.BBox.Y1, 1.5, r.BBox.Y1+r.CellSize) var wg sync.WaitGroup cnts := make(chan int) doContours := func() { defer wg.Done() for hIdx := range cnts { if c := tracer.Contours(heights[hIdx]); len(c) > 0 { areas[hIdx] = buildMultipolygon(c, reprojX, reprojY) } } } for n := runtime.NumCPU(); n >= 1; n-- { wg.Add(1) go doContours() } for i := range heights { cnts <- i } close(cnts) wg.Wait() log.Infof("Tracing areas took %v\n", time.Since(start)) return areas } type contour []contourmap.Point type bboxNode struct { box Box2D cnt contour children []*bboxNode } func (cnt contour) contains(o contour) bool { return contains(cnt, o[0].X, o[0].Y) || contains(cnt, o[len(o)/2].X, o[len(o)/2].Y) } func (cnt contour) length() int { return len(cnt) } func (cnt contour) point(i int) (float64, float64) { return cnt[i].X, cnt[i].Y } func (cnt contour) bbox() 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} } func (bn *bboxNode) insert(cnt contour, box Box2D) { // check if children are inside new var nr *bboxNode for i, r := range bn.children { if r.box.Inside(box) && cnt.contains(r.cnt) { 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 } // check if new is inside an old for _, r := range bn.children { if box.Inside(r.box) && r.cnt.contains(cnt) { r.insert(cnt, box) return } } // its a new child node. nr = &bboxNode{box: box, cnt: cnt} bn.children = append(bn.children, nr) } func (bn *bboxNode) insertRoot(cnt contour) { bn.insert(cnt, cnt.bbox()) } type bboxOutFunc func(contour, []contour) func (bn *bboxNode) generate(out bboxOutFunc) { var grands []*bboxNode holes := make([]contour, len(bn.children)) for i, ch := range bn.children { holes[i] = ch.cnt grands = append(grands, ch.children...) } out(bn.cnt, holes) // the grand children are new polygons. for _, grand := range grands { grand.generate(out) } } func (bn *bboxNode) generateRoot(out bboxOutFunc) { for _, r := range bn.children { r.generate(out) } } func buildMultipolygon( cnts []contourmap.Contour, reprojX, reprojY func(float64) float64, ) wkb.MultiPolygonGeom { var forest bboxNode for _, cnt := range cnts { forest.insertRoot(contour(cnt)) } //log.Debugf("cnts: %d roots: %d\n", len(cnts), len(bf.roots)) var mp wkb.MultiPolygonGeom out := func(sh contour, hls []contour) { polygon := make(wkb.PolygonGeom, 1+len(hls)) // Handle shell shell := make(wkb.LinearRingGeom, len(sh)) for j, pt := range sh { shell[j] = wkb.PointGeom{ X: reprojX(pt.X), Y: reprojY(pt.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 { hole[j] = wkb.PointGeom{ X: reprojX(pt.X), Y: reprojY(pt.Y), } } if !hole.CCW() { hole.Reverse() } polygon[1+i] = hole } mp = append(mp, polygon) } forest.generateRoot(out) return mp }