Mercurial > gemma
view pkg/octree/triangulation.go @ 3534:034657d6604f waterlevel-in-crossprofile
client: fairway profiles: improved legend
author | Markus Kottlaender <markus@intevation.de> |
---|---|
date | Wed, 29 May 2019 18:20:53 +0200 |
parents | 4fa92d468164 |
children | 6248a4bc10c7 |
line wrap: on
line source
// Copyright (C) 2018 Michael Fogleman // // Permission is hereby granted, free of charge, to any person obtaining // a copy of this software and associated documentation files (the "Software"), // to deal in the Software without restriction, including without limitation // the rights to use, copy, modify, merge, publish, distribute, sublicense, // and/or sell copies of the Software, and to permit persons to whom the // Software is furnished to do so, subject to the following conditions: // // The above copyright notice and this permission notice shall be included // in all copies or substantial portions of the Software. // // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, // OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. package octree import ( "fmt" "math" ) type Triangulation struct { Points []Vertex ConvexHull []Vertex Triangles []int32 Halfedges []int32 } // Triangulate returns a Delaunay triangulation of the provided points. func Triangulate(points []Vertex) (*Triangulation, error) { t := newTriangulator(points) err := t.triangulate() return &Triangulation{points, t.convexHull(), t.triangles, t.halfedges}, err } func (t *Triangulation) TriangleSlices() [][]int32 { sl := make([][]int32, len(t.Triangles)/3) var j int for i := range sl { sl[i] = t.Triangles[j : j+3] j += 3 } return sl } func (t *Triangulation) Tin() *Tin { min := Vertex{math.MaxFloat64, math.MaxFloat64, math.MaxFloat64} max := Vertex{-math.MaxFloat64, -math.MaxFloat64, -math.MaxFloat64} vertices := t.Points for _, v := range vertices { min.Minimize(v) max.Maximize(v) } return &Tin{ Vertices: vertices, Triangles: t.TriangleSlices(), Min: min, Max: max, } } func (t *Triangulation) area() float64 { var result float64 points := t.Points ts := t.Triangles for i := 0; i < len(ts); i += 3 { p0 := points[ts[i+0]] p1 := points[ts[i+1]] p2 := points[ts[i+2]] result += area(p0, p1, p2) } return result / 2 } // Validate performs several sanity checks on the Triangulation to check for // potential errors. Returns nil if no issues were found. You normally // shouldn't need to call this function but it can be useful for debugging. func (t *Triangulation) Validate() error { // verify halfedges for i1, i2 := range t.Halfedges { if i1 != -1 && t.Halfedges[i1] != i2 { return fmt.Errorf("invalid halfedge connection") } if i2 != -1 && t.Halfedges[i2] != int32(i1) { return fmt.Errorf("invalid halfedge connection") } } // verify convex hull area vs sum of triangle areas hull1 := t.ConvexHull hull2 := ConvexHull(t.Points) area1 := polygonArea(hull1) area2 := polygonArea(hull2) area3 := t.area() if math.Abs(area1-area2) > 1e-9 || math.Abs(area1-area3) > 1e-9 { return fmt.Errorf("hull areas disagree: %f, %f, %f", area1, area2, area3) } // verify convex hull perimeter perimeter1 := polygonPerimeter(hull1) perimeter2 := polygonPerimeter(hull2) if math.Abs(perimeter1-perimeter2) > 1e-9 { return fmt.Errorf("hull perimeters disagree: %f, %f", perimeter1, perimeter2) } return nil }