Mercurial > gemma

comparison pkg/octree/simplify.go @ 3775:33fa76994b8a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fixed simplification a bit.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Tue, 02 Jul 2019 11:43:12 +0200
parents 545304d3ff93
children
comparison
equal deleted inserted replaced
3774:bb62c98fcf05 3775:33fa76994b8a
15 15
16 import ( 16 import (
17 "math" 17 "math"
18 ) 18 )
19 19
20 func (points MultiPointZ) Simplify(tolerance float64) MultiPointZ { 20 func handleTriangle(
21 21 t *Triangle,
22 if len(points) < 2 { 22 maxDist, tolerance float64,
23 return points 23 maxIdx int,
24 } 24 points MultiPointZ,
25 25 result *MultiPointZ,
26 if tolerance < 0 { 26 ) bool {
27 tolerance = -tolerance 27 if maxDist <= tolerance {
28 } 28 return false
29 29 }
30 min := Vertex{X: math.MaxFloat64, Y: math.MaxFloat64, Z: math.MaxFloat64} 30
31 max := Vertex{X: -math.MaxFloat64, Y: -math.MaxFloat64, Z: -math.MaxFloat64} 31 if len(points) == 1 {
32 32 *result = append(*result, points[0])
33 var maxIdx int 33 return true
34 34 }
35 for i, v := range points { 35
36 min.Minimize(v) 36 var (
37 37 tris [3]Triangle
38 if v.X < min.X { 38 planes [3]Plane3D
39 min.X = v.X 39 maxDists [3]float64
40 } 40 maxIdxs [3]int
41 if v.X > max.X { 41 parts [3]MultiPointZ
42 max.X = v.X 42 )
43 }
44 if v.Y < min.Y {
45 min.Y = v.Y
46 }
47 if v.Y > max.Y {
48 max.Y = v.Y
49 }
50 if v.Z < min.Z {
51 min.Z = v.Z
52 }
53 if v.Z > max.Z {
54 max.Z = v.Z
55 maxIdx = i
56 }
57
58 max.Maximize(v)
59 }
60
61 /*
62 log.Printf("(%.5f, %.5f, %.5f) - (%.5f, %.5f, %.5f)\n",
63 min.X, min.Y, min.Z,
64 max.X, max.Y, max.Z)
65 */
66
67 below := min.Z - 3*tolerance
68 xMin := min.X - tolerance
69 xMax := max.X + tolerance
70 yMin := min.Y - tolerance
71 yMax := max.Y + tolerance
72
73 corners := []Vertex{
74 {xMin, yMin, below},
75 {xMax, yMin, below},
76 {xMax, yMax, below},
77 {xMin, yMax, below},
78 }
79 43
80 top := points[maxIdx] 44 top := points[maxIdx]
81 45 for i := 0; i < 3; i++ {
82 tris := make([]Triangle, len(corners)) 46 tris[i] = Triangle{t[i], t[(i+1)%3], top}
83 planes := make([]Plane3D, len(corners))
84
85 for i, v1 := range corners {
86 v2 := corners[(i+1)%len(corners)]
87 tris[i] = Triangle{v1, v2, top}
88 planes[i] = tris[i].Plane3D() 47 planes[i] = tris[i].Plane3D()
89 } 48 }
90
91 parts := make([][]Vertex, len(tris))
92
93 maxDists := make([]float64, len(planes))
94 maxIdxs := make([]int, len(planes))
95 49
96 nextPoint: 50 nextPoint:
97 for i, v := range points { 51 for i, v := range points {
98 if i == maxIdx { 52 if i == maxIdx {
99 continue 53 continue
109 continue nextPoint 63 continue nextPoint
110 } 64 }
111 } 65 }
112 } 66 }
113 67
114 result := make([]Vertex, 0, len(points))
115
116 var handleTriangle func(*Triangle, float64, int, []Vertex) bool
117
118 handleTriangle = func(
119 t *Triangle,
120 maxDist float64, maxIdx int,
121 points []Vertex,
122 ) bool {
123 if maxDist <= tolerance {
124 return false
125 }
126
127 if len(points) == 1 {
128 result = append(result, points[0])
129 }
130
131 var (
132 tris [3]Triangle
133 planes [3]Plane3D
134 maxDists [3]float64
135 maxIdxs [3]int
136 parts [3][]Vertex
137 )
138
139 top := points[maxIdx]
140 for i := 0; i < 3; i++ {
141 tris[i] = Triangle{t[i], t[(i+1)%3], top}
142 planes[i] = tris[i].Plane3D()
143 }
144
145 nextPoint:
146 for i, v := range points {
147 if i == maxIdx {
148 continue
149 }
150
151 for j := range tris {
152 if tris[j].Contains(v.X, v.Y) {
153 if dist := math.Abs(planes[j].Eval(v)); dist > maxDists[j] {
154 maxDists[j] = dist
155 maxIdxs[j] = len(parts[j])
156 }
157 parts[j] = append(parts[j], v)
158 continue nextPoint
159 }
160 }
161 }
162
163 var found bool
164 for i, part := range parts {
165 if len(part) > 0 && handleTriangle(&tris[i], maxDists[i], maxIdxs[i], part) {
166 found = true
167 }
168 }
169
170 if found {
171 result = append(result, top)
172 }
173
174 return found
175 }
176
177 var found bool 68 var found bool
178 for i, part := range parts { 69 for i, part := range parts {
179 if len(part) > 0 && handleTriangle(&tris[i], maxDists[i], maxIdxs[i], part) { 70 if len(part) > 0 && handleTriangle(
71 &tris[i],
72 maxDists[i], tolerance,
73 maxIdxs[i],
74 part,
75 result,
76 ) {
77 found = true
78 }
79 }
80
81 if found {
82 *result = append(*result, top)
83 }
84
85 return found
86 }
87
88 func (points MultiPointZ) Simplify(tolerance float64) MultiPointZ {
89
90 if len(points) < 2 {
91 return points
92 }
93
94 if tolerance < 0 {
95 tolerance = -tolerance
96 }
97
98 min := Vertex{X: math.MaxFloat64, Y: math.MaxFloat64, Z: math.MaxFloat64}
99 max := Vertex{X: -math.MaxFloat64, Y: -math.MaxFloat64, Z: -math.MaxFloat64}
100
101 var maxIdx int
102
103 for i, v := range points {
104 min.Minimize(v)
105
106 if v.X < min.X {
107 min.X = v.X
108 }
109 if v.X > max.X {
110 max.X = v.X
111 }
112 if v.Y < min.Y {
113 min.Y = v.Y
114 }
115 if v.Y > max.Y {
116 max.Y = v.Y
117 }
118 if v.Z < min.Z {
119 min.Z = v.Z
120 }
121 if v.Z > max.Z {
122 max.Z = v.Z
123 maxIdx = i
124 }
125
126 max.Maximize(v)
127 }
128
129 /*
130 log.Printf("(%.5f, %.5f, %.5f) - (%.5f, %.5f, %.5f)\n",
131 min.X, min.Y, min.Z,
132 max.X, max.Y, max.Z)
133 */
134
135 below := min.Z - 3*tolerance
136 xMin := min.X - tolerance
137 xMax := max.X + tolerance
138 yMin := min.Y - tolerance
139 yMax := max.Y + tolerance
140
141 corners := []Vertex{
142 {xMin, yMin, below},
143 {xMax, yMin, below},
144 {xMax, yMax, below},
145 {xMin, yMax, below},
146 }
147
148 top := points[maxIdx]
149
150 tris := make([]Triangle, len(corners))
151 planes := make([]Plane3D, len(corners))
152
153 for i, v1 := range corners {
154 v2 := corners[(i+1)%len(corners)]
155 tris[i] = Triangle{v1, v2, top}
156 planes[i] = tris[i].Plane3D()
157 }
158
159 parts := make([][]Vertex, len(tris))
160
161 maxDists := make([]float64, len(planes))
162 maxIdxs := make([]int, len(planes))
163
164 nextPoint:
165 for i, v := range points {
166 if i == maxIdx {
167 continue
168 }
169
170 for j := range tris {
171 if tris[j].Contains(v.X, v.Y) {
172 if dist := math.Abs(planes[j].Eval(v)); dist > maxDists[j] {
173 maxDists[j] = dist
174 maxIdxs[j] = len(parts[j])
175 }
176 parts[j] = append(parts[j], v)
177 continue nextPoint
178 }
179 }
180 }
181
182 result := make(MultiPointZ, 0, len(points))
183
184 var found bool
185 for i, part := range parts {
186 if len(part) > 0 && handleTriangle(
187 &tris[i],
188 maxDists[i], tolerance,
189 maxIdxs[i],
190 part,
191 &result,
192 ) {
180 found = true 193 found = true
181 } 194 }
182 } 195 }
183 196
184 if found { 197 if found {