Calculate alpha bounding polygon from Delaunay triangulation

Question

Given the many points on the plane, the concept of an alpha form for a given positive number of alpha is determined by finding the Delaunay triangulation and removing any triangles for which at least one edge exceeds alpha in length. Here is an example using d3:

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

The problem is that when there are thousands of points, just drawing all the inner triangles is too slow for interactive rendering, so I would just like to find bounding polygons. This is not so simple, because, as you can see from this example, sometimes there can be two such polygons.

As a simplification, suppose some clustering has been performed so that there is guaranteed to be a single bounding polygon for each triangulation. What is the best way to find this bounding polygon? In particular, the edges should be ordered sequentially and should support the possibility of “holes” (think about the shape of a torus or donut - this can be expressed in GeoJSON).

Answer 1

• Create a graph in which the nodes correspond to Delaunay triangles and in which there is a graph face between two triangles if and only if they separate two vertices.

• Compute related graph components.

• For each connected component, find all nodes that have less than three neighboring nodes (that is, those that have degrees 0, 1, or 2). They correspond to boundary triangles . We define the edges of a boundary triangle that are not shared with another triangle, as the boundary edges of this boundary triangle.

As an example, I highlighted these boundary triangles in your “question mark” Delan triangulation example:

By definition, each boundary triangle is adjacent to no more than two other boundary triangles. The boundaries of the boundaries of the boundary triangles form cycles. You can simply go through these loops to determine the polygon shapes for the borders. This will also work for holes polygons if you consider them in your implementation.

Answer 2

Here is the Python code that does what you need. I changed the calculation of the alpha shape (concave body) here so that it does not insert inner edges ( only_outer parameter). I also made it self-sufficient so that it does not depend on an external 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 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 

If you run it with the following test code, you will get the following digit :

 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() 

Answer 3

It turns out that TopoJSON has a merge algorithm that performs exactly this task: https://github.com/mbostock/topojson/wiki/API-Reference#merge

There even an example shows it in action: http://bl.ocks.org/mbostock/9927735

In my case, it was easy for me to generate TopoJSON data, and this library function did an excellent job of this.

Answer 4

Based on @Timothy's answer, I used the following algorithm to calculate the Delaunay triangulation boundary rings.

 from matplotlib.tri import Triangulation import numpy as np def get_boundary_rings(x, y, elements): mpl_tri = Triangulation(x, y, elements) idxs = np.vstack(list(np.where(mpl_tri.neighbors == -1))).T unique_edges = list() for i, j in idxs: unique_edges.append((mpl_tri.triangles[i, j], mpl_tri.triangles[i, (j+1) % 3])) unique_edges = np.asarray(unique_edges) ring_collection = list() initial_idx = 0 for i in range(1, len(unique_edges)-1): if unique_edges[i-1, 1] != unique_edges[i, 0]: try: idx = np.where( unique_edges[i-1, 1] == unique_edges[i:, 0])[0][0] unique_edges[[i, idx+i]] = unique_edges[[idx+i, i]] except IndexError: ring_collection.append(unique_edges[initial_idx:i, :]) initial_idx = i continue # if there is just one ring, the exception is never reached, # so populate ring_collection before returning. if len(ring_collection) == 0: ring_collection.append(np.asarray(unique_edges)) return ring_collection 

Answer 5

Reconsider Hanniel's answer a bit for the 3-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 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 

Answer 6

Alpha forms are defined as delaunay triangulation with no borders exceeding alpha. First remove all inner triangles, and then all edges greater than alpha.