我是 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
最佳答案
正如 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);
现在只选择一个西格玛。如果您可以向量化底层的 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))
关于performance - 向量化三个 for 循环,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/21805714/