Kdtree求点云数据边缘点

1、MATLAB下Kdtree求点云数据边缘点。 2、dianyun.txt为点云文件，需要更换点云文件时直接替换就行。 3、直接运行dianyun.m文件即可求出点云边缘。

clc
clear

%获取点云并显示
load（‘dianyun.txt‘‘-ascii‘）%用ASCII码打开dianyun.txt点云文档
P=dianyun;%点云数据保存在P里
x=dianyun（:1）;%建立x坐标
y=dianyun（:2）;%建立y坐标
z=dianyun（:3）;%建立z坐标
scatter3（xyz）%以x，y，z描绘点云三维图像
[Eindex]=sortrows（P1）;%对点云数据进行矩阵排列
m=size（x）;%确定点的个数
x1=E（:1）;%确定x1坐标
y1=E（:2）;%确定y1坐标
z1=E（:3）;%确定z1坐标
w=[x1y1z1];%把取得的点保存在w中
mm=0;

for i=1:m
   
    Q=[];
    L=[];
    K=[];
    T=[];
    TT=[];
    k=0;
    k3=0;
    m3=0;
    
    %待判断点
    x2=w（i1）;
    y2=w（i2）;
    z2=w（i3）;
    
        %求邻域K
    for j=1:m
        x3=w（j1）;
        y3=w（j2）;
        z3=w（j3）;
        d=sqrt（（x3-x2）^2+（y3-y2）^2+（z3-z2）^2）;%计算两点距离
        if d<200&&d~=0      %根据确定邻域范围大小
            k=k+1;
            KK=[x3y3z3];
            K（k:）=KK;
        end
	end
	m2=size（K（:1））;%邻域内点的个数
    
    %邻域内点的坐标
    x4=K（:1）;
    y4=K（:2）;
    z4=K（:3）;

    %拟合平面
    P=[x4 y4 ones（length（x4）1）]\z4 ;    %最小二乘法拟合函数  
    A=P（1）;
    B=P（2）;
    C=-1;
    D=P（3）;

N=[ABC];%平面法向量
N=N/norm（N）;%向量单位化

%把点投影到拟合平面上并求出坐标
kk1=-（A*x2+B*y2+C*z2+D）;
kk2=（A^2+B^2+C^2）;
kk=-（A*x2+B*y2+C*z2+D）/（A^2+B^2+C^2）;

%求出的投影坐标
x5=x4+kk*A;
y5=y4+kk*B;
z5=z4+kk*C;
Q=[x5y5z5];

%待判断点投影坐标
x6=x2+kk*A;
y6=y2+kk*B;
z6=z2+kk*C;



%投影面上，其他投影点到待判断投影点的向量
        for k2=1:m2
        x7=Q（k21）-x6;
        y7=Q（k22）-y6;
        z7=Q（k23）-z6;
        
            if x7~=0||y7~=0||z7~=0
            	Qqq=[x7y7z7];
                k3=k3+1;
                Qq（k3:）=Qqq;
                m3=m3+1;
            end
        end
        
        %求向量的旋转角
        
        for k3=1:m3
            q1=Qq（1:）;
            q2=Qq（k3:）;
            ff1=dot（q1q2）;
            
            ff2=cross（q1q2）;            
            ff3=dot（ff2N）;
            theta=atan2（ff3ff1）;
            if theta<0
                theta=2*pi+theta;
            end
            T（k3）=theta;         
        end
        T=sort（T‘descend‘）; %从大到小排列旋转角
        
        for k4=1:m3-1
          L（k4）=T（k4）-T（k4+1）;
         
        end  
        ma=max（L）;
        mi=min（L）;
        if ma>=pi/2
                mm=mm+1;
                M（mm:）=[x2y2z2];
        end


end

x0=M（:1）;
y0=M（:2）;
z0=M（:3）;

%显示边缘点结果
figure（2）
scatter3（x0y0z0）

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     目录           0  2014-09-17 16:32  kt\
     文件        2552  2014-09-17 16:32  kt\dianyun.m
     文件        8400  2007-09-30 00:55  kt\dianyun.txt
     文件        1362  2014-05-13 20:59  kt\kdtree_ball_query.m
     文件       13824  2014-05-13 20:59  kt\kdtree_ball_query.mexw32
     文件         850  2014-05-13 20:59  kt\kdtree_build.m
     文件       28672  2014-05-13 20:59  kt\kdtree_build.mexw32
     文件         462  2014-05-13 20:59  kt\kdtree_delete.m
     文件        8704  2014-05-13 21:45  kt\kdtree_delete.mexw32
     文件           0  2014-05-13 20:59  kt\kdtree_delete.tmp1
     文件        1338  2014-05-13 20:59  kt\kdtree_k_nearest_neighbors.m
     文件       19456  2014-05-13 20:59  kt\kdtree_k_nearest_neighbors.mexw32