How to reduce for loops for Hartley transform?

Question

I am currently working on 2D Hartley transform. The code is shown below:

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

It has 4 for loops and it takes a very long time to execute this. Is there any method to make it more efficient by reducing the amount of for loops? Anything regarding this would be very helpful.

The formula used for this 2-D Hartley transform is shown below:

Reference:Separable two-dimensional discrete Hartley transform by Andrew B. Watson and Allen Poirson.

score 3 · Answer 1 · answered Dec 14 '15 at 21:14

If you can fit into memory, you can use bsxfun with some additional singleton dimensions:

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

Admittedly this is quite brute-force, one could probably come up with a more memory-efficient solution. The solution heavily utilizes that every array implicitly possesses an infinite number of trailing singleton dimensions, which I tried to note in the comments.

Comparison to your original looping version with 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

which means they give the same solution within machine precision.

rayryeng · Answer 2 · 2015-12-14T21:27:56.520

Andras and Divakar pretty much suggested the same method that I was in the middle of writing. As such, this answer will simply provide you with a method for speeding up using your canonical implementation.

You can eliminate the pair of nested for loops that use x and y as variables by manually specifying a ndgrid of spatial coordinates that span the extent of the image. Also, the transform in your equation starts indexing from 0, yet your code is starting from 1. You'll need to start indexing at 1 with MATLAB, but the compute the terms in the equation, you need to start from 0. As such, when you compute a and b, you'll have to subtract 1 from x, y, u and v.

In any case, you can then compute element-wise products and sum all of the values together. It's also a good idea to pre-allocate the output for efficiency:

%// 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

You will hopefully get better performance improvements. Compare this one with Andras's or Divakar's solution and see which one you're comfortable in using.

@Divakar because rayryeng fixed the semantic problem: the formula sums from 0, while OP sums from 1. So our solution gives identical results, rayryeng's solution probably gives a correct one;) — Andras Deak -- Слава Україні, Dec 14 '15 at 21:47
@AndrasDeak Ah yes, didn't notice that in the formula, thanks! :) — Divakar, Dec 14 '15 at 21:48
@Divakar yup OP needs to sum from 0 not 1. Thanks for checking! — rayryeng, Dec 14 '15 at 21:52

score 3 · Answer 3 · edited May 23 '17 at 12:32

3

Here's one using bsxfun and fast matrix-multiplication for a vectorized solution -

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

edited May 23 '17 at 12:32

Community

1
1

answered Dec 14 '15 at 21:17

Divakar

218,885
19
262
358

I knew it!;) Efficient as always. – Andras Deak -- Слава Україні Dec 14 '15 at 21:18
@AndrasDeak Well the bottleneck seems to be the cosine and sine calculations, too bad really. – Divakar Dec 14 '15 at 22:25
@AndrasDeak .. and the *co-partner* – Divakar Dec 14 '15 at 22:51
All the codes provided work very well on a small set of values, but when I try to apply on an image of size 480 x 640, bsxfun gives me an error which says 'array exceeds maximum array size preference'. How do I go about with this? – shreyas kamath Dec 15 '15 at 02:12
1

@shreyaskamath For memory exceed error, you might want to go with a partially vectorized solution, something like suggested in [`Ray's solution`](http://stackoverflow.com/a/34276759/3293881). – Divakar Dec 15 '15 at 06:14

How to reduce for loops for Hartley transform?

3 Answers3