双三次插值(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)最近的16个像素点作为计算目标图像B(X,Y)处像素值的参数,利用BiCubic基函数求出16个像素点的权重,图B像素(x,y)的值就等于16个像素点的加权叠加。

 

 根据比例关系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;
						}
				}
			}
		}
	}
}

【转载】https://blog.csdn.net/qq_34885184/article/details/79163991

              https://blog.csdn.net/u010555688/article/details/24352343

已标记关键词 清除标记
相关推荐
©️2020 CSDN 皮肤主题: 技术黑板 设计师:CSDN官方博客 返回首页