function [z,h]=kdest2d(x,y,A,k,h)

% KDEST2D     Bi-variate Kernel Density Estimator with generalised smoothing.
%
%             Examples: load kddata
%                       z=kdest2d(0,0,lukesh1(:,1:2),'knorm',0) 
% 
%             KDEST2D(X,Y,A,K,H) evaluates a bi-variate Kernel Density 
%             Estimate with kernel given by the function in K.
%             'KNORM' is the multivariate normal distribution.
%
%             (X,Y) contains the point(s) for which the density 
%             estimate is required.  X and Y can be scalars, vectors 
%             or matrices.
%
%             The bi-variate data points are stored row by row in A.
%
%             The parameter H represents the window width:
%             (1) If a single non-zero value of H is given then smoothing is the 
%             same in all directions.
%             (2) If H=[H1, H2] then smoothing is by H1 in the x direction and H2
%             in the y direction.  If H=0 is specified, automatic 'optimal' values 
%             are used. 
%             (3) If H=[H1, H2, H3] then the smoothing can be in orientations
%             different from those of the co-ordinate directions.  Note: the 
%             size of H3 must be less than H1*H2. 
%
%             Christian C. Beardah 1994

if nargin==4,
  h=0;
end;

temp=h;

[n d]=size(A); % Establish dimensions of A, 
               % n=number of data points,
               % d=dimension of the data. 

if nargin<4 | nargin>5,
  error('Incorrect number of input arguments');
end;

if min(h)==0,  % If no value of h has been specified by the user then 
               % calculate optimal values.

  h=[hns(A(:,1),k), hns(A(:,2),k)];

%  h=[hdpi(A(:,1),k,2), hdpi(A(:,2),k,2)];

%  h=[hste(A(:,1),k), hste(A(:,2),k)];

  if length(temp)==2,
    h=temp;
    if temp(1)==0,
      h(1)=hx;
    end;
    if temp(2)==0,
      h(2)=hy;
    end;
  end;
end;

if length(h)==1,
  hx=h;
  hy=h;
  else if length(h)==2,
         hx=h(1);
         hy=h(2);    
         else
           if abs(h(3))>=h(1)*h(2),
             error('Orientation parameter too large');
             else, 
               H=[h(1)^2 h(3); h(3) h(2)^2];
           end; 
       end;
end;

if length(h)<3, 

  sum=feval(k,(x-A(1,1)*ones(size(x)))/hx,(y-A(1,2)*ones(size(y)))/hy);

  for i=2:n,     % Sum over all data points
    sum=sum+feval(k,(x-A(i,1)*ones(size(x)))/hx,(y-A(i,2)*ones(size(y)))/hy);
  end;

  z=sum/(n*hx*hy); % z holds the kernel density estimates for the points in (x,y)

  else, % If h has three elements (fully general)

    H1=inv(sqrtm(H));

    v1=x-A(1,1)*ones(size(x));
    v2=y-A(1,2)*ones(size(y));

    [p,q]=size(v1);

    w=H1*[v1(:)';v2(:)'];
    w1=reshape(w(1,:),p,q);
    w2=reshape(w(2,:),p,q);
    
    sum=feval(k,w1,w2);

    for i=2:n,     % Sum over all data points
      v1=x-A(i,1)*ones(size(x));
      v2=y-A(i,2)*ones(size(y));

      w=H1*[v1(:)';v2(:)'];
      w1=reshape(w(1,:),p,q);
      w2=reshape(w(2,:),p,q);

      sum=sum+feval(k,w1,w2);
    end;

    z=sum/(n*sqrt(det(H))); % z holds the kernel density estimates for the points in (x,y)

end;

