calculating a 2D Savitzky-Golay kernel

I've been using Savitzky-Golay (SG) filters for quite a while, but always on 1D data, using an implementation based on the algorithm from the Numerical Recipes to calculate the kernel (mask, coefficients, however you wish to call them).

Apparently this filter is also useful for image processing (something I do not yet have a lot of experience with), so I thought I'd make myself a Core Image Unit to see what one can do with this filter in this new context. I haven't found a good, versatile enough C or C++ implementation to generate a 2D SG kernel allowing to specify width, order and derivative at will, so I'm using my trusty old function to generate a 1D kernel, and then make a 2D kernel out of that.My simple (naive?) idea was to rotate the 1D kernel to generate a 2D one. Then I realised that the kernel isn't necessarily symmetrical around its centre point, so the way and direction I rotate it might well have a meaning.What I create now is a 2D array in which the centre row corresponds to the 1D kernel (i.e. from left to right), and the rest is generated by rotating that vector counter-clockwise (over just less than 180°), presuming that this would correspond to a space in which the positive quadrant is the upper-right one. Thus, the centre column also corresponds to the 1D kernel, but now reading from bottom to top.

Is this approach correct? I calculate the 1D kernel to fill an array of a length that corresponds to the 2D kernel's diagonal. During the "rotation", I simply calculate an offset in the 1D kernel by applying Pythagoras to the 2D co-ordinates and rounding off to the nearest integer value. Should I instead do an interpolation, e.g. using a piece-wise linear fit?

% Map the local coordinates (0, 0) of "D" to an index. This index % corresponds to the row in "X" which evaluates the polynomial at (0, 0). index = sub2ind(size(D), yextent(1)+1, xextent(1)+1);

% Instead of evaluating the polynomial for all window pixels, we only % evaluate it at the origin. We also exploit "x=0" and "y=0" to simplify % the calculation. % Y = X * A; Y = Y(index); Y = X(index,1) * A(1);

Thanks - I'd seen that reference, but I'm afraid I cannot make much of it (my algebra has gone more than rusty, and my Matlab experience is mostly limited to Simulink )Even if I presume that the filter calculated is for a 0-order derivative (pure smoothing), I do not understand why there are multiple matrices being returned.The implementation from the Numerical Recipes uses a single (1D) vector for smoothing 1D data through a simple convolution, while the Microsoft reference gives at least 3 2D matrices.

About your implementation, is C the coefficient matrix (kernel) I was expecting to obtain? Guess I'll have to see if my Octave install is still recent enough

The basic premise of the Savitzky-Golay Filter is to perform local polynomial regression to smooth out a data set. So for some data point x, we look at its neighbors and fit a polynomial to the local data points. Then we just evaluate the polynomial at x. It's easy to see that we can control the order of the polynomial and the extent of the neighborhood (window size).http://en.wikipedia.org/wiki/Savitzky%E2%80%93Golay_filter_for_smoothing_and_differentiation

Fitting a polynomial involves finding the coefficients which allow the polynomial to best fit our specified data points. This can be performed using Linear Least Squares -- Using a little Linear Algebra, we can "easily" find our coefficients with a simple matrix equation (the Pseudo-Inverse). Also note that finding the derivative of a polynomial is pretty simple, so the Savitzky-Golay Filter is also a nice way to evaluating numerical derivatives.http://en.wikipedia.org/wiki/Linear_least_squares_%28mathematics%29

Note that there is no explicit restriction on the dimensionality of our data set. If we have a 1-dimensional data set, we need 1-dimensional polynomials. If we have a 2-dimensonal data set, we need 2-dimensional polynomials. If we have an n-dimensional data set, we need an n-dimensional polynomial. However, there is a restriction that the data points be regularly spaced.

http://research.microsoft.com/en-us/um/people/jckrumm/SavGol/SavGol.htmThe author of the Microsoft reference extends Savitzky-Golay filtering into two dimensions (i.e. images) using a 2-dimensional polynomial. He begins his formulation by serializing the window data points (2-dimensional) into a 1-dimensional vector. So for a 5x5 window D, we get a vector d of 25 data points.

The next step is formulating Eq(1) which states that multiplying the polynomial matrix X with the coefficients a should give us our data points d. We want to solve this equation for the coefficients a. In general, this equation is not solvable because there may not be a perfect polynomial fit -- So we use Linear Least Squares to find the best fit using the Pseudo-Inverse C.

The inner product of the i'th row in C with the vector d will give us the i'th coefficient. This gives us a. Now we can plug a into Eq(1) to evaluate the polynomial for all points in d.

Remember that we serialized D into d so we having a mapping from d back to D. We only need the polynomial evaluated at (0, 0) in D so we can use this mapping to see find the k'th value in d which maps to (0, 0).

At this point, we're pretty much done, but we can try to be a little more clever. We notice that (x=0, y=0) implies that we only need the coefficient a00. To see this, substitute the (x=0, y=0) back into the polynomial. In fact, a00 is also the result of evaluating our polynomial at (x=0, y=0). So we really only need to find the inner product of the k'th row in C with the vector d.

We can use the mapping of d back to D to reshape each row of C into matrices C00, C10, C01, etc. This is nice because we can see that C00 is a convolution kernel which gives us a00 for each data point -- We can perform 2D Savitzky-Golay filtering on an image by using C00 as a convolution kernel!

The basic premise of the Savitzky-Golay Filter is to perform local polynomial regression to smooth out a data set. So for some data point x, we look at its neighbors and fit a polynomial to the local data points. Then we just evaluate the polynomial at x. It's easy to see that we can control the order of the polynomial and the extent of the neighborhood (window size).

Yeah, I know, it's a bit like fitting a cubic spline, but without the constraint of going through all the original points.

Quote:

The author of the Microsoft reference extends Savitzky-Golay filtering into two dimensions (i.e. images) using a 2-dimensional polynomial. He begins his formulation by serializing the window data points (2-dimensional) into a 1-dimensional vector. So for a 5x5 window D, we get a vector d of 25 data points.

Does he say anything about the boundary problem that arises (or rather, gets suppressed)? Points (5,1) and (1,2) become neighbours in that approach, while they're not supposed to be.

Quote:

We can perform 2D Savitzky-Golay filtering on an image by using C00 as a convolution kernel!

Hmmm, this is what I thought I'd understood ... but if I look at the values in the X0Y0 arrays in SavGol.cpp, they don't correspond to the values calculated by the 1D Numerical Recipe. Should your implementation give the 1D coefficients when one feeds it a 1D image and height=1 ?

Does he say anything about the boundary problem that arises (or rather, gets suppressed)? Points (5,1) and (1,2) become neighbours in that approach, while they're not supposed to be.

The points are neighbors in the d vector, but this doesn't matter. It only matters that we correctly and consistently map d to X and D. The polynomial is fitted using the local coordinates of each element in d. So as long as our local coordinates are correct, the arrangement of d will not matter.

Quote:

Quote:

We can perform 2D Savitzky-Golay filtering on an image by using C00 as a convolution kernel!

Hmmm, this is what I thought I'd understood ... but if I look at the values in the X0Y0 arrays in SavGol.cpp, they don't correspond to the values calculated by the 1D Numerical Recipe. Should your implementation give the 1D coefficients when one feeds it a 1D image and height=1 ?

So yes, I believe my implementation will give the 1-dimensional coefficients if you provide a 1-dimensional "image".

Krumm's (Microsoft) coefficients will not match the 1-dimensional coefficients. His coefficients are fitted to data points from square windows with sizes greater than 1. This illustrates why we can't just "rotate" 1-dimensional coefficients to get 2-dimensional coefficients.

Krumm's coefficients also will not match my implementation. This is because he builds his polynomial to exclude some of the coupling terms while my implementation includes all coupling terms. Basically, for the same order, we are using different polynomials. We could modify my polynomials to match Krumm.

Pity Krumm only considers the "0th derivative" case ... though I guess that it must be quite straightforward to include support for arbitrary derivatives given a sufficient grasp of the math involved...

Pity Krumm only considers the "0th derivative" case ... though I guess that it must be quite straightforward to include support for arbitrary derivatives given a sufficient grasp of the math involved...

So to evaluate the k'th partial derivative with respect to x for an image, we can simply use Ck0 as a convolution kernel then multiply everything by k! (factorial). The same process follows for the k'th partial derivative with respect to y, except we use C0k.

We just need to be careful to select the correct Ck0 or C0k -- This will depend on how we ordered our polynomial terms. For example, in my example implementation, a 2'nd order polynomial is arranged as:

Code:

1, y, y^2, x, x*y, x*y^2, x^2, x^2*y, x^2*y^2

So the 1'st derivative with respect to y requires C01 which corresponds to the 2'nd term:

Code:

C01 = C(:,:,2);

And the 1'st derivative with respect to x requires C10 which corresponds to the 4'th term:

Heh, thanks for more math lessons (Now I really should find some to time sit down as long as required to wrap my head around it again ... but holidays first I'm afraid ...)

Am I right looking at your latest post that it won't be that easy to implement support for generic derivatives in your function? And that applying the filter means convoluting with Ck0 and then C0k (order shouldn't matter)?

Am I right looking at your latest post that it won't be that easy to implement support for generic derivatives in your function? And that applying the filter means convoluting with Ck0 and then C0k (order shouldn't matter)?

Hmm -- I think my posts may have been a little confusing. Some clarifications...

The example implementation should demonstrate that C can be constructed solely from the specified order and window. Similar to Krumm, the example implementation should also demonstrate that C can be rearranged to give us C00, C01, C10, etc. The rearrangement should be performed carefully. Remember that C** is the convolution kernel which gives us a**(i,j) when applied to data point I(i,j).

Previous posts (and Krumm) should show that:- The filtered value for data point I(i,j) is simply a00(i,j).- The k'th partial derivative w.r.t. x for data point I(i,j) is simply k!*ak0(i,j).- The k'th partial derivative w.r.t. y for data point I(i,j) is simply k!*a0k(i,j).And for application:- For a given polynomial order and window, we can easily construct C00, C01, C10, etc. (convolution kernels).- We can evaluate the k'th partial derivative w.r.t. x by convolving image I with Ck0 and multiplying by k! (factorial).- We can evaluate the k'th partial derivative w.r.t. y by convolving image I with C0k and multiplying by k! (factorial).- We can evaluate the filter as a 0'th partial derivative.