function h=kdint2(A,x,y,z)

% KDINT2      Calculate kde height using linear interpolation.
%
%             Example: load kddata
%                      [x,y,z]=kdfft2(mask4,'knorm',64,0);
%                      h=kdint2(mask4,x,y,z);
%                      [mask4,h] 
%
%             KDINT2(A,X,Y,Z)
%
%             Christian C. Beardah 1995 

z=z';

[n,d]=size(A);

h=zeros(n,1);

deltax=x(2)-x(1);
deltay=y(2)-y(1);

binx=ceil((A(:,1)-x(1))/deltax);
biny=ceil((A(:,2)-y(1))/deltay);

f1=find(binx==0);
for j=1:length(f1),
  binx(f1(j))=1;
end; 

f2=find(biny==0);
for j=1:length(f2),
  biny(f2(j))=1;
end; 

f5=find(binx>=1 & binx<length(x) & biny>=1 & biny<length(y))';

w=deltax*deltay;

w11=zeros(n,1);
w12=w11;
w21=w11;
w22=w11;

w11(f5)=(x(binx(f5)+1)'-A(f5,1)).*(y(biny(f5)+1)'-A(f5,2));
w12(f5)=(A(f5,1)-x(binx(f5))').*(y(biny(f5)+1)'-A(f5,2));  
w21(f5)=(x(binx(f5)+1)'-A(f5,1)).*(A(f5,2)-y(biny(f5))');
w22(f5)=(A(f5,1)-x(binx(f5))').*(A(f5,2)-y(biny(f5))');

z1=zeros(n,1);
z2=z1;
z3=z1;
z4=z1;

for i=1:length(f5),
  z1(f5(i),1)=z(binx(f5(i)),biny(f5(i)));
  z2(f5(i),1)=z(binx(f5(i))+1,biny(f5(i)));
  z3(f5(i),1)=z(binx(f5(i)),biny(f5(i))+1);
  z4(f5(i),1)=z(binx(f5(i))+1,biny(f5(i))+1);
end;

h=(w11.*z1+w12.*z2+w21.*z3+w22.*z4)/w;