Mercurial > hg > octave-image

changeset 227:1a4b84759d2c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
New functions and an update of imresize to use imremap
author hauberg
date Sat, 13 Jan 2007 14:43:38 +0000
parents 015626391be7
children 5dc5f8566898
files INDEX inst/imperspectivewarp.m inst/imremap.m inst/imresize.m
diffstat 4 files changed, 445 insertions(+), 58 deletions(-) [+]
line wrap: on
line diff
--- a/INDEX
+++ b/INDEX
@@ -10,7 +10,9 @@
  imginfo
  bmpwrite jpgread jpgwrite pngread pngwrite
 Reshape
- imcrop 
+ imcrop
+ imremap
+ imperspectivewarp 
  imresize 
  imrotate
  imrotate_Fourier
new file mode 100644
--- /dev/null
+++ b/inst/imperspectivewarp.m
@@ -0,0 +1,161 @@
+## Copyright (C) 2006  Søren Hauberg
+## 
+## 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, 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 file.  If not, write to the Free Software Foundation,
+## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} @var{warped} = imperspectivewarp(@var{im}, @var{P}, @var{interp}, @var{bbox}, @var{extrapval})
+## @deftypefnx{Function File} [@var{warped}, @var{valid}] = imperspectivewarp(...)
+## Applies the spatial perspective homogeneous transformation @var{P} to the image @var{im}.
+## The transformation matrix @var{P} must be a 3x3 homogeneous matrix, or 2x2 or 2x3
+## affine transformation matrix.
+## 
+## The resulting image @var{warped} is computed using an interpolation method that
+## can be selected through the @var{interp} argument. This must be one
+## of the following strings
+## @table @code
+## @item "nearest"
+## Nearest neighbor interpolation.
+## @item "linear"
+## @itemx "bilinear"
+## Bilinear interpolation. This is the default behavior.
+## @item "cubic"
+## @itemx "bicubic"
+## Bicubic interpolation.
+## @end table
+##
+## By default the resulting image contains the entire warped image. In some situation
+## you only parts of the warped image. The argument @var{bbox} controls this, and can
+## be one of the following strings
+## @table @code
+## @item "loose"
+## The entire warped result is returned. This is the default behavior.
+## @item "crop"
+## The central part of the image of the same size as the input image is returned.
+## @end table
+##
+## All values of the result that fall outside the original image will
+## be set to @var{extrapval}. For images of class @code{double} @var{extrapval}
+## defaults to @code{NA} and for other classes it defaults to 0.
+##
+## The optional output @var{valid} is a matrix of the same size as @var{warped}
+## that contains the value 1 in pixels where @var{warped} contains an interpolated
+## value, and 0 in pixels where @var{warped} contains an extrapolated value.
+## @seealso{imremap, imrotate, imresize, imshear, interp2}
+## @end deftypefn
+
+function [warped, valid] = imperspectivewarp(im, P, interp = "bilinear", bbox = "loose", extrapolation_value = NA)
+  ## Check input
+  if (nargin < 2)
+    print_usage();
+  endif
+  
+  if (!(isgray(im) || isrgb(im)))
+    error("imperspectivewarp: first input argument must be an image");
+  endif
+  
+  if (ismatrix(P) && ndims(P) == 2)
+    if (issquare(P) && rows(P) == 3) # 3x3 matrix
+      if (P(3,3) != 0)
+        P /= P(3,3);
+      else
+        error("imperspectivewarp: P(3,3) must be non-zero");
+      endif
+    elseif (rows(P) == 2 && (columns(P) == 2 || columns(P) == 3)) # 2x2 or 2x3 matrix
+      P(3,3) = 1;
+    else # unsupported matrix size
+      error("imperspectivewarp: transformation matrix must be 2x2, 2x3, or 3x3");
+    endif
+  else
+    error("imperspectivewarp: transformation matrix not valid");
+  endif
+  
+  if (!any(strcmpi(interp, {"nearest", "linear", "bilinear", "cubic", "bicubic", "pchip"})))
+    error("imperspectivewarp: unsupported interpolation method");
+  endif
+  if (any(strcmpi(interp, {"bilinear", "bicubic"})))
+    interp = interp(3:end); # Remove "bi"
+  endif
+  interp = lower(interp);
+  
+  if (!any(strcmpi(bbox, {"loose", "crop"})))
+    error("imperspectivewarp: bounding box must be either 'loose' or 'crop'");
+  endif
+  
+  if (!isscalar(extrapolation_value))
+    error("imperspective: extrapolation value must be a scalar");
+  endif
+  
+  ## Do the transformation
+  [y, x, tmp] = size(im);
+  ## Transform corners
+  corners = [1, 1, 1;
+             1, y, 1;
+             x, 1, 1;
+             x, y, 1]';
+  Tcorners = P*corners;
+  Tx = Tcorners(1,:)./Tcorners(3,:);
+  Ty = Tcorners(2,:)./Tcorners(3,:);
+
+  ## Do cropping?
+  x1 = round(min(Tx)); x2 = round(max(Tx));
+  y1 = round(min(Ty)); y2 = round(max(Ty));
+  # FIXME: This seems to work fine for rotations, but
+  # somebody who knows computational geometry should
+  # be able to come up with a better algorithm.
+  if (strcmpi(bbox, "crop"))
+    xl = x2 - x1 + 1;
+    yl = y2 - y1 + 1;
+    xd = (xl - x)/2;
+    yd = (yl - y)/2;
+    x1 += xd; x2 -= xd;
+    y1 += yd; y2 -= yd;
+  endif
+ 
+  ## Transform coordinates
+  [X, Y] = meshgrid(x1:x2, y1:y2);
+  [sy, sx] = size(X);
+  D = [X(:), Y(:), ones(sx*sy, 1)]';
+  PD = inv(P)*D;
+  XI = PD(1,:)./PD(3,:);
+  YI = PD(2,:)./PD(3,:);
+  XI = reshape(XI, sy, sx);
+  YI = reshape(YI, sy, sx);
+  
+  ## Interpolate
+  [warped, valid] = imremap(im, XI, YI, interp, extrapolation_value);
+
+endfunction
+
+%!demo
+%! ## Generate a synthetic image and show it
+%! I = tril(ones(100)) + abs(rand(100)); I(I>1) = 1;
+%! I(20:30, 20:30) = !I(20:30, 20:30);
+%! I(70:80, 70:80) = !I(70:80, 70:80);
+%! imshow(I);
+%! ## Resize the image to the double size and show it
+%! P = diag([1, 1, 0.5]);
+%! warped = imperspectivewarp(I, P);
+%! imshow(warped);
+
+%!demo
+%! ## Generate a synthetic image and show it
+%! I = tril(ones(100)) + abs(rand(100)); I(I>1) = 1;
+%! I(20:30, 20:30) = !I(20:30, 20:30);
+%! I(70:80, 70:80) = !I(70:80, 70:80);
+%! imshow(I);
+%! ## Rotate the image around (0, 0) by -0.4 radians and show it
+%! R = [cos(-0.4) sin(-0.4); -sin(-0.4) cos(-0.4)];
+%! warped = imperspectivewarp(I, R);
+%! imshow(warped);
new file mode 100644
--- /dev/null
+++ b/inst/imremap.m
@@ -0,0 +1,223 @@
+## Copyright (C) 2006  Søren Hauberg
+## 
+## 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, 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 file.  If not, write to the Free Software Foundation,
+## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} @var{warped} = imremap(@var{im}, @var{XI}, @var{YI})
+## @deftypefnx{Function File} @var{warped} = imremap(@var{im}, @var{XI}, @var{YI}, @var{interp}, @var{extrapval})
+## @deftypefnx{Function File} [@var{warped}, @var{valid} ] = imremap(...)
+## Applies any geometric transformation to the image @var{im}.
+##
+## The arguments @var{XI} and @var{YI} are lookup tables that define the resulting
+## image
+## @example
+## @var{warped}(y,x) = @var{im}(@var{YI}(y,x), @var{XI}(y,x))
+## @end example
+## where @var{im} is assumed to be a continuous function, which is achieved
+## by interpolation. Note that the image @var{im} is expressed in a (X, Y)-coordinate
+## system and not a (row, column) system.
+##
+## The argument @var{interp} selects the used interpolation method, and most be one
+## of the following strings
+## @table @code
+## @item "nearest"
+## Nearest neighbor interpolation.
+## @item "linear"
+## @itemx "bilinear"
+## Bilinear interpolation. This is the default behavior.
+## @item "cubic"
+## @itemx "bicubic"
+## Bicubic interpolation.
+## @end table
+##
+## All values of the result that fall outside the original image will
+## be set to @var{extrapval}. For images of class @code{double} @var{extrapval}
+## defaults to @code{NA} and for other classes it defaults to 0.
+##
+## The optional output @var{valid} is a matrix of the same size as @var{warped}
+## that contains the value 1 in pixels where @var{warped} contains an interpolated
+## value, and 0 in pixels where @var{warped} contains an extrapolated value.
+## @seealso{imperspectivewarp, imrotate, imresize, imshear, interp2}
+## @end deftypefn
+
+function [warped, valid] = imremap(im, XI, YI, interp = "bilinear", extrapval = NA)
+  ## Check input
+  if (nargin < 3)
+    print_usage();
+  endif
+  
+  if (!(isgray(im) || isrgb(im)))
+    error("imremap: first input argument must be an image");
+  endif
+  
+  if (!size_equal(XI, YI) || !ismatrix(XI) || ndims(XI) != 2)
+    error("imremap: XI and YI must be matrices of the same size");
+  endif
+  
+  if (!any(strcmpi(interp, {"nearest", "linear", "bilinear", "cubic", "bicubic"})))
+    error("imremap: unsupported interpolation method");
+  endif
+  if (any(strcmpi(interp, {"bilinear", "bicubic"})))
+    interp = interp(3:end); # Remove "bi"
+  endif
+  interp = lower(interp);
+  
+  if (!isscalar(extrapval))
+    error("imremap: extrapolation value must be a scalar");
+  endif
+  
+  ## Interpolate
+  if (isgray(im))
+    warped = grayinterp(im, XI, YI, interp, NA);
+  else # rgb image
+    for i = 3:-1:1
+      warped(:,:,i) = grayinterp(im(:,:,i), XI, YI, interp, NA);
+    endfor
+  endif
+  valid = !isna(warped);
+  warped(!valid) = extrapval;
+
+  ## Change the class of the results according to the class of the image
+  c = class(im);
+  if (strcmpi(c, "uint8"))
+    warped = uint8(warped);
+  elseif (strcmpi(c, "uint16"))
+    warped = uint16(warped);
+  endif
+
+endfunction
+
+function [warped, valid] = grayinterp(im, XI, YI, interp, extrapval)
+  if (strcmp(interp, "cubic"))
+    warped = graybicubic(double(im), XI, YI, NA);
+  else
+    warped = interp2(double(im), XI, YI, interp, NA);
+  endif
+  valid = !isna(warped);
+  warped(!valid) = extrapval;
+endfunction
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{zi}=} bicubic (@var{x}, @var{y}, @var{z}, @var{xi}, @var{yi})
+## Reference:
+## Image Processing, Analysis, and Machine Vision, 2nd Ed.
+## Sonka et.al.
+## Brooks/Cole Publishing Company
+## ISBN: 0-534-95393-X
+## @seealso{interp2}
+## @end deftypefn
+
+function ZI = graybicubic (Z, XI, YI, extrapval = NA)
+  
+  ## Allocate output
+  [X, Y] = meshgrid(1:columns(Z), 1:rows(Z));
+  [Zr, Zc] = size(XI);
+  ZI = zeros(Zr, Zc);
+  
+  ## Find inliers
+  inside = !( XI < X(1) | XI > X(end) | YI < Y(1) | YI > Y(end) );
+  
+  ## Scale XI and YI to match indices of Z (not needed when interpolating images)
+  #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(inside) = AY_2 .* AX_2 .* Z(sym_sub2ind(sz, L+2, K+2)) ...
+  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;
+
+endfunction
+
+## Checks if data is meshgrided
+function b = isgriddata(X)
+  D = X - repmat(X(1,:), rows(X), 1);
+  b = all(D(:) == 0);
+endfunction
+
+## Checks if data is equally spaced (assumes data is meshgrided)
+function b = isequallyspaced(X)
+  Dx = gradient(X(1,:));
+  b = all(Dx == Dx(1));
+endfunction
+
+## Computes the 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);
+  Y(Y>sz(1)) = 2*sz(1) - Y(Y>sz(1));
+  X(X<1) = 1 - X(X<1);
+  X(X>sz(2)) = 2*sz(2) - X(X>sz(2));
+  ind = sub2ind(sz, Y, X);
+endfunction
+
+%!demo
+%! ## Generate a synthetic image and show it
+%! I = tril(ones(100)) + abs(rand(100)); I(I>1) = 1;
+%! I(20:30, 20:30) = !I(20:30, 20:30);
+%! I(70:80, 70:80) = !I(70:80, 70:80);
+%! imshow(I);
+%! ## Resize the image to the double size and show it
+%! [XI, YI] = meshgrid(linspace(1, 100, 200));
+%! warped = imremap(I, XI, YI);
+%! imshow(warped);
+
+%!demo
+%! ## Generate a synthetic image and show it
+%! I = tril(ones(100)) + abs(rand(100)); I(I>1) = 1;
+%! I(20:30, 20:30) = !I(20:30, 20:30);
+%! I(70:80, 70:80) = !I(70:80, 70:80);
+%! imshow(I);
+%! ## Rotate the image around (0, 0) by -0.4 radians and show it
+%! [XI, YI] = meshgrid(1:100);
+%! R = [cos(-0.4) sin(-0.4); -sin(-0.4) cos(-0.4)];
+%! RXY = [XI(:), YI(:)] * R;
+%! XI = reshape(RXY(:,1), [100, 100]); YI = reshape(RXY(:,2), [100, 100]);
+%! warped = imremap(I, XI, YI);
+%! imshow(warped);
--- a/inst/imresize.m
+++ b/inst/imresize.m
@@ -1,4 +1,4 @@
-## Copyright (C) 2005 Søren Hauberg
+## Copyright (C) 2005 Søren Hauberg
 ## 
 ## 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
@@ -15,47 +15,47 @@
 ## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 
 ## -*- texinfo -*-
-## @deftypefn  {Function File} @var{B}= imresize (@var{A}, @var{m})
-## @deftypefnx {Function File} @var{B}= imresize (@var{A}, @var{m}, @var{method})
+## @deftypefn {Function File} @var{B}= imresize (@var{A}, @var{m})
+## Scales the image @var{A} by a factor @var{m} using nearest neighbour
+## interpolation. If @var{m} is less than 1 the image size will be reduced,
+## and if @var{m} is greater than 1 the image will be enlarged. If the image
+## is being enlarged the it will be convolved with a 11x11 Gaussian FIR filter
+## to reduce aliasing. See below on how to alter this behavior.
+##
+## @deftypefnx {Function File} @var{B}= imresize (@var{A}, @var{m}, @var{interp})
+## Same as above except @var{interp} interpolation is performed instead of
+## using nearest neighbour. @var{interp} can be any interpolation method supported by interp2.
+##
 ## @deftypefnx {Function File} @var{B}= imresize (@var{A}, [@var{mrow} @var{mcol}])
-## @deftypefnx {Function File} @var{B}= imresize (@var{A}, [@var{mrow} @var{mcol}], @var{method})
-## @deftypefnx {Function File} @var{B}= imresize (..., @var{method}, @var{fsize})
-## @deftypefnx {Function File} @var{B}= imresize (..., @var{method}, @var{filter})
-## Resizes the image @var{A} to a given size using any interpolation
-## method supported by @code{interp2}.
-##
-## If the second argument is a scalar @var{m} the image will be scaled
-## by a factor @var{m}. If the second argument is a two-vector 
-## [@var{mrow}, @var{mcol}] the resulting image will be resized to have
-## this size.
+## Scales the image @var{A} to be of size @var{mrow}x@var{mcol} using nearest
+## neighbour interpolation. If the image is being enlarged it will be convolved
+## with a lowpass FIR filter as described above.
 ##
-## The third argument controls the kind of interpolation used to perform
-## the resizing. This can be any method supported by @code{interp2}, and
-## defaults to nearest neighbour interpolation.
+## @deftypefnx {Function File} @var{B}= imresize (@var{A}, [@var{mrow} @var{mcol}], @var{interp})
+## Same as above except @var{interp} interpolation is performed instead of using
+## nearest neighbour. @var{interp} can be any interpolation method supported by interp2.
 ##
-## The fourth argument controls if the image is filtered before resizing
-## to prevent aliasing. If the fourth argument is a two-vector 
-## [@code{frow}, @code{fcol}], the image will be convolved with a
-## @code{frow}x@code{fcol} gaussian filter. If it is a matrix it will
-## be used as a filter. The default is not to filter the image if it
-## is reduced in size, and to filter it by an 11x11 gaussian filter if
-## the image is enlarged. Setting the filter to 0 will turn of any
-## filtering.
-## @seealso{interp2}
+## @deftypefnx {Function File} @var{B}= imresize (..., @var{interp}, @var{fsize})
+## If the image the image is being reduced it will usually be convolved with
+## a 11x11 Gaussian FIR filter. By setting @var{fsize} to 0 this will be turned
+## off, and if @var{fsize} > 0 the image will be convolved with a @var{fsize}
+## by @var{fsize} Gaussian FIR filter.
+##
+## @deftypefnx {Function File} @var{B}= imresize (..., @var{interp}, @var{filter})
+## If the image size is being reduced and the @var{filter} argument is passed to
+## imresize the image will be convolved with @var{filter} before the resizing
+## takes place.
+##
+## @seealso{imremap, imrotate, interp2}
 ## @end deftypefn
 
-## Author: Søren Hauberg <hauberg at gmail dot com>
-## 
-## 2005-04-14 Søren Hauberg <hauberg at gmail dot com>
-## * Initial revision
+function ret = imresize(im, m, interp = "nearest", filter = 11)
+  if (!isgray(im) && !isrgb(im))
+    error("imresize: the first argument has must be an image");
+  endif
+  [row, col, num_planes] = size(im);
 
-function [ ret ] = imresize (im, m, method, filter)
-  if (!isgray(im))
-    error("The first argument has to be a gray-scale image.");
-  endif
-  [row col] = size(im);
-
-  # Handle the argument that describes the size of the result
+  ## Handle the argument that describes the size of the result
   if (length(m) == 1)
     new_row = round(row*m);
     new_col = round(col*m);
@@ -64,36 +64,37 @@
     new_col = m(2);
     m = min( new_row/row, new_col/col );
   else
-    error("Bad second argument");
+    error("imresize: second argument mest be a scalar or a 2-vector");
   end
 
-  # Handle the method argument
-  if (nargin < 3)
-    method = "nearest";
+  ## Handle the method argument
+  if (!any(strcmpi(interp, {"nearest", "linear", "bilinear", "cubic", "bicubic", "pchip"})))
+    error("imresize: unsupported interpolation method");
   endif
 
-  # Handle the filterargument
-  if (!strcmp(method, "nearest"))
-    if (nargin < 4)
-      filter = 11;
-    endif
-    if (m > 1 & filter > 0)
-      # If the image is being enlarged and filter > 0 then
-      # convolve the image with a filter*filter gaussian.
-      mu = round(filter/2);
+  ## Handle the filter argument
+  if (!strcmp(interp, "nearest") && m < 1)
+    if (isscalar(filter) && filter > 0)
+      ## If the image is being reduced and filter > 0 then
+      ## convolve the image with a filter*filter gaussian.
+      mu = (filter/2);
       sigma = mu/3;
       x = 1:filter;
       gauss = 1/sqrt(2*pi*sigma^2) * exp( (-(x-mu).^2)/(2*sigma^2) );
-      im = conv2(gauss, gauss, im, "same");
-    elseif (m < 1 & nargin == 4)
-      # If the image size is being reduced and a fourth argument
-      # is given, use it as a FIR filter.
-      im = conv2(im, filter, "same");
+      for i = 1:num_planes
+          im(:,:,i) = conv2(gauss, gauss, im(:,:,i), "same");
+      endfor
+    elseif (nargin >= 4 && ismatrix(filter) && ndims(filter) == 2)
+      ## If the image size is being reduced and a fourth argument
+      ## is given, use it as a FIR filter.
+      for i = 1:num_planes
+        im(:,:,i) = conv2(im(:,:,i), filter, "same");
+      endfor
     endif
   endif
   
-  # Perform the actual resizing
-  [XI YI] = meshgrid( linspace(1,col,new_col), linspace(1,row,new_row) );
-  ret = interp2(im, XI, YI, method);
+  ## Perform the actual resizing
+  [XI, YI] = meshgrid( linspace(1,col,new_col), linspace(1,row,new_row) );
+  ret = imremap(im, XI, YI, interp);
 
 endfunction