Reputation: 446
I have a list of 2D points x,y. And I need to find a smooth curve for the upper and lower edges (red and blue curves, correspondingly). See the picture below:
Here I've found a good example, where the outer edge of x,y points is detected. Using these I have work I have:
# https://stackoverflow.com/a/50714300/7200745
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 an edge 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
def find_edges_with(i, edge_set):
i_first = [j for (x,j) in edge_set if x==i]
i_second = [j for (j,x) in edge_set if x==i]
return i_first,i_second
def stitch_boundaries(edges):
edge_set = edges.copy()
boundary_lst = []
while len(edge_set) > 0:
boundary = []
edge0 = edge_set.pop()
boundary.append(edge0)
last_edge = edge0
while len(edge_set) > 0:
i,j = last_edge
j_first, j_second = find_edges_with(j, edge_set)
if j_first:
edge_set.remove((j, j_first[0]))
edge_with_j = (j, j_first[0])
boundary.append(edge_with_j)
last_edge = edge_with_j
elif j_second:
edge_set.remove((j_second[0], j))
edge_with_j = (j, j_second[0]) # flip edge rep
boundary.append(edge_with_j)
last_edge = edge_with_j
if edge0[0] == last_edge[1]:
break
boundary_lst.append(boundary)
return boundary_lst[0]
#generating of random points
N = 1000
r = 1 - 2*np.random.random((N,2))
r_norm = np.linalg.norm(r, axis=1)
points = r[r_norm <= 1]
plt.figure(figsize=(10,10))
plt.scatter(points[:,0], points[:,1], color='k', s=1)
# Computing the alpha shape
edges = alpha_shape(points, alpha=1, only_outer=True)
#order edges
edges = stitch_boundaries(edges)
plt.axis('equal')
edge_points = np.zeros((len(edges),2))
k=0
for i, j in edges:
edge_points[k,:] = points[[i, j], 0][0] , points[[i, j], 1][0]
k += 1
plt.plot(edge_points[:,0],edge_points[:,1])
#theoretical/expected edges
# xx = np.linspace(-1,1, 100)
# yy_upper = np.sqrt(1 - xx**2)
# yy_lower = -np.sqrt(1 - xx**2)
# plt.plot(xx, yy_upper, 'r:')
# plt.plot(xx, yy_lower, 'b:')
Right now the cloud of points is black. Blue line is obtained from the algorithm above.
UPDATE: the starting and final points can be chosen the most left point (OR by hand it is no problem) I am expecting the following result:
Upvotes: 3
Views: 226
Reputation: 2056
You can use the fact (see the scipy.spatial.Delaunay
documentation) that
"for 2-D, the triangle points are oriented counterclockwise".
Therefore, the outer edge points constructed in the Alpha shape will always be oriented counterclockwise i.e., the inner side of the shape will be to their left (if this was not certified in the documentation we could have checked it ourselves and flipped if necessary, but here there is no need).
This means that the points along the polygon between the leftmost point and the rightmost point are the lower hull, and the points between the rightmost and the leftmost are the upper hull. The code below implements this idea.
min_x_ind = np.argmin(edge_points[:, 0])
max_x_ind = np.argmax(edge_points[:, 0])
if min_x_ind < max_x_ind:
lower_hull = edge_points[min_x_ind:max_x_ind+1, :]
upper_hull = np.concatenate([edge_points[max_x_ind:, :], edge_points[:min_x_ind+1, :]])
else:
upper_hull = edge_points[max_x_ind:min_x_ind+1, :]
lower_hull = np.concatenate([edge_points[min_x_ind:, :], edge_points[:max_x_ind+1, :]])
and the result can be visualized for example with:
plt.plot(upper_hull[:,0],upper_hull[:,1], "r", lw=3)
plt.plot(lower_hull[:,0],lower_hull[:,1], "b", lw=3)
Upvotes: 1