Add bicubic interpolation to image.scale.

* Added 'bicubic' method to image.scale.
* Added unit tests of scaling.

Does interpolation for both upsizing and downsizing and does not preserve the
average pixel intensity. Note this is different from 'bilinear' which preserves
the average pixel intensity on downsizing but not on upsizing.

At edges the missing data is added by extending the first/last line segment.

4 files changed, 202 insertions, 29 deletions

diff --git a/README.md b/README.md
index 9c71d7e..3912801 100644
--- a/README.md
+++ b/README.md
@@ -87,7 +87,8 @@ Rescale the height and width of image `src` to have
 width `width` and height `height`.  Variable `mode` specifies 
 type of interpolation to be used. Valid values include 
 [bilinear](https://en.wikipedia.org/wiki/Bilinear_interpolation)
-(the default) or *simple* interpolation. Returns a new `res` Tensor.
+(the default), [bicubic](https://en.wikipedia.org/wiki/Bicubic_interpolation),
+or *simple* interpolation. Returns a new `res` Tensor.
 
 ### [res] image.scale(src, size, [mode]) ###
 Rescale the height and width of image `src`. 
diff --git a/generic/image.c b/generic/image.c
index 27469ab..496b1a9 100755
--- a/generic/image.c
+++ b/generic/image.c
@@ -47,14 +47,14 @@ static long image_(Main_op_depth)( THTensor *T){
   return 1; /* greyscale */
 }
 
-static void image_(Main_scale_rowcol)(THTensor *Tsrc,
-                                      THTensor *Tdst,
-                                      long src_start,
-                                      long dst_start,
-                                      long src_stride,
-                                      long dst_stride,
-                                      long src_len,
-                                      long dst_len ) {
+static void image_(Main_scaleLinear_rowcol)(THTensor *Tsrc,
+                                            THTensor *Tdst,
+                                            long src_start,
+                                            long dst_start,
+                                            long src_stride,
+                                            long dst_stride,
+                                            long src_len,
+                                            long dst_len ) {
 
   real *src= THTensor_(data)(Tsrc);
   real *dst= THTensor_(data)(Tdst);
@@ -110,6 +110,61 @@ static void image_(Main_scale_rowcol)(THTensor *Tsrc,
   }
 }
 
+static void image_(Main_scaleCubic_rowcol)(THTensor *Tsrc,
+                                           THTensor *Tdst,
+                                           long src_start,
+                                           long dst_start,
+                                           long src_stride,
+                                           long dst_stride,
+                                           long src_len,
+                                           long dst_len ) {
+
+  real *src= THTensor_(data)(Tsrc);
+  real *dst= THTensor_(data)(Tdst);
+
+  if ( dst_len == src_len ){
+    long i;
+    for( i = 0; i < dst_len; i++ )
+      dst[ dst_start + i*dst_stride ] = src[ src_start + i*src_stride ];
+  } else {
+    long di;
+    float si_f;
+    long si_i;
+    float scale = (float)(src_len - 1) / (dst_len - 1);
+
+    for( di = 0; di < dst_len - 1; di++ ) {
+      long dst_pos = dst_start + di*dst_stride;
+      si_f = di * scale; si_i = (long)si_f; si_f -= si_i;
+
+      real p0;
+      real p1 = src[ src_start + si_i * src_stride ];
+      real p2 = src[ src_start + (si_i + 1) * src_stride ];
+      real p3;
+      if (si_i > 0) {
+        p0 = src[ src_start + (si_i - 1) * src_stride ];
+      } else {
+        p0 = 2*p1 - p2;
+      }
+      if (si_i + 2 < src_len) {
+        p3 = src[ src_start + (si_i + 2) * src_stride ];
+      } else {
+        p3 = 2*p2 - p1;
+      }
+
+      real a0 = p1;
+      real a1 = -(real)1/(real)2*p0 + (real)1/(real)2*p2;
+      real a2 = p0 - (real)5/(real)2*p1 + (real)2*p2 - (real)1/(real)2*p3;
+      real a3 = -(real)1/(real)2*p0 + (real)3/(real)2*p1 - (real)3/(real)2*p2 +
+                (real)1/(real)2*p3;
+
+      dst[dst_pos] = a0 + si_f * (a1 + si_f * (a2 + a3 * si_f));
+    }
+
+    dst[ dst_start + (dst_len - 1) * dst_stride ] =
+      src[ src_start + (src_len - 1) * src_stride ];
+  }
+}
+
 static int image_(Main_scaleBilinear)(lua_State *L) {
 
   THTensor *Tsrc = luaT_checkudata(L, 1, torch_Tensor);
@@ -147,27 +202,90 @@ static int image_(Main_scaleBilinear)(lua_State *L) {
   for(k=0;k<image_(Main_op_depth)(Tsrc);k++) {
     /* compress/expand rows first */
     for(j = 0; j < src_height; j++) {
-      image_(Main_scale_rowcol)(Tsrc,
-				Ttmp,
-				0*src_stride2+j*src_stride1+k*src_stride0,
-				0*tmp_stride2+j*tmp_stride1+k*tmp_stride0,
-				src_stride2,
-				tmp_stride2,
-				src_width,
-				tmp_width );
+      image_(Main_scaleLinear_rowcol)(Tsrc,
+                                      Ttmp,
+                                      0*src_stride2+j*src_stride1+k*src_stride0,
+                                      0*tmp_stride2+j*tmp_stride1+k*tmp_stride0,
+                                      src_stride2,
+                                      tmp_stride2,
+                                      src_width,
+                                      tmp_width );
+
+    }
+
+    /* then columns */
+    for(i = 0; i < dst_width; i++) {
+      image_(Main_scaleLinear_rowcol)(Ttmp,
+                                      Tdst,
+                                      i*tmp_stride2+0*tmp_stride1+k*tmp_stride0,
+                                      i*dst_stride2+0*dst_stride1+k*dst_stride0,
+                                      tmp_stride1,
+                                      dst_stride1,
+                                      tmp_height,
+                                      dst_height );
+    }
+  }
+  THTensor_(free)(Ttmp);
+  return 0;
+}
+
+static int image_(Main_scaleBicubic)(lua_State *L) {
+
+  THTensor *Tsrc = luaT_checkudata(L, 1, torch_Tensor);
+  THTensor *Tdst = luaT_checkudata(L, 2, torch_Tensor);
+  THTensor *Ttmp;
+  long dst_stride0, dst_stride1, dst_stride2, dst_width, dst_height;
+  long src_stride0, src_stride1, src_stride2, src_width, src_height;
+  long tmp_stride0, tmp_stride1, tmp_stride2, tmp_width, tmp_height;
+  long i, j, k;
+
+  image_(Main_op_validate)(L, Tsrc,Tdst);
+
+  int ndims;
+  if (Tdst->nDimension == 3) ndims = 3;
+  else ndims = 2;
+
+  Ttmp = THTensor_(newWithSize2d)(Tsrc->size[ndims-2], Tdst->size[ndims-1]);
 
+  dst_stride0= image_(Main_op_stride)(Tdst,0);
+  dst_stride1= image_(Main_op_stride)(Tdst,1);
+  dst_stride2= image_(Main_op_stride)(Tdst,2);
+  src_stride0= image_(Main_op_stride)(Tsrc,0);
+  src_stride1= image_(Main_op_stride)(Tsrc,1);
+  src_stride2= image_(Main_op_stride)(Tsrc,2);
+  tmp_stride0= image_(Main_op_stride)(Ttmp,0);
+  tmp_stride1= image_(Main_op_stride)(Ttmp,1);
+  tmp_stride2= image_(Main_op_stride)(Ttmp,2);
+  dst_width=   Tdst->size[ndims-1];
+  dst_height=  Tdst->size[ndims-2];
+  src_width=   Tsrc->size[ndims-1];
+  src_height=  Tsrc->size[ndims-2];
+  tmp_width=   Ttmp->size[1];
+  tmp_height=  Ttmp->size[0];
+
+  for(k=0;k<image_(Main_op_depth)(Tsrc);k++) {
+    /* compress/expand rows first */
+    for(j = 0; j < src_height; j++) {
+      image_(Main_scaleCubic_rowcol)(Tsrc,
+                                     Ttmp,
+                                     0*src_stride2+j*src_stride1+k*src_stride0,
+                                     0*tmp_stride2+j*tmp_stride1+k*tmp_stride0,
+                                     src_stride2,
+                                     tmp_stride2,
+                                     src_width,
+                                     tmp_width );
     }
 
     /* then columns */
     for(i = 0; i < dst_width; i++) {
-      image_(Main_scale_rowcol)(Ttmp,
-				Tdst,
-				i*tmp_stride2+0*tmp_stride1+k*tmp_stride0,
-				i*dst_stride2+0*dst_stride1+k*dst_stride0,
-				tmp_stride1,
-				dst_stride1,
-				tmp_height,
-				dst_height );
+      image_(Main_scaleCubic_rowcol)(Ttmp,
+                                     Tdst,
+                                     i*tmp_stride2+0*tmp_stride1+k*tmp_stride0,
+                                     i*dst_stride2+0*dst_stride1+k*dst_stride0,
+                                     tmp_stride1,
+                                     dst_stride1,
+                                     tmp_height,
+                                     dst_height );
     }
   }
   THTensor_(free)(Ttmp);
@@ -1876,6 +1994,7 @@ int image_(Main_rgb2y)(lua_State *L) {
 static const struct luaL_Reg image_(Main__) [] = {
   {"scaleSimple", image_(Main_scaleSimple)},
   {"scaleBilinear", image_(Main_scaleBilinear)},
+  {"scaleBicubic", image_(Main_scaleBicubic)},
   {"rotate", image_(Main_rotate)},
   {"rotateBilinear", image_(Main_rotateBilinear)},
   {"polar", image_(Main_polar)},
diff --git a/init.lua b/init.lua
index 0819ec0..7f32f98 100644
--- a/init.lua
+++ b/init.lua
@@ -500,15 +500,15 @@ local function scale(...)
                        {type='torch.Tensor', help='input image', req=true},
                        {type='number', help='destination width', req=true},
                        {type='number', help='destination height', req=true},
-                       {type='string', help='mode: bilinear | simple', default='bilinear'},
+                       {type='string', help='mode: bilinear | bicubic |simple', default='bilinear'},
                        '',
                        {type='torch.Tensor', help='input image', req=true},
                        {type='string | number', help='destination size: "WxH" or "MAX" or "^MIN" or MAX', req=true},
-                       {type='string', help='mode: bilinear | simple', default='bilinear'},
+                       {type='string', help='mode: bilinear | bicubic | simple', default='bilinear'},
                        '',
                        {type='torch.Tensor', help='destination image', req=true},
                        {type='torch.Tensor', help='input image', req=true},
-                       {type='string', help='mode: bilinear | simple', default='bilinear'}))
+                       {type='string', help='mode: bilinear | bicubic | simple', default='bilinear'}))
       dok.error('incorrect arguments', 'image.scale')
    end
    if size then
@@ -548,10 +548,12 @@ local function scale(...)
    mode = mode or 'bilinear'
    if mode=='bilinear' then
       src.image.scaleBilinear(src,dst)
+   elseif mode=='bicubic' then
+      src.image.scaleBicubic(src,dst)
    elseif mode=='simple' then
       src.image.scaleSimple(src,dst)
    else
-      dok.error('mode must be one of: simple | bilinear', 'image.scale')
+      dok.error('mode must be one of: simple | bicubic | bilinear', 'image.scale')
    end
    return dst
 end
diff --git a/test/test_scale.lua b/test/test_scale.lua
new file mode 100644
index 0000000..f2225a7
--- /dev/null
+++ b/test/test_scale.lua
@@ -0,0 +1,51 @@
+require 'image'
+require 'torch'
+
+
+torch.setdefaulttensortype('torch.FloatTensor')
+
+local tester = torch.Tester()
+local tests = {}
+
+
+local function outerProduct(x)
+  x = torch.Tensor(x)
+  return torch.ger(x, x)
+end
+
+
+function tests.bilinearUpscale()
+  local im = outerProduct{1, 2, 4, 2}
+  local expected = outerProduct{1, 1.5, 2, 3, 4, 3, 2}
+  local actual = image.scale(im, expected:size(1), expected:size(2), 'bilinear')
+  tester:assertTensorEq(actual, expected, 1e-5)
+end
+
+
+function tests.bilinearDownscale()
+  local im = outerProduct{1, 2, 4, 2}
+  local expected = outerProduct{1.25, 3, 2.5}
+  local actual = image.scale(im, expected:size(1), expected:size(2), 'bilinear')
+  tester:assertTensorEq(actual, expected, 1e-5)
+end
+
+
+function tests.bicubicUpscale()
+  local im = outerProduct{1, 2, 4, 2}
+  local expected = outerProduct{1, 1.4375, 2, 3.1875, 4, 3.25, 2}
+  local actual = image.scale(im, expected:size(1), expected:size(2), 'bicubic')
+  tester:assertTensorEq(actual, expected, 1e-5)
+end
+
+
+function tests.bicubicDownscale()
+  local im = outerProduct{1, 2, 4, 2}
+  local expected = outerProduct{1, 3.1875, 2}
+  local actual = image.scale(im, expected:size(1), expected:size(2), 'bicubic')
+  tester:assertTensorEq(actual, expected, 1e-5)
+end
+
+
+tester:add(tests)
+tester:run()
+