图像重采样算法之bicubic双立方插值算法

概述

最邻近插值算法的目标像素值由源图上单个像素决定，双线性插值算法由源像素某点周围4个像素点按一定权重获得，而双立方插值算法更进一步参考了源像素某点周围4*4个像素来获得。

数学原理

1. 基础知识

如果已知一个函数f(x)以及它在x=0,x=1处的导数，那么函数可以在[0,1]之间插值，这个约束条件很重要。当函数表达为三次多项式时我们称之为立方插值。

我们用三次多项式来表示一条曲线

对这个函数求导得

那么

因此我们可以算出a，b，c，d，也就得到这个f(x)函数，代入我们在[0,1]区间的点，就可以求出函数值。

2. 近似求导数

假设我们有四个值，分别p0，p1，p2，p3代表x=-1，x=0，x=1和x=2位置的函数值，我们需要求x在[0,1]范围内的某个值对应的y值。但是我们仅有这四个点，不知道函数在0和1处的导数。怎么办呢？为了得到平滑的曲线，我们采用近似法：我们通过当前点相邻的间隔的两点连线，用这条直线的斜率来近似替代这个位置导数。

得斜率a为

因为两点x轴距离为2

得

从而可以算出a，b，c，d值

3. 扩展到二维上

假设我们有16个点,i,j分别表示行位置和列位置，我们可以将行方向和列方向分解，比如先在行方向做三次二项式

从4x4输入像素值进行双立方插值以计算一个输出像素值。

1. 输出像素必须位于相对于输入像素（黑点）的灰色区域。
2. 通过从四个输入像素的每一行进行水平插值来独立计算四个中间值（红色点）。
3. 使用相同的插值方法，从四个中间值的列中垂直插值输出值（蓝点）。

每个插值总计为四个值的加权和，权重由灰色正方形内水平或垂直位置的三次多项式给出。双线性插值的工作原理类似，但是仅以灰色正方形的四个角的像素（2x2像素）作为输入，并且多项式是线性的。

关键代码

我在iOS平台上，使用Objective-C语言实现这个算法。首先需要获取到图片的RGBA数据，也叫RAW数据，这是一个一维数组，但是在实际处理中我们需要以二维的思路来遍历。

1. 获取图片RGBA数据

UInt32* pixelData = (UInt32 *)calloc(width * height, sizeof(UInt32));
CGColorSpaceRef colorSpace = CGColorSpaceCreateDeviceRGB();
CGContextRef context = CGBitmapContextCreate(pixelData,
                                                width,
                                                height,
                                                bitsPerComponent,
                                                bytesPerRow,
                                                colorSpace,
                                                kCGImageAlphaPremultipliedLast | kCGBitmapByteOrder32Big);

CGContextDrawImage(context, CGRectMake(0, 0, width, height), cgOriginalImage);

2. 计算宽、高的缩放常数

float rowRatio = ((float)sourceHeight) / ((float)desHeight);
float colRatio = ((float)sourceWidth) / ((float)desWidth);

3. 根据缩放常数，找到当前遍历位置所在的源图的位置，这个位置是有小数的，这里我们直接取整数部分，j代表在原图的行，k代表在原图的列，小数部分就是我们要求的u和t，在灰色区块中的位置。

double srcRow = ((float)row) * rowRatio;
double j = floor(srcRow);
double u = srcRow - j;

double srcCol = ((float)col) * colRatio;
double k = floor(srcCol);
double t = srcCol - k;

4. 这里直接把16个值代入公式，直接求得所有的多项式系数

static void updateCoefficients(UInt32 p[4][4]) {
    a00 = p[1][1];
    a01 = -.5*p[1][0] + .5*p[1][2];
    a02 = p[1][0] - 2.5*p[1][1] + 2*p[1][2] - .5*p[1][3];
    a03 = -.5*p[1][0] + 1.5*p[1][1] - 1.5*p[1][2] + .5*p[1][3];
    a10 = -.5*p[0][1] + .5*p[2][1];
    a11 = .25*p[0][0] - .25*p[0][2] - .25*p[2][0] + .25*p[2][2];
    a12 = -.5*p[0][0] + 1.25*p[0][1] - p[0][2] + .25*p[0][3] + .5*p[2][0] - 1.25*p[2][1] + p[2][2] - .25*p[2][3];
    a13 = .25*p[0][0] - .75*p[0][1] + .75*p[0][2] - .25*p[0][3] - .25*p[2][0] + .75*p[2][1] - .75*p[2][2] + .25*p[2][3];
    a20 = p[0][1] - 2.5*p[1][1] + 2*p[2][1] - .5*p[3][1];
    a21 = -.5*p[0][0] + .5*p[0][2] + 1.25*p[1][0] - 1.25*p[1][2] - p[2][0] + p[2][2] + .25*p[3][0] - .25*p[3][2];
    a22 = p[0][0] - 2.5*p[0][1] + 2*p[0][2] - .5*p[0][3] - 2.5*p[1][0] + 6.25*p[1][1] - 5*p[1][2] + 1.25*p[1][3] + 2*p[2][0] - 5*p[2][1] + 4*p[2][2] - p[2][3] - .5*p[3][0] + 1.25*p[3][1] - p[3][2] + .25*p[3][3];
    a23 = -.5*p[0][0] + 1.5*p[0][1] - 1.5*p[0][2] + .5*p[0][3] + 1.25*p[1][0] - 3.75*p[1][1] + 3.75*p[1][2] - 1.25*p[1][3] - p[2][0] + 3*p[2][1] - 3*p[2][2] + p[2][3] + .25*p[3][0] - .75*p[3][1] + .75*p[3][2] - .25*p[3][3];
    a30 = -.5*p[0][1] + 1.5*p[1][1] - 1.5*p[2][1] + .5*p[3][1];
    a31 = .25*p[0][0] - .25*p[0][2] - .75*p[1][0] + .75*p[1][2] + .75*p[2][0] - .75*p[2][2] - .25*p[3][0] + .25*p[3][2];
    a32 = -.5*p[0][0] + 1.25*p[0][1] - p[0][2] + .25*p[0][3] + 1.5*p[1][0] - 3.75*p[1][1] + 3*p[1][2] - .75*p[1][3] - 1.5*p[2][0] + 3.75*p[2][1] - 3*p[2][2] + .75*p[2][3] + .5*p[3][0] - 1.25*p[3][1] + p[3][2] - .25*p[3][3];
    a33 = .25*p[0][0] - .75*p[0][1] + .75*p[0][2] - .25*p[0][3] - .75*p[1][0] + 2.25*p[1][1] - 2.25*p[1][2] + .75*p[1][3] + .75*p[2][0] - 2.25*p[2][1] + 2.25*p[2][2] - .75*p[2][3] - .25*p[3][0] + .75*p[3][1] - .75*p[3][2] + .25*p[3][3];
}

5. 通过坐标和上面计算出来的多项式系数，计算像素值。这里的坐标就是像素位置在上面中间灰色区块中的坐标位置，范围在[0,1]之间，体现在上面计算的t和u值。

static double getValue(double x, double y) {
    double x2 = x * x;
    double x3 = x2 * x;
    double y2 = y * y;
    double y3 = y2 * y;

    return (a00 + a01 * y + a02 * y2 + a03 * y3) +
           (a10 + a11 * y + a12 * y2 + a13 * y3) * x +
           (a20 + a21 * y + a22 * y2 + a23 * y3) * x2 +
           (a30 + a31 * y + a32 * y2 + a33 * y3) * x3;
}

最终效果

Source Code

GitHub

参考链接

• 图像处理（一）bicubic解释推导
• 图像放缩之双立方插值
• How does nearest neighbour, bilinear and cubic interpolation work in images?
• Image Processing in iOS Part 1: Raw Bitmap Modification
• UIImage resize with hard edges
• Capturing Uncompressed Image Data
• Image Processing in C
• Converting RGB data into a bitmap in Objective-C++ Cocoa
• Converting Array of Pixel Brightness Values Back into UIImage?
• 常见的图像重采样算法