Fast average calculation of submatrices in large matrix

Hi,

I have a matrix M=1024x512 elements. What I would like to do is to calculate the average of every submatrix of a size of 4x4 and assign that average to all the elements in that submatrix. Then move on to the next submatrix and do the routine all over again. I have managed to do this by two for loops, but the problem is that I have lots of matrices of this size and this routine is very slow. How can I speed up this calculations? I have the impression that the use of for loops isn't that quick at all :(.

This is the code I used:

for i = 1:4:(1024-4)
    for j = 1:4:(512-4)
        average = mean(mean(M(i:(i+4-1),j:(j+4-1))));
        M(i:(i+4-1),j:(j+4-1)) = average;
    end
end

Your help would greatly be apreciated!

Answer 1

[Andrei Bobrov]

PART 1

use function blockproc from Image Processing Toolbox

a1 = blockproc(M,[4,4],@(x)mean2(x.data))

or without blockproc, but with conv2

a1 = conv2(M,ones(4)/16,'valid');
a1 = a1(1:4:end,1:4:end);  %CORRECT

or with cellfun

a1 = cellfun(@mean2,mat2cell(M,4*ones(size(M,1)/4,1),4*ones(size(M,2)/4,1)));

average = kron(a1,ones(4));

or

average = cell2mat(arrayfun(@(x)x*ones(4),a1,'un',0));

ADD

average = conv2(M,ones(4)/16,'valid');

[Matt J]

Stijn commented:

Thanks Andrei for your answer! However this will return a matrix that is 4 times smaller than the original. How do I assign the first value of 'average' to all the elements in the 4x4 submatrix of the original matrix and the second value of 'average' to all the elements the second submatris etc...?

[Andrei Bobrov]

See ADD part in my answer or Matt J's answer.

[Stijn]

I've tested this with this matrix

M= [
1 2 3 0 1 2 5 8
1 2 3 0 1 2 5 8
1 2 3 0 1 2 5 8
1 2 3 0 1 2 5 8
1 2 3 0 1 2 5 8
1 2 3 0 1 2 5 8
1 2 3 0 1 2 5 8
1 2 3 0 1 2 5 8]

The average is calculated correctly:

average = [
1.5 4
1.5 4]

however, the substitution into M does not go well. It returns me this:

    1.5000    1.5000    1.5000    2.0000    4.0000
    1.5000    1.5000    1.5000    2.0000    4.0000
    1.5000    1.5000    1.5000    2.0000    4.0000
    1.5000    1.5000    1.5000    2.0000    4.0000
    1.5000    1.5000    1.5000    2.0000    4.0000

this should be an 8x8 matrix only with the values 1.5 and 4

[Andrei Bobrov]

see PART 1 in my answer

[Stijn]

I think I'm missing something, because it still doesn't work. Sorry for the misunderstanding.

edit: I got it working, thanks! Now check whether it is faster ;)

[Andrei Bobrov]

correct

[Stijn]

can I use blockproc to shift one element in stead of a whole block?

Answer 2

[Matt J]

This will be faster than BLOCKPROC,

    kernel=ones(1,4)/4;

    M=conv2(kernel, kernel, M,'valid');

Answer 3

[Azzi Abdelmalek]

To improve the speed of your code two times:

for i = 1:4:(1024-4)
    for j = 1:4:(512-4)
        v=M1(i:(i+4-1),j:(j+4-1));
        M1(i:(i+4-1),j:(j+4-1)) =mean(v(:));
    end
end

Answer 4

[Stijn]

I want to expand this question a bit. It is in the line of the previous problem, only I would like to do the following:

Suppose I have the following matrix:

m=[1  2  3  4  5
   6  7  8  9  10
   11 12 13 14 15
   16 17 18 19 20
   21 22 23 24 25]

What i would like to establish here is a some kind of moving average. Instead of jumping ahead with a submatrix of 3x3, I want to move only 1 element forward at the time and calculate the average again and store it a different matrix. To be more clear, suppose you pick an element (say m(2,2)), I want to get an average around that element by averaging the elements m(1:3,1:3), thus go 1 element to each side to get a average vlue of the 3x3 submatrix. (If you use 2 elements to each side you will get an average of a 5x5 submatrix). Only elements on the edges of the matrix are calculated differently.

The matrix would than be obtained like this (I hope this image clarifies it a bit)

https://dl.dropboxusercontent.com/u/23491985/movingaverage.png

[Andrei Bobrov]

m = [1    2    3    4    5
    6    7    8    9   10
   11   12   13   14   15
   16   17   18   19   20
   21   22   23   24   25];

[a b] = size(m);
c = true([a b]);
c(2:end-1,2:end-1) = false;
ii = bsxfun(@plus,[0 (b-1)*a],[1;a]);
d = c*6;
d(ii) = 4;
out = conv2(m,ones(3)/9);
out(c) = out(c)*9./d(d>0);

Answer 5

[Jan]

With the C-Mex FEX: BlockMean:

   D = rand(1024, 512);
   M = BlockMean(M, 4, 4);
   M = kron(M, ones(4, 4));