2D sliding/moving/running average window

Hello, I want to perform a matrix 'patch' moving average but I am not sure how to. So, for example, defining my moving window as 3 x 3. I would like, for each element in the matrix, to perform the average of the elements as in the sequence below (excerpt) for an exemplary 6 x 10 matrix. I know this can be done more efficiently through the use of conv2 or in frequency domain with fft2, using the window kernel, but I still want to do it with a for loop, as this will be later translated to another real time language code. In addition, I would like to be able to select the elements as a 1D array (a reshape of the n_rows, n_cols matrix). At the edges of the matrix I select the next row, so it probably makes more sense to do it as a 1D array. Finally, when reaching the end of the matrix I can only select a few elements in the corner, but in this case I will do the average just of those elements.

thanks!!

[Matt J]

I know this can be done more efficiently through the use of conv2 or in frequency domain with fft2, using the window kernel, but I still want to do it with a for loop

It sounds like you already know how you're going to do it, so there doesn't seem to be a question here.

[Albert Zurita]

Yes, but I don't know how to select the elements, and especially given the 'corners' effect problem...

Answer 1

[Bruno Luong]

Just some hint if your want to implement efficiently using loop :

1. To do mean, you might do sliding sum onf the data, then divide by sliding sum on the array of size data but values are replaced with 1s (ones(size(data)).
2. The sliding 2D sum is "separable" meaning you just need to do sliding sum along one dimension, and apply the sliding sum along other dimension.
3. The sliding sum in 1D can be donne efficiently by simply adding sum of the previous step with the new entry element and substract the one that exits the window. If ther data is "quantified" such as multiples of the 2^power something, this method does not have cumulative error, otherwise you might have a small cumulative error compared to standard sum on the whole windows.

[Albert Zurita]

Thanks for the tips! I'm still struggling with how to select the elements. The reply from David is fine, but not directly used for a different window size.

[Bruno Luong]

I don't see yet you showing any effort to make your code/modification after many guidances, wonder how you translate the code to "realtime" language. But here we go :

% Create test data
A = peaks;
A = A + 0.1*randn(size(A));
win = [7 7];
B = slidingmean(A, win);

subplot(1,2,1)
surf(A)
title('Original')
subplot(1,2,2)
surf(B)
title('Sliding mean')

%%
function B = slidingmean(A, win)
B = slidingsum2(A, win) ./ slidingsum2(ones(size(A)), win);
end

%%
function B = slidingsum2(A, win)
if isscalar(win)
    win = [win,win];
end
B = slidingsum(A, 1, win(1));
B = slidingsum(B, 2, win(2));
end

%%
function B = slidingsum(A, dim, win)
m = ndims(A);
szA = size(A);
nA = szA(dim);
wh = floor(win/2);
wt = win-wh-1;
B = zeros(size(A));
idx = repmat({':'},1,m);
szS = szA;
szS(dim) = 1;
S = zeros(szS);
for k=1-wh:nA+wt
    in = k+wh;
    if in <= nA 
        idx{dim} = in;
        S = S + A(idx{:});
    end
    out = k-wt-1;
    if out >= 1
        idx{dim} = out;
        S = S - A(idx{:});
    end
    if k>=1 && k<=nA
        idx{dim} = k;
        B(idx{:}) = S;
    end
end
end

[Albert Zurita]

Thanks a lot! I have been trying different options offline, adding NaN rows and columns and other strategies without much success, didn't want to post here draft code only, apologies.

Your code looks good to me but it doesn't seem to handle well the edges exactly as in the sequence. When an edge is found it simply averages the values it can. However it does not jump to the next line as in the sequence. I think we need to add the modulus approach as in David's reply.

Thanks again

[Bruno Luong]

I have no idea what is "edges", "sequence" and "jump". Somthing probably related to your image in the topic that I admid don't take a close look.

My code simply does a plain sliding averaging on the 2D array, period.

[Albert Zurita]

Hi Bruno, and I truly apreciate your help. Your code is very good.

My point is on this case below. You average 29+30+39+40+49+50, which makes total sense for a 2D averaging. But I wanted to include as well in the average the continuation of the stream, so adding as well 31+41+51. I know it sounds a bit weird, but there is some sort of correlation between the data.

[Bruno Luong]

I see now that you want to do the wrap around with increment+1 in the row.

Instead of doing some kind of weird algorithm you could simply apply the normal algorithm on pad array on the right side by the left side:

All the specific wrap around is done in the way you pad your data, and leave the algorithm untouched.

A=zeros(4,6);
A(1:numel(A))=1:numel(A);
A

A = 4×6

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

win=3;
Apad=[A, [A(2:end,1:win-1); nan(1,win-1)]],

Apad = 4×8

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

% Doing average on Apad then remove the exceed
Bpad = slidingmean(Apad, win);
B = Bpad(:,1:size(A,2)),

B = 4×6

3.5000 5.5000 9.5000 13.5000 17.5000 13.8333 4.0000 6.0000 10.0000 14.0000 18.0000 14.3333 5.0000 7.0000 11.0000 15.0000 19.0000 16.6250 5.5000 7.5000 11.5000 15.5000 19.5000 18.0000

Note that if you want to discard NaN in the average you need to change the normalization from

... ./ slidingsum2(ones(size(A)), win)

to

... ./ slidingsum2(isfinite(A), win)

%%

function B = slidingmean(A, win)
valid = ~isnan(A);
A(~valid) = 0;
B = slidingsum2(A, win) ./ slidingsum2(valid, win);
end
%%
function B = slidingsum2(A, win)
if isscalar(win)
    win = [win,win];
end
B = slidingsum(A, 1, win(1));
B = slidingsum(B, 2, win(2));
end
%%
function B = slidingsum(A, dim, win)
m = ndims(A);
szA = size(A);
nA = szA(dim);
wh = floor(win/2);
wt = win-wh-1;
B = zeros(size(A));
idx = repmat({':'},1,m);
szS = szA;
szS(dim) = 1;
S = zeros(szS);
for k=1-wh:nA+wt
    in = k+wh;
    if in <= nA 
        idx{dim} = in;
        S = S + A(idx{:});
    end
    out = k-wt-1;
    if out >= 1
        idx{dim} = out;
        S = S - A(idx{:});
    end
    if k>=1 && k<=nA
        idx{dim} = k;
        B(idx{:}) = S;
    end
end
end

[Albert Zurita]

Great, now that case is handled perfectly. I made the modifications. There are still however cases not fully covered. For instance when we need to go 'one row up'. For example those two cases should give 41 and 48, respectively, but instead I get 41.5 and 46.5 when using your code.

[Bruno Luong]

Cheese, I hope you should know by now what you want with this weird wrapping and make the Padding exactly like you want it to be. This is the LAST time I adapt the code for you.

I warn you, next question along special case on top level, this answer will be deleted, because I don't have any patient to follow you, sorry.

A=(1:10)+(0:10:50)'

A = 6×10

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 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

win = 3;

Apad = [[nan(1,win-1); A(1:end,end-win+2:end)], [A; nan(1,size(A,2))], [A(2:end,1:win-1); nan(2,win-1)]]

Apad = 7×14

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

Bpad = slidingmean(Apad, win);
B = Bpad(1:end-1,win:size(A,2)+win-1)

B = 6×10

7.2000 7.0000 8.0000 9.0000 10.0000 11.0000 12.0000 13.0000 14.0000 15.0000 12.3750 12.0000 13.0000 14.0000 15.0000 16.0000 17.0000 18.0000 19.0000 20.0000 21.0000 22.0000 23.0000 24.0000 25.0000 26.0000 27.0000 28.0000 29.0000 30.0000 31.0000 32.0000 33.0000 34.0000 35.0000 36.0000 37.0000 38.0000 39.0000 40.0000 41.0000 42.0000 43.0000 44.0000 45.0000 46.0000 47.0000 48.0000 49.0000 48.6250 48.0000 47.0000 48.0000 49.0000 50.0000 51.0000 52.0000 53.0000 54.0000 53.8000

function B = slidingmean(A, win)
valid = ~isnan(A);
A(~valid) = 0;
B = slidingsum2(A, win) ./ slidingsum2(valid, win);
end
%%
function B = slidingsum2(A, win)
if isscalar(win)
    win = [win,win];
end
B = slidingsum(A, 1, win(1));
B = slidingsum(B, 2, win(2));
end
%%
function B = slidingsum(A, dim, win)
m = ndims(A);
szA = size(A);
nA = szA(dim);
wh = floor(win/2);
wt = win-wh-1;
B = zeros(size(A));
idx = repmat({':'},1,m);
szS = szA;
szS(dim) = 1;
S = zeros(szS);
for k=1-wh:nA+wt
    in = k+wh;
    if in <= nA 
        idx{dim} = in;
        S = S + A(idx{:});
    end
    out = k-wt-1;
    if out >= 1
        idx{dim} = out;
        S = S - A(idx{:});
    end
    if k>=1 && k<=nA
        idx{dim} = k;
        B(idx{:}) = S;
    end
end
end

[Albert Zurita]

Excellent Bruno. That works perfectly. Note that this behaviour was already anticipated in my original post. That's why I find it difficult to code. And honestly I still have some troubles in fully following your answer. If you get annoyed with my enquiries I apologize for that, but please don't reply and forget about the post in such case. For relatively newbies like me it's not always a clear path. In any case, thank you for your time and efforts.

[Bruno Luong]

Sorry I was a little bit nerveous yesterday.

[Albert Zurita]

No problem. I am now parametrizing for a non-square window. Playing with the padding step.

[Albert Zurita]

...and still not clear how to deal with a case where win is not a single value (describing a square matrix) but defined by two values:

Instead of:

win = 3;

I would use:

win = [3 5];

I presume the deal is only in padding and unpadding operations for Apad and Bpad

[Bruno Luong]

Just folow the simple rule:

• Working along the first dimension (rows) => use win(1)
• Working along the second dimension (columns) => use win(2)

[Albert Zurita]

Indeed, yet the padding was not that evident though ;)

[Bruno Luong]

Well you start with my original code, then make these modifications:

• when you see win in the first array index, replace it with win(1),
• when you see win in the second array index, replace it with win(2).

[Albert Zurita]

For:

win = [nrows ncols];

I had to do these modifications, which are slightly different:

Apad = [[nan(1,win(2)-1); A(1:end,end-win(2)+2:end)], [A; nan(1,size(A,2))], [A(2:end,1:win(1)-1); nan(2,win(1)-1)]]

And for Bpad I had to use this:

B = Bpad(1:end-1,win(1):size(A,2)+win(1)-1);

This seems to be the only way to make it work...

[Bruno Luong]

Obviously you don't follow the rule I state.

[Albert Zurita]

I known. If I understood correctly you proposed:

Apad = [[nan(1,win(1)-1); A(1:end,end-win(1)+2:end)], [A; nan(1,size(A,2))], [A(2:end,1:win(2)-1); nan(2,win(2)-1)]];
B = Bpad(1:end-1,win(2):size(A,2)+win(2)-1);

But then the Apad matrix does not seem to pad the number of columns correctly

[Bruno Luong]

"Apad matrix does not seem to pad the number of columns correctly."

Sorry but You still get it wrong again. The correct is

Apad = [[nan(1,win(2)-1); ...
         A(1:end,end-win(2)+2:end)], ...
         [A; nan(1,size(A,2))], ...
       [A(2:end,1:win(2)-1); ...
        nan(2,win(2)-1)]];
Bpad = slidingmean(Apad, win);
B = Bpad(1:end-1,win(2):size(A,2)+win(2)-1);

Up to know you only specify the wrap around on the horizontal direction, so only win(2) matter.

You might expect something different now, and I again insist that you must pad according whatever the weird wrap around you want. It's up to you to change the pad if you change the wrap shiftting that beside you nobody know.

[Bruno Luong]

Also you might change some of the hard-code first index

nan(1, ...
A(2:...
Bpad(1:end-1 ...

with expression using win(1)

Answer 2

[David Hill]

a=randi(100,15);%whatever initial size of your matrix
b=size(a,1);
m=zeros(b);
for x=1:b^2
  idx=[x-b-1,x-1,x+b-1;x-b,x,x+b;x-b+1,x+1,x+b+1];
  if idx(2,3)>b^2
      idx(:,3)=[];
  end
  if mod(idx(3,2),b)==1
      idx(3,:)=[];
  end
  if idx(2,1)<1
      idx(:,1)=[];
  end
  if mod(idx(1,2),b)==0
      idx(1,:)=[];
  end
  m(x)=mean(a(idx),'all');
end

[Albert Zurita]

That looks good! But how can I parametrize the size of the running window to another value instead of 3 x 3 (e.g. 3 x 5)?

[Image Analyst]

Albert, did you even see my answer below? Please read it. It does what you asked. You can specify the image file and window size and is 100% manual (non-vectorized). At least give it a try even if you prefer Bruno's answer (which has vectorized indexing).

[Albert Zurita]

Yes, please see my comment to your reply below

Answer 3

[Image Analyst]

See my attached "manual" convolution demo. It lets you pick a standard demo image (converts to gray scale if necessary) and then asks you for the number of rows and columns in the scanning filter window (kernel). Then it has a 4-nested for loop where it sums the image pixel times the kernel value. It also counts the number of pixels in the kernel that overlap the image so in essence it shrinks the window when it "leaves the image" by ignoring those pixels. So if the 3x3 kernel is centered over (1,1) only 4 pixels are considered rather than the full 9 because 5 pixels are off the image and don't overlap. The input image and output image are displayed side by side. It does not use any MATLAB convolution functions or vectorized indexing to get any of the pixels in the window so it's completely 100% manual.

[Albert Zurita]

Hi, indeed your proposal looks very good. The only issue is the 'circular smooth' that Bruno has implemented as per my original post. Since there is some sort of similarity between the end of one row and the start of the next one Bruno creates that padded matrix. Perhaps your code could be adapted to work with the padded matrix, but the indexing should be changed.