function Z = hist2D(X,Y,Mx,My,bw,k)

[nx ny]= size(Mx);
N = length(X);
if nargin<4, error('hist2D takes 4 arguments'), end
if nargin<5, bw = 1; end

Z = zeros(nx,ny);
if (nargin<6) | (k==0)
    dx = abs(Mx(1,2)-Mx(1,1));
    dy = abs(My(2,1)-My(1,1));
    px = bw*[-dx/2 -dx/2 dx/2  dx/2];
    py = bw*[-dy/2  dy/2 dy/2 -dy/2];
    for i=1:nx
        for j=1:ny
            cx = Mx(i,j)+px;
            cy = My(i,j)+py;
            [ic,oc] = inpolygon(X,Y,cx,cy);
            %if any(ic), error(' '), end
            if any(oc),
                ic = single(ic);
                ic(find(oc))=.5; 
            end
            Z(i,j) = sum(ic);
        end
    end
end

Z = Z/(dx*dy*sum(sum(Z)));