Mercurial > gemma

comparison cmd/srsimplify/main.go @ 3772:545304d3ff93 simplify-sounding-results

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Moved simplification algorithm to octree package.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Mon, 01 Jul 2019 11:08:19 +0200
parents b7530ed07561
children
comparison
equal deleted inserted replaced
3771:b7530ed07561 3772:545304d3ff93
17 "bufio" 17 "bufio"
18 "flag" 18 "flag"
19 "fmt" 19 "fmt"
20 "io" 20 "io"
21 "log" 21 "log"
22 "math"
23 "os" 22 "os"
24 "strconv" 23 "strconv"
25 "strings" 24 "strings"
26 25
27 "gemma.intevation.de/gemma/pkg/octree" 26 "gemma.intevation.de/gemma/pkg/octree"
68 points[i].X, points[i].Y, points[i].Z) 67 points[i].X, points[i].Y, points[i].Z)
69 } 68 }
70 return out.Flush() 69 return out.Flush()
71 } 70 }
72 71
73 func simplify(points octree.MultiPointZ, tolerance float64) octree.MultiPointZ {
74
75 if len(points) < 2 {
76 return points
77 }
78
79 if tolerance < 0 {
80 tolerance = -tolerance
81 }
82
83 min := octree.Vertex{X: math.MaxFloat64, Y: math.MaxFloat64, Z: math.MaxFloat64}
84 max := octree.Vertex{X: -math.MaxFloat64, Y: -math.MaxFloat64, Z: -math.MaxFloat64}
85
86 var maxIdx int
87
88 for i, v := range points {
89 min.Minimize(v)
90
91 if v.X < min.X {
92 min.X = v.X
93 }
94 if v.X > max.X {
95 max.X = v.X
96 }
97 if v.Y < min.Y {
98 min.Y = v.Y
99 }
100 if v.Y > max.Y {
101 max.Y = v.Y
102 }
103 if v.Z < min.Z {
104 min.Z = v.Z
105 }
106 if v.Z > max.Z {
107 max.Z = v.Z
108 maxIdx = i
109 }
110
111 max.Maximize(v)
112 }
113
114 log.Printf("(%.5f, %.5f, %.5f) - (%.5f, %.5f, %.5f)\n",
115 min.X, min.Y, min.Z,
116 max.X, max.Y, max.Z)
117
118 below := min.Z - 3*tolerance
119 xMin := min.X - tolerance
120 xMax := max.X + tolerance
121 yMin := min.Y - tolerance
122 yMax := max.Y + tolerance
123
124 corners := []octree.Vertex{
125 {xMin, yMin, below},
126 {xMax, yMin, below},
127 {xMax, yMax, below},
128 {xMin, yMax, below},
129 }
130
131 top := points[maxIdx]
132
133 tris := make([]octree.Triangle, len(corners))
134 planes := make([]octree.Plane3D, len(corners))
135
136 for i, v1 := range corners {
137 v2 := corners[(i+1)%len(corners)]
138 tris[i] = octree.Triangle{v1, v2, top}
139 planes[i] = tris[i].Plane3D()
140 }
141
142 parts := make([][]octree.Vertex, len(tris))
143
144 maxDists := make([]float64, len(planes))
145 maxIdxs := make([]int, len(planes))
146
147 nextPoint:
148 for i, v := range points {
149 if i == maxIdx {
150 continue
151 }
152
153 for j := range tris {
154 if tris[j].Contains(v.X, v.Y) {
155 if dist := math.Abs(planes[j].Eval(v)); dist > maxDists[j] {
156 maxDists[j] = dist
157 maxIdxs[j] = len(parts[j])
158 }
159 parts[j] = append(parts[j], v)
160 continue nextPoint
161 }
162 }
163 }
164
165 result := make([]octree.Vertex, 0, len(points))
166
167 var handleTriangle func(*octree.Triangle, float64, int, []octree.Vertex) bool
168
169 handleTriangle = func(
170 t *octree.Triangle,
171 maxDist float64, maxIdx int,
172 points []octree.Vertex,
173 ) bool {
174 if maxDist <= tolerance {
175 return false
176 }
177
178 if len(points) == 1 {
179 result = append(result, points[0])
180 }
181
182 var (
183 tris [3]octree.Triangle
184 planes [3]octree.Plane3D
185 maxDists [3]float64
186 maxIdxs [3]int
187 parts [3][]octree.Vertex
188 )
189
190 top := points[maxIdx]
191 for i := 0; i < 3; i++ {
192 tris[i] = octree.Triangle{t[i], t[(i+1)%3], top}
193 planes[i] = tris[i].Plane3D()
194 }
195
196 nextPoint:
197 for i, v := range points {
198 if i == maxIdx {
199 continue
200 }
201
202 for j := range tris {
203 if tris[j].Contains(v.X, v.Y) {
204 if dist := math.Abs(planes[j].Eval(v)); dist > maxDists[j] {
205 maxDists[j] = dist
206 maxIdxs[j] = len(parts[j])
207 }
208 parts[j] = append(parts[j], v)
209 continue nextPoint
210 }
211 }
212 }
213
214 var found bool
215 for i, part := range parts {
216 if len(part) > 0 && handleTriangle(&tris[i], maxDists[i], maxIdxs[i], part) {
217 found = true
218 }
219 }
220
221 if found {
222 result = append(result, top)
223 }
224
225 return found
226 }
227
228 var found bool
229 for i, part := range parts {
230 if len(part) > 0 && handleTriangle(&tris[i], maxDists[i], maxIdxs[i], part) {
231 found = true
232 }
233 }
234
235 if found {
236 result = append(result, top)
237 }
238
239 return result
240 }
241
242 func main() { 72 func main() {
243 73
244 var tolerance float64 74 var tolerance float64
245 75
246 flag.Float64Var(&tolerance, "t", 0.1, "accepted tolerance (shorthand)") 76 flag.Float64Var(&tolerance, "t", 0.1, "accepted tolerance (shorthand)")
251 points, err := loadXYZ(os.Stdin) 81 points, err := loadXYZ(os.Stdin)
252 if err != nil { 82 if err != nil {
253 log.Fatalf("err: %v\n", err) 83 log.Fatalf("err: %v\n", err)
254 } 84 }
255 85
256 points = simplify(points, tolerance) 86 points = points.Simplify(tolerance)
257 87
258 if err := storeXYZ(points, os.Stdout); err != nil { 88 if err := storeXYZ(points, os.Stdout); err != nil {
259 log.Fatalf("err: %v\n", err) 89 log.Fatalf("err: %v\n", err)
260 } 90 }
261 } 91 }