Calculate bounding polygon of alpha shape from the Delaunay triangulation

Question

Given a set of points in the plane, a notion of alpha-shape, for a given positive number alpha, is defined by finding the Delaunay triangulation and deleting any triangles for which at least one edge exceeds alpha in length. Here's an example using d3:

http://bl.ocks.org/gka/1552725

The problem is that, when there are thousands of points, simply drawing all the interior triangles is too slow for an interactive visualization, so I'd like to just find the bounding polygons. This isn't so simple, because as you can see from that example sometimes there might be two such polygons.

As a simplification, suppose some clustering has been performed so that there's guaranteed to be a unique bounding polygon for each triangulation. What's the best way to find this bounding polygon? In particular, the edges have to be ordered consistently and it must support the possibility of "holes" (think a torus or donut shape--this is expressible in GeoJSON).

Solution 1

1. Create a graph in which the nodes correspond to Delaunay triangles and in which there is a graph edge between two triangles if and only if they share two vertices.

2. Compute the connected components of the graph.

3. For each connected component, find all of the nodes that have less than three adjacent nodes (that is, those that have degree 0, 1, or 2). These correspond to the boundary triangles. We define the edges of a boundary triangle that are not shared with another triangle to be the boundary edges of that boundary triangle.

As an example, I have highlighted these boundary triangles in your example "question mark" Delaunay triangulation:

By definition, every boundary triangle is adjacent to at most two other boundary triangles. The boundary edges of boundary triangles form cycles. You can simply traverse those cycles to determine polygon shapes for the boundaries. This will also work for polygons with holes if you keep them in mind in your implementation.

Solution 2

Here is some Python code that does what you need. I modified the alpha-shape (concave hull) computation from here so that it doesn't insert inner edges (the only_outer parameter). I also made it self-contained so it doesn't depend on an outside library.

from scipy.spatial import Delaunay
import numpy as np


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 is 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.simplices:
        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

If you run it with the following test code you will get this figure:

from matplotlib.pyplot import *

# Constructing the input point data
np.random.seed(0)
x = 3.0 * np.random.rand(2000)
y = 2.0 * np.random.rand(2000) - 1.0
inside = ((x ** 2 + y ** 2 > 1.0) & ((x - 3) ** 2 + y ** 2 > 1.0) & ((x - 1.5) ** 2 + y ** 2 > 0.09))
points = np.vstack([x[inside], y[inside]]).T

# Computing the alpha shape
edges = alpha_shape(points, alpha=0.25, only_outer=True)

# Plotting the output
figure()
axis('equal')
plot(points[:, 0], points[:, 1], '.')
for i, j in edges:
    plot(points[[i, j], 0], points[[i, j], 1])
show()

Solution 3

There now exists a python package alphashape which is extremely easy to use, and can be installed by pip or conda.

The main function has similar inputs to the answer given by @Iddo Hanniel, adjusting the second positional argument would give you the desired outline. Alternatively, you could leave the seconda positional argument blank and the function would optimize that parameter for you to give you the best concave hull. Beware, the computational time is increased greatly if you let the function optimize the value.

Solution 4

I know it is a delayed answer, but the methods posted here did not work for me for various reasons.

The package Alphashape that is mentioned is generally good but its downside is that it uses Shapely's cascade_union and triangulation methods to give you the final concave hull which is super slow for large datasets and low alpha values (high precision). In this case the code posted by Iddo Hanniel (and the revision by Harold) will keep a great number of edges on the interior and the only way to dissolve them is to use the aforementioned cascade_union and triangulation from Shapely.

I generally work with complex forms and the code below works fine and it is fast (2 seconds for 100K 2D points):

import numpy as np
from shapely.geometry import MultiLineString
from shapely.ops import unary_union, polygonize
from scipy.spatial import Delaunay
from collections import Counter
import itertools


def concave_hull(coords, alpha):  # coords is a 2D numpy array

    # i removed the Qbb option from the scipy defaults.
    # it is much faster and equally precise without it.
    # unless your coords are integers.
    # see http://www.qhull.org/html/qh-optq.htm
    tri = Delaunay(coords, qhull_options="Qc Qz Q12").vertices

    ia, ib, ic = (
        tri[:, 0],
        tri[:, 1],
        tri[:, 2],
    )  # indices of each of the triangles' points
    pa, pb, pc = (
        coords[ia],
        coords[ib],
        coords[ic],
    )  # coordinates of each of the triangles' points

    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) * 0.5  # Semi-perimeter of triangle

    area = np.sqrt(
        s * (s - a) * (s - b) * (s - c)
    )  # Area of triangle by Heron's formula

    filter = (
        a * b * c / (4.0 * area) < 1.0 / alpha
    )  # Radius Filter based on alpha value

    # Filter the edges
    edges = tri[filter]

    # now a main difference with the aforementioned approaches is that we dont
    # use a Set() because this eliminates duplicate edges. in the list below
    # both (i, j) and (j, i) pairs are counted. The reasoning is that boundary
    # edges appear only once while interior edges twice
    edges = [
        tuple(sorted(combo)) for e in edges for combo in itertools.combinations(e, 2)
    ]

    count = Counter(edges)  # count occurrences of each edge

    # keep only edges that appear one time (concave hull edges)
    edges = [e for e, c in count.items() if c == 1]

    # these are the coordinates of the edges that comprise the concave hull
    edges = [(coords[e[0]], coords[e[1]]) for e in edges]

    # use this only if you need to return your hull points in "order" (i think
    # its CCW)
    ml = MultiLineString(edges)
    poly = polygonize(ml)
    hull = unary_union(list(poly))
    hull_vertices = hull.exterior.coords.xy

    return edges, hull_vertices

Solution 5

Slightly revise Hanniel's answer for 3d point case (tetrahedron).

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, 3) 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 set already
        """
        if (i, j) in edges or (j, i) in edges:
            # already added
            if only_outer:
                # if both neighboring triangles are in shape, it's not a boundary edge
                if (j, i) in edges:
                    edges.remove((j, i))
            return
        edges.add((i, j))

    tri = Delaunay(points)
    edges = set()
    # Loop over triangles:
    # ia, ib, ic, id = indices of corner points of the tetrahedron
    print(tri.vertices.shape)
    for ia, ib, ic, id in tri.vertices:
        pa = points[ia]
        pb = points[ib]
        pc = points[ic]
        pd = points[id]

        # Computing radius of tetrahedron Circumsphere
        # http://mathworld.wolfram.com/Circumsphere.html

        pa2 = np.dot(pa, pa)
        pb2 = np.dot(pb, pb)
        pc2 = np.dot(pc, pc)
        pd2 = np.dot(pd, pd)

        a = np.linalg.det(np.array([np.append(pa, 1), np.append(pb, 1), np.append(pc, 1), np.append(pd, 1)]))

        Dx = np.linalg.det(np.array([np.array([pa2, pa[1], pa[2], 1]),
                                     np.array([pb2, pb[1], pb[2], 1]),
                                     np.array([pc2, pc[1], pc[2], 1]),
                                     np.array([pd2, pd[1], pd[2], 1])]))

        Dy = - np.linalg.det(np.array([np.array([pa2, pa[0], pa[2], 1]),
                                       np.array([pb2, pb[0], pb[2], 1]),
                                       np.array([pc2, pc[0], pc[2], 1]),
                                       np.array([pd2, pd[0], pd[2], 1])]))

        Dz = np.linalg.det(np.array([np.array([pa2, pa[0], pa[1], 1]),
                                     np.array([pb2, pb[0], pb[1], 1]),
                                     np.array([pc2, pc[0], pc[1], 1]),
                                     np.array([pd2, pd[0], pd[1], 1])]))

        c = np.linalg.det(np.array([np.array([pa2, pa[0], pa[1], pa[2]]),
                                    np.array([pb2, pb[0], pb[1], pb[2]]),
                                    np.array([pc2, pc[0], pc[1], pc[2]]),
                                    np.array([pd2, pd[0], pd[1], pd[2]])]))

        circum_r = math.sqrt(math.pow(Dx, 2) + math.pow(Dy, 2) + math.pow(Dz, 2) - 4 * a * c) / (2 * abs(a))
        if circum_r < alpha:
            add_edge(edges, ia, ib)
            add_edge(edges, ib, ic)
            add_edge(edges, ic, id)
            add_edge(edges, id, ia)
            add_edge(edges, ia, ic)
            add_edge(edges, ib, id)
    return edges

View more solutions

Comments

[Timothy Shields]

This is wrong. Look at the picture in my answer. There are many boundary edges that are longer than interior edges.

[Timothy Shields]

Your answer seems to suggest that you can just start repeatedly deleting the longest Delaunay edge until you're left with a bunch of polygons. That won't work. The "question mark" shape has plenty of longer-than-most edges around its boundary, yet those should not be deleted. Additionally there are shorter-than-most edges in the interior of the shapes that should be deleted. - Maybe your answer is trying to say something different? You could add more explanation.

[juniper-]

Wouldn't it be easier to just remove all shared edges and reconnect the remaining ones to form the bounding polygons?

[Timothy Shields]

@juniper- how does that differ from what I have described? Keep in mind that the approach I have described allows the algorithm to retain the topology of the boundaries (e.g. A bubble letter O has two boundaries, one inside of the other).

[Justin Meiners]

The hard part is doing step 1 in a non-convex manner.

[Ta946]

your solution was the fastest for me to get the edges without using shapely (i keep gettings errors trying to run it), The itertools combinations was slow for me so i ended up replacing it with numpy slicing and reduced the time by almost 50%

[A. Hendry]

Can you plot a shape with this? This algorithm is giving me trouble.

[Nick Fernandez]

Thanks for the code. How would we break down the edges into separate shapes/polygons if our alpha shape returns disconnected regions?

[Iddo Hanniel]

In another answer, I added code that stitches the edges. Look at the end of this answer stackoverflow.com/a/50714300/9702190 , I believe it does what you want.