Triangulate a set of points with a concave domain

Question

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.

Answer 1

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.

Answer 2

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.