Mercurial > gemma

comparison pkg/octree/tree.go @ 2466:a1e751c08c56 octree-diff

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Calculate difference on single core.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Mon, 25 Feb 2019 18:21:33 +0100
parents fe1aa62195c2
children c85b16db8a02
comparison
equal deleted inserted replaced
2465:86c7a023400e 2466:a1e751c08c56
30 Min Vertex 30 Min Vertex
31 // Max is the upper right corner of the bbox. 31 // Max is the upper right corner of the bbox.
32 Max Vertex 32 Max Vertex
33 } 33 }
34 34
35 func (ot *Tree) Vertices() []Vertex {
36 return ot.vertices
37 }
38
35 var scale = [4][4]float64{ 39 var scale = [4][4]float64{
36 {0.0, 0.0, 0.5, 0.5}, 40 {0.0, 0.0, 0.5, 0.5},
37 {0.5, 0.0, 1.0, 0.5}, 41 {0.5, 0.0, 1.0, 0.5},
38 {0.0, 0.5, 0.5, 1.0}, 42 {0.0, 0.5, 0.5, 1.0},
39 {0.5, 0.5, 1.0, 1.0}, 43 {0.5, 0.5, 1.0, 1.0},
40 } 44 }
41 45
42 // Vertical does a vertical cross cut from (x1, y1) to (x2, y2). 46 func (ot *Tree) Value(x, y float64) (float64, bool) {
43 func (ot *Tree) Vertical(x1, y1, x2, y2 float64, fn func(*Triangle)) {
44
45 box := Box2D{
46 X1: math.Min(x1, x2),
47 Y1: math.Min(y1, y2),
48 X2: math.Max(x1, x2),
49 Y2: math.Max(y1, y2),
50 }
51 47
52 // out of bounding box 48 // out of bounding box
53 if box.X2 < ot.Min.X || ot.Max.X < box.X1 || 49 if x < ot.Min.X || ot.Max.X < x ||
54 box.Y2 < ot.Min.Y || ot.Max.Y < box.Y1 { 50 y < ot.Min.Y || ot.Max.Y < y {
55 return 51 return 0, false
56 } 52 }
57
58 line := NewPlane2D(x1, y1, x2, y2)
59 53
60 type frame struct { 54 type frame struct {
61 pos int32 55 pos int32
62 Box2D 56 Box2D
63 } 57 }
64 58
65 dupes := map[int32]struct{}{}
66
67 all := Box2D{ot.Min.X, ot.Min.Y, ot.Max.X, ot.Max.Y} 59 all := Box2D{ot.Min.X, ot.Min.Y, ot.Max.X, ot.Max.Y}
68 //log.Printf("area: %f\n", (box.x2-box.x1)*(box.y2-box.y1))
69 //log.Printf("all: %f\n", (all.x2-all.x1)*(all.y2-all.y1))
70 60
71 stack := []frame{{1, all}} 61 stack := []frame{{1, all}}
72 62
73 for len(stack) > 0 { 63 for len(stack) > 0 {
74 top := stack[len(stack)-1] 64 top := stack[len(stack)-1]
88 dx*scale[i][0] + top.X1, 78 dx*scale[i][0] + top.X1,
89 dy*scale[i][1] + top.Y1, 79 dy*scale[i][1] + top.Y1,
90 dx*scale[i][2] + top.X1, 80 dx*scale[i][2] + top.X1,
91 dy*scale[i][3] + top.Y1, 81 dy*scale[i][3] + top.Y1,
92 } 82 }
83 if nbox.Contains(x, y) {
84 if a != 0 {
85 stack = append(stack, frame{a, nbox})
86 }
87 if b != 0 {
88 stack = append(stack, frame{b, nbox})
89 }
90 break
91 }
92 }
93 }
94 } else { // leaf
95 pos := -top.pos - 1
96 n := ot.index[pos]
97 indices := ot.index[pos+1 : pos+1+n]
98
99 for _, idx := range indices {
100 tri := ot.triangles[idx]
101 t := Triangle{
102 ot.vertices[tri[0]],
103 ot.vertices[tri[1]],
104 ot.vertices[tri[2]],
105 }
106 if t.Contains(x, y) {
107 return t.Plane3D().Z(x, y), true
108 }
109 }
110 }
111 }
112
113 return 0, false
114 }
115
116 // Vertical does a vertical cross cut from (x1, y1) to (x2, y2).
117 func (ot *Tree) Vertical(x1, y1, x2, y2 float64, fn func(*Triangle)) {
118
119 box := Box2D{
120 X1: math.Min(x1, x2),
121 Y1: math.Min(y1, y2),
122 X2: math.Max(x1, x2),
123 Y2: math.Max(y1, y2),
124 }
125
126 // out of bounding box
127 if box.X2 < ot.Min.X || ot.Max.X < box.X1 ||
128 box.Y2 < ot.Min.Y || ot.Max.Y < box.Y1 {
129 return
130 }
131
132 line := NewPlane2D(x1, y1, x2, y2)
133
134 type frame struct {
135 pos int32
136 Box2D
137 }
138
139 dupes := map[int32]struct{}{}
140
141 all := Box2D{ot.Min.X, ot.Min.Y, ot.Max.X, ot.Max.Y}
142 //log.Printf("area: %f\n", (box.x2-box.x1)*(box.y2-box.y1))
143 //log.Printf("all: %f\n", (all.x2-all.x1)*(all.y2-all.y1))
144
145 stack := []frame{{1, all}}
146
147 for len(stack) > 0 {
148 top := stack[len(stack)-1]
149 stack = stack[:len(stack)-1]
150
151 if top.pos > 0 { // node
152 if index := ot.index[top.pos:]; len(index) > 7 {
153 for i := 0; i < 4; i++ {
154 a := index[i]
155 b := index[i+4]
156 if a == 0 && b == 0 {
157 continue
158 }
159 dx := top.X2 - top.X1
160 dy := top.Y2 - top.Y1
161 nbox := Box2D{
162 dx*scale[i][0] + top.X1,
163 dy*scale[i][1] + top.Y1,
164 dx*scale[i][2] + top.X1,
165 dy*scale[i][3] + top.Y1,
166 }
93 if nbox.Intersects(box) && nbox.IntersectsPlane(line) { 167 if nbox.Intersects(box) && nbox.IntersectsPlane(line) {
94 if a != 0 { 168 if a != 0 {
95 stack = append(stack, frame{a, nbox}) 169 stack = append(stack, frame{a, nbox})
96 } 170 }
97 if b != 0 { 171 if b != 0 {