Mercurial > gemma
comparison pkg/mesh/triangulator.go @ 4827:f4abfd0ee8ad remove-octree-debris
Renamed octree package to mesh as there is no octree any more.
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Tue, 05 Nov 2019 14:30:22 +0100 |
parents | pkg/octree/triangulator.go@d12c2f4d3483 |
children |
comparison
equal
deleted
inserted
replaced
4826:ec5afada70ec | 4827:f4abfd0ee8ad |
---|---|
1 // Copyright (C) 2018 Michael Fogleman | |
2 // | |
3 // Permission is hereby granted, free of charge, to any person obtaining | |
4 // a copy of this software and associated documentation files (the "Software"), | |
5 // to deal in the Software without restriction, including without limitation | |
6 // the rights to use, copy, modify, merge, publish, distribute, sublicense, | |
7 // and/or sell copies of the Software, and to permit persons to whom the | |
8 // Software is furnished to do so, subject to the following conditions: | |
9 // | |
10 // The above copyright notice and this permission notice shall be included | |
11 // in all copies or substantial portions of the Software. | |
12 // | |
13 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS | |
14 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, | |
15 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL | |
16 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER | |
17 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, | |
18 // OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. | |
19 | |
20 package mesh | |
21 | |
22 import ( | |
23 "errors" | |
24 "math" | |
25 "sort" | |
26 ) | |
27 | |
28 type triangulator struct { | |
29 points []Vertex | |
30 squaredDistances []float64 | |
31 ids []int32 | |
32 center Vertex | |
33 triangles []int32 | |
34 halfedges []int32 | |
35 trianglesLen int | |
36 hull *node | |
37 hash []*node | |
38 } | |
39 | |
40 func newTriangulator(points []Vertex) *triangulator { | |
41 return &triangulator{points: points} | |
42 } | |
43 | |
44 // sorting a triangulator sorts the `ids` such that the referenced points | |
45 // are in order by their distance to `center` | |
46 func (tri *triangulator) Len() int { | |
47 return len(tri.points) | |
48 } | |
49 | |
50 func (tri *triangulator) Swap(i, j int) { | |
51 tri.ids[i], tri.ids[j] = tri.ids[j], tri.ids[i] | |
52 } | |
53 | |
54 func (tri *triangulator) Less(i, j int) bool { | |
55 d1 := tri.squaredDistances[tri.ids[i]] | |
56 d2 := tri.squaredDistances[tri.ids[j]] | |
57 if d1 != d2 { | |
58 return d1 < d2 | |
59 } | |
60 p1 := tri.points[tri.ids[i]] | |
61 p2 := tri.points[tri.ids[j]] | |
62 if p1.X != p2.X { | |
63 return p1.X < p2.X | |
64 } | |
65 return p1.Y < p2.Y | |
66 } | |
67 | |
68 func (tri *triangulator) triangulate() error { | |
69 points := tri.points | |
70 | |
71 n := len(points) | |
72 if n == 0 { | |
73 return nil | |
74 } | |
75 | |
76 tri.ids = make([]int32, n) | |
77 | |
78 // compute bounds | |
79 x0 := points[0].X | |
80 y0 := points[0].Y | |
81 x1 := points[0].X | |
82 y1 := points[0].Y | |
83 for i, p := range points { | |
84 if p.X < x0 { | |
85 x0 = p.X | |
86 } | |
87 if p.X > x1 { | |
88 x1 = p.X | |
89 } | |
90 if p.Y < y0 { | |
91 y0 = p.Y | |
92 } | |
93 if p.Y > y1 { | |
94 y1 = p.Y | |
95 } | |
96 tri.ids[i] = int32(i) | |
97 } | |
98 | |
99 var i0, i1, i2 int32 | |
100 | |
101 // pick a seed point close to midpoint | |
102 m := Vertex{X: (x0 + x1) / 2, Y: (y0 + y1) / 2} | |
103 minDist := infinity | |
104 for i, p := range points { | |
105 d := p.SquaredDistance2D(m) | |
106 if d < minDist { | |
107 i0 = int32(i) | |
108 minDist = d | |
109 } | |
110 } | |
111 | |
112 // find point closest to seed point | |
113 minDist = infinity | |
114 for i, p := range points { | |
115 if int32(i) == i0 { | |
116 continue | |
117 } | |
118 d := p.SquaredDistance2D(points[i0]) | |
119 if d > 0 && d < minDist { | |
120 i1 = int32(i) | |
121 minDist = d | |
122 } | |
123 } | |
124 | |
125 // find the third point which forms the smallest circumcircle | |
126 minRadius := infinity | |
127 for i, p := range points { | |
128 if int32(i) == i0 || int32(i) == i1 { | |
129 continue | |
130 } | |
131 r := circumradius(points[i0], points[i1], p) | |
132 if r < minRadius { | |
133 i2 = int32(i) | |
134 minRadius = r | |
135 } | |
136 } | |
137 if minRadius == infinity { | |
138 return errors.New("no delaunay triangulation exists for this input") | |
139 } | |
140 | |
141 // swap the order of the seed points for counter-clockwise orientation | |
142 if area(points[i0], points[i1], points[i2]) < 0 { | |
143 i1, i2 = i2, i1 | |
144 } | |
145 | |
146 tri.center = circumcenter(points[i0], points[i1], points[i2]) | |
147 | |
148 // sort the points by distance from the seed triangle circumcenter | |
149 tri.squaredDistances = make([]float64, n) | |
150 for i, p := range points { | |
151 tri.squaredDistances[i] = p.SquaredDistance2D(tri.center) | |
152 } | |
153 sort.Sort(tri) | |
154 | |
155 // initialize a hash table for storing edges of the advancing convex hull | |
156 hashSize := int(math.Ceil(math.Sqrt(float64(n)))) | |
157 tri.hash = make([]*node, hashSize) | |
158 | |
159 // initialize a circular doubly-linked list that will hold an advancing convex hull | |
160 nodes := make([]node, n) | |
161 | |
162 e := newNode(nodes, i0, nil) | |
163 e.t = 0 | |
164 tri.hashEdge(e) | |
165 | |
166 e = newNode(nodes, i1, e) | |
167 e.t = 1 | |
168 tri.hashEdge(e) | |
169 | |
170 e = newNode(nodes, i2, e) | |
171 e.t = 2 | |
172 tri.hashEdge(e) | |
173 | |
174 tri.hull = e | |
175 | |
176 maxTriangles := 2*n - 5 | |
177 tri.triangles = make([]int32, maxTriangles*3) | |
178 tri.halfedges = make([]int32, maxTriangles*3) | |
179 | |
180 tri.addTriangle(i0, i1, i2, -1, -1, -1) | |
181 | |
182 pp := Vertex{X: infinity, Y: infinity} | |
183 for k := 0; k < n; k++ { | |
184 i := tri.ids[k] | |
185 p := points[i] | |
186 | |
187 // skip nearly-duplicate points | |
188 if p.SquaredDistance2D(pp) < eps { | |
189 continue | |
190 } | |
191 pp = p | |
192 | |
193 // skip seed triangle points | |
194 if i == int32(i0) || i == int32(i1) || i == int32(i2) { | |
195 continue | |
196 } | |
197 | |
198 // find a visible edge on the convex hull using edge hash | |
199 var start *node | |
200 key := tri.hashKey(p) | |
201 for j := 0; j < len(tri.hash); j++ { | |
202 start = tri.hash[key] | |
203 if start != nil && start.i >= 0 { | |
204 break | |
205 } | |
206 key++ | |
207 if key >= len(tri.hash) { | |
208 key = 0 | |
209 } | |
210 } | |
211 start = start.prev | |
212 | |
213 e := start | |
214 for area(p, points[e.i], points[e.next.i]) >= 0 { | |
215 e = e.next | |
216 if e == start { | |
217 e = nil | |
218 break | |
219 } | |
220 } | |
221 if e == nil { | |
222 // likely a near-duplicate point; skip it | |
223 continue | |
224 } | |
225 walkBack := e == start | |
226 | |
227 // add the first triangle from the point | |
228 t := tri.addTriangle(e.i, i, e.next.i, -1, -1, e.t) | |
229 e.t = t // keep track of boundary triangles on the hull | |
230 e = newNode(nodes, i, e) | |
231 | |
232 // recursively flip triangles from the point until they satisfy the Delaunay condition | |
233 e.t = tri.legalize(t + 2) | |
234 | |
235 // walk forward through the hull, adding more triangles and flipping recursively | |
236 q := e.next | |
237 for area(p, points[q.i], points[q.next.i]) < 0 { | |
238 t = tri.addTriangle(q.i, i, q.next.i, q.prev.t, -1, q.t) | |
239 q.prev.t = tri.legalize(t + 2) | |
240 tri.hull = q.remove() | |
241 q = q.next | |
242 } | |
243 | |
244 if walkBack { | |
245 // walk backward from the other side, adding more triangles and flipping | |
246 q := e.prev | |
247 for area(p, points[q.prev.i], points[q.i]) < 0 { | |
248 t = tri.addTriangle(q.prev.i, i, q.i, -1, q.t, q.prev.t) | |
249 tri.legalize(t + 2) | |
250 q.prev.t = t | |
251 tri.hull = q.remove() | |
252 q = q.prev | |
253 } | |
254 } | |
255 | |
256 // save the two new edges in the hash table | |
257 tri.hashEdge(e) | |
258 tri.hashEdge(e.prev) | |
259 } | |
260 | |
261 tri.triangles = tri.triangles[:tri.trianglesLen] | |
262 tri.halfedges = tri.halfedges[:tri.trianglesLen] | |
263 | |
264 return nil | |
265 } | |
266 | |
267 func (tri *triangulator) hashKey(point Vertex) int { | |
268 d := point.Sub2D(tri.center) | |
269 return int(pseudoAngle(d.X, d.Y) * float64(len(tri.hash))) | |
270 } | |
271 | |
272 func (tri *triangulator) hashEdge(e *node) { | |
273 tri.hash[tri.hashKey(tri.points[e.i])] = e | |
274 } | |
275 | |
276 func (tri *triangulator) addTriangle(i0, i1, i2, a, b, c int32) int32 { | |
277 i := int32(tri.trianglesLen) | |
278 tri.triangles[i] = i0 | |
279 tri.triangles[i+1] = i1 | |
280 tri.triangles[i+2] = i2 | |
281 tri.link(i, a) | |
282 tri.link(i+1, b) | |
283 tri.link(i+2, c) | |
284 tri.trianglesLen += 3 | |
285 return i | |
286 } | |
287 | |
288 func (tri *triangulator) link(a, b int32) { | |
289 tri.halfedges[a] = b | |
290 if b >= 0 { | |
291 tri.halfedges[b] = a | |
292 } | |
293 } | |
294 | |
295 func (tri *triangulator) legalize(a int32) int32 { | |
296 // if the pair of triangles doesn't satisfy the Delaunay condition | |
297 // (p1 is inside the circumcircle of [p0, pl, pr]), flip them, | |
298 // then do the same check/flip recursively for the new pair of triangles | |
299 // | |
300 // pl pl | |
301 // /||\ / \ | |
302 // al/ || \bl al/ \a | |
303 // / || \ / \ | |
304 // / a||b \ flip /___ar___\ | |
305 // p0\ || /p1 => p0\---bl---/p1 | |
306 // \ || / \ / | |
307 // ar\ || /br b\ /br | |
308 // \||/ \ / | |
309 // pr pr | |
310 | |
311 b := tri.halfedges[a] | |
312 | |
313 a0 := a - a%3 | |
314 b0 := b - b%3 | |
315 | |
316 al := a0 + (a+1)%3 | |
317 ar := a0 + (a+2)%3 | |
318 bl := b0 + (b+2)%3 | |
319 | |
320 if b < 0 { | |
321 return ar | |
322 } | |
323 | |
324 p0 := tri.triangles[ar] | |
325 pr := tri.triangles[a] | |
326 pl := tri.triangles[al] | |
327 p1 := tri.triangles[bl] | |
328 | |
329 illegal := inCircle(tri.points[p0], tri.points[pr], tri.points[pl], tri.points[p1]) | |
330 | |
331 if illegal { | |
332 tri.triangles[a] = p1 | |
333 tri.triangles[b] = p0 | |
334 | |
335 // edge swapped on the other side of the hull (rare) | |
336 // fix the halfedge reference | |
337 if tri.halfedges[bl] == -1 { | |
338 e := tri.hull | |
339 for { | |
340 if e.t == bl { | |
341 e.t = a | |
342 break | |
343 } | |
344 e = e.next | |
345 if e == tri.hull { | |
346 break | |
347 } | |
348 } | |
349 } | |
350 | |
351 tri.link(a, tri.halfedges[bl]) | |
352 tri.link(b, tri.halfedges[ar]) | |
353 tri.link(ar, bl) | |
354 | |
355 br := b0 + (b+1)%3 | |
356 | |
357 tri.legalize(a) | |
358 return tri.legalize(br) | |
359 } | |
360 | |
361 return ar | |
362 } | |
363 | |
364 func (tri *triangulator) convexHull() []Vertex { | |
365 var result []Vertex | |
366 e := tri.hull | |
367 for e != nil { | |
368 result = append(result, tri.points[e.i]) | |
369 e = e.prev | |
370 if e == tri.hull { | |
371 break | |
372 } | |
373 } | |
374 return result | |
375 } |