Mercurial > gemma

comparison cmd/srsimplify/main.go @ 3771:b7530ed07561 simplify-sounding-results

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Partition recursively. Does not produce good triangles.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Wed, 26 Jun 2019 12:06:28 +0200
parents 71164b817d6e
children 545304d3ff93
comparison
equal deleted inserted replaced
3770:71164b817d6e 3771:b7530ed07561
62 } 62 }
63 63
64 func storeXYZ(points octree.MultiPointZ, w io.Writer) error { 64 func storeXYZ(points octree.MultiPointZ, w io.Writer) error {
65 out := bufio.NewWriter(w) 65 out := bufio.NewWriter(w)
66 for i := range points { 66 for i := range points {
67 fmt.Fprintf(out, "%.5f %.5f %.5f\n", 67 fmt.Fprintf(out, "%.5f,%.5f,%.5f\n",
68 points[i].X, points[i].Y, points[i].Z) 68 points[i].X, points[i].Y, points[i].Z)
69 } 69 }
70 return out.Flush() 70 return out.Flush()
71 } 71 }
72 72
73 func simplify(points octree.MultiPointZ, tolerance float64) octree.MultiPointZ { 73 func simplify(points octree.MultiPointZ, tolerance float64) octree.MultiPointZ {
74 74
75 if len(points) < 2 { 75 if len(points) < 2 {
76 return points 76 return points
77 }
78
79 if tolerance < 0 {
80 tolerance = -tolerance
77 } 81 }
78 82
79 min := octree.Vertex{X: math.MaxFloat64, Y: math.MaxFloat64, Z: math.MaxFloat64} 83 min := octree.Vertex{X: math.MaxFloat64, Y: math.MaxFloat64, Z: math.MaxFloat64}
80 max := octree.Vertex{X: -math.MaxFloat64, Y: -math.MaxFloat64, Z: -math.MaxFloat64} 84 max := octree.Vertex{X: -math.MaxFloat64, Y: -math.MaxFloat64, Z: -math.MaxFloat64}
81 85
109 113
110 log.Printf("(%.5f, %.5f, %.5f) - (%.5f, %.5f, %.5f)\n", 114 log.Printf("(%.5f, %.5f, %.5f) - (%.5f, %.5f, %.5f)\n",
111 min.X, min.Y, min.Z, 115 min.X, min.Y, min.Z,
112 max.X, max.Y, max.Z) 116 max.X, max.Y, max.Z)
113 117
114 below := min.Z - 2*tolerance 118 below := min.Z - 3*tolerance
115 xMin := min.X - 1 119 xMin := min.X - tolerance
116 xMax := max.X + 1 120 xMax := max.X + tolerance
117 yMin := min.Y - 1 121 yMin := min.Y - tolerance
118 yMax := max.Y + 1 122 yMax := max.Y + tolerance
119 123
120 corners := []octree.Vertex{ 124 corners := []octree.Vertex{
121 {xMin, yMin, below}, 125 {xMin, yMin, below},
122 {xMax, yMin, below}, 126 {xMax, yMin, below},
123 {xMax, yMax, below}, 127 {xMax, yMax, below},
134 tris[i] = octree.Triangle{v1, v2, top} 138 tris[i] = octree.Triangle{v1, v2, top}
135 planes[i] = tris[i].Plane3D() 139 planes[i] = tris[i].Plane3D()
136 } 140 }
137 141
138 parts := make([][]octree.Vertex, len(tris)) 142 parts := make([][]octree.Vertex, len(tris))
139 var rest int
140 143
141 maxDists := make([]float64, len(planes)) 144 maxDists := make([]float64, len(planes))
142 maxIdxs := make([]int, len(planes)) 145 maxIdxs := make([]int, len(planes))
143 146
144 nextPoint: 147 nextPoint:
155 } 158 }
156 parts[j] = append(parts[j], v) 159 parts[j] = append(parts[j], v)
157 continue nextPoint 160 continue nextPoint
158 } 161 }
159 } 162 }
160 rest++ 163 }
161 } 164
162 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
163 for i, part := range parts { 229 for i, part := range parts {
164 log.Printf("%d: %d %.5f\n", i, len(part), maxDists[i]) 230 if len(part) > 0 && handleTriangle(&tris[i], maxDists[i], maxIdxs[i], part) {
165 } 231 found = true
166 log.Printf("rest: %d\n", rest) 232 }
167 233 }
168 // TODO: Implement me! 234
169 return points 235 if found {
236 result = append(result, top)
237 }
238
239 return result
170 } 240 }
171 241
172 func main() { 242 func main() {
173 243
174 var tolerance float64 244 var tolerance float64