How to calculate the mean of all neighbouring elements in a matrix

Hello, I have a m x n matrix and I want tot calculate the mean for every element and all its neighbouring values (horizontal, vertical, diagonal). Here I also want to consider the number of neighbours an element has dependent on its position. With matrix M = [3,4,8,6;1,6,9,0;2,5,7,1], M(1,1) should only have 3 neighbours, M(1,2) 5 neigbours and M(2,2) 8 neighbours.

How can I implement a function to calculate the mean of all neighbouring elements with using the matrix positions i,j which represents the rows and column values of the matrix. So M(i,j-1) should be the left element, M(i,j+1) the right element, M(i+1,j) the upper element, M(i-1) the lower element, M(i-1,j-1) the south-west element and so on.

The result should be saved in a new matrix. Should I use a nested for-loop to achieve this?

Answer 1

[John D'Errico]

This is actually trivially easy to solve, using a nice trick that is well worth remembering. (Note that this is a great trick that has worked for me in multiple places. For example, you can use it to write a simple code for the game of life.) I'll use my own matrix for the example, so we can see how it works.

A = magic(7)

A = 7×7

30 39 48 1 10 19 28 38 47 7 9 18 27 29 46 6 8 17 26 35 37 5 14 16 25 34 36 45 13 15 24 33 42 44 4 21 23 32 41 43 3 12 22 31 40 49 2 11 20

Now we wish to compute the mean of 3x3 blocks of the matrix. But near the edges, we have fewer neighbors. In a corner, even fewer neighbors yet. So what can we do? First, we will use convolution, but we do so in a tricky way. I'll do so here for the 2-dimensional problem, but higher dimensional problems would work just as nicely, because we can use convn.

K = ones(3,3); % the convolution kernel, indicating a 3x3 moving mean
Amean = conv2(A,K,'same')./conv2(ones(size(A)),K,'same')

Amean = 7×7

38.5000 34.8333 25.1667 15.5000 14.0000 21.8333 25.7500 34.3333 29.8889 20.2222 16.0000 18.0000 25.4444 29.1667 26.0000 20.7778 16.5556 17.7778 25.2222 31.8889 34.8333 16.5000 16.3333 17.5556 25.0000 32.4444 33.6667 33.5000 15.1667 18.1111 24.7778 32.2222 33.4444 29.2222 24.0000 20.8333 24.5556 32.0000 34.0000 29.7778 20.1111 15.6667 24.2500 28.1667 36.0000 34.5000 24.8333 15.1667 11.5000

Why did this work? Perhaps we need to look at the pieces. What did this do?

conv2(A,K,'same')

ans = 7×7

154 209 151 93 84 131 103 206 269 182 144 162 229 175 156 187 149 160 227 287 209 99 147 158 225 292 303 201 91 163 223 290 301 263 144 125 221 288 306 268 181 94 97 169 216 207 149 91 46

So that just adds up all of the elements in the vicinity of any original element in A.

conv2(ones(size(A)),K,'same')

ans = 7×7

4 6 6 6 6 6 4 6 9 9 9 9 9 6 6 9 9 9 9 9 6 6 9 9 9 9 9 6 6 9 9 9 9 9 6 6 9 9 9 9 9 6 4 6 6 6 6 6 4

Do you see how this works? It counted the number of terms in the sum. Essentially, it counts up the number of elements that were added to form the sum in the first call to conv2. And what is a mean, except the sum of a list of elements, divided by the number of elements in the list?

Now do you see how to make this work for a higher dimensional matrix? Easy, peasy.

A = randi(9,[4,4,4])

A = 
A(:,:,1) =

     5     2     2     9
     7     9     4     3
     8     7     2     2
     3     3     1     6


A(:,:,2) =

     9     7     3     9
     5     2     9     3
     3     6     7     9
     3     9     1     4


A(:,:,3) =

     3     8     8     9
     6     1     9     1
     1     8     6     4
     6     1     3     2


A(:,:,4) =

     8     9     2     5
     9     7     6     3
     8     9     2     6
     1     9     7     6

K = ones(3,3,3); % the convolution kernel, indicating a 3x3x3 moving mean
Amean = convn(A,K,'same')./convn(ones(size(A)),K,'same')

Amean = 
Amean(:,:,1) =

    5.7500    5.3333    5.1667    5.2500
    5.8333    5.3889    5.2778    5.1667
    5.4167    4.9444    4.8333    4.2500
    5.2500    4.4167    4.7500    4.0000


Amean(:,:,2) =

    5.3333    5.5000    5.4444    5.7500
    5.3889    5.4444    5.5185    5.5000
    4.8889    4.8148    4.5185    4.2222
    4.8333    4.3333    4.5000    3.9167


Amean(:,:,3) =

    6.1667    6.1667    5.6111    5.5833
    6.0556    5.9630    5.8519    5.6111
    5.2222    5.3333    5.1852    4.8889
    5.3333    5.0000    5.5000    4.7500


Amean(:,:,4) =

    6.3750    6.3333    5.6667    5.3750
    6.4167    6.1111    5.7222    5.0833
    5.5000    5.5000    5.0000    4.5833
    5.3750    5.0833    5.2500    4.5000

Answer 2

[Jan]

What about:

M = [3,4,8,6; 1,6,9,0; 2,5,7,1];
Result = movmean(movmean(M, 3, 1), 3, 2);

Or:

n = 3;  % Window size
[s1, s2]  = size(M);
D         = repmat(9, size(M));
D(1, 1)   = 4;
D(1, s2)  = 4;
D(s1, 1)  = 4;
D(s1, s2) = 4;
D(1,  2:s2-1) = 6;
D(s1, 2:s2-1) = 6;
D(2:s1-1, 1)  = 6;
D(2:s1-1, s2) = 6;

Result = conv2(M, ones(n), 'same') ./ D;

Or nicer:

Result = conv2(M, ones(n), 'same') ./ conv2(ones(size(M)), ones(n), 'same')

The approach with nested loops is ugly and slower.

[Radde1683]

Thank you really much for your answer! The first approach works fine but the second doesn't because n is not defined. But unfortunately I also have to do it with a for loop as a task.

[Jan]

n is the width of the window, in your case: n=3

You really need loops? Then this sounds like a homework question. It would not be useful, if I solve your homework. So post, what you have tried so far and ask a specific question, if you have one.

Start with 2 loops over the rows and columns of the input matrix. Get the submatrices using min() and max(), to crop the elements at the borders:

[s1, s2] = size(M)
...
    SubM = M(max(1, i1-w):min(s1, i1+w), max(1, i2-w):min(s2, i2+w));

where i1 and i2 are the counter over the elements of M. Define w such, that 2*w+1 is your window size.

If you have such a submatrix, calculate the mean. Note that sum(X(:))/numel(X) is much faster than mean().

[Márton Tamás Birosz]

Hey, is it possible to use 'movmean' the same way for a 3D (nxnxn) matrix?

[John D'Errico]

@Márton Tamás Birosz : See my answer.