## Copyright (C) 2008, 2009 Moreno Marzolla, INFN Padova
## (moreno.marzolla(at)pd.infn.it)
##
## 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{ti} =} griddata_cgs(@var{x}, @var{y}, @var{t}, @var{xx}, @var{yy}, @var{method}, @var{c})
##
## Perform interpolation of scattered data @code{[ @var{x}(i) @var{y}(i)
## @var{t}(i) ]} using continuous global surfaces, as described in
## Stephen D. Billings, Rick K. Beatson, and Garry N. Newsam,
## "Interpolation of geophysical data using continuous global surfaces".
## GEOPHYSICS, VOL. 67, NO. 6 (NOVEMBER-DECEMBER 2002); P. 1810--1822.
## This implementation supports thin-plate splines, multiquadric,
## gaussian and inverse_multiquadric radial base functions. This
## implementation performs an interpolation considering an additional
## linear polynomial.
##
## If there are @code{k} data points @code{[@var{x}(i) @var{y}(i) @var{t}(i)]},
## interpolation is done according to the function:
##
## @iftex
## @tex $$f(x,y) = \sum_{i=1}^k \lambda_i \phi( d_i ) + a_0 + a_1 x + a_2 y $$
## @end tex
## @end iftex
## @ifinfo
## @example
## @group
##           -- k
## f(x,y) =  >      lambda_i phi( d_i ) + a_0 + a_1 x + a_2 y 
##           -- i=1
## @end group
## @end example
## @end ifinfo
##
## where @code{d_i} is the euclidean distance of the point
## @code{[@var{x} @var{y}]} from the point @code{[@var{x}(i)
## @var{y}(i)]} and @code{phi} is an appropriate radial base function.
## Thie function computes the values for the unknown coefficients
## @code{lambda_i} and @code{a_0, a_1, a_2}.
##
## The interpolation method can be @code{"linear"} (in which case the
## @code{griddata()} octave function is used), @code{"thin_plate"}
## (which uses the base function @code{phi(d_i) = d_i^2 log(d+i +
## epsilon)}, with @code{epsilon = 1e-5}), @code{"multiquadric"} (base
## function @code{phi(d_i) = (d_i^2 + c^2) ^ 0.5}
## @code{"inverse_multiquadric"} (base function @code{(d_i^2 + c^2) ^
## -0.5}) and @code{"gaussian"} (base function @code{exp(-c * d_i^2)}).
##
## @var{x}, @var{y} and @var{t} must be (1xn) vectors;
##
## @var{xx} and @var{yy} are (m*n) matrices where @code{[@var{xx}(i,j)
## @var{yy}(i,j)]} is the coordinate of the (i,j) point in the grid. The
## result @code{@var{ti}(i,j)} is the interpolated value at the
## coordinates @code{[@var{xx}(i,j) @var{yy}(i,j)]}. The parameters are
## used as in the @code{griddata()} octave function.
##
## @var{c} is the radial function parameter (if any). Not used for
## linear interpolation and thin-plate splines; for other radial base
## functions, if omitted a default value is used.
##
## @seealso{griddata}
##
## @end deftypefn
##
## Author: Moreno Marzolla ( moreno.marzolla (at) pd.infn.it ) Created:
## 2008-07-19 Version: 0.1 Revision: $Revision: 71 $
function ti = griddata_cgs( x, y, t, xx, yy, method = "linear", c )
  if ( nargin < 5 )
    print_usage();
    return;
  endif
  isvector(x) || error( "x must be a vector" );
  isvector(y) || error( "y must be a vector" );
  isvector(t) || error( "t must be a vector" );
  size(xx) == size(yy) || error( "xx and yy must be of equal size" );
  ischar(method) || error( "method must be a string" );
  
  if ( strcmp(method, "linear") )
    ti = griddata( x, y, t, xx, yy );
    return;
  elseif ( strcmp(method, "thin_plate") )
    phi = @thin_plate;
  elseif ( strcmp(method, "multiquadric") )
    if ( nargin < 7 )
      c = 0.2;
    endif
    phi = @(x) multiquadric(x,c);
  elseif ( strcmp(method, "gaussian") )
    if ( nargin < 7 )
      c = 0.2;
    endif
    phi = @(x) gaussian(x,c);
  elseif ( strcmp(method, "inverse_multiquadric") )
    if ( nargin < 7 )
      c = 0.2;
    endif
    phi = @(x) inverse_multiquadric(x,c);
  elseif ( strcmp(method, "tension_spline") )
    if ( nargin < 7 )
      c = 1.0;
    endif
    phi = @(x) tension_spline(x,c);
  else
    error( [ "Unrecognized interpolation method: " method ] );
  endif
  n = max(size(x));
  x = reshape( x, [1 n] );
  y = reshape( y, [1 n] );
  t = reshape( t, [1 n] );
  A = zeros( n+3, n+3 );
  x_i = dmult(x,ones(n,n));
  x_j = x_i';
  y_i = dmult(y,ones(n,n));
  y_j = y_i';

  d_2 = (x_i-x_j).^2 + (y_i-y_j).^2; # d_2(i,j) = (x(i)-x(j)^2 + (y(i)-y(j))^2
  A(1:n,1:n) = phi( d_2 );
  A(n+1,1:n) = 1;
  A(1:n,n+1) = 1;
  A(n+2,1:n) = x;
  A(n+3,1:n) = y;
  A(1:n,n+2) = x';
  A(1:n,n+3) = y';
  t_prime = [t 0 0 0];
  lambda = t_prime / A;
  ti = zeros( size(xx) );
  F = zeros(1, n+3 );
  F(n+1) = 1;
  for i=1:size(xx,1)
    for j=1:size(xx,2)
      ## D_2(k) = (xx(i,j)-x(k))^2 + (yy(i,j)-y(k))^2;
      D_2 = (x - xx(i,j)).^2 + (y - yy(i,j)).^2;
      F(1:n) = phi( D_2 );
      F(n+2) = xx(i,j);
      F(n+3) = yy(i,j);
      ti(i,j) = dot(lambda, F);
    endfor
  endfor
endfunction
%!demo
%! x=4*rand(100,1)-2;
%! y=4*rand(size(x))-2;
%! z=sinc(2*(x.^2+y.^2));
%! [xx,yy]=meshgrid(linspace(-2,2,32));
%! zz = griddata_cgs(x,y,z,xx,yy);
%! zz2 = sinc(2*(xx.^2+yy.^2));
%! mesh(xx,yy,zz2);
%! title('Original profile');
%! disp("Press Enter"); pause;
%! mesh(xx,yy,zz);
%! title('nonuniform grid sampled at 100 points, linear');
%! disp("Press Enter"); pause;
%! tic(); zz = griddata_cgs(x,y,z,xx,yy,"thin_plate"); toc();
%! mesh(xx,yy,zz);
%! title('nonuniform grid sampled at 100 points, thin-plate spline');
%! disp("Press Enter"); pause;
%! tic(); zz = griddata_cgs(x,y,z,xx,yy,"multiquadric",0.4); toc();
%! mesh(xx,yy,zz);
%! title('nonuniform grid sampled at 100 points, multiquadric (c=0.4)');

%!demo
%! x=4*rand(500,1)-2;
%! y=4*rand(size(x))-2;
%! z=sinc(2*(x.^2+y.^2));
%! [xx,yy]=meshgrid(linspace(-2,2,32));
%! zz = griddata_cgs(x,y,z,xx,yy);
%! zz2 = sinc(2*(xx.^2+yy.^2));
%! mesh(xx,yy,zz2);
%! title('Original profile');
%! disp("Press Enter"); pause;
%! mesh(xx,yy,zz);
%! title('nonuniform grid sampled at 500 points, linear');
%! disp("Press Enter"); pause;
%! tic(); zz = griddata_cgs(x,y,z,xx,yy,"thin_plate"); toc();
%! mesh(xx,yy,zz);
%! title('nonuniform grid sampled at 500 points, thin-plate spline');
%! disp("Press Enter"); pause;
%! tic(); zz = griddata_cgs(x,y,z,xx,yy,"multiquadric",0.4); toc();
%! mesh(xx,yy,zz);
%! title('nonuniform grid sampled at 500 points, multiquadric (c=0.4)');

## Thin-plate base function
function result = thin_plate( D_2 )
  epsilon = 1e-5;
  result = D_2 .* log( D_2.^0.5 + epsilon );
endfunction

## Gaussian base function
function result = gaussian( D_2, c )
  result = exp( -c * D_2 );
endfunction

## Multiquadric base function
function result = multiquadric( D_2, c )
  result = (D_2 + c^2) .^ 0.5;
endfunction

## Inverse multiquadric base function
function result = inverse_multiquadric( D_2, c )
  result = (D_2 + c^2) .^ -0.5;
endfunction

## Tension spline base function
function result = tension_spline( D_2, delta )
  error( "tension spline interpolation not working yet" );
  ## epsilon = 1e-6;
  ## result = 2*log( delta * D_2.^.5 / 2 + epsilon ) + expint_Ei( delta^2 * D_2 / 4) + e(1);
endfunction