双三次插值 - 插值图像任意位置亚像素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:16Value = Value + pV(i)*W(i,1)*W(i,2);end
endfunction 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;elseW(1) = 0;enddistY = 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;elseW(2) = 0;end
end