function [x, y, zave, zmedian, zmin, zmax, zstd, n] = hist2d( xi, yi, zi, ni, mi )

% HIST2D populates 2d histogram
%   [x, y, zave, zmedian, zmin, zmax, zstd, n] = hist2d( xi, yi, zi, ni, mi )
%   ni & mi either number of divisions in xi and yi or
%   [xmin xmax xdelta] and [ymin ymax ydelta] respectively.
%
% J. G. Bellingham
% Feb 18, 2014

ok = ~isnan(xi.*yi.*zi);
xi = xi(ok);
yi = yi(ok);
zi = zi(ok);

% Set up binning space for histogram
if size(ni) == 1
    temp = sort(xi);
    ntemp = length(xi);
    xmin = temp(round(ntemp*0.05));
    xmax = temp(round(ntemp*0.95));
    xdelta = (xmax - xmin)/(ni);
else
    xmin = ni(1);
    xdelta = ni(2);
    xmax = ni(3);
end;

if size(mi) == 1
    temp = sort(yi);
    ntemp = length(yi);
    ymin = temp(round(ntemp*0.05));
    ymax = temp(round(ntemp*0.95));
    ydelta = (ymax - ymin)/(mi);    
else
    ymin = mi(1);
    ydelta = mi(2);
    ymax = mi(3);
end;

% Only consider data 'in bounds'
ok = xi <= xmax & xi >= xmin & yi <= ymax & yi >= ymin;
xi = xi(ok);
yi = yi(ok);
zi = zi(ok);

% Make sure we know number of x and y bins...
ni = ceil((xmax- xmin)/xdelta+1);
mi = ceil((ymax- ymin)/ydelta+1);
xmax = (ni-1)*xdelta + xmin; % In case xmax not an integral number of dx over xmin
ymax = (mi-1)*ydelta + ymin; % In case ymax not an integral number of dy over ymin

% Now make easily searchable 
ix = floor((xi - xmin)/xdelta)+1;
iy = floor((yi - ymin)/ydelta);
ii = round(ix+ni*iy);

n(1:ni, 1:mi) = 0;
zave(1:ni, 1:mi) = nan;
zmax = zave;
zmin = zave;
zstd = zave;
zmedian = zave;

% Compute histogram & keep statistics
for i=1:(ni)*(mi)
    zbin = zi(ii == i);
    n(i) = length(zbin);
    if n(i) > 0
        zave(i) =  sum(zbin)/n(i);
        zmedian(i) = median(zbin);
        zsort = sort(zbin);
        zmax(i) = zsort(round(n(i)*0.95));
        zmin(i) = zsort(round(n(i)*0.05)+1);
        % zmax(i) = max(zbin);
        % zmin(i) = min(zbin);
        zstd(i) = std(zbin);
    end;
end;

% Make this the right orientation...
n = n';
zave = zave';
zmedian = zmedian';
zmin = zmin';
zmax = zmax';
zstd = zstd';

% Creat grid to match...
x = xmin:xdelta:(xmax);
y = ymin:ydelta:(ymax);
[x, y] = meshgrid(x, y);
