2

我想加快我的代码。我总是使用矢量化。但在这段代码中，我不知道如何避免for-loop. 我真的很感激提示如何进行。

非常感谢你的时间。

close all
clear
clc  

 % generating sample data
 x = linspace(10,130,33);
 y = linspace(20,100,22);

 [xx, yy] = ndgrid(x,y);

 k = 2*pi/50;
 s = [sin(k*xx+k*yy)];    

 % generating query points
xi  = 10:5:130;
yi  = 20:5:100;

 [xxi, yyi] = ndgrid(xi,yi);

P   = [xxi(:), yyi(:)];

% interpolation algorithm
 dx = x(2) - x(1);
dy = y(2) - y(1);
x_  = [x(1)-dx  x x(end)+dx x(end)+2*dx];
y_  = [y(1)-dy  y y(end)+dy y(end)+2*dy];   

s_  = [s(1)        s(1,:)    s(1,end)   s(1,end)
       s(:,1)      s         s(:,end)   s(:,end)
       s(end,1)    s(end,:)  s(end,end) s(end,end) 
       s(end,1)    s(end,:)  s(end,end) s(end,end)];


si   = P(:,1)*0;

M = 1/6*[-1     3   -3      1
         3 -6 3 0
     -3 0 3 0
     1 4 1 0];

tic
 for nn = 1:numel(P(:,1))
     u              = mod(P(nn,1)- x_(1), dx)/dx;     
     jj             = floor((P(nn,1) - x_(1))/dx) + 1;

     v              = mod(P(nn,2)- y_(1), dy)/dy;     
     ii             = floor((P(nn,2) - y_(1))/dy) + 1;

     D              = [s_(jj-1,ii-1)     s_(jj-1,ii)  s_(jj-1,ii+1)    s_(jj-1,ii+2)
                       s_(jj,ii-1)       s_(jj,ii)    s_(jj,ii+1)      s_(jj,ii+2)
                       s_(jj+1,ii-1)     s_(jj+1,ii)  s_(jj+1,ii+1)    s_(jj+1,ii+2)
                       s_(jj+2,ii-1)     s_(jj+2,ii)  s_(jj+2,ii+1)    s_(jj+2,ii+2)];

     U              = [u.^3 u.^2 u 1];
     V              = [v.^3 v.^2 v 1];

     si(nn)         = U*M*D*M'*V';

 end
 toc    

    scatter3(P(:,1), P(:,2), si)
    hold on
    mesh(xx,yy,s)

这是完整的示例，是二维空间中的三次 B 样条曲面插值算法。

4

0 回答 0