Mercurial > gemma

comparison pkg/octree/triangulation.go @ 2483:620038ade708 octree-diff

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Incorporated fogleman's fast Delaunay triangulation adjuted to our vertex model. License: MIT Home: https://github.com/fogleman/delaunay
author Sascha L. Teichmann <sascha.teichmann@intevation.de>
date Fri, 01 Mar 2019 15:33:27 +0100
parents
children 4fa92d468164
comparison
equal deleted inserted replaced
2481:3cf5d27a6c8b 2483:620038ade708
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 octree
21
22 import (
23 "fmt"
24 "math"
25 )
26
27 type Triangulation struct {
28 Points []Vertex
29 ConvexHull []Vertex
30 Triangles []int
31 Halfedges []int
32 }
33
34 // Triangulate returns a Delaunay triangulation of the provided points.
35 func Triangulate(points []Vertex) (*Triangulation, error) {
36 t := newTriangulator(points)
37 err := t.triangulate()
38 return &Triangulation{points, t.convexHull(), t.triangles, t.halfedges}, err
39 }
40
41 func (t *Triangulation) area() float64 {
42 var result float64
43 points := t.Points
44 ts := t.Triangles
45 for i := 0; i < len(ts); i += 3 {
46 p0 := points[ts[i+0]]
47 p1 := points[ts[i+1]]
48 p2 := points[ts[i+2]]
49 result += area(p0, p1, p2)
50 }
51 return result / 2
52 }
53
54 // Validate performs several sanity checks on the Triangulation to check for
55 // potential errors. Returns nil if no issues were found. You normally
56 // shouldn't need to call this function but it can be useful for debugging.
57 func (t *Triangulation) Validate() error {
58 // verify halfedges
59 for i1, i2 := range t.Halfedges {
60 if i1 != -1 && t.Halfedges[i1] != i2 {
61 return fmt.Errorf("invalid halfedge connection")
62 }
63 if i2 != -1 && t.Halfedges[i2] != i1 {
64 return fmt.Errorf("invalid halfedge connection")
65 }
66 }
67
68 // verify convex hull area vs sum of triangle areas
69 hull1 := t.ConvexHull
70 hull2 := ConvexHull(t.Points)
71 area1 := polygonArea(hull1)
72 area2 := polygonArea(hull2)
73 area3 := t.area()
74 if math.Abs(area1-area2) > 1e-9 || math.Abs(area1-area3) > 1e-9 {
75 return fmt.Errorf("hull areas disagree: %f, %f, %f", area1, area2, area3)
76 }
77
78 // verify convex hull perimeter
79 perimeter1 := polygonPerimeter(hull1)
80 perimeter2 := polygonPerimeter(hull2)
81 if math.Abs(perimeter1-perimeter2) > 1e-9 {
82 return fmt.Errorf("hull perimeters disagree: %f, %f", perimeter1, perimeter2)
83 }
84
85 return nil
86 }