MATLAB Alpha Shape算法检测边缘点

α-shapes的具体算法流程如下：

①从点集PS中任意选择一点P1开始，将与之距离小于2×α的点组成的新的点集PS2 如示意图4.3左图中的所有点，从PS2中任意选取一点P2 ，利用公式（4.1）和（4.2）求出过P1、P2 点的圆心P12。

②遍历点集PS2，依次求出（P1、P2 除外）其它点到P12的距离D：

a．如果求出的所有D均大于半径α，如示意图4.3（a）中的上方两圆，其内没有其他点，则可以判断出点P1和P2 是边缘轮廓点，P1 P2 是边界线段；

b．如果当求出的D中有小于半径α的值时，如图4.3（a）中的下方两圆，在圆内有其他点，则可以判断不是边缘轮廓点，停止遍历，转到③。

③选择PS2中的下一个点按步骤①、②进行判断，直到PS2中的全部点全部判断结束。

④取PS 中下一个点按步骤①、②、③进行判断，直到PS中的全部点判断结束。

当上述步骤运算结束后，离散点集PS的边缘轮廓点就被检测出来了，如图4.3（b）可以看出用α-shapes方法不光可以检测出凸多边形的边缘轮廓点，还可以检测凹多边形的边缘轮廓点，对有空洞的建筑物也可以进行边缘轮廓点的检测。

在我的程序中Y表示横坐标，Z表示纵坐标。有兴趣的专家可以稍作修改后参考看看。

% Alpha Shape算法检测边缘点
data=load('matlab.mat');
% data=LASreadAll(points);%点云读取
alpha=0.05;%输入alpha值
Y=(data.Pp(1,:));
Z=(data.Pp(2,:));
S1=[Y;Z];%点集S1
n=length(Y);
d=zeros(n);
ys=zeros(n);
zs=zeros(n);
for i=1:n
    for j=1:n
        p=S1(:,i)-S1(:,j);
        d(i,j)=sqrt(p(1)^2+p(2)^2);
        if d(i,j)<2*alpha
            ys(i,j)=Y(j);
            zs(i,j)=Z(j);
        end
    end
end
clear i j
a=cell(1,n);b=cell(1,n);e=[];f=[];
for i=1:n
    a{i}=(0);
    b{i}=(0);
    e=ys(i,:);
    f=zs(i,:);
    e(e==0)=[];
    f(f==0)=[];
    a{i}=a{i}+e;
    b{i}=b{i}+f;
end
y=[];z=[];D=[];y0=0;z0=0;
for i=1:n
    y=a{i};
    z=b{i};
    S2=[y;z];
    m1=length(y);
    for j1=1:m1
        if S2(:,j1)==[Y(i);Z(i)]
            S2(:,j1)=[];
            clear y z
            y=S2(1,:);
            z=S2(2,:);
            m2=length(y);
            break;
        end
    end
    for j2=1:m2
        ss=(Y(i)-y(j2))^2+(Z(i)-z(j2))^2;
        H=sqrt(alpha*alpha/ss-0.25);
        y0=Y(i)+0.5*(y(j2)-Y(i))+H*(z(j2)-Z(i));
        z0=Z(i)+0.5*(z(j2)-Z(i))+H*(Y(i)-y(j2));
        D=zeros(1,m2);
        for  k=1:m2
            D(k)=sqrt((y(k)-y0)^2+(z(k)-z0)^2);
        end
        if min(D)>=alpha
            Y(i)=smooth(Y(i));
            Z(i)=smooth(Z(i));
            plot(Y(i),Z(i),'g','LineWidth',3)
            hold on
            break;
%         else plot(Y(i),Z(i),'b.')
        end
    end
end
axis equal;