Mercurial > gemma
diff pkg/octree/raster.go @ 4768:a2f16bbcc846 direct-diff
Morph differences: Directly raster A and subtract B as a raster.
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Mon, 21 Oct 2019 02:01:56 +0200 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pkg/octree/raster.go Mon Oct 21 02:01:56 2019 +0200 @@ -0,0 +1,404 @@ +// 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 octree + +import ( + "log" + "math" + "runtime" + "sync" + "time" + + "gemma.intevation.de/gemma/pkg/common" + "gemma.intevation.de/gemma/pkg/wkb" + "github.com/fogleman/contourmap" +) + +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.Printf("info raster extent: %.2f / %.2f", width, height) + + xCells := int(math.Ceil(width / cellSize)) + yCells := int(math.Ceil(height / cellSize)) + + log.Printf("info 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) +} + +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) +} + +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 []float64) []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.Printf("info: 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.Printf("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 +}