Mercurial > gemma
changeset 4726:c91e759007da stack-polygons
Started with stacking shells and holes of generated polygons correctly.
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Thu, 17 Oct 2019 18:33:46 +0200 |
parents | 462d8f71da62 |
children | 1d48be66ce69 |
files | pkg/controllers/diff.go pkg/octree/areas.go pkg/octree/vertex.go |
diffstat | 3 files changed, 206 insertions(+), 31 deletions(-) [+] |
line wrap: on
line diff
--- a/pkg/controllers/diff.go Thu Oct 17 16:36:58 2019 +0200 +++ b/pkg/controllers/diff.go Thu Oct 17 18:33:46 2019 +0200 @@ -84,15 +84,8 @@ $2, ST_Transform( ST_Multi( - ST_CollectionExtract( - ST_MakeValid( - ST_Multi( - ST_Collectionextract( - ST_SimplifyPreserveTopology(ST_GeomFromWKB($4, $3::integer), $5), - 3 - ) - ) - ), + ST_Collectionextract( + ST_SimplifyPreserveTopology(ST_GeomFromWKB($4, $3::integer), $5), 3 ) ),
--- a/pkg/octree/areas.go Thu Oct 17 16:36:58 2019 +0200 +++ b/pkg/octree/areas.go Thu Oct 17 18:33:46 2019 +0200 @@ -189,6 +189,10 @@ 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() { @@ -200,28 +204,7 @@ continue } - // We need to bring it back to the - // none raster coordinate system. - a := make(wkb.MultiPolygonGeom, len(c)) - - for i, pl := range c { - shell := make(wkb.LinearRingGeom, len(pl)) - for j, pt := range pl { - shell[j] = wkb.PointGeom{ - X: reprojX(pt.X), - Y: reprojY(pt.Y), - } - } - /* - if !shell.CCW() { - log.Println("not ccw") - shell.Reverse() - } - */ - a[i] = wkb.PolygonGeom{shell} - } - - areas[hIdx] = a + areas[hIdx] = buildMultipolygon(c, reproj) } } @@ -240,3 +223,197 @@ 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 { + y, x := 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 { + y, x := 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 +}
--- a/pkg/octree/vertex.go Thu Oct 17 16:36:58 2019 +0200 +++ b/pkg/octree/vertex.go Thu Oct 17 18:33:46 2019 +0200 @@ -108,6 +108,11 @@ } } +func (a Box2D) Inside(b Box2D) bool { + return a.X1 >= b.X1 && a.X2 <= b.X2 && + a.Y1 >= b.Y1 && a.Y2 <= b.Y2 +} + func (p Plane3D) Z(x, y float64) float64 { // p.A*x + p.B*y + p.C*z + p.D = 0 return -(p.A*x + p.B*y + p.D) / p.C