# 2D gaussian filter with a variable sigma

112 views (last 30 days)
Chad Greene on 1 Apr 2019
Commented: Image Analyst on 4 Apr 2019
I have a large gridded dataset I'd like to lowpass filter. The catch is, need to specify a different sigma value for each pixel of the grid.
Using a single value of sigma for the entire grid is easy with imgaussfilt, so for example, I can filter the grid I using a sigma value of 3 like this:
% A grid of data:
I = 10*rand(100);
% The grid, filtered to sigma=3:
If = imgaussfilt(I,3);
But I don't want tuse a single value of sigma for the entire grid. Instead, I want to specify a different sigma for every pixel. For this 100x100 grid it's easy to loop through each row and column, filtering the grid to specified values of sigma like this:
% A grid of sigma values corresponding to the grid of data:
sigma = abs(peaks(100))+0.1;
% Preallocate:
If2 = NaN(size(I));
% Loop through each row and column:
for row = 1:size(I,1)
for col = 1:size(I,2)
% Perform Gaussian filter on the entire grid using sigma(row,col):
tmp = imgaussfilt(I,sigma(row,col));
% Only save the value corresponding to this row and col:
If2(row,col) = tmp(row,col);
end
end
The nested loop approach gives the answer I want, but it's inelegant, and slow for very large grids. Is there a better way to define a locally-varying filter?

Wolfgang Schwanghart on 3 Apr 2019
Edited: Wolfgang Schwanghart on 3 Apr 2019
below is some code that does the trick using nlfilter. To be more efficient, I created a look-up table with a finite set of sliding windows.
% Variable standard deviations
st = 1:100;
% Size of sliding window
% Should be an odd number
m = 101;
kcenter = ceil(m/2);
% Cell array of filters
st = cellfun(@(x) fspecial('gaussian',[m m],x),num2cell(st),'Unif',false);
% Image to be smoothed
A = rand(200);
% Image with variable standard deviations
STSIZE = randi(20,size(A))+5;
% Image with linear indices
IX = reshape(1:numel(A),size(A));
B = nlfilter(IX,[m m],...
@(ix) fun(ix,A,STSIZE,st,kcenter));
imagesc(B)
function out = fun(ix,A,STSIZE,st,kcenter)
if any(ix(:) == 0)
% return nan
out = nan;
else
% get window
ixc = ix(kcenter,kcenter);
stsize = STSIZE(ixc);
f = st{stsize};
a = A(ix);
out = f(:)'*a(:);
end
end