1\

1. Delaunay三角网
   参考blog:
   Scipy
   笔记
   [Geometry] Alpha Shapes - 原理及我的实现
2. Alpha Shape
   Widyaningrum E , Peters R Y , Lindenbergh R C . Building outline extraction from als point clouds using medial axis transform descriptors[J]. Pattern Recognition, 2020:107447.
   

3. Codes

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3Dfrom scipy.spatial import Delaunay
from sklearn.neighbors import KDTreedef plot_circle(centers,rs,ax):N = centers.shape[0]for i in range(N):theta = np.arange(0, 2*np.pi, 0.01)x = centers[i,0] + rs[i] * np.cos(theta)y = centers[i,1] + rs[i] * np.sin(theta)ax.plot(x, y, 'b-', alpha=0.1)def edge_check_vaild(e,tree,r,err):xp = e[0]xq = e[1]L = np.sqrt(np.dot(xq-xp,xq-xp))if L > 2*r:return False, -1vec = (xq-xp)/L # the vector from p to qnormal = np.array([vec[1],-vec[0]])c1 = (xp + xq) / 2 + normal * np.sqrt(r**2-(L/2)**2)c2 = (xp + xq) / 2 - normal * np.sqrt(r**2-(L/2)**2)c = np.array([c1,c2])count = tree.query_radius(c,r=r+err,return_distance=False,count_only=True,sort_results=False)if count[0]<=2:return True, c[0]elif count[1]<=2:return True, c[1]else:return False, -1def boundary_extract(points,alpha,err=10e-3):"""Here, parameter err was place, because there are errors when calculating distancemeanwhile, this err was different for different scaling 2D point cloudso, a parameter was placed here to considering the calculation errors"""R = 1 / alphapts = np.copy(points)tree = KDTree(pts, leaf_size=2)tri = Delaunay(pts)s = tri.simplicesN = s.shape[0]i = 0edges = []centers = []while i <= N - 1:if s[i,0]==-1:i = i + 1continuep3 = s[i]e1 = np.array([points[p3[0],:],points[p3[1],:]])e2 = np.array([points[p3[1],:],points[p3[2],:]])e3 = np.array([points[p3[0],:],points[p3[2],:]])e = [e1,e2,e3]for j in range(3):flag, center = edge_check_vaild(e[j],tree,R,err)if flag:edges.append(e[j])centers.append(center)nb = tri.neighbors[i]nb_valid = nb[nb!=-1]#nb_valid_num = nb_valid.shape[0]#s[nb_valid] = -1i = i + 1return edges, centersdef show_edge(edges,points,circle=None,r=None):fig = plt.figure()ax = fig.add_subplot(111)ax.scatter(*zip(*points),s=4,c='k')for i in range(len(edges)):ax.plot(*zip(*edges[i]),'-r')if circle is not None:plot_circle(circle,r,ax)plt.show()if __name__ == "__main__":pts = np.random.rand(200, 2) # 随机生成10个2维点alpha = 6edges, centers = boundary_extract(pts,alpha,err=10e-5)show_edge(edges,pts,circle=np.array(centers),r=np.ones(len(centers))/alpha)print("over!!!")

4. Results