MeanShift濾波演算法與實現

本文將嘗試使用MeanShift濾波來做磨皮演算法；

MeanShift即均值漂移，最早由Fukunage在1975年提出，論文名字為：The Estimation of the Gradient of a density function.

MeanShift一般是指一個迭代的步驟，即先算出當前點的偏移均值，然後以此為新的起始點，繼續移動，直到滿足一定的結束條件；MeanShift廣泛應用於圖像聚類、平滑、分割和跟蹤方面，本文主要講的是圖像的平滑濾波，嘗試應用於人像的磨皮演算法中；

我們使用一張圖來講解MeanShift的演算法原理(此圖來自網路)：

Fig.1基本MeanShift演算法示意圖

我們假設起始位置的濾波半徑為Radius，也就是圖a中的藍色圓形區域半徑，圖a為起始位置，假設紅色點為目標像素，每個目標像素包含位置特徵和像素RGB特徵；

1，計算圖a起始位置處，半徑Radius內目標像素的位置特徵和像素RGB特徵的均值M，如圖c所示；

2，將起始位置的初始特徵(位置特徵和RGB特徵)更新為特徵M；

3，計算M處半徑Radius區域內，目標像素的均值特徵M；

4，按照1-3的過程進行迭代，直到滿足一定的迭代次數和限制條件；

5，圖a中起始位置的RGB特徵值即為迭代完成時M的RGB特徵值，如圖f所示；

整個過程也叫均值漂移，實際上不是位置從圖a起始值漂移到了f圖中的位置，而是圖a和圖f處的特徵值歸為了一類，當然這裡指的是RGB像素值；

這裡我們只講最基本的MeanShift平滑濾波演算法，對於改進的MeanShift演算法不做講解；

演算法流程如下：

1，假設當前像素點P(i,j)，濾波半徑為R，迭代次數閾值為maxIter，像素差值閾值為threshold；

2，計算以P為中心，R為半徑的圓形區域S內目標像素的均值特徵，包含像素rgb的均值特徵和位置的均值特徵(質心)，計算公式如下：

其中K為核函數，這裡取得是|x-y|；

1，將P的特徵值M更新為2中計算的新特徵值；

2，按照2-3的步驟進行迭代，直到滿足迭代次數閾值maxIter停止，P處的像素值即迭代終結時的rgb特徵值；

上述即為MeanShift平滑濾波演算法的流程，該演算法最大缺點為速度慢，本文用它來嘗試磨皮效果，採用YCbCr顏色空間，僅對Y通道處理，以此加速；

效果圖如下所示：

完整C代碼如下：

[cpp] view plain copy

1. #include "string.h"
2. #include "stdio.h"
3. #include "stdlib.h"
4. #include "math.h"
5. #include"f_MeanShiftFilter.h"
6. #include"TRGB2YCbCr.h"
7. #define MIN2(a, b) ((a) < (b) ? (a) : (b))
8. #define MAX2(a, b) ((a) > (b) ? (a) : (b))
9. #define CLIP3(x, a, b) MIN2(MAX2(a,x), b)

12. int

    MeanShiftOneChannel(unsigned

    char

    * srcData,

    int

    width ,

    int

    height,

    int

    radius,

    int

    threshold,

    int

    maxIter)
13. {

14. int

    len =

    sizeof

    (unsigned

    long

    ) * width * height;

15. int

    i, j;

16. int

    gray = 0, sum = 0, srcGray = 0, count = 0;
17. unsigned

    char

    * tempData = (unsigned

    char

    *) malloc(

    sizeof

    (unsigned

    char

    ) * height * width);
18. memcpy(tempData, srcData,

    sizeof

    (unsigned

    char

    ) * height * width);

19. for

    (j = 0; j < height; j++ )

20. {

21. for

    (i = 0; i < width; i++)
22. {
23. len = i + j * width;

24. int

    nIter = 0, cx = 0, cy = 0, sumx = 0, sumy = 0;
25. srcGray = tempData[len];
26. cx = i;
27. cy = j;

29. while

    (nIter < maxIter)
30. {
31. sum = 0;
32. sumx = 0;
33. sumy = 0;
34. count = 0;

35. for

    (

    int

    y = cy - radius; y <= cy + radius; y++)
36. {

37. for

    (

    int

    x = cx - radius; x <= cx + radius; x++)
38. {

39. int

    px = CLIP3(x, 0, width - 1);

40. int

    py = CLIP3(y, 0, height - 1);
41. len = px + py * width;
42. gray = tempData[len];

43. if

    (abs(gray - srcGray) < threshold)

44. {
45. count++;
46. sum += gray;
47. sumx += x;
48. sumy += y;
49. }
50. }
51. }

52. if

    (count == 0)

53. break

    ;
54. srcGray = sum / count;
55. cx = sumx / count;
56. cy = sumy / count;
57. nIter++;
58. }
59. srcData[i + j * width] = srcGray;
60. }
61. }
62. free(tempData);

63. return

    0;
64. };

66. void

    f_MeanShiftFilter(unsigned

    char

    * srcData,

    int

    nWidth,

    int

    nHeight,

    int

    nStride,

    int

    radius,

    int

    threshold,

    int

    maxIter)
67. {

68. if

    (srcData == NULL)
69. {

70. return

    ;
71. }

72. if

    (radius == 0 || threshold == 0)

73. return

    ;
74. unsigned

    char

    * yData = (unsigned

    char

    *)malloc(

    sizeof

    (unsigned

    char

    ) * nWidth * nHeight);
75. unsigned

    char

    * cbData = (unsigned

    char

    *)malloc(

    sizeof

    (unsigned

    char

    ) * nWidth * nHeight);
76. unsigned

    char

    * crData = (unsigned

    char

    *)malloc(

    sizeof

    (unsigned

    char

    ) * nWidth * nHeight);
77. unsigned

    char

    * pSrc = srcData;

78. int

    Y, CB, CR;
79. unsigned

    char

    * pY = yData;

80. unsigned

    char

    * pCb = cbData;
81. unsigned

    char

    * pCr = crData;

82. for

    (

    int

    j = 0; j < nHeight; j++)
83. {

84. for

    (

    int

    i = 0; i < nWidth; i++)
85. {
87. RGBToYCbCr(pSrc[2],pSrc[1],pSrc[0],&Y,&CB,&CR);
88. *pY = Y;
89. *pCb = CB;
90. *pCr = CR;
91. pY++;
92. pCb++;
93. pCr++;
94. pSrc += 4;
95. }
96. }
97. MeanShiftOneChannel(yData, nWidth, nHeight, radius, threshold, maxIter);
98. pSrc = srcData;
99. pY = yData;
100. pCb = cbData;
101. pCr = crData;

102. int

     R, G, B;

103. for

     (

     int

     j = 0; j < nHeight; j++)
104. {

105. for

     (

     int

     i = 0; i < nWidth; i++)
106. {
107. YCbCrToRGB(*pY, *pCb, *pCr, &R, &G, &B);
108. pSrc[0] = B;
109. pSrc[1] = G;
110. pSrc[2] = R;
111. pY++;
112. pCb++;
113. pCr++;
114. pSrc += 4;
115. }
116. }
117. free(yData);
118. free(cbData);
119. free(crData);
120. }

代碼已經貼出來了，這裡就不給DEMO了，大家可以直接使用代碼進行測試，代碼中YCbCr轉換函數前文博客中有，或者大家自己實現，都可以。

喜歡這篇文章嗎？立刻分享出去讓更多人知道吧！