Mercurial > gemma
changeset 3604:6248a4bc10c7
Moved triangle elimination to triangulation code.
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Wed, 05 Jun 2019 11:39:25 +0200 |
parents | 2cf7e167dad9 |
children | d02d4e31491b |
files | pkg/imports/sr.go pkg/octree/triangulation.go |
diffstat | 2 files changed, 81 insertions(+), 32 deletions(-) [+] |
line wrap: on
line diff
--- a/pkg/imports/sr.go Wed Jun 05 09:38:59 2019 +0200 +++ b/pkg/imports/sr.go Wed Jun 05 11:39:25 2019 +0200 @@ -28,13 +28,11 @@ "os" "path" "path/filepath" - "sort" "strconv" "strings" "time" shp "github.com/jonas-p/go-shp" - "gonum.org/v1/gonum/stat" "gemma.intevation.de/gemma/pkg/common" "gemma.intevation.de/gemma/pkg/misc" @@ -366,36 +364,7 @@ feedback.Info("No boundary given.") feedback.Info("Eliminate triangles with long edges.") - lens := make([]float64, len(tri.Triangles)/3) - for i := 0; i < len(tri.Triangles); i += 3 { - vs := tri.Triangles[i : i+3] - var max float64 - for j := 0; j < 3; j++ { - k := (j + 1) % 3 - if l := tri.Points[vs[j]].SquaredDistance2D(tri.Points[vs[k]]); l > max { - max = l - } - } - lens[i/3] = max - } - qs := make([]float64, len(lens)) - copy(qs, lens) - sort.Float64s(qs) - q := stat.Quantile(0.98, stat.LinInterp, qs, nil) - feedback.Info("98%% quantile of the edge lengths: %.2f\n", math.Sqrt(q)) - - ts := tri.Triangles - - for i := len(lens) - 1; i >= 0; i-- { - if lens[i] >= q { - if i != len(lens)-1 { - copy(ts[i*3:], ts[(i+1)*3:]) - } - ts = ts[:len(ts)-3] - } - } - feedback.Info("Remaining triangles: %d (%.2f%%)\n", - len(ts)/3, float64(len(ts))/float64(len(tri.Triangles))*100.0) + tri.ConcaveHull(0.98) } // TODO: Implement me!
--- a/pkg/octree/triangulation.go Wed Jun 05 09:38:59 2019 +0200 +++ b/pkg/octree/triangulation.go Wed Jun 05 11:39:25 2019 +0200 @@ -21,7 +21,11 @@ import ( "fmt" + "log" "math" + "sort" + + "gonum.org/v1/gonum/stat" ) type Triangulation struct { @@ -38,6 +42,82 @@ return &Triangulation{points, t.convexHull(), t.triangles, t.halfedges}, err } +func (t *Triangulation) ConcaveHull(quantile float64) []int32 { + + numTris := len(t.Triangles) / 3 + + // Calculate the max edge lengths per triangle. + lengths := make([]float64, numTris) + for i := range lengths { + 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 { + max = l + } + } + lengths[i] = max + } + + qs := make([]float64, len(lengths)) + copy(qs, lengths) + sort.Float64s(qs) + q := stat.Quantile(quantile, stat.LinInterp, qs, nil) + log.Printf("info: %.2f%% quantile: %.2f\n", quantile, math.Sqrt(q)) + + removed := map[int32]bool{} + + isBorder := func(n int32) bool { + n *= 3 + for i := int32(0); i < 3; i++ { + e := n + i + if t.Halfedges[e] == -1 || removed[e] { + return true + } + } + return false + } + + oldCandidates := make([]int32, 0, int(math.Ceil(float64(numTris)*(1.0-quantile)))) + log.Printf("info: cap candidates: %d\n", cap(oldCandidates)) + + for i, length := range lengths { + if length >= q { + oldCandidates = append(oldCandidates, int32(i)) + } + } + + var newCandidates []int32 + + for len(oldCandidates) > 0 { + log.Printf("candidates: %d\n", len(oldCandidates)) + + oldRemoved := len(removed) + + for _, i := range oldCandidates { + + if isBorder(i) { + removed[i] = true + } else { + newCandidates = append(newCandidates, i) + } + } + + if oldRemoved == len(removed) { + break + } + + oldCandidates = newCandidates + newCandidates = newCandidates[:0] + } + + log.Printf("info: triangles to remove: %d\n", len(removed)) + + return nil +} + func (t *Triangulation) TriangleSlices() [][]int32 { sl := make([][]int32, len(t.Triangles)/3) var j int