function [Pxy, X, Y]=bs_sample_pdf_2d(x,y)

% Prepare X-axis
mmin=min(x(:));
mmax=max(x(:));
dX=(mmax-mmin)/60;
Xmedian=mean(x);
Xplus= Xmedian:dX:mmax+dX;
Xminus= Xmedian:-dX:mmin-dX;
Xminus=Xminus(end:-1:2);
X=[Xminus Xplus];

% Prepare Y-axis
mmin=min(y(:));
mmax=max(y(:));
dY=(mmax-mmin)/60;
Ymedian=mean(y);
Yplus= Ymedian:dY:mmax+dY;
Yminus= Ymedian:-dY:mmin-dY;
Yminus=Yminus(end:-1:2);
Y=[Yminus Yplus];

% produce the 2D axis
[X, Y]=meshgrid(X,Y);

[I,J]=size(X);
jpdf=zeros(I,J);

% count the elements falling in each bin
for i=1:I
  for j=1:J
     in=find ( x >= X(i,j)-dX/2 & x < X(i,j)+dX/2   & ...
               y >= Y(i,j)-dY/2 & y < Y(i,j)+dY/2)  ;
     jpdf(i,j) = length(in);
  end
end

% normalize
Pxy=jpdf/sum(jpdf(:).*dX*dY);