如何减少 Hartley 变换的 for 循环？

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 循环的数量来提高效率？任何与此相关的事情都会非常有帮助。

用于此二维 Hartley 变换的公式如下所示：

参考：Andrew B. Watson 和 Allen Poirson 的可分离二维离散 Hartley 变换。

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));

诚然，这是相当蛮力的，人们可能会想出一种更节省内存的解决方案。该解决方案大量利用每个数组隐含地拥有无限数量的尾随单例维度，我试图在 cmets 中注明。

与 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

Andras 和 Divakar 几乎提出了与我在写作过程中相同的方法。因此，此答案将为您提供一种加速使用规范实现的方法。

您可以通过手动指定跨越图像范围的空间坐标ndgrid 来消除使用x 和y 作为变量的嵌套for 循环对。此外，方程中的变换从 0 开始索引，但您的代码从 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);