Removing Gaussian noise from an image (method comparison)

Can someone give me a steer on the best way to remove random speckle from a monochrome image?

I have a series of 16-bit monochrome images which show the pattern of light uniformity on a flat surface at different wavelengths. The images are stored in a 3-d array such that the first two dimensions (row, column) are the image pixel values and the third dimension ('page') is the wavelength value.

The images are 'noisy'. That is, if you chose a page of the array and plotted its 2-d surface, it would not be smooth but show small peaks and troughs either side of a mean value at every pixel location. I want to smooth out these peaks and troughs.

The imgaussfilt function applies a Gaussian smoothing filter to images. There is also a method using the polyfitn function on file exchange. What is the difference in outcome between these approaches?

I tried comparing the two using the code below (my histogram comparison method is crude!). Visually, method 2 looks best but I don't know whether it's a fair test. Also, I can select the polyfitn polynomial order that works best only when I know the original noise-less signal. How should I select the order when the signal is unknown?

%Make up 3 dummy test images 64 x 80
rows = 64;
columns = 80;
y=1:rows;
x=1:columns;
x1=(2-exp(x./1000)).*50000;
y1= 1 - (y.*0.008-0.4).^2; % profile at wavelength 1
y2= 1 - (y.*0.009-0.4).^2; % profile at wavelength 2
y3= 1 - (y.*0.01-0.4).^2; % profile at wavelength 3
y1=y1';y2=y2';y3=y3';
zSignal= cat (3, y1*x1, y2*x1, y3*x1);

% Add some noise to images
z= zSignal + 350.*rand(rows,columns,3); % randn more realistic?
z= round(z); %

% SMOOTHING METHOD 1 (Try blurring the image on page 1)
image1 = uint16(squeeze(z(:,:,1))); % create the dummy image
zBlur = double(imgaussfilt(image1,2));
figure;
histogram (zBlur-zSignal(:,:,1)); % check for absolute error #1

%SMOOTHING METHOD 2 (Create a model zPred for each page)
% (with thanks to ImageAnalyst)
zPred = zeros (64,80,3);
[X,Y] = meshgrid(1:columns, 1:rows);
x1d = reshape(X, numel(X), 1);
y1d = reshape(Y, numel(Y), 1);
zPred = zeros (64,80,3);
for page = 1:3
    image = uint16(squeeze(z(:,:,page))); % create the dummy image
    z1d = double(reshape(image, numel(image), 1));
    polynomialOrder = 4;
	p = polyfitn([x1d y1d], z1d, polynomialOrder);
    zg = polyvaln(p, [x1d y1d]);
	zPred(:,:,page) = reshape(zg, [rows columns]);
end
figure;
histogram (zPred(:,:,1)-zSignal(:,:,1)); % check for error #2
return;

[Mathieu NOE]

as a very general remark, if you only smooth data you take the easiest way but you may not get the best solution.

If you bring more info's and knowledge in the process and decide to fit instead an (appropriate) model , you should normally get a (much) better result. The difficulty here is to be sure to have a model that ensures a good fit to your data. if you are sure about this assumption, then I would personnaly stick to the "model fit" approach rather than just smoothing the data.

But there is always the risk that one day your model is making nonsense vs new data and you will have to find another model. that's life...better results comes at the price of greater attention to all those details .

a simple remark about your code above , another basic suggestion : I can also get an idea about how good is the result by taking the std of the error surface (instead of plotting the histogram) , but both approaches give here the same answer : model fit is better by far vs smoothing :

std1 = 35.6091 % smoothing

std2 = 4.9295 % polyfitn

%Make up 3 dummy test images 64 x 80
rows = 64;
columns = 80;
y=1:rows;
x=1:columns;
x1=(2-exp(x./1000)).*50000;
y1= 1 - (y.*0.008-0.4).^2; % profile at wavelength 1
y2= 1 - (y.*0.009-0.4).^2; % profile at wavelength 2
y3= 1 - (y.*0.01-0.4).^2; % profile at wavelength 3
y1=y1';y2=y2';y3=y3';
zSignal= cat (3, y1*x1, y2*x1, y3*x1);
% Add some noise to images
z= zSignal + 350.*rand(rows,columns,3); % randn more realistic?
z= round(z); %
% SMOOTHING METHOD 1 (Try blurring the image on page 1)
image1 = uint16(squeeze(z(:,:,1))); % create the dummy image
zBlur = double(imgaussfilt(image1,2));
figure;
histogram (zBlur-zSignal(:,:,1)); % check for absolute error #1

std1 = std(zBlur-zSignal(:,:,1),1,'all') % check for absolute error #1

%SMOOTHING METHOD 2 (Create a model zPred for each page)
% (with thanks to ImageAnalyst)
zPred = zeros (64,80,3);
[X,Y] = meshgrid(1:columns, 1:rows);
x1d = reshape(X, numel(X), 1);
y1d = reshape(Y, numel(Y), 1);
zPred = zeros (64,80,3);
for page = 1:3
    image = uint16(squeeze(z(:,:,page))); % create the dummy image
    z1d = double(reshape(image, numel(image), 1));
    polynomialOrder = 4;
	p = polyfitn([x1d y1d], z1d, polynomialOrder);
    zg = polyvaln(p, [x1d y1d]);
	zPred(:,:,page) = reshape(zg, [rows columns]);
end
figure;
histogram (zPred(:,:,1)-zSignal(:,:,1)); % check for error #2

std2 = std(zPred(:,:,1)-zSignal(:,:,1),1,'all') % check for absolute error #2

[Steve Francis]

Thank you very much, @Mathieu NOE, for your comment and excellent suggestion regarding standard deviation as a comparator. You're right; empirically, I have found that model fitting gives much better results, particularly if I can use a lower polynomial order. 3 or 4 seems to work well for my images. I'm deducing this using a 'before' and 'after' visual comparison so it's somewhat subjective.