Mercurial > gemma
changeset 3733:ec86a7155377 concave-hull
Estimated too large triangles as triangles which have an edge which is at least 3.5 times as long as the standard dev of the longest egde per inner triangle.
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Mon, 24 Jun 2019 11:39:09 +0200 |
parents | 9fa108794e21 |
children | 879c297c47e9 |
files | pkg/imports/sr.go pkg/octree/triangulation.go |
diffstat | 2 files changed, 42 insertions(+), 5 deletions(-) [+] |
line wrap: on
line diff
--- a/pkg/imports/sr.go Thu Jun 20 19:43:25 2019 +0200 +++ b/pkg/imports/sr.go Mon Jun 24 11:39:09 2019 +0200 @@ -68,7 +68,6 @@ ) const ( - tooLongEdge = 50.0 pointsPerSquareMeter = 2 ) @@ -398,6 +397,7 @@ if boundary == nil { feedback.Info("No boundary given. Calulate from XYZ data.") + tooLongEdge := tri.EstimateTooLong() feedback.Info("Eliminate triangles with edges longer than %.2f meters.", tooLongEdge) var polygon octree.LineStringZ
--- a/pkg/octree/triangulation.go Thu Jun 20 19:43:25 2019 +0200 +++ b/pkg/octree/triangulation.go Mon Jun 24 11:39:09 2019 +0200 @@ -23,6 +23,8 @@ "fmt" "log" "math" + + "gonum.org/v1/gonum/stat" ) type Triangulation struct { @@ -39,19 +41,54 @@ 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 := t.Points[vj].SquaredDistance2D(t.Points[vs[k]]); l > max { + if l := points[vj].SquaredDistance2D(points[vs[k]]); l > max { max = l } } @@ -188,12 +225,12 @@ curr := rings[0] polygon := LineStringZ{ - t.Points[curr.a], - t.Points[curr.b], + points[curr.a], + points[curr.b], } for curr = curr.next; curr != rings[0]; curr = curr.next { - polygon = append(polygon, t.Points[curr.b]) + polygon = append(polygon, points[curr.b]) } polygon = append(polygon, t.Points[rings[0].a])