--- a/scripts/ChangeLog	Thu Oct 22 10:12:54 2009 +0200
+++ b/scripts/ChangeLog	Tue Oct 20 20:22:10 2009 +0200
@@ -1,3 +1,8 @@
+2009-10-20  Soren Hauberg  <hauberg@gmail.com>
+
+	* general/interp2.m: improved error checking and support for bicubic
+	interpolation when X and Y are meshgrid format.
+
 2009-10-22  Jaroslav Hajek  <highegg@gmail.com>
 
 	* general/interp1.m: Perform optimizations, improve code (use switch

--- a/scripts/general/interp2.m	Thu Oct 22 10:12:54 2009 +0200
+++ b/scripts/general/interp2.m	Tue Oct 20 20:22:10 2009 +0200
@@ -57,7 +57,7 @@
 ## @item 'linear'
 ## Linear interpolation from nearest neighbors.
 ## @item 'pchip'
-## Piece-wise cubic hermite interpolating polynomial (not implemented yet).
+## Piece-wise cubic hermite interpolating polynomial.
 ## @item 'cubic'
 ## Cubic interpolation from four nearest neighbors.
 ## @item 'spline'
@@ -344,11 +344,64 @@
 
     ## FIXME bicubic/__splinen__ don't handle arbitrary XI, YI
     if (strcmp (method, "cubic"))
-      ZI = bicubic (X, Y, Z, XI(1,:), YI(:,1), extrapval);
+      if (isgriddata (XI) && isgriddata (YI'))
+        ZI = bicubic (X, Y, Z, XI (1, :), YI (:, 1), extrapval);
+      elseif (isgriddata (X) && isgriddata (Y'))
+        ## Allocate output
+        ZI = zeros (size (X));
+  
+        ## Find inliers
+        inside = !(XI < X (1) | XI > X (end) | YI < Y (1) | YI > Y (end));
+  
+        ## Scale XI and YI to match indices of Z
+        XI = (columns (Z) - 1) * (XI - X (1)) / (X (end) - X (1)) + 1;
+        YI = (rows (Z) - 1) * (YI - Y (1)) / (Y (end) - Y (1)) + 1;
+  
+        ## Start the real work
+        K = floor (XI);
+        L = floor (YI);
+
+        ## Coefficients
+        AY1  = bc ((YI - L + 1));
+        AX1  = bc ((XI - K + 1));
+        AY0  = bc ((YI - L + 0));
+        AX0  = bc ((XI - K + 0));
+        AY_1 = bc ((YI - L - 1));
+        AX_1 = bc ((XI - K - 1));
+        AY_2 = bc ((YI - L - 2));
+        AX_2 = bc ((XI - K - 2));
+
+        ## Perform interpolation
+        sz = size(Z);
+        ZI = AY_2 .* AX_2 .* Z (sym_sub2ind (sz, L+2, K+2)) ...
+           + AY_2 .* AX_1 .* Z (sym_sub2ind (sz, L+2, K+1)) ...
+           + AY_2 .* AX0  .* Z (sym_sub2ind (sz, L+2, K))   ...
+           + AY_2 .* AX1  .* Z (sym_sub2ind (sz, L+2, K-1)) ...
+           + AY_1 .* AX_2 .* Z (sym_sub2ind (sz, L+1, K+2)) ...
+           + AY_1 .* AX_1 .* Z (sym_sub2ind (sz, L+1, K+1)) ...
+           + AY_1 .* AX0  .* Z (sym_sub2ind (sz, L+1, K))   ...
+           + AY_1 .* AX1  .* Z (sym_sub2ind (sz, L+1, K-1)) ...
+           + AY0  .* AX_2 .* Z (sym_sub2ind (sz, L,   K+2)) ...
+           + AY0  .* AX_1 .* Z (sym_sub2ind (sz, L,   K+1)) ...
+           + AY0  .* AX0  .* Z (sym_sub2ind (sz, L,   K))   ...
+           + AY0  .* AX1  .* Z (sym_sub2ind (sz, L,   K-1)) ...
+           + AY1  .* AX_2 .* Z (sym_sub2ind (sz, L-1, K+2)) ...
+           + AY1  .* AX_1 .* Z (sym_sub2ind (sz, L-1, K+1)) ...
+           + AY1  .* AX0  .* Z (sym_sub2ind (sz, L-1, K))   ...
+           + AY1  .* AX1  .* Z (sym_sub2ind (sz, L-1, K-1));
+        ZI (!inside) = extrapval;
+      
+      else
+        error ("interp2: input data must have `meshgrid' format");
+      endif
 
     elseif (strcmp (method, "spline"))
-      ZI = __splinen__ ({Y(:,1).', X(1,:)}, Z, {YI(:,1), XI(1,:)}, extrapval, 
+      if (isgriddata (XI) && isgriddata (YI'))
+        ZI = __splinen__ ({Y(:,1).', X(1,:)}, Z, {YI(:,1), XI(1,:)}, extrapval, 
 			"spline");
+      else
+        error ("interp2: input data must have `meshgrid' format");
+      endif
     else
       error ("interpolation method not recognized");
     endif
@@ -356,6 +409,38 @@
   endif
 endfunction
 
+function b = isgriddata (X)
+  d1 = diff (X, 1, 1);
+  d2 = diff (X, 1, 2);
+  b = all (d1 (:) == 0) & all (d2 (:) == d2 (1));
+endfunction
+
+## Compute the bicubic interpolation coefficients
+function o = bc(x)
+  x = abs(x);
+  o = zeros(size(x));
+  idx1 = (x < 1);
+  idx2 = !idx1 & (x < 2);
+  o(idx1) = 1 - 2.*x(idx1).^2 + x(idx1).^3;
+  o(idx2) = 4 - 8.*x(idx2) + 5.*x(idx2).^2 - x(idx2).^3;
+endfunction
+
+## This version of sub2ind behaves as if the data was symmetrically padded
+function ind = sym_sub2ind(sz, Y, X)
+  Y (Y < 1) = 1 - Y (Y < 1);
+  while (any (Y (:) > 2 * sz (1)))
+    Y (Y > 2 * sz (1)) = round (Y (Y > 2 * sz (1)) / 2);
+  endwhile
+  Y (Y > sz (1)) = 1 + 2 * sz (1) - Y (Y > sz (1));
+  X (X < 1) = 1 - X (X < 1);
+  while (any (X (:) > 2 * sz (2)))
+    X (X > 2 * sz (2)) = round (X (X > 2 * sz (2)) / 2);
+  endwhile
+  X (X > sz (2)) = 1 + 2 * sz (2) - X (X > sz (2));
+  ind = sub2ind(sz, Y, X);
+endfunction
+
+
 %!demo
 %! A=[13,-1,12;5,4,3;1,6,2];
 %! x=[0,1,4]; y=[10,11,12];