Mercurial > gemma
diff 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 diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pkg/mesh/triangulation.go Tue Nov 05 14:30:22 2019 +0100 @@ -0,0 +1,318 @@ +// 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 +}