阅读 133

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

概述

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

数学原理

  1. 基础知识

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

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

f(x)=ax^3+bx^2+cx+d

对这个函数求导得

f'(x)=3ax^2+2bx+c

那么

f(0)=d
f(1)=a+b+c+d
f'(0)=c
f'(1)=3a+2b+c

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

  1. 近似求导数

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

y=ax+b

得斜率a为

a=\frac{y_{1}-y_{2}}{x_{1}-x_{2}}

因为两点x轴距离为2

x_{1}-x_{2}=2

f'(0)= \frac{p_{2}-p_{0}}{2}
f'(1)= \frac{p_{3}-p_{1}}{2}
f(0)=p_{1}
f(1)=p_{2}

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

  1. 扩展到二维上

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

p(x,y)=\sum_{i=0}^{3}\sum_{j=0}^{3}a_{ij}x^iy^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);
复制代码
  1. 计算宽、高的缩放常数
float rowRatio = ((float)sourceHeight) / ((float)desHeight);
float colRatio = ((float)sourceWidth) / ((float)desWidth);
复制代码
  1. 根据缩放常数,找到当前遍历位置所在的源图的位置,这个位置是有小数的,这里我们直接取整数部分,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;
复制代码
  1. 这里直接把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];
}
复制代码
  1. 通过坐标和上面计算出来的多项式系数,计算像素值。这里的坐标就是像素位置在上面中间灰色区块中的坐标位置,范围在[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

参考链接

关注下面的标签,发现更多相似文章
评论