局部最小二乘曲面拟合循环中出现的问题

yanshutianshi · 发表于 2008-3-5 12:55

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？我要加入

x

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=[ch' a' b' xy' x2' y2']; %拟合
         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

eight · 发表于 2008-3-5 15:49

原帖由 yanshutianshi 于 2008-3-5 12:55 发表

 登录/注册后可看大图

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楼吧

yanshutianshi · 发表于 2008-4-10 10:45

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=[ch' a' b' xy' x2' y2']; %平面拟合
         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

经过修改后的效率要稍好些完全可以运行

账号		自动登录	找回密码
密码			我要加入

[综合讨论] 局部最小二乘曲面拟合循环中出现的问题

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

本帖被以下淘专辑推荐:

评分