h=sgdf_2d(x,y,nx,ny,flag_coupling)

function h=sgdf_2d(x,y,nx,ny,flag_coupling)
%2-D Savitzky-Golay Differentiation Filter
%Here, the filter coefficients for the central point and the first order
%derivative (differentiation) is taken into account.
%Usage:
%       h=sgdf_2d(x,y,nx,ny,flag_coupling)
%       x    = x data point, e.g., -3:3
%       y    = y data point, e.g., -2:2 
%       nx    =x polynomial order       default=1              
%       ny    =y polynomial order       default=1
%       flag_coupling  = with or without the consideration of the coupling terms, between x and y. default=0
%Example:
%       sgdf_2d(-3:3,-4:4,2,2)
%       sgdf_2d(-3:3,-4:4,2,3,1)     
%Author:
%       Jianwen Luo <[e-mail], [e-mail]> 2004-10-31
%       Department of Biomedical Engineering
%       Tsinghua University, Beijing 100084, P. R. China  
%Reference
%[1]A. Savitzky and M. J. E. Golay, "Smoothing and Differentiation of Data by Simplified Least Squares Procedures," 
%   Analytical Chemistry, vol. 36, pp. 1627-1639, 1964.
%[2]K. L. Ratzlaff and J. T. Johnson, "Computation of Two-Dimensional Polynomial Least-Squares Convolution Smoothing Integers," 
%   Analytical Chemistry, vol. 61, pp. 1303-1305, 1989.
%[3]J. E. Kuo, H. Wang, and S. Pickup, "Multidimensional Least-Squares Smoothing Using Orthogonal Polynomials," 
%   Analytical Chemistry, vol. 63, pp. 630-635, 1991.
%[4]http://research.microsoft.com/users/jckrumm/SavGol/SavGol.htm
%[5] J. W. Luo, K. Ying, P. He and J. Bai, 
% Properties of Savitzky-Golay digital differentiators, 
%  Digital Signal Processing, in press, 
%  Available at http://www.sciencedirect.com/science?_ob=IssueURL&_tockey=%23TOC%236768%239999%23999999999%2399999%23FLP%23Articles_in_Press&_auth=y&view=c&_acct=C000053663&_version=1&_urlVersion=0&_userid=1553430&md5=6a3653781dc7d2ec0b942b6efd816a32


if nargin<5
    flag_coupling=0;
end

if nargin<4
    ny=1;
end

if nargin<3
    nx=1;
end

[x,y]=meshgrid(x,y);
[ly,lx]=size(x);

if(nx>lx-1) || (ny>ly-1)
    error('polynomial order too large!');
end

x=x(:);
y=y(:);

A=[];

if flag_coupling
    for kx=nx:-1:0 %nx,nx-1,...,x,1 (nx+1)
        for ky=ny:-1:0 %ny,ny-1,...,y,1 (ny+1)
            A=[A x.^kx.*y.^ky];
        end
    end     
else
    for k=nx:-1:1
        A=[A x.^k];
    end
    for k=ny:-1:0
        A=[A y.^k];
    end  
end

h=inv(A'*A)*A';
if flag_coupling
    h=h(nx*(ny+1),:);% ?*x 
else    
    h=h(nx,:);% ?*x 
end
h=reshape(h,ly,lx);