2D gaussian filter with a variable sigma

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?

Answer 1

[Wolfgang Schwanghart]

Hi Chad,

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)
    % image is zero padded 
    % 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

[Chad Greene]

Oh, interesting, I was not aware of nlfilter. However for a gaussian filter, the size of the window should vary depending on the sigma of the filter. Is there any way to do this without keeping m and n fixed?

[Wolfgang Schwanghart]

True, the size of the window is fixed, but the standard deviation of the Gaussian Kernel varies. I have now updated the code. While m and n remain fixed (required by nlfilter), the standard deviation of the kernel varies freely, yet it never extends the size of the kernel.

[Image Analyst]

I second this answer - I was going to suggest nlfilter() also, until I saw Wolfgang beat me to it.