Mercurial > gemma
changeset 4737:e8df3aeb0a32
Merged stack-polygons branch back into default.
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Fri, 18 Oct 2019 11:20:33 +0200 |
parents | 7a40a39853a9 (current diff) fcbd585ae308 (diff) |
children | a0a4298a4756 |
files | |
diffstat | 5 files changed, 278 insertions(+), 67 deletions(-) [+] |
line wrap: on
line diff
--- a/pkg/controllers/diff.go Fri Oct 18 10:45:25 2019 +0200 +++ b/pkg/controllers/diff.go Fri Oct 18 11:20:33 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/imports/sr.go Fri Oct 18 10:45:25 2019 +0200 +++ b/pkg/imports/sr.go Fri Oct 18 11:20:33 2019 +0200 @@ -174,15 +174,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 Fri Oct 18 10:45:25 2019 +0200 +++ b/pkg/octree/areas.go Fri Oct 18 11:20:33 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,220 @@ return areas } + +type contour []contourmap.Point + +type bboxNode struct { + box Box2D + cnt contour + children []*bboxNode +} + +type bboxForest struct { + roots []*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) closed() bool { + return len(cnt) >= 3 +} + +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 (bf *bboxForest) insert(cnt contour) { + box := cnt.bbox() + + // check if roots are inside new + var nr *bboxNode + + for i, r := range bf.roots { + if r.box.Inside(box) && cnt.contains(r.cnt) { + 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 + } + + // check if new is inside an old + for _, r := range bf.roots { + if box.Inside(r.box) && r.cnt.contains(cnt) { + r.insert(cnt, box) + return + } + } + + // its a new root node. + nr = &bboxNode{box: box, cnt: cnt} + bf.roots = append(bf.roots, nr) +} + +type bboxOutFunc func(contour, []contour) + +func (bn *bboxNode) generate(out bboxOutFunc) { + if len(bn.children) == 0 { + out(bn.cnt, nil) + return + } + + 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) + for _, grand := range grands { + grand.generate(out) + } +} + +func (bf *bboxForest) generate(out bboxOutFunc) { + for _, r := range bf.roots { + r.generate(out) + } +} + +func buildMultipolygon( + cnts []contourmap.Contour, + reproj func(float64, float64) (float64, float64), +) wkb.MultiPolygonGeom { + + var bf bboxForest + + for _, cnt := range cnts { + bf.insert(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 { + x, y := 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 { + x, y := 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) + + return mp +}
--- a/pkg/octree/polygon.go Fri Oct 18 10:45:25 2019 +0200 +++ b/pkg/octree/polygon.go Fri Oct 18 11:20:33 2019 +0200 @@ -281,19 +281,18 @@ // No intersection -> check inside or outside // if an abritrary point is inside or not. - point := []float64{box.X1, box.Y1} // Check holes first: inside a hole means outside. if len(p.rings) > 1 { for _, hole := range p.rings[1:] { - if hole.contains(point) { + if contains(hole, box.X1, box.Y1) { return IntersectionOutSide } } } // Check shell - if p.rings[0].contains(point) { + if contains(p.rings[0], box.X1, box.Y1) { return IntersectionInside } return IntersectionOutSide @@ -324,37 +323,59 @@ } // No intersection -> check inside or outside // if an abritrary point is inside or not. - point := []float64{t[0].X, t[0].Y} + pX, pY := t[0].X, t[0].Y // Check holes first: inside a hole means outside. if len(p.rings) > 1 { for _, hole := range p.rings[1:] { - if hole.contains(point) { + if contains(hole, pX, pY) { return IntersectionOutSide } } } // Check shell - if p.rings[0].contains(point) { + if contains(p.rings[0], pX, pY) { return IntersectionInside } return IntersectionOutSide } -func (rng ring) isClosed() bool { return (len(rng) / 2) >= 3 } +func (rng ring) closed() bool { + return (len(rng) / 2) >= 3 +} + +func (rng ring) length() int { + return len(rng) / 2 +} -func (rng ring) contains(point []float64) bool { - if !rng.isClosed() { +func (rng ring) point(i int) (float64, float64) { + i *= 2 + return rng[i], rng[i+1] +} + +type segments interface { + closed() bool + length() int + point(int) (float64, float64) +} + +func contains(s segments, pX, pY float64) bool { + if !s.closed() { return false } - end := len(rng)/2 - 1 + n := s.length() + + sX, sY := s.point(0) + eX, eY := s.point(n - 1) - contains := intersectsWithRaycast(point, rng[:2], rng[end*2:end*2+2]) + contains := intersectsWithRaycast(pX, pY, sX, sY, eX, eY) - for i := 2; i < len(rng); i += 2 { - if intersectsWithRaycast(point, rng[i-2:i], rng[i:i+2]) { + for i := 1; i < n; i++ { + sX, sY := s.point(i - 1) + eX, eY := s.point(i) + if intersectsWithRaycast(pX, pY, sX, sY, eX, eY) { contains = !contains } } @@ -365,46 +386,45 @@ // Using the raycast algorithm, this returns whether or not the passed in point // Intersects with the edge drawn by the passed in start and end points. // Original implementation: http://rosettacode.org/wiki/Ray-casting_algorithm#Go -func intersectsWithRaycast(point, start, end []float64) bool { +func intersectsWithRaycast(pX, pY, sX, sY, eX, eY float64) bool { // Always ensure that the the first point // has a y coordinate that is less than the second point - if start[1] > end[1] { + if sY > eY { // Switch the points if otherwise. - start, end = end, start + sX, sY, eX, eY = eX, eY, sX, sY } // Move the point's y coordinate // outside of the bounds of the testing region // so we can start drawing a ray - for point[1] == start[1] || point[1] == end[1] { - y := math.Nextafter(point[1], math.Inf(1)) - point = []float64{point[0], y} + for pY == sY || pY == eY { + pY = math.Nextafter(pY, math.Inf(1)) } // If we are outside of the polygon, indicate so. - if point[1] < start[1] || point[1] > end[1] { + if pY < sY || pY > eY { return false } - if start[0] > end[0] { - if point[0] > start[0] { + if sX > eX { + if pX > sX { return false } - if point[0] < end[0] { + if pX < eX { return true } } else { - if point[0] > end[0] { + if pX > eX { return false } - if point[0] < start[0] { + if pX < sX { return true } } - raySlope := (point[1] - start[1]) / (point[0] - start[0]) - diagSlope := (end[1] - start[1]) / (end[0] - start[0]) + raySlope := (pY - sY) / (pX - sX) + diagSlope := (eY - sY) / (eX - sX) return raySlope >= diagSlope }
--- a/pkg/octree/vertex.go Fri Oct 18 10:45:25 2019 +0200 +++ b/pkg/octree/vertex.go Fri Oct 18 11:20:33 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