双三次插值 - 插值图像任意位置亚像素C++

双三次插值 - 插值图像任意位置亚像素C++

一、概念

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

加权算法（a可以不取-0.5）：点开链接

根据比例关系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)像素值为：

二、代码实现C++(可配合读取图片程序使用)

#include <iostream>
#include <tchar.h>
#include "ImageBuffer.h"

using namespace CamCtrl;

void test()
{
	CImageBuffer mCImageBuffer;
	mCImageBuffer.LoadFromBmp(_T("1.bmp"));
}

int main()
{
	test();
	return 0;
}

//功能：双三边插值
//image：输入图片
//iImgWidth：图像宽度
//iImgHeight：图像高度
//dx，dy：亚像素点与整像素点posy在下x，y上的距离
//ipos = (posy - 1)*(iImgWidth)+(posx - 1);  //得到ipos
double CubicPolynomial(unsigned char * image, int iImgWidth, int iImgHeight, double dx, double dy, int ipos)
{
	// image为图像各个位置的像素值， iImgWidth为图像的宽度, iImgHeight为图像的高度, ipos为当前插值的位置。

	// 先计算padding 图像的大小。
	int padding_size = (iImgHeight + 2) * (iImgWidth + 2);
	double * padding_img = new double[padding_size];

	// 接着对 padding_img 进行赋值。
	// 将原图像的值赋给padding_img的(2,2)-(iImgHeight+1, iImgWidth+1)
	// 该部分验证无误。
	for (int i = 0; i < iImgHeight; i++)
	{
		for (int j = 0; j < iImgWidth; j++)
		{
			int ipos1 = (i + 1) * (iImgWidth + 2) + (j + 1);

			int temp = i * iImgWidth + j;

			padding_img[ipos1] = image[temp];
		}
	}

	// 第一行的值： 第二行、第三行、第四行的值的线性组合.
	// 该部分验证无误。

	for (int i = 1; i < iImgWidth + 1; i++)
	{
		int ipos1 = i;

		double temp1 = padding_img[ipos1 + iImgWidth + 2];
		double temp2 = padding_img[ipos1 + 2 * (iImgWidth + 2)];
		double temp3 = padding_img[ipos1 + 3 * (iImgWidth + 2)];

		padding_img[ipos1] = 3 * temp1 - 3 * temp2 + temp3;
	}

	// 最后一行的值：倒数第二、第三、第四行的值的线性组合。
	// 该部分验证无误。

	for (int i = 1; i < iImgWidth + 1; i++)
	{
		int ipos1 = i + (iImgHeight + 1)*(iImgWidth + 2);

		double temp1 = padding_img[ipos1 - (iImgWidth + 2)];
		double temp2 = padding_img[ipos1 - 2 * (iImgWidth + 2)];
		double temp3 = padding_img[ipos1 - 3 * (iImgWidth + 2)];

		padding_img[ipos1] = 3 * temp1 - 3 * temp2 + temp3;
	}

	// 第一列的值：第二列、第三列、第四列的值的线性组合。
	// 该部分验证无误。
	for (int i = 1; i < iImgHeight + 1; i++)
	{
		int ipos1 = i * (iImgWidth + 2);

		double temp1 = padding_img[ipos1 + 1];
		double temp2 = padding_img[ipos1 + 2];
		double temp3 = padding_img[ipos1 + 3];

		padding_img[ipos1] = 3 * temp1 - 3 * temp2 + temp3;
	}

	// 最后一列的值：倒数第二、第三、第四列的值的线性组合。
	// 该部分验证无误。
	for (int i = 1; i < iImgHeight + 1; i++)
	{
		int ipos1 = (i + 1)*(iImgWidth + 2) - 1;

		double temp1 = padding_img[ipos1 - 1];
		double temp2 = padding_img[ipos1 - 2];
		double temp3 = padding_img[ipos1 - 3];

		padding_img[ipos1] = 3 * temp1 - 3 * temp2 + temp3;
	}

	// 后面4个点的组合方向均为列的方向。即ipos+1, ipos+2, ipos+3，或ipos-1, ipos-2, ipos-3.

	// 左上角的点。
	padding_img[0] = 3 * padding_img[1] - 3 * padding_img[2] + padding_img[3];

	// 右上角的点。
	padding_img[iImgWidth + 1] = 3 * padding_img[iImgWidth] - 3 * padding_img[iImgWidth - 1] + padding_img[iImgWidth - 2];

	// 左下角的点。

	int pos_left_bottom = (iImgHeight + 1)*(iImgWidth + 2);
	padding_img[pos_left_bottom] = 3 * padding_img[pos_left_bottom + 1] - 3 * padding_img[pos_left_bottom + 2] + padding_img[pos_left_bottom + 3];

	// 右下角的点。

	int pos_right_bottom = (iImgHeight + 2)*(iImgWidth + 2) - 1;
	padding_img[pos_right_bottom] = 3 * padding_img[pos_right_bottom - 1] - 3 * padding_img[pos_right_bottom - 2] + padding_img[pos_right_bottom - 3];


	// padding结束之后，进行插值操作。首先要把原图的ipos对应到新图的ipos。

	int height = (ipos + 1) / iImgWidth;

	int width = (ipos + 1) - height * iImgWidth;

	// 在新图中的位置。

	int new_ipos = width + (height + 1) * (iImgWidth + 2);


	// x方向卷积
	double dx_2 = dx * dx;
	double dx_3 = dx * dx*dx;
	double a4 = (-dx_3 + 2 * dx_2 - dx) / 2;
	double a5 = (3 * dx_3 - 5 * dx_2 + 2) / 2;
	double a6 = (-3 * dx_3 + 4 * dx_2 + dx) / 2;
	double a7 = (dx_3 - dx_2) / 2;



	double a8 = a4 * padding_img[new_ipos - 1 - (iImgWidth + 2)] + a5 * padding_img[new_ipos - (iImgWidth + 2)] + a6 * padding_img[new_ipos + 1 - (iImgWidth + 2)] + a7 * padding_img[new_ipos + 2 - (iImgWidth + 2)];
	double a9 = a4 * padding_img[new_ipos - 1] + a5 * padding_img[new_ipos] + a6 * padding_img[new_ipos + 1] + a7 * padding_img[new_ipos + 2];
	double a10 = a4 * padding_img[new_ipos - 1 + (iImgWidth + 2)] + a5 * padding_img[new_ipos + (iImgWidth + 2)] + a6 * padding_img[new_ipos + 1 + (iImgWidth + 2)] + a7 * padding_img[new_ipos + 2 + (iImgWidth + 2)];
	double a11 = a4 * padding_img[new_ipos - 1 + 2 * (iImgWidth + 2)] + a5 * padding_img[new_ipos + 2 * (iImgWidth + 2)] + a6 * padding_img[new_ipos + 1 + 2 * (iImgWidth + 2)] + a7 * padding_img[new_ipos + 2 + 2 * (iImgWidth + 2)];



	// y方向卷积
	double dy_2 = dy * dy;
	double dy_3 = dy * dy*dy;
	double a12 = (-dy_3 + 2 * dy_2 - dy) / 2;
	double a13 = (3 * dy_3 - 5 * dy_2 + 2) / 2;
	double a14 = (-3 * dy_3 + 4 * dy_2 + dy) / 2;
	double a15 = (dy_3 - dy_2) / 2;

	double gray_value;
	gray_value = a8 * a12 + a9 * a13 + a10 * a14 + a11 * a15;


	delete padding_img;

	return gray_value;


}

三、C++代码（基于opencv）

#include <vector>
#include "opencv2/core.hpp"
#include "opencv2/imgproc.hpp"
#include "opencv2/imgcodecs.hpp"
#include <opencv2/opencv.hpp>
#include "bmptools.h"
#include <cstdlib>
using namespace std;

using namespace cv;



double CubicPolynomial(cv::Mat pGray, double dx, double dy, int posx, int posy)
{	
	int rows = pGray.rows;
	int cols = pGray.cols;

	// 创建padding一圈的图像，用来做卷积
	cv::Mat padded_image = cv::Mat::zeros(rows + 2, cols + 2, CV_64FC1);

	// 将原始图像放置padding图像的中间
	for (int i = 1; i < padded_image.rows-1; i++)
	{	
		for (int j = 1; j < padded_image.cols-1; j++)
		{
			double temp = pGray.at<uchar>(i - 1, j - 1);
			padded_image.at<double>(i, j) = temp;
		}
	}

	// 分别计算周围padding一圈的像素值
	for (int j = 1; j < padded_image.cols-1; j++)
	{
		// 第一行
		padded_image.at<double>(0, j) = 3 * padded_image.at<double>(1, j) - 3 * padded_image.at<double>(2, j) + padded_image.at<double>(3, j);
		//最后一行
		padded_image.at<double>(rows + 1, j) = 3 * padded_image.at<double>(rows, j) - 3 * padded_image.at<double>(rows - 1, j) + padded_image.at<double>(rows - 2, j);
	}

	for (int i = 1; i < padded_image.rows-1; i++)
	{
		// 第一列
		padded_image.at<double>(i, 0) = 3 * padded_image.at<double>(i, 1) - 3 * padded_image.at<double>(i, 2) + padded_image.at<double>(i, 3);
		// 最后一列
		padded_image.at<double>(i, cols + 1) = 3 * padded_image.at<double>(i, cols) - 3 * padded_image.at<double>(i, cols - 1) + padded_image.at<double>(i, cols - 2);
	}


	padded_image.at<double>(0, 0) = 3 * padded_image.at<double>(0, 1) - 3 * padded_image.at<double>(0, 2) + padded_image.at<double>(0, 3);

	padded_image.at<double>(0, cols + 1) = 3 * padded_image.at<double>(0, cols) - 3 * padded_image.at<double>(0, cols - 1) + padded_image.at<double>(0, cols - 2);

	padded_image.at<double>(rows + 1, 0) = 3 * padded_image.at<double>(rows + 1, 1) - 3 * padded_image.at<double>(rows + 1, 2) + padded_image.at<double>(rows + 1, 3);

	padded_image.at<double>(rows + 1, cols + 1) = 3 * padded_image.at<double>(rows + 1, cols) - 3 * padded_image.at<double>(rows + 1, cols - 1) + padded_image.at<double>(rows + 1, cols - 2);

	// 卷积的过程

	// x方向卷积
	double dx_2 = dx*dx;
	double dx_3 = dx*dx*dx;
	double a4 = (-dx_3 + 2 * dx_2 - dx) / 2;
	double a5 = (3 * dx_3 - 5 * dx_2 + 2) / 2;
	double a6 = (-3 * dx_3 + 4 * dx_2 + dx) / 2;
	double a7 = (dx_3 - dx_2) / 2;

	double a8 = a4*padded_image.at<double>(posy - 1, posx - 1) + a5*padded_image.at<double>(posy - 1, posx) + a6*padded_image.at<double>(posy - 1, posx + 1) + a7*padded_image.at<double>(posy - 1, posx + 2);
	double a9 = a4*padded_image.at<double>(posy, posx - 1) + a5*padded_image.at<double>(posy, posx) + a6*padded_image.at<double>(posy, posx + 1) + a7*padded_image.at<double>(posy, posx + 2);
	double a10 = a4*padded_image.at<double>(posy + 1, posx - 1) + a5*padded_image.at<double>(posy + 1, posx) + a6*padded_image.at<double>(posy + 1, posx + 1) + a7*padded_image.at<double>(posy + 1, posx + 2);
	double a11 = a4*padded_image.at<double>(posy + 2, posx - 1) + a5*padded_image.at<double>(posy + 2, posx) + a6*padded_image.at<double>(posy + 2, posx + 1) + a7*padded_image.at<double>(posy + 2, posx + 2);

	// y方向卷积
	double dy_2 = dy*dy;
	double dy_3 = dy*dy*dy;
	double a12 = (-dy_3 + 2 * dy_2 - dy) / 2;
	double a13 = (3 * dy_3 - 5 * dy_2 + 2) / 2;
	double a14 = (-3 * dy_3 + 4 * dy_2 + dy) / 2;
	double a15 = (dy_3 - dy_2) / 2;

	double gray_value;
	gray_value = a8*a12 + a9*a13 + a10*a14 + a11*a15;

	return gray_value;
}


// 这里的ipos是原图像中的位置，要把它对应到padding后图像的ipos，才能进行插值计算。

double CubicPolynomial(unsigned char * image, int iImgWidth, int iImgHeight, double dx, double dy, int ipos)
{
	// image为图像各个位置的像素值， iImgWidth为图像的宽度, iImgHeight为图像的高度, ipos为当前插值的位置。

	// 先计算padding 图像的大小。
	int padding_size = (iImgHeight + 2) * (iImgWidth + 2);
	double * padding_img = new double[padding_size];

	// 接着对 padding_img 进行赋值。
	// 将原图像的值赋给padding_img的(2,2)-(iImgHeight+1, iImgWidth+1)
	// 该部分验证无误。
	for (int i = 0; i < iImgHeight; i++)
	{
		for (int j = 0; j < iImgWidth; j++)
		{
			int ipos1 = (i+1) * (iImgWidth+2) + (j+1);

			int temp = i*iImgWidth + j;

			padding_img[ipos1] = image[temp];
		}
	}


	// 第一行的值： 第二行、第三行、第四行的值的线性组合.
	// 该部分验证无误。

	for (int i = 1; i < iImgWidth+1; i++)
	{
		int ipos1 = i;

		double temp1 = padding_img[ipos1 + iImgWidth + 2];
		double temp2 = padding_img[ipos1 + 2 * (iImgWidth + 2)];
		double temp3 = padding_img[ipos1 + 3 * (iImgWidth + 2)];

		padding_img[ipos1] = 3 * temp1 - 3 * temp2 + temp3;
	}


	// 最后一行的值：倒数第二、第三、第四行的值的线性组合。
	// 该部分验证无误。

	for (int i = 1; i < iImgWidth+1; i++)
	{
		int ipos1 = i + (iImgHeight + 1)*(iImgWidth + 2);

		double temp1 = padding_img[ipos1 - (iImgWidth + 2)];
		double temp2 = padding_img[ipos1 - 2 * (iImgWidth + 2)];
		double temp3 = padding_img[ipos1 - 3 * (iImgWidth + 2)];

		padding_img[ipos1] = 3 * temp1 - 3 * temp2 + temp3;
	}

	// 第一列的值：第二列、第三列、第四列的值的线性组合。
	// 该部分验证无误。
	for (int i = 1; i < iImgHeight+1; i++)
	{
		int ipos1 = i * (iImgWidth + 2);
		
		double temp1 = padding_img[ipos1 + 1];
		double temp2 = padding_img[ipos1 + 2];
		double temp3 = padding_img[ipos1 + 3];

		padding_img[ipos1] = 3 * temp1 - 3 * temp2 + temp3;
	}

	// 最后一列的值：倒数第二、第三、第四列的值的线性组合。
	// 该部分验证无误。
	for (int i = 1; i < iImgHeight+1; i++)
	{
		int ipos1 = (i + 1)*(iImgWidth + 2) - 1;

		double temp1 = padding_img[ipos1 - 1];
		double temp2 = padding_img[ipos1 - 2];
		double temp3 = padding_img[ipos1 - 3];

		padding_img[ipos1] = 3 * temp1 - 3 * temp2 + temp3;
	}

	// 后面4个点的组合方向均为列的方向。即ipos+1, ipos+2, ipos+3，或ipos-1, ipos-2, ipos-3.

	// 左上角的点。
	padding_img[0] = 3 * padding_img[1] - 3 * padding_img[2] + padding_img[3];

	// 右上角的点。
	padding_img[iImgWidth + 1] = 3 * padding_img[iImgWidth] - 3 * padding_img[iImgWidth - 1] + padding_img[iImgWidth - 2];

	// 左下角的点。

	int pos_left_bottom = (iImgHeight + 1)*(iImgWidth + 2);
	padding_img[pos_left_bottom] = 3 * padding_img[pos_left_bottom + 1] - 3 * padding_img[pos_left_bottom + 2] + padding_img[pos_left_bottom + 3];

	// 右下角的点。

	int pos_right_bottom = (iImgHeight + 2)*(iImgWidth + 2) - 1;
	padding_img[pos_right_bottom] = 3 * padding_img[pos_right_bottom - 1] - 3 * padding_img[pos_right_bottom - 2] + padding_img[pos_right_bottom - 3];

	
	// padding结束之后，进行插值操作。首先要把原图的ipos对应到新图的ipos。

	int height = (ipos + 1) / iImgWidth;

	int width = (ipos + 1) - height*iImgWidth;

	// 在新图中的位置。

	int new_ipos = width + (height+1) * (iImgWidth + 2);


	// x方向卷积
	double dx_2 = dx*dx;
	double dx_3 = dx*dx*dx;
	double a4 = (-dx_3 + 2 * dx_2 - dx) / 2;
	double a5 = (3 * dx_3 - 5 * dx_2 + 2) / 2;
	double a6 = (-3 * dx_3 + 4 * dx_2 + dx) / 2;
	double a7 = (dx_3 - dx_2) / 2;



	double a8 = a4*padding_img[new_ipos - 1 - (iImgWidth + 2)] + a5*padding_img[new_ipos - (iImgWidth + 2)] + a6 * padding_img[new_ipos + 1 - (iImgWidth + 2)] + a7 * padding_img[new_ipos + 2 - (iImgWidth + 2)];
	double a9 = a4*padding_img[new_ipos - 1] + a5*padding_img[new_ipos] + a6 * padding_img[new_ipos + 1] + a7 * padding_img[new_ipos + 2];
	double a10 = a4*padding_img[new_ipos - 1 + (iImgWidth + 2)] + a5*padding_img[new_ipos + (iImgWidth + 2)] + a6 * padding_img[new_ipos + 1 + (iImgWidth + 2)] + a7 * padding_img[new_ipos + 2 + (iImgWidth + 2)];
	double a11 = a4*padding_img[new_ipos - 1 + 2 * (iImgWidth + 2)] + a5*padding_img[new_ipos + 2 * (iImgWidth + 2)] + a6 * padding_img[new_ipos + 1 + 2 * (iImgWidth + 2)] + a7 * padding_img[new_ipos + 2 + 2 * (iImgWidth + 2)];
	


	// y方向卷积
	double dy_2 = dy*dy;
	double dy_3 = dy*dy*dy;
	double a12 = (-dy_3 + 2 * dy_2 - dy) / 2;
	double a13 = (3 * dy_3 - 5 * dy_2 + 2) / 2;
	double a14 = (-3 * dy_3 + 4 * dy_2 + dy) / 2;
	double a15 = (dy_3 - dy_2) / 2;

	double gray_value;
	gray_value = a8*a12 + a9*a13 + a10*a14 + a11*a15;

	return gray_value;


}


int main(int argc, char **argv)
{

	// imageData为图像数据，unsigned char * 型

	Mat image = imread("F:\\img\\image\\result1.jpg", 0);
	int image_size = image.cols * image.rows;
	unsigned char* imageData = new unsigned char[image_size];

	int a = 0;
	for (int i = 0; i<image.rows; i++)
	{
		for (int j = 0; j<image.cols; j++)
		{
			imageData[a] = image.at<uchar>(i, j);
			a++;
		}
	}



	int iImgHeight = image.rows;
	int iImgWidth = image.cols;

	// 进行多次随机试验。

	for (int i = 0; i < 1000; i++)
	{

		// x是列方向，y是行方向，在cv::Mat做输入的函数中，起始点为(1,1)点。



		// 生成随机位置。
		int x = rand() % (iImgWidth - 1) + 1;
		int y = rand() % (iImgHeight - 1) + 1;

		// 对应回原图像下的ipos。
		int ipos = (y - 1)*(iImgWidth)+(x - 1);


		// 生成随机的dx, dy.
		int N = 10000;
		double dx = rand() % (N + 1) / (float)(N + 1);;
		double dy = rand() % (N + 1) / (float)(N + 1);;


		// 基于unsigned char * 数据的插值结果。
		double temp1 = CubicPolynomial(imageData, iImgWidth, iImgHeight, dx, dy, ipos);

		// 基于cv::Mat 数据的插值结果。
		double temp2 = CubicPolynomial(image, dx, dy, x, y);

		double difference = abs(temp2 - temp1);

		cout << "The difference is : " << difference << endl;


	}


	system("pause");
	return 0;
}

四、Matlab代码

function Value = CubicInterpoly(y,x,img)
    uw = floor(x);
    vw = floor(y);
    pV = zeros(16,1);
    W = zeros(16,2);
    
    pV(1) = img(vw-1,uw-1);
    W(1,:) = BiCubic(y,x,vw-1,uw-1);
    pV(2) = img(vw-1,uw);
    W(2,:) = BiCubic(y,x,vw-1,uw);
    pV(3) = img(vw-1,uw+1);
    W(3,:) = BiCubic(y,x,vw-1,uw+1);
    pV(4) = img(vw-1,uw+2);
    W(4,:) = BiCubic(y,x,vw-1,uw+2);

    pV(5) = img(vw,uw-1);
    W(5,:) = BiCubic(y,x,vw,uw-1);
    pV(6) = img(vw,uw);
    W(6,:) = BiCubic(y,x,vw,uw);
    pV(7) = img(vw,uw+1);
    W(7,:) = BiCubic(y,x,vw,uw+1);
    pV(8) = img(vw,uw+2);
    W(8,:) = BiCubic(y,x,vw,uw+2);
    
    pV(9) = img(vw+1,uw-1);
    W(9,:) = BiCubic(y,x,vw+1,uw-1);
    pV(10) = img(vw+1,uw);
    W(10,:) = BiCubic(y,x,vw+1,uw);
    pV(11) = img(vw+1,uw+1);
    W(11,:) = BiCubic(y,x,vw+1,uw+1);
    pV(12) = img(vw+1,uw+2);
    W(12,:) = BiCubic(y,x,vw+1,uw+2);
    
    pV(13) = img(vw+2,uw-1);
    W(13,:) = BiCubic(y,x,vw+2,uw-1);
    pV(14) = img(vw+2,uw);
    W(14,:) = BiCubic(y,x,vw+2,uw);
    pV(15) = img(vw+2,uw+1);
    W(15,:) = BiCubic(y,x,vw+2,uw+1);
    pV(16) = img(vw+2,uw+2);
    W(16,:) = BiCubic(y,x,vw+2,uw+2);
    
    Value = 0;
    for i=1:16
        Value = Value + pV(i)*W(i,1)*W(i,2);
    end
end

function W = BiCubic(y,x,yi,xi)
    W = zeros(1,2);
    a = -0.5;
    distX = abs(x-xi);
    if(distX <=1)
        W(1) = (a+2)*distX^3-(a+3)*distX^2+1;
    elseif(distX > 1 && distX <2)
        W(1) = a*distX^3-5*a*distX^2+8*a*distX-4*a;
    else
        W(1) = 0;
    end
    distY = abs(y-yi);
    if(distY <=1)
        W(2) = (a+2)*distY^3-(a+3)*distY^2+1;
    elseif(distY > 1 && distY <2)
        W(2) = a*distY^3-5*a*distY^2+8*a*distY-4*a;
    else
        W(2) = 0;
    end
end