2

I want to speed up my code. I always use vectorization. But in this code I have no idea how to avoid the for-loop. I would really appreciate a hint how to proceed.

thank u so much for your time.

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)

This is the full example and is a cubic B-spline surface interpolation algorithm in 2D space.

BhushanK
  • 1,205
  • 6
  • 23
  • 39
PHS
  • 21
  • 2
  • you can probably do this with `bsxfun` \ `im2col` \ `permute` etc but that doesn't mean that the for-loop solution is slower. – bla Dec 17 '15 at 20:11
  • Don't know what you mean exactly. I expanded the U*M*D*M'*V' term. Then it's easy to do it vectorized. But for small sets of interpolation points its slower, than the loop. For large sets (10k) it gets faster than the loop. Unfortunately, I have a lot smaller sets. – PHS Dec 18 '15 at 08:57

0 Answers0