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转换函数前文博客中有，或者大家自己实现，都可以。