## Copyright (C) 2015 Beke J.
## 
## This program is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
## 
## This program 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 General Public License for more details.
## 
## You should have received a copy of the GNU General Public License
## along with this program.  If not, see <http://www.gnu.org/licenses/>.

## -*- texinfo -*- 
## @deftypefn {Function File} {@var{ZZ} =} griddataIDW (@var{x} ,@var{y} ,@var{z} ,@var{XX} ,@var{YY} ,@var{p} ,@var{wfunc} )
##
##  This function implements to interpolate irregular 2D data to a 
##  regular 2D grid using inverse distance weighting.
##
##  Each new point is calculated as a weighted average of the inverse distance (to power p)
##  of the known function values. The larger p, the less influence from far points.
##  A custom weighting function can be provided. Distance^p will be given as input.
##
##  @var{x} vector with known x values
##
##  @var{y} vector with known y values
##
##  @var{z} vector with known z values
##
##  @var{XX}, @var{YY} grid to interpolate values
##
##  
##  @var{p} (optional, default p=4) power of the distance.
##
##  @var{wfunc} (optional, default wfunc = @@(x) 1./x with x the distance^p)
##               Optional weighting function which will get distance^p as input.
##               The function must be an element-by-element function.
##
## @seealso{griddata, meshgrid}
## @seealso{@url{http://en.wikipedia.org/wiki/Inverse_distance_weighting}}
##
## @end deftypefn

## Author: Beke J.
## Created: 2015-01-06

function [ZZ] = griddataIDW (x, y, z, XX, YY, p, wfunc)
  #http://en.wikipedia.org/wiki/Inverse_distance_weighting
  
  # check input
  if (nargin < 5 || nargin > 7)
    error("usage: [ZZ] = griddataIDW (x, y, z, XX, YY, p, wfunc)");
  end
  if (nargin < 6)
    p=4;
  end
  if (nargin < 7)
    wfunc = @(x) 1./x;
  end
  if (!isvector(x) || !isvector(y) || !isvector(z))
    error("x, y and z must be vectors");
  end
  if(!isscalar(p) || iscomplex(p))
    error("p must be scalar");
  end
  if (!ismatrix(XX) || !ismatrix(YY))
    error("XX, YY must be matrices");
  end
  
  #check sizes
  if (!( (size(x)==size(y)) && (size(x)==size(z)) ))
    error("x, y and z must be vectors of the same size");
  end
  if (!(size(XX)==size(YY)))
    error("XX and YY must be matrices of the same size");
  end
  
  # Actual function:
  for row=1:rows(XX)
    for col=1:columns(XX)
    
      xi = XX(row,col);
      yi = YY(row,col);
      
      distances = ((x.-xi).^2 .+(y.-yi).^2).^(p/2);
      w = wfunc(distances);
      ZZ(row,col) = sum(z.*w)/sum(w);
      
    end
   end
  
end


%!demo
%!    #normal test function
%!	[XX,YY]=meshgrid(0:0.1:1,0:0.1:1);
%!	
%!	func = @(x,y) 0.75./exp((x*.5).^2.*(y.*5).^2);
%!	
%!	func = @(x,y) 0.75./exp((x*.5).^2.*(y.*5).^2);
%!	
%!	ZZ = func(XX,YY);
%!	
%!	#random points of test function
%!	x=rand(1,100);
%!	y=rand(1,100);
%!	z=func(x,y);
%!	
%!	newplot();
%!	hold on;
%!	subplot (1, 2, 1);
%!	[c,h]=contourf(XX,YY,ZZ,10);
%!	axis equal;
%!	#colorbar ();
%!	clabel (c, h);
%!	
%!	
%!	p=4;
%!  [ZZi] = griddataIDW (x, y, z, XX, YY, p);
%!	subplot (1, 2, 2);
%!	hold on;
%!	[c,h]=contourf(XX,YY,ZZi,10);
%!	plot(x,y,"+")
%!	axis equal;
%!	#colorbar ();
%!	clabel (c, h);
%!

%!test
%!	x = 1.0;
%!	y = 4.0;
%!	
%!	x1 = 0;
%!	y1 = 0;
%!	z1 = 5.0;
%!	d1 =sqrt((x-x1)^2+(y-y1)^2);
%!	
%!	x2 = 1;
%!	y2 = 0;
%!	z2 = 6.0;
%!	d2 =sqrt((x-x2)^2+(y-y2)^2);
%!	
%!	x3 = 0;
%!	y3 = 1;
%!	z3 = 7.0;
%!	d3 =sqrt((x-x3)^2+(y-y3)^2);
%!
%!	p = 4.0;
%!	wfunc = @(x) 1./x;
%!	z = (z1*wfunc(d1^p)+z2*wfunc(d2^p)+z3*wfunc(d3^p))/(wfunc(d1^p)+wfunc(d2^p)+wfunc(d3^p));
%!	ZZ = griddataIDW([x1,x2,x3], [y1,y2,y3], [z1,z2,z3], x, y);
%!	assert(abs(z-ZZ)<1.0e-10)
%!
%!	p = 1.0;
%!	wfunc = @(x) exp(-x);
%!	z = (z1*wfunc(d1^p)+z2*wfunc(d2^p)+z3*wfunc(d3^p))/(wfunc(d1^p)+wfunc(d2^p)+wfunc(d3^p));
%!	ZZ = griddataIDW([x1,x2,x3], [y1,y2,y3], [z1,z2,z3], x, y, p, wfunc);
%!	assert(abs(z-ZZ)<1.0e-10)
%!

%!test
%!	p = 1.0;
%!	wfunc = @(x) 1./x;
%!	z = (0.5/0.75+1/0.25)/(1/0.25+1/0.75);
%!	ZZ = griddataIDW([0,1], [0, 0], [0.5,1], 0.75, 0, p, wfunc);
%!	assert(abs(z-ZZ)<1.0e-10)