function [x,y,z] = xyzbin(varargin)
%Bins scatter points according to a grid specified by a grid
%
% Usage:
%   [x,y,z] = xyzbin(...)
%   [x,y,z] = xyzbin('infile',infile,'outfile',outfile,'dx',dx)
%   [x,y,z] = xyzbin('x',xin,'y',yin,'z',zin,'dx',dx)
%   [x,y,z] = xyzbin('x',xin,'y',yin,'z',zin,'dx',dx,'dy',dy)
%   [x,y,z] = xyzbin('infile',infile,'dx',dx,'x0',x0,'y0',y0,'theta',theta)
%
% Description:
%   Averages the x, y, and z values of points within the same bin.
%   The grid used to define the bins may be user specified are calculated
%   based on the extent of the scatter points. 
%
%written by Alex Sanchez, USACE-CHL 2010

if(nargin==0)
    xyzbindemo    
    return    
end

theta = 0;
for i=1:2:length(varargin)
    eval([varargin{i},'=varargin{',num2str(i+1),'};'])
end

if(exist('infile','var'))
    xyz = load(infile);
    x = xyz(:,1);
    y = xyz(:,2);
    z = xyz(:,3);
end
if(~exist('x','var') || ~exist('y','var') || ~exist('z','var'))
    if(~exist('xyz','var'))
        error('Invalid input. Missing input file name or x,y,z values.')
    else
        x = xyz(:,1);
        y = xyz(:,2);
        z = xyz(:,3);
    end
end
ind = isfinite(z);
xi = x(ind);
yi = y(ind);
zi = z(ind);

if(~exist('x0','var'))
    x0 = min(xi);
end
if(~exist('y0','var'))
    y0 = min(yi);
end
if(~exist('dx','var'))
    dx = mean(diff(sort(x)));
end
if(~exist('dy','var'))
    dy = dx;
end
if(~exist('nx','var'))
    nx = (max(xi)-x0)/dx+1;
else
    nx= 100;
end
if(~exist('ny','var'))
    ny = (max(yi)-y0)/dy+1;
else
    ny = 100;
end

%Translate and rotate interpolation coordinates
theta = theta*pi/180;
co = cos(-theta); %note negative angle
si = sin(-theta);
x2 = yi*si + xi*co - x0;
yi = yi*co - xi*si - y0;
xi = x2;

%Index of (xi,yi) to grid
i = xi/dx+1;
j = yi/dy+1;

%Points outside grid
out = (i>nx | i<1 | j>ny | j<1);

%Points on boundary
bi = (i==nx);
bj = (j==ny);
i(bi) = nx-1;
j(bj) = ny-1;
i(out) = [];
j(out) = [];

%Clip out bad points
x = xi(~out);
y = yi(~out);
z = zi(~out);

%Grid index, Note: Points with the same k are in the same grid cell
k = floor(j) + floor(i-1)*ny;
[k,ik] = sort(k);

%Sort points
x = x(ik);
y = y(ik);
z = z(ik);

%Index of points in same cells
ind = [0; (diff(round(k))==0); 0];

%Start and end of binned cells
s = find(ind(1:end-1)==0 & ind(2:end)==1);
e = find(ind(1:end-1)==1 & ind(2:end)==0);

%Points merged
n = sum(ind)+1;
if(n==0)
    disp('No points binned')
    x = xi;
    y = yi;
    z = zi;
    return
end
for i=1:length(s)
    z(e(i)) = mean(z(s(i):e(i)));
    x(e(i)) = mean(x(s(i):e(i)));
    y(e(i)) = mean(y(s(i):e(i)));
end
disp([num2str(n),' points binned'])
x = x(~ind(2:end));
y = y(~ind(2:end));
z = z(~ind(2:end));
x2 = x0 - y*si + x*co;
y = y0 + y*co + x*si;
x = x2;

if(exist('outfile','var'))
   fid = fopen(outfile,'w');
   fprintf(fid,'%f %f %f\n',[x,y,z]');
   fclose(fid);
end

return

%----------------------------------------------------------
function xyzbindemo
n = 100;
xi = abs(rand(n,1));
yi = abs(rand(n,1));
zi = 5-peaks(xi,yi) + 0.2*abs(rand(n,1));
dx = 0.1;
clc
[x,y] = xyzbin('x',xi,'y',yi,'z',zi,'dx',dx,'outfile','test.xyz');
clf
plot(xi,yi,'.b',x,y,'or')
return



