向量化三个 for 循环

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

正如 Dennis 所指出的，这是很多代码，将其减少到探查器给出的速度较慢的最小值会有所帮助。我不确定我的代码是否与您的代码等效，您可以尝试一下并对其进行概要分析吗？ 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);

现在只选择一个 sigma。如果您可以对底层的 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))