# Nonlinear operations using imfilter 9

Posted by **Steve Eddins**,

Some nonlinear image processing operations can be expressed in terms of linear filtering. When this is true, it often provides a recipe for a speedy MATLAB implementation. Today I'll show two examples: Computing the local standard deviation, and computing the local geometric mean.

The local standard deviation operator is often used as a measure of "busy-ness" throughout the image. For each pixel, the standard deviation of that pixel's neighbors is computed:

sqrt( sum ( (x - xbar).^2 ) )

(Scale factors omitted.)

Mathematically, the summation can be rewritten as:

sum(x.^2) - sum(x).^2/N

Both of these local sums can be computed by using `imfilter` with an all-ones filter.

```
I = im2double(imread('cameraman.tif'));
imshow(I)
```

*Original image credit: Massachusetts Institute of Technology*

To compute the sum(x.^2) term, we square (elementwise) the input to `imfilter`.

h = ones(5,5); term1 = imfilter(I.^2, h); imshow(term1, [])

Notice the dark band around the edge. This is because `imfilter` zero-pads by default. We might want to use the 'symmetric' option instead.

```
term1 = imfilter(I.^2, h, 'symmetric');
imshow(term1, [])
```

To compute the sum(x).^2 term, we square the output of `imfilter`.

```
term2 = imfilter(I, h, 'symmetric').^2 / numel(h);
imshow(term2, [])
```

Then we subtract the second term from the first and take the square root.

```
local_std = sqrt(term1 - term2); % scale factor omitted
imshow(local_std, [])
```

*Cautionary notes*

- The procedure shown here is not always a numerically sound way to compute the standard deviation, because it can suffer from both overflow (squared terms) and underflow (cancellation in the subtraction of large numbers). For typical image pixel values and window sizes, though, it works reasonably well.

- Round-off error in the computation of term1 and term2 can sometimes make (term1 - term2) go slightly negative, resulting in complex outputs from square root operator. Avoid this problem by using this code:

local_std = sqrt(max(term1 - term2, 0));

Recent releases of the Image Processing Toolbox include the function `stdfilt`, which does all this work for you.

The local geometric mean filter multiplies together all the pixel values in the neighborhood and then takes the N-th root, where N is the number of pixels in the neighborhood. The geometric mean filter is said to be slightly better than the arithmetic mean at preserving image detail.

Use the old logarithm trick to express the geometric mean in terms of a summation. Then `imfilter` can be used to compute the neighborhood summation, like this.

```
geo_mean = imfilter(log(I), h, 'replicate');
geo_mean = exp(geo_mean);
geo_mean = geo_mean .^ (1/numel(h));
imshow(geo_mean, [])
```

Does anyone have other creative uses of `imfilter` to share?

Copyright 2008 The MathWorks, Inc.

Get
the MATLAB code

Published with MATLAB® 7.6

## 9 CommentsOldest to Newest

This very same “sum of squares minus sqaure of sum” consideration drove me to re-implement the Kuwahara nonlinear filter in a faster way, though I did not use imfilter but just straightforward convolution, as I did not want the Image Processing Toolbox to be needed:

http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=15027&objectType=FILE

Nice example by the way Steve! And thanks for pointing to the stdfilt function.

Luca—Thanks for submitting your code to the File Exchange!

Hi Steve,

You mention that the local standard deviation technique measures the ‘business’ of an image. How is this different, if at all, from edge detection? It looks like your first example does a really amazing job at finding edges, albeit they are rather thick.

thanks

-francisco

Francisco—At last count, there were 13,273,474 journal and conference papers describing different methods of edge detection. :-)

Think about regions of differing texture. The local standard deviation operator would have a relatively higher response in regions of “busier” texture, and it would have a relatively lower response in “smoother” regions. This is different than what most people think an edge detector does.

I am aware this post is old, but just to note:

At this line:

geo_mean = imfilter(log(I), h, ‘replicate’);

You are calculating the log() of an array of doubles where 0 <= x <= 1.

That will produce an array of negative values, and possibly -Inf.

also replicating to the image as black squares.

Diego—It doesn’t matter; the procedure still gives the correct answer. Note that the geometric mean of a set of values containing a 0 is 0. Try this:

Hi Steve,

I realize this is old but I have a question about it. Why is there such a significant difference between your imfilter method and stdfilt with numbers in the range 0-255. Is it just numerical round off? The thing that irks me is that the greatest difference between the two is in the middle where there is a significant standard deviation and that the standard deviation from your imfilter example doesn’t seem reasonable.

Thanks!

-Sean

%%%

I = double(imread(‘cameraman.tif’));

h = ones(11,11);

Istd1 = sqrt(imfilter(I.^2, h, ‘symmetric’) – imfilter(I, h, ‘symmetric’).^2 / numel(h));

Istd2 = stdfilt(I,h);

imshow(Istd1-Istd2,[])

norm(Istd1(:)-Istd2(:))

%{

ans = 75950

%}

max([Istd1(:), Istd2(:)])

%{

ans = 1092.4 99.719

%}

Okay, I guess by opening stdfilt you can see the difference between the two. So why the discrepancy?

Sean—A lot of the mathematical details are omitted in this post. See the places, for example, where I include a note that says “scale factors omitted.” There is a significant difference between stdfilt and the code in this post because, well, they don’t compute the same thing.