雙三次插值(BiCubic插值)
雙三次插值(BiCubic插值 )
雙三次插值又稱立方卷積插值。三次卷積插值是一種更加複雜的插值方式。該演算法利用待取樣點周圍16個點的灰度值作三次插值,不僅考慮到4 個直接相鄰點的灰度影響,而且考慮到各鄰點間灰度值變化率的影響。三次運算可以得到更接近高解析度影象的放大效果,但也導致了運算量的急劇增加。這種演算法需要選取插值基函式來擬合數據,其最常用的插值基函式如圖1所示,本次實驗採用如圖所示函式作為基函式。
影象放大並進行BiCubic插值 Matlab/C++程式碼點選開啟連結
假設源影象A大小為m*n,縮放K倍後的目標影象B的大小為M*N,即K=M/m。A的每一個畫素點是已知的,B是未知的,我們想要求出目標影象B中每一畫素點(X,Y)的值,必須先找出畫素(X,Y)在源影象A中對應的畫素(x,y),再根據源影象A距離畫素(x,y)最近
根據比例關係x/X=m/M=1/K,我們可以得到B(X,Y)在A上的對應座標為A(x,y)=A(X*(m/M),Y*(n/N))=A(X/K,Y/K)。如圖所示P點就是目標影象B在(X,Y)處對應於源影象A中的位置,P的座標位置會出現小數部分,所以我們假設 P的座標為P(x+u,y+v),其中x,y分別表示整數部分,u,v分別表示小數部分(藍點到a11方格中紅點的距離)。那麼我們就可以得到如圖所示的最近16個畫素的位置,在這裡用a(i,j)(i,j=0,1,2,3)來表示,如上圖。
我們要做的就是求出BiCubic函式中的引數x,從而獲得上面所說的16個畫素所對應的權重W(x)。BiCubic基函式是一維的,而畫素是二維的,所以我們將畫素點的行與列分開計算。BiCubic函式中的引數x表示該畫素點到P點的距離,例如a00距離P(x+u,y+v)的距離為(1+u,1+v),因此a00的橫座標權重i_0=W(1+u),縱座標權重j_0=W(1+v),a00對B(X,Y)的貢獻值為:(a00畫素值)* i_0* j_0。因此,a0X的橫座標權重分別為W(1+u),W(u),W(1-u),W(2-u);ay0的縱座標權重分別為W(1+v),W(v),W(1-v),W(2-v);B(X,Y)畫素值為:
上述可以用矩陣形式表示,下面圖片來源:點選開啟連結
加權演算法(a可以不取-0.5):
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 <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 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<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;
}
}
}
}
}
}