function [NamesCoefs, NamesTerms, XPow, YPow, SG] = SavGol(nOrder,nSize)
% Compute Savitzky-Golay coefficients
% John Krumm, Microsoft Research, August 2001
% Requires MatLabSymbolic Math Toolbox
% On return:
% NamesCoefs(i,:) gives the name of coefficient i, e.g. a23
% NamesTerms(i,:) gives the name of the polynomial term i, e.g. (x^2)(y^3)
% XPow(i) and YPow(i) give exponents on x and y for coefficient i
% SG(:,:,i) gives the nSize x nSize filter for computing coefficient i

% Set up polynomial terms for a given order
Terms = [];
NamesCoefs = [];
NamesTerms = [];
XPow = [];
YPow = [];
syms x y real;
for j=0:nOrder
for i=0:nOrder-j
% NamesTerms and NamesCoefs each must have strings all the same length
% Each string in NamesCoefs is 3 characters long
% Each string in NamesTerms is 10 characters long
% There will be a problem if nOrder >= 10
NamesCoefs = [NamesCoefs; sprintf('a%1d%1d',i,j)]; % must all be same length
if (i>0 & j>0)
NamesTerms = [NamesTerms; sprintf('(x^%d)(y^%d)',i,j)];
end
if (i>0 & j==0)
NamesTerms = [NamesTerms; sprintf('(x^%d) ',i)];
end
if (i==0 & j>0)
NamesTerms = [NamesTerms; sprintf('(y^%d) ',j)];
end
if (i==0 & j==0)
NamesTerms = [NamesTerms; sprintf('1 ')];
end
Terms = [Terms; (x^i)*(y^j)];
XPow = [XPow; i];
YPow = [YPow; j];
end
end

% Compute A matrix for a nSize x nSize window
A = [];
for y = -(nSize-1)/2:(nSize-1)/2 % important to loop through in same scan order as image patch pixels
for x = -(nSize-1)/2:(nSize-1)/2
%sprintf ('%f %f',x,y)
A = [A; subs(Terms')];
end
end

% Compute coefficient matrix
C = inv(A'*A)*A';

% Pull out coefficients
SG = [];
[nTerms, nDum] = size(Terms);
for i=1:nTerms
SG(:,:,i) = reshape(C(i,:),[nSize,nSize]);
end

% Print
%for i=1:nTerms
%     NamesCoefs(i,:)
%     NamesTerms(i,:)
%     SG(:,:,i)
%end