FIXED (again): INTERP2 problem

[Pierre Baldensperger]

Thanks Paul,

Actually I'm pretty much a beginner with Octave / Matlab. At first, I didn't
manage to keep the vectorized syntax because expressions like a(yidx,xidx)
weren't interpreted the way I expected them to be. Actually I expected them
to return the vector of elements a(yidx(i),xidx(i)) for i=1:size(xidx), but
something else shows up that I don't really understand. Example:

a = [ 1 2 3 4 5 6 ; 7 8 9 10 11 12 ; 13 14 15 16 17 18 ; 19 20 21 22 23 24 ]
xidx = floor(rand(40,1)*6)+1
yidx = floor(rand(40,1)*4)+1
result = a(yidx,xidx)
size(result)

ans =
 40  40

What's this beast? If somebody has a teacher's soul and could explain to me
what this "something else" is, what this syntax actually does and what 
it can
be used for: your clarifications are welcome (RTFM links if applicable :-)!

So I reverted to the "for" loops... and indeed noticed the very poor (even
unusable) performance, which eventually made me give up Octave and try 
Scilab
instead (which looks very impressive BTW, and worked flawlessly after 
getting
used to a few syntax differences).

I was expecting somebody else to optimize that "interp2.m" for me, but since
you hurt my pride by pointing out my lousy programming ;-), I went back 
to it
and came up with a better (considerably faster) version using 1D-indexing
through the "sub2ind" function. Proceeding with the above example:

idx = sub2ind(size(a),yidx,xidx)
result = a(idx)

As mentioned in the previous posts: the behaviour of this version of INTERP2
with respect to XI and YI has been changed!! So if your code was 
assuming the
former behaviour (imposed meshgridded XI and YI), you may need to modify the
way you call the function. Also, this version may still be sub-optimal and
badly syntaxed, so before it goes into the repository, I'd prefer some 
"guru"
out there (Paul, hear my call?) to have a look at it and see :

- whether the initial 2D-indexing syntax could be restored while keeping the
 new behaviour with XI/YI
- whether it would be (performance-wise,...) meaningful to do so

The m-file is attached.

Regards,
-Pierre.





Paul Kienzle wrote:

> Thanks for your efforts.  As this is a community project it only
> improves when people write code to fix things and bother to tell
> others about it.
>
> I'm not going to post your new version to the repository
> because the for loops will kill performance on a large mesh.
> I don't see anything that would prohibit replacing them with
> the usual vectorization techniques, but I didn't look very closely.
>
> - Paul
>
> On Nov 24, 2005, at 4:58 PM, Pierre Baldensperger wrote:
>
> > Hello again,
> >
> > Just in case anybody runs into the same INTERP2 problem as mine (see 
> > previous
> > posts): the "interp2.m" routine (from octave-forge-2005-06-13) 
> > assumes that
> > interpolating points XI and YI are in "meshgrid" format (and if they are
> > vectors, the "meshgridding" is forced at the beginning of the 
> > routine). This
> > behaviour is _not_consistent_ with INTERP2's help message, nor with the
> > behaviour of its namesake routine in Matlab.
> >
> > To solve my problem, I finally modified the routine myself to make its
> > behaviour more consistent with what it should (IMO) be. I modified 
> > both the
> > "linear" and "nearest" methods. You will find a copy of the modified 
> > routine
> > attached to this post. I don't know how to submit it to the 
> > maintainer and
> > don't want to bother with that. So hopefully, if this fix is useful, 
> > somebody
> > will find it here, adapt it and transfer it to the main octave 
> > repository.
> >
> > Please be advised that this is a quick & dirty fix and probably not 
> > very smart
> > programming: use at your own risks, but at least it works for me. In 
> > case you
> > want to try it, remember that the behaviour of this routine wrt XI 
> > and YI has
> > changed: you may need to modify code that assumes the former behaviour.
> >
> > Regards,
> > -Pierre.
>
>
>
>
> -------------------------------------------------------------
> Octave is freely available under the terms of the GNU GPL.
>
> Octave's home on the web:  http://www.octave.org
> How to fund new projects:  http://www.octave.org/funding.html
> Subscription information:  http://www.octave.org/archive.html
> -------------------------------------------------------------

## Copyright (C) 2000  Kai Habel
##
## 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 2 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, write to the Free Software
## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

## -*- texinfo -*-
## @deftypefn {Function File} address@hidden interp2 (@var{x}, @var{y}, 
@var{z}, @var{xi}, @var{yi})
## @deftypefnx {Function File} address@hidden interp2 (@var{Z}, @var{xi}, 
@var{yi})
## @deftypefnx {Function File} address@hidden interp2 (@var{Z}, @var{n})
## @deftypefnx {Function File} address@hidden interp2 (... , '@var{method}')
## Two-dimensional interpolation. @var{X},@var{Y} and @var{Z} describe a
## surface function. If @var{X} and @var{Y} are vectors their length
## must correspondent to the size of @var{Z}. @var{X} and @var{Y} must be
## monotonic. If they are matrices they  must have the 'meshgrid' format. 
##
## ZI = interp2 (X, Y, Z, XI, YI, ...) returns a matrix corresponding
## to the points described by the matrices @var{XI}, @var{YI}. 
##
## If the last argument is a string, the interpolation method can
## be specified. At the moment only 'linear' and 'nearest' methods
## are provided. If it is omitted 'linear' interpolation  is 
## assumed.
##
## ZI = interp2 (Z, XI, YI) assumes X = 1:rows(Z) and Y = 1:columns(Z)
## 
## ZI = interp2 (Z, n) interleaves the Matrix Z n-times. If n is ommited
## n = 1 is assumed
##
## @seealso{interp1}
## @end deftypefn

## Author:      Kai Habel <address@hidden>
## 2005-03-02 Thomas Weber <address@hidden> 
##     * Add test cases
## 2005-03-02 Paul Kienzle <address@hidden>
##     * Simplify
## 2005-04-23 Dmitri A. Sergatskov <address@hidden>
##     * Modified demo and test for new gnuplot interface
## 2005-11-24 Pierre Baldensperger <address@hidden>
##     * Rather big modification (XI,YI no longer need to be
##       "meshgridded") to be consistent with the help message
##       above and also with the Matlab namesake.
## 2005-11-26 Pierre Baldensperger <address@hidden>
##     * Replaced performance-killing "for" loops (introduced
##       in the previous modification) with vector syntax.

function ZI = interp2 (varargin)
  Z = X = Y = xi = yi = [];
  n = 1;
  method = "linear";

  switch nargin
  case 1
    Z = varargin{1};
  case 2
    if isstr(varargin{2})
      [Z,method] = deal(varargin{:});
    else
      [Z,n] = deal(varargin{:});
    endif
  case 3
    if isstr(varargin{3})
      [Z,n,method] = deal(varargin{:});
    else
      [Z,XI,YI] = deal(varargin{:});
    endif
  case 4
    [Z,XI,YI,method] = deal(varargin{:});
  case 5
    [X,Y,Z,XI,YI] = deal(varargin{:});
  case 6 
    [X,Y,Z,XI,YI,method] = deal(varargin{:});
  otherwise
    usage ("interp2 (X, Y, Z, XI, YI, method)");
  endswitch

  % Type checking.
  if !ismatrix(Z), error("interp2 expected matrix Z"); endif
  if !isscalar(n), error("interp2 expected scalar n"); endif
  if !isnumeric(X) || !isnumeric(Y), error("interp2 expected numeric X,Y"); 
endif
  if !isnumeric(XI) || !isnumeric(YI), error("interp2 expected numeric XI,YI"); 
endif
  if !isstr(method), error("interp2 expected string 'method'"); endif

  % Define X,Y,XI,YI if needed
  [zr, zc] = size (Z);
  if isempty(X),  X=[1:zc]; Y=[1:zr]; endif
  if isempty(XI), p=2^n; XI=[p:p*zc]/p; YI=[p:p*zr]/p; endif
  
  % Are X and Y vectors? => produce a grid from them
  if isvector (X) && isvector (Y)
    [X, Y] = meshgrid (X, Y);
  elseif ! all(size (X) == size (Y))
    error ("X and Y must be matrices of same size");
  endif
  if any(size (X) != size (Z))
    error ("X and Y size must match Z dimensions");
  endif
  
  if any(size(XI) != size(YI))
    error ("XI and YI must be matrices of same size");
  endif
 
  shape = size(XI);
  XI = reshape(XI,1,prod(shape));
  YI = reshape(YI,1,prod(shape));

  xidx = lookup(X(1, 2:end-1), XI) + 1;
  yidx = lookup(Y(2:end-1, 1), YI) + 1;

  if strcmp (method, "linear")
    ## each quad satisfies the equation z(x,y)=a+b*x+c*y+d*xy
    ##
    ## a-b
    ## | |
    ## c-d
    a = Z(1:(zr - 1), 1:(zc - 1));
    b = Z(1:(zr - 1), 2:zc) - a;
    c = Z(2:zr, 1:(zc - 1)) - a;
    d = Z(2:zr, 2:zc) - a - b - c;

    idx = sub2ind(size(a),yidx,xidx);

    ## scale XI,YI values to a 1-spaced grid
    Xsc = (XI - X(1, xidx)) ./ (X(1, xidx + 1) - X(1, xidx));
    Ysc = (YI - Y(yidx, 1)') ./ (Y(yidx + 1, 1) - Y(yidx, 1))';

    ## apply plane equation
    ZI = a(idx) + b(idx).*Xsc + c(idx).*Ysc + d(idx).*Xsc.*Ysc;

  elseif strcmp (method, "nearest")
    xtable = X(1, :);
    ytable = Y(:, 1)';
    ii = (XI - xtable(xidx) > xtable(xidx + 1) - XI);
    jj = (YI - ytable(yidx) > ytable(yidx + 1) - YI);
    idx = sub2ind(size(Z),yidx+jj,xidx+ii);
    ZI = Z(idx);
  else
    error ("interpolation method not (yet) supported");
  endif

  ## set points outside the table to NaN

  ZI( XI < X(1,1) | XI > X(1,end) | YI < Y(1,1) | YI > Y(end,1) ) = NaN;

  ZI = reshape(ZI,shape(1),shape(2));


endfunction

%!demo
%! A=[13,-1,12;5,4,3;1,6,2];
%! x=[0,1,4]; y=[10,11,12];
%! xi=linspace(min(x),max(x),17);
%! yi=linspace(min(y),max(y),26);
%! mesh(xi,yi,interp2(x,y,A,xi,yi,'linear'));
%! [x,y] = meshgrid(x,y); 
%! __gnuplot_raw__ ("set nohidden3d;\n")
%! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;

%!demo
%! A=[13,-1,12;5,4,3;1,6,2];
%! x=[0,1,4]; y=[10,11,12];
%! xi=linspace(min(x),max(x),17);
%! yi=linspace(min(y),max(y),26);
%! mesh(xi,yi,interp2(x,y,A,xi,yi,'nearest'));
%! [x,y] = meshgrid(x,y); 
%! __gnuplot_raw__ ("set nohidden3d;\n")
%! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;

%!#demo
%! A=[13,-1,12;5,4,3;1,6,2];
%! x=[0,1,2]; y=[10,11,12];
%! xi=linspace(min(x),max(x),17);
%! yi=linspace(min(y),max(y),26);
%! mesh(xi,yi,interp2(x,y,A,xi,yi,'cubic'));
%! [x,y] = meshgrid(x,y); 
%! __gnuplot_raw__ ("set nohidden3d;\n")
%! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;

%!test
%!      % simple test
%!  x = [1,2,3];
%!  y = [4,5,6,7];
%!  [X, Y] = meshgrid(x,y);
%!  Orig = X.^2 + Y.^3;
%!  xi = [1.2,2, 1.5];
%!  yi = [6.2, 4.0, 5.0]';
%!
%!  Expected = \
%!    [243,   245.4,  243.9;
%!      65.6,  68,     66.5;
%!     126.6, 129,    127.5];
%!  Result = interp2(x,y,Orig, xi, yi);
%!
%!  assert(Result, Expected, 1000*eps);

%!test
%!      % test for values outside of boundaries
%!  x = [1,2,3];
%!  y = [4,5,6,7];
%!  [X, Y] = meshgrid(x,y);
%!  Orig = X.^2 + Y.^3;
%!  xi = [0,4];
%!  yi = [3, 8]';
%!  Expected = [nan,nan; nan,nan];
%!  Result = interp2(x,y,Orig, xi, yi);
%!  
%!  assert(Result, Expected);

%!test
%!      % test for values at boundaries
%!  A=[1,2;3,4];
%!  x=[0,1]; 
%!  y=[2,3];
%!  xi=x;
%!  yi=y;
%!  Expected = A;
%!  
%!  Result = interp2(x,y,A,xi,yi,'linear');
%!  assert(Result, Expected);
%! 
%!  Result = interp2(x,y,A,xi,yi,'nearest');
%!  assert(Result, Expected);