# HG changeset patch
# User Christof Kaufmann <christofkaufmann.dev@gmail.com>
# Date 1643304199 -3600
#      Thu Jan 27 18:23:19 2022 +0100
# Node ID e238a14a9148d48213cae1e6b07f2f70dc50229f
# Parent  232da3e2938032c6e454a5ad611785b799be8209
imresize: Implement anti-aliasing (bug #61760)

diff --git a/inst/imresize.m b/inst/imresize.m
--- a/inst/imresize.m
+++ b/inst/imresize.m
@@ -18,6 +18,7 @@
 ## @deftypefn  {Function File} {} imresize (@var{im}, @var{scale})
 ## @deftypefnx {Function File} {} imresize (@var{im}, [@var{M} @var{N}])
 ## @deftypefnx {Function File} {} imresize (@dots{}, @var{method})
+## @deftypefnx {Function File} {} imresize (@dots{}, @dots{}, @var{property}, @var{value}, @dots{})
 ## Resize image with interpolation
 ##
 ## Scales the image @var{im} by a factor @var{scale} or into the size @var{M}
@@ -37,53 +38,175 @@
 ## as to preserve aspect ratio.
 ##
 ## The optional argument @var{method} defines the interpolation method to be
-## used.  All methods supported by @code{interp2} can be used.  By default, the
-## @code{cubic} method is used.
+## used.  The following methods are available, see below for custom methods.
+## @table @code
+## @item "nearest", "box"
+## Nearest neighbor method.  Only the nearest pixel is used.  This gives hard
+## edges.
 ##
-## For @sc{matlab} compatibility, the methods @code{bicubic} (same as
-## @code{cubic}), @code{bilinear} and @code{triangle} (both the same as
-## @code{linear}) are also supported.
+## @item "linear", "bilinear", "triangle"
+## Bilinear interpolation method using the four neighbor pixels.
+##
+## @item "cubic", "bicubic" (default)
+## Bicubic interpolation method using the 16 neighbor pixels.  At the borders
+## symmetric padding is used.  This is the default method.
+## @end table
+## By default, the @code{cubic} method is used.
 ##
-## @table @asis
-## @item bicubic (default)
-## Same as @code{cubic}.
+## For custom interpolation kernels, specify a two-element cell array for
+## @var{method}: {@var{kernel}, @var{size}}.  @var{kernel} must be an
+## interpolation kernel function that can handle vector input.  It must be zero
+## outside -@var{size}/2 <= x <= @var{size}/2.  The following example does a
+## bilinear interpolation, just as using "bilinear":
+## @example
+## @group
+## im = magic(6);
+## lin = @(x) (1 - abs(x)) .* (abs(x) < 1);
+## size = 2;
+## imresize(im, 0.5, {lin, size});
+## @end group
+## @end example
+## Note that also for custom kernels anti-aliasing is applied by default.
 ##
-## @item bilinear
-## Same as @code{linear}.
+## Additionally the following optional property-value-pairs can be used:
+## @table @code
+## @item "Antialiasing"
+## If this is set to @code{true} and the scale factor in horizontal or vertical
+## direction is less than 1, anti-aliasing will used for that direction.  This
+## means the interpolation kernel is broadened by 1/scale to reduce the
+## frequency components that cause aliasing effects.  Hence more neighbors than
+## described above are used, e. g. "bilinear" with a scale of 0.5 in both
+## directions uses 16 neighbors.  The default value is @code{true}, except for
+## the method "nearest" / "box".
 ##
-## @item triangle
-## Same as @code{linear}.
+## @item "Method"
+## The interpolation method as string or a custom interpolation kernel, see
+## above.
+##
+## @item "OutputSize"
+## Specify the output size @var{M} rows by @var{N} columns as vector
+## [@var{M} @var{N}], see above.
+##
+## @item "Scale"
+## Either a scalar to use the same scale factor for both direction or a vector
+## [@var{scale_rows} @var{scale_columns}] to use different scaling factors.
 ## @end table
 ##
-## Note: Antialiasing is not implemented yet.
+## Note: Currently there is no special support for categorical images or images
+## with indexed colors.
 ##
 ## @seealso{imremap, imrotate, interp2}
 ## @end deftypefn
 
-function im = imresize (im, scale, method = "cubic")
+function im = imresize (im, varargin)
+
+  if (nargin < 2 || nargin > 7)
+    print_usage ();
+  endif
+
+  antialiasing = [];
+  scale_or_M_N = [];
+  scale = []; # for the property, that can specify [scale_rows, scale_cols]
+  method = [];
+
+  if (nargin/2 == round (nargin/2))  # even number of inputs (2, 4 or 6)
+    ## imresize (im, scale_or_M_N)
+    ## imresize (im, scale_or_M_N, property1, value1)
+    ## imresize (im, scale_or_M_N, property1, value1, property2, value2)
+    scale_or_M_N = varargin{1};
+    n_start = 2;
+  else                              # odd number of inputs (3, 5 or 7)
+    ## imresize (im, property1, value1) # property1 must be Scale or OutputSize
+    ## imresize (im, scale_or_M_N, method)
+    ## imresize (im, scale_or_M_N, method, property1, value1)
+    ## imresize (im, scale_or_M_N, method, property1, value1, property2, value2)
+    if isnumeric(varargin{1})
+      scale_or_M_N = varargin{1};
+      method = varargin{2};
+      n_start = 3;
+    else
+      n_start = 1;
+    endif
+  endif
 
-  if (nargin < 2 || nargin > 3)
-    print_usage ();
+  for n = n_start:2:(nargin-1)      # process parameter-values pairs
+    if (strcmpi (varargin{n}, "Antialiasing"))
+      antialiasing = varargin{n+1};
+    elseif (strcmpi (varargin{n}, "Method"))
+      method = varargin{n+1};
+    elseif (strcmpi (varargin{n}, "OutputSize"))
+      if (! isempty (scale_or_M_N))
+        error ("imresize: OutputSize must not be specified, when SCALE or [M N] is also specified.")
+      endif
+      scale_or_M_N = varargin{n+1};
+    elseif (strcmpi (varargin{n}, "Scale"))
+      scale = varargin{n+1};
+    else
+      error ("imresize: invalid PROPERTY given")
+    endif
+  endfor
+
+  ## defaults
+  if (isempty (method))
+    method = "cubic";
+  elseif (ischar (method))
+    ## convert to lower case. Replace "box" by "nearest", replace "bicubic" by "cubic" and replace "bilinear" and "triangle" by "linear".
+    method = interp_method (method);
+  endif
+
+  if (isempty (antialiasing))
+    if (strcmpi (method, "nearest"))
+      antialiasing = false;
+    else
+      antialiasing = true;
+    endif
+  endif
+
+  ## check input arguments
+  if (isempty (scale_or_M_N) && isempty (scale))
+    error ("imresize: Scale or output size must be specified.");
+  elseif (! isempty(scale_or_M_N) && ! isempty (scale))
+    error ("imresize: Scale and OutputSize must not be specified both.")
+  elseif (! (islogical (antialiasing) || isnumeric (antialiasing)) || ! isscalar (antialiasing))
+    ## accept also numbers as logical values (even complex, which Matlab does not accept)
+    error ("imresize: Antialiasing must be true, false or a number (0 interpreted as false, otherwise as true).")
   elseif (! ((isnumeric (im) || islogical (im)) && ! issparse (im) && ! isempty (im)))
     error ("imresize: IM must be an image")
-  elseif (! isnumeric (scale) || any (scale <= 0) || all (isnan (scale)))
+  elseif (! isempty (scale_or_M_N) && (! isnumeric (scale_or_M_N) || any (scale_or_M_N <= 0)))
     error ("imresize: SCALE or [M N] must be numeric positive values")
-  elseif (! ischar (method))
-    error ("imresize: METHOD must be a string with interpolation method")
+  elseif (! isempty (scale_or_M_N) && numel (scale_or_M_N) > 2)
+    error ("imresize: SCALE or [M N] argument must be a scalar or a 2 element vector");
+  elseif (! isempty (scale_or_M_N) && all (isnan (scale_or_M_N)))
+    error ("imresize: In [M N] only one value may be NaN to maintain aspect ratio.")
+  elseif (! isempty (scale) && (! isnumeric (scale) || any (scale <= 0) || any (isnan (scale)) || numel (scale) > 2))
+    error ("imresize: Scale must be one or two numeric positive values")
+  elseif (! ischar (method) && ! (iscell (method) && length (method) == 2))
+    error ("imresize: METHOD must be a string with the interpolation method or a two-element cell array with a custom kernel and size.")
   endif
-  method = interp_method (method);
+
 
   in_rows = rows (im);
   in_cols = columns (im);
 
-  if (isscalar (scale))
-    out_rows = ceil (in_rows * scale);
-    out_cols = ceil (in_cols * scale);
-    scale_rows = scale;
-    scale_cols = scale;
-  elseif (numel (scale) == 2)
-    out_rows = scale(1);
-    out_cols = scale(2);
+  if (isscalar (scale_or_M_N))
+    scale = scale_or_M_N;
+    scale_or_M_N = [];
+  endif
+
+  if (isempty (scale_or_M_N))
+    if (isscalar (scale))
+      scale_rows = scale;
+      scale_cols = scale;
+    else
+      scale_rows = scale(1);
+      scale_cols = scale(2);
+    endif
+    out_rows = ceil (in_rows * scale_rows);
+    out_cols = ceil (in_cols * scale_cols);
+  else
+    # scale_or_M_N contains output size
+    out_rows = scale_or_M_N(1);
+    out_cols = scale_or_M_N(2);
 
     ## maintain aspect ratio if requested
     if (isnan (out_rows))
@@ -96,9 +219,7 @@ function im = imresize (im, scale, metho
     scale_cols = out_cols / in_cols;
     out_rows = ceil (out_rows);
     out_cols = ceil (out_cols);
-  else
-    error ("imresize: SCALE argument must be a scalar or a 2 element vector");
-  end
+  endif
 
   ## calculate the new pixel indices in terms of the old pixel indices
   off_rows = 1 / scale_rows / 2;
@@ -106,28 +227,164 @@ function im = imresize (im, scale, metho
   idx_rows = 0.5 + off_rows + (0:out_rows-1) / scale_rows;
   idx_cols = 0.5 + off_cols + (0:out_cols-1) / scale_cols;
 
-  ## Simple replicate padding for nearest neighbor and bilinear interpolation. Padding for bicubic interpolation is handled below.
-  ## Limits indices to the bounds 1 and end.
-  if (strcmp(method, "linear"))
+  ## trivial cases
+  if (scale_rows == 1 && scale_cols == 1)
+    ## no resizing to do
+    return
+  elseif (ischar (method) && strcmp (method, "nearest") ...
+          && (! antialiasing || (scale_rows >= 1 && scale_cols >= 1)))
     idx_rows = max (min (idx_rows, in_rows), 1);
     idx_cols = max (min (idx_cols, in_cols), 1);
+    im = im(round (idx_rows), round (idx_cols), :);
+    return
   endif
 
-  ## perform the actual resizing
-  if (in_rows == out_rows && in_cols == out_cols)
-    ## no resizing to do
-  elseif (strcmp (method, "nearest"))
-    im = im(round (idx_rows), round (idx_cols), :);
-  elseif (strcmp (method, "cubic"))
-    ## When rounding the output size up, it can happen that some interpolation points lap out at the right or the bottom of the input image.
-    ## In these cases we do not want to set the outlying points to 0 afterwards, but use the padded values. "symmetric" is also used in Matlab.
+  ## cast to floating point for accuracy
+  orig_class = class (im);
+  switch orig_class
+    case {"int8", "uint8", "int16", "uint16"}
+      im = single (im);
+      idx_cols = single (idx_cols);
+      idx_rows = single (idx_rows);
+    case {"int32", "uint32", "int64", "uint64"}
+      im = double (im);
+  endswitch
+
+  ## actual interpolation
+  for scale_and_idx = {scale_cols, scale_rows;
+                       idx_cols,   idx_rows';
+                       2,          1}
+    [scale, idx, axis] = scale_and_idx{:};
+    if scale == 1
+      continue;
+    endif
+
+    if (iscell (method))
+      kernel_size = method{2};
+      if (scale < 1 && antialiasing)
+        kernel = @(h) scale * method{1} (scale * h);
+        kernel_size /= scale;
+      else
+        kernel = method{1};
+      endif
+    elseif (strcmp (method, "nearest"))
+      if (scale < 1)
+        kernel = @(h) scale * box (scale * h);
+        kernel_size = 1 / scale;
+      else
+        if axis == 1
+          im = im(round (idx), :, :);
+        else
+          im = im(:, round (idx), :);
+        endif
+        continue;
+      endif
+    elseif (strcmp (method, "linear"))
+      kernel_size = 2;
+      if (scale < 1 && antialiasing)
+        kernel = @(h) scale * triangle (scale * h);
+        kernel_size /= scale;
+      else
+        kernel = @triangle;
+      endif
+    elseif (strcmp (method, "cubic"))
+      kernel_size = 4;
+      if (scale < 1 && antialiasing)
+        kernel = @(h) scale * cubic (scale * h);
+        kernel_size /= scale;
+      else
+        kernel = @cubic;
+      endif
+    else
+      error ("imresize: Interpolation method not supported");
+    endif
+
+    ## When rounding the output size up, it can happen that some interpolation
+    ## points lie out at the right or the bottom of the input image. For these
+    ## points symmetric padding is very important and also used in Matlab.
     ## For an example, see tests below with scale 1/3.
-    im = imremap (im, idx_cols, idx_rows, method, "symmetric");
+    im = conv_interp_vec (im, idx, kernel, kernel_size, axis, "symmetric");
+  endfor
+
+  ## we return image on same class as input
+  im = cast (im, orig_class);
+endfunction
+
+
+## box / nearest neighbor interpolation kernel.
+function w = box (d)
+  w = -0.5 < d & d <= 0.5;
+endfunction
+
+## linear interpolation kernel.
+function w = triangle (d)
+  absd = abs (d);
+  absd01 = absd <= 1;
+  w = (1 - absd) .* absd01;  ## for |d| <= 1
+endfunction
+
+## cubic interpolation kernel with a = -0.5 for MATLAB compatibility.
+function w = cubic (h)
+  absh = abs (h);
+  absh01 = absh <= 1;
+  absh12 = absh <= 2 & ~absh01;
+  absh_sqr = absh .* absh;
+  absh_cube = absh_sqr .* absh;
+  w = (1.5 * absh_cube - 2.5 * absh_sqr + 1)             .* absh01 ...  ## for |h| <= 1
+    + (-0.5 * absh_cube + 2.5 * absh_sqr - 4 * absh + 2) .* absh12;     ## for 1 < |h| <= 2
+end
+
+
+## padding by changing indices. Cannot mimic constant value padding, like zero padding
+function idx = pad_indices (i, sz, method = "symmetric")
+  if strcmp (method, "replicate")
+    idx = max (min (i, sz), 1);
+  elseif strcmp (method, "symmetric")
+    idx = i - 1;
+    m = mod (idx, sz);
+    odd = mod (floor (idx / sz), 2) == 1;
+    idx(odd) = sz - m(odd);
+    idx(!odd) = m(!odd) + 1;
+  elseif strcmp (method, "reflect")
+    idx = i - 1;
+    while (any (idx(:) < 0 | idx(:) >= sz))
+      idx(idx < 0) = -idx(idx < 0);
+      idx(idx >= sz) = 2*sz - 2 - idx(idx >= sz);
+    endwhile
+    idx += 1;
+  elseif strcmp (method, "circular")
+    idx = mod (i - 1, sz) + 1;
   else
-    im = imremap (im, idx_cols, idx_rows, method);
+    error (['imresize: Invalid argument for PADDING. Valid are "replicate", "symmetric", "reflect", "circular". You gave "', method, '"'])
   endif
 endfunction
 
+## interpolation using convolution kernel
+function out = conv_interp_vec (img, ZI, kernel, kernel_size, axis, padding = "symmetric")
+  ## get indexes and distances
+  idx = floor (ZI);
+  DZ = ZI - idx;
+  pad_size = ceil (kernel_size / 2);
+  pad_border = size (img, axis);
+
+  ## allocate output
+  out_shape = size (img);
+  out_shape(axis) = length (ZI);
+  out = zeros (out_shape);
+
+  ## interpolate
+  for shift = 1-pad_size : pad_size
+    h = shift - DZ;
+    idx_padded = pad_indices (idx + shift, pad_border, padding);
+
+    if axis == 1
+      out += img(idx_padded, :, :) .* kernel (h);
+    else
+      out += img(:, idx_padded, :) .* kernel (h);
+    endif
+  endfor
+endfunction
+
 ## Test basic features.
 %!test
 %!
@@ -190,14 +447,17 @@ endfunction
 %! assert (imresize (uint8 (in), [14 10], "nearest"), uint8 (out))
 %!
 ## Test equivalence of different input writing styles:
+%! assert (imresize (uint8 (in), 1.5, "box"), imresize (uint8 (in), 1.5, "MeTHoD", "nearest"))
+%! assert (imresize (uint8 (in), "Scale", 1.5, "Method", "box"), imresize (uint8 (in), 1.5, {@(h) -0.5 < h & h <= 0.5, 1}))
 %! assert (imresize (uint8 (in), 1.5, "bicubic"), imresize (uint8 (in), 1.5, "cubic"))
 %! assert (imresize (uint8 (in), [NaN, size(in,2)*1.5], "bicubic"), imresize (uint8 (in), 1.5, "cubic"))
 %! assert (imresize (uint8 (in), [size(in,1)*1.5, NaN], "bicubic"), imresize (uint8 (in), 1.5, "cubic"))
+%! assert (imresize (uint8 (in), "outputsize", [size(in,1)*1.5, NaN], "method", "bicubic"), imresize (uint8 (in), 1.5, "cubic"))
 %! assert (imresize (uint8 (in), 1.5, "linear"), imresize (uint8 (in), 1.5, "LIneAR"))
 %! assert (imresize (uint8 (in), 1.5, "linear"), imresize (uint8 (in), 1.5, "triangle"))
-%!
-## nearest neighbour test, check that a pixel from the neighbour locations
-## gets picked. The distance is the same for all neighbours here.
+
+## nearest neighbour test. The distance is the same for all neighbours here.
+%!test
 %! in = [116  227  153   69  146  194   59  130  139  106
 %!         2   47  137  249   90   75   16   24  158   44
 %!       155   68   46   84  166  156   69  204   32  152
@@ -207,14 +467,48 @@ endfunction
 %!         9   47  133  135  130  207  236   43   19   73
 %!       129   60   59  243   64  181  249   56   32   86];
 %!
-%! out = imresize (in, 0.5, "nearest");
+## Check that a pixel from the neighbour locations gets picked.
+%! out = imresize (in, 0.5, "nearest", "Antialiasing", false);
+%! for x = 1:columns (out)
+%!   for y = 1:rows (out)
+%!     x_in = 2 * (x-1) + 1;
+%!     y_in = 2 * (y-1) + 1;
+%!     sub = in(y_in:y_in+1, x_in:x_in+1);
+%!     assert (any (any (sub == out(y, x))))
+%!   endfor
+%! endfor
+%!
+## Check that with anti-aliasing the mean of the neighbour pixels is used.
+%! out = imresize (in, 0.5, "nearest", "Antialiasing", true);
 %! for x = 1:columns (out)
-%!     for y = 1:rows (out)
-%!         x_in = 2 * (x-1) + 1;
-%!         y_in = 2 * (y-1) + 1;
-%!         sub = in(y_in:y_in+1, x_in:x_in+1);
-%!         assert (any (any (sub == out(y,x))))
-%!     endfor
+%!   for y = 1:rows (out)
+%!     x_in = 2 * (x-1) + 1;
+%!     y_in = 2 * (y-1) + 1;
+%!     val = mean (mean (in(y_in:y_in+1, x_in:x_in+1)));
+%!     assert (val, out(y, x))
+%!   endfor
+%! endfor
+%!
+## Check that anti-aliasing also works in only y direction.
+%! out = imresize (in,  "Scale", [0.5, 2],  "Method", "nearest",  "Antialiasing", true);
+%! for x = 1:columns (out)
+%!   for y = 1:rows (out)
+%!     x_in = floor (0.5 * (x-1) + 1);
+%!     y_in = 2 * (y-1) + 1;
+%!     val = mean (in(y_in:y_in+1, x_in));
+%!     assert (val, out(y, x))
+%!   endfor
+%! endfor
+%!
+## Check that anti-aliasing also works in only x direction.
+%! out = imresize (in,  "Scale", [2, 0.5],  "Method", "nearest",  "Antialiasing", true);
+%! for x = 1:columns (out)
+%!   for y = 1:rows (out)
+%!     x_in = 2 * (x-1) + 1;
+%!     y_in = floor (0.5 * (y-1) + 1);
+%!     val = mean (in(y_in, x_in:x_in+1));
+%!     assert (val, out(y, x))
+%!   endfor
 %! endfor
 
 
@@ -232,21 +526,49 @@ endfunction
 %!test
 %!
 %! for channels = 1:3
-%!     in = rand (5, 4, channels);
-%!     for method = {"nearest", "bilinear", "bicubic"}
-%!         method = method{1};
-%!         out = imresize (in, 2, method);
-%!         for i = 1:size (in, 3)
-%!             assert (out(:,:,i), imresize (in(:,:,i), 2, method))
-%!         endfor
+%!   in = rand (5, 4, channels);
+%!   for method = {"nearest", "bilinear", "bicubic"}
+%!     out = imresize (in, 2, method{1});
+%!     for i = 1:size (in, 3)
+%!       assert (out(:, :, i), imresize (in(:, :, i), 2, method{1}))
 %!     endfor
+%!   endfor
+%! endfor
+
+## Test scaling down to a single row
+%!test
+%!
+%! for channels = 1:3
+%!   in = rand (5, 4, channels);
+%!   out = imresize (in, [1, columns(in)], "nearest", "Antialiasing", true);
+%!   for i = 1:columns (in)
+%!     assert (out(1, i, :), mean (in(:, i, :), 1), 10*eps)
+%!   endfor
 %! endfor
 
-## Test linear interpolations against some reference results from matlab
-## without anti-aliasing with
-## imresize (uint8 (in), <FACTOR>, 'bilinear', 'Antialiasing', false)
-## The floating point error is less than 1e-13, but it affects the rounding
-## to int. So a difference of 1 is acceptable.
+## Test scaling down to a single column
+%!test
+%!
+%! for channels = 1:3
+%!   in = rand (5, 4, channels);
+%!   out = imresize (in, [rows(in), 1], "nearest", "Antialiasing", true);
+%!   for i = 1:rows (in)
+%!     assert (out(i, 1, :), mean (in(i, :, :), 2), 10*eps)
+%!   endfor
+%! endfor
+
+## Test scaling down to a single pixel
+%!test
+%!
+%! for channels = 1:3
+%!   in = rand (5, 4, channels);
+%!   out = imresize (in, [1, 1], "nearest", "Antialiasing", true);
+%!   assert (out(1, 1, :), mean (mean (in(:, :, :))), 10*eps)
+%! endfor
+
+## Test linear interpolations against some reference results from matlab.
+## The floating point error is less than 1e-13, but for int matlab uses an
+## optimized algorithm. So a difference of 1 is acceptable.
 %!test
 %!
 %! in = [116  227  153   69  146  194   59  130  139  106
@@ -258,6 +580,17 @@ endfunction
 %!         9   47  133  135  130  207  236   43   19   73
 %!       129   60   59  243   64  181  249   56   32   86];
 %!
+## Factor 0.91 yields same output size, but interpolation must not be skipped
+%! out = [115  208  134  100  163  117  101  136  109  103
+%!         26   61  149  182   95   53   41  116   73   60
+%!        133  101   82  140  167  125  152   71  126  144
+%!         88  184  137  164  142  145  110   81  104  108
+%!         86  146  109   55   73  110  111  156  153  150
+%!         33  104  131  100  130  184  147   97  133  142
+%!         84   59  114  164  133  219  120   33   72   81
+%!        126   60  104  181  116  218  125   38   77   86];
+%! assert (imresize (uint8 (in), 0.91, "bilinear", "Antialiasing", false), uint8 (out), 1)
+%!
 ## Factor 1.5, gives an output of size 12 x 15 (without requiring to round the size)
 %! out = [116  172  215  165  111   82  133  170  171   81   95  132  138  123  106
 %!         59   98  138  144  152  152  125  127  119   54   58   89  137  112   75
@@ -274,11 +607,17 @@ endfunction
 %! assert (imresize (uint8 (in), 1.5, "bilinear"), uint8 (out), 1)
 %!
 ## Factor 0.5, gives an output of size 4 x 5 (without requiring to round the size)
-%! out = [ 98  152  126   57  112
-%!        129  124  171  145   83
-%!        116   95   80  116  164
-%!         61  143  146  146   53];
-%! assert (imresize (uint8 (in), 0.5, "bilinear"), uint8 (out), 1)
+%! out = [ 98  152  126   58  112
+%!        129  125  172  146   83
+%!        116   96   80  116  164
+%!         62  143  146  147   53];
+%! assert (imresize (uint8 (in), 0.5, "bilinear", "Antialiasing", false), uint8 (out), 1)
+%!
+%! out = [108  136  125   89  107
+%!        111  132  143  114   99
+%!        106  110  106  127  136
+%!         75  124  154  142   75];
+%! assert (imresize (uint8 (in), 0.5, "bilinear", "Antialiasing", true), uint8 (out), 1)
 %!
 ## Factor 4/3, gives an output of size 10.6667 x 13.3333 rounded up to 11 x 14
 %! out = [116  185  199  143   80  117  164  177   76  103  133  135  110  106
@@ -294,15 +633,23 @@ endfunction
 %!        129   86   60   82  220  131  108  190  241  128   47   39   79   86];
 %! assert (imresize (uint8 (in), 4/3, "bilinear"), uint8 (out), 1)
 %!
+## Define custom bilinear interpolation kernel
+%! lin = @(x) (1 - abs(x)) .* (abs(x) < 1);
+%!
 ## Factor 1/3, gives an output of size 2.6667 x 3.3333 rounded up to 3 x 4
 %! out = [ 47   90   24   44
 %!        143   51   21   90
 %!         60   64   56   86];
-%! assert (imresize (uint8 (in), 1/3, "bilinear"), uint8 (out), 1)
+%! assert (imresize (uint8 (in), 1/3, "bilinear", "Antialiasing", false), uint8 (out), 1)
+%! assert (imresize (uint8 (in), 1/3, {lin, 2}, "Antialiasing", false), uint8 (out), 1)
+%!
+%! out = [115  131  101  102
+%!        114  117  120  121
+%!         91  147  116   76];
+%! assert (imresize (uint8 (in), 1/3, "bilinear", "Antialiasing", true), uint8 (out), 1)
+%! assert (imresize (uint8 (in), 1/3, {lin, 2}, "Antialiasing", true), uint8 (out), 1)
 
-## Test bicubic interpolations against some reference results from matlab
-## without anti-aliasing with
-## imresize (uint8 (in), <FACTOR>, 'bicubic', 'Antialiasing', false)
+## Test bicubic interpolations against some reference results from matlab.
 ## The floating point error is less than 4e-13, but with integer Matlab rounds
 ## wrong. So a difference of 1 is acceptable.
 ## Factor 1.5, gives an output of size 12 x 15 (without requiring to round the size)
@@ -331,11 +678,17 @@ endfunction
 %! assert (imresize (uint8 (in), 1.5, "bicubic"), uint8 (out), 1)
 %!
 ## Factor 0.5, gives an output of size 4 x 5 (without requiring to round the size)
-%! out = [ 92   164   123    38   118
-%!        139   116   188   167    69
-%!        121    87    67   108   180
-%!         54   153   141   149    42];
-%! assert (imresize (uint8 (in), 0.5, "bicubic"), uint8 (out), 1)
+%! out = [ 92  164  123   38  118
+%!        139  116  188  167   69
+%!        121   87   67  108  180
+%!         54  153  141  149   42];
+%! assert (imresize (uint8 (in), 0.5, "bicubic", "Antialiasing", false), uint8 (out), 1)
+%!
+%! out = [105  140  126   81  109
+%!        110  134  153  114   93
+%!        108  108   94  127  146
+%!         67  126  162  149   62];
+%! assert (imresize (uint8 (in), 0.5, "bicubic", "Antialiasing", true), uint8 (out), 1)
 %!
 ## Factor 4/3, gives an output of size 10.6667 x 13.3333 rounded up to 11 x 14
 %! out = [116   203   221   141    62   110   180   191    70   104   143   136   111   106
@@ -352,10 +705,15 @@ endfunction
 %! assert (imresize (uint8 (in), 4/3, "bicubic"), uint8 (out), 1)
 %!
 ## Factor 1/3, gives an output of size 2.6667 x 3.3333 rounded up to 3 x 4
-%! out = [ 47    90    24    44
-%!        143    51    21    90
-%!         60    64    56    86];
-%! assert (imresize (uint8 (in), 1/3, "bicubic"), uint8 (out), 1)
+%! out = [ 47   90   24   44
+%!        143   51   21   90
+%!         60   64   56   86];
+%! assert (imresize (uint8 (in), 1/3, "bicubic", "Antialiasing", false), uint8 (out), 1)
+%!
+%! out = [115  135   97  101
+%!        113  119  124  125
+%!         81  157  118   64];
+%! assert (imresize (uint8 (in), 1/3, "bicubic", "Antialiasing", true), uint8 (out), 1)
 
 ## Reduce the size of a checkerboard, such that it turns into gray
 %!test
@@ -363,7 +721,7 @@ endfunction
 %! in = checkerboard (1, [2 2]);
 %! out = [0.5  0.35
 %!        0.5  0.35];
-%! assert ( imresize (in, 0.5, 'bilinear'), out)
+%! assert ( imresize (in, 0.5, "bilinear", "Antialiasing", false), out)
 
 ## check complex inputs (Matlab allows them)
 %!test
@@ -425,3 +783,15 @@ endfunction
 %! out_expected = cat (3, plane_expected_cubic, plane_expected_cubic, plane_expected_cubic);
 %! out_cubic = imresize (in_col, [3, 3], "cubic");
 %! assert (out_cubic, out_expected, 1e-13);
+
+## performance tests for scale down and scale up (commented out by default)
+%!#test
+%! r_big_int16 = randi(16000, 8000, 8000, 'int16');
+%! r_big_double = randi(16000, 8000, 8000);
+%! tic; imresize (r_big_int16,  100/8000); toc;
+%! tic; imresize (r_big_double, 100/8000); toc;
+%!
+%! r_small_int16 = randi(16000, 100, 100, 'int16');
+%! r_small_double = randi(16000, 100, 100);
+%! tic; imresize (r_small_int16,  8000/100); toc;
+%! tic; imresize (r_small_double, 8000/100); toc;
diff --git a/inst/private/interp_method.m b/inst/private/interp_method.m
--- a/inst/private/interp_method.m
+++ b/inst/private/interp_method.m
@@ -26,5 +26,6 @@ function method = interp_method (method)
     case "bicubic",   method = "cubic";
     case "bilinear",  method = "linear";
     case "triangle",  method = "linear";  # interpolation kernel
+    case "box",       method = "nearest";
   endswitch
 endfunction
