首先讲插值问题：
假设存在这样一个3x3的矩阵：
1 3 5
7 9 11
13 15 17
如果我们想把它放大为4x4的矩阵：
？ ？ ？ ？
？ ？ ？ ？
？ ？ ？ ？
？ ？ ？ ？
那么这些问号该怎么填值呢？
简单讲一些基础的方法。
邻近插值：
在上面的两个矩阵中，假设定义左上角为0点，每个点可以用二维(x,y)表示，按惯例，x表示横轴，y表示纵轴。
并设源3x3矩阵的每个点用(x,y)表示、宽高用w和h表示，目标4x4矩阵的每个点用(x',y')表示、宽高用w'和h'表示，可以推出下面的等式：
x' = x * (w' / w)
y' = y * (h' / h)
x = x' * (w / w')
y = y' * (h / h')
现在需要解决的问题是，已知x'和y'，需要求出对应的x和y，以便填写这些问号。
我们先套用上面的这两个公式：
x = x' * (w / w')
y = y' * (h / h')
求取4x4矩阵，首先计算(0,0)，得：
x = 0 * (3/4) = 0
y = 0 * (3/4) = 0
即对应的源坐标为(0,0)。
同样的方法，计算(1,0)，得：
x = 1 * (3/4) = 0.75
y = 0 * (3/4) = 0
这里对应的源坐标为(0.75,0)，坐标怎么能为小数呢，只能取整，如果按四舍五入取整，得到源坐标为(1,0)。
如法计算，得到整幅4x4的源坐标对应关系，如下：
(0,0) (1,0) (2,0) (2,0)
(0,1) (1,1) (2,1) (2,1)
(0,2) (1,2) (2,2) (2,2)
(0,2) (1,2) (2,2) (2,2)
将值替换，即：
1 3 5 5
7 9 11 11
13 15 17 17
13 15 17 17
这种放大图像的方法叫做最临近插值算法，这是一种最基本、最简单的图像缩放算法，效果也是最不好的，放大后的图像有很严重的马赛克，缩小后的图像有很严重的失真，效果不好的根源就是其简单的最临近插值方法引入了严重的图像失真。比如，当由目标图的坐标反推得到的源图的的坐标是一个浮点数的时候，采用了四舍五入的方法，直接采用了和这个浮点数最接近的象素的值，这种方法是很不科学的，当推得坐标值为 0.75的时候，不应该就简单的取为1，既然是0.75，比1要小0.25 ，比0要大0.75，那么目标象素值其实应该根据这个源图中虚拟的点四周的四个真实的点来按照一定的规律计算出来的，这样才能达到更好的缩放效果。
附邻近插值实现：
namespace dakuang
{
// 邻近插值
void insertByNear(const cv::Mat& matSrc, cv::Mat& matDst)
{
for (int y = 0; y < matDst.rows; ++y)
{
for (int x = 0; x < matDst.cols; ++x)
{
int nSrcY = std::round( y * ((double)matSrc.rows / matDst.rows) );
int nSrcX = std::round( x * ((double)matSrc.cols / matDst.cols) );
nSrcY = std::min(nSrcY, matSrc.rows - 1);
nSrcX = std::min(nSrcX, matSrc.cols - 1);
matDst.at<uchar>(y, x) = matSrc.at<uchar>(nSrcY, nSrcX);
}
}
}
}
也可以使用opencv的resize()方法，最后一个参数取INTER_NEAREST。
双线性插值：
鉴于临近插值的的缺点，尤其是在图像缩放场合，需要一种更贴近真实感的算法。
在讲双线性插值方法前，先看看反向距离权重IDW的定义：
反向距离权重 (IDW) 插值显式假设：彼此距离较近的事物要比彼此距离较远的事物更相似。当为任何未测量的位置预测值时，反距离权重法会采用预测位置周围的测量值。与距离预测位置较远的测量值相比，距离预测位置最近的测量值对预测值的影响更大。反距离权重法假定每个测量点都有一种局部影响，而这种影响会随着距离的增大而减小。由于这种方法为距离预测位置最近的点分配的权重较大，而权重却作为距离的函数而减小，因此称之为反距离权重法。
IDW具体的计算方法为，假若存在两个点，其值分别为a和b，ab之间的距离为l，则ab之间离a相距m的值为：
a((l-m)/l) + b(m/l)
令l为1，简化公式为：
a(1-m) + bm
回到双线性插值的定义：
利用源图中虚拟点四周的四个真实存在的像素值，使用IDW来共同计算目标图中的一个像素值。由于图像属于2维数据，对虚拟点取值的计算，需要进行两轮运算。借用网络上的图（坐标系与图片的不一致的）：
第1轮：分别计算P点的两条水平IDW插值。
即计算Q11和Q21的IDW插值R1，计算Q12和Q22的IDW插值R2。
第2轮：计算两个水平IDW插值间的垂直IDW插值。
即计算R1和R2的IDW插值P。
假设Q11、Q21、Q12、Q22、P的坐标分别为：
(x1,y1) (x2,y1) (x1,y2) (x2,y2) (x,y)
则R1的值为：
f(x,y1) = ((x2-x)/(x2-x1)) * f(Q11) + ((x-x1)/(x2-x1)) * f(Q21)
R2的值为：
f(x,y2) = ((x2-x)/(x2-x1)) * f(Q12) + ((x-x1)/(x2-x1)) * f(Q22)
P的值为：
f(x,y) = ((y2-y)/(y2-y1)) * f(x,y1) + ((y-y1)/(y2-y1)) * f(x,y2)
在处理图像缩放的时候，如果求取的点在源矩阵中不存在，使用双线性插值法，首先应该找到该虚拟点周围的四个像素，然后使用上面的公式计算插值。对前面的3x3矩阵采用双线性插值法，计算出的4x4矩阵值为：
1 2 4 5
5 6 8 9
10 11 13 14
13 14 16 17
附双线性插值实现：
namespace dakuang
{
// 双线性插值
void insertByLiner(const cv::Mat& matSrc, cv::Mat& matDst)
{
for (int y = 0; y < matDst.rows; ++y)
{
for (int x = 0; x < matDst.cols; ++x)
{
double fSrcY = y * ((double)matSrc.rows / matDst.rows);
double fSrcX = x * ((double)matSrc.cols / matDst.cols);
int nIntY = (int)fSrcY;
int nIntX = (int)fSrcX;
double fSmallY = fSrcY - nIntY;
double fSmallX = fSrcX - nIntX;
int nYPlus1 = std::min(nIntY + 1, matSrc.rows - 1);
int nXPlus1 = std::min(nIntX + 1, matSrc.cols - 1);
// 提取4个点的值
int nYX = matSrc.at<uchar>(nIntY, nIntX);
int nYX1 = matSrc.at<uchar>(nIntY, nXPlus1);
int nY1X = matSrc.at<uchar>(nYPlus1, nIntX);
int nY1X1 = matSrc.at<uchar>(nYPlus1, nXPlus1);
// x轴的两条线性融合
// (x,y) 与 (x+1,y) 的线性融合
int nBlendX1 = (1.0 - fSmallX) * nYX + fSmallX * nYX1;
// (x,y+1) 与 (x+1,y+1) 的线性融合
int nBlendX2 = (1.0 - fSmallX) * nY1X + fSmallX * nY1X1;
// y轴的一条线性融合
int nBlendY = (1.0 - fSmallY) * nBlendX1 + fSmallY * nBlendX2;
matDst.at<uchar>(y, x) = nBlendY;
}
}
}
}
也可以使用opencv的resize()方法，最后一个参数取INTER_LINEAR。
散列数据插值：
前面提到的源矩阵，当需要插值时，由于其数据是规则的，很容易找到四周的邻近数据，无论是进行就近取值，还是IDW插值。
有一些情况下，比如地质勘探场合，或者是一些立体图像扫描仪场景，需要根据事先测量的少量的位置不规则的点，求取整个曲面。
源图：
目标图：
有两种方法可用来求取这个问题：
IDW加权求和：
设空间待插点为P(Xp, Yp, Zp)，P点邻域内已知散乱点Qi(Xi, Yi, Zi)，i=1,2,...,n，利用反向距离加权对P点的属性Zp进行插值。
其插值原理是待插点的属性值是待插点邻域内已知散乱点属性值的加权平均，权的大小与待插点邻域内混乱点之间的距离有关，
是距离k(0<=k<=2)(k一般取2)次方的倒数。即
其中di为待插点与其邻域内第i个点之间的距离。
克里金法：
该方法基于包含自相关（即，测量点之间的统计关系）的统计模型。因此，地统计方法不仅具有产生预测表面的功能，而且能够对预测的确定性或准确性提供某种度量。克里金法假定采样点之间的距离或方向可以反映可用于说明表面变化的空间相关性。克里金法工具可将数学函数与指定数量的点或指定半径内的所有点进行拟合以确定每个位置的输出值。克里金法是一个多步过程；它包括数据的探索性统计分析、变异函数建模和创建表面，还包括研究方差表面。当您了解数据中存在空间相关距离或方向偏差后，便会认为克里金法是最适合的方法。该方法通常用在土壤科学和地质中。
由于克里金法支持不同的样本模型预测，经多年的理论和实践证明，克里金法的实际效果要优于IDW方法。鉴于克里金法的复杂性，本文暂时不予发挥。
下面我将以IDW加权求和法，根据前面的公式，编写插值函数，代码如下：
namespace dakuang
{
// 使用IDW搜索法进行插值
// 要求源和目标矩阵格式都为CV_32FC1
// 要求源矩阵待填值的位置初使值为0
void fillByIDWSearch(const cv::Mat& matSrc, cv::Mat& matDst)
{
if (matSrc.type() != CV_32FC1)
return;
matDst = matSrc.clone();
for (int y = 0; y < matSrc.rows; ++y)
{
for (int x = 0; x < matSrc.cols; ++x)
{
float fValSrc = matSrc.at<float>(y, x);
if (fValSrc > 0)
continue;
// 分子、分母积分
float fSumTop = 0.0, fSumBottom = 0.0;
// 积分的样本数
int nSampleCount = 0;
// 逐渐扩大邻域搜索距离
for (int d = 1; d < std::max(matSrc.rows, matSrc.cols); ++d)
{
// 搜索距离为d时的矩形区域
int sx1 = x - d;
int sx2 = x + d;
int sy1 = y - d;
int sy2 = y + d;
// 上邻域
if (sy1 >= 0)
{
for (int sx = sx1; sx <= sx2; ++sx)
{
if (sx >= 0 && sx < matSrc.cols)
{
// 寻找有效点
float fVal = matSrc.at<float>(sy1, sx);
if (fVal > 0)
{
fSumTop += fVal / (d * d);
fSumBottom += 1.0 / (d * d);
nSampleCount++;
}
}
}
}
// 下邻域
if (sy2 < matSrc.rows)
{
for (int sx = sx1; sx <= sx2; ++sx)
{
if (sx >= 0 && sx < matSrc.cols)
{
// 寻找有效点
float fVal = matSrc.at<float>(sy2, sx);
if (fVal > 0)
{
fSumTop += fVal / (d * d);
fSumBottom += 1.0 / (d * d);
nSampleCount++;
}
}
}
}
// 左邻域
if (sx1 >= 0)
{
for (int sy = sy1 + 1; sy <= sy2 - 1; ++sy)
{
if (sy >= 0 && sy < matSrc.rows)
{
// 寻找有效点
float fVal = matSrc.at<float>(sy, sx1);
if (fVal > 0)
{
fSumTop += fVal / (d * d);
fSumBottom += 1.0 / (d * d);
nSampleCount++;
}
}
}
}
// 右邻域
if (sx2 < matSrc.cols)
{
for (int sy = sy1 + 1; sy <= sy2 - 1; ++sy)
{
if (sy >= 0 && sy < matSrc.rows)
{
// 寻找有效点
float fVal = matSrc.at<float>(sy, sx2);
if (fVal > 0)
{
fSumTop += fVal / (d * d);
fSumBottom += 1.0 / (d * d);
nSampleCount++;
}
}
}
}
// 最多参照不远大于8的样本
if (nSampleCount >= 8)
break;
}
if (fSumBottom != 0.0)
{
matDst.at<float>(y, x) = fSumTop / fSumBottom;
}
}
}
}
}
编写测试代码：
cv::Mat mat1(8, 8, CV_8UC1, cv::Scalar(0));
mat1.at<uchar>(0, 3) = 10;
mat1.at<uchar>(0, 6) = 30;
mat1.at<uchar>(2, 2) = 50;
mat1.at<uchar>(2, 4) = 70;
mat1.at<uchar>(3, 6) = 90;
mat1.at<uchar>(4, 3) = 110;
mat1.at<uchar>(5, 1) = 130;
mat1.at<uchar>(5, 5) = 150;
mat1.at<uchar>(7, 2) = 170;
mat1.convertTo(mat1, CV_32FC1, 1.0 / 255);
cv::Mat mat2;
dakuang::fillByIDWSearch(mat1, mat2);
mat2.convertTo(mat2, CV_8UC1, 255);
这里构造了一个8x8的矩阵，并且模拟填充一些采样点，然后进行插值填充，效果如下：
源图：
填充效果图：
