function [f,h]=kdest3d(x1,x2,x3,A,k,h)

% KDEST3D     Tri-variate Kernel Density Estimator with generalised smoothing.
%
%             Examples: load kddata
%                       f=kdest3d(0,0,0,lukesh1(:,1:3),'knorm',0) 
% 
%             KDEST3D(X1,X2,X3,A,K,H) evaluates a tri-variate Kernel Density 
%             Estimate with kernel given by the function in K.
%             'KNORM' is the multivariate normal distribution.
%
%             (X1,X2,X3) contains the point(s) for which the density 
%             estimate is required.  X1, X2 and X3 can be scalars, vectors 
%             or matrices.
%
%             The tri-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, H3] then smoothing is by HI in the xi direction.  
%             If H=0 is specified, automatic 'optimal' values are used. 
%
%             Christian C. Beardah 1995

if nargin==5 | h==0,
  h=[0 0 0];
end;

[r,c]=size(h);

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

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

if r==1 & min(h)==0,  % If no value of h has been specified by the user then 
                      % calculate optimal values.
  if h(1)==0,
    h(1)=hns(A(:,1),k);
  end;
  if h(2)==0,
    h(2)=hns(A(:,2),k);
  end;
  if h(3)==0,
    h(3)=hns(A(:,3),k);
  end;
  
end;

if r==1, 

  sum=feval(k,(x1-A(1,1)*ones(size(x1)))/h(1),(x2-A(1,2)*ones(size(x2)))/h(2),(x3-A(1,3)*ones(size(x3)))/h(3));

  for i=2:n,     % Sum over all data points
    sum=sum+feval(k,(x1-A(i,1)*ones(size(x1)))/h(1),(x2-A(i,2)*ones(size(x2)))/h(2),(x3-A(i,3)*ones(size(x3)))/h(3));
  end;

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

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

    H1=inv(sqrtm(h));
    
    v1=x1-A(1,1)*ones(size(x1));
    v2=x2-A(1,2)*ones(size(x2));
    v3=x3-A(1,3)*ones(size(x3));

    [p,q,s]=size(v1);

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

    for i=2:n,     % Sum over all data points
      v1=x1-A(i,1)*ones(size(x1));
      v2=x2-A(i,2)*ones(size(x2));
      v3=x3-A(i,3)*ones(size(x3));
      
      w=H1*[v1(:)';v2(:)';v3(:)'];
      w1=reshape(w(1,:),p,q,s);
      w2=reshape(w(2,:),p,q,s);
      w3=reshape(w(3,:),p,q,s);

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

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

end;