如何減少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);