## Copyright (C) 2008, 2009, 2012 Moreno Marzolla, Universita' di Bologna
## <marzolla (at) cs.unibo.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})
##
## Interpolate 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
## uses an additional linear polynomial, see the reference for details.
##
## If there are @code{n} 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}^n \lambda_i \phi( d_i ) + a_0 + a_1 x + a_2 y $$
## @end tex
## @end iftex
## @ifinfo
## @example
## @group
##           -- n
## 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 point
## @code{[@var{x} @var{y}]} from point @code{[@var{x}(i)
## @var{y}(i)]}, and @code{phi} is an appropriate radial base function.
##
## 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)}).
## Note that the @var{method} parameter string is case sensitive.
##
## @var{x}, @var{y} and @var{t} must be vectors of length @math{n}.
##
## @var{xx} and @var{yy} are rectangular matrices of the same size, where @code{[@var{xx}(i,j)
## @var{yy}(i,j)]} is the coordinate of the @math{(i,j)}-th 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 ( marzolla (at) cs.unibo.it )
## Last updated: 2012/09/29
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" );
  ( length(x) == length(y) && length(y) == length(t) ) || ...
      error("x, y, t must have the same length");
  x = x(:)';
  y = y(:)';
  t = t(:)';
  ismatrix(xx) || error( "xx must be a matric" );
  ismatrix(yy) || error( "yy must be a matrix" );
  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 = length(x);
  A = zeros( n+3, n+3 );
  x_i = diag(x)*ones(n); # x_i(i,j) = x(i)
  x_j = x_i';            # x_j(i,j) = x(j)
  y_i = diag(y)*ones(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 implemented yet" );
  ## epsilon = 1e-6;
  ## result = 2*log( delta * D_2.^.5 / 2 + epsilon ) + expint_Ei( delta^2 * D_2 / 4) + e(1);
endfunction
