BiCubic插值原理：

双三次插值又称立方卷积插值。三次卷积插值是一种更加复杂的插值方式。该算法利用待采样点周围16个点的灰度值作三次插值，不仅考虑到4 个直接相邻点的灰度影响，而且考虑到各邻点间灰度值变化率的影响。三次运算可以得到更接近高分辨率图像的放大效果，但也导致了运算量的急剧增加。这种算法需要选取插值基函数来拟合数据，其最常用的插值基函数如图1所示，本次实验采用如图所示函数作为基函数。

构造BiCubic函数：

其中，a取-0.5.

BiCubic函数具有如下形状：

[source:  R. Keys, (1981). "Cubic convolution interpolation for digital image processing". IEEE Transactions on Signal Processing, Acoustics, Speech, and Signal Processing 29 (6): 1153–1160.]

对待插值的像素点(x,y)(x和y可以为浮点数)，取其附近的4x4邻域点(xi,yj), i,j = 0,1,2,3。按如下公式进行插值计算：

Matlab实现代码：

%双三次插值具体实现

clc,clear;

fff=imread('E:\Documents\BUPT\DIP\图片\lena.bmp');

ff =rgb2gray(fff);%转化为灰度图像

[mm,nn]=size(ff); %将图像隔行隔列抽取元素，得到缩小的图像f

m=mm/2;

n=nn/2;

f =zeros(m,n);

for i=1:m

for j=1:n

f(i,j)=ff(2*i,2*j);

end

end

k=5; %设置放大倍数

bijiao1 =imresize(f,k,'bilinear');%双线性插值结果比较

bijiao =uint8(bijiao1);

a=f(1,:);

c=f(m,:); %将待插值图像矩阵前后各扩展两行两列,共扩展四行四列

b=[f(1,1),f(1,1),f(:,1)',f(m,1),f(m,1)];

d=[f(1,n),f(1,n),f(:,n)',f(m,n),f(m,n)];

a1=[a;a;f;c;c];

b1=[b;b;a1';d;d];

ffff=b1';

f1=double(ffff);

g1 =zeros(k*m,k*n);

fori=1:k*m %利用双三次插值公式对新图象所有像素赋值

u=rem(i,k)/k;

i1=floor(i/k)+2;

A=[sw(1+u) sw(u) sw(1-u) sw(2-u)];

for j=1:k*n

v=rem(j,k)/k;

j1=floor(j/k)+2;

C=[sw(1+v);sw(v);sw(1-v);sw(2-v)];

B=[f1(i1-1,j1-1) f1(i1-1,j1) f1(i1-1,j1+1)f1(i1-1,j1+2)

f1(i1,j1-1) f1(i1,j1) f1(i1,j1+1) f1(i1,j1+2)

f1(i1+1,j1-1) f1(i1+1,j1) f1(i1+1,j1+1) f1(i1+1,j1+2)

f1(i1+2,j1-1) f1(i1+2,j1) f1(i1+2,j1+1)f1(i1+2,j1+2)];

g1(i,j)=(A*B*C);

end

end

g=uint8(g1);

imshow(uint8(f));title('缩小的图像'); %显示缩小的图像

figure,imshow(ff);title('原图'); %显示原图像

figure,imshow(g);title('双三次插值放大的图像'); %显示插值后的图像

figure,imshow(bijiao);title('双线性插值放大结果'); %显示插值后的图像

mse=0;

ff=double(ff);

g=double(g);

ff2=fftshift(fft2(ff)); %计算原图像和插值图像的傅立叶幅度谱

g2=fftshift(fft2(g));

figure,subplot(1,2,1),imshow(log(abs(ff2)),[8,10]);title('原图像的傅立叶幅度谱');

subplot(1,2,2),imshow(log(abs(g2)),[8,10]);title('双三次插值图像的傅立叶幅度谱');

基函数代码：

functionA=sw(w1)

w=abs(w1);

ifw<1&&w>=0

A=1-2*w^2+w^3;

elseifw>=1&&w<2

A=4-8*w+5*w^2-w^3;

else

A=0;

end

C++实现代码：

#include "opencv2/imgproc/imgproc.hpp"

#include "opencv2/highgui/highgui.hpp"

#include

#include

#include

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 image

char* imageName = "images/Lenna_256.png";

Mat image;

image = imread(imageName,1);if(!image.data)

{

cout << "No image data" << endl;

return -1;

}

// show image

namedWindow("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;

}

else

return 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 image

float 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 image

dst.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(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(y0);

q1 = src.ptr(y1);

q2 = src.ptr(y2);

q3 = src.ptr(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;

}

}

}

}

}

}