很长时间不敲代码，感觉一闲下来就忘了很多。想着把一些图像算法自己实现一遍也好，一方面加深算法的学习和理解，另一方面又可以练练编码能力。对于我这个非科班出身的，也挺有好处的。

不管那么多，先把这个坑挖了。闲着想敲代码了就慢慢填。

下面先把图像处理中常见的三种插值算法实现了。先解释一下什么叫插值。老样子先看看维基百科怎么解释。

数学的数值分析领域中，内插或称插值（英语：interpolation)是一种通过已知的、离散的数据点，在范围内推求新数据点的过程或方法。

450px-Splined_epitrochoid.svg.png

一组离散数据点在一个外延的插值。曲线中实际已知数据点是红色的；连接它们的蓝色曲线即为插值。

再举个例子：

x1 = 1, y1 = 3
x2 = 3, y2 = 7
x3 = 5, y3 = 24

求x = 4时， y = ？

所以可以看到，可以这么理解插值，在一个函数里面，自变量是离散有间隔的，插值就是往自变量的间隔之间插入新的自变量，然后求解新的自变量函数值。这有什么作用呢，从上图可以看到，散点图是可以利用插值来拟合曲线的，蓝色的就是插入的密密麻麻的点。

然而在图像处理上的应用可以体现在图像的缩放上面。如放大一张图片，在像素点的层面上其实就是往像素点之间插入新的像素点从而增大图像的分辨率。

插值算法有很多种，具体可以参考维基百科插值

这里就实现图像处理比较常用的三种：最近邻域差值、双线性插值、双三次插值

最近邻域插值

最近邻域内插法(Nearest Neighbor interpolation)

最近邻域是三种插值之中最简单的一种，原理就是选取距离插入的像素点（x+u, y+v)【注：x,y为整数， u，v为小数】最近的一个像素点，用它的像素点的灰度值代替插入的像素点。说到距离，顺便复习一下像素点之间的三种空间距离。

对于像素p、q和z，分别具有坐标(x,y)，(s,t)和(u,v)，如果

(1) D(p,q) ≥ 0 (当且仅当p=q时，D(p,q)=0)

(2) D(p,q) = D(q,p)

(3) D(p,z) ≤ D(p,q) + D(q,z)

则称D是距离函数或度量

欧氏距离.gif

D4距离（城市距离）

D4距离.gif

D8距离（棋盘距离）

D8距离.gif

这里代码实现采用的是欧式距离

//最近邻域插值
//根据目标图像的像素点（浮点坐标）找到原始图像中的4个像素点，取距离该像素点最近的一个原始像素值作为该点的值。

#include<opencv2/opencv.hpp>
#include<cassert>

namespace mycv {
    void nearestIntertoplation(cv::Mat& src, cv::Mat& dst, const int rows, const int cols);
}//mycv

double distance(const double x1, const double y1, const double x2, const double y2);//两点之间距离，这里用欧式距离

int main()
{
    cv::Mat img = cv::imread("lena.jpg", 0);
    if (img.empty()) return -1;
    cv::Mat dst;
    mycv::nearestIntertoplation(img, dst, 600, 600);
    cv::imshow("img", img);
    cv::imshow("dst", dst);
    cv::waitKey(0);
    return 0;
    return 0;
}//main

void mycv::nearestIntertoplation(cv::Mat& src, cv::Mat& dst, const int rows, const int cols)
{
    //比例尺
    const double scale_row = static_cast<double>(src.rows) / rows;
    const double scale_col = static_cast<double>(src.rows) / cols;

    //扩展src到dst
    dst = cv::Mat(rows, cols, src.type());
    assert(src.channels() == 1 && dst.channels() == 1);

    for (int i = 0; i < rows; ++i)//dst的行
        for (int j = 0; j < cols; ++j)//dst的列
        {
            //求插值的四个点
            double y = (i + 0.5) * scale_row + 0.5;
            double x = (j + 0.5) * scale_col + 0.5;
            int x1 = static_cast<int>(x);//col对应x
            if (x1 >= (src.cols - 2)) x1 = src.cols - 2;//防止越界
            int x2 = x1 + 1;
            int y1 = static_cast<int>(y);//row对应y
            if (y1 >= (src.rows - 2))  y1 = src.rows - 2;
            int y2 = y1 + 1;
            //根据目标图像的像素点（浮点坐标）找到原始图像中的4个像素点，取距离该像素点最近的一个原始像素值作为该点的值。
            assert(0 < x2 && x2 < src.cols && 0 < y2 &&  y2 < src.rows);
            std::vector<double> dist(4);
            dist[0] = distance(x, y, x1, y1);
            dist[1] = distance(x, y, x2, y1);
            dist[2] = distance(x, y, x1, y2);
            dist[3] = distance(x, y, x2, y2);

            int min_val = dist[0];
            int min_index = 0;
            for (int i = 1; i < dist.size(); ++i)
                if (min_val > dist[i])
                {
                    min_val = dist[i];
                    min_index = i;
                }

            switch (min_index)
            {
            case 0:
                dst.at<uchar>(i, j) = src.at<uchar>(y1, x1);
                break;
            case 1:
                dst.at<uchar>(i, j) = src.at<uchar>(y1, x2);
                break;
            case 2:
                dst.at<uchar>(i, j) = src.at<uchar>(y2, x1);
                break;
            case 3:
                dst.at<uchar>(i, j) = src.at<uchar>(y2, x2);
                break;
            default:
                assert(false);
            }           
        }
}

double distance(const double x1, const double y1, const double x2, const double y2)//两点之间距离，这里用欧式距离
{
    return (x1 - x2)*(x1 - x2) + (y1 - y2)*(y1 - y2);//只需比较大小，返回距离平方即可
}

双线性插值

双线性插值，顾名思义，在像素点矩阵上面，x和y两个方向的线性插值所得的结果。那么先看看线性插值

Linear_interpolation.png

image.png

image.png

以上是一维的，接下来看看二维中的双线性插值

Bilinear_interpolation.png

首先在x方向上面线性插值，得到R2、R1

image.png

然后以R2,R1在y方向上面再次线性插值

image.png

如果选择一个坐标系统使得 f 的四个已知点坐标分别为 (0, 0)、(0, 1)、(1, 0) 和 (1, 1)，那么插值公式就可以化简为

image.png

用矩阵表示

image.png

下面贴上代码实现

//线性内插（双线性插值)
#include<opencv2/opencv.hpp>
#include<opencv2/core/matx.hpp>
#include<cassert>

namespace mycv {
    void bilinearIntertpolatioin(cv::Mat& src, cv::Mat& dst, const int rows, const int cols);
}//mycv

int main()
{
    cv::Mat img = cv::imread("lena.jpg", 0);
    if (img.empty()) return -1;
    cv::Mat dst;
    mycv::bilinearIntertpolatioin(img, dst, 600, 600);
    cv::imshow("img", img);
    cv::imshow("dst", dst);
    cv::waitKey(0);
    return 0;
}//main

void mycv::bilinearIntertpolatioin(cv::Mat& src, cv::Mat& dst, const int rows, const int cols)
{
    //比例尺
    const double scale_row = static_cast<double>(src.rows) / rows;
    const double scale_col = static_cast<double>(src.rows) / cols;

    //扩展src到dst
    dst = cv::Mat(rows, cols, src.type());
    assert(src.channels() == 1 && dst.channels() == 1);
    
    for(int i = 0; i < rows; ++i)//dst的行
        for (int j = 0; j < cols; ++j)//dst的列
        {
            //求插值的四个点
            double y = (i + 0.5) * scale_row + 0.5;
            double x = (j + 0.5) * scale_col + 0.5;
            int x1 = static_cast<int>(x);//col对应x
            if (x1 >= (src.cols - 2)) x1 = src.cols - 2;//防止越界
            int x2 = x1 + 1;
            int y1 = static_cast<int>(y);//row对应y
            if (y1 >= (src.rows - 2))  y1 = src.rows - 2;
            int y2 = y1 + 1;

            assert(0 < x2 && x2 < src.cols && 0 < y2 &&  y2 < src.rows);
            //插值公式,参考维基百科矩阵相乘的公式https://zh.wikipedia.org/wiki/%E5%8F%8C%E7%BA%BF%E6%80%A7%E6%8F%92%E5%80%BC
    
            cv::Matx12d matx = { x2 - x, x - x1 };
            cv::Matx22d matf = { static_cast<double>(src.at<uchar>(y1, x1)), static_cast<double>(src.at<uchar>(y2, x1)),
                                 static_cast<double>(src.at<uchar>(y1, x2)), static_cast<double>(src.at<uchar>(y2, x2)) };
            cv::Matx21d maty = {
                y2 - y,
                y - y1
            };

            auto  val = (matx * matf * maty);
            dst.at<uchar>(i, j) = val(0,0);
        }

}

双三次插值

双三次插值就稍有点复杂了，先看维基百科

在数值分析这个数学分支中，双三次插值（英语：Bicubic interpolation）是二维空间中最常用的插值方法。在这种方法中，函数 f 在点 (x, y) 的值可以通过矩形网格中最近的十六个采样点的加权平均得到，在这里需要使用两个多项式插值三次函数，每个方向使用一个

双三次插值计算公式

image.png

即

image.png

那么这个a(i, j)便是介绍里面所说的加权系数了，所以关键是要把它求解出来。

求解加权系数的公式如下

image.png

求解加强系数代码如下

std::vector<double> mycv::getW(double coor, double a)
{
    std::vector<double> w(4);
    int base = static_cast<int>(coor);//取整作为基准
    double e = coor - static_cast<double>(base);//多出基准的小数部分
    std::vector<double> tmp(4);//存放公式中 |x| <= 1,  1 < |x| < 2四个值
    // 4 x 4的16个点，所以tmp[0]和tmp[4]距离较远，值在[1, 2]区间
    tmp[0] = 1.0 + e;// 1 < x < 2
    tmp[1] = e;//x <= 1
    tmp[2] = 1.0 - e;// x <= 1
    tmp[3] = 2.0 - e;// 1 < x < 2

    //按照bicubic公式计算系数w
    // x <= 1
    w[1] = (a + 2.0) * std::abs(std::pow(tmp[1], 3)) - (a + 3.0) * std::abs(std::pow(tmp[1], 2)) + 1;
    w[2] = (a + 2.0) * std::abs(std::pow(tmp[2], 3)) - (a + 3.0) * std::abs(std::pow(tmp[2], 2)) + 1;
    // 1 < x < 2
    w[0] = a * std::abs(std::pow(tmp[0], 3)) - 5.0 * a * std::abs(std::pow(tmp[0], 2)) + 8.0*a*std::abs(tmp[0]) - 4.0*a;
    w[3] = a * std::abs(std::pow(tmp[3], 3)) - 5.0 * a * std::abs(std::pow(tmp[3], 2)) + 8.0*a*std::abs(tmp[3]) - 4.0*a;

    return w;
}

求解出来之后，wx和wy就是4×4组成的16个系数，与插值所需要用到4×4的16个点相匹配。插值步骤的关键代码

                //4x4数量的点(rr, cc) -> (y, x)
                std::vector<std::vector<int> > src_arr = {
                    { src.at<uchar>(rr - 1, cc - 1), src.at<uchar>(rr, cc - 1), src.at<uchar>(rr + 1, cc - 1), src.at<uchar>(rr + 2, cc - 1)},
                    { src.at<uchar>(rr - 1, cc), src.at<uchar>(rr, cc), src.at<uchar>(rr + 1, cc), src.at<uchar>(rr + 2, cc)},
                    { src.at<uchar>(rr - 1, cc + 1), src.at<uchar>(rr, cc + 1), src.at<uchar>(rr + 1, cc + 1), src.at<uchar>(rr + 2, cc + 1)},
                    { src.at<uchar>(rr - 1, cc + 2), src.at<uchar>(rr, cc + 2), src.at<uchar>(rr + 1, cc + 2), src.at<uchar>(rr + 2, cc + 2)}
                };
                for(int p = 0; p < 3; ++p)
                    for (int q = 0; q < 3; ++q)
                    {
                        //val(p, q) = w(p,q) * src(p, q)
                        val += wr[p] * wc[q] * static_cast<double>(src_arr[p][q]);
                    }
                assert(i < dst.rows && j < dst.cols);
                dst.at<uchar>(i, j) = static_cast<int>(val);

下面放出完整代码

//双三次插值

#include<opencv2/opencv.hpp>
#include<iostream>
#include<vector>
#include<cassert>

namespace mycv {
    void bicubicInsterpolation(cv::Mat& src, cv::Mat& dst, const int rows, const int cols);
    std::vector<double> getW(double coor, double a = -0.5);//a默认-0.5
}//mycv

int main(void)
{
    cv::Mat img = cv::imread("lena.jpg", 0);
    if (img.empty()) return -1;
    cv::Mat dst;
    mycv::bicubicInsterpolation(img, dst, 600, 600);
    cv::imshow("img", img);
    cv::imshow("dst", dst);
    cv::waitKey(0);
    return 0;
}//main

void mycv::bicubicInsterpolation(cv::Mat& src, cv::Mat& dst, const int rows, const int cols)
{
    dst = cv::Mat(rows, cols, src.type(), cv::Scalar::all(0));//初始化dst
    //比例尺
    double row_scale = static_cast<double>(src.rows) / rows;
    double col_scale = static_cast<double>(src.cols) / cols;

    switch (src.channels())
    {

    case 1://灰度
        for(int i = 2; i < dst.rows - 2; ++i)
            for (int j = 2; j < dst.cols - 2; ++j)
            {
                //计算系数w
                double r = static_cast<double>(i *  row_scale);
                double c = static_cast<double>(j *  col_scale);
                
                //防止越界
                if (r < 1.0) r += 1.0;
                if (c < 1.0) c += 1.0;
            

                std::vector<double> wr = mycv::getW( r);
                std::vector<double> wc = mycv::getW( c);
                
                //最后计算插值得到的灰度值
                double val = 0;
                int cc = static_cast<int>(c);
                int rr = static_cast<int>(r);
                //防止越界
                if (cc > src.cols - 3)
                {
                    cc = src.cols - 3;
                    
                }
                if (rr > src.rows - 3) rr = src.rows - 3;

                assert(0 <= (rr - 1) && 0 <= (cc - 1) && (rr + 2) < src.rows && (cc + 2) < src.cols);

                //4x4数量的点(rr, cc) -> (y, x)
                std::vector<std::vector<int> > src_arr = {
                    { src.at<uchar>(rr - 1, cc - 1), src.at<uchar>(rr, cc - 1), src.at<uchar>(rr + 1, cc - 1), src.at<uchar>(rr + 2, cc - 1)},
                    { src.at<uchar>(rr - 1, cc), src.at<uchar>(rr, cc), src.at<uchar>(rr + 1, cc), src.at<uchar>(rr + 2, cc)},
                    { src.at<uchar>(rr - 1, cc + 1), src.at<uchar>(rr, cc + 1), src.at<uchar>(rr + 1, cc + 1), src.at<uchar>(rr + 2, cc + 1)},
                    { src.at<uchar>(rr - 1, cc + 2), src.at<uchar>(rr, cc + 2), src.at<uchar>(rr + 1, cc + 2), src.at<uchar>(rr + 2, cc + 2)}
                };
                for(int p = 0; p < 3; ++p)
                    for (int q = 0; q < 3; ++q)
                    {
                        //val(p, q) = w(p,q) * src(p, q)
                        val += wr[p] * wc[q] * static_cast<double>(src_arr[p][q]);
                    }
                assert(i < dst.rows && j < dst.cols);
                dst.at<uchar>(i, j) = static_cast<int>(val);
                
            }
        break;
    case 3://彩色(原理一样多了两个通道而已)
        break;
    default:
        break;
    }
}

std::vector<double> mycv::getW(double coor, double a)
{
    std::vector<double> w(4);
    int base = static_cast<int>(coor);//取整作为基准
    double e = coor - static_cast<double>(base);//多出基准的小数部分
    std::vector<double> tmp(4);//存放公式中 |x| <= 1,  1 < |x| < 2四个值
    // 4 x 4的16个点，所以tmp[0]和tmp[4]距离较远，值在[1, 2]区间
    tmp[0] = 1.0 + e;// 1 < x < 2
    tmp[1] = e;//x <= 1
    tmp[2] = 1.0 - e;// x <= 1
    tmp[3] = 2.0 - e;// 1 < x < 2

    //按照bicubic公式计算系数w
    // x <= 1
    w[1] = (a + 2.0) * std::abs(std::pow(tmp[1], 3)) - (a + 3.0) * std::abs(std::pow(tmp[1], 2)) + 1;
    w[2] = (a + 2.0) * std::abs(std::pow(tmp[2], 3)) - (a + 3.0) * std::abs(std::pow(tmp[2], 2)) + 1;
    // 1 < x < 2
    w[0] = a * std::abs(std::pow(tmp[0], 3)) - 5.0 * a * std::abs(std::pow(tmp[0], 2)) + 8.0*a*std::abs(tmp[0]) - 4.0*a;
    w[3] = a * std::abs(std::pow(tmp[3], 3)) - 5.0 * a * std::abs(std::pow(tmp[3], 2)) + 8.0*a*std::abs(tmp[3]) - 4.0*a;

    return w;
}

三次插值的代码计算量特别大，而且代码没什么优化，速度是慢的可以。但感觉可读性还是ok得，便于理解吧。