Mercurial > gemma
view pkg/octree/triangulator.go @ 4488:bff6c5c1db4f
client: pdf-gen: improve adding bottleneck info to pdf
* Check if the bottleneck is in the current view to add its info to the exported pdf and the pdf filename, this avoid wrong filename and wrong info in pdf in case view has been changed to another location.
* Set the bottleneck to print after moving to it in map.
author | Fadi Abbud <fadi.abbud@intevation.de> |
---|---|
date | Fri, 27 Sep 2019 11:15:02 +0200 |
parents | d12c2f4d3483 |
children |
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 octree import ( "errors" "math" "sort" ) type triangulator struct { points []Vertex squaredDistances []float64 ids []int32 center Vertex triangles []int32 halfedges []int32 trianglesLen int hull *node hash []*node } func newTriangulator(points []Vertex) *triangulator { return &triangulator{points: points} } // sorting a triangulator sorts the `ids` such that the referenced points // are in order by their distance to `center` func (tri *triangulator) Len() int { return len(tri.points) } func (tri *triangulator) Swap(i, j int) { tri.ids[i], tri.ids[j] = tri.ids[j], tri.ids[i] } func (tri *triangulator) Less(i, j int) bool { d1 := tri.squaredDistances[tri.ids[i]] d2 := tri.squaredDistances[tri.ids[j]] if d1 != d2 { return d1 < d2 } p1 := tri.points[tri.ids[i]] p2 := tri.points[tri.ids[j]] if p1.X != p2.X { return p1.X < p2.X } return p1.Y < p2.Y } func (tri *triangulator) triangulate() error { points := tri.points n := len(points) if n == 0 { return nil } tri.ids = make([]int32, n) // compute bounds x0 := points[0].X y0 := points[0].Y x1 := points[0].X y1 := points[0].Y for i, p := range points { if p.X < x0 { x0 = p.X } if p.X > x1 { x1 = p.X } if p.Y < y0 { y0 = p.Y } if p.Y > y1 { y1 = p.Y } tri.ids[i] = int32(i) } var i0, i1, i2 int32 // pick a seed point close to midpoint m := Vertex{X: (x0 + x1) / 2, Y: (y0 + y1) / 2} minDist := infinity for i, p := range points { d := p.SquaredDistance2D(m) if d < minDist { i0 = int32(i) minDist = d } } // find point closest to seed point minDist = infinity for i, p := range points { if int32(i) == i0 { continue } d := p.SquaredDistance2D(points[i0]) if d > 0 && d < minDist { i1 = int32(i) minDist = d } } // find the third point which forms the smallest circumcircle minRadius := infinity for i, p := range points { if int32(i) == i0 || int32(i) == i1 { continue } r := circumradius(points[i0], points[i1], p) if r < minRadius { i2 = int32(i) minRadius = r } } if minRadius == infinity { return errors.New("no delaunay triangulation exists for this input") } // swap the order of the seed points for counter-clockwise orientation if area(points[i0], points[i1], points[i2]) < 0 { i1, i2 = i2, i1 } tri.center = circumcenter(points[i0], points[i1], points[i2]) // sort the points by distance from the seed triangle circumcenter tri.squaredDistances = make([]float64, n) for i, p := range points { tri.squaredDistances[i] = p.SquaredDistance2D(tri.center) } sort.Sort(tri) // initialize a hash table for storing edges of the advancing convex hull hashSize := int(math.Ceil(math.Sqrt(float64(n)))) tri.hash = make([]*node, hashSize) // initialize a circular doubly-linked list that will hold an advancing convex hull nodes := make([]node, n) e := newNode(nodes, i0, nil) e.t = 0 tri.hashEdge(e) e = newNode(nodes, i1, e) e.t = 1 tri.hashEdge(e) e = newNode(nodes, i2, e) e.t = 2 tri.hashEdge(e) tri.hull = e maxTriangles := 2*n - 5 tri.triangles = make([]int32, maxTriangles*3) tri.halfedges = make([]int32, maxTriangles*3) tri.addTriangle(i0, i1, i2, -1, -1, -1) pp := Vertex{X: infinity, Y: infinity} for k := 0; k < n; k++ { i := tri.ids[k] p := points[i] // skip nearly-duplicate points if p.SquaredDistance2D(pp) < eps { continue } pp = p // skip seed triangle points if i == int32(i0) || i == int32(i1) || i == int32(i2) { continue } // find a visible edge on the convex hull using edge hash var start *node key := tri.hashKey(p) for j := 0; j < len(tri.hash); j++ { start = tri.hash[key] if start != nil && start.i >= 0 { break } key++ if key >= len(tri.hash) { key = 0 } } start = start.prev e := start for area(p, points[e.i], points[e.next.i]) >= 0 { e = e.next if e == start { e = nil break } } if e == nil { // likely a near-duplicate point; skip it continue } walkBack := e == start // add the first triangle from the point t := tri.addTriangle(e.i, i, e.next.i, -1, -1, e.t) e.t = t // keep track of boundary triangles on the hull e = newNode(nodes, i, e) // recursively flip triangles from the point until they satisfy the Delaunay condition e.t = tri.legalize(t + 2) // walk forward through the hull, adding more triangles and flipping recursively q := e.next for area(p, points[q.i], points[q.next.i]) < 0 { t = tri.addTriangle(q.i, i, q.next.i, q.prev.t, -1, q.t) q.prev.t = tri.legalize(t + 2) tri.hull = q.remove() q = q.next } if walkBack { // walk backward from the other side, adding more triangles and flipping q := e.prev for area(p, points[q.prev.i], points[q.i]) < 0 { t = tri.addTriangle(q.prev.i, i, q.i, -1, q.t, q.prev.t) tri.legalize(t + 2) q.prev.t = t tri.hull = q.remove() q = q.prev } } // save the two new edges in the hash table tri.hashEdge(e) tri.hashEdge(e.prev) } tri.triangles = tri.triangles[:tri.trianglesLen] tri.halfedges = tri.halfedges[:tri.trianglesLen] return nil } func (tri *triangulator) hashKey(point Vertex) int { d := point.Sub2D(tri.center) return int(pseudoAngle(d.X, d.Y) * float64(len(tri.hash))) } func (tri *triangulator) hashEdge(e *node) { tri.hash[tri.hashKey(tri.points[e.i])] = e } func (tri *triangulator) addTriangle(i0, i1, i2, a, b, c int32) int32 { i := int32(tri.trianglesLen) tri.triangles[i] = i0 tri.triangles[i+1] = i1 tri.triangles[i+2] = i2 tri.link(i, a) tri.link(i+1, b) tri.link(i+2, c) tri.trianglesLen += 3 return i } func (tri *triangulator) link(a, b int32) { tri.halfedges[a] = b if b >= 0 { tri.halfedges[b] = a } } func (tri *triangulator) legalize(a int32) int32 { // if the pair of triangles doesn't satisfy the Delaunay condition // (p1 is inside the circumcircle of [p0, pl, pr]), flip them, // then do the same check/flip recursively for the new pair of triangles // // pl pl // /||\ / \ // al/ || \bl al/ \a // / || \ / \ // / a||b \ flip /___ar___\ // p0\ || /p1 => p0\---bl---/p1 // \ || / \ / // ar\ || /br b\ /br // \||/ \ / // pr pr b := tri.halfedges[a] a0 := a - a%3 b0 := b - b%3 al := a0 + (a+1)%3 ar := a0 + (a+2)%3 bl := b0 + (b+2)%3 if b < 0 { return ar } p0 := tri.triangles[ar] pr := tri.triangles[a] pl := tri.triangles[al] p1 := tri.triangles[bl] illegal := inCircle(tri.points[p0], tri.points[pr], tri.points[pl], tri.points[p1]) if illegal { tri.triangles[a] = p1 tri.triangles[b] = p0 // edge swapped on the other side of the hull (rare) // fix the halfedge reference if tri.halfedges[bl] == -1 { e := tri.hull for { if e.t == bl { e.t = a break } e = e.next if e == tri.hull { break } } } tri.link(a, tri.halfedges[bl]) tri.link(b, tri.halfedges[ar]) tri.link(ar, bl) br := b0 + (b+1)%3 tri.legalize(a) return tri.legalize(br) } return ar } func (tri *triangulator) convexHull() []Vertex { var result []Vertex e := tri.hull for e != nil { result = append(result, tri.points[e.i]) e = e.prev if e == tri.hull { break } } return result }