图像缩放算法 最近邻插值算法 & 双线性插值算法

图像缩放：

resize函数是在OpenCV中经常使用的函数，功能是将一副加载到Mat中的图像调整尺寸或者按照比例进行缩放。

其中有两种简单又常用的插值算法用来实现图像缩放，分别是最近邻插值算法和双线性插值算法。

最近邻插值算法:

最近邻插值算法的思想十分简单

设原始图像src的高为h，宽为w，图像上像素值为(x,y)。

设目标图像dst的高H，宽为W，图像上的像素值为(X,Y)。

那么根据比例缩放的思想有$\frac{x}{X}=\frac{w}{W}$
同理，纵向上的像素对应比值为$\frac{y}{Y}=\frac{h}{H}$

那么在考虑目标图像dst上面的像素点(X,Y)对应到原始图像src上面像素点的位置为

$(x,y)=(\frac{w}{W}X,\frac{h}{H}Y)$

这里要对$\frac{w}{W}X$与$\frac{h}{H}Y$两个数值进行取整，在OpenCV源码当中，是用cvFloor函数取整，即舍去小数部分。

代码：

void InterNearestResize(const Mat &src, Mat &dst, Size &dsize, double fx = 0.0, double fy = 0.0)//原始图形以及缩放比例
{

    Size ssize = src.size();//获取矩阵大小
    CV_Assert(ssize.area() > 0);//保证矩阵的长宽都大于0

    if (!dsize.area())//如果dsize为(0,0)
    {
        dsize = Size(saturate_cast<int>(src.cols*fx),//satureate_cast防止数据溢出
            saturate_cast<int>(src.rows*fy));

        CV_Assert(dsize.area());
    }
    else
    {
        fx = (double)dsize.width / src.cols;//Size中的宽高和mat中的行列是相反的
        fy = (double)dsize.height / src.rows;
    }

    dst.create(Size(src.cols*fy,src.rows*fx), src.type());//将dst矩阵大小变为想要的尺寸

    fx = 1.0 / fx;
    fy = 1.0 / fy;
    uchar *ps = src.data;
    uchar *pd = dst.data;
    int channels = src.channels(),x,y,out,in;

    for (int row = 0; row < dst.rows; row++)
    {
        x = cvFloor(row * fx);
        for (int col = 0; col < dst.cols; col++)
        {
            y = cvFloor(col * fy);
            for (int c = 0; c < channels; c++)
            {
                out = (row * dst.cols + col) * channels + c;
                in = (x * src.cols + y) * channels + c;
                pd[out] = ps[in];
            }

            //dst.at<Vec3b>(row, col) = src.at<Vec3b>(x, y);
        }
    }
}

双线性插值算法:

如果原始图像src的大小是3×3，目标图像dst的大小是4×4，考虑dst中(1,1)点像素对应原始图像像素点的位置为(0.75,0.75)，如果使用最近邻算法来计算，原始图像的位置在浮点数取整后为坐标(0,0)。

上面这样粗暴的计算会丢失很多信息，考虑(0.75,0.75)这个信息，它表示在原始图像中的坐标位置，相比较取(0,0)点，(0.75,0.75)貌似更接近(1,1)点，那如果将最近邻算法中的取整方式改为cvRound（四舍五入）的方式取(1,1)点，同样会有丢的信息，即丢失了“0.25”部分的(0,0)点、(1,0)点和(0,1)点。

可以看到，dst图像上(X,Y)对应到src图像上的点，最好是根据计算出的浮点数坐标，按照百分比各取四周的像素点的部分。
如下图：

双线性插值的原理相类似，这里不写双线性插值计算点坐标的方法，容易把思路带跑偏，直接就按照比率权重的思想考虑。

将$(\frac{w}{W}X,\frac{h}{H}Y)$写成$(x'+u,y'+v)$的形式，表示将$x$与$y$中的整数和小数分开表示$uv$分别代表小数部分。

这样，根据权重比率的思想得到计算公式

$(X,Y)=(1-u)·(1-v)·(x,y)+(u-1)·v·(x,y+1)+u·(v-1)·(x+1,y)+(u·v)·(x,y)$

在实际的代码编写中，会有两个问题，一个是图像会发生偏移，另一个是效率问题。

几何中心对齐：

由于计算的图像是离散坐标系，如果使用$(\frac{w}{W}X,\frac{h}{H}Y)$公式来计算，得到的(X,Y)值是错误的比率计算而来的，即(x+1,y)、(x,y+1)、(x+1,y+1)这三组点中，有可能有几个没有参与到比率运算当中，或者这个插值的比率直接是错误的。
例如，src图像大小是$a×a$，dst图像的大小是$0.5a×0.5a$。

根据原始公式计算

$(\frac{w}{W},\frac{h}{H})$得到$(2,2)$（注意这不是表示点坐标，而是x和y对应的比率）

如果要计算dst点(0,0)对应的插值结果，由于$(2,2)$是整数，没有小数，所以最后得到dst点在$(0,0)$点的像素值就是src图像上在$(0,0)$点的值。然而，我们想要的dst在$(0,0)$上的结果是应该是有(0,0)(1,0)(0,1)(1,1)这四个点各自按照0.5×0.5的比率加权的结果。
所以我们要将dst上面的点，按照比率$(\frac{w}{W},\frac{h}{H})$向右下方向平移0.5个单位。
公式如下：

$$(x,y)=(X\frac{w}{W}+0.5(\frac{w}{W}-1),Y\frac{h}{H}+0.5(\frac{h}{H}-1))$$

运算优化：

由计算公式可以得知，在计算每一个dst图像中的像素值时会涉及到大量的浮点数运算，性能不佳。

可以考虑将浮点数变换成一个整数，即扩大一定的倍数，运算得到的结果再除以这个倍数。

举一个简单的例子，计算0.25×0.75，可以将0.25和0.75都乘上8，得到2×6=12，结果再除以$8^2$，这样运算的结果与直接计算浮点数没有差别。

在程序中，没有办法取得一个标准的整数，使得两个相互运算的浮点数都变成类似“2”和”6“一样的标准整数，只能取一个适当的值来尽量的减少误差，在源码当中取值为$2^{11}$=2048，即2的固定幂数，最后结果可以通过用位移来表示除以一个2整次幂数，计算速度会有很大的提高。

代码：

void Inter_Linear(const Mat &src, Mat &dst, Size &dsize ,double fx = 0.0, double fy = 0.0)//双线性插值
{
    Size ssize = src.size();//获取矩阵大小
    CV_Assert(ssize.area() > 0);//保证矩阵的长宽都大于0

    if (!dsize.area())//如果dsize为(0,0)
    {
        dsize = Size(saturate_cast<int>(src.cols*fx),//satureate_cast防止数据溢出
            saturate_cast<int>(src.rows*fy));

        CV_Assert(dsize.area());
    }
    else
    {
        fx = (double)dsize.width / src.cols;//Size中的宽高和mat中的行列是相反的
        fy = (double)dsize.height / src.rows;
    }

    dst.create(dsize, src.type());

    double ifx = 1. / fx;
    double ify = 1. / fy;

    uchar* dp = dst.data;
    uchar* sp = src.data;

    int iWidthSrc = src.cols;//宽(列数）
    int iHiehgtSrc = src.rows;//高(行数)
    int channels = src.channels();
    short cbufy[2];
    short cbufx[2];

    for (int row = 0; row < dst.rows; row++)
    {
        float fy = (float)((row + 0.5) * ify - 0.5);
        int sy = cvFloor(fy);//整数部分
        fy -= sy;//小数部分
        sy = std::min(sy, iHiehgtSrc - 2);
        sy = std::max(0, sy);

        cbufy[0] = cv::saturate_cast<short>((1.f - fy) * 2048);
        cbufy[1] = 2048 - cbufy[0];

        for (int col = 0; col < dst.cols; col++)
        {
            float fx = (float)((col + 0.5) * ifx - 0.5);
            int sx = cvFloor(fx);
            fx -= sx;

            if (sx < 0) 
            {
                fx = 0, sx = 0;
            }
            if (sx >= iWidthSrc - 1) 
            {
                fx = 0, sx = iWidthSrc - 2;
            }

            cbufx[0] = cv::saturate_cast<short>((1.f - fx) * 2048);
            cbufx[1] = 2048 - cbufx[0];

            for (int k = 0; k < src.channels(); ++k)
            {
                dp[(row * dst.cols + col) * channels + k] =
                    (
                    sp[ ( sy * src.cols + sx ) * channels + k] * cbufx[0] * cbufy[0] +
                    sp[((sy + 1) * src.cols + sx) * channels + k] * cbufx[0] * cbufy[1] +
                    sp[(sy * src.cols + (sx + 1)) * channels + k] * cbufx[1] * cbufy[0] +
                    sp[((sy + 1) * src.cols + (sx + 1)) * channels + k] * cbufx[1] * cbufy[1]
                    ) >> 22;
            }
        }
    }

}

参考博客

对应OpenCV源码部分为
opencv-2.4.9\modules\imgproc\src\imgwarp 部分