前言

图像处理中有三种常用的插值算法：

---

最邻近插值

双线性插值

双立方（三次卷积）插值

---

其中效果最好的是双立方（三次卷积）插值，本文介绍它的原理以及使用

如果想先看效果和源码，可以拉到最底部

本文的契机是某次基于canvas做图像处理时，发现canvas自带的缩放功能不尽人意，于是重温了下几种图像插值算法，并整理出来。

为何要进行双立方插值

---

• 对图像进行插值的目的是为了获取缩小或放大后的图片
• 常用的插值算法中，双立方插值效果最好
• 本文中介绍双立方插值的一些数学理论以及实现

---

双立方和三次卷积只是这个插值算法的两种不同叫法而已，可以自行推导，会发现最终可以将求值转化为卷积公式

另外，像Photoshop等图像处理软件中也有这三种算法的实现

数学理论

双立方插值计算涉及到16个像素点，如下图

简单分析如下：

• 其中P00代表目标插值图中的某像素点(x, y)在原图中最接近的映射点

• 譬如映射到原图中的坐标为(1.1, 1.1)，那么P00就是(1, 1)

• 而最终插值后的图像中的(x, y)处的值即为以上16个像素点的权重卷积之和

---

下图进一步分析

如下是对图的一些简单分析

• 譬如计算插值图中(distI, distJ)处像素的值

• 首先计算它映射到原图中的坐标(i + v, j + u)

• 也就是说，卷积计算时，p00点对应(i, j)坐标

• 最终，插值后的图中(distI, distJ)坐标点对应的值是原图中(i, j)处邻近16个像素点的权重卷积之和

• i, j的范围是[i - 1, i + 2]，[j - 1, j + 2]

卷积公式

• 设采样公式为S(x)
• 原图中每一个(i, j)坐标点的值得表达式为f(i, j)
• 插值后对应坐标的值为F(i + v, j + u)（这个值会作为(distI, distJ)坐标点的值）
  那公式为
  
  

采样公式

在卷积公式中有一个S(x)，它就是关键的卷积插值公式

不同的公式，插值效果会有所差异（会导致加权值不一样）

本文中采用WIKI-Bicubic interpolation中给出的插值公式：

公式中的特点是：

• S(0) = 1

• S(n) = 0(当n为整数时)

• 当x超出范围时，S(x)为0

• 当a取不同值时可以用来逼近不同的样条函数（常用值-0.5, -0.75）

当a取值为-1
公式如下：

此时，逼近的函数是y = sin(xPI)/(xPI)，如图

当a取值为-0.5
公式如下：

此时对应三次Hermite样条

不同a的简单对比

代码实现

 
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"
#include <iostream>
#include <cmath>
#include <fstream>
using namespace cv;
using namespace std;
#define PI 3.14159265
float BiCubicPoly(float x);
void MyScaleBiCubicInter(Mat& src, Mat& dst, float TransMat[3][3]);
/*** @function main*/
int main( int argc, char** argv )
{// load imagechar* imageName = "images/Lenna_256.png";Mat image;image = imread(imageName,1);if(!image.data){cout << "No image data" << endl;return -1;}// show imagenamedWindow("image", CV_WINDOW_AUTOSIZE);imshow("image", image);  Mat dst;float transMat[3][3] = { {2.0, 0, 0}, {0, 2.0, 0}, {0, 0, 1} };MyScaleBiCubicInter(image, dst, transMat);namedWindow("out_image", CV_WINDOW_AUTOSIZE);imshow("out_image", dst);imwrite("Lenna_scale_biCubic2.jpg", dst);waitKey(0);return 0;
}
float BiCubicPoly(float x)
{float abs_x = abs(x);float a = -0.5;if( abs_x <= 1.0 ){return (a+2)*pow(abs_x,3) - (a+3)*pow(abs_x,2) + 1;}else if( abs_x < 2.0 ){return a*pow(abs_x,3) - 5*a*pow(abs_x,2) + 8*a*abs_x - 4*a;}elsereturn 0.0;
}void MyScaleBiCubicInter(Mat& src, Mat& dst, float TransMat[3][3])
{CV_Assert(src.data);CV_Assert(src.depth() != sizeof(uchar));// calculate margin point of dst imagefloat left =  0;float right =  0;float top =  0;float down =  0;float x = src.cols * 1.0f;float y = 0.0f;float u1 = x * TransMat[0][0] + y * TransMat[0][1];float v1 = x * TransMat[1][0] + y * TransMat[1][1];x = src.cols * 1.0f;y = src.rows * 1.0f;float u2 = x * TransMat[0][0] + y * TransMat[0][1];float v2 = x * TransMat[1][0] + y * TransMat[1][1];x = 0.0f;y = src.rows * 1.0f;float u3 = x * TransMat[0][0] + y * TransMat[0][1];float v3 = x * TransMat[1][0] + y * TransMat[1][1];left =  min( min( min(0.0f,u1), u2 ), u3);right =  max( max( max(0.0f,u1), u2 ), u3);top =  min( min( min(0.0f,v1), v2 ), v3);down =  max( max( max(0.0f,v1), v2 ), v3);// create dst imagedst.create(int(abs(right-left)), int(abs(down-top)), src.type());	CV_Assert( dst.channels() == src.channels() );int channels = dst.channels();int i,j;uchar* p;uchar* q0;uchar* q1;uchar* q2;uchar* q3;for( i = 0; i < dst.rows; ++i){p = dst.ptr<uchar>(i);for ( j = 0; j < dst.cols; ++j){// x = (j+left)/TransMat[0][0]  ; y = (i+top)/TransMat[1][1] ;int x0 = int(x) - 1;int y0 = int(y) - 1;int x1 = int(x);int y1 = int(y);int x2 = int(x) + 1;int y2 = int(y) + 1;int x3 = int(x) + 2;int y3 = int(y) + 2;if( (x0 >= 0) && (x3 < src.cols) && (y0 >= 0) && (y3 < src.rows) ) {q0 = src.ptr<uchar>(y0);q1 = src.ptr<uchar>(y1);q2 = src.ptr<uchar>(y2);q3 = src.ptr<uchar>(y3);float dist_x0 = BiCubicPoly(x-x0);float dist_x1 = BiCubicPoly(x-x1);float dist_x2 = BiCubicPoly(x-x2);float dist_x3 = BiCubicPoly(x-x3);float dist_y0 = BiCubicPoly(y-y0);float dist_y1 = BiCubicPoly(y-y1);float dist_y2 = BiCubicPoly(y-y2);float dist_y3 = BiCubicPoly(y-y3);float dist_x0y0 = dist_x0 * dist_y0;float dist_x0y1 = dist_x0 * dist_y1;float dist_x0y2 = dist_x0 * dist_y2;float dist_x0y3 = dist_x0 * dist_y3;float dist_x1y0 = dist_x1 * dist_y0;float dist_x1y1 = dist_x1 * dist_y1;float dist_x1y2 = dist_x1 * dist_y2;float dist_x1y3 = dist_x1 * dist_y3;float dist_x2y0 = dist_x2 * dist_y0;float dist_x2y1 = dist_x2 * dist_y1;float dist_x2y2 = dist_x2 * dist_y2;float dist_x2y3 = dist_x2 * dist_y3;float dist_x3y0 = dist_x3 * dist_y0;float dist_x3y1 = dist_x3 * dist_y1;float dist_x3y2 = dist_x3 * dist_y2;float dist_x3y3 = dist_x3 * dist_y3;switch(channels){case 1:{break;}case 3:{p[3*j] =    (uchar)(q0[3*x0] * dist_x0y0 +q1[3*x0] * dist_x0y1 +q2[3*x0] * dist_x0y2 +q3[3*x0] * dist_x0y3 +q0[3*x1] * dist_x1y0 +q1[3*x1] * dist_x1y1 +q2[3*x1] * dist_x1y2 +q3[3*x1] * dist_x1y3 +q0[3*x2] * dist_x2y0 +q1[3*x2] * dist_x2y1 +q2[3*x2] * dist_x2y2 +q3[3*x2] * dist_x2y3 +q0[3*x3] * dist_x3y0 +q1[3*x3] * dist_x3y1 +q2[3*x3] * dist_x3y2 +q3[3*x3] * dist_x3y3 ) ;p[3*j+1] =  (uchar)(q0[3*x0+1] * dist_x0y0 +q1[3*x0+1] * dist_x0y1 +q2[3*x0+1] * dist_x0y2 +q3[3*x0+1] * dist_x0y3 +q0[3*x1+1] * dist_x1y0 +q1[3*x1+1] * dist_x1y1 +q2[3*x1+1] * dist_x1y2 +q3[3*x1+1] * dist_x1y3 +q0[3*x2+1] * dist_x2y0 +q1[3*x2+1] * dist_x2y1 +q2[3*x2+1] * dist_x2y2 +q3[3*x2+1] * dist_x2y3 +q0[3*x3+1] * dist_x3y0 +q1[3*x3+1] * dist_x3y1 +q2[3*x3+1] * dist_x3y2 +q3[3*x3+1] * dist_x3y3 ) ;p[3*j+2] =  (uchar)(q0[3*x0+2] * dist_x0y0 +q1[3*x0+2] * dist_x0y1 +q2[3*x0+2] * dist_x0y2 +q3[3*x0+2] * dist_x0y3 +q0[3*x1+2] * dist_x1y0 +q1[3*x1+2] * dist_x1y1 +q2[3*x1+2] * dist_x1y2 +q3[3*x1+2] * dist_x1y3 +q0[3*x2+2] * dist_x2y0 +q1[3*x2+2] * dist_x2y1 +q2[3*x2+2] * dist_x2y2 +q3[3*x2+2] * dist_x2y3 +q0[3*x3+2] * dist_x3y0 +q1[3*x3+2] * dist_x3y1 +q2[3*x3+2] * dist_x3y2 +q3[3*x3+2] * dist_x3y3 ) ;float thre = 198.0f;if( (abs(p[3*j]-q1[3*x1]) > thre) || (abs(p[3*j+1]-q1[3*x1+1]) > thre) ||(abs(p[3*j+2]-q1[3*x1+2]) > thre) ){p[3*j] = q1[3*x1];p[3*j+1] = q1[3*x1+1];p[3*j+2] = q1[3*x1+2];}					break;}}}}}
}