Mercurial > gemma

comparison pkg/mesh/triangulation.go @ 4828:39ee00d06309

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merged remove-octree-debris branch back into default.
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Tue, 05 Nov 2019 14:31:22 +0100
parents f4abfd0ee8ad
children 3f0382e9f302
comparison
equal deleted inserted replaced
4825:5eb05714353a 4828:39ee00d06309
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 }