clc
clear

%獲取點云并顯示
load（‘dianyun.txt‘‘-ascii‘）%用ASCII碼打開dianyun.txt點云文檔
P=dianyun;%點云數(shù)據(jù)保存在P里
x=dianyun（:1）;%建立x坐標
y=dianyun（:2）;%建立y坐標
z=dianyun（:3）;%建立z坐標
scatter3（xyz）%以x，y，z描繪點云三維圖像
[Eindex]=sortrows（P1）;%對點云數(shù)據(jù)進行矩陣排列
m=size（x）;%確定點的個數(shù)
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??????%根據(jù)確定鄰域范圍大小
????????????k=k+1;
????????????KK=[x3y3z3];
????????????K（k:）=KK;
????????end
	end
	m2=size（K（:1））;%鄰域內(nèi)點的個數(shù)
????
????%鄰域內(nèi)點的坐標
????x4=K（:1）;
????y4=K（:2）;
????z4=K（:3）;

????%擬合平面
????P=[x4?y4?ones（length（x4）1）]\z4?;????%最小二乘法擬合函數(shù)??
????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