Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

How to eliminate for loop in energy calculation?

I'm making an algorithm which extracts information from images - most of these informations are energy-like data. This is basically done by running through the image with a kernel (size is given as parameter), and get the squared sums of values in that kernel.

This is done in three scales, where the three kernel (patch) sizes are (at this point): smallestPatchSize, smallestPatchSize * 3, smallestPatchSize * 9 (with overlapping kernels in the second and third case). This is done for several color channels, gradient filters, etc (17 of these altogether).

My question is whether it is possible to vectorize the code below; obviously this part of the code requires much more time to run than any other. I'm quite a beginner in Matlab and still try to get the hang of vectorization but this one beats me :)

for dim = 1:17

for i = 1:smallestPatchSize:(size(image,1) - currentPatchSize)
    for j = 1:smallestPatchSize:(size(image,2) - currentPatchSize)

        % calculate the position in the energy matrix
        % which has different dimensions than the input pictures
        iPosition = (i - 1 + smallestPatchSize) / smallestPatchSize;
        jPosition = (j - 1 + smallestPatchSize) / smallestPatchSize;

        % calculate the energy values and save them into the energy matrix
        energy(iPosition, jPosition, dim) = sum(sum(abs(...
            filters(i:i+currentPatchSize, j:j+currentPatchSize,dim)))) ^ 2;
    end
end
end

Thanks in advance - this is my first question @ StackOverflow :)

like image 846
Oszkar Avatar asked Oct 12 '22 11:10

Oszkar


2 Answers

You're taking local sums-of-values of blocks. A sum-filter is separable, i.e. if you want to do a sum over m-by-n, you can break it down into taking first the m-by-1 sum, and then the 1-by-n sum on the result. Note that with a full convolution you'll take the sum a little bit more often, but it's still faster than using blkproc or the more recent blockproc.

Thus, for even smallestPatchSize, you can write your inner loop as:

tmp = conv2(...
            conv2(abs(filter),ones(currentPatchSize,1),'same'),...
      ones(1,currentPatchSize,'same').^2;
%# tmp contains the energy for a sliding window filter.
energy = tmp((currentPatchSize/2+1):(currentPatchSize+1):size(image,1)-(currentPatchSize/2+1),...
      (currentPatchSize/2+1):(currentPatchSize+1):size(image,2)-(currentPatchSize/2+1));
like image 63
Jonas Avatar answered Oct 25 '22 05:10

Jonas


If you have access to the image processing toolkit you can use blkproc:

B = blkproc(A,[m n],fun)

In your case:

[m,n] = size(filters) ./ patchSize;
fun = @(x) sum( sum(x^2) ); % or whatever

B should end up being the right size. You can also use overlapping regions if you want.

like image 39
Alex Avatar answered Oct 25 '22 05:10

Alex