function ZI = griddata_average(X,Y,Z,XI,YI)
%GRIDDATA_AVERAGE maps (bins) non-uniform data to regular grid by averaging
%   
%   A main advantage over griddata is that data is not interpolated, and
%   thus there is no 'fantasy' data in the output. If there is no data in a
%   gridcell, the output value is simply NaN. If there are multiple points
%   in a gridcel, the output is calculated by averaging (no weighting) of
%   the datapoints within that cell.
%   
%   The function only works when XI and YI are generated by meshgrid or
%   similar functions.
%    
%   Syntax:
%   ZI = griddata_averaging(X,Y,Z,XI,YI)
%
%   Input:
%   varargin  =
%
%   Output:
%   varargout =
%
%   Example 1: 
%     binning of high res scattered data (as e.g. LIDAR)
%
%     x = (rand(100)-.5).*6;
%     y = (rand(100)-.5).*6;
%     z = peaks(x,y)+(rand(100)-.5)*7;
%     nn = 18;
%     [XI,YI] = meshgrid(linspace(-4,4,nn),linspace(-4,4,nn));
%
%     subplot(2,1,1)
%     ZI = griddata_average(x,y,z,XI,YI);
%     surf(XI,YI,ZI)
%     title 'griddata average'
%
%     subplot(2,1,2)
%     ZI = griddata(x,y,z,XI,YI); 
%     surf(XI,YI,ZI)
%     title griddata
%     hold off
%
%   See also: GRIDDATA, GRIDDATA_NEAREST, GRIDDATA_AVERAGE, GRIDDATE_REMAP,
%   INTERP2, BIN2 

%% Copyright notice
%   --------------------------------------------------------------------
%   Copyright (C) 2010 <COMPANY>
%       tda
%
%       <EMAIL>	
%
%       <ADDRESS>
%
%   This library is free software: you can redistribute it and/or
%   modify it under the terms of the GNU Lesser General Public
%   License as published by the Free Software Foundation, either
%   version 2.1 of the License, or (at your option) any later version.
%
%   This library is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
%   Lesser General Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with this library. If not, see <http://www.gnu.org/licenses/>.
%   --------------------------------------------------------------------

% This tool is part of <a href="http://OpenEarth.nl">OpenEarthTools</a>.
% OpenEarthTools is an online collaboration to share and manage data and 
% programming tools in an open source, version controlled environment.
% Sign up to recieve regular updates of this function, and to contribute 
% your own tools.

%% Version <http://svnbook.red-bean.com/en/1.5/svn.advanced.props.special.keywords.html>
% Created: 26 Mar 2010
% Created with Matlab version: 7.10.0.499 (R2010a)

% $Id$
% $Date$
% $Author$
% $Revision$
% $HeadURL$
% $Keywords: $

%% vectorize data and sort along x (improves speed within loops)
x=X(:);
y=Y(:);
z=Z(:);

[x ind]=sort(x,'ascend');
y=y(ind);
z=z(ind);

%% remap x and y to rounded gridpoints
unique_XI=unique(XI)';
unique_YI=unique(YI)';

unique_XI = [2*unique_XI(1) - unique_XI(2) - 2*eps(unique_XI(1)) unique_XI 2*unique_XI(end) - unique_XI(end-1) + 2*eps(unique_XI(end))];
unique_YI = [2*unique_YI(1) - unique_YI(2) - 2*eps(unique_YI(1)) unique_YI 2*unique_YI(end) - unique_YI(end-1) + 2*eps(unique_YI(end))];

for ii = find(~ismember(x,unique_XI))' 
    [dummy,nn] = min(abs(x(ii)-unique_XI));
    x(ii) = unique_XI(nn);
end

for ii = find(~ismember(y,unique_YI))' 
    [dummy,nn] = min(abs(y(ii)-unique_YI));
    y(ii) = unique_YI(nn);
end

%% remove NaN's in input data
temp=any([isnan(x) isnan(y) isnan(z)],2);
x(temp)=[];
y(temp)=[];
z(temp)=[];

%% map points of z to full grid points of zn
ZI = nan(length(unique_YI)-2,length(unique_XI)-2);
unique_x = unique(x)';

for ix=2:length(unique_XI)-1 % skip first and last ... then step through each x
    if ismember(unique_XI(ix),unique_x) % returns 1 when to-grid x-value is present in the from-x-grid
        f=find(unique_XI(ix)==x,1,'first');
        l=find(unique_XI(ix)==x,1,'last');
        x_temp=x(f:l); % generate temporary, limited vectors to improve search speed of data
        y_temp=y(f:l); 
        z_temp=z(f:l);
        
        for iy=2:length(unique_YI)-1 % skip first and last
            temp=(y_temp==unique_YI(iy)&x_temp==unique_XI(ix));
            if any(temp)
                temp2=z_temp(temp);
                ZI(iy-1,ix-1)=sum(temp2)/length(temp2); % map the data, use average if more than one datapoint falls within grid
            end
        end
    end
end
