Merge pull request #86 from jagapiou/bicubic

Add bicubic interpolation to image.scale.

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 bfbf29a..6d54691 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()
+