diff --git a/inst/imresize.m b/inst/imresize.m
--- a/inst/imresize.m
+++ b/inst/imresize.m
@@ -71,83 +71,61 @@ function im = imresize (im, scale, metho
   endif
   method = interp_method (method);
 
-  inRows = rows (im);
-  inCols = columns (im);
+  in_rows = rows (im);
+  in_cols = columns (im);
 
-  ## we may be able to use clever indexing instead of interpolation
-  int_row_scale = false;
-  int_col_scale = false;
   if (isscalar (scale))
-    outRows = ceil (rows (im) * scale);
-    outCols = ceil (columns (im) * scale);
-
-    ## check if we can use clever indexing
-    scale_rows = scale_cols = scale;
-    int_row_scale = int_col_scale = is_for_integer (scale);
+    out_rows = ceil (in_rows * scale);
+    out_cols = ceil (in_cols * scale);
+    scale_rows = scale;
+    scale_cols = scale;
   elseif (numel (scale) == 2)
-    outRows = scale(1);
-    outCols = scale(2);
+    out_rows = scale(1);
+    out_cols = scale(2);
+
     ## maintain aspect ratio if requested
-    if (isnan (outRows))
-      outRows = inRows * (outCols / inCols);
-    elseif (isnan (outCols))
-      outCols = inCols * (outRows / inRows);
+    if (isnan (out_rows))
+      out_rows = in_rows * (out_cols / in_cols);
+    elseif (isnan (out_cols))
+      out_cols = in_cols * (out_rows / in_rows);
     endif
-    outRows = ceil (outRows);
-    outCols = ceil (outCols);
 
-    ## we will need this to use clever indexing. In this case, we will also need
-    ## to check that we are changing the rows and columns of the image in the
-    ## same direction
-    scale_rows = (outRows/inRows);
-    scale_cols = (outCols/inCols);
-    int_row_scale = is_for_integer (scale_rows);
-    int_col_scale = is_for_integer (scale_cols);
+    scale_rows = out_rows / in_rows;
+    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
 
-  ## Perform the actual resizing
-  if (inRows == outRows && inCols == outCols)
+  ## calculate the new pixel indices in terms of the old pixel indices
+  off_rows = 1 / scale_rows / 2;
+  off_cols = 1 / scale_cols / 2;
+  idx_rows = 0.5 + off_rows + (0:out_rows-1) / scale_rows;
+  idx_cols = 0.5 + off_cols + (0:out_cols-1) / scale_cols;
+
+  ## work around for imremap, valid for linear interpolation, probably bad
+  ## for cubic interpolation. Limits indices to the bounds 1 and end.
+  idx_rows = max (min (idx_rows, in_rows), 1);
+  idx_cols = max (min (idx_cols, in_cols), 1);
+
+  ## perform the actual resizing
+  if (in_rows == out_rows && in_cols == out_cols)
     ## no resizing to do
-  elseif (strcmpi (method, "nearest") && all ([int_row_scale int_col_scale]))
-    ## we are matlab incompatible here on purpose. We can the stuff here in 2
-    ## ways. With interp2 or by clever indexing. Indexing is much much faster
-    ## than interp2 but they return different results (the way we are doing it
-    ## at least). Matlab does the same as we are doing if the both columns and
-    ## rows go the same direction but if they increase one and decrease the
-    ## other, then they return the same as if we were using interp2. We are
-    ## smarter and use indexing even in that case but then the results differ
-    if (int_row_scale == 1)
-      row_idx = (1:rows (im))(ones (1, scale_rows), :);
-    elseif (int_row_scale == -1)
-      row_idx = ceil (linspace (floor (1/(scale_rows * 2)) + 1, inRows, outRows));
-    endif
-    if (int_col_scale == 1)
-      col_idx = (1:columns (im))(ones (scale_cols, 1), :);
-    elseif (int_col_scale == -1)
-      col_idx = ceil (linspace (floor (1/(scale_cols * 2)) + 1, inCols, outCols));
-    endif
-    im = im(row_idx, col_idx);
+  elseif (strcmpi (method, "nearest"))
+    im = im(round(idx_rows), round(idx_cols), :);
 
   else
-    [XI, YI] = meshgrid (linspace (1, inCols, outCols), linspace (1, inRows, outRows));
-    im = imremap (im, XI, YI, method);
+    idx_cols = repmat (idx_cols, [out_rows 1]);
+    idx_rows = repmat (idx_rows', [1 out_cols]);
+    im = imremap (im, idx_cols, idx_rows, method);
   endif
 endfunction
 
-function retval = is_for_integer (scale)
-  retval = false;
-  if (fix (scale) == scale)
-    retval = 1;
-  elseif (fix (1/scale) == (1/scale))
-    ## if scale/n is an integer then we are resizing to one half, one third, etc
-    ## and we can also use clever indexing
-    retval = -1;
-  endif
-endfunction
-
+## Test basic features.
 %!test
+%!
+## Test scaling with 1:
 %! 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
@@ -158,6 +136,7 @@ endfunction
 %! assert (imresize (uint8 (in), 1, "nearest"), uint8 (in))
 %! assert (imresize (uint8 (in), 1, "bicubic"), uint8 (in))
 %!
+## Test nearest neighbour with a scale factor of 2, also by aspect ratio preservation:
 %! out = [116  116  227  227  153  153   69   69  146  146  194  194   59   59  130  130  139  139  106  106
 %!        116  116  227  227  153  153   69   69  146  146  194  194   59   59  130  130  139  139  106  106
 %!          2    2   47   47  137  137  249  249   90   90   75   75   16   16   24   24  158  158   44   44
@@ -177,6 +156,7 @@ endfunction
 %! assert (imresize (uint8 (in), [14 NaN], "nearest"), uint8 (out))
 %! assert (imresize (uint8 (in), [NaN 20], "nearest"), uint8 (out))
 %!
+## Test nearest neighbour with a scaling of 2 for x and 1 for y:
 %! out = [116  116  227  227  153  153   69   69  146  146  194  194   59   59  130  130  139  139  106  106
 %!          2    2   47   47  137  137  249  249   90   90   75   75   16   16   24   24  158  158   44   44
 %!        155  155   68   68   46   46   84   84  166  166  156  156   69   69  204  204   32   32  152  152
@@ -186,96 +166,169 @@ endfunction
 %!          9    9   47   47  133  133  135  135  130  130  207  207  236  236   43   43   19   19   73   73];
 %! assert (imresize (uint8 (in), [7 20], "nearest"), uint8 (out))
 %!
+## Test nearest neighbour with a scaling of 1 for x and 2 for y:
+%! out = [116  227  153   69  146  194   59  130  139  106
+%!        116  227  153   69  146  194   59  130  139  106
+%!          2   47  137  249   90   75   16   24  158   44
+%!          2   47  137  249   90   75   16   24  158   44
+%!        155   68   46   84  166  156   69  204   32  152
+%!        155   68   46   84  166  156   69  204   32  152
+%!         71  221  137  230  210  153  192  115   30  118
+%!         71  221  137  230  210  153  192  115   30  118
+%!        107  143  108   52   51   73  101   21  175   90
+%!        107  143  108   52   51   73  101   21  175   90
+%!         54  158  143   77   26  168  113  229  165  225
+%!         54  158  143   77   26  168  113  229  165  225
+%!          9   47  133  135  130  207  236   43   19   73
+%!          9   47  133  135  130  207  236   43   19   73];
+%! assert (imresize (uint8 (in), [14 10], "nearest"), uint8 (out))
+%!
+## Test equivalence of different input writing styles:
 %! 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), 1.5, "linear"), imresize (uint8 (in), 1.5, "LIneAR"))
 %! assert (imresize (uint8 (in), 1.5, "linear"), imresize (uint8 (in), 1.5, "triangle"))
 %!
-%! out = [ 47  249   75   24   44
-%!        221  230  153  115  118
-%!        158   77  168  229  225
-%!        47   135  207   43   73];
-%! assert (imresize (uint8 (in), 0.5, "nearest"), uint8 (out))
+## nearest neighbour test, check that a pixel from the neighbour locations
+## gets picked. The distance is the same for all neighbours here.
+%! 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
+%!        71  221  137  230  210  153  192  115   30  118
+%!       107  143  108   52   51   73  101   21  175   90
+%!        54  158  143   77   26  168  113  229  165  225
+%!         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");
+%! 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
 
-## Do not enforce floating point images to be in the [0 1] range (bug #43846)
-%!assert (imresize (repmat (5, [3 3]), 2), repmat (5, [6 6]), eps*100)
 
+## Test floating point range and and scaling of multi-channel images.
+%!test
+%!
+## Do not enforce floating point images to be in the [0 1] range (bug #43846):
+%! assert (imresize (repmat (5, [3 3]), 2), repmat (5, [6 6]), eps*100)
+%!
 ## Similarly, do not enforce images to have specific dimensions and only
-## expand on the first 2 dimensions.
-%!assert (imresize (repmat (5, [3 3 2]), 2), repmat (5, [6 6 2]), eps*100)
+## expand on the first 2 dimensions:
+%! assert (imresize (repmat (5, [3 3 2]), 2), repmat (5, [6 6 2]), eps*100)
 
-## The following are the matlab results. We have slighlty different
-## results but not by much. If there's would be any fixes, they
-## would have to be on interp2 or maybe in imremap.
+## Test that scaling a multi-channel image is equivalent to scaling its channels.
+%!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
+%!     endfor
+%! endfor
 
-%!shared in
+## 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
+%!
 %! 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
 %!        71  221  137  230  210  153  192  115   30  118
 %!       107  143  108   52   51   73  101   21  175   90
 %!        54  158  143   77   26  168  113  229  165  225
-%!         9   47  133  135  130  207  236   43   19   73];
+%!         9   47  133  135  130  207  236   43   19   73
+%!       129   60   59  243   64  181  249   56   32   86];
+%!
+## 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
+%!         27   39   62  110  172  202  123   96   78   36   40   68  123  100   62
+%!        129   97   64   62   87  119  146  148  128   74  117  154   73   94  134
+%!        113  129  136  101  125  162  183  172  151  135  146  139   53   83  135
+%!         77  143  195  145  166  197  186  162  146  171  138   92   62   84  113
+%!        101  129  149  120   98   81   78   82   91  111   77   56  132  123   95
+%!         81  116  147  130   96   61   43   80  119  109  116  132  162  164  158
+%!         46   93  139  141  114   80   50  109  168  141  166  189  151  171  200
+%!         16   41   77  123  130  123  115  157  204  214  145   69   48   71   98
+%!         69   62   61   89  143  174  112  146  202  235  147   46   30   53   80
+%!        129   95   60   59  151  213   94  123  192  238  153   52   36   59   86];
+%! assert (max (max (abs (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 (max (max (abs (imresize (uint8 (in), 0.5, "bilinear") - 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
+%!         45   89  126  148  177  138  114  109   43   52   97  141   78   67
+%!         59   57   73  114  177  145  114   96   45   71   99  108   88   85
+%!        145  109   76   63   96  146  166  147   93  152  133   47  134  148
+%!         82  157  174  137  201  208  186  156  174  145   90   42  111  122
+%!         94  143  152  119  119  114  108  107  131   86   80  119  104  101
+%!         87  126  139  114   69   49   67  109  106  102  126  167  145  141
+%!         48  108  143  135   91   56   89  167  134  177  184  154  199  206
+%!         15   44   88  133  129  121  149  204  219  124   55   44   85   92
+%!         84   66   67  102  189  132  127  198  237  123   42   34   74   81
+%!        129   86   60   82  220  131  108  190  241  128   47   39   79   86];
+%! assert (max (max (abs (imresize (uint8 (in), 4/3, "bilinear") - 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 (max (max (abs (imresize (uint8 (in), 1/3, "bilinear") - uint8 (out)))) <= 1)
 
+
+## Comparison to the reference result from matlab with cubic interpolation
+## fails. This has to be fixed in imremap and here (remove the work around
+## above).
 %!xtest
-%! out = [116  185  235  171   96   64  134  189  186   74   90   141  140  124  108
-%!         44   92  143  149  164  163  119  123  118   44   38    80  151  118   63
-%!         14   21   47  107  195  228  115   81   70   24   19    56  137  105   49
-%!        145   98   49   49   71  107  148  159  132   58  124   176   61   85  145
-%!        118  139  144   92  116  168  201  188  159  140  167   158   27   69  152
-%!         62  151  218  145  174  219  201  164  146  187  148    84   48   76  115
-%!        102  132  151  119   90   72   72   72   83  114   60    31  144  130   81
-%!         82  121  154  133   87   41   19   67  116   95  108   140  183  180  164
-%!         40   96  152  149  117   74   34  108  179  131  175   215  153  177  219
-%!         11   33   73  127  137  125  113  158  212  229  148    55   35   63   96
-%!          4   17   53  121  141  138  133  171  220  253  141    16    7   36   67];
+%!
+%! 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
+%!        71  221  137  230  210  153  192  115   30  118
+%!       107  143  108   52   51   73  101   21  175   90
+%!        54  158  143   77   26  168  113  229  165  225
+%!         9   47  133  135  130  207  236   43   19   73
+%!       129   60   59  243   64  181  249   56   32   86];
+%!
+%! out = [116  187  237  171   94   61  135  191  187   75   91  142  140  124  108
+%!         43   92  143  149  164  163  119  123  118   44   38   80  151  118   62
+%!         13   21   47  107  195  228  115   81   70   24   19   56  137  105   48
+%!        146   98   49   49   71  107  148  159  132   58  124  176   61   85  146
+%!        118  139  144   92  116  168  201  188  159  140  167  158   27   69  153
+%!         61  151  218  145  174  219  201  164  146  187  148   84   48   76  115
+%!        102  132  151  119   90   72   72   72   83  114   60   31  144  130   80
+%!         81  121  154  133   87   41   19   67  116   95  108  140  183  180  163
+%!         37   95  152  150  117   73   35  108  179  130  174  214  153  176  219
+%!          3   29   73  131  136  120  116  162  214  229  147   54   35   62   96
+%!         67   54   51   83  153  187  111  141  210  255  149   22   13   42   74
+%!        142   99   53   43  164  237   77  103  197  254  159   42   31   59   91];
 %! assert (imresize (uint8 (in), 1.5, "bicubic"), uint8 (out))
 
-%!xtest
-%! 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
-%!         27   39   62  110  172  202  123   96   78   36   40    68  123  100   62
-%!        129   97   64   62   87  119  146  148  128   74  117   154   73   94  134
-%!        113  129  136  101  125  162  183  172  151  135  146   139   53   83  135
-%!         77  143  195  145  166  197  186  162  146  171  138    92   62   84  113
-%!        101  129  149  120   98   81   78   82   91  111   77    56  132  123   95
-%!         81  116  147  130   96   61   43   80  119  109  116   132  162  164  158
-%!         46   93  139  141  114   80   50  109  168  141  166   189  151  171  200
-%!         16   41   77  123  130  123  115  157  204  214  145    69   48   71   98
-%!          9   28   61  119  134  134  131  169  212  231  140    39   23   46   73];
-%! assert (imresize (uint8 (in), 1.5, "bilinear"), uint8 (out))
 
-%!xtest
-%! out = [108  136  125   89  107
-%!        111  132  143  114   99
-%!        106  110  106  127  136
-%!         47  121  163  138   68];
-%! assert (imresize (uint8 (in), 0.5, "bilinear"), uint8 (out))
-
-%!xtest
-%! out = [103  141  124   78  110
-%!        111  134  153  114   91
-%!        115  108   93  128  146
-%!         38  124  175  143   54];
-%! assert (imresize (uint8 (in), 0.5, "bicubic"), uint8 (out))
-
-%!xtest
-%! ## bug #55202
-%! assert (imresize (zeros (1, 20), [1 30], "nearest"), zeros (1, 30))
-
+## Reduce the size of a checkerboard, such that it turns into gray
 %!test
-%! ## a test for bug #55780
-%! a = zeros (7, "uint8");
-%! a (4,4) = 255;
-%! c = [   4   13   24   31   31   24   13    4;
-%!        13   44   79  103  103   79   44   13;
-%!        24   79  140  182  182  140   79   24;
-%!        31  103  182  238  238  182  103   31;
-%!        31  103  182  238  238  182  103   31;
-%!        24   79  140  182  182  140   79   24;
-%!        13   44   79  103  103   79   44   13;
-%!         4   13   24   31   31   24   13    4];
-%! b = imresize (a, 4, "bicubic");
-%! assert (b(11:18,11:18), c, eps)
-
+%!
+%! in = checkerboard (1, [2 2]);
+%! out = [0.5  0.35
+%!        0.5  0.35];
+%! assert ( imresize (in, 0.5, 'bilinear'), out)
