function [counts,lims] = hist2(vx,vy,bins,contourLevels)

% Usage: [counts,lims] = hist2(vx,vy,bins,contourLevels=15)
% Alternative: hist2(vx,vy,bins) just draws the histogram
% vx,vy: equal-length vectors
% bins: number of intervals
% lims: first row x, second row y
% count: bins x bins, x is horiz dim min at left, y is vert dim min at top
% contourLevels: number of levels in contour plot (optional)

% Created by: Kristian Kirk
% Last modified: 26 Aug 2004

min_x = min(vx);
max_x = max(vx);
min_y = min(vy);
max_y = max(vy);
range_x = max_x-min_x;
range_y = max_y-min_y;

nvx = (vx-min_x)/range_x;
nvy = (vy-min_y)/range_y;

lims_x = linspace(min_x,max_x,bins);
lims_y = linspace(min_y,max_y,bins);
lims = [lims_x; lims_y];
counts = zeros(bins);

N = length(vx);

for n = 1 : N
    binx = max(1,ceil(bins*nvx(n)));
    biny = max(1,ceil(bins*nvy(n)));
    counts(biny,binx) = counts(biny,binx) + 1;
end


if nargout == 0
    if nargin < 4
        contourLevels = 15;
    end
%    mesh(counts)
    %counts = min(prctile(counts(:),99),counts);
    [xgrid,ygrid] = meshgrid(lims_x,lims_y);
    levels = linspace(min(counts(:)),max(counts(:)),contourLevels);
%    contour(xgrid,ygrid,log(counts),contourLevels)
    contour(xgrid,ygrid,counts,contourLevels)
    %imagesc([min_x max_x]+range_x/bins/2*[1 -1], [min_y max_y]+range_y/bins/2*[1 -1], counts)
    axis xy
    %colorbar
    hold on
    xlabel('x');
    ylabel('y');
    % division grid:
% lims = [linspace(min_x,max_x,bins+1); linspace(min_y,max_y,bins+1)];
%     for i = 1:length(lims)
%         plot(lims(1,i)*[1 1],lims(2,[1 end]),'k')
%         plot(lims(1,[1 end]),lims(2,i)*[1 1],'k')
%     end
end