Asked by Chad Greene
on 1 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?

Answer by Wolfgang Schwanghart
on 3 Apr 2019

Edited by Wolfgang Schwanghart
on 3 Apr 2019

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
on 3 Apr 2019

Wolfgang Schwanghart
on 3 Apr 2019

Image Analyst
on 4 Apr 2019

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

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.