Mercurial > gemma
view pkg/octree/tree.go @ 2624:9dbaf69c7a66
Improve geoserver config to better calculate bounding boxes
* Disable the use of estimated extents for the postgis storage
configuration for geoserver, which is set via the gemma middleware.
This way we are able to get better bounding boxes for many layers
where the postgis function `ST_EstimatedExtent()` would be far off.
author | Bernhard Reiter <bernhard@intevation.de> |
---|---|
date | Wed, 13 Mar 2019 16:18:39 +0100 |
parents | 7686c7c23506 |
children | a6c671abbc35 |
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 ( "math" ) // Tree is an Octree holding triangles. type Tree struct { // EPSG is the projection. EPSG uint32 vertices []Vertex triangles [][]int32 index []int32 // Min is the lower left corner of the bbox. Min Vertex // Max is the upper right corner of the bbox. Max Vertex } type boxFrame struct { pos int32 Box2D } func (ot *Tree) Vertices() []Vertex { return ot.vertices } var scale = [4][4]float64{ {0.0, 0.0, 0.5, 0.5}, {0.5, 0.0, 1.0, 0.5}, {0.0, 0.5, 0.5, 1.0}, {0.5, 0.5, 1.0, 1.0}, } func (ot *Tree) Value(x, y float64) (float64, bool) { // out of bounding box if x < ot.Min.X || ot.Max.X < x || y < ot.Min.Y || ot.Max.Y < y { return 0, false } all := Box2D{ot.Min.X, ot.Min.Y, ot.Max.X, ot.Max.Y} stack := []boxFrame{{1, all}} for len(stack) > 0 { top := stack[len(stack)-1] stack = stack[:len(stack)-1] if top.pos > 0 { // node if index := ot.index[top.pos:]; len(index) > 7 { for i := 0; i < 4; i++ { a := index[i] b := index[i+4] if a == 0 && b == 0 { continue } dx := top.X2 - top.X1 dy := top.Y2 - top.Y1 nbox := Box2D{ dx*scale[i][0] + top.X1, dy*scale[i][1] + top.Y1, dx*scale[i][2] + top.X1, dy*scale[i][3] + top.Y1, } if nbox.Contains(x, y) { if a != 0 { stack = append(stack, boxFrame{a, nbox}) } if b != 0 { stack = append(stack, boxFrame{b, nbox}) } break } } } } else { // leaf pos := -top.pos - 1 n := ot.index[pos] indices := ot.index[pos+1 : pos+1+n] for _, idx := range indices { tri := ot.triangles[idx] t := Triangle{ ot.vertices[tri[0]], ot.vertices[tri[1]], ot.vertices[tri[2]], } if t.Contains(x, y) { return t.Plane3D().Z(x, y), true } } } } return 0, false } // Vertical does a vertical cross cut from (x1, y1) to (x2, y2). func (ot *Tree) Vertical(x1, y1, x2, y2 float64, fn func(*Triangle)) { box := Box2D{ X1: math.Min(x1, x2), Y1: math.Min(y1, y2), X2: math.Max(x1, x2), Y2: math.Max(y1, y2), } // out of bounding box if box.X2 < ot.Min.X || ot.Max.X < box.X1 || box.Y2 < ot.Min.Y || ot.Max.Y < box.Y1 { return } line := NewPlane2D(x1, y1, x2, y2) dupes := map[int32]struct{}{} all := Box2D{ot.Min.X, ot.Min.Y, ot.Max.X, ot.Max.Y} //log.Printf("area: %f\n", (box.x2-box.x1)*(box.y2-box.y1)) //log.Printf("all: %f\n", (all.x2-all.x1)*(all.y2-all.y1)) stack := []boxFrame{{1, all}} for len(stack) > 0 { top := stack[len(stack)-1] stack = stack[:len(stack)-1] if top.pos > 0 { // node if index := ot.index[top.pos:]; len(index) > 7 { for i := 0; i < 4; i++ { a := index[i] b := index[i+4] if a == 0 && b == 0 { continue } dx := top.X2 - top.X1 dy := top.Y2 - top.Y1 nbox := Box2D{ dx*scale[i][0] + top.X1, dy*scale[i][1] + top.Y1, dx*scale[i][2] + top.X1, dy*scale[i][3] + top.Y1, } if nbox.Intersects(box) && nbox.IntersectsPlane(line) { if a != 0 { stack = append(stack, boxFrame{a, nbox}) } if b != 0 { stack = append(stack, boxFrame{b, nbox}) } } } } } else { // leaf pos := -top.pos - 1 n := ot.index[pos] indices := ot.index[pos+1 : pos+1+n] for _, idx := range indices { if _, found := dupes[idx]; found { continue } tri := ot.triangles[idx] t := Triangle{ ot.vertices[tri[0]], ot.vertices[tri[1]], ot.vertices[tri[2]], } v0 := line.Eval(t[0].X, t[0].Y) v1 := line.Eval(t[1].X, t[1].Y) v2 := line.Eval(t[2].X, t[2].Y) if onPlane(v0) || onPlane(v1) || onPlane(v2) || sides(sides(sides(0, v0), v1), v2) == 3 { fn(&t) } dupes[idx] = struct{}{} } } } } // Horizontal does a horizontal cross cut. func (ot *Tree) Horizontal(h float64, fn func(*Triangle)) { if h < ot.Min.Z || ot.Max.Z < h { return } type frame struct { pos int32 min float64 max float64 } dupes := map[int32]struct{}{} stack := []frame{{1, ot.Min.Z, ot.Max.Z}} for len(stack) > 0 { top := stack[len(stack)-1] stack = stack[:len(stack)-1] pos := top.pos if pos == 0 { continue } min, max := top.min, top.max if pos > 0 { // node if mid := (max-min)*0.5 + min; h >= mid { pos += 4 // nodes with z-bit set min = mid } else { max = mid } if pos < int32(len(ot.index)) { if index := ot.index[pos:]; len(index) > 3 { stack = append(stack, frame{index[0], min, max}, frame{index[1], min, max}, frame{index[2], min, max}, frame{index[3], min, max}) } } } else { // leaf pos = -pos - 1 n := ot.index[pos] //log.Printf("%d %d %d\n", pos, n, len(ot.index)) indices := ot.index[pos+1 : pos+1+n] for _, idx := range indices { if _, found := dupes[idx]; found { continue } tri := ot.triangles[idx] t := Triangle{ ot.vertices[tri[0]], ot.vertices[tri[1]], ot.vertices[tri[2]], } if !(math.Min(t[0].Z, math.Min(t[1].Z, t[2].Z)) > h || math.Max(t[0].Z, math.Max(t[1].Z, t[2].Z)) < h) { dupes[idx] = struct{}{} fn(&t) } } } } } 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 }