Delaunay三角網 alpha shape 2D點集邊緣線提取

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 Axes3D

from scipy.spatial import Delaunay
from sklearn.neighbors import KDTree


def 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, -1
    vec = (xq-xp)/L # the vector from p to q
    normal = 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, -1


def boundary_extract(points,alpha,err=10e-3):
    """
    Here, parameter err was place, because there are errors when calculating distance
    meanwhile, this err was different for different scaling 2D point cloud
    so, a parameter was placed here to considering the calculation errors
    """
    R = 1 / alpha
    pts = np.copy(points)
    tree = KDTree(pts, leaf_size=2)
    tri = Delaunay(pts)
    s = tri.simplices
    N = s.shape[0]
    i = 0
    edges = []
    centers = []
    while i <= N - 1:
        if s[i,0]==-1:
            i = i + 1
            continue
        p3 = 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] = -1
        i = i + 1
    return edges, centers


def 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 = 6
    edges, 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