用凹域对一组点进行三角剖分 [英] Triangulate a set of points with a concave domain

问题描述

设置

给出凸包内的一组节点，假设该域包含一个或多个凹区：

其中蓝色点是点，黑线表示域。假设这些点被保存为2D数组点，长度为 n ，其中 n 是点对的数量。

让我们使用scipy.spatial的Delaunay方法对点进行三角剖分：

如您所见，可能会遇到穿过三角形的创建

问题

有什么好的算法来删除跨域的三角形？理想情况下但并非一定如此，在此情况下，单纯形边仍保留域形状（即，去除三角形的位置没有大的间隙）。

---

由于我的问题似乎仍在继续进行，因此我想跟进我当前正在使用的应用程序。

假设您定义边界后，您可以使用

在上图所示的输入中（取自我之前的 answer ） c $ c> is_inside（1.5，0.0，points，edge）将返回 True ，而 is_inside（1.5，3.0 ，点，边）将返回 False 。

请注意，<$ c上面的$ c> is_inside 函数无法处理退化的情况。我添加了两个断言来检测这种情况（您可以定义适合您的应用程序的任何epsilon值）。在许多应用程序中，这已经足够了，但是如果没有，并且遇到这些最终情况，则需要将它们分开处理。
例如，请参见此处有关实施时的鲁棒性和精度问题几何算法。

Setup

Given some set of nodes within a convex hull, assume the domain contains one or more concave areas:

where blue dots are points, and the black line illustrates the domain. Assume the points are held as a 2D array points of length n, where n is the number of point-pairs.

Let us then triangulate the points, using something like the Delaunay method from scipy.spatial:

As you can see, one may experience the creation of triangles crossing through the domain.

Question

What is a good algorithmic approach to removing any triangles that span outside of the domain? Ideally but not necessarily, where the simplex edges still preserve the domain shape (i.e., no major gaps where triangles are removed).

---

Since my question is seeming to continue to get a decent amount of activity, I wanted to follow up with the application that I'm currently using.

Assuming that you have your boundary defined, you can use a ray casting algorithm to determine whether or not the polygon is inside the domain.

To do this:

1. Take the centroid of each polygon as C_i = (x_i,y_i).
2. Then, imagine a line L = [C_i,(+inf,y_i)]: that is, a line that spans east past the end of your domain.
3. For each boundary segment s_i in boundary S, check for intersections with L. If yes, add +1 to an internal counter intersection_count; else, add nothing.

4. After the count of all intersections between L and s_i for i=1..N are calculated:

   if intersection_count % 2 == 0:
       return True # triangle outside convex hull
   else:
       return False # triangle inside convex hull
   

If your boundary is not explicitly defined, I find it helpful to 'map' the shape onto an boolean array and use a neighbor tracing algorithm to define it. Note that this approach assumes a solid domain and you will need to use a more complex algorithm for domains with 'holes' in them.

解决方案

Here is some Python code that does what you want.

First, building the alpha shape (see my previous answer):

def alpha_shape(points, alpha, only_outer=True):
    """
    Compute the alpha shape (concave hull) of a set of points.
    :param points: np.array of shape (n,2) points.
    :param alpha: alpha value.
    :param only_outer: boolean value to specify if we keep only the outer border or also inner edges.
    :return: set of (i,j) pairs representing edges of the alpha-shape. (i,j) are the indices in the points array.
    """
    assert points.shape[0] > 3, "Need at least four points"

    def add_edge(edges, i, j):
        """
        Add a line between the i-th and j-th points,
        if not in the list already
        """
        if (i, j) in edges or (j, i) in edges:
            # already added
            assert (j, i) in edges, "Can't go twice over same directed edge right?"
            if only_outer:
                # if both neighboring triangles are in shape, it's not a boundary edge
                edges.remove((j, i))
            return
        edges.add((i, j))

    tri = Delaunay(points)
    edges = set()
    # Loop over triangles:
    # ia, ib, ic = indices of corner points of the triangle
    for ia, ib, ic in tri.vertices:
        pa = points[ia]
        pb = points[ib]
        pc = points[ic]
        # Computing radius of triangle circumcircle
        # www.mathalino.com/reviewer/derivation-of-formulas/derivation-of-formula-for-radius-of-circumcircle
        a = np.sqrt((pa[0] - pb[0]) ** 2 + (pa[1] - pb[1]) ** 2)
        b = np.sqrt((pb[0] - pc[0]) ** 2 + (pb[1] - pc[1]) ** 2)
        c = np.sqrt((pc[0] - pa[0]) ** 2 + (pc[1] - pa[1]) ** 2)
        s = (a + b + c) / 2.0
        area = np.sqrt(s * (s - a) * (s - b) * (s - c))
        circum_r = a * b * c / (4.0 * area)
        if circum_r < alpha:
            add_edge(edges, ia, ib)
            add_edge(edges, ib, ic)
            add_edge(edges, ic, ia)
    return edges

To compute the edges of the outer boundary of the alpha shape use the following example call:

edges = alpha_shape(points, alpha=alpha_value, only_outer=True)

Now, after the edges of the outer boundary of the alpha-shape of points have been computed, the following function will determine whether a point (x,y) is inside the outer boundary.

def is_inside(x, y, points, edges, eps=1.0e-10):
    intersection_counter = 0
    for i, j in edges:
        assert abs((points[i,1]-y)*(points[j,1]-y)) > eps, 'Need to handle these end cases separately'
        y_in_edge_domain = ((points[i,1]-y)*(points[j,1]-y) < 0)
        if y_in_edge_domain:
            upper_ind, lower_ind = (i,j) if (points[i,1]-y) > 0 else (j,i)
            upper_x = points[upper_ind, 0] 
            upper_y = points[upper_ind, 1]
            lower_x = points[lower_ind, 0] 
            lower_y = points[lower_ind, 1]

            # is_left_turn predicate is evaluated with: sign(cross_product(upper-lower, p-lower))
            cross_prod = (upper_x - lower_x)*(y-lower_y) - (upper_y - lower_y)*(x-lower_x)
            assert abs(cross_prod) > eps, 'Need to handle these end cases separately'
            point_is_left_of_segment = (cross_prod > 0.0)
            if point_is_left_of_segment:
                intersection_counter = intersection_counter + 1
    return (intersection_counter % 2) != 0

On the input shown in the above figure (taken from my previous answer) the call is_inside(1.5, 0.0, points, edges) will return True, whereas is_inside(1.5, 3.0, points, edges) will return False.

Note that the is_inside function above does not handle degenerate cases. I added two assertions to detect such cases (you can define any epsilon value that fits your application). In many applications this is sufficient, but if not and you encounter these end cases, they need to be handled separately. See, for example, here on robustness and precision issues when implementing geometric algorithms.

这篇关于用凹域对一组点进行三角剖分的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！