向量化三个循环

Question

我对Matlab很陌生，我需要帮助来加快我的代码的一部分。我正在写一个执行3D矩阵卷积的Matlab应用程序，但与标准卷积不同，内核不是恒定的，它需要为图像的每个像素计算。

到目前为止，我已经结束了一个工作代码，但慢得令人难以置信：

function result = calculateFilteredImages(images, T) 

% images - matrix [480,360,10] of 10 grayscale images of height=480 and width=360 
% reprezented as a value in a range [0..1] 
% i.e. images(10,20,5) = 0.1231; 

% T - some matrix [480,360,10, 3,3] of double values, calculated earlier 

    kerN = 5;    %kernel size 
    mid=floor(kerN/2);  %half the kernel size 
    offset=mid+1;   %kernel offset 
    [h,w,n] = size(images); 
    %add padding so as not to get IndexOutOfBoundsEx during summation: 
    %[i.e. changes [1 2 3...10] to [0 0 1 2 ... 10 0 0]] 
    images = padarray(images,[mid, mid, mid]); 

    result(h,w,n)=0;   %preallocate, faster than zeros(h,w,n) 
    kernel(kerN,kerN,kerN)=0; %preallocate 

    % the three parameters below are not important in this problem 
    % (are used to calculate sigma in x,y,z direction inside the loop) 
    sigMin=0.5; 
    sigMax=3; 
    d = 3; 

    for a=1:n; 
     tic; 
     for b=1:w; 
      for c=1:h; 
       M(:,:)=T(c,b,a,:,:); % M is now a 3x3 matrix 
       [R D] = eig(M); %get eigenvectors and eigenvalues - R and D are now 3x3 matrices  

       % eigenvalues 
       l1 = D(1,1); 
       l2 = D(2,2); 
       l3 = D(3,3); 

       sig1=sig(l1 , sigMin, sigMax, d); 
       sig2=sig(l2 , sigMin, sigMax, d); 
       sig3=sig(l3 , sigMin, sigMax, d); 

       % calculate kernel 
       for i=-mid:mid 
        for j=-mid:mid 
         for k=-mid:mid 
          x_new = [i,j,k] * R; %calculate new [i,j,k] 
          kernel(offset+i, offset+j, offset+k) = exp(- (((x_new(1))^2)/(sig1^2) + ((x_new(2))^2)/(sig2^2) + ((x_new(3))^2)/(sig3^2)) /2); 
         end 
        end 
       end 
       % normalize 
       kernel=kernel/sum(kernel(:)); 

       %perform summation 
       xm_sum=0; 
       for i=-mid:mid 
        for j=-mid:mid 
         for k=-mid:mid 
          xm_sum = xm_sum + kernel(offset+i, offset+j, offset+k) * images(c+mid+i, b+mid+j, a+mid+k); 
         end 
        end 
       end 
       result(c,b,a)=xm_sum; 

      end 
     end 
     toc; 
    end 
end 

我试着用

sigma=[sig1 sig2 sig3] 
[x,y,z] = ndgrid(-mid:mid,-mid:mid,-mid:mid); 
k2 = arrayfun(@(x, y, z) exp(-(norm([x,y,z]*R./sigma)^2)/2), x,y,z); 

取代“计算内核”的一部分，但它竟然是甚至比循环更慢。我浏览了几篇有关矢量化的文章和教程，但我对此一直很困惑。 它可以被矢量化或以某种方式加速使用别的东西？ 我是Matlab新手，也许有一些内置函数可以帮助在这种情况下？

更新 剖析结果：
T.mat
grayImages.mat

Answer 1

如丹尼斯指出的，这是一个很大的代码，切割它：其分析期间使用

样本数据下降到分析器给出的最低速度将会有所帮助。我不确定我的代码是否与您的代码相同，您是否可以尝试并分析它？ Matlab矢量化的“窍门”是使用。*和。^，它们逐个操作而不必使用循环。 http://www.mathworks.com/help/matlab/ref/power.html

把你重写部分：

sigma=[sig1 sig2 sig3] 
[x,y,z] = ndgrid(-mid:mid,-mid:mid,-mid:mid); 
k2 = arrayfun(@(x, y, z) exp(-(norm([x,y,z]*R./sigma)^2)/2), x,y,z); 

而且随便挑一个西格马现在。如果可以矢量化底层k2公式，则循环3个不同的sigma不是性能问题。

编辑：更改matrix_to_norm代码为x(:)，并没有逗号。见Generate all possible combinations of the elements of some vectors (Cartesian product)

然后尝试：

% R & mid my test variables 
R = [1 2 3; 4 5 6; 7 8 9]; 
mid = 5; 
[x,y,z] = ndgrid(-mid:mid,-mid:mid,-mid:mid); 
% meshgrid is also a possibility, check that you are getting the order you want 
% Going to break the equation apart for now for clarity 
% Matrix operation, should already be fast. 
matrix_to_norm = [x(:) y(:) z(:)]*R/sig1 
% Ditto 
matrix_normed = norm(matrix_to_norm) 
% Note the .^ - I believe you want element-by-element exponentiation, this will 
% vectorize it. 
k2 = exp(-0.5*(matrix_normed.^2))