如何减少Hartley变换的循环？

Question

我目前正在研究2D Hartley变换。代码如下所示：

for u=1:size(img,1) 
    for v=1:size(img,2) 
     for x=1:size(img,1) 
      for y=1:size(img,2) 
       a = 2*pi*u*x/size(img,1); 
       b = 2*pi*v*y/size(img,2); 
       temp= img(x,y)*(cos(a+b) + sin(a+b)); 
       X(u,v)= X(u,v)+temp; 
      end 
     end 
    end 
end 

它有4个for循环，它需要很长的时间来执行这一点。有没有什么方法可以通过减少for循环的数量来提高效率？任何关于此的将会非常有帮助。

用于该2-d哈特利的公式变换如下：

参考：可分二维离散Hartley由Andrew B.沃森和Allen Poirson变换。

Answer 1

如果能够装入内存，则可以使用bsxfun一些额外的单维度：

N = size(img,1); 
M = size(img,2); 
x = [1:N].'; %' vector of size N x 1 (x 1 x 1) 
y = 1:M;  % vector of size 1 x M (x 1 x 1) 
u = permute(1:N,[1 3 2]); %vector of size 1 x 1 x N (x 1) 
v = permute(1:M,[1 3 4 2]); %vector of size 1 x 1 x 1 x M 

a = 2*pi/N*bsxfun(@times,u,x); % N x 1 x N x 1 
b = 2*pi/M*bsxfun(@times,v,y); % 1 x M x 1 x M 
apb = bsxfun(@plus,a,b); % N x M x N x M 
%img is N x M (x 1 x 1) 

X2 = squeeze(sum(sum(bsxfun(@times,img,cos(apb)+sin(apb)),1),2)); 

诚然，这是一个相当强力，一个也许可以想出一个内存更有效的解决方案。该解决方案很大程度上利用了每个数组隐含地拥有无限数量的尾随单身维度，我在评论中试图说明这一点。

比较你原来的循环版本N=20; M=30; img=rand(N,M);：

>> max(max(abs(X-X2))) 
ans = 
    1.023181539494544e-12 
>> max(max(abs(X))) 
ans = 
    3.091143465722029e+02 

，这意味着他们给机器精度内相同的解决方案。

Answer 2

安德拉斯和Divakar几乎认为我是在写作中同样的方法。因此，这个答案只会为您提供一种加速使用规范实现的方法。

可以消除一对通过手动指定跨越图像的程度的空间坐标的一个ndgrid使用x和y作为变量的嵌套for循环。此外，等式中的变换从开始索引，但是您的代码从1开始。您需要使用MATLAB开始索引1，但计算方程中的项需要从0开始。因此，当你计算a和b，你必须从x，y，u和v减去1。

在任何情况下，你就可以计算元素方面的产品，总结所有的值一起。这也是一个好主意，预先分配输出效率：

%// Change - preallocate 
X = zeros(size(img)); 
%// New - define spatial coordinates 
[x,y] = ndgrid(0:size(img,1)-1, 0:size(img,2)-1); 
for u=1:size(img,1) 
    for v=1:size(img,2) 
     %// Change 
     a = 2*pi*(u-1)*x/size(img,1); 
     b = 2*pi*(v-1)*y/size(img,2); 
     temp = img.*(cos(a+b) + sin(a+b)); 
     %// Change 
     X(u,v) = sum(temp(:));   
    end 
end 

你会希望得到更好的性能改进。将这一个与Andras或Divakar的解决方案进行比较，看看您使用哪一种解决方案。

Answer 3

下面是使用一个bsxfun和fast matrix-multiplication了量化的解决方案 -

%// Store size parameters 
[m,n] = size(img); 

%// Get vectorized versions of a and b 
A = 2*pi*(1:m).'*(1:m)/m; 
B = 2*pi*(1:n).'*(1:n)/n; 

%// Add vectorized a and b's to form a 4D array. Get cos + sin version. 
AB = bsxfun(@plus,permute(A,[1 3 2 4]),permute(B,[3 1 4 2])); 
cosAB = cos(AB) + sin(AB); 

%// Finally bring in magic of matrix-multiplication for final output 
Xout = reshape(img(:).'*reshape(cosAB,m*n,[]),m,n);