上一篇開頭提到了一些可用於磨皮的去噪演算法，下面我們實現這些演算法並且觀察效果，咱不考慮實時性的問題

該演算法利用影象區域性統計特性進行濾波處理，例如NXM畫素的灰度圖，首先計算點(i,j)所在視窗內（大小為(2n+1)(2m+1))的平均值m(i,j)

以及均方差：

得到加性去噪後的結果為：

其中：

1.根據原文提出的優化方法，首先是建立兩個積分圖，如圖所示，點4的積分即為Sum(Ra)+Sum(Rb)+Sum(Rc)+Sum(Rd)。積分圖的建立演算法可以參考這篇文章進行簡單優化，然後即可根據積分圖計算公式中的m值和v值。

例如半徑為r的視窗的m(i,j)為Integral(i+r,j+r) + Integral(i-r-1,j-r-1)-Integral(i+r,j-r-1)-Integral(i-r-1,j+r)。程式碼如下，分別求1次方和平方的積分圖。

1. void MagicBeauty::initIntegral(uint8_t* inputMatrix){  
2.     LOGE("initIntegral start");  
3.     if(mIntegralMatrix == NULL)  
4.         mIntegralMatrix = new uint64_t[mImageWidth * mImageHeight];  
5.     if(mIntegralMatrixSqr == NULL)  
6.         mIntegralMatrixSqr = new uint64_t[mImageWidth * mImageHeight];  
7.     uint64_t *columnSum = new uint64_t[mImageWidth];  
8.     uint64_t *columnSumSqr = new uint64_t[mImageWidth];  
9.     columnSum[0] = inputMatrix[0];  
10.     columnSumSqr[0] = inputMatrix[0] * inputMatrix[0];  
11.     mIntegralMatrix[0] = columnSum[0];  
12.     mIntegralMatrixSqr[0] = columnSumSqr[0];  
13.     for(int i = 1;i < mImageWidth;i++){  
14.         columnSum[i] = inputMatrix[i];  
15.         columnSumSqr[i] = inputMatrix[i] * inputMatrix[i];  
16.         mIntegralMatrix[i] = columnSum[i];  
17.         mIntegralMatrix[i] += mIntegralMatrix[i-1];  
18.         mIntegralMatrixSqr[i] = columnSumSqr[i];  
19.         mIntegralMatrixSqr[i] += mIntegralMatrixSqr[i-1];  
20.     }  
21.     for (int i = 1;i < mImageHeight; i++){  
22.         int offset = i * mImageWidth;  
23.         columnSum[0] += inputMatrix[offset];  
24.         columnSumSqr[0] += inputMatrix[offset] * inputMatrix[offset];  
25.         mIntegralMatrix[offset] = columnSum[0];  
26.         mIntegralMatrixSqr[offset] = columnSumSqr[0];  
27.          // other columns
28.         for(int j = 1; j < mImageWidth; j++){  
29.             columnSum[j] += inputMatrix[offset+j];  
30.             columnSumSqr[j] += inputMatrix[offset+j] * inputMatrix[offset+j];  
31.             mIntegralMatrix[offset+j] = mIntegralMatrix[offset+j-1] + columnSum[j];  
32.             mIntegralMatrixSqr[offset+j] = mIntegralMatrixSqr[offset+j-1] + columnSumSqr[j];  
33.         }  
34.     }  
35.     delete[] columnSum;  
36.     delete[] columnSumSqr;  
37.     LOGE("initIntegral end");  
38. }  

2.根據網上抄來的RGB膚色檢測計算膚色區域

1. void MagicBeauty::initSkinMatrix(){  
2.     LOGE("start - initSkinMatrix");  
3.     if(mSkinMatrix == NULL)  
4.         mSkinMatrix = new uint8_t[mImageWidth * mImageHeight];  
5.     for(int i = 0; i < mImageHeight; i++){  
6.         for(int j = 0; j < mImageWidth; j++){  
7.             int offset = i*mImageWidth+j;  
8.             ARGB RGB;  
9.             BitmapOperation::convertIntToArgb(mImageData_rgb[offset],&RGB);  
10.             if ((RGB.blue>95 && RGB.green>40 && RGB.red>20 &&  
11.                     RGB.blue-RGB.red>15 && RGB.blue-RGB.green>15)||//uniform illumination
12.                     (RGB.blue>200 && RGB.green>210 && RGB.red>170 &&  
13.                     abs(RGB.blue-RGB.red)<=15 && RGB.blue>RGB.red&& RGB.green>RGB.red))//lateral illumination
14.                 mSkinMatrix[offset] = 255;  
15.             else
16.                 mSkinMatrix[offset] = 0;  
17.         }  
18.     }  
19.     LOGE("end - initSkinMatrix");  
20. }  

3.根據公式對RGB通道或者將RGB通道轉化為YCbCr格式單獨對Y通道進行濾波

1. void MagicBeauty::startLocalStatisticsSmooth(float sigema){  
2.     if(mIntegralMatrix == NULL || mIntegralMatrixSqr == NULL ||  
3.             mImageData_yuv_y == NULL || mSkinMatrix == NULL || mImageData_yuv == NULL){  
4.         LOGE("not init correctly");  
5.         return;  
6.     }  
7.     int radius = mImageWidth > mImageHeight ? mImageWidth * 0.02 : mImageHeight * 0.02;  
8.     LOGE("startSmooth");  
9.     for(int i = 1; i < mImageHeight; i++){  
10.         for(int j = 1; j < mImageWidth; j++){  
11.             int offset = i * mImageWidth + j;  
12.             if(mSkinMatrix[offset] == 255){  
13.                 int iMax = i + radius >= mImageHeight-1 ? mImageHeight-1 : i + radius;  
14.                 int jMax = j + radius >= mImageWidth-1 ? mImageWidth-1 :j + radius;  
15.                 int iMin = i - radius <= 1 ? 1 : i - radius;  
16.                 int jMin = j - radius <= 1 ? 1 : j - radius;  
17.                 int squar = (iMax - iMin + 1)*(jMax - jMin + 1);  
18.                 int i4 = iMax*mImageWidth+jMax;  
19.                 int i3 = (iMin-1)*mImageWidth+(jMin-1);  
20.                 int i2 = iMax*mImageWidth+(jMin-1);  
21.                 int i1 = (iMin-1)*mImageWidth+jMax;  
22.                 float m = (mIntegralMatrix[i4]  
23.                         + mIntegralMatrix[i3]  
24.                         - mIntegralMatrix[i2]  
25.                         - mIntegralMatrix[i1]) / squar;  
26.                 float v = (mIntegralMatrixSqr[i4]  
27.                         + mIntegralMatrixSqr[i3]  
28.                         - mIntegralMatrixSqr[i2]  
29.                         - mIntegralMatrixSqr[i1]) / squar - m*m;  
30.                 float k = v / (v + sigema);  
31.                     mImageData_yuv[offset*3] = m - k * m + k * mImageData_yuv_y[offset];</span>  
32.             }  
33.         }  
34.     }  
35.     endLocalStatisticsSmooth();  
36. }