詳解opencv resize之INTER_LINEAR和INTER

一。先簡單介紹一下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的右邊像素反映射到目標圖上的浮點坐標，記為fx1fx = fx <= 0 ? 0.f : fx - cvFloor(fx); //如果fx1在dx+1之前，說明dx其實也映射到sx+1上了，則會取兩個像素進行插值。否則直接取sx的像素值}if( sx < ksize2-1 ){xmin = dx+1; // 雙線性插值邏輯中其實沒用到xminif( 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 = 3print([(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;
}