]> Creatis software - CreaPhase.git/blobdiff - octave_packages/plot-1.1.0/hist2d.m
Add a useful package (from Source forge) for octave
[CreaPhase.git] / octave_packages / plot-1.1.0 / hist2d.m
diff --git a/octave_packages/plot-1.1.0/hist2d.m b/octave_packages/plot-1.1.0/hist2d.m
new file mode 100644 (file)
index 0000000..bc7e33b
--- /dev/null
@@ -0,0 +1,50 @@
+## Copyright (C) 2006 Paul Kienzle <pkienzle@users.sf.net>
+## This program is in the public domain
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{counts}, @var{xbins}, @var{ybins}] =} hist2d ([@var{x}, @var{y}], @var{xbins}, @var{ybins}, @var{norm})
+## Produce a 2D histogram.
+##
+## Points xi,yi are stored in a 2-column array.
+## If ybins is missing, use xbins.
+## If bins is a scalar, use that many bins.
+## If bins is a vector, it represents bin edges.
+##
+## @end deftypefn
+
+function [ret_counts, xbins, ybins] = hist2d(M,xbins,ybins)
+
+  if nargin < 1 && nargin > 3
+    print_usage
+  end
+
+  lo = min(M);
+  hi = max(M);
+  if nargin == 1
+    ybins = xbins = 10;
+  elseif nargin == 2
+    ybins = xbins;
+  endif
+
+  # If n bins, find centers based on n+1 bin edges
+  if isscalar(xbins) 
+    xbins = linspace(lo(1),hi(1),xbins+1);
+    xbins = (xbins(1:end-1)+xbins(2:end))/2;
+  end
+  if isscalar(ybins)
+    ybins = linspace(lo(2),hi(2),ybins+1); 
+    ybins = (ybins(1:end-1)+ybins(2:end))/2;
+  end
+
+  xcut = (xbins(1:end-1)+xbins(2:end))/2;
+  ycut = (ybins(1:end-1)+ybins(2:end))/2;
+  xidx = lookup(xcut,M(:,1))+1;
+  yidx = lookup(ycut,M(:,2))+1;
+  counts = sparse(xidx,yidx,1,length(xbins),length(ybins),'sum');
+
+  if nargout
+    ret_counts = full(counts');
+  else
+    mesh(xbins,ybins,full(counts'));
+  end
+endfunction