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)];
  hx = h(1);
  hy = h(2);

%  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;
elseif 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;

if length(h)<3, 

   my_sum = 0;
   for i=1:n;
      if strcmp(k,'knorm')==1;
         temp = knorm((x-A(i,1)*ones(size(x)))/hx,(y-A(i,2)*ones(size(y)))/hy);
      elseif strcmp(k,'kepan')==1;
         temp = kepan((x-A(i,1)*ones(size(x)))/hx,(y-A(i,2)*ones(size(y)))/hy);
      elseif strcmp(k,'krect')==1;
         temp = krect((x-A(i,1)*ones(size(x)))/hx,(y-A(i,2)*ones(size(y)))/hy);
      elseif strcmp(k,'ktria')==1;
         temp = ktria((x-A(i,1)*ones(size(x)))/hx,(y-A(i,2)*ones(size(y)))/hy);
      elseif strcmp(k,'kbiwe')==1;
         temp = kbiwe((x-A(i,1)*ones(size(x)))/hx,(y-A(i,2)*ones(size(y)))/hy);
      elseif strcmp(k,'ktriw')==1;
         temp = ktriw((x-A(i,1)*ones(size(x)))/hx,(y-A(i,2)*ones(size(y)))/hy);
      end;
      my_sum = my_sum+temp;
   end;

%  my_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
%    my_sum=my_sum+feval(k,(x-A(i,1)*ones(size(x)))/hx,(y-A(i,2)*ones(size(y)))/hy);
%  end;

  z=my_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 = sqrtm(H);
   H1= eye(size(H1,1))/H1;
   my_sum = 0;
   for i=1:n;

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

      if i==1;
         [p,q]=size(v1);
      end;

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


      if strcmp(k,'knorm')==1;
         temp = knorm(w1,w2);
      elseif strcmp(k,'kepan')==1;
         temp = kepan(w1,w2);
      elseif strcmp(k,'krect')==1;
         temp = krect(w1,w2);
      elseif strcmp(k,'ktria')==1;
         temp = ktria(w1,w2);
      elseif strcmp(k,'kbiwe')==1;
         temp = kbiwe(w1,w2);
      elseif strcmp(k,'ktriw')==1;
         temp = ktriw(w1,w2);
      end;
      my_sum = my_sum+temp;
   end;




%    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);
    
%    my_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);

%      my_sum=my_sum+feval(k,w1,w2);
%    end;

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

end;

