removing up and down spikes from a surface

I have a 2-d array of data that represents a surface. The surface was scanned with a laser, and the laser data has both upward and downward spikes. However, the surface also has legitimate topography (it is smooth, but with transitions). I would like to preserve the topography, but get rid of the spikes.

I have had success using a 7,7 median filter to remove the spikes. However, the median filter also removes a lot of the "valid" detail. So, I use the MAX function to compare the original data with the filtered data. That gives me a noise-free top surface, while preserving the data at the bottom. I also use the MIN function to get a noise-free bottom surface. So the max data now has downward spikes and the min data has upward spikes.

My challenge: Is there a way to combine these two surfaces, possibly including the original data, to get rid of both sets of spikes, while preserving the underlying detail? I think I almost want the opposite of a median of the three arrays, where if two of the three corresponding elements are the same, I want to keep the third one.

Answer 1

[DGM]

It would help to have an example of the data to know what the ranges look like.

At first glance, I might try to throw an adaptive median noise removal filter at it (or maybe a fixed median NR filter with appropriately tailored thresholds)

https://www.mathworks.com/matlabcentral/answers/22368-adaptive-median-filter#answer_1143972

% some fake data (a smooth surface)
sz = 50;
Z = peaks(sz);

% add some impulse noise
noise = 2*randn(sz); imrange(noise)
noise(abs(noise)<5) = 0;
Z = Z + noise;

surf(Z)

% adaptive median noise removal filter
% these parameters specify an 11x11 maximum window size
% and a more generous relative order statistic for identifying noise
Zf = amedfilt(Z,5,0.5); % MIMT
surf(Zf)

That would only remove local outliers, but using a plain median filter without any noise discrimination (e.g. medfilt2()) will remove more than the spikes, and may change the apparent shape of the surface.

% a plain spatial filter with no discrimination
Zf = medfilt2(Z,[7 7]);
surf(Zf)

amedfilt() is part of MIMT. It's not exactly meant for this purpose, but depending on the relative scale of noise vs non-noise features, it might work adequately.

EDIT:

There's also rmoutliers(), though I admit that I'm not as familiar with it:

https://www.mathworks.com/help/matlab/ref/rmoutliers.html

[Jim McIntyre]

Thank you for this answer. It provided some good ideas, and I tried the AMEDFILT option. However, the noise in my data appears to be more diffuse than salt-and-pepper noise, so it is difficult to separate the noise from real features. AMEDFILT, even trying several options for radius, did not work nearly as well as MEDFILT2, followed by min(Raw, Filtered).

Here is an example of some raw data:

This is what Medfilt2 does:

Here is Min(Raw, Filt):

And here is Max(Raw, Filt):

Here are raw and max views from below the surface: