Savitzky-Golay Filters for 2D Images

John Krumm
Microsoft Research
Microsoft Corporation
Redmond, WA 98052

August 2001

Table of Contents
1. Introduction
2. Example of How to Compute Two-Dimensional Savitzky-Golay Filters
3. Example of How to Use Two-Dimensional Savitzky-Golay Filters
4. MatLab Routine for Computing Savitzky-Golay Filters
5. MatLab Routine for Checking Savitzky-Golay Filters
6. C Program with Savitzky-Golay Filters

1. Introduction

Savitzky-Golay filters have been popularized by the Numerical Recipes books. The one-dimensional filters presented in the book and the original papers are used for smoothing one-dimensional, tabulated data and also for computing numerical derivatives. The fundamental idea is to fit a different polynomial to the data surrounding each data point. The smoothed points are computed by replacing each data point with the value of its fitted polynomial. Numerical derivatives come from computing the derivative of each fitted polynomial at each data point.

While fitting polynomials for these purposes is obvious, the surprising part is that the polynomial coefficients can be computed with a linear filter. For smoothing, only one coefficient of the polynomial is needed, so the whole process of least squares fitting at every point becomes a simple process of applying the appropriate linear filter at every point. The Numerical Recipes books give a description of the filters and source code (FORTRAN and C)  for computing these one-dimensional filters for both smoothing and numerical derivatives.

The Savitzky-Golay filters are just as useful for image processing, where the idea is to fit a two-dimensional polynomial to a two-dimensional subsection of the image for smoothing and numerical derivatives. While the Numerical Recipes books give a way of computing the filters for one-dimensional data, this article shows how to compute two-dimensional filters for image data. It gives a MatLab routine for computing the filters and C source code for the actual filter coefficients for different image patch sizes and different polynomial orders.

back to top

2. Example of How to Compute Two-Dimensional Savitzky-Golay Filters

As an example, suppose we want to smooth or differentiate an image based on 5x5 image patches. The image patch is laid out as follows:

  xi
-2 -1 0 1 2
yi -2 d(0) d(1) d(2) d(3) d(4)
-1 d(5) d(6) d(7) d(8) d(9)
0 d(10) d(11) d(12) d(13) d(14)
1 d(15) d(16) d(17) d(18) d(19)
2 d(20) d(21) d(22) d(23) d(24)

Table 1: Patch of image data d(i) with its local coordinate system.

d(i) is the pixel value, and the column vector d represents all the image data, i.e.

d = (d(0) d(1) ... d(24))T

We want to fit a 3rd order, two-dimensional polynomial to this data. The polynomial is

d(i) f(xi,yi) = a00 + a10xi + a01yi + a20xi2+ a11xiyi + a02yi2 + a30xi3 + a21xi2yi + a12xiyi2 + a03yi3

Note that the coefficient of xiyj is aij. (xi,yi) is the pixel coordinate of d(i).

To compute the coefficients from the data we set up a matrix equation:

Xa = d

(1)

where

X = 1 x0 y0 x02 x0y0 y02 x03 x02y0 x0y02 y03
1 x1 y1 x12 x1y1 y12 x13 x12y1 x1y12 y13
: : : : : : : : : :
1 x24 y24 x242 x24y24 y242 x243 x242y24 x24y242 y243

and a is the vector of polynomial coefficients:

a = (a00 a10 a01 a20 a11 a02 a30 a21 a12 a03)T

Equation (1) simply reproduces the polynomial for each pixel in the image patch. We solve for the polynomial coefficients using least squares:

a = (XTX)-1XTd

C = (XTX)-1XT is the pseudo-inverse of X, and it is independent of the image data. Each polynomial coefficient is computed as the inner product of one row of C and the column of pixel values d. This is the surprising part about Savitzky-Golay filters: the polynomial coefficients are computed using a linear filter on the data. Just as one can reassemble d back into a rectangular patch of pixels, one can also assemble each row of C into the same size rectangle to get a traditional-looking image filter.

To complete this example, here are the filters for the first three polynomial coefficients. We will use the naming convention that coefficient aij is computed from rectangular filter Cij. The first three filters are:

C00 = -0.0743 0.0114 0.0400 0.0114 -0.0743
0.0114 0.0971 0.1257 0.0971 0.0114
0.0400 0.1257 0.1543 0.1257 0.0400
0.0114 0.0971 0.1257 0.0971 0.0114
-0.0743 0.0114 0.0400 0.0114 -0.0743

 

C10 = 0.0738 -0.0119 -0.0405 -0.0119 0.0738
-0.1048 -0.1476 -0.1619 -0.1476 -0.1048
0 0 0 0 0
0.1048 0.1476 0.1619 0.1476 0.1048
-0.0738 0.0119 0.0405 0.0119 -0.0738

 

C01 = 0.0738 -0.1048 0 0.1048 -0.0738
-0.0119 -0.1476 0 0.1476 0.0119
-0.0405 -0.1619 0 0.1619 0.0405
-0.0119 -0.1476 0 0.1476 0.0119
0.0738 -0.1048 0 0.1048 -0.0738

back to top

3. Example of How to Use Two-Dimensional Savitzky-Golay Filters

Suppose we want to smooth an image. Using a Savitzky-Golay filters, we are conceptually fitting a two-dimensional polynomial to the image patch surrounding each pixel and then evaluating this polynomial. The local coordinate system that we use for the image patch (see Table 1) has (x,y) = (0,0) at the pixel of interest in the middle of the patch. Thus to compute the smoothed value of the pixel, we just evaluate the polynomial at (x,y) = (0,0). This turns out to be merely a00, which we can compute by applying filter C00 to the image patch.

Suppose we want to compute partial derivatives on the patch. The two partial derivatives of the fitted polynomial are

fx(xi,yi) = a10 + 2a20xi + a11yi + 3a30xi2 + 2a21xiyi + a12yi2

fy(xi,yi) = a01+ a11xi + 2a02yi + a21xi2 + 2a12xiyi + 3a03yi2

Evaluating at (x,y) = (0,0), the results are simply fx(0,0) = a10 and fy(0,0) = a01, which are computed with filters C10 and C01 above.

back to top

4. MatLab Routine for Computing Savitzky-Golay Filters

Generalizing the example of Section 2, the following MatLab function called SavGol.m, generates the Savitzky-Golay filters for a user-specified square patch size and polynomial order. This function requires the MatLab Symbolic Math Toolbox. Note that this routine does not check to make sure the patch size is large enough to fit the requested polynomial order. Note also that for square patches with even dimensions, the origin of the local patch coordinates is at the center of the center four pixels.

SavGol.m

back to top

5. MatLab Routine for Checking Savitzky-Golay Filters

The following MatLab routine, called SavGolTest.m, takes the output of SavGol and tests it on an image patch made from a polynomial whose coefficients are all one.

SavGolTest.m

back to top

6. C Program with Savitzky-Golay Filters

The following C code gives the two-dimensional Savitzky-Golay filters for different patch sizes and different polynomial orders. The filters are stored in one-dimensional arrays in row-major order. The code was generated using the SavGol() MatLab routine and another MatLab routine to format the numbers into array declarations.

SavGol.cpp

back to top