function [NN,map1,map2]=hist2d(X,m1spec,m2spec,varargin)
% hist2d - 2D histogram
%
% hist2d :: 
%    [[N,2]]	~'N rows of 2D data',
%    (K1:natural | dmap(K1)) ~'num bins for dim 1, or bins',
%    (K2:natural | dmap(K2)) ~'num bins for dim 1, or bins'
% -> [[K1,K2]->natural] ~'bin counts',
%    dmap(K1),
%    dmap(K2).
%
% if second map spec is empty, it means use the same map
% for both dimensions.

mins=min(X); mm=max(X); maxs=mm+eps(mm);
map1=mkdmap(mins(1),maxs(1),m1spec,[]);
map2=mkdmap(mins(2),maxs(2),m2spec,m1spec);

N=accumhist([map1(X(:,1)) map2(X(:,2))],1,[cardr(map1) cardr(map2)]);
if nargout==0,
	opts=prefs('cmap',@(t)log(1+t),varargin{:});
	imagexy(centres(map2),centres(map1),opts.cmap(N));
else
	NN=N;
end


function M=mkdmap(min,max,spec,def)
	if isempty(spec), spec=def; end
	if isa(spec,'dmap') % is alread a dmap
		M=spec;
	elseif length(spec)==1, % spec=number of bins
		M=linmap(min,max,spec);
	else
		error('invalid discretisation map');
	end

