1 雙邊濾波簡介
雙邊濾波（Bilateral filter）是一種非線性的濾波方法，是結合圖像的空間鄰近度和像素值相似度的一種折衷處理，同時考慮空域信息和灰度相似性，達到保邊去噪的目的。具有簡單、非迭代、局部的特點。
雙邊濾波器的好處是可以做邊緣保存（edge preserving），一般過去用的維納濾波或者高斯濾波去降噪，都會較明顯地模糊邊緣，對於高頻細節的保護效果並不明顯。雙邊濾波器顧名思義比高斯濾波多了一個高斯方差sigma－d，它是基於空間分布的高斯濾波函數，所以在邊緣附近，離的較遠的像素不會太多影響到邊緣上的像素值，這樣就保證了邊緣附近像素值的保存。但是由於保存了過多的高頻信息，對於彩色圖像里的高頻噪聲，雙邊濾波器不能夠干凈的濾掉，只能夠對於低頻信息進行較好的濾波。
2 雙邊濾波原理
濾波算法中，目標點上的像素值通常是由其所在位置上的周圍的一個小局部鄰居像素的值所決定。在2D高斯濾波中的具體實現就是對周圍的一定范圍內的像素值分別賦以不同的高斯權重值，並在加權平均后得到當前點的最終結果。而這里的高斯權重因子是利用兩個像素之間的空間距離（在圖像中為2D）關系來生成。通過高斯分布的曲線可以發現，離目標像素越近的點對最終結果的貢獻越大，反之則越小。其公式化的描述一般如下所述：
其中的c即為基於空間距離的高斯權重，而用來對結果進行單位化。
高斯濾波在低通濾波算法中有不錯的表現，但是其卻有另外一個問題，那就是只考慮了像素間的空間位置上的關系，因此濾波的結果會丟失邊緣的信息。這里的邊緣主要是指圖像中主要的不同顏色區域（比如藍色的天空，黑色的頭發等），而Bilateral就是在Gaussian blur中加入了另外的一個權重分部來解決這一問題。Bilateral濾波中對於邊緣的保持通過下述表達式來實現：
其中的s為基於像素間相似程度的高斯權重，同樣用來對結果進行單位化。對兩者進行結合即可以得到基於空間距離、相似程度綜合考量的Bilateral濾波：
上式中的單位化分部綜合了兩種高斯權重於一起而得到，其中的c與s計算可以詳細描述如下：
且有
且有
上述給出的表達式均是在空間上的無限積分，而在像素化的圖像中當然無法這么做，而且也沒必要如此做，因而在使用前需要對其進行離散化。而且也不需要對於每個局部像素從整張圖像上進行加權操作，距離超過一定程度的像素實際上對當前的目標像素影響很小，可以忽略的。限定局部子區域后的離散化公就可以簡化為如下形式：
上述理論公式就構成了Bilateral濾波實現的基礎。為了直觀地了解高斯濾波與雙邊濾波的區別，我們可以從下列圖示中看出依據。假設目標源圖像為下述左右區域分明的帶有噪聲的圖像（由程序自動生成），藍色框的中心即為目標像素所在的位置，那么當前像素處所對應的高斯權重與雙邊權重因子3D可視化后的形狀如后邊兩圖所示：
左圖為原始的噪聲圖像；中間為高斯采樣的權重；右圖為Bilateral采樣的權重。從圖中可以看出Bilateral加入了相似程度分部以后可以將源圖像左側那些跟當前像素差值過大的點給濾去，這樣就很好地保持了邊緣。為了更加形象地觀察兩者間的區別，使用Matlab將該圖在兩種不同方式下的高度圖3D繪制出來，如下：
上述三圖從左到右依次為：雙邊濾波，原始圖像，高斯濾波。從高度圖中可以明顯看出Bilateral和Gaussian兩種方法的區別，前者較好地保持了邊緣處的梯度，而在高斯濾波中，由於其在邊緣處的變化是線性的，因而就使用連累的梯度呈現出漸變的狀態，而這表現在圖像中的話就是邊界的丟失（圖像的示例可見於后述）。
3 雙邊濾波實現步驟
有了上述理論以后實現Bilateral Filter就比較簡單了，其實它也與普通的Gaussian Blur沒有太大的區別。這里主要包括3部分的操作: 基於空間距離的權重因子生成；基於相似度的權重因子的生成；最終filter顏色的計算。
3.1Spatial Weight
這就是通常的Gaussian Blur中使用的計算高斯權重的方法，其主要通過兩個pixel之間的距離並使用如下公式計算而來：
3.1Similarity Weight
與基於距離的高斯權重計算類似，只不過此處不再根據兩個pixel之間的空間距離，而是根據其相似程度（或者兩個pixel的值之間的距離）。
其中的表示兩個像素值之間的距離，可以直接使用其灰度值之間的差值或者RGB向量之間的歐氏距離。
3.4 Color Filtering
其中的相似度分部的權重s主要根據兩個Pixel之間的顏色差值計算面來。對於灰度圖而言，這個差值的范圍是可以預知的，即[-255, 255]，因而為了提高計算的效率我們可以將該部分權重因子預計算生成並存表，在使用時快速查詢即可。使用上述實現的算法對幾張帶有噪聲的圖像進行濾波后的結果如下所示：
3 雙邊濾波源碼實現
1 #include <time.h> 2 #include <stdio.h> 3 #include <stdlib.h> 4 #include <math.h> 5 #include "bf.h" 6 7 8 //-------------------------------------------------- 9 /** 10 雙邊濾波器圖像去噪方法 11 \param pDest 去噪處理后圖像buffer指針 12 \param pSrc 原始圖像buffer指針 13 \paran nImgWid 圖像寬 14 \param nImgHei 圖像高 15 \param nSearchR 搜索區域半徑 16 \param nSigma 噪聲方差 17 \param nCoff DCT變換系數個數 18 19 return void 20 */ 21 //-------------------------------------------------- 22 void BilateralFilterRGB(float *pDest, BYTE *pSrc,int nImgWid, int nImgHei, int nSearchR, int nSigma, int nCoff) 23 { 24 25 BYTE *pSrcR = new BYTE[nImgHei * nImgWid]; // 新建圖像緩沖區R 26 BYTE *pSrcG = new BYTE[nImgHei * nImgWid]; // 新建圖像緩沖區G 27 BYTE *pSrcB = new BYTE[nImgHei * nImgWid]; // 新建圖像緩沖區B 28 29 30 for(int i = 0;i < nImgHei; i++) 31 { 32 for (int j = 0; j < nImgWid; j++) 33 { 34 pSrcR[i * nImgWid + j] = pSrc[i * nImgWid * 3 + j * 3 + 0]; 35 pSrcG[i * nImgWid + j] = pSrc[i * nImgWid * 3 + j * 3 + 1]; 36 pSrcB[i * nImgWid + j] = pSrc[i * nImgWid * 3 + j * 3 + 2]; 37 } 38 } 39 40 41 float* pDestR = new float[nImgHei * nImgWid]; // 去噪后的圖像數據R通道 42 float* pDestG = new float[nImgHei * nImgWid]; // 去噪后的圖像數據G通道 43 float* pDestB = new float[nImgHei * nImgWid]; // 去噪后的圖像數據B通道 44 memset(pDestR,255,nImgHei * nImgWid * sizeof(float)); // 傳遞變量初始化 45 memset(pDestG,255,nImgHei * nImgWid * sizeof(float)); 46 memset(pDestB,255,nImgHei * nImgWid * sizeof(float)); 47 48 49 float *pII = new float[nImgHei * nImgWid ]; // 積分圖 50 float *pWeights = new float[nImgHei * nImgWid ]; // 每一個窗的權重之和，用於歸一化 51 52 53 float mDctc[MAX_RANGE]; // DCT變換系數 54 float mGker[MAX_RANGE]; // 高斯核函數的值 55 56 for (int i = 0; i< MAX_RANGE; ++i) 57 { 58 mGker[i] = exp(-0.5 * pow(i / nSigma, 2.0)); // 首先計算0-255的灰度值在高斯變換下的取值，然后存在mGker數組中 59 } 60 61 CvMat G = cvMat(MIN_RANGE, MAX_RANGE, CV_32F, mGker); // 轉換成opencv中的向量 62 CvMat D = cvMat(MIN_RANGE, MAX_RANGE, CV_32F, mDctc); // 把mDctc系數數組轉換成opencv中的向量形式 63 64 cvDCT(&G,&D,CV_DXT_ROWS); // 然后將高斯核進行離散的dct變換 65 mDctc[0] /= sqrt(2.0); // 離散余弦變換的第一個系數是1/1.414; 66 67 68 bf(pSrcR, pDestR, mDctc, pII, pWeights, nImgHei, nImgWid, nCoff, nSearchR); 69 bf(pSrcG, pDestG, mDctc, pII, pWeights, nImgHei, nImgWid, nCoff, nSearchR); 70 bf(pSrcB, pDestB, mDctc, pII, pWeights, nImgHei, nImgWid, nCoff, nSearchR); 71 72 73 for(int i=0; i < nImgHei;i++) 74 { 75 for (int j=0; j < nImgWid; j++) 76 { 77 pDest[i * nImgWid * 3 + j * 3 + 0] = pDestR[i * nImgWid + j]; 78 pDest[i * nImgWid * 3 + j * 3 + 1] = pDestG[i * nImgWid + j]; 79 pDest[i * nImgWid * 3 + j * 3 + 2] = pDestB[i * nImgWid + j]; 80 } 81 } 82 83 84 free(pSrcR); // 釋放臨時緩沖區 85 free(pSrcG); 86 free(pSrcB); 87 free(pDestR); 88 free(pDestG); 89 free(pDestB); 90 91 }
1 //-------------------------------------------------- 2 /** 3 雙邊濾波器圖像去噪方法 4 \param pDest 去噪處理后圖像buffer指針 5 \param pSrc 原始圖像buffer指針 6 \paran nImgWid 圖像寬 7 \param nImgHei 圖像高 8 \param pII 積分圖緩沖區 9 \param pWeights 歸一化權重系數 10 \param nCoff DCT變換系數 11 \param nPatchWid 圖像塊寬度 12 13 return void 14 */ 15 //-------------------------------------------------- 16 void bf (BYTE *pSrc, float *pDest, float *pDctc, float *pII, float *pWeights, int nImgHei, int nImgWid, int nCoff, int nPatchWid) 17 { 18 if (pDest == NULL || pSrc ==NULL) 19 { 20 return; 21 } 22 23 float *cR = (float*) malloc((nCoff - 1) * MAX_RANGE * sizeof(float)); 24 float *sR = (float*) malloc((nCoff - 1) * MAX_RANGE * sizeof(float)); 25 float *dcR = (float*) malloc((nCoff - 1) * MAX_RANGE * sizeof(float)); 26 float *dsR = (float*) malloc((nCoff - 1) * MAX_RANGE * sizeof(float)); 27 28 29 // 余弦變換查找表 30 int fx = 0; 31 int ck = 0; 32 int r2 = pow(2.0*nPatchWid + 1.0, 2.0); // 這個是圖像塊的像素點的總個數 33 34 float c0 = pDctc[0]; // 這是DCT的第一個系數 35 float c0r2 = pDctc[0] * r2; // 這個用於初始化歸一化權重之和 36 37 38 for (ck = 1; ck < nCoff; ck++) // 注意到這里面的系數並不包含C0，所以后面才有把C0額外的加上 39 { // 一個矩陣，用於存放各種系數和數值，便於查找 40 int ckr = (ck - 1) * MAX_RANGE; // 數組是從0開始的，所以這里減1 41 42 for (fx = 0; fx<MAX_RANGE; fx++) // fx就是圖像的灰度值 43 { 44 float tmpfx = PI * float(fx) * float(ck) / MAX_RANGE; // ck其實就相當於那個余弦函數里面的t，fx相當於u，PI/MAX相當於前面的那個系數，這都是余弦變換相關的東西 45 46 cR[ckr + fx] = cos(tmpfx); // 存儲余弦變換，這個用在空間域上 47 sR[ckr + fx] = sin(tmpfx); // 存儲正弦變換 48 dcR[ckr + fx] = pDctc[ck] * cos(tmpfx); // 存儲余弦變換和系數的乘積，這個則用在強度范圍上 49 dsR[ckr + fx] = pDctc[ck] * sin(tmpfx); // 存儲正弦變換和系數的乘積 50 } 51 } 52 53 54 float *pw = pWeights; // 新建一個歸一化權重的中間變量進行數據的初始化 55 float *pwe = pWeights + nImgWid * nImgHei; // 限定最大位置 56 57 while (pw < pwe) 58 { 59 *pw++ = c0r2; // 賦初值，讓它們都等於第一個DCT系數乘以圖像塊中總的像素點的個數 60 } 61 62 for (int ck = 1; ck < nCoff; ck++) 63 { 64 int ckr = (ck-1)*MAX_RANGE; // 數組是從0開始的，所以這里減1 65 66 add_lut(pWeights, pSrc, pII,cR + ckr, dcR + ckr, nImgHei, nImgWid, nPatchWid, nPatchWid); // add cos term to pWeights 67 add_lut(pWeights, pSrc, pII,sR + ckr, dsR + ckr, nImgHei, nImgWid, nPatchWid, nPatchWid); // add cos term to pWeights 68 69 } 70 71 ImgRectSumO(pDest, pSrc, pII, c0, nImgHei, nImgWid, nPatchWid, nPatchWid); //加上C0的變換之后，再初始化濾波后的函數 72 73 for (int ck = 1; ck < nCoff; ck++) 74 { 75 int ckr = (ck-1)*MAX_RANGE; 76 77 add_f_lut(pDest, pSrc, pII, cR + ckr, dcR + ckr, nImgHei, nImgWid, nPatchWid, nPatchWid); // add cosine term to dataf 78 add_f_lut(pDest, pSrc, pII, sR + ckr, dsR + ckr, nImgHei, nImgWid, nPatchWid, nPatchWid); // add sine term to dataf 79 80 } 81 82 83 float *pd = pDest + nPatchWid * nImgWid; 84 pw = pWeights + nPatchWid * nImgWid; 85 86 float *pdie = pd + nImgWid - nPatchWid; 87 float *pdend = pDest + (nImgHei - nPatchWid) * nImgWid; 88 89 while (pd < pdend) 90 { 91 pd += nPatchWid; //把邊都給去掉了 92 pw += nPatchWid; //這個也是把邊去掉了 93 94 while (pd < pdie) 95 { 96 *pd++ /= ( (*pw++)*255); // 之所以要除以255，是為了把它化到[0,1]之間，便於顯示 97 } 98 99 pd += nPatchWid; // 把邊都給去掉了 100 pw += nPatchWid; // 這個也是把邊去掉了 101 pdie += nImgWid; // 過渡到下一行 102 } 103 104 free(cR); 105 free(sR); 106 free(dcR); 107 free(dsR); // 釋放緩沖區 108 }
//-------------------------------------------------- /** 余弦函數首系數積分圖 \param pDest 去噪處理后圖像buffer指針 \param pSrc 原始圖像buffer指針 \paran nImgWid 圖像寬 \param nImgHei 圖像高 \param pII 積分圖緩沖區 \param c0 DCT變換第一個系數 \param nPatchWid 圖像塊寬度 \param nPatchHei 圖像塊高度 return void */ //-------------------------------------------------- void ImgRectSumO (float* pDest, const BYTE* pSrc, float* pII, float c0,const int nImgHei, const int nImgWid, const int nPatchWid, const int nPatchHei) { if (pDest == NULL || pSrc ==NULL) { return; } int ri1 = nPatchWid + 1; // 中心像素點 int rj1 = nPatchHei + 1; int ri21 = 2 * nPatchWid + 1; // 圖像塊的寬度 int rj21 = 2 * nPatchHei + 1; const BYTE *pSrcImg = pSrc; // 原始圖像數據 const BYTE *pSrcImgEnd = pSrc + nImgHei * nImgWid; // 最大的像素點位置 const BYTE *pSrcImgWid = pSrcImg + nImgWid; // 第一行的最大的像素點位置，下面循環用的的變量 float *pIITmp = pII; // 積分圖數據 float *ppII = pIITmp; // 積分圖中間變量 *pIITmp++ = c0 * float( *pSrcImg++ ); // 第一個系數乘以圖像數據 while (pSrcImg < pSrcImgWid) // 遍歷第一行進行求積分圖，其實這個是圖像數據的積分圖 { *pIITmp++ = (*ppII++) + c0 * float(*pSrcImg++); } while (pSrcImg < pSrcImgEnd) // 遍歷所有的像素點的變量 { pSrcImgWid += nImgWid; // 進行第二行的循環 ppII = pIITmp; // 積分圖中間變量 *pIITmp++ = c0 * float(*pSrcImg++); // 第二行第一個像素點的處理 while (pSrcImg < pSrcImgWid) { (*pIITmp++) = (*ppII++) + c0 * float(*pSrcImg++); // 求第二行的積分圖 } float *pIIWid = pIITmp; // 當前位置像素點 pIITmp = pIITmp - nImgWid; // 上一行的起始點位置 float *pIIup = pIITmp - nImgWid; // 上上一行的起始點位置 while (pIITmp < pIIWid) // 遍歷整個行的每一個數據 { (*pIITmp++) = (*pIITmp) + (*pIIup++); // 行與行之間的積分圖 } } float *pres = pDest + ri1 * nImgWid; // 最小的行位置 float *pend = pDest + (nImgHei - nPatchWid) * nImgWid; // 最大的行位置 float *pii1 = NULL; float *pii2 = NULL; float *pii3 = NULL; float *pii4 = NULL; // 求積分圖的四個點 pii1 = pII + ri21 * nImgWid + rj21; // 右下角 pii2 = pII + ri21 * nImgWid; // 左下角 pii3 = pII + rj21; // 右上角 pii4 = pII; // 左上角 while (pres < pend) { float *pe = pres + nImgWid - nPatchHei; pres += rj1; while (pres < pe) // 這個只是求圖像數據的積分圖 { (*pres++) = (*pii1++) - (*pii2++) - (*pii3++) + (*pii4++); } pres += nPatchHei; pii1 += rj21; pii2 += rj21; pii3 += rj21; pii4 += rj21; } }
//-------------------------------------------------- /** 余弦積分圖加函數 \param pNomalWts 像素點的歸一化權重矩陣 \param pSrc 原始圖像buffer指針 \param pII 積分圖 \paran nImgWid 圖像寬 \param nImgHei 圖像高 \param pII 積分圖緩沖區 \param pCosMtx 余弦函數矩陣 \param pCosDctcMtx 余弦函數與DCT變換系數乘積矩陣 \param nPatchWid 圖像塊寬度 \param nPatchHei 圖像塊高度 return void */ //-------------------------------------------------- void add_lut (float* pNomalWts, const BYTE* pSrc, float* pII, float* pCosMtx, float* pCosDctcMtx, const int nImgHei, const int nImgWid, const int nPatchWid, const int nPatchHei) { int ri1 = nPatchWid + 1; // kernel相關 int rj1 = nPatchHei + 1; int ri21 = 2 * nPatchWid + 1; // 這是kernel的寬度和長度 int rj21 = 2 * nPatchHei + 1; const BYTE *pi = pSrc; // 這個是圖像數據 float *pii = pII; // 這個是中間的緩沖變量，積分圖的緩沖區 const BYTE *piw = pi + nImgWid; // 也是賦初值的時候使用的中間變量 float *pii_p = pii; // 直線積分圖指針的指針 *pii++ = pCosMtx[*pi++]; //圖像數據值所對應的查找表中的余弦或者正弦變換 while (pi < piw) { *pii++ = (*pii_p++) + pCosMtx[*pi++]; // 第一行的積分圖 } const BYTE *piend = pSrc + nImgHei * nImgWid; // 限定循環位置 while (pi < piend) // 遍歷整個圖像數據 { piw += nImgWid; // 第二行 pii_p = pii; // 指向積分圖指針的指針 *pii++ = pCosMtx[*pi++]; // 獲取第一個圖像數據值所對應的查找表中的余弦或者正弦變換 while (pi < piw) { (*pii++) = (*pii_p++) + pCosMtx[*pi++]; // 求這一行的積分圖 } float *piiw = pii; pii = pii - nImgWid; float *pii_p1 = pii - nImgWid; while (pii < piiw) { (*pii++) = (*pii) + (*pii_p1++); // 行與行之間求積分圖 } } float *pres = pNomalWts + ri1 * nImgWid; // 定位要處理點的像素的那一行 float *pend = pNomalWts + (nImgHei - nPatchWid) * nImgWid; // 限定了邊界 float *pii1 = NULL; float *pii2 = NULL; float *pii3 = NULL; float *pii4 = NULL; pii1 = pII + ri21 * nImgWid + rj21; // 獲得積分圖中那個最大的數據右下角 pii2 = pII + ri21 * nImgWid; // 積分圖左下角 pii3 = pII + rj21; // 積分圖右上角 pii4 = pII; // 積分圖左上角 pi = pSrc + ri1 * nImgWid; // 定位要處理的像素的位置 while (pres < pend) // 設定高度做循環 { float *pe = pres + nImgWid - nPatchHei; // 限定了寬度的范圍 pres = pres + rj1; // 定位了要處理的像素點的歸一化系數矩陣的位置 pi = pi + rj1; // 定位了要處理的像素點 while (pres < pe) // 設定寬度做循環 { (*pres++) = (*pres) + pCosDctcMtx[*pi++] * ( (*pii1++) - (*pii2++) - (*pii3++) + (*pii4++) ); // 得到的就是歸一化系數矩陣之和 } pres += nPatchHei; pi += nPatchHei; pii1 += rj21; pii2 += rj21; pii3 += rj21; pii4 += rj21; // 略過不處理的邊界區域 } }
//-------------------------------------------------- /** 積分圖 \param pDest 像素點的歸一化權重矩陣 \param pSrc 原始圖像buffer指針 \param pII 積分圖 \paran nImgWid 圖像寬 \param nImgHei 圖像高 \param pII 積分圖緩沖區 \param pCosMtx 余弦函數矩陣 \param pCosDctcMtx 余弦函數與DCT變換系數乘積矩陣 \param nPatchWid 圖像塊寬度 \param nPatchHei 圖像塊高度 return void */ //-------------------------------------------------- void add_f_lut (float* pDest,const BYTE* pSrc,float* pII,float* pCosMtx,float* pCosDctcMtx,const int nImgHei, const int nImgWid,const int nPatchWid, const int nPatchHei) { int ri1 = nPatchWid + 1; // kernel相關，目標是指向中心像素點 int rj1 = nPatchHei + 1; // 目標是指向中心像素點 int ri21 = 2 * nPatchWid + 1; // 這是kernel的寬度和長度 int rj21 = 2 * nPatchHei + 1; // 其實就是指向搜索區域的右下角，也就是最大位置的那個像素點，就是邊界了 const BYTE *pi = pSrc; // 這個是圖像數據,圖像的灰度值 float *pii = pII; // 這是積分圖 const BYTE *piw = pi + nImgWid; // 指向第一行的末尾位置，用於循環控制 float *pii_p = pii; // 指向積分圖指針的指針 *pii++ = float(*pi) * pCosMtx[*pi]; // 灰度值乘以余弦系數 ++pi; while (pi < piw) // 第一行遍歷 { *pii++ = (*pii_p++) + float(*pi) * pCosMtx[*pi]; ++pi; } const BYTE *piend = pSrc + nImgHei * nImgWid; while (pi < piend) { piw += nImgWid; pii_p = pii; *pii++ = float(*pi) * pCosMtx[*pi]; ++pi; while (pi < piw) { (*pii++) = (*pii_p++) + float(*pi) * pCosMtx[*pi]; ++pi; } float *piiw = pii; pii = pii - nImgWid; // 上一行起點 float *pii_p1 = pii - nImgWid; // 上上一行的起始點 while (pii < piiw) // 循環一行 { (*pii++) += (*pii_p1++); //其實就是在列的位置上加上上一行 } } float *pres = pDest + ri1*nImgWid; float *pend = pDest + (nImgHei - nPatchWid) * nImgWid; float *pii1 = NULL; float *pii2 = NULL; float *pii3 = NULL; float *pii4 = NULL; pii1 = pII + ri21 * nImgWid + rj21; // 積分圖搜索區域最大位置的元素，堪稱右下角的元素 pii2 = pII + ri21 * nImgWid; // 可以看成搜索區域中左下角的像素位置 pii3 = pII + rj21; // 搜索區域右上角像素位置 pii4 = pII; // 搜索區域左上角像素位置 pi = pSrc + ri1 * nImgWid; // 定位要處理的像素的位置的那一行 while (pres < pend) { float *pe = pres + nImgWid - nPatchHei; //限定了寬度的范圍 pres = pres + rj1; //定位了要處理的像素點的歸一化系數矩陣的位置 pi = pi + rj1; //定位了要處理的像素點 while (pres < pe) //遍歷整個行 { (*pres++) = (*pres) + pCosDctcMtx[*pi++] * ( (*pii1++) - (*pii2++) - (*pii3++) + (*pii4++) ); //這個其實是計算整個像素塊 } pres += nPatchHei; pi += nPatchHei; pii1 += rj21; pii2 += rj21; pii3 += rj21; pii4 += rj21; } }