# HG changeset patch # User Sascha L. Teichmann # Date 1551620741 -3600 # Node ID cb55d7eaaa3679afd23cb54574ae8dd9e8fc4996 # Parent bd46ffbb944e1c04eb246aa2335d4a8996395087 Started with the idea to clip an octree by an bounding polygon. diff -r bd46ffbb944e -r cb55d7eaaa36 cmd/octreediff/main.go --- a/cmd/octreediff/main.go Sun Mar 03 13:00:50 2019 +0100 +++ b/cmd/octreediff/main.go Sun Mar 03 14:45:41 2019 +0100 @@ -65,6 +65,23 @@ 4326 ) ` + clippingPolygonSQL = ` +WITH joined AS ( + SELECT + sr.area AS area, + sr.date_info AS date_info + FROM waterway.sounding_results sr JOIN + waterway.bottlenecks bn ON sr.bottleneck_id = bn.id + WHERE bn.bottleneck_id = $1 +) +SELECT ST_AsBinary( + ST_Buffer( + ST_intersection( + (SELECT ST_Transform(area::geometry, $2::int) FROM joined WHERE date_info = $3::date), + (SELECT ST_Transform(area::geometry, $2::int) FROM joined WHERE date_info = $4::date) + ), + 0.001)) AS area +` insertContourSQLClipped = ` WITH joined AS ( SELECT @@ -347,6 +364,21 @@ last = now log.Printf("num points: %d\n", len(result)) + var clip []byte + + if err := tx.QueryRowContext( + ctx, clippingPolygonSQL, + bottleneck, + first.EPSG, + firstDate, secondDate, + ).Scan(&clip); err != nil { + return err + } + + now = time.Now() + log.Printf("loading clipping polygon took %v\n", now.Sub(last)) + last = now + tri, err := result.triangulate() if err != nil { return err diff -r bd46ffbb944e -r cb55d7eaaa36 pkg/octree/polygon.go --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pkg/octree/polygon.go Sun Mar 03 14:45:41 2019 +0100 @@ -0,0 +1,36 @@ +// 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 + +package octree + +type Polygon struct { + // TODO: Implement me! +} + +type IntersectionType byte + +const ( + IntersectionInside IntersectionType = iota + IntersectionOutSide + IntersectionOverlaps +) + +func (p *Polygon) IntersectionBox2D(box Box2D) IntersectionType { + // TODO: Implement me + return IntersectionOutSide +} + +func (p *Polygon) IntersectionWithTriangle(t *Triangle) IntersectionType { + // TODO: Implement me + return IntersectionOutSide +} diff -r bd46ffbb944e -r cb55d7eaaa36 pkg/octree/tree.go --- a/pkg/octree/tree.go Sun Mar 03 13:00:50 2019 +0100 +++ b/pkg/octree/tree.go Sun Mar 03 14:45:41 2019 +0100 @@ -32,6 +32,11 @@ Max Vertex } +type boxFrame struct { + pos int32 + Box2D +} + func (ot *Tree) Vertices() []Vertex { return ot.vertices } @@ -43,6 +48,105 @@ {0.5, 0.5, 1.0, 1.0}, } +func (ot *Tree) Clip(p *Polygon) { + + all := Box2D{ot.Min.X, ot.Min.Y, ot.Max.X, ot.Max.Y} + + stack := []boxFrame{{1, all}} + + triChecks := make(map[int32]IntersectionType) + +frames: + for len(stack) > 0 { + top := stack[len(stack)-1] + stack = stack[:len(stack)-1] + + if top.pos > 0 { // node + switch p.IntersectionBox2D(top.Box2D) { + case IntersectionInside: + // all inside so nothing to clip. + continue frames + case IntersectionOutSide: + // all outside -> clip from tree. + if index := ot.index[top.pos:]; len(index) > 7 { + for i := range index { + index[i] = 0 + } + } + continue frames + default: // Overlaps + if index := ot.index[top.pos:]; len(index) > 7 { + children: + 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, + } + switch p.IntersectionBox2D(nbox) { + case IntersectionInside: + // all inside so nothing to clip. + continue children + case IntersectionOutSide: + // all are ouside -> clip from tree. + index[i] = 0 + index[i+4] = 0 + continue + default: // Overlaps + 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] + tris: + for i := len(indices) - 1; i >= 0; i-- { + triIndex := indices[i] + what, found := triChecks[triIndex] + if !found { + tri := ot.triangles[triIndex] + t := Triangle{ + ot.vertices[tri[0]], + ot.vertices[tri[1]], + ot.vertices[tri[2]], + } + what = p.IntersectionWithTriangle(&t) + triChecks[triIndex] = what + } + switch what { + case IntersectionInside: + // triangle inside -> stay. + continue tris + default: + // outside or not fully covered -> remove. + if i < len(indices)-1 { + copy(indices[i:], indices[i+1:]) + } + indices[len(indices)-1] = 0 + indices = indices[:len(indices)-1] + } + } + ot.index[pos] = int32(len(indices)) + } + } +} + func (ot *Tree) Value(x, y float64) (float64, bool) { // out of bounding box @@ -51,14 +155,9 @@ return 0, false } - type frame struct { - pos int32 - Box2D - } - all := Box2D{ot.Min.X, ot.Min.Y, ot.Max.X, ot.Max.Y} - stack := []frame{{1, all}} + stack := []boxFrame{{1, all}} for len(stack) > 0 { top := stack[len(stack)-1] @@ -82,10 +181,10 @@ } if nbox.Contains(x, y) { if a != 0 { - stack = append(stack, frame{a, nbox}) + stack = append(stack, boxFrame{a, nbox}) } if b != 0 { - stack = append(stack, frame{b, nbox}) + stack = append(stack, boxFrame{b, nbox}) } break } @@ -131,18 +230,13 @@ line := NewPlane2D(x1, y1, x2, y2) - type frame struct { - pos int32 - Box2D - } - 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 := []frame{{1, all}} + stack := []boxFrame{{1, all}} for len(stack) > 0 { top := stack[len(stack)-1] @@ -166,10 +260,10 @@ } if nbox.Intersects(box) && nbox.IntersectsPlane(line) { if a != 0 { - stack = append(stack, frame{a, nbox}) + stack = append(stack, boxFrame{a, nbox}) } if b != 0 { - stack = append(stack, frame{b, nbox}) + stack = append(stack, boxFrame{b, nbox}) } } }