Never mind, the error is in the linspace range. Here's the correct code to show 2d hist.
function hist3_test
close all; clc;
events = 10000;
mytest = 0;
if (mytest)
xrange=0:.1:1;
nxs=length(xrange);
x=[]; y=[];
for xx=0:.1:1
xs=xx+sqrt(0.005)*randn(round(events/(nxs-1)),1);
if (xx<.3)
ys=cos(xs)+ sqrt(0.001)*randn(round(events/(nxs-1)),1);
else
ys=cos(xs)+ sqrt(0.1)*randn(round(events/(nxs-1)),1);
end
x = [x; xs];
y = [y; ys];
end
else
events = 1000000;
x1 = sqrt(0.05)*randn(events,1)-0.5; x2 = sqrt(0.05)*randn(events,1)+0.5;
y1 = sqrt(0.05)*randn(events,1)+0.5; y2 = sqrt(0.05)*randn(events,1)-0.5;
x= [x1;x2]; y = [y1;y2];
end
% For linearly spaced edges:
figure;
subplot(2,1,1);
plot(x,y,'.');
% correct the linspace range
ax=axis;
xedges = linspace(ax(1),ax(2),64); yedges = linspace(ax(3),ax(4),64);
histmat = hist2(x, y, xedges, yedges);
subplot(2,1,2);
% pcolor(xedges,yedges,histmat'); colorbar ; axis square tight ;
izero = histmat==0;
histmat(izero) = nan;
h=pcolor(xedges,yedges,histmat); colorbar ; axis square tight ;
set(h,'edgecolor','none');
end
function histmat = hist2(x, y, xedges, yedges)
% function histmat = hist2(x, y, xedges, yedges)
%
% Extract 2D histogram data containing the number of events
% of [x , y] pairs that fall in each bin of the grid defined by
% xedges and yedges. The edges are vectors with monotonically
% non-decreasing values.
%
%EXAMPLE
%
% events = 1000000;
% x1 = sqrt(0.05)*randn(events,1)-0.5; x2 = sqrt(0.05)*randn(events,1)+0.5;
% y1 = sqrt(0.05)*randn(events,1)+0.5; y2 = sqrt(0.05)*randn(events,1)-0.5;
% x= [x1;x2]; y = [y1;y2];
%
%For linearly spaced edges:
% xedges = linspace(-1,1,64); yedges = linspace(-1,1,64);
% histmat = hist2(x, y, xedges, yedges);
% figure; pcolor(xedges,yedges,histmat'); colorbar ; axis square tight ;
%
%For nonlinearly spaced edges:
% xedges_ = logspace(0,log10(3),64)-2; yedges_ = linspace(-1,1,64);
% histmat_ = hist2(x, y, xedges_, yedges_);
% figure; pcolor(xedges_,yedges_,histmat_'); colorbar ; axis square tight ;
% University of Debrecen, PET Center/Laszlo Balkay 2006
% email: balkay@pet.dote.hu
if nargin ~= 4
error ('The four input arguments are required!');
return;
end
if any(size(x) ~= size(y))
error ('The size of the two first input vectors should be same!');
return;
end
[xn, xbin] = histc(x,xedges);
[yn, ybin] = histc(y,yedges);
%xbin, ybin zero for out of range values
% (see the help of histc) force this event to the
% first bins
xbin((xbin == 0)) = 1;
ybin((ybin == 0)) = 1;
xnbin = length(xedges);
ynbin = length(yedges);
if xnbin >= ynbin
xy = ybin*(xnbin) + xbin;
indexshift = xnbin;
else
xy = xbin*(ynbin) + ybin;
indexshift = xnbin;
end
%[xyuni, m, n] = unique(xy);
xyuni = unique(xy);
hstres = histc(xy,xyuni);
histmat = zeros(xnbin,ynbin);
histmat(xyuni-indexshift) = hstres;
histmat = histmat';
end
