avfilter/scale_cuda: add bicubic interpolation

5 files changed, 368 insertions, 36 deletions

diff --git a/libavfilter/Makefile b/libavfilter/Makefile
index 1e60c55f6f..65d03f9191 100644
--- a/libavfilter/Makefile
+++ b/libavfilter/Makefile
@@ -374,7 +374,8 @@ OBJS-$(CONFIG_ROBERTS_OPENCL_FILTER)         += vf_convolution_opencl.o opencl.o
 OBJS-$(CONFIG_ROTATE_FILTER)                 += vf_rotate.o
 OBJS-$(CONFIG_SAB_FILTER)                    += vf_sab.o
 OBJS-$(CONFIG_SCALE_FILTER)                  += vf_scale.o scale_eval.o
-OBJS-$(CONFIG_SCALE_CUDA_FILTER)             += vf_scale_cuda.o vf_scale_cuda.ptx.o scale_eval.o
+OBJS-$(CONFIG_SCALE_CUDA_FILTER)             += vf_scale_cuda.o scale_eval.o \
+                                                vf_scale_cuda.ptx.o vf_scale_cuda_bicubic.ptx.o
 OBJS-$(CONFIG_SCALE_NPP_FILTER)              += vf_scale_npp.o scale_eval.o
 OBJS-$(CONFIG_SCALE_QSV_FILTER)              += vf_scale_qsv.o
 OBJS-$(CONFIG_SCALE_VAAPI_FILTER)            += vf_scale_vaapi.o scale_eval.o vaapi_vpp.o
diff --git a/libavfilter/cuda/vector_helpers.cuh b/libavfilter/cuda/vector_helpers.cuh
new file mode 100644
index 0000000000..67332ef030
--- /dev/null
+++ b/libavfilter/cuda/vector_helpers.cuh
@@ -0,0 +1,112 @@
+/*
+ * This file is part of FFmpeg.
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included in
+ * all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ */
+
+#ifndef AVFILTER_CUDA_VECTORHELPERS_H
+#define AVFILTER_CUDA_VECTORHELPERS_H
+
+typedef unsigned char uchar;
+typedef unsigned short ushort;
+
+template<typename T> struct vector_helper { };
+template<> struct vector_helper<uchar>   { typedef float  ftype; typedef int  itype; };
+template<> struct vector_helper<uchar2>  { typedef float2 ftype; typedef int2 itype; };
+template<> struct vector_helper<uchar4>  { typedef float4 ftype; typedef int4 itype; };
+template<> struct vector_helper<ushort>  { typedef float  ftype; typedef int  itype; };
+template<> struct vector_helper<ushort2> { typedef float2 ftype; typedef int2 itype; };
+template<> struct vector_helper<ushort4> { typedef float4 ftype; typedef int4 itype; };
+template<> struct vector_helper<int>     { typedef float  ftype; typedef int  itype; };
+template<> struct vector_helper<int2>    { typedef float2 ftype; typedef int2 itype; };
+template<> struct vector_helper<int4>    { typedef float4 ftype; typedef int4 itype; };
+
+#define floatT typename vector_helper<T>::ftype
+#define intT typename vector_helper<T>::itype
+
+template<typename T, typename V> inline __device__ V to_floatN(const T &a) { return (V)a; }
+template<typename T, typename V> inline __device__ T from_floatN(const V &a) { return (T)a; }
+
+#define OPERATORS2(T) \
+    template<typename V> inline __device__ T operator+(const T &a, const V &b) { return make_ ## T (a.x + b.x, a.y + b.y); } \
+    template<typename V> inline __device__ T operator-(const T &a, const V &b) { return make_ ## T (a.x - b.x, a.y - b.y); } \
+    template<typename V> inline __device__ T operator*(const T &a, V b) { return make_ ## T (a.x * b, a.y * b); } \
+    template<typename V> inline __device__ T operator/(const T &a, V b) { return make_ ## T (a.x / b, a.y / b); } \
+    template<typename V> inline __device__ T operator>>(const T &a, V b) { return make_ ## T (a.x >> b, a.y >> b); } \
+    template<typename V> inline __device__ T operator<<(const T &a, V b) { return make_ ## T (a.x << b, a.y << b); } \
+    template<typename V> inline __device__ T &operator+=(T &a, const V &b) { a.x += b.x; a.y += b.y; return a; } \
+    template<typename V> inline __device__ void vec_set(T &a, const V &b) { a.x = b.x; a.y = b.y; } \
+    template<typename V> inline __device__ void vec_set_scalar(T &a, V b) { a.x = b; a.y = b; } \
+    template<> inline __device__ float2 to_floatN<T, float2>(const T &a) { return make_float2(a.x, a.y); } \
+    template<> inline __device__ T from_floatN<T, float2>(const float2 &a) { return make_ ## T(a.x, a.y); }
+#define OPERATORS4(T) \
+    template<typename V> inline __device__ T operator+(const T &a, const V &b) { return make_ ## T (a.x + b.x, a.y + b.y, a.z + b.z, a.w + b.w); } \
+    template<typename V> inline __device__ T operator-(const T &a, const V &b) { return make_ ## T (a.x - b.x, a.y - b.y, a.z - b.z, a.w - b.w); } \
+    template<typename V> inline __device__ T operator*(const T &a, V b) { return make_ ## T (a.x * b, a.y * b, a.z * b, a.w * b); } \
+    template<typename V> inline __device__ T operator/(const T &a, V b) { return make_ ## T (a.x / b, a.y / b, a.z / b, a.w / b); } \
+    template<typename V> inline __device__ T operator>>(const T &a, V b) { return make_ ## T (a.x >> b, a.y >> b, a.z >> b, a.w >> b); } \
+    template<typename V> inline __device__ T operator<<(const T &a, V b) { return make_ ## T (a.x << b, a.y << b, a.z << b, a.w << b); } \
+    template<typename V> inline __device__ T &operator+=(T &a, const V &b) { a.x += b.x; a.y += b.y; a.z += b.z; a.w += b.w; return a; } \
+    template<typename V> inline __device__ void vec_set(T &a, const V &b) { a.x = b.x; a.y = b.y; a.z = b.z; a.w = b.w; } \
+    template<typename V> inline __device__ void vec_set_scalar(T &a, V b) { a.x = b; a.y = b; a.z = b; a.w = b; } \
+    template<> inline __device__ float4 to_floatN<T, float4>(const T &a) { return make_float4(a.x, a.y, a.z, a.w); } \
+    template<> inline __device__ T from_floatN<T, float4>(const float4 &a) { return make_ ## T(a.x, a.y, a.z, a.w); }
+
+OPERATORS2(int2)
+OPERATORS2(uchar2)
+OPERATORS2(ushort2)
+OPERATORS2(float2)
+OPERATORS4(int4)
+OPERATORS4(uchar4)
+OPERATORS4(ushort4)
+OPERATORS4(float4)
+
+template<typename V> inline __device__ void vec_set(int &a, V b) { a = b; }
+template<typename V> inline __device__ void vec_set(float &a, V b) { a = b; }
+template<typename V> inline __device__ void vec_set(uchar &a, V b) { a = b; }
+template<typename V> inline __device__ void vec_set(ushort &a, V b) { a = b; }
+template<typename V> inline __device__ void vec_set_scalar(int &a, V b) { a = b; }
+template<typename V> inline __device__ void vec_set_scalar(float &a, V b) { a = b; }
+template<typename V> inline __device__ void vec_set_scalar(uchar &a, V b) { a = b; }
+template<typename V> inline __device__ void vec_set_scalar(ushort &a, V b) { a = b; }
+
+template<typename T>
+inline __device__ T lerp_scalar(T v0, T v1, float t) {
+    return t*v1 + (1.0f - t)*v0;
+}
+
+template<>
+inline __device__ float2 lerp_scalar<float2>(float2 v0, float2 v1, float t) {
+    return make_float2(
+        lerp_scalar(v0.x, v1.x, t),
+        lerp_scalar(v0.y, v1.y, t)
+    );
+}
+
+template<>
+inline __device__ float4 lerp_scalar<float4>(float4 v0, float4 v1, float t) {
+    return make_float4(
+        lerp_scalar(v0.x, v1.x, t),
+        lerp_scalar(v0.y, v1.y, t),
+        lerp_scalar(v0.z, v1.z, t),
+        lerp_scalar(v0.w, v1.w, t)
+    );
+}
+
+#endif
diff --git a/libavfilter/version.h b/libavfilter/version.h
index b8ba489da7..2db35f85af 100644
--- a/libavfilter/version.h
+++ b/libavfilter/version.h
@@ -31,7 +31,7 @@
 
 #define LIBAVFILTER_VERSION_MAJOR   7
 #define LIBAVFILTER_VERSION_MINOR  88
-#define LIBAVFILTER_VERSION_MICRO 100
+#define LIBAVFILTER_VERSION_MICRO 101
 
 
 #define LIBAVFILTER_VERSION_INT AV_VERSION_INT(LIBAVFILTER_VERSION_MAJOR, \
diff --git a/libavfilter/vf_scale_cuda.c b/libavfilter/vf_scale_cuda.c
index 9d59ed4863..b287bd8c12 100644
--- a/libavfilter/vf_scale_cuda.c
+++ b/libavfilter/vf_scale_cuda.c
@@ -55,6 +55,15 @@ static const enum AVPixelFormat supported_formats[] = {
 
 #define CHECK_CU(x) FF_CUDA_CHECK_DL(ctx, s->hwctx->internal->cuda_dl, x)
 
+enum {
+    INTERP_ALGO_DEFAULT,
+
+    INTERP_ALGO_BILINEAR,
+    INTERP_ALGO_BICUBIC,
+
+    INTERP_ALGO_COUNT
+};
+
 typedef struct CUDAScaleContext {
     const AVClass *class;
 
@@ -98,6 +107,9 @@ typedef struct CUDAScaleContext {
     CUdeviceptr srcBuffer;
     CUdeviceptr dstBuffer;
     int         tex_alignment;
+
+    int interp_algo;
+    int interp_use_linear;
 } CUDAScaleContext;
 
 static av_cold int cudascale_init(AVFilterContext *ctx)
@@ -269,10 +281,32 @@ static av_cold int cudascale_config_props(AVFilterLink *outlink)
     AVCUDADeviceContext *device_hwctx = frames_ctx->device_ctx->hwctx;
     CUcontext dummy, cuda_ctx = device_hwctx->cuda_ctx;
     CudaFunctions *cu = device_hwctx->internal->cuda_dl;
+    char buf[64];
     int w, h;
     int ret;
 
+    char *scaler_ptx;
+    const char *function_infix = "";
+
     extern char vf_scale_cuda_ptx[];
+    extern char vf_scale_cuda_bicubic_ptx[];
+
+    switch(s->interp_algo) {
+    case INTERP_ALGO_BILINEAR:
+        scaler_ptx = vf_scale_cuda_ptx;
+        function_infix = "_Bilinear";
+        s->interp_use_linear = 1;
+        break;
+    case INTERP_ALGO_DEFAULT:
+    case INTERP_ALGO_BICUBIC:
+        scaler_ptx = vf_scale_cuda_bicubic_ptx;
+        function_infix = "_Bicubic";
+        s->interp_use_linear = 0;
+        break;
+    default:
+        av_log(ctx, AV_LOG_ERROR, "Unknown interpolation algorithm\n");
+        return AVERROR_BUG;
+    }
 
     s->hwctx = device_hwctx;
     s->cu_stream = s->hwctx->stream;
@@ -281,31 +315,37 @@ static av_cold int cudascale_config_props(AVFilterLink *outlink)
     if (ret < 0)
         goto fail;
 
-    ret = CHECK_CU(cu->cuModuleLoadData(&s->cu_module, vf_scale_cuda_ptx));
+    ret = CHECK_CU(cu->cuModuleLoadData(&s->cu_module, scaler_ptx));
     if (ret < 0)
         goto fail;
 
-    CHECK_CU(cu->cuModuleGetFunction(&s->cu_func_uchar, s->cu_module, "Subsample_Bilinear_uchar"));
+    snprintf(buf, sizeof(buf), "Subsample%s_uchar", function_infix);
+    CHECK_CU(cu->cuModuleGetFunction(&s->cu_func_uchar, s->cu_module, buf));
     if (ret < 0)
         goto fail;
 
-    CHECK_CU(cu->cuModuleGetFunction(&s->cu_func_uchar2, s->cu_module, "Subsample_Bilinear_uchar2"));
+    snprintf(buf, sizeof(buf), "Subsample%s_uchar2", function_infix);
+    CHECK_CU(cu->cuModuleGetFunction(&s->cu_func_uchar2, s->cu_module, buf));
     if (ret < 0)
         goto fail;
 
-    CHECK_CU(cu->cuModuleGetFunction(&s->cu_func_uchar4, s->cu_module, "Subsample_Bilinear_uchar4"));
+    snprintf(buf, sizeof(buf), "Subsample%s_uchar4", function_infix);
+    CHECK_CU(cu->cuModuleGetFunction(&s->cu_func_uchar4, s->cu_module, buf));
     if (ret < 0)
         goto fail;
 
-    CHECK_CU(cu->cuModuleGetFunction(&s->cu_func_ushort, s->cu_module, "Subsample_Bilinear_ushort"));
+    snprintf(buf, sizeof(buf), "Subsample%s_ushort", function_infix);
+    CHECK_CU(cu->cuModuleGetFunction(&s->cu_func_ushort, s->cu_module, buf));
     if (ret < 0)
         goto fail;
 
-    CHECK_CU(cu->cuModuleGetFunction(&s->cu_func_ushort2, s->cu_module, "Subsample_Bilinear_ushort2"));
+    snprintf(buf, sizeof(buf), "Subsample%s_ushort2", function_infix);
+    CHECK_CU(cu->cuModuleGetFunction(&s->cu_func_ushort2, s->cu_module, buf));
     if (ret < 0)
         goto fail;
 
-    CHECK_CU(cu->cuModuleGetFunction(&s->cu_func_ushort4, s->cu_module, "Subsample_Bilinear_ushort4"));
+    snprintf(buf, sizeof(buf), "Subsample%s_ushort4", function_infix);
+    CHECK_CU(cu->cuModuleGetFunction(&s->cu_func_ushort4, s->cu_module, buf));
     if (ret < 0)
         goto fail;
 
@@ -352,17 +392,19 @@ fail:
 static int call_resize_kernel(AVFilterContext *ctx, CUfunction func, int channels,
                               uint8_t *src_dptr, int src_width, int src_height, int src_pitch,
                               uint8_t *dst_dptr, int dst_width, int dst_height, int dst_pitch,
-                              int pixel_size)
+                              int pixel_size, int bit_depth)
 {
     CUDAScaleContext *s = ctx->priv;
     CudaFunctions *cu = s->hwctx->internal->cuda_dl;
     CUdeviceptr dst_devptr = (CUdeviceptr)dst_dptr;
     CUtexObject tex = 0;
-    void *args_uchar[] = { &tex, &dst_devptr, &dst_width, &dst_height, &dst_pitch, &src_width, &src_height };
+    void *args_uchar[] = { &tex, &dst_devptr, &dst_width, &dst_height, &dst_pitch, &src_width, &src_height, &src_pitch, &bit_depth };
     int ret;
 
     CUDA_TEXTURE_DESC tex_desc = {
-        .filterMode = CU_TR_FILTER_MODE_LINEAR,
+        .filterMode = s->interp_use_linear ?
+                      CU_TR_FILTER_MODE_LINEAR :
+                      CU_TR_FILTER_MODE_POINT,
         .flags = CU_TRSF_READ_AS_INTEGER,
     };
 
@@ -404,73 +446,73 @@ static int scalecuda_resize(AVFilterContext *ctx,
         call_resize_kernel(ctx, s->cu_func_uchar, 1,
                            in->data[0], in->width, in->height, in->linesize[0],
                            out->data[0], out->width, out->height, out->linesize[0],
-                           1);
+                           1, 8);
         call_resize_kernel(ctx, s->cu_func_uchar, 1,
-                           in->data[1], in->width/2, in->height/2, in->linesize[0]/2,
-                           out->data[1], out->width/2, out->height/2, out->linesize[0]/2,
-                           1);
+                           in->data[1], in->width / 2, in->height / 2, in->linesize[0] / 2,
+                           out->data[1], out->width / 2, out->height / 2, out->linesize[0] / 2,
+                           1, 8);
         call_resize_kernel(ctx, s->cu_func_uchar, 1,
-                           in->data[2], in->width/2, in->height/2, in->linesize[0]/2,
-                           out->data[2], out->width/2, out->height/2, out->linesize[0]/2,
-                           1);
+                           in->data[2], in->width / 2, in->height / 2, in->linesize[0] / 2,
+                           out->data[2], out->width / 2, out->height / 2, out->linesize[0] / 2,
+                           1, 8);
         break;
     case AV_PIX_FMT_YUV444P:
         call_resize_kernel(ctx, s->cu_func_uchar, 1,
                            in->data[0], in->width, in->height, in->linesize[0],
                            out->data[0], out->width, out->height, out->linesize[0],
-                           1);
+                           1, 8);
         call_resize_kernel(ctx, s->cu_func_uchar, 1,
                            in->data[1], in->width, in->height, in->linesize[0],
                            out->data[1], out->width, out->height, out->linesize[0],
-                           1);
+                           1, 8);
         call_resize_kernel(ctx, s->cu_func_uchar, 1,
                            in->data[2], in->width, in->height, in->linesize[0],
                            out->data[2], out->width, out->height, out->linesize[0],
-                           1);
+                           1, 8);
         break;
     case AV_PIX_FMT_YUV444P16:
         call_resize_kernel(ctx, s->cu_func_ushort, 1,
                            in->data[0], in->width, in->height, in->linesize[0] / 2,
                            out->data[0], out->width, out->height, out->linesize[0] / 2,
-                           2);
+                           2, 16);
         call_resize_kernel(ctx, s->cu_func_ushort, 1,
                            in->data[1], in->width, in->height, in->linesize[1] / 2,
                            out->data[1], out->width, out->height, out->linesize[1] / 2,
-                           2);
+                           2, 16);
         call_resize_kernel(ctx, s->cu_func_ushort, 1,
                            in->data[2], in->width, in->height, in->linesize[2] / 2,
                            out->data[2], out->width, out->height, out->linesize[2] / 2,
-                           2);
+                           2, 16);
         break;
     case AV_PIX_FMT_NV12:
         call_resize_kernel(ctx, s->cu_func_uchar, 1,
                            in->data[0], in->width, in->height, in->linesize[0],
                            out->data[0], out->width, out->height, out->linesize[0],
-                           1);
+                           1, 8);
         call_resize_kernel(ctx, s->cu_func_uchar2, 2,
-                           in->data[1], in->width/2, in->height/2, in->linesize[1],
-                           out->data[1], out->width/2, out->height/2, out->linesize[1]/2,
-                           1);
+                           in->data[1], in->width / 2, in->height / 2, in->linesize[1],
+                           out->data[1], out->width / 2, out->height / 2, out->linesize[1] / 2,
+                           1, 8);
         break;
     case AV_PIX_FMT_P010LE:
         call_resize_kernel(ctx, s->cu_func_ushort, 1,
-                           in->data[0], in->width, in->height, in->linesize[0]/2,
-                           out->data[0], out->width, out->height, out->linesize[0]/2,
-                           2);
+                           in->data[0], in->width, in->height, in->linesize[0] / 2,
+                           out->data[0], out->width, out->height, out->linesize[0] / 2,
+                           2, 10);
         call_resize_kernel(ctx, s->cu_func_ushort2, 2,
-                           in->data[1], in->width / 2, in->height / 2, in->linesize[1]/2,
+                           in->data[1], in->width / 2, in->height / 2, in->linesize[1] / 2,
                            out->data[1], out->width / 2, out->height / 2, out->linesize[1] / 4,
-                           2);
+                           2, 10);
         break;
     case AV_PIX_FMT_P016LE:
         call_resize_kernel(ctx, s->cu_func_ushort, 1,
                            in->data[0], in->width, in->height, in->linesize[0] / 2,
                            out->data[0], out->width, out->height, out->linesize[0] / 2,
-                           2);
+                           2, 16);
         call_resize_kernel(ctx, s->cu_func_ushort2, 2,
                            in->data[1], in->width / 2, in->height / 2, in->linesize[1] / 2,
                            out->data[1], out->width / 2, out->height / 2, out->linesize[1] / 4,
-                           2);
+                           2, 16);
         break;
     default:
         return AVERROR_BUG;
@@ -552,6 +594,9 @@ fail:
 static const AVOption options[] = {
     { "w",      "Output video width",  OFFSET(w_expr),     AV_OPT_TYPE_STRING, { .str = "iw"   }, .flags = FLAGS },
     { "h",      "Output video height", OFFSET(h_expr),     AV_OPT_TYPE_STRING, { .str = "ih"   }, .flags = FLAGS },
+    { "interp_algo", "Interpolation algorithm used for resizing", OFFSET(interp_algo), AV_OPT_TYPE_INT, { .i64 = INTERP_ALGO_DEFAULT }, 0, INTERP_ALGO_COUNT - 1, FLAGS, "interp_algo" },
+        { "bilinear",    "bilinear",     0, AV_OPT_TYPE_CONST, { .i64 = INTERP_ALGO_BILINEAR }, 0, 0, FLAGS, "interp_algo" },
+        { "bicubic",     "bicubic",      0, AV_OPT_TYPE_CONST, { .i64 = INTERP_ALGO_BICUBIC  }, 0, 0, FLAGS, "interp_algo" },
     { "force_original_aspect_ratio", "decrease or increase w/h if necessary to keep the original AR", OFFSET(force_original_aspect_ratio), AV_OPT_TYPE_INT, { .i64 = 0}, 0, 2, FLAGS, "force_oar" },
     { "disable",  NULL, 0, AV_OPT_TYPE_CONST, {.i64 = 0 }, 0, 0, FLAGS, "force_oar" },
     { "decrease", NULL, 0, AV_OPT_TYPE_CONST, {.i64 = 1 }, 0, 0, FLAGS, "force_oar" },
diff --git a/libavfilter/vf_scale_cuda_bicubic.cu b/libavfilter/vf_scale_cuda_bicubic.cu
new file mode 100644
index 0000000000..8a27927e60
--- /dev/null
+++ b/libavfilter/vf_scale_cuda_bicubic.cu
@@ -0,0 +1,174 @@
+/*
+ * This file is part of FFmpeg.
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included in
+ * all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ */
+
+#include "cuda/vector_helpers.cuh"
+
+__device__ inline float4 bicubic_coeffs(float x)
+{
+    const float A = -0.75f;
+
+    float4 res;
+    res.x = ((A * (x + 1) - 5 * A) * (x + 1) + 8 * A) * (x + 1) - 4 * A;
+    res.y = ((A + 2) * x - (A + 3)) * x * x + 1;
+    res.z = ((A + 2) * (1 - x) - (A + 3)) * (1 - x) * (1 - x) + 1;
+    res.w = 1.0f - res.x - res.y - res.z;
+
+    return res;
+}
+
+__device__ inline void bicubic_fast_coeffs(float x, float *h0, float *h1, float *s)
+{
+    float4 coeffs = bicubic_coeffs(x);
+
+    float g0 = coeffs.x + coeffs.y;
+    float g1 = coeffs.z + coeffs.w;
+
+    *h0 = coeffs.y / g0 - 0.5f;
+    *h1 = coeffs.w / g1 + 1.5f;
+    *s  = g0 / (g0 + g1);
+}
+
+template<typename V>
+__device__ inline V bicubic_filter(float4 coeffs, V c0, V c1, V c2, V c3)
+{
+    V res = c0 * coeffs.x;
+    res  += c1 * coeffs.y;
+    res  += c2 * coeffs.z;
+    res  += c3 * coeffs.w;
+
+    return res;
+}
+
+template<typename T>
+__device__ inline void Subsample_Bicubic(cudaTextureObject_t src_tex,
+                                         T *dst,
+                                         int dst_width, int dst_height, int dst_pitch,
+                                         int src_width, int src_height,
+                                         int bit_depth)
+{
+    int xo = blockIdx.x * blockDim.x + threadIdx.x;
+    int yo = blockIdx.y * blockDim.y + threadIdx.y;
+
+    if (yo < dst_height && xo < dst_width)
+    {
+        float hscale = (float)src_width / (float)dst_width;
+        float vscale = (float)src_height / (float)dst_height;
+        float xi = (xo + 0.5f) * hscale - 0.5f;
+        float yi = (yo + 0.5f) * vscale - 0.5f;
+        float px = floor(xi);
+        float py = floor(yi);
+        float fx = xi - px;
+        float fy = yi - py;
+
+        float factor = bit_depth > 8 ? 0xFFFF : 0xFF;
+
+        float4 coeffsX = bicubic_coeffs(fx);
+        float4 coeffsY = bicubic_coeffs(fy);
+
+#define PIX(x, y) tex2D<floatT>(src_tex, (x), (y))
+
+        dst[yo * dst_pitch + xo] = from_floatN<T, floatT>(
+            bicubic_filter<floatT>(coeffsY,
+                bicubic_filter<floatT>(coeffsX, PIX(px - 1, py - 1), PIX(px, py - 1), PIX(px + 1, py - 1), PIX(px + 2, py - 1)),
+                bicubic_filter<floatT>(coeffsX, PIX(px - 1, py    ), PIX(px, py    ), PIX(px + 1, py    ), PIX(px + 2, py    )),
+                bicubic_filter<floatT>(coeffsX, PIX(px - 1, py + 1), PIX(px, py + 1), PIX(px + 1, py + 1), PIX(px + 2, py + 1)),
+                bicubic_filter<floatT>(coeffsX, PIX(px - 1, py + 2), PIX(px, py + 2), PIX(px + 1, py + 2), PIX(px + 2, py + 2))
+            ) * factor
+        );
+
+#undef PIX
+    }
+}
+
+/* This does not yield correct results. Most likely because of low internal precision in tex2D linear interpolation */
+template<typename T>
+__device__ inline void Subsample_FastBicubic(cudaTextureObject_t src_tex,
+                                             T *dst,
+                                             int dst_width, int dst_height, int dst_pitch,
+                                             int src_width, int src_height,
+                                             int bit_depth)
+{
+    int xo = blockIdx.x * blockDim.x + threadIdx.x;
+    int yo = blockIdx.y * blockDim.y + threadIdx.y;
+
+    if (yo < dst_height && xo < dst_width)
+    {
+        float hscale = (float)src_width / (float)dst_width;
+        float vscale = (float)src_height / (float)dst_height;
+        float xi = (xo + 0.5f) * hscale - 0.5f;
+        float yi = (yo + 0.5f) * vscale - 0.5f;
+        float px = floor(xi);
+        float py = floor(yi);
+        float fx = xi - px;
+        float fy = yi - py;
+
+        float factor = bit_depth > 8 ? 0xFFFF : 0xFF;
+
+        float h0x, h1x, sx;
+        float h0y, h1y, sy;
+        bicubic_fast_coeffs(fx, &h0x, &h1x, &sx);
+        bicubic_fast_coeffs(fy, &h0y, &h1y, &sy);
+
+#define PIX(x, y) tex2D<floatT>(src_tex, (x), (y))
+
+        floatT pix[4] = {
+            PIX(px + h0x, py + h0y),
+            PIX(px + h1x, py + h0y),
+            PIX(px + h0x, py + h1y),
+            PIX(px + h1x, py + h1y)
+        };
+
+#undef PIX
+
+        dst[yo * dst_pitch + xo] = from_floatN<T, floatT>(
+            lerp_scalar(
+                lerp_scalar(pix[3], pix[2], sx),
+                lerp_scalar(pix[1], pix[0], sx),
+                sy) * factor
+        );
+    }
+}
+
+extern "C" {
+
+#define BICUBIC_KERNEL(T) \
+    __global__ void Subsample_Bicubic_ ## T(cudaTextureObject_t src_tex,                  \
+                                            T *dst,                                       \
+                                            int dst_width, int dst_height, int dst_pitch, \
+                                            int src_width, int src_height,                \
+                                            int bit_depth)                                \
+    {                                                                                     \
+        Subsample_Bicubic<T>(src_tex, dst,                                                \
+                             dst_width, dst_height, dst_pitch,                            \
+                             src_width, src_height,                                       \
+                             bit_depth);                                                  \
+    }
+
+BICUBIC_KERNEL(uchar)
+BICUBIC_KERNEL(uchar2)
+BICUBIC_KERNEL(uchar4)
+
+BICUBIC_KERNEL(ushort)
+BICUBIC_KERNEL(ushort2)
+BICUBIC_KERNEL(ushort4)
+
+}