How to apply a moving median filter on a time series of 2D scans in Matlab?

Question

I have a huge set of data of a timelapse of 2D laser scans of waves running up and down stairs (see fig.1fig.2fig.3). There is a lot of noise in the scans, since the water splashes a lot. Now I want to smoothen the scans.

I have 2 questions:

1. How do I apply a moving median filter (as recommended by another study dealing with a similar problem)? I can only find instructions for single e.g. (x,y) or (t,y) plots but not for x and y values that vary over time. Maybe an average filter would do it as well, but I do not have a clue on that either.

2. The scanner is at a fixed point (222m) so all the data spikes point towards that point at the ceiling. Is it possible or necessary to include this into the smoothing process?

This is the part of the code (I hope it's enough to get it):

% Plot data as real time profile
x1=data.x;y1=data.y;
t=data.t;
% add moving median filter here?
h1=plot(x1(1,:),y1(1,:));
axis([210 235 3 9]) 
ht=title('Scanner data');
for i=1:1:length(t);    
set(h1,'XData',x1(i,:),'YData',y1(i,:));set(ht,'String',sprintf('t = %5.2f 
s',data.t(i)));pause(.01);end

The data.x values are stored in a (mxn) matrix in which the change in time is arranged vertically and the x values i.e. "laser points" of the scanner are horizontally arranged. The data.y is stored in the same way. The data.t values are stored in a (mx1) matrix.

I hope I explained everything clearly and that somebody can help me. I am already pretty desperate about it... If there is anything missing or confusing, please let me know.

Answer 1

If you're trying to apply a median filter in the x-y plane, then consider using medfilt2 from the Image Processing Toolbox. Note that this function only accepts 2-D inputs, so you'll have to loop over the third dimension.

Also note that medfilt2 assumes that the x and y data are uniformly spaced, so if your x and y data don't fall onto a uniformly spaced grid you may have to manually loop over indices, extract the corresponding patches, and compute the median.

If you can/want to apply an averaging filter instead of a median filter, and if you have uniformly spaced data, then you can use convn to compute a k x k moving average by doing:

y = convn(x, ones(k,k)/(k*k), 'same');

Note that you'll get some bias on the boundaries because you're technically trying to compute an average of k^2 pixels when you have less than that number of values available.

Alternatively, you can use nested calls to movmean since the averaging operation is separable:

y = movmean(movmean(x, k, 2), k, 1);

If your grid is separable, but not uniform, you can still use movmean, just use the SamplePoints name-value pair:

y = movmean(movmean(x, k, 2, 'SamplePoints', yv), k, 1, 'SamplePoints', xv);

You can also control the endpoint handling in movmean with the Endpoints name-value pair.