## Copyright (C) 2008  Carlo de Falco, Massimiliano Culpo
##
## This file is part of 
##
##                   PTC - FEM Patches Toolbox for Octave
## 
##   PTC 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 2 of the License, or
##   (at your option) any later version.
## 
##   PTC 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 PTC; if not, see <http://www.gnu.org/licenses/>.
##
##   MAIN AUTHOR:
##   Culpo Massimiliano
##   Bergische Universität Wuppertal
##   Fachbereich C - Mathematik und Naturwissenschaften
##   Arbeitsgruppe für Angewandte Mathematik
##   D-42119 Wuppertal, Germany

## -*- texinfo -*-
## @deftypefn {Function File} {@var{zi} =} PTC2griddata (@var{x}, @
## @var{y}, @var{z}, @var{xi}, @var{yi}, @var{method}, @var{tri})
## 
## Generate a regular mesh from irregular data using interpolation.
## The function is defined by @code{@var{z} = f (@var{x}, @var{y})}.
## The interpolation points are all @code{(@var{xi}, @var{yi})}.  If
## @var{xi}, @var{yi} are vectors then they are made into a 2D mesh,
## unless a triangulation @var{tri} has been passed to the function.
##
## The interpolation method can be @code{"nearest"}, or
## @code{"linear"}. If method is omitted it defaults to @code{"linear"}.
## @seealso{delaunay}
## @end deftypefn
##
## This function has in major part been copied from octave "griddata"
## Original Author:	Kai Habel <kai.habel@@gmx.de>
## Adapted-by:  Alexander Barth <barth.alexander@@gmail.com>
##              xi and yi are not "meshgridded" if both are vectors 
##              of the same size (for compatibility) 
  

function [rz] = PTC2griddata (x, y, z, xi, yi, method, tri)

  if (ischar (method))
    method = tolower (method);
  endif
  
  if (! all (size (x) == size (y) & size (x) == size (z)))
    error ("PTC2griddata: x, y, and z must be vectors of same length");
  endif

  [nr, nc] = size (xi);
  
  zi = nan (size (xi));

  if (strcmp (method, "nearest"))
    ## search index of nearest point
    idx = dsearch (x, y, tri, xi, yi);
    valid = !isnan (idx);
    zi(valid) = z(idx(valid));

  elseif (strcmp (method, "linear"))
    ## search for every point the enclosing triangle
    tri_list= tsearch (x, y, tri, xi(:), yi(:));
    
    ## only keep the points within triangles.
    valid = !isnan (reshape (tri_list, size (xi)));
    tri_list = tri_list(!isnan (tri_list));
    nr_t = rows (tri_list);

    if nr_t##if one point is within trgs
      ## assign x,y,z for each point of triangle
      x1 = x(tri(tri_list,1));
      x2 = x(tri(tri_list,2));
      x3 = x(tri(tri_list,3));
      
      y1 = y(tri(tri_list,1));
      y2 = y(tri(tri_list,2));
      y3 = y(tri(tri_list,3));

      z1 = z(tri(tri_list,1));
      z2 = z(tri(tri_list,2));
      z3 = z(tri(tri_list,3));

      ## calculate norm vector
      N = cross ([x2-x1, y2-y1, z2-z1], [x3-x1, y3-y1, z3-z1]);
      N_norm = sqrt (sumsq (N, 2));
      N = N ./ N_norm(:,[1,1,1]);
    
      ## calculate D of plane equation
      ## Ax+By+Cz+D=0;
      D = -(N(:,1) .* x1 + N(:,2) .* y1 + N(:,3) .* z1);
    
      ## calculate zi by solving plane equation for xi,yi
      zi(valid) = -(N(:,1).*xi(valid) + N(:,2).*yi(valid) + D ) ./ N(:,3);
    endif
    
  else
    error ("PTC2griddata: unknown interpolation method");
  endif
  
  rz = zi;

endfunction

%!test
%xp  = [.25 .25 .25 .275 .275 .275 .3 .3 .3];
%yp  = [.25 .275 .3 .25 .275 .3 .25 .275 .3];
%zp  = [1 0 0 0 0 0 0 0 0];
%xc  = [0 0 0 .5 .5 .5 1 1 1];
%yc  = [0 .5 1 0 .5 1 0 .5 1];
%tri = [ 1   4   5\
%        2   5   6 \
%        4   7   8 \
%        5   8   9 \
%        1   5   2 \
%        2   6   3 \
%        4   8   5 \
%        5   9   6];
%rz  = PTC2griddata(xp,yp,zp,xc,yc,"linear",tri);
%assert(rz,nan(size(xc)));