从凹多边形生成三角形网格的高效算法

问题描述

我正在从事一个项目,该项目涉及从潜在的凹多边形生成三角形网格。我对高性能解决方案比最佳解决方案更感兴趣。我遇到的大部分内容都涉及使用 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