局部最小二乘曲面拟合循环中出现的问题
function MRsosuo()%
%
load xyz10.dat
PT=xyz10;
R=10;
x=PT(:,1);
y=PT(:,2);
z=PT(:,3);
minx=min(x);
miny=min(y);
maxx=max(x);
maxy=max(y);
dx=maxx-minx
dy=maxy-miny
m=fix(dx/R)+1
n=fix(dy/R)+1
for l=0:n
for r=0:m
k=1;
q=1;
h=1;
for i=1:length(PT)
rr(k)=sqrt(((minx+r*R)-x(i))^2+((miny+l*R)-y(i))^2);
if rr(k)<R
a(q)=x(i);
b(q)=y(i);
c(q)=z(i);
ch(q)=i/i;
q=q+1;
else
e(h)=x(i);
f(h)=y(i);
g(h)=z(i);
h=h+1
end
k=k+1;
end
rr
xy=a.*b;%---------问题出现的地方 当rr(k)<R对于所有点都不成立时 a,b,c,ch在前次循环结束时就给释放了
x2=a.^2;%--------此时程序若不需要继续此次循环直接到下一个循环，如何进行？？
y2=b.^2;
N=; %拟合
N1=N'*N;
N2=pinv(N1);
M=N2*N'*c' ;
ZN=M'*N' ;
for i=1:length(a)
CC(i)=ZN(i)-c(i);
end
CC; % 残差＝拟合的高程－观测值
%----------------------------------------------------------------------
%分离出地面点和地物点
%判断若点在拟合面上，则认为是地物点，否则为地面点
j=1;
p=1;
dwPT(j)=0;
dmPT(p)=0;
for i=1:length(CC)
if CC(i)<0
dwPT(j)=i; %地物点所在的行号
j=j+1;
else
dmPT(p)=i; %地面点所在的行号
p=p+1;
end
end
dwPT
dmPT
%挑出的地物点
if dwPT~=0
for i=1:length(dwPT)
dw_m=dwPT(i);
dwPTx(i)=a(dw_m);
dwPTy(i)=b(dw_m);
dwPTz(i)=c(dw_m);
CCdw(i)=CC(dw_m);
end
else
disp('dwPT=NULL');
end
%挑出的地面点
if dmPT~=0
for i=1:length(dmPT)
dm_n=dmPT(i);
dmPTx(i)=a(dm_n);
dmPTy(i)=b(dm_n);
dmPTz(i)=c(dm_n);
CCdm(i)=CC(dm_n);
end
else
disp('dmPT=NULL');
end
% 打开文件将拟合的数据写入
fid = fopen('out.dat','wt');
for i=1:length(a)
fprintf(fid,'%10.3f\t %10.3f\t %10.3f\t %10.3f\t %10.3f\t\n',...
a(i),b(i),c(i),ZN(i),CC(i));
end
fprintf(fid,'\n');
fclose(fid);
%将地物点数据存入dwPT.dat文件中
fid = fopen('dwPT.dat','wt');
if dwPT~=0
for i=1:length(dwPT)
fprintf(fid,'%10.3f\t %10.3f\t %10.3f\t %10.3f\t\n',...
dwPTx(i),dwPTy(i),dwPTz(i),CCdw(i));
end
fprintf(fid,'\n');
fclose(fid);
else
fprintf(fid,'dwPT file is NULL');
end
%将地面点数据存入dmPT.dat文件中
fid = fopen('dmPT.dat','wt');
if dmPT~=0
for i=1:length(dmPT)
fprintf(fid,'%10.3f\t %10.3f\t %10.3f\t %10.3f\t\n',...
dmPTx(i),dmPTy(i),dmPTz(i),CCdm(i));
end
fprintf(fid,'\n');
fclose(fid);
else
fprintf(fid,'dmPT file is NULL');
end
clear a;
clear b;
clear c;
clear ch;
clear e;
clear f;
clear g;
clear CC;
clear rr;
clear dwPT;
clear dmPT;
end
end
end 原帖由 yanshutianshi 于 2008-3-5 12:55 发表 http://chinavib.com/forum/images/common/back.gif
function MRsosuo()
%
%
load xyz10.dat
PT=xyz10;
R=10;
x=PT(:,1);
y=PT(:,2);
z=PT(:,3);
minx=min(x);
miny=min(y);
maxx=max(x);
maxy=max(y);
dx=maxx-minx
...
作个条件判断不行吗？你的代码这么长，除非有专稿这个的版友路过，否则根本没有人愿意看的。最好换个简单容易理解的例子，看看会员守则的2楼吧 function MRsosuo1()
%
%
load FZ.dat
PT=FZ;
R=50;
x=PT(:,1);
y=PT(:,2);
z=PT(:,3);
minx=min(x);
miny=min(y);
maxx=max(x);
maxy=max(y);
dx=maxx-minx;
dy=maxy-miny;
m=fix(dx/R)+1
n=fix(dy/R)+1
for l=0:n
for r=0:m
k=1;
q=1;
a=0;
b=0;
c=0;
ch=0;
for i=1:length(PT)
rr(k)=sqrt(((minx+r*R)-x(i))^2+((miny+l*R)-y(i))^2);
if rr(k)<R
a(q)=x(i);
b(q)=y(i);
c(q)=z(i);
ch(q)=1;
q=q+1;
end
k=k+1;
end
xy=a.*b;
x2=a.^2;
y2=b.^2;
N=; %平面拟合
N1=N'*N;
N2=pinv(N1); %
M=N2*N'*c' ; % ZN=M'*N' ;
for i=1:length(a)
CC(i)=ZN(i)-c(i);
end
CC; % 残差＝拟合的高程－观测值
%----------------------------------------------------------------------
%分离出地面点和地物点
%判断若点在拟合面上，则认为是地物点，否则为地面点
j=1;
p=1;
dwPT(j)=0;
dmPT(p)=0;
for i=1:length(CC)
if CC(i)<0
dwPT(j)=i; %地物点所在的行号
j=j+1;
else
dmPT(p)=i; %地面点所在的行号
p=p+1;
end
end
%挑出的地物点
if dwPT~=0
for i=1:length(dwPT)
dw_m=dwPT(i);
dwPTx(i)=a(dw_m);
dwPTy(i)=b(dw_m);
dwPTz(i)=c(dw_m);
CCdw(i)=CC(dw_m);
end
end
%挑出的地面点
if dmPT~=0
for i=1:length(dmPT)
dm_n=dmPT(i);
dmPTx(i)=a(dm_n);
dmPTy(i)=b(dm_n);
dmPTz(i)=c(dm_n);
CCdm(i)=CC(dm_n);
end
end
% 打开文件将拟合的数据写入
fid = fopen('out.dat','at');
for i=1:length(a)
fprintf(fid,'%10.3f\t %10.3f\t %10.3f\t %10.3f\t %10.3f\t\n',...
a(i),b(i),c(i),ZN(i),CC(i));
end
fprintf(fid,'\n');
fclose(fid);
%将地物点数据存入dwPT.txt文件中
fid = fopen('dwPT.dat','at');
if dwPT~=0
for i=1:length(dwPT)
fprintf(fid,'%10.3f\t %10.3f\t %10.3f\t %10.3f\t\n',...
dwPTx(i),dwPTy(i),dwPTz(i),CCdw(i));
end
fprintf(fid,'\n');
fclose(fid);
end
%将地面点数据存入dmPT.txt文件中
fid = fopen('dmPT.dat','at');
if dmPT~=0
for i=1:length(dmPT)
fprintf(fid,'%10.3f\t %10.3f\t %10.3f\t %10.3f\t\n',...
dmPTx(i),dmPTy(i),dmPTz(i),CCdm(i));
end
fprintf(fid,'\n');
fclose(fid);
end
clear a;
clear b;
clear c;
clear ch;
clear CC;
clear rr;
clear dwPT;
clear dmPT;
end
end
end
经过修改后的 效率要稍好些 完全可以运行
页:
[1]