Designing and Implementing Linear Filters in the Spatial Domain

Overview

Filtering is a technique for modifying or enhancing an image. For example, you can filter an image to emphasize certain features or remove other features. Image processing operations implemented with filtering include smoothing, sharpening, and edge enhancement. Filtering is a neighborhood operation, in which the value of any given pixel in the output image is determined by applying some algorithm to the values of the pixels in the neighborhood of the corresponding input pixel. A pixel's neighborhood is some set of pixels, defined by their locations relative to that pixel.  Linear filtering is filtering in which the value of an output pixel is a linear combination of the values of the pixels in the input pixel's neighborhood.

Convolution

Linear filtering of an image is accomplished through an operation called convolution. Convolution is a neighborhood operation in which each output pixel is the weighted sum of neighboring input pixels. The matrix of weights is called the convolution kernel, also known as the filter. A convolution kernel is a correlation kernel that has been rotated 180 degrees.

For example, suppose the image is

A = [17  24   1   8  15

     23   5   7  14  16

      4   6  13  20  22

     10  12  19  21   3

     11  18  25   2   9]

and the convolution kernel is

h = [8   1   6

     3   5   7

     4   9   2]

The following figure shows how to compute the (2,4) output pixel using these steps:

1. Rotate the convolution kernel 180 degrees about its center element.
2. Slide the center element of the convolution kernel so that it lies on top of the (2,4) element of A.
3. Multiply each weight in the rotated convolution kernel by the pixel of A underneath.
4. Sum the individual products from step 3.

Performing Linear Filtering of Images Using imfilter

Filtering of images, either by correlation or convolution, can be performed using the toolbox function imfilter. This example filters an image with a 5-by-5 filter containing equal weights. Such a filter is often called an averaging filter.

I = imread('coins.png');

h = ones(5,5) / 25;

I2 = imfilter(I,h);

imshow(I), title('Original Image');

figure, imshow(I2), title('Filtered Image')

Data Types

The imfilter function handles data types similarly to the way the image arithmetic functions do The output image has the same data type, or numeric class, as the input image. The imfilter function computes the value of each output pixel using double-precision, floating-point arithmetic. If the result exceeds the range of the data type, the imfilter function truncates the result to that data type's allowed range. If it is an integer data type, imfilter rounds fractional values.

Because of the truncation behavior, you might sometimes want to consider converting your image to a different data type before calling imfilter. In this example, the output of imfilter has negative values when the input is of class double.

A = magic(5)

A =

    17    24     1     8    15

    23     5     7    14    16

     4     6    13    20    22

    10    12    19    21     3

    11    18    25     2     9

h = [-1 0 1]

h =

    -1     0     1

imfilter(A,h)

ans =

    24   -16   -16    14    -8

     5   -16     9     9   -14

     6     9    14     9   -20

    12     9     9   -16   -21

    18    14   -16   -16    -2

Notice that the result has negative values. Now suppose A is of class uint8, instead of double.

A = uint8(magic(5));

imfilter(A,h)

ans =

    24     0     0    14     0

     5     0     9     9     0

     6     9    14     9     0

    12     9     9     0     0

    18    14     0     0     0

Since the input to imfilter is of class uint8, the output also is of class uint8, and so the negative values are truncated to 0. In such cases, it might be appropriate to convert the image to another type, such as a signed integer type, single, or double, before calling imfilter.

However, if you want to perform filtering using convolution instead, you can pass the string 'conv' as an optional input argument to imfilter. For example:

A = magic(5);

h = [-1 0 1]

imfilter(A,h)   % filter using correlation

ans =

    24   -16   -16    14    -8

     5   -16     9     9   -14

     6     9    14     9   -20

    12     9     9   -16   -21

    18    14   -16   -16    -2

imfilter(A,h,'conv')   % filter using convolution

ans =

   -24    16    16   -14     8

    -5    16    -9    -9    14

    -6    -9   -14    -9    20

   -12    -9    -9    16    21

   -18   -14    16    16     2

Boundary Padding Options

When computing an output pixel at the boundary of an image, a portion of the convolution or correlation kernel is usually off the edge of the image

When the Values of the Kernel Fall Outside the Image

The imfilter function normally fills in these off-the-edge image pixels by assuming that they are 0. This is called zero padding and is illustrated in the following figure.

When you filter an image, zero padding can result in a dark band around the edge of the image, as shown in this example.

I = imread('eight.tif');

h = ones(5,5) / 25;

I2 = imfilter(I,h);

imshow(I), title('Original Image');

figure, imshow(I2), title('Filtered Image with Black Border')

To eliminate the zero-padding artifacts around the edge of the image, imfilter offers an alternative boundary padding method called border replication. In border replication, the value of any pixel outside the image is determined by replicating the value from the nearest border pixel. 

Replicated Boundary Pixels

To filter using border replication, pass the additional optional argument 'replicate' to imfilter.

I3 = imfilter(I,h,'replicate');

figure, imshow(I3);

title('Filtered Image with Border Replication')

The imfilter function supports other boundary padding options, such as 'circular' and 'symmetric'. See the reference page for imfilter for details.

Multidimensional Filtering

The imfilter function can handle both multidimensional images and multidimensional filters. A convenient property of filtering is that filtering a three-dimensional image with a two-dimensional filter is equivalent to filtering each plane of the three-dimensional image individually with the same two-dimensional filter. This example shows how easy it is to filter each color plane of a truecolor image with the same filter:

1. Read in an RGB image and display it.

2.  rgb = imread('peppers.png');

imshow(rgb);

3. Filter the image and display it.

4.  h = ones(5,5)/25;

5.  rgb2 = imfilter(rgb,h);

figure, imshow(rgb2)

Relationship to Other Filtering Functions

MATLAB^® has several two-dimensional and multidimensional filtering functions. The function filter2 performs two-dimensional correlation, conv2 performs two-dimensional convolution, and convn performs multidimensional convolution. Each of these filtering functions always converts the input to double, and the output is always double. These other filtering functions always assume the input is zero padded, and they do not support other padding options.

Filtering an Image with Predefined Filter Types

The fspecial function produces several kinds of predefined filters, in the form of correlation kernels. After creating a filter with fspecial, you can apply it directly to your image data using imfilter. This example illustrates applying an unsharp masking filter to a grayscale image. The unsharp masking filter has the effect of making edges and fine detail in the image more crisp.

I = imread('moon.tif');

h = fspecial('unsharp');

I2 = imfilter(I,h);

imshow(I), title('Original Image')

figure, imshow(I2), title('Filtered Image')

Designing Linear Filters in the Frequency Domain

FIR Filters

The Image Processing Toolbox™ software supports one class of linear filter: the two-dimensional finite impulse response (FIR) filter. FIR filters have a finite extent to a single point, or impulse. All the Image Processing Toolbox filter design functions return FIR filters.

FIR filters have several characteristics that make them ideal for image processing in the MATLAB^® environment:

• FIR filters are easy to represent as matrices of coefficients.
• Two-dimensional FIR filters are natural extensions of one-dimensional FIR filters.
• There are several well-known, reliable methods for FIR filter design.
• FIR filters are easy to implement.
• FIR filters can be designed to have linear phase, which helps prevent distortion.

Another class of filter, the infinite impulse response (IIR) filter, is not as suitable for image processing applications. It lacks the inherent stability and ease of design and implementation of the FIR filter. Therefore, this toolbox does not provide IIR filter support.

Frequency Transformation Method

The frequency transformation method transforms a one-dimensional FIR filter into a two-dimensional FIR filter. The frequency transformation method preserves most of the characteristics of the one-dimensional filter, particularly the transition bandwidth and ripple characteristics. This method uses a transformation matrix, a set of elements that defines the frequency transformation.

The toolbox function ftrans2 implements the frequency transformation method. This function's default transformation matrix produces filters with nearly circular symmetry. By defining your own transformation matrix, you can obtain different symmetries. (See Jae S. Lim, Two-Dimensional Signal and Image Processing, 1990, for details.)

The frequency transformation method generally produces very good results, as it is easier to design a one-dimensional filter with particular characteristics than a corresponding two-dimensional filter. For instance, the next example designs an optimal equiripple one-dimensional FIR filter and uses it to create a two-dimensional filter with similar characteristics. The shape of the one-dimensional frequency response is clearly evident in the two-dimensional response.

b = remez(10,[0 0.4 0.6 1],[1 1 0 0]);

h = ftrans2(b);

[H,w] = freqz(b,1,64,'whole');

colormap(jet(64))

plot(w/pi-1,fftshift(abs(H)))

figure, freqz2(h,[32 32])

Frequency Sampling Method

The frequency sampling method creates a filter based on a desired frequency response. Given a matrix of points that define the shape of the frequency response, this method creates a filter whose frequency response passes through those points. Frequency sampling places no constraints on the behavior of the frequency response between the given points; usually, the response ripples in these areas. (Ripples are oscillations around a constant value. The frequency response of a practical filter often has ripples where the frequency response of an ideal filter is flat.)

The toolbox function fsamp2 implements frequency sampling design for two-dimensional FIR filters. fsamp2 returns a filter h with a frequency response that passes through the points in the input matrix Hd. The example below creates an 11-by-11 filter using fsamp2 and plots the frequency response of the resulting filter. (The freqz2 function in this example calculates the two-dimensional frequency response of a filter.

Hd = zeros(11,11); Hd(4:8,4:8) = 1;

[f1,f2] = freqspace(11,'meshgrid');

mesh(f1,f2,Hd), axis([-1 1 -1 1 0 1.2]), colormap(jet(64))

h = fsamp2(Hd);

figure, freqz2(h,[32 32]), axis([-1 1 -1 1 0 1.2])

Notice the ripples in the actual frequency response, compared to the desired frequency response. These ripples are a fundamental problem with the frequency sampling design method. They occur wherever there are sharp transitions in the desired response.

You can reduce the spatial extent of the ripples by using a larger filter. However, a larger filter does not reduce the height of the ripples, and requires more computation time for filtering. To achieve a smoother approximation to the desired frequency response, consider using the frequency transformation method or the windowing method.

Windowing Method

The windowing method involves multiplying the ideal impulse response with a window function to generate a corresponding filter, which tapers the ideal impulse response. Like the frequency sampling method, the windowing method produces a filter whose frequency response approximates a desired frequency response. The windowing method, however, tends to produce better results than the frequency sampling method.

The toolbox provides two functions for window-based filter design, fwind1 and fwind2. fwind1 designs a two-dimensional filter by using a two-dimensional window that it creates from one or two one-dimensional windows that you specify. fwind2 designs a two-dimensional filter by using a specified two-dimensional window directly.

fwind1 supports two different methods for making the two-dimensional windows it uses:

• Transforming a single one-dimensional window to create a two-dimensional window that is nearly circularly symmetric, by using a process similar to rotation
• Creating a rectangular, separable window from two one-dimensional windows, by computing their outer product

The example below uses fwind1 to create an 11-by-11 filter from the desired frequency response Hd. The example uses the Signal Processing Toolbox hamming function to create a one-dimensional window, which fwind1 then extends to a two-dimensional window.

Hd = zeros(11,11); Hd(4:8,4:8) = 1;

[f1,f2] = freqspace(11,'meshgrid');

mesh(f1,f2,Hd), axis([-1 1 -1 1 0 1.2]), colormap(jet(64))

h = fwind1(Hd,hamming(11));

figure, freqz2(h,[32 32]), axis([-1 1 -1 1 0 1.2])

Creating the Desired Frequency Response Matrix

The filter design functions fsamp2, fwind2, and fwind2 all create filters based on a desired frequency response magnitude matrix. Frequency response is a mathematical function describing the gain of a filter in response to different input frequencies.

You can create an appropriate desired frequency response matrix using the freqspace function. freqspace returns correct, evenly spaced frequency values for any size response. If you create a desired frequency response matrix using frequency points other than those returned by freqspace, you might get unexpected results, such as nonlinear phase.

For example, to create a circular ideal lowpass frequency response with cutoff at 0.5, use

[f1,f2] = freqspace(25,'meshgrid');

Hd = zeros(25,25); d = sqrt(f1.^2 + f2.^2) < 0.5;

Hd(d) = 1;

mesh(f1,f2,Hd)

Ideal Circular Lowpass Frequency Response

Note that for this frequency response, the filters produced by fsamp2, fwind1, and fwind2 are real. This result is desirable for most image processing applications. To achieve this in general, the desired frequency response should be symmetric about the frequency origin (f1 = 0, f2 = 0).

Computing the Frequency Response of a Filter

The freqz2 function computes the frequency response for a two-dimensional filter. With no output arguments, freqz2 creates a mesh plot of the frequency response. For example, consider this FIR filter,

h =[0.1667    0.6667    0.1667

    0.6667   -3.3333    0.6667

    0.1667    0.6667    0.1667];

This command computes and displays the 64-by-64 point frequency response of h.

freqz2(h)

Frequency Response of a Two-Dimensional Filter

To obtain the frequency response matrix H and the frequency point vectors f1 and f2, use output arguments

[H,f1,f2] = freqz2(h);

freqz2 normalizes the frequencies f1 and f2 so that the value 1.0 corresponds to half the sampling frequency, or π radians.

For a simple m-by-n response, as shown above, freqz2 uses the two-dimensional fast Fourier transform function fft2. You can also specify vectors of arbitrary frequency points, but in this case freqz2 uses a slower algorithm.