Mercurial > gemma
comparison pkg/mesh/triangulation.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/triangulation.go@ec86a7155377 |
children | 3f0382e9f302 |
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 "fmt" | |
24 "log" | |
25 "math" | |
26 | |
27 "gonum.org/v1/gonum/stat" | |
28 ) | |
29 | |
30 type Triangulation struct { | |
31 Points []Vertex | |
32 ConvexHull []Vertex | |
33 Triangles []int32 | |
34 Halfedges []int32 | |
35 } | |
36 | |
37 // Triangulate returns a Delaunay triangulation of the provided points. | |
38 func Triangulate(points []Vertex) (*Triangulation, error) { | |
39 t := newTriangulator(points) | |
40 err := t.triangulate() | |
41 return &Triangulation{points, t.convexHull(), t.triangles, t.halfedges}, err | |
42 } | |
43 | |
44 func (t *Triangulation) EstimateTooLong() float64 { | |
45 | |
46 num := len(t.Triangles) / 3 | |
47 | |
48 lengths := make([]float64, 0, num) | |
49 | |
50 points := t.Points | |
51 | |
52 tris: | |
53 for i := 0; i < num; i++ { | |
54 idx := i * 3 | |
55 var max float64 | |
56 vs := t.Triangles[idx : idx+3] | |
57 for j, vj := range vs { | |
58 if t.Halfedges[idx+j] < 0 { | |
59 continue tris | |
60 } | |
61 k := (j + 1) % 3 | |
62 if l := points[vj].Distance2D(points[vs[k]]); l > max { | |
63 max = l | |
64 } | |
65 } | |
66 lengths = append(lengths, max) | |
67 } | |
68 | |
69 std := stat.StdDev(lengths, nil) | |
70 return 3.5 * std | |
71 } | |
72 | |
73 func (t *Triangulation) ConcaveHull(tooLong float64) (LineStringZ, map[int32]struct{}) { | |
74 | |
75 if tooLong <= 0 { | |
76 tooLong = t.EstimateTooLong() | |
77 } | |
78 | |
79 tooLong *= tooLong | |
80 | |
81 var candidates []int32 | |
82 | |
83 points := t.Points | |
84 | |
85 for i, num := 0, len(t.Triangles)/3; i < num; i++ { | |
86 idx := i * 3 | |
87 var max float64 | |
88 vs := t.Triangles[idx : idx+3] | |
89 for j, vj := range vs { | |
90 k := (j + 1) % 3 | |
91 if l := points[vj].SquaredDistance2D(points[vs[k]]); l > max { | |
92 max = l | |
93 } | |
94 } | |
95 if max > tooLong { | |
96 candidates = append(candidates, int32(i)) | |
97 } | |
98 } | |
99 | |
100 removed := map[int32]struct{}{} | |
101 | |
102 isBorder := func(n int32) bool { | |
103 n *= 3 | |
104 for i := int32(0); i < 3; i++ { | |
105 e := n + i | |
106 o := t.Halfedges[e] | |
107 if o < 0 { | |
108 return true | |
109 } | |
110 if _, found := removed[o/3]; found { | |
111 return true | |
112 } | |
113 } | |
114 return false | |
115 } | |
116 | |
117 var newCandidates []int32 | |
118 | |
119 log.Printf("info: candidates: %d\n", len(candidates)) | |
120 for len(candidates) > 0 { | |
121 | |
122 oldRemoved := len(removed) | |
123 | |
124 for _, i := range candidates { | |
125 | |
126 if isBorder(i) { | |
127 removed[i] = struct{}{} | |
128 } else { | |
129 newCandidates = append(newCandidates, i) | |
130 } | |
131 } | |
132 | |
133 if oldRemoved == len(removed) { | |
134 break | |
135 } | |
136 | |
137 candidates = newCandidates | |
138 newCandidates = newCandidates[:0] | |
139 } | |
140 | |
141 log.Printf("info: candidates left: %d\n", len(candidates)) | |
142 log.Printf("info: triangles: %d\n", len(t.Triangles)/3) | |
143 log.Printf("info: triangles to remove: %d\n", len(removed)) | |
144 | |
145 type edge struct { | |
146 a, b int32 | |
147 prev, next *edge | |
148 } | |
149 | |
150 isClosed := func(e *edge) bool { | |
151 for curr := e.next; curr != nil; curr = curr.next { | |
152 if curr == e { | |
153 return true | |
154 } | |
155 } | |
156 return false | |
157 } | |
158 | |
159 open := map[int32]*edge{} | |
160 var rings []*edge | |
161 | |
162 for i, num := int32(0), int32(len(t.Triangles)/3); i < num; i++ { | |
163 if _, found := removed[i]; found { | |
164 continue | |
165 } | |
166 n := i * 3 | |
167 for j := int32(0); j < 3; j++ { | |
168 e := n + j | |
169 f := t.Halfedges[e] | |
170 if f >= 0 { | |
171 if _, found := removed[f/3]; !found { | |
172 continue | |
173 } | |
174 } | |
175 a := t.Triangles[e] | |
176 b := t.Triangles[n+(j+1)%3] | |
177 | |
178 curr := &edge{a: a, b: b} | |
179 | |
180 if old := open[a]; old != nil { | |
181 delete(open, a) | |
182 if old.a == a { | |
183 old.prev = curr | |
184 curr.next = old | |
185 } else { | |
186 old.next = curr | |
187 curr.prev = old | |
188 } | |
189 | |
190 if isClosed(old) { | |
191 rings = append(rings, old) | |
192 } | |
193 } else { | |
194 open[a] = curr | |
195 } | |
196 | |
197 if old := open[b]; old != nil { | |
198 delete(open, b) | |
199 if old.b == b { | |
200 old.next = curr | |
201 curr.prev = old | |
202 } else { | |
203 old.prev = curr | |
204 curr.next = old | |
205 } | |
206 | |
207 if isClosed(old) { | |
208 rings = append(rings, old) | |
209 } | |
210 } else { | |
211 open[b] = curr | |
212 } | |
213 } | |
214 } | |
215 | |
216 if len(open) > 0 { | |
217 log.Printf("warn: open vertices left: %d\n", len(open)) | |
218 } | |
219 | |
220 if len(rings) == 0 { | |
221 log.Println("warn: no ring found") | |
222 return nil, removed | |
223 } | |
224 | |
225 curr := rings[0] | |
226 | |
227 polygon := LineStringZ{ | |
228 points[curr.a], | |
229 points[curr.b], | |
230 } | |
231 | |
232 for curr = curr.next; curr != rings[0]; curr = curr.next { | |
233 polygon = append(polygon, points[curr.b]) | |
234 } | |
235 | |
236 polygon = append(polygon, t.Points[rings[0].a]) | |
237 | |
238 log.Printf("length of boundary: %d\n", len(polygon)) | |
239 | |
240 return polygon, removed | |
241 } | |
242 | |
243 func (t *Triangulation) TriangleSlices() [][]int32 { | |
244 sl := make([][]int32, len(t.Triangles)/3) | |
245 var j int | |
246 for i := range sl { | |
247 sl[i] = t.Triangles[j : j+3] | |
248 j += 3 | |
249 } | |
250 return sl | |
251 } | |
252 | |
253 func (t *Triangulation) Tin() *Tin { | |
254 | |
255 min := Vertex{math.MaxFloat64, math.MaxFloat64, math.MaxFloat64} | |
256 max := Vertex{-math.MaxFloat64, -math.MaxFloat64, -math.MaxFloat64} | |
257 | |
258 vertices := t.Points | |
259 | |
260 for _, v := range vertices { | |
261 min.Minimize(v) | |
262 max.Maximize(v) | |
263 } | |
264 | |
265 return &Tin{ | |
266 Vertices: vertices, | |
267 Triangles: t.TriangleSlices(), | |
268 Min: min, | |
269 Max: max, | |
270 } | |
271 } | |
272 | |
273 func (t *Triangulation) area() float64 { | |
274 var result float64 | |
275 points := t.Points | |
276 ts := t.Triangles | |
277 for i := 0; i < len(ts); i += 3 { | |
278 p0 := points[ts[i+0]] | |
279 p1 := points[ts[i+1]] | |
280 p2 := points[ts[i+2]] | |
281 result += area(p0, p1, p2) | |
282 } | |
283 return result / 2 | |
284 } | |
285 | |
286 // Validate performs several sanity checks on the Triangulation to check for | |
287 // potential errors. Returns nil if no issues were found. You normally | |
288 // shouldn't need to call this function but it can be useful for debugging. | |
289 func (t *Triangulation) Validate() error { | |
290 // verify halfedges | |
291 for i1, i2 := range t.Halfedges { | |
292 if i1 != -1 && t.Halfedges[i1] != i2 { | |
293 return fmt.Errorf("invalid halfedge connection") | |
294 } | |
295 if i2 != -1 && t.Halfedges[i2] != int32(i1) { | |
296 return fmt.Errorf("invalid halfedge connection") | |
297 } | |
298 } | |
299 | |
300 // verify convex hull area vs sum of triangle areas | |
301 hull1 := t.ConvexHull | |
302 hull2 := ConvexHull(t.Points) | |
303 area1 := polygonArea(hull1) | |
304 area2 := polygonArea(hull2) | |
305 area3 := t.area() | |
306 if math.Abs(area1-area2) > 1e-9 || math.Abs(area1-area3) > 1e-9 { | |
307 return fmt.Errorf("hull areas disagree: %f, %f, %f", area1, area2, area3) | |
308 } | |
309 | |
310 // verify convex hull perimeter | |
311 perimeter1 := polygonPerimeter(hull1) | |
312 perimeter2 := polygonPerimeter(hull2) | |
313 if math.Abs(perimeter1-perimeter2) > 1e-9 { | |
314 return fmt.Errorf("hull perimeters disagree: %f, %f", perimeter1, perimeter2) | |
315 } | |
316 | |
317 return nil | |
318 } |