Mercurial > gemma
view pkg/octree/simplify.go @ 4316:3d6a2c6b436c
shape upload stretch import: Store features in database. Still broken [WIP].
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Tue, 03 Sep 2019 19:01:51 +0200 |
parents | 33fa76994b8a |
children |
line wrap: on
line source
// This is Free Software under GNU Affero General Public License v >= 3.0 // without warranty, see README.md and license for details. // // SPDX-License-Identifier: AGPL-3.0-or-later // License-Filename: LICENSES/AGPL-3.0.txt // // Copyright (C) 2019 by via donau // – Österreichische Wasserstraßen-Gesellschaft mbH // Software engineering by Intevation GmbH // // Author(s): // * Sascha L. Teichmann <sascha.teichmann@intevation.de> package octree import ( "math" ) func handleTriangle( t *Triangle, maxDist, tolerance float64, maxIdx int, points MultiPointZ, result *MultiPointZ, ) bool { if maxDist <= tolerance { return false } if len(points) == 1 { *result = append(*result, points[0]) return true } var ( tris [3]Triangle planes [3]Plane3D maxDists [3]float64 maxIdxs [3]int parts [3]MultiPointZ ) top := points[maxIdx] for i := 0; i < 3; i++ { tris[i] = Triangle{t[i], t[(i+1)%3], top} planes[i] = tris[i].Plane3D() } nextPoint: for i, v := range points { if i == maxIdx { continue } for j := range tris { if tris[j].Contains(v.X, v.Y) { if dist := math.Abs(planes[j].Eval(v)); dist > maxDists[j] { maxDists[j] = dist maxIdxs[j] = len(parts[j]) } parts[j] = append(parts[j], v) continue nextPoint } } } var found bool for i, part := range parts { if len(part) > 0 && handleTriangle( &tris[i], maxDists[i], tolerance, maxIdxs[i], part, result, ) { found = true } } if found { *result = append(*result, top) } return found } func (points MultiPointZ) Simplify(tolerance float64) MultiPointZ { if len(points) < 2 { return points } if tolerance < 0 { tolerance = -tolerance } min := Vertex{X: math.MaxFloat64, Y: math.MaxFloat64, Z: math.MaxFloat64} max := Vertex{X: -math.MaxFloat64, Y: -math.MaxFloat64, Z: -math.MaxFloat64} var maxIdx int for i, v := range points { min.Minimize(v) if v.X < min.X { min.X = v.X } if v.X > max.X { max.X = v.X } if v.Y < min.Y { min.Y = v.Y } if v.Y > max.Y { max.Y = v.Y } if v.Z < min.Z { min.Z = v.Z } if v.Z > max.Z { max.Z = v.Z maxIdx = i } max.Maximize(v) } /* log.Printf("(%.5f, %.5f, %.5f) - (%.5f, %.5f, %.5f)\n", min.X, min.Y, min.Z, max.X, max.Y, max.Z) */ below := min.Z - 3*tolerance xMin := min.X - tolerance xMax := max.X + tolerance yMin := min.Y - tolerance yMax := max.Y + tolerance corners := []Vertex{ {xMin, yMin, below}, {xMax, yMin, below}, {xMax, yMax, below}, {xMin, yMax, below}, } top := points[maxIdx] tris := make([]Triangle, len(corners)) planes := make([]Plane3D, len(corners)) for i, v1 := range corners { v2 := corners[(i+1)%len(corners)] tris[i] = Triangle{v1, v2, top} planes[i] = tris[i].Plane3D() } parts := make([][]Vertex, len(tris)) maxDists := make([]float64, len(planes)) maxIdxs := make([]int, len(planes)) nextPoint: for i, v := range points { if i == maxIdx { continue } for j := range tris { if tris[j].Contains(v.X, v.Y) { if dist := math.Abs(planes[j].Eval(v)); dist > maxDists[j] { maxDists[j] = dist maxIdxs[j] = len(parts[j]) } parts[j] = append(parts[j], v) continue nextPoint } } } result := make(MultiPointZ, 0, len(points)) var found bool for i, part := range parts { if len(part) > 0 && handleTriangle( &tris[i], maxDists[i], tolerance, maxIdxs[i], part, &result, ) { found = true } } if found { result = append(result, top) } return result }