详解opencv resize之INTER_LINEAR和INTER_AREA

一。先简单介绍一下resize的用法

src：输入图，

dst：输出图

dsize：输出图的宽高，如果dsize不为空（即宽高都不是0），则以dsize为准进行resize。

fx, fy是放大缩小的比例，是输出图宽(高)与输入图的宽(高)的比值，这两个都有默认值0，如果想直接指定fx, fy，则可以把dsize置空。

interpolation就是插值的方法了，插值的方法有很多种，本贴只分析INTER_AREA与INTER_LINEAR。

例1：通过dsize指定目标图大小

void resize_demo()
{
    cv::Mat src = cv::imread("../data/car.jpg");
    int height   = src.rows;
    int width    = src.cols;
    cv::Mat dest;
    cv::resize(src, dest, cv::Size(width / 2, height / 4), 0, 0, cv::INTER_LINEAR);
    cv::imwrite("../data/cat_resize.jpg", dest);
}

左边为原图，右边为目标图 

例2：通过fx, fy指定缩放大小，这次都指定为0.5

void resize_demo1()
{
    cv::Mat src = cv::imread("../data/car.jpg");
    int height   = src.rows;
    int width    = src.cols;
    cv::Mat dest;
    cv::resize(src, dest, cv::Size(), 0.5, 0.5, cv::INTER_LINEAR);
    cv::imwrite("../data/cat_resize.jpg", dest);
}

那么目标图上的像素值是怎么算出来的呢？

比如一个5乘5的原图，放大到7乘7目标图上，目标图上的(2,2)位置的像素，映射回原图，x,y坐标都是2*5/7，约为(1.4, 1.4)。如果说正好是(1,1)，那就取原图上(1,1)处的像素值，但现在不是个整数，到底应该怎么取呢？这就涉及到不同的插值算法了。如果我直接取个最近位置的像素值，那就是最近邻算法了(INTER_NEAREST)，比如离(1.4, 1.4)最近的就是(1,1)。

二。INTER_LINEAR

1.插值逻辑

接上面的问题，如果取(1.4, 1.4)周围的4个像素的值，并且按照远近的反比来做一个加权平均，就是双线性插了。

4个像素值分别为v11, v12, v21, v22（名字就跟坐标对应了），中间那个黑点就是v1.4，4个区域的面积分别为a00, a01, a10, a11，它们的和肯定是1。则v1.4=v11*a11+v12*a10+v21*a01+v22*a00。每个像素乘的都是它斜对面的那个面积，其实直观理解也说的通，黑点离a11最近，那a11是不是应该分配一个最大的权重。

具体来分析的话，首先考虑一维场景，在x为1和2的两个点(v1和v2)中间有一个x=1.4的点，如果用v1和v2来计算v1.4，那按照远近关系的反比来求加权平均，就直接明显了，肯定是v1.4=v1*a1+v2*a0

通过x轴方向上的线性插值，就能求出下图两个紫色点的像素值

然后再用它俩在y轴方向上进行线性插值，就能求出最终的v1.4了，所以叫双线性插值。

但是讲到这里还没有结束，刚才并没有严谨地说明坐标值到底是什么含义。一个整数坐标就是一个完整的格子吗？那肯定不是，不然的话0和1之间岂不是没有任何值了。上面在说具体的插值公式时我直接重新画了一个图2，暂时回避了这个问题。

现在再看看图1，v11,v12,v21,v22对应的到底是哪4个点呢？下面的红点分别是原图上的(1.4, 1.4)和目标图上的(2,2)。而4个黑点就是v11,v12,v21,v22

具体的坐标约定参照下图，此图来源于博客：

https://medium.com/@wenrudong/what-is-opencvs-inter-area-actually-doing-282a626a09b3

具体来说，每个像素格子的正中央就是整数坐标点，从(0,0)开始，下图是一维的场景，二维也是一样的道理。那么格子的边框自然就是x.5了，比如0的左边是-0.5，右边是0.5（matlab中的约定跟opencv还不一样，可能matlab中的整数坐标是从1开始的？暂时没去了解matlab）

由上面的约定，我们就能确定目标图上的点到原图上的具体映射公式了，x坐标和y坐标的映射公式是一样的，所以只说x。假设目标图上的像素对应坐标为dx，求其映射回原图的坐标sx。原图的宽度除以目标图的宽度的比值为scale，则：

sx=(dx+0.5)*scale-0.5    ------------------（式1）

这个式子是怎么来的？目标图上的每个点都在原图上有对应点

这是一个线性映射，故可以设sx=a*dx+b，再代入两组具体的点就可以求解此方程：

1.dx=-0.5, sx=-0.5

2.dx=-0.5+1, sx=-0.5+scale （其实也可以不要第2组点，因为a就是scale！）

最终就能求出式1。平时我们直接用cv::resize的时候，可能并不需要知道这个公式，甚至在进行物体检测后处理的时候，我们用目标图上的物体框坐标，直接乘scale得到原图物体框就行了，可能都不用严格按照上式（尤其是scale不太大的时候）。但是如果需要我们自己去实现插值逻辑，比如想在cuda中完成resize，以提升前处理的速度，那就需要用到这个式子了。

2.边界值

现在考虑一下边界情况，比如对于目标图上的(0,0)点，计算出来的原图坐标为(-0.14, -0.14)，那它插值所需要的4个点就变成了(-1,-1),(-1,0), (0,-1),(0,0)，问题是前三个点对应的像素根本就不存在啊。所以直接就用(0,0)这个点的像素值就行了。代码实现上其实仍然是进行插值的，只不过插值前取点的时候，如果发现x坐标小于0了，修正为0，y坐标小于0了，也修正为0。

要说具体是怎么处理的，那肯定是直接上代码最直接，下面截取了一段cv::hal::resize中的代码，这段代码就是在设置目标像素与原图像素的映射关系与插值权重，这是x方向的。y方向的逻辑基本一样。

另外说明一下，这段代码中也包含了INTER_AREA的插值处理，不过是在scale小于1的时候，即放大图片的时候，因为在放大图片的时候，INTER_AREA最终调用的插值代码跟INTER_LINEAR是一样的，只不过它们的插值权重分配有区别。在“三。INTER_AREA”中会详述

    for( dx = 0; dx < dsize.width; dx++ )
    {
        if( !area_mode ) // 非INTER_AREA，比如INTER_LINEAR, INTER_CUBIC，INTER_LANCZOS4
        {
            fx = (float)((dx+0.5)*scale_x - 0.5); // 目标图像素坐标映射到原图的浮点坐标
            sx = cvFloor(fx);
            fx -= sx;
        }
        else
        {
            sx = cvFloor(dx*scale_x); // INTER_AREA的像素坐标映射，向下取整
            fx = (float)((dx+1) - (sx+1)*inv_scale_x); //此式的后半部分是sx的右边像素反映射到目标图上的浮点坐标，记为fx1
            fx = fx <= 0 ? 0.f : fx - cvFloor(fx); //如果fx1在dx+1之前，说明dx其实也映射到sx+1上了，则会取两个像素进行插值。否则直接取sx的像素值
        }

        if( sx < ksize2-1 )
        {
            xmin = dx+1; // 双线性插值逻辑中其实没用到xmin
            if( sx < 0 && (interpolation != INTER_CUBIC && interpolation != INTER_LANCZOS4))
                fx = 0, sx = 0; // 1-fx是左边元素权重，所以fx为0就是只取左边元素，而sx就是左边元素的坐标，所以就完成了左边越界的处理。
        }

        if( sx + ksize2 >= src_width )
        {
            xmax = std::min( xmax, dx ); // 大于等于xmax的dx，只取左边元素
            if( sx >= src_width-1 && (interpolation != INTER_CUBIC && interpolation != INTER_LANCZOS4))
                fx = 0, sx = src_width-1; // 原理其实跟上面一样，右边越界时，还是只取左边元素。注意:sx本身没有可能越右边的界哦。
        }

        for( k = 0, sx *= cn; k < cn; k++ )
            xofs[dx*cn + k] = sx + k; // 目标图像素与原图像素映射关系（对于双线性插值来说就是左边的像素）

        if( interpolation == INTER_CUBIC )
            interpolateCubic( fx, cbuf );
        else if( interpolation == INTER_LANCZOS4 )
            interpolateLanczos4( fx, cbuf );
        else // INTER_LINEAR、INTER_AREA
        {
            cbuf[0] = 1.f - fx; // x轴方向左边元素权重
            cbuf[1] = fx; // 右边元素权重
        }

        // alpha就是x轴方向上每个目标图像素对应到原图的插值权重，对于INTER_LINEAR和INTER_AREA来说都是两个像素，这边只要权重，坐标的映射在前面的xofs里
        if( fixpt ) // 如果depth为CV_8U，则权重用的是short，而非float，最终插值的时候采用整数运算
        {
            for( k = 0; k < ksize; k++ ) // 先设第0个通道的权重
                ialpha[dx*cn*ksize + k] = saturate_cast<short>(cbuf[k]*INTER_RESIZE_COEF_SCALE); 
            for( ; k < cn*ksize; k++ ) //其它通道的权重分配直接抄第0个通道就行了
                ialpha[dx*cn*ksize + k] = ialpha[dx*cn*ksize + k - ksize];
        }
        else
        {
            for( k = 0; k < ksize; k++ )
                alpha[dx*cn*ksize + k] = cbuf[k];
            for( ; k < cn*ksize; k++ )
                alpha[dx*cn*ksize + k] = alpha[dx*cn*ksize + k - ksize];
        }
    }

3.适用场景

平时就用INTER_LINEAR就好啦，那为什么下面还要介绍INTER_AREA呢，其实它两适用的场景是不一样的，放大图片的时候用INTER_LINEAR比较好（如果想更好可以用INTER_CUBIC，但是计算量更大），但是缩小图片就不一定了。

比如看下面这个例子，6乘6的图，缩小为2乘2，scale为3

    scale = 3
    print([(dx+0.5)*scale-0.5 for dx in range(2)])

    img = np.array([[86, 85, 83, 92, 104, 127],
                    [113, 105, 92, 103, 118, 135],
                    [156, 137, 105, 121, 141, 147],
                    [144, 138, 129, 148, 162, 130],
                    [126, 135, 149, 168, 176, 108],
                    [133, 132, 130, 138, 140, 85]], dtype=np.uint8)
    print(cv2.resize(img, (2, 2), interpolation=cv2.INTER_LINEAR))

输出结果如下：

[1.0, 4.0]   ----这就是由dx为0、1时计算出的sx，得到两个整数，所以不会用到4个像素的值
[[105 118]       ----分别对应原图的(1,1), (1, 4)
 [135 176]]     -----分别对应原图的(4,1), (4, 4)

映射图如下，目标图上每个像素都只用到原图上一个像素，效果自然不会太好。

如果说刚才只是一个巧合，其实整体来说，缩小图片时，如果缩小的比例比较大（即原图尺度除以目标图的比例），用INTER_LINEAR的效果就会打折扣。

比如现在是8乘8的图，缩放到2乘2，如果仍然用INTER_LINEAR，则4个点映射关系如下，原图和目标图的对应点分别用了4种颜色，每个点都处在4个像素的边框界限上，故其4周的4个像素正好可以取平均值，也就是说是正常插值的。但是再想想看，目标图是2乘2，每一个像素其实对应在原图上的区域是4乘4的区域（即16个像素），但插值只用了4个像素。所以说没有充分利用到原图上的像素。

那我怎样才能用到原图上的所有像素呢？就像上面的4个颜色的框子一样？那就是下面要介绍的INTER_AREA了。

三。INTER_AREA

 INTER_AREA的映射逻辑跟INTER_LINEAR不一样，INTER_LINEAR是点坐标映射，而 INTER_AREA更像是像素映射。并且INTER_AREA的代码分了几种情况进行处理。

1.整数倍缩放

其实图9就是一个整数倍缩放的例子，比如原宽高是目标图的iscale_x和iscale_y倍，则每目标图中每一个像素必然对应原图中的iscale_x*iscale_y个像素，并且目标图上不同像素在原图中的映射不会重叠，且总共会覆盖到整个原图。这很好理解，代码处理起来也很简单。最终把iscale_x*iscale_y个像素取个平均值就是缩放后的像素值，这就是INTER_AREA的fast流程。

这里再提一个特例，就是如果iscale_x和iscale_y都是2的话，opencv会调用fast的fast流程~~，因为4个像素取平均值通过移位操作就可以完成了，并且还用到了SIMD（单指令多数据），这个我暂时没研究~~。另外，如果是INTER_LINEAR，并且也满足iscale_x和iscale_y都是2的条件，那么实际上它的插值逻辑等阶于INTER_AREA，所以它也会直接走INTER_AREA的fast的fast流程。

这里你可能有个疑问，我怎么通过代码判断是不是整数倍缩放呢？代码如下，scale_x, scale_y是原图宽高除以目标图的宽高，类型为double，saturate_cast即cvRound，应该就类似于round()。

    int iscale_x = saturate_cast<int>(scale_x);
    int iscale_y = saturate_cast<int>(scale_y);

    bool is_area_fast = std::abs(scale_x - iscale_x) < DBL_EPSILON &&
            std::abs(scale_y - iscale_y) < DBL_EPSILON;

而DBL_EPSILON是double数的最小间隔数值，在我的gcc头文件里看到的值为2.22044604925031308084726333618164062e-16L，注意，它可不是最小的double正值，它是不带指数位的情况下可表示的两个连续的double数之间的最小差值（或者说如下图的注释所说，是比1大的最小数与1的间隔）。要理解这个，需要先了解IEEE浮点标准。

 所以DBL_EPSILON就是2的负52次方,即2.22E-16。那此处我们为什么通过它来判断是否是整数呢？直接用DBL_MIN行不行呢？我也不能完全说上来【狗头】。这个问题未必那么重要，如果scale和真正的整数已经是这么小的差异了，当成整数完全是没有问题的，因为在”2.非整数倍缩放“中你会发现，这么小的差异已经被忽略了。

不过此处还是插入一段浮点数精确表示整数的分析（不感兴趣的童鞋请跳过）。用下面这段代码来看看float32是不是能表示所有的INT32。

void print_float(float a)
{
    unsigned int *ai = (unsigned int*)&a;
    unsigned int a_sign = *ai >> 31;
    unsigned int a_exp = (*ai >> 23) & ((1 << 8) - 1);
    unsigned int a_base = *ai & ((1 << 23) - 1);
    printf("a: %f, a_sign: %X, a_exp: %X, a_base: %X\n", a, a_sign, a_exp, a_base);

}


int main(int argc, char* argv[])
{
    float a;

    int error_count = 0;
    for (int i = 0; i < INT_MAX; i++)
    {
        a = i;
        if (std::abs(a - i) > __FLT_EPSILON__ || std::abs(i - (int)a > 0)) // 为什么多加了一个条件，因为根本不存在16777217这个浮点数
        // if (std::abs(a - i) > FLT_MIN || std::abs(i - (int)a > 0)) // 其实是一样的结果
        {
            printf("i: %d, ", i);
            print_float(a);

            if (error_count > 10)
            {
                break;
            }
            error_count++;
        }
    }

    return 0;
}

运行结果如下： 

可以看到，从16777217开始，就不都能由float32精确表示了，其实原因很简单 16777216的时候，由上图可以看到，指数位为0x97，即151，减去指数偏移值127，得到24。尾数位为0，故最终的浮点值为1.0*(2^24)=16777216，其实2的24次方就是16777216。（为什么是1.0乘？因为IEEE标准默认尾数位是从1.X开始表示的，即小数点前面一定是1，这个1不用占用尾数位）。此时再让它增加一个最小值，只能让尾数最右位加1，最终整个浮点数加的值为：(2^-23)*(2^24)=2。所以只能得到16777218，跳过了16777217这个数。

所以就明白为什么16777216之前的整数都能精确表示出来了，因为在此之前，指数位m小于24，我始终能让下一个数比上一个数大1，比如我在尾数位第n位加1（从左往右数的第n位），那么最终增加的值就是(2^-n)*(2^m)，我只要让n==m，加的值就是1！，而只要m小于24，我都能找到这个n。

不过有一点要注意的是，并不是说因为float32能表示16777216之前的所有整数，你就能用float32精确计算出16777216之前的所有scale，因为分子会比16777216更大啦。但是实际情况下往往也没这么大的图，也没这么大的scale，所以float压力都不大，double更加不是问题啦！

2.非整数倍缩放

非整数倍的缩放其实也是很简单的，比如3乘3缩到2乘2，那目标图上每个像素在原图上会分到3/2乘3/2个像素，仍然是瓜分全图的。不管是什么比例都能瓜分全图，比如m乘m缩放到n乘n，那每个分m/n乘m/n，n*n*(m/n)*(m/n)=m*m。至于原图各像素所占的权重，其实就是映射到它上面的面积。

可以看一个跨更多像素的例子，这里是12乘12的图缩放到5乘5，缩放比例2.4，则目标图上的(1,1)像素，对应原图范围为(2.4,2.4)到(4.8,4.8)的矩形范围，一共占了9个像素，蓝色的那个是完整覆盖，其它的8个都只覆盖了部分面积，覆盖的9个像素的总面积为A（即2.4乘2.4），则计算目标像素时，每个原像素所占的权重就是它被覆盖到的面积除以总面积A。

这里可能有个疑问，2.4为什么是下图中的橙色点，而不是红色点的位置？

目标图上的(0,0)对应的肯定是原图上的(0,0)到(2.4,2.4)的范围吧，2.4就是scale，这样就想明白了，所以说其实对于INTER_AREA来说，如果考虑浮点坐标，则左上角的起始坐标为(0,0)，右下角的坐标为(12,12)，前闭后开，(0,0)坐标能取到，(12,12)的坐标实际取不到。

3.放大

其实放大的插值原理跟“2.非整数倍缩放”是一样的，但是既然是放大，目标图在原图上的覆盖面积必然小于1，它有可能只会覆盖到一个像素（覆盖面在一个像素内），最多覆盖4个像素（覆盖面在4个像素交界处），这似乎跟双线性插值有一点像，都是4个像素按权重求平均麻，所以说INTER_AREA的放大流程最终跟INTER_LINEAR调的是相同的函数，只是权重分配逻辑有区别，关于权重分配和像素映射的代码参考“二。INTER_LINEAR”中的那段。

我们只把最关键的那段代码再抄一遍：

sx = cvFloor(dx*scale_x); // INTER_AREA的像素坐标映射，向下取整
            
//此式的后半部分是sx的右边像素反映射到目标图上的浮点坐标，记为fx1
fx = (float)((dx+1) - (sx+1)*inv_scale_x); 

//如果fx1在dx+1之前，说明dx其实也映射到sx+1上了，则会取两个像素进行插值。否则直接取sx的像素值
fx = fx <= 0 ? 0.f : fx - cvFloor(fx); 

比如现在是5乘5的图，放大到12乘12，scale为5/12，inv_scale就是12/5。为了方便展示，只画x轴方向的图。下图演示的是目标图dx=2像素，映射回原图覆盖的区域，映射回去的浮点值是0.83，向下取整得到sx为0，然后把sx+1反映射到目标图上（就是绿色那条线），得到2.4，然后用dx+1(即2+1)减2.4，得到最终的fx=0.6，它就是原图上sx+1所分配的权重，sx分配的权重自然就是1-fx=0.4。也就是说opencv如上的代码，它其实是通过目标图上的AB、BC段的比例来计算权重。不过其实AB和BC的比例与ab和bc的比例是一样的，似乎在原图上直接计算权重也没问题。

 当然喽，不是所有目标像素在原图上都能覆盖两个像素（仅考虑x轴方向最多是2个像素，如果是二维的就是最多4个像素啦），比如考虑dx=1的像素，会发现C跑到B前面去了，那C-B就得到负值，按如上代码的逻辑，fx取0，所以sx+1的权重为0。直观理解也很明显，ac段就是在原图上覆盖的区域，只覆盖到了sx=0，自然只会取该像素值。

 再回顾一下刚才说直接在原图上通过ab、bc的比例计算权重似乎也行，真的是这样吗？其实未必，因为在目标图上，直接计算BC就能得出fx，不需要再除什么东西，因为AC的长度就是一个像素的长度，分母是1。但如果在原图上计算，ac必然不是1，因为它就是scale啊，而放大的时候，scale必然小于1，所以得用bc/scale，bc经过了浮点计算，再拿它除浮点数，有可能造成误差放大。

下面结合代码再看一下，fx_opencv就是opencv中在目标图上计算权重的方法，fx_my就是在原图上计算权重，可以看到它最后要多除一个scale_x。

float fx_opencv(int dx, double scale_x, double inv_scale_x)
{
    int sx = cvFloor(dx * scale_x);
    float fx = (float)((dx + 1) - (sx + 1) * inv_scale_x);
    fx = fx <= 0 ? 0.f : fx - cvFloor(fx);
    return fx;
}

float fx_my(int dx, double scale_x)
{
    float fx = dx * scale_x;
    int sx = cvFloor(fx);
    fx = fx + scale_x - (sx + 1);
    if (fx < 0)
    {
        fx = 0.f;
    }
    else
    {
        fx /= scale_x;
    }
    return fx;
}

int main(int argc, char* argv[])
{
    int src_width = 5, dest_width = 12;
    double scale = src_width / (double)dest_width;
    double inv_scale = 1 / scale;
    for (int i = 0; i < dest_width; i++)
    {
        printf("i: %d, fx_opencv: %f, fx_my: %f\n", i, fx_opencv(i, scale, inv_scale), fx_my(i, scale));
    }
    return 0;
}

这次的结果倒是一样，不过如果也许除法不如乘法快？暂时就不研究那么多了~~

另外如果是一切更特殊的数字，比如13放大到17，那结果还是有点差别的，不过也不大