Gaussian filter with fspecial versus imgaussfilt

I have two questions. First, I want to filter a function with a Gaussian that has a standard deviation with dimensions of length; units are microns. The dx spacing in the code below has units of microns. If the filters in Matlab require their standard deviation provided in 'nodes', is the way to associate the two as such:

x = linspace(-5,5,200); % Units in microns
dx = x(end) - x(end-1); % Spacing in microns
stdev = 0.5; % Desired Gaussian filter standard deviation. Units in microns.
nodesPerSigma = round(stdev/dx,0); % Is this the way to associate nodal sigma with a dimensional sigma?
% Filter with IMGAUSSFILT
zf = imgaussfilt(Z,nodesPerSigma);

My second question is: in the code below, why do zf1 and zf2 produce different results if zf2 matches the result from imgauss filt? The built-in function imgaussfilt selects the filter size based on the sigma provided. It uses 2*ceil(2*sigma)+1 to determine the filter size. This is smaller than the 'rule of thumb' (see link below) of 2*3*sigma. Does this mean that imgaussfilt is not using the entire Gaussian??

clear;
x = linspace(-5,5,200); % Units in microns
[X,Y] = meshgrid(x,x);

Z = X .* exp(-X.^2 - Y.^2);
%figure; surf(X,Y,Z)
%figure; plot(X(100,:),Z(100,:))

dx = x(end) - x(end-1); % Spacing in microns
stdev = 0.5; % Desired Gaussian filter standard deviation. Units in microns.

nodesPerSigma = round(stdev/dx,0); % Is this the way to associate nodal sigma with a dimensional sigma?

% Filter with IMGAUSSFILT
zf = imgaussfilt(Z,nodesPerSigma);

filtSize = 2*3*nodesPerSigma;
if mod(filtSize,2) == 0
    filtSize = filtSize + 1;
end
% Filter with fspecial's Gaussian using the 'rule of thumb' for the filter
% size. See link below.
gfilter = fspecial('gaussian',[filtSize filtSize],nodesPerSigma);
zf1 = imfilter(Z,gfilter,'same');

% Filter with fspecial's Gaussian but using the filter size in documentation.
filtSize2 = 2*ceil(2*nodesPerSigma)+1; % This will produce a result matching the imgaussfilt.
gfilter2 = fspecial('gaussian',[filtSize2 filtSize2],nodesPerSigma);
zf2 = imfilter(Z,gfilter2,'same');

figure; plot(X(100,:),Z(100,:),'k',X(100,:),zf(100,:),'r',X(100,:),zf1(100,:),'g',X(100,:),zf2(100,:),'*b')
legend('Original Function','Using imgaussfilt','Using fspecial Gaussian','Using fspecial Gaussian Again')

Link containing info of a rule of thumb for sizing the Gaussian filter. In essence the goal is to include all of the Gaussian curve and not much more, for efficiency.

https://www.mathworks.com/matlabcentral/answers/231351-what-is-the-correct-number-of-pixels-for-a-gaussian-with-a-given-standard-deviation?s_tid=srchtitle

Answer 1

[Maneet Kaur Bagga]

Hi Paul,

As the answer to your first question:

The formula you've used is:

nodesPerSigma = round(stdev/dx,0);

The "stdev" is the desired standard deviation of the Gaussian filter in "microns", and "dx" is te spacing between the adjacent points in the grid, which is also in "microns". The conversion is necessary because the filtering functions expect the standard deviation in terms of pixels (or nodes), and not in physical unit like microns.

The workaround for the second question is:

The difference seen id due to the filter size, "imgaussfilt" automatically calculates the filter size, but when using "fspecial" you have control over the filter size.

• When using "fspecial", explicitly calculate the filter size as done in your code. For applications requiring a closer match to "imgaussfilt", use the same formula for filter size calculation "(2*ceil(2*sigma)+1)".
• For more extensive smoothing, as attempted with "zf1", your approach of using a larger filter size based on the "rule of thumb" is valid.

Please refer to the code below for better understanding:

% Assuming stdev, dx, and Z are defined as per your code

% Correct conversion from microns to nodes
nodesPerSigma = round(stdev/dx,0); 

% Using imgaussfilt directly with nodesPerSigma
zf = imgaussfilt(Z,nodesPerSigma);

% For fspecial, matching imgaussfilt's behavior
filtSize2 = 2*ceil(2*nodesPerSigma)+1; % Ensures matching result with imgaussfilt
gfilter2 = fspecial('gaussian',[filtSize2 filtSize2],nodesPerSigma);
zf2 = imfilter(Z,gfilter2,'same');

% Plotting to compare
figure; plot(X(100,:),Z(100,:),'k',X(100,:),zf(100,:),'r',X(100,:),zf2(100,:),'*b')
legend('Original Function','Using imgaussfilt','Using fspecial Gaussian Again')

Below I have attached MathWorks Documentation for further understanding:

imgaussfilt: https://www.mathworks.com/help/images/ref/imgaussfilt.html

fspecial: https://www.mathworks.com/help/images/ref/fspecial.html

Hope this helps!