问题描述
我正在从事一个项目,该项目涉及从潜在的凹多边形生成三角形网格。我对高性能解决方案比最佳解决方案更感兴趣。我遇到的大部分内容都涉及使用 OpenGL(不是一个选项)或专注于优化或近似优化分解为凸多边形,使用 O(n3) 或 O(n4)。剪耳似乎是最直接的解决方案,应该是 O(n2):
for /* ever */ {
for each(vertex) {
candidate := Triangle{prevIoUsvertex,currentVertex,nextVertex}
if !candidate.IsInterior() || candidate.ContainsAnyOtherVertex() {
continue
}
result.Push(candidate)
vertex = vertex.Except(currentVertex)
}
if len(vertex) == 3 {
result.Push(Triangle(vertex))
break
}
}
我不明白这怎么可能是 O(n2)。在我看来,这必须涉及三个嵌套循环,每个循环都与 N 成正比:上面的两个显式循环和 candidate.ContainsAnyOtherVertex()
。我错过了什么吗?是否有某种方法可以检测候选三角形是否包含任何其他剩余顶点,而无需遍历所有剩余顶点?
或者,我应该使用其他算法吗?我对分解为凸多边形的方法没问题,但 A) 我不在乎解决方案是否最优,B) 我对阅读 Cgal 的源代码或学术文献不感兴趣。
解决方法
如果你这样做的话,剪耳是 O(n2):
- 首先,对于每个角度小于 180 度的顶点,计算其角度内的顶点数。 (O(n2))
- 您可以使用 0 计数剪切 + 删除任何此类顶点。您将执行此 O(n) 次。当您移除 vetex 时:
- 首先将它从它在 (O(n)) 中的任何三角形的计数中删除;和
- 计算您通过裁剪形成的最多 2 个新三角形中的其他顶点(也是 O(n))。
对于其他遇到此问题的人,这里是我使用的精简版本:
package mesh
import "sort"
type Vec []float32
type Face []Vec
type FaceTri struct {
face Face
index [3]int
}
func triangulate(face Face,offset int) [][3]int {
normal := face.Normal()
index := make([]int,len(face))
convex := map[int]bool{}
reflex := map[int]bool{}
ear := map[int]bool{}
// Mark convex and reflex
for i := range index {
index[i] = i
tri := face.tri(i,i+1,i+2)
// Skip colinear vertices
if tri.Area().Len() == 0 {
continue
}
if tri.Area().Dot(normal) > 0 {
convex[tri.index[1]] = true
} else {
reflex[tri.index[1]] = true
}
}
// Mark ears
for i := range convex {
if isEar(face.tri(i-1,i,i+1),reflex) {
ear[i] = true
}
}
var tris [][3]int
for len(reflex) > 0 || len(index) > 3 {
ni := len(index)
// Next ear
var j int
for j = range ear {
break
}
// Find J in the list of unclipped vertices and get 2 neighbors to each side
x := sort.SearchInts(index,j)
h,k,l := (ni+x-2)%ni,(ni+x-1)%ni,(x+1)%ni,(x+2)%ni
h,l = index[h],index[i],index[k],index[l]
// Clip (I,J,K)
index = append(index[:x],index[x+1:]...)
tris = append(tris,[3]int{offset + i,offset + j,offset + k})
delete(ear,j)
delete(convex,j)
// Update prior vertex
update(face.tri(h,k),normal,reflex,convex,ear)
// Update later vertex
if h != l {
update(face.tri(i,l),ear)
}
}
tris = append(tris,[3]int{
offset + index[0],offset + index[1],offset + index[2],})
return tris
}
func update(tri *FaceTri,faceNormal Vec,ear map[int]bool) {
idx := tri.index[1]
wasConvex,wasEar := convex[idx],ear[idx]
isConvex := wasConvex || tri.Area().Dot(faceNormal) > 0
if !wasConvex && isConvex {
convex[idx] = true
delete(reflex,idx)
}
if !wasEar && isConvex && isEar(tri,reflex) {
ear[idx] = true
}
}
func isEar(tri *FaceTri,reflex map[int]bool) bool {
// It is sufficient to only check reflex vertices - a convex vertex
// cannot be contained without a reflex vertex also being contained
for j := range reflex {
if tri.HasIndex(j) {
continue
}
if tri.ContainsPoint(tri.face[j]) {
return false
}
}
return true
}
// Euclidean length of V
func (v Vec) Len() float32
// Dot product of A and B
func (a Vec) Dot(b Vec) float32
// Calculates the normal vector to the face - for concave polygons,this
// is the summation of the normals of each vertex (normalized)
func (Face) Normal() Vec
// Tri returns a FaceTri for I,K,modulo len(f)
func (Face) tri(i,j,k int) *FaceTri
// Area returns the cross product of the vector from v1 to v0 and the vector
// from v1 to v2
func (*FaceTri) Area() Vec
// Returns true if v is in f.index
func (f *FaceTri) HasIndex(v int) bool
// Returns true if P is contained within the described triangle
func (*FaceTri) ContainsPoint(p Vec) bool
这是基于 https://www.geometrictools.com/Documentation/TriangulationByEarClipping.pdf 的描述。我认为这与马特描述的基本相同。
为了清晰和简洁起见,我省略了一些方法的内容,以及以降低可读性的方式重用结果的情况。 IMO 我遗漏的唯一有趣的部分是 ContainsPoint
。在 2D 中,确定三角形 ABC 是否包含点 P 是直截了当的。在 3D 中,不那么重要,因为 P 不一定与 ABC 共面:
V = vector from B to A
U = vector from B to C
Q = vector from B to P
M = MatrixFromColumns(V,U,V x U)
X = (inverse of M) * Q
Q === (V * X[0]) + (U * X[1]) + ((V x U) * X[2])
If X[2] is zero,then Q has a component out of the plane of the triangle
P is in ABC if and only if X[0] and X[1] are positive and X[0] + X[1] <= 1