Mercurial > gemma
view 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 |
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 mesh import ( "fmt" "log" "math" "gonum.org/v1/gonum/stat" ) 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) EstimateTooLong() float64 { num := len(t.Triangles) / 3 lengths := make([]float64, 0, num) points := t.Points tris: for i := 0; i < num; i++ { idx := i * 3 var max float64 vs := t.Triangles[idx : idx+3] for j, vj := range vs { if t.Halfedges[idx+j] < 0 { continue tris } k := (j + 1) % 3 if l := points[vj].Distance2D(points[vs[k]]); l > max { max = l } } lengths = append(lengths, max) } std := stat.StdDev(lengths, nil) return 3.5 * std } func (t *Triangulation) ConcaveHull(tooLong float64) (LineStringZ, map[int32]struct{}) { if tooLong <= 0 { tooLong = t.EstimateTooLong() } tooLong *= tooLong var candidates []int32 points := t.Points for i, num := 0, len(t.Triangles)/3; i < num; i++ { idx := i * 3 var max float64 vs := t.Triangles[idx : idx+3] for j, vj := range vs { k := (j + 1) % 3 if l := points[vj].SquaredDistance2D(points[vs[k]]); l > max { max = l } } if max > tooLong { candidates = append(candidates, int32(i)) } } removed := map[int32]struct{}{} isBorder := func(n int32) bool { n *= 3 for i := int32(0); i < 3; i++ { e := n + i o := t.Halfedges[e] if o < 0 { return true } if _, found := removed[o/3]; found { return true } } return false } var newCandidates []int32 log.Printf("info: candidates: %d\n", len(candidates)) for len(candidates) > 0 { oldRemoved := len(removed) for _, i := range candidates { if isBorder(i) { removed[i] = struct{}{} } else { newCandidates = append(newCandidates, i) } } if oldRemoved == len(removed) { break } candidates = newCandidates newCandidates = newCandidates[:0] } log.Printf("info: candidates left: %d\n", len(candidates)) log.Printf("info: triangles: %d\n", len(t.Triangles)/3) log.Printf("info: triangles to remove: %d\n", len(removed)) type edge struct { a, b int32 prev, next *edge } isClosed := func(e *edge) bool { for curr := e.next; curr != nil; curr = curr.next { if curr == e { return true } } return false } open := map[int32]*edge{} var rings []*edge for i, num := int32(0), int32(len(t.Triangles)/3); i < num; i++ { if _, found := removed[i]; found { continue } n := i * 3 for j := int32(0); j < 3; j++ { e := n + j f := t.Halfedges[e] if f >= 0 { if _, found := removed[f/3]; !found { continue } } a := t.Triangles[e] b := t.Triangles[n+(j+1)%3] curr := &edge{a: a, b: b} if old := open[a]; old != nil { delete(open, a) if old.a == a { old.prev = curr curr.next = old } else { old.next = curr curr.prev = old } if isClosed(old) { rings = append(rings, old) } } else { open[a] = curr } if old := open[b]; old != nil { delete(open, b) if old.b == b { old.next = curr curr.prev = old } else { old.prev = curr curr.next = old } if isClosed(old) { rings = append(rings, old) } } else { open[b] = curr } } } if len(open) > 0 { log.Printf("warn: open vertices left: %d\n", len(open)) } if len(rings) == 0 { log.Println("warn: no ring found") return nil, removed } curr := rings[0] polygon := LineStringZ{ points[curr.a], points[curr.b], } for curr = curr.next; curr != rings[0]; curr = curr.next { polygon = append(polygon, points[curr.b]) } polygon = append(polygon, t.Points[rings[0].a]) log.Printf("length of boundary: %d\n", len(polygon)) return polygon, removed } 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 }