You don't need a bivariate histogram to fit the bivariate normal--just use the sample means and covariance matrix. Here's an example:
% Let's say your data are in an n,2 matrix called xy.
% Here is one randomly generated to use in the example.
muXY = [100, 200];
sigmaXY = [15^2, 5^2; 5^2, 20^2];
xy = mvnrnd(muXY,sigmaXY,10000);
% Here is your bivariate histogram:
figure; histogram2(xy(:,1),xy(:,2));
% Now estimate the parameters of the best-fitting Gaussian:
xybar = mean(xy);
xycovar=cov(xy);
% Plot the best-fitting bivariate pdf:
xsteps = min(xy(:,1)):1:max(xy(:,1)); % Adjust with step sizes appropriate for your
ysteps = min(xy(:,2)):1:max(xy(:,2)); % x and y values.
[X,Y] = meshgrid(xsteps,ysteps);
F = mvnpdf([X(:) Y(:)],xybar,xycovar); % Note that xybar and xycovar are used here.
F = reshape(F,length(ysteps),length(xsteps));
figure; surf(xsteps,ysteps,F);
caxis([min(F(:))-.5*range(F(:)),max(F(:))]);
xlabel('x'); ylabel('y'); zlabel('Probability Density');
