为给定的点云拆分上边缘线和下边缘线？
Split the upper edge line and lower edge line for a given cloud of points?
我有一个二维点 x,y 的列表。我需要为上边缘和下边缘找到一条平滑的曲线（相应的红色和蓝色曲线）。
见下图：
我找到了一个很好的例子，其中检测到 x,y 点的外边缘。
使用这些我有我的工作：
#
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:')
现在点云是黑色的。蓝线是从上面的算法中得到的。
更新：
起点和终点可以选择最左边的点（或者手动也没有问题）
我期待以下结果：
您可以使用以下事实（参见 scipy.spatial.Delaunay documentation）
“对于二维，三角形点逆时针方向”。
因此，在 Alpha 形状中构造的外边缘点将始终逆时针方向，即形状的内侧将在它们的左侧（如果文档中没有证明这一点，我们可以自己检查并在必要时翻转，但是这里不需要）。
这意味着最左边点和最右边点之间的多边形上的点是下层船体，最右边和最左边之间的点是上层船体。下面的代码实现了这个想法。
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, :]])
结果可以可视化，例如：
plt.plot(upper_hull[:,0],upper_hull[:,1], "r", lw=3)
plt.plot(lower_hull[:,0],lower_hull[:,1], "b", lw=3)
我有一个二维点 x,y 的列表。我需要为上边缘和下边缘找到一条平滑的曲线（相应的红色和蓝色曲线）。
见下图：
#
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:')
现在点云是黑色的。蓝线是从上面的算法中得到的。
更新：
起点和终点可以选择最左边的点（或者手动也没有问题）
我期待以下结果：
您可以使用以下事实（参见 scipy.spatial.Delaunay documentation）
“对于二维，三角形点逆时针方向”。
因此，在 Alpha 形状中构造的外边缘点将始终逆时针方向，即形状的内侧将在它们的左侧（如果文档中没有证明这一点，我们可以自己检查并在必要时翻转，但是这里不需要）。
这意味着最左边点和最右边点之间的多边形上的点是下层船体，最右边和最左边之间的点是上层船体。下面的代码实现了这个想法。
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, :]])
结果可以可视化，例如：
plt.plot(upper_hull[:,0],upper_hull[:,1], "r", lw=3)
plt.plot(lower_hull[:,0],lower_hull[:,1], "b", lw=3)