Code covered by the BSD License  

Highlights from
histograms for ND data

image thumbnail

histograms for ND data

by

 

Create histograms for ND data using ND bins

histnd(varargin)
function [histmat,nOOF,OOFidx]  = histnd(varargin)
% histnd Histogram count of ND data with ND bins
%
% histmat = histnd(x, y, ..., xedges, yedges, ...)
% Extract ND histogram data containing the number of events
% of [x, y, ...] tuples that fall in each bin of the ND-grid defined by
% xedges, yedges, .... The edges are passed to histc internally and should
% therefore conform to histc's input restrictions: the edge-vectors should
% be monotonically non-decreasing.
%
% [histmat, nOOF, OOFidx] = histnd(x, y, ..., xedges, yedges, ...)
% If any values are outside of the range of the edges, they are not
% counted. The number of those cases and their linear index in the input
% data is however returned in the second and third output arguments.
%
%EXAMPLES
%
% 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 = histnd(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_ = histnd(x, y, xedges_, yedges_);
% figure; pcolor(xedges_,yedges_,histmat_'); colorbar ; axis square tight ;
%
%3D data
% x = 3.*randn(640000,1);
% y = 1.*randn(640000,1);
% z = 1.*randn(640000,1);
% histmat = histnd(x,y,z,linspace(min(x),max(x),20),linspace(min(y),max(y),20),linspace(min(z),max(z),20));
% % make 3D hist, color of points indicates count
% [xp,yp,zp] = meshgrid(linspace(min(x),max(x),20),linspace(min(y),max(y),20),linspace(min(z),max(z),20));
% % cut away histogram positions where count is 0
% qzero = histmat==0;
% histmat(qzero) = [];
% xp(qzero) = [];
% yp(qzero) = [];
% zp(qzero) = [];
% % draw points
% figure;%('Renderer','OpenGL') % might need the openGL renderer to handle so many points
% ax = scatter3(xp(:),yp(:),zp(:),'.');
% % color them according to count
% cdata = histmat(:)./max(histmat);
% set(ax,'CData',cdata);
% xlabel('X'), ylabel('Y'), zlabel('Z')

% Diederick C. Niehorster, University of Hong Kong, department of
% Psychology. 2011. email: dcnieho@gmail.com
%
% based on hist2 by:
% University of Debrecen, PET Center/Laszlo Balkay 2006
% email: balkay@pet.dote.hu

if mod(nargin,2)~=0
    error ('An even number of input arguments is required (you need to specify edges for each dimension!');
end
data = varargin(1:nargin/2);
edges= varargin(nargin/2+1:end);
% ensure data is column vectors
data = cellfun(@(x) x(:),data,'uni',false);

if any(diff(cellfun(@length,data)))
    error ('The sizes of all data input vectors should be same!');
end

% determine number of bins in each dimension
sz = cellfun(@length,edges);

% for each dimension, determine which bin the point falls in
binIdx = zeros(length(data{1}),length(data));
for p=1:length(data)
    [~,binIdx(:,p)] = histc(data{p},edges{p});
end

% binIdx zero for out of range values (see the help of histc)
% count these events and remove them
qOutOfRange = ~all(binIdx,2);
nOOF        = sum(qOutOfRange);
OOFidx      = find(qOutOfRange);
binIdx(qOutOfRange,:) = [];

% now determine which ND-bin each point falls in
binIdx   = num2cell(binIdx, 1);
ndBinIdx = sub2ind(sz, binIdx{:});

% count number per bin
histmat = histc(ndBinIdx, 1:prod(sz));

% create output
histmat = reshape(histmat, sz);

Contact us