Mercurial > gemma
changeset 2516:1ec4c5633eb6 octree-diff
Clip STR tree and not Octree.
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Tue, 05 Mar 2019 17:08:16 +0100 |
parents | 6bcaa8bf2603 |
children | ac61c250337b |
files | cmd/octreediff/main.go pkg/imports/sr.go pkg/octree/builder.go pkg/octree/strtree.go pkg/octree/tree.go |
diffstat | 5 files changed, 97 insertions(+), 346 deletions(-) [+] |
line wrap: on
line diff
--- a/cmd/octreediff/main.go Tue Mar 05 15:55:14 2019 +0100 +++ b/cmd/octreediff/main.go Tue Mar 05 17:08:16 2019 +0100 @@ -395,26 +395,29 @@ var str octree.STRTree - str.Build(tri) + tin := tri.Tin() + + str.Build(tin) now = time.Now() log.Printf("building STR tree took %v\n", now.Sub(last)) last = now - builder := octree.NewBuilder(tri.Tin()) - builder.Build() + removed := str.Clip(&clippingPolygon) + now = time.Now() + log.Printf("clipping STR tree took %v\n", now.Sub(last)) + last = now + + log.Printf("Number of triangles to clip: %d\n", len(removed)) + + builder := octree.NewBuilder(tin) + builder.Build(removed) now = time.Now() log.Printf("building octree took %v\n", now.Sub(last)) last = now tree := builder.Tree() - tree.Clip(&clippingPolygon) - //tree.Dump() - - now = time.Now() - log.Printf("clipping octree took %v\n", now.Sub(last)) - last = now log.Printf("min/max: %f %f\n", tree.Min.Z, tree.Max.Z)
--- a/pkg/imports/sr.go Tue Mar 05 15:55:14 2019 +0100 +++ b/pkg/imports/sr.go Tue Mar 05 17:08:16 2019 +0100 @@ -265,7 +265,7 @@ feedback.Info("Building octree...") builder := octree.NewBuilder(tin) start = time.Now() - builder.Build() + builder.Build(nil) octreeIndex, err := builder.Bytes() tin = nil // not needed from now on feedback.Info("building octree took %s", time.Since(start))
--- a/pkg/octree/builder.go Tue Mar 05 15:55:14 2019 +0100 +++ b/pkg/octree/builder.go Tue Mar 05 17:08:16 2019 +0100 @@ -71,11 +71,24 @@ } // Build builds the Octree. -func (tb *Builder) Build() { +func (tb *Builder) Build(removed map[int32]struct{}) { + + var triangles []int32 - triangles := make([]int32, len(tb.t.Triangles)) - for i := range triangles { - triangles[i] = int32(i) + if len(removed) > 0 { + triangles = make([]int32, len(tb.t.Triangles)-len(removed)) + idx := 0 + for i := range tb.t.Triangles { + if _, found := removed[int32(i)]; !found { + triangles[idx] = int32(i) + idx++ + } + } + } else { + triangles = make([]int32, len(tb.t.Triangles)) + for i := range triangles { + triangles[i] = int32(i) + } } n := runtime.NumCPU()
--- a/pkg/octree/strtree.go Tue Mar 05 15:55:14 2019 +0100 +++ b/pkg/octree/strtree.go Tue Mar 05 17:08:16 2019 +0100 @@ -21,19 +21,16 @@ const numEntries = 8 type STRTree struct { - tri *Triangulation - index []int32 - triangles [][]int32 - bboxes []Box2D + tin *Tin + index []int32 + bboxes []Box2D } -func (s *STRTree) Build(t *Triangulation) { - - s.tri = t +func (s *STRTree) Build(t *Tin) { - s.triangles = t.TriangleSlices() + s.tin = t - all := make([]int32, len(s.triangles)) + all := make([]int32, len(t.Triangles)) for i := range all { all[i] = int32(i) @@ -46,11 +43,63 @@ s.index[0] = root } +func (s *STRTree) Clip(p *Polygon) map[int32]struct{} { + + removed := make(map[int32]struct{}) + + stack := []int32{s.index[0]} + + for len(stack) > 0 { + top := stack[len(stack)-1] + stack = stack[:len(stack)-1] + + if top > 0 { // node + switch p.IntersectionBox2D(s.bbox(top)) { + case IntersectionInside: + // all triangles are inside polygon + case IntersectionOutSide: + // all triangles are outside polygon + s.allTriangles(top, removed) + default: + // mixed bag + for i, n := int32(0), s.index[top+1]; i < n; i++ { + stack = append(stack, s.index[top+2+i]) + } + } + } else { // leaf + top = -top - 1 + for i, n := int32(0), s.index[top+1]; i < n; i++ { + removed[s.index[top+2+i]] = struct{}{} + } + } + } + + return removed +} + +func (s *STRTree) allTriangles(pos int32, tris map[int32]struct{}) { + stack := []int32{pos} + for len(stack) > 0 { + top := stack[len(stack)-1] + stack = stack[:len(stack)-1] + if top > 0 { // node + for i, n := int32(0), s.index[top+1]; i < n; i++ { + stack = append(stack, s.index[top+2+i]) + } + } else { // leaf + top = -top - 1 + for i, n := int32(0), s.index[top+1]; i < n; i++ { + tris[s.index[top+2+i]] = struct{}{} + } + } + } +} + func (s *STRTree) build(items []int32) int32 { sort.Slice(items, func(i, j int) bool { - ti := s.triangles[items[i]] - tj := s.triangles[items[j]] - return s.tri.Points[ti[0]].X < s.tri.Points[tj[0]].X + ti := s.tin.Triangles[items[i]] + tj := s.tin.Triangles[items[j]] + return s.tin.Vertices[ti[0]].X < s.tin.Vertices[tj[0]].X }) P := int(math.Ceil(float64(len(items)) / float64(numEntries))) @@ -75,9 +124,9 @@ for _, slice := range slices { sort.Slice(slice, func(i, j int) bool { - ti := s.triangles[slice[i]] - tj := s.triangles[slice[j]] - return s.tri.Points[ti[0]].Y < s.tri.Points[tj[0]].Y + ti := s.tin.Triangles[slice[i]] + tj := s.tin.Triangles[slice[j]] + return s.tin.Vertices[ti[0]].Y < s.tin.Vertices[tj[0]].Y }) for len(slice) > 0 { @@ -174,8 +223,8 @@ s.index = append(s.index, 0, int32(len(items))) s.index = append(s.index, items...) if len(items) > 0 { - vertices := s.tri.Points - ti := s.triangles[items[0]] + vertices := s.tin.Vertices + ti := s.tin.Triangles[items[0]] t := Triangle{ vertices[ti[0]], vertices[ti[1]], @@ -184,7 +233,7 @@ box := t.BBox() for i := 1; i < len(items); i++ { it := items[i] - ti := s.triangles[it] + ti := s.tin.Triangles[it] t := Triangle{ vertices[ti[0]], vertices[ti[1]],
--- a/pkg/octree/tree.go Tue Mar 05 15:55:14 2019 +0100 +++ b/pkg/octree/tree.go Tue Mar 05 17:08:16 2019 +0100 @@ -14,10 +14,7 @@ package octree import ( - "fmt" - "log" "math" - "os" ) // Tree is an Octree holding triangles. @@ -51,317 +48,6 @@ {0.5, 0.5, 1.0, 1.0}, } -func (ot *Tree) Dump() error { - fake, err := os.Create("/tmp/dump.geojson") - if err != nil { - return err - } - defer fake.Close() - - fmt.Fprintln(fake, - `{ -"crs": { - "type": "name", - "properties": { - "name": "urn:ogc:def:crs:EPSG::32633" - } - }, - - "type": "FeatureCollection", - "features": [`) - - first := true - - already := make(map[int32]struct{}) - - stack := []int32{1} - for len(stack) > 0 { - pos := stack[len(stack)-1] - stack = stack[:len(stack)-1] - - if pos > 0 { // node - index := ot.index[pos:] - if len(index) > 8 { - index = index[:8] - } - for i := range index { - if index[i] != 0 { - stack = append(stack, index[i]) - } - } - } else { - pos = -pos - 1 - n := ot.index[pos] - indices := ot.index[pos+1 : pos+1+n] - for _, index := range indices { - if _, found := already[index]; found { - continue - } - already[index] = struct{}{} - tri := ot.triangles[index] - t := Triangle{ - ot.vertices[tri[0]], - ot.vertices[tri[1]], - ot.vertices[tri[2]], - } - if first { - first = false - } else { - fmt.Fprintln(fake, ",") - } - fmt.Fprintf(fake, - `{ "type": "Feature", - "properties": { - "id": %d - }, - "geometry": { - "type": "Polygon", - "coordinates": [[ - [ %.6f, %.6f ], - [ %.6f, %.6f ], - [ %.6f, %.6f ], - [ %.6f, %.6f ]] - ] - } -}`, index, - t[0].X, t[0].Y, - t[1].X, t[1].Y, - t[2].X, t[2].Y, - t[0].X, t[0].Y) - } - } - } - - fmt.Fprintln(fake, - `] -}`) - return nil -} - -func (ot *Tree) eliminate(removed map[int32]IntersectionType) { - - var bad int - - stack := []int32{1} - for len(stack) > 0 { - pos := stack[len(stack)-1] - stack = stack[:len(stack)-1] - - if pos > 0 { // node - index := ot.index[pos:] - if len(index) > 8 { - index = index[:8] - } - for i := range index { - if index[i] != 0 { - stack = append(stack, index[i]) - } - } - } else { - pos = -pos - 1 - n := ot.index[pos] - indices := ot.index[pos+1 : pos+1+n] - - for i := len(indices) - 1; i >= 0; i-- { - index := indices[i] - if t, found := removed[index]; found && t != IntersectionInside { - if i < len(indices)-1 { - copy(indices[i:], indices[i+1:]) - } - indices[len(indices)-1] = 0 - indices = indices[:len(indices)-1] - bad++ - } - } - ot.index[pos] = int32(len(indices)) - } - } - - log.Printf("bad triangles: %d\n", bad) -} - -func (ot *Tree) Clip(p *Polygon) { - - /* - fake, _ := os.Create("/tmp/removed.geojson") - defer fake.Close() - - fmt.Fprintln(fake, - `{ - "crs": { - "type": "name", - "properties": { - "name": "urn:ogc:def:crs:EPSG::32633" - } - }, - - "type": "FeatureCollection", - "features": [`) - - first := true - */ - - log.Printf("info: num triangles: %d\n", len(ot.triangles)) - - all := Box2D{ot.Min.X, ot.Min.Y, ot.Max.X, ot.Max.Y} - - stack := []boxFrame{{1, all}} - - triChecks := make(map[int32]IntersectionType) - - var nodeTests int - var nodesClipped int - var trianglesClipped int - var nodesAllInside int - -frames: - for len(stack) > 0 { - top := stack[len(stack)-1] - stack = stack[:len(stack)-1] - - if top.pos > 0 { // node - nodeTests++ - switch p.IntersectionBox2D(top.Box2D) { - case IntersectionInside: - // all inside so nothing to clip. - nodesAllInside++ - continue frames - case IntersectionOutSide: - // all outside -> clip from tree. - index := ot.index[top.pos:] - if len(index) > 8 { - index = index[:8] - } - for i := range index { - if index[i] != 0 { - index[i] = 0 - nodesClipped++ - } - } - 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. - nodesAllInside++ - continue children - case IntersectionOutSide: - // all are ouside -> clip from tree. - if a != 0 { - index[i] = 0 - nodesClipped++ - } - if b != 0 { - index[i+4] = 0 - nodesClipped++ - } - continue children - 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 - - /* - - if what != IntersectionInside { - if first { - first = false - } else { - fmt.Fprintln(fake, ",") - } - fmt.Fprintf(fake, - `{ "type": "Feature", - "properties": { - "id": %d - }, - "geometry": { - "type": "Polygon", - "coordinates": [[ - [ %.6f, %.6f ], - [ %.6f, %.6f ], - [ %.6f, %.6f ], - [ %.6f, %.6f ]] - ] - } - }`, triIndex, - t[0].X, t[0].Y, - t[1].X, t[1].Y, - t[2].X, t[2].Y, - t[0].X, t[0].Y) - } - */ - } - switch what { - case IntersectionInside: - // triangle inside -> stay. - continue tris - default: - trianglesClipped++ - // 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)) - } - } - /* - fmt.Fprintln(fake, - `] - }`) - */ - log.Printf("info: node tests: %d\n", nodeTests) - log.Printf("info: nodes clipped: %d\n", nodesClipped) - log.Printf("info: nodes all inside: %d\n", nodesAllInside) - log.Printf("info: triangle clipped: %d\n", trianglesClipped) - log.Printf("info: tested triangles: %d\n", len(triChecks)) - - ot.eliminate(triChecks) -} - func (ot *Tree) Value(x, y float64) (float64, bool) { // out of bounding box