function hist = kerneldensity2(i_edges, j_edges, ij_pts, kernel)
	% hist = kerneldensity2(i_edges, j_edges, ij_pts, kernel)
	%
	% kernel density estimation of 2D histogram
	%
	% Input:
	%	i_edges	edges of bins for first dimension of histogram
	%	j_edges	edges of bins for second dimension of histogram
	%	ij_pts	2 x N array of data points to plot on histogram
	%	kernel	function handle to kernel density function
	%
	% Output:
	%	hist		Ni x Nj histogram smoothed with kernel density function
	%
	% kernel should have the following properties:
	%	kernel(0,0) is a local maximum
	%  int(int(kernel(x,y),x)y) = 1
	
	%% get bins
	i_bins		= (i_edges(1:end-1) + i_edges(2:end))/2;
	j_bins		= (j_edges(1:end-1) + j_edges(2:end))/2;
	Ni				= length(i_bins);
	Nj				= length(j_bins);
	[ib, jb]		= ndgrid(i_bins, j_bins);
	
	%% accumulate data
	assert(size(ij_pts,1) == 2);
	n_data		= size(ij_pts,2);
	
	hist_ary		= arrayfun(@(i, j) kernel(ib-i, jb-j), ...
									ij_pts(1,:), ij_pts(2,:), 'UniformOutput', false);
	hist_ary		= reshape(hist_ary, [1 1 n_data]);
	hist			= sum(cell2mat(hist_ary), 3);
end