Translate an image by fractional number of pixels using bilinear interpolation.

translate_bilinear_interp.py

#!/usr/bin/env python

"""
Translate an image by fractional number of pixels using bilinear interpolation.

Based on various snippets of code on StackOverflow.
"""

import numpy as np
import skimage.draw as skdraw

import pycuda.autoinit
import pycuda.gpuarray as gpuarray
from pycuda.compiler import SourceModule

# dest - destination image
# src - source image
# x_shift, y_shift - translation in number of pixels (may be fractional)
# height, width - source image height and width in pixels
translate_bilinear_interp_kernel = """
__global__ void kernel(float *dest, const float *src,
                       const float x_shift, const float y_shift,
                       const int height, const int width) {
    // Integer x and y coordinates of translated pixel to interpolate:
    const unsigned int x = blockIdx.x*blockDim.x+threadIdx.x;
    const unsigned int y = blockIdx.y*blockDim.y+threadIdx.y;

    // Serial index of pixel to interpolate:
    const unsigned int idx = x+y*width;

    // Fractional x and y coordinates of source pixel:
    const float src_x = x-x_shift;
    const float src_y = y-y_shift;

    if (src_x >= 0 && src_x < width && src_y >= 0 && src_y < height) {

        // Integer coordinates of pixel to the left/right of the fractional source pixel:
        int src_x0 = (int)floor(src_x);
        int src_x1 = src_x0+1;
        int src_y0 = (int)floor(src_y);
        int src_y1 = src_y0+1;

        // Serial coordinates of pixels to the left/right/above/below 
        // of the fractional source pixel:
        int idx_src00 = min(max(0, src_x0+src_y0*width), width*height-1);
        int idx_src10 = min(max(0, src_x1+src_y0*width), width*height-1);
        int idx_src01 = min(max(0, src_x0+src_y1*width), width*height-1);
        int idx_src11 = min(max(0, src_x1+src_y1*width), width*height-1);

        // Interpolate source pixel:
        float sx = src_x-src_x0;
        float sy = src_y-src_y0;
        dest[idx]  = (1.0f-sx)*(1.0f-sy)*src[idx_src00];
        dest[idx] += (     sx)*(1.0f-sy)*src[idx_src10];
        dest[idx] += (1.0f-sx)*(     sy)*src[idx_src01];
        dest[idx] += (     sx)*(     sy)*src[idx_src11];
    } else {
        // Zero out regions shifted from outside of the image bounds:
        dest[idx] = 0.0f;
    }
}
"""

# Image to translate:
width = height = 200
src = np.zeros((height, width), np.float32)
src[skdraw.circle(height/2, width/2, 30, (height, width))] = 1.0
src[skdraw.circle(height/2, width/2, 15, (height, width))] = 0.0

# Desired translation:
x_shift = -5.3
y_shift = 7.25

# Perform translation:
src_gpu = gpuarray.to_gpu(src)
dest_gpu = gpuarray.zeros_like(src_gpu)

mod = SourceModule(translate_bilinear_interp_kernel)
kernel = mod.get_function('kernel')
block = (32, 32, 1)
grid = (int(np.ceil(float(height)/block[0])), int(np.ceil(float(width)/block[1])))
kernel(dest_gpu, src_gpu, np.float32(x_shift), np.float32(y_shift), np.int32(height), np.int32(width),
       block=block, grid=grid)

dest = dest_gpu.get()