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:

Take the centroid of each polygon as C_i = (x_i,y_i).
Then, imagine a line L = [C_i,(+inf,y_i)]: that is, a line that spans east past the end of your domain.
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.

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.

Do you have access to the domain, or do you have to find out what it is? — Mauricio, Jul 21 '17 at 00:11
Are you familiar with alpha hulls / alpha shapes? https://en.wikipedia.org/wiki/Alpha_shape — Rethunk, Jul 21 '17 at 01:57
I would use the marching square algorithm for any convex polygon. It was made to find borders easily — asdfasdf, Feb 06 '18 at 19:02

score 5 · Accepted Answer · answered Jun 28 '18 at 11:41

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.

score 1 · Answer 2 · answered Jul 23 '17 at 11:11

One of Classic DT algorithms generates first a bounding triangle, then adds all new triangles sorted by x, then prunes out all triangles having a vertex in the supertriangle.

At least from the provided image one can derive the heuristics of pruning out also some triangles having all vertices on the concave hull. Without a proof, the triangles to be pruned out have a negative area when their vertices are sorted in the same order as the concave hull is defined.

This may need the concave hull to be inserted as well, and to be pruned out.

score 1 · Answer 3 · answered Jun 14 '18 at 16:08

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:

Take the centroid of each polygon as C_i = (x_i,y_i).
Then, imagine a line L = [C_i,(+inf,y_i)]: that is, a line that spans east past the end of your domain.
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.

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.

score 0 · Answer 4 · edited Feb 20 '18 at 15:54

0

You can try a constrained delaunay algorithm for example with sloan algoritm or cgal library.

[1] A Brute-Force Constrained Delaunay Triangulation?

edited Feb 20 '18 at 15:54

Ripi2

6,098
1
12
28

answered Jul 23 '17 at 10:13

Gigamegs

12,342
7
31
71

score 0 · Answer 5 · answered Jun 05 '20 at 15:52

A simple but elegant way is to loop over the triangels and check wether they are within our domain or not. The shapely package could do the trick for you.

for more on this please check the following post: https://gis.stackexchange.com/a/352442 Note that triangulation in shapely is also implemented, even for MultiPoin objects.

I used it, the performance was amazing and the code was only like five lines.

score -1 · Answer 6 · answered Jul 21 '17 at 06:43

-1

Compute the triangles centroid an check if it's inside the polygon using this algorithm.

answered Jul 21 '17 at 06:43

abenci

7,370
13
61
124

Not sure that would work; consider a "pointy bit" on the polygon that might happen to contain the centroid, but not the rest of the triangle – meowgoesthedog Jul 21 '17 at 08:12
The final DT doesn't contain such a triangle, even though the intermediate tessellation can. – Aki Suihkonen Jul 23 '17 at 10:54

Triangulate a set of points with a concave domain

Setup

Question

6 Answers6

Linked