Mercurial > gemma
view pkg/mesh/triangulator.go @ 5418:e89ff1894bb4 marking-single-beam
Don't clip the points against the border polygon. The implemented way was wrong.
author | Sascha L. Teichmann <sascha.teichmann@intevation.de> |
---|---|
date | Wed, 07 Jul 2021 18:27:40 +0200 |
parents | f4abfd0ee8ad |
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 mesh 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 }