前言
图像处理中有三种常用的插值算法：
最邻近插值
双线性插值
双立方（三次卷积）插值
其中效果最好的是双立方（三次卷积）插值，本文介绍它的原理以及使用
如果想先看效果和源码，可以拉到最底部
本文的契机是某次基于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;}}}}}
}
