Mercurial > gemma

changeset 2516:1ec4c5633eb6 octree-diff

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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