学习OpenCV2——MeanShift之图形分割

1. 原理

    用meanshift做图像平滑和分割，其实是一回事。其本质是经过迭代，将收敛点的像素值代替原来的像素值，从而去除了局部相似的纹理，同时保留了边缘等差异较大的特征。

        OpenCV中自带有基于meanshift的分割方法pyrMeanShiftFiltering()。由函数名pyrMeanShiftFiltering可知，这里是将meanshift算法和图像金字塔相结合用来分割的。

 

1. <span style="font-size:18px;">void PyrMeanShiftFiltering( const CvArr* srcarr, //输入图像

2. CvArr* dstarr, //输出图像

3. double sp, //颜色域半径

4. double sr, //空间域半径

5. int max_level, //金字塔最大层数

6. CvTermCriteria termcrit ) //迭代终止条件</span>

 

    要求输入和输出图像都是CV_8UC3类型，而且两者尺寸一样。实际上并不需要去先定义dstarr，因为程序里会将srcarr的格式赋值给dstarr。

    termcrit有三种情况，迭代次数、迭代精度和两者同时满足。默认为迭代次数为5同时迭代精度为1。termcrit是个结构体，其结构如下

 

1. <span style="font-size:18px;">typedef struct CvTermCriteria

2. {

3. int type; /*CV_TERMCRIT_ITER或CV_TERMCRIT_EPS 或二者都是*/

4. int max_iter; /* 最大迭代次数 */

5. double epsilon; /* 结果的精确性 */

6. }

7. CvTermCriteria;</span>

     使用pyrMeanShiftFiltering()进行图像分割非常简单，只需要定义sp0,sr,max_level和termrit，然后调用pyrMeanShiftFiltering()就行了。

 

    在实际操作时，为了使分割的结果显示得更明显，经常用floodFill( )将不同连通域涂上不同的颜色。具体情况参看下 面的实例。

 

2. 程序实例

    来看看OpenCV自带的一个用meanshift进行分割的例子

    原程序见   “  .\OpenCV249\sources\samples\cpp\meanshift_segmentation.cpp”

 

1. <span style="font-size:18px;">#include "opencv2/highgui/highgui.hpp"

2. #include "opencv2/core/core.hpp"

3. #include "opencv2/imgproc/imgproc.hpp"

4.  

5. #include <iostream>

6.  

7. using namespace cv;

8. using namespace std;

9.  

10. static void help(char** argv)

11. {

12. cout << "\nDemonstrate mean-shift based color segmentation in spatial pyramid.\n"

13. << "Call:\n " << argv[0] << " image\n"

14. << "This program allows you to set the spatial and color radius\n"

15. << "of the mean shift window as well as the number of pyramid reduction levels explored\n"

16. << endl;

17. }

18.  

19. //This colors the segmentations

20. static void floodFillPostprocess( Mat& img, const Scalar& colorDiff=Scalar::all(1) )

21. {

22. CV_Assert( !img.empty() );

23. RNG rng = theRNG();

24. Mat mask( img.rows+2, img.cols+2, CV_8UC1, Scalar::all(0) );

25. for( int y = 0; y < img.rows; y++ )

26. {

27. for( int x = 0; x < img.cols; x++ )

28. {

29. if( mask.at<uchar>(y+1, x+1) == 0 )

30. {

31. Scalar newVal( rng(256), rng(256), rng(256) );

32. floodFill( img, mask, Point(x,y), newVal, 0, colorDiff, colorDiff );

33. }

34. }

35. }

36. }

37.  

38. string winName = "meanshift";

39. int spatialRad, colorRad, maxPyrLevel;

40. Mat img, res;

41.  

42. static void meanShiftSegmentation( int, void* )

43. {

44. cout << "spatialRad=" << spatialRad << "; "

45. << "colorRad=" << colorRad << "; "

46. << "maxPyrLevel=" << maxPyrLevel << endl;

47. pyrMeanShiftFiltering( img, res, spatialRad, colorRad, maxPyrLevel );

48. //Mat imgGray;

49. //cvtColor(res,imgGray,CV_RGB2GRAY);

50. //imshow("res",res);

51. floodFillPostprocess( res, Scalar::all(2) );

52. imshow( winName, res );

53. }

54.  

55. int main(int argc, char** argv)

56. {

57. img = imread("rubberwhale1.png");

58. //img = imread("pic2.png");

59.  
60.  

61. if( img.empty() )

62. return -1;

63.  

64. spatialRad = 10;

65. colorRad = 10;

66. maxPyrLevel = 1;

67.  

68. namedWindow( winName, WINDOW_AUTOSIZE );

69. //imshow("img",img);

70.  
71.  

72. createTrackbar( "spatialRad", winName, &spatialRad, 80, meanShiftSegmentation );

73. createTrackbar( "colorRad", winName, &colorRad, 60, meanShiftSegmentation );

74. createTrackbar( "maxPyrLevel", winName, &maxPyrLevel, 5, meanShiftSegmentation );

75.  

76. meanShiftSegmentation(0, 0);

77. //floodFillPostprocess( img, Scalar::all(2) );

78. //imshow("img2",img);

79. waitKey();

80. return 0;

81. }</span>

 

 

程序很简单，来看看floodFill()函数，有两种形式
    int floodFill( InputOutputArray image, Point seedPoint, Scalar newVal, CV_OUT Rect* rect=0, Scalar loDiff=Scalar(), Scalar upDiff=Scalar(), int flags=4 );
    int floodFill( InputOutputArray image,  InputOutputArray mask, Point seedPoint,  Scalar newVal,  CV_OUT Rect* rect=0,  Scalar loDiff=Scalar(),  Scalar upDiff=Scalar(),  int flags=4 );

     InputOutputArray image    输入输出图像，要求格式为1通道或3通道，8位或浮点

     InputOutputArray mask   掩膜，比image的宽和高各大两像素点

     Point seedPoint    填充的起始点

    Scalar newVal   像素点被染色的值

    CV_OUT Rect* rect=0  可选参数，设置floodFill()要重绘区域的最小边界矩形区域

    Scalar loDiff=Scalar()  定义当前像素值与起始点像素值的亮度或颜色负差的最大值

    Scalar upDiff=Scalar()  定义当前像素值与起始点像素值的亮度或颜色正差的最大值

    flags 操作标志符    

 

程序结果

    处理后一些细小的纹理都平滑掉了，例如图中绿色线条所指示的区域。未填充时，很多地方看得并不明显，填充后就能明显看出差别来了。填充后的图很好地体现了meanshift聚类的思想！

    

    再来看一组更“夸张”的效果图

    使用meanshift方法进行处理后，原来的三个矩形区域消失了！平滑掉了！    

   

    meanshift算法的两个关键参数是空间域半径sr和颜色域半径sp，别说max_level，那是构建图像金字塔的参数好吧。最后，我们来看看sr和sp对结果的影响。

       显然颜色域半径sp对结果的影响比空间域半径sr对结果的影响大。sp和sr越小，细节保留得越多，sp和sr越大，平滑力度越大。边缘和颜色突变的区域的特征保留的较好。因为meanshift要对每个像素点进行操作，所以算法的时间花销很大。

 

3. 深入代码

 

1. <span style="font-size:14px;">/****************************************************************************************\

2. * Meanshift *

3. \****************************************************************************************/

4.  

5. CV_IMPL void

6. cvPyrMeanShiftFiltering( const CvArr* srcarr, CvArr* dstarr,

7. double sp0, double sr, int max_level,

8. CvTermCriteria termcrit )

9. {

10. const int cn = 3;

11. const int MAX_LEVELS = 8;

12.  

13. if( (unsigned)max_level > (unsigned)MAX_LEVELS )

14. CV_Error( CV_StsOutOfRange, "The number of pyramid levels is too large or negative" ); //限定max_level不超过8

15.  

16. std::vector<cv::Mat> src_pyramid(max_level+1); //+1是因为原始图和最终图都定义为图像金字塔的第0层

17. std::vector<cv::Mat> dst_pyramid(max_level+1);

18. cv::Mat mask0;

19. int i, j, level;

20. //uchar* submask = 0;

21.  

22. #define cdiff(ofs0) (tab[c0-dptr[ofs0]+255] + \

23. tab[c1-dptr[(ofs0)+1]+255] + tab[c2-dptr[(ofs0)+2]+255] >= isr22)

24.  

25. double sr2 = sr * sr;

26. int isr2 = cvRound(sr2), isr22 = MAX(isr2,16);

27. int tab[768];

28. cv::Mat src0 = cv::cvarrToMat(srcarr); //arr转Mat

29. cv::Mat dst0 = cv::cvarrToMat(dstarr);

30.  

31. //确保src和dst都是CV_8UC3，且同尺寸

32. if( src0.type() != CV_8UC3 )

33. CV_Error( CV_StsUnsupportedFormat, "Only 8-bit, 3-channel images are supported" );

34. if( src0.type() != dst0.type() )

35. CV_Error( CV_StsUnmatchedFormats, "The input and output images must have the same type" );

36. if( src0.size() != dst0.size() )

37. CV_Error( CV_StsUnmatchedSizes, "The input and output images must have the same size" );

38.  

39. //确保迭代次数在1到100次，默认则为5；迭代精度默认为1.

40. if( !(termcrit.type & CV_TERMCRIT_ITER) )

41. termcrit.max_iter = 5;

42. termcrit.max_iter = MAX(termcrit.max_iter,1);

43. termcrit.max_iter = MIN(termcrit.max_iter,100);

44. if( !(termcrit.type & CV_TERMCRIT_EPS) )

45. termcrit.epsilon = 1.f;

46. termcrit.epsilon = MAX(termcrit.epsilon, 0.f);

47.  

48. for( i = 0; i < 768; i++ )

49. tab[i] = (i - 255)*(i - 255); //tab[]存的是（-255）^2到512^2

50.  

51. // 1. 构造金字塔

52. src_pyramid[0] = src0;

53. dst_pyramid[0] = dst0;

54. for( level = 1; level <= max_level; level++ )

55. {

56. //src_pyramid和dst_pyramid尺寸一样，下一层是上一层尺寸的一半

57. src_pyramid[level].create( (src_pyramid[level-1].rows+1)/2,

58. (src_pyramid[level-1].cols+1)/2, src_pyramid[level-1].type() );

59. dst_pyramid[level].create( src_pyramid[level].rows,

60. src_pyramid[level].cols, src_pyramid[level].type() );

61. //对src_pyramid[level-1]下采样，结果存入src_pyramid[level]

62. cv::pyrDown( src_pyramid[level-1], src_pyramid[level], src_pyramid[level].size() );

63. //CV_CALL( cvResize( src_pyramid[level-1], src_pyramid[level], CV_INTER_AREA ));

64. }

65.  

66. mask0.create(src0.rows, src0.cols, CV_8UC1);

67. //CV_CALL( submask = (uchar*)cvAlloc( (sp+2)*(sp+2) ));

68.  

69. // 2. 从顶层（最小层）开始应用meanshift算法。

70. for( level = max_level; level >= 0; level-- )

71. {

72. cv::Mat src = src_pyramid[level];

73. cv::Size size = src.size();

74. uchar* sptr = src.data; //sptr指向图像矩阵的起始地址，也就是第一行的起始地址

75. int sstep = (int)src.step; //sstep是图像矩阵每一行的长度（以字节为单位），以便后面计算地址

76. uchar* mask = 0;

77. int mstep = 0;

78. uchar* dptr;

79. int dstep;

80. float sp = (float)(sp0 / (1 << level));

81. sp = MAX( sp, 1 ); //这里保证了sp≥1，那么窗口最小是3×3

82.  

83. //这段语句主要作用1、通过上采样得到dst_pyramid[level]；2、得到掩码mask

84. if( level < max_level )

85. {

86. cv::Size size1 = dst_pyramid[level+1].size();

87. cv::Mat m( size.height, size.width, CV_8UC1, mask0.data );

88. dstep = (int)dst_pyramid[level+1].step;

89. dptr = dst_pyramid[level+1].data + dstep + cn;

90. mstep = (int)m.step;

91. mask = m.data + mstep;

92. //cvResize( dst_pyramid[level+1], dst_pyramid[level], CV_INTER_CUBIC );

93. cv::pyrUp( dst_pyramid[level+1], dst_pyramid[level], dst_pyramid[level].size() ); //上采样

94. m.setTo(cv::Scalar::all(0));

95.  

96. for( i = 1; i < size1.height-1; i++, dptr += dstep - (size1.width-2)*3, mask += mstep*2 )

97. {

98. for( j = 1; j < size1.width-1; j++, dptr += cn )

99. {

100. int c0 = dptr[0], c1 = dptr[1], c2 = dptr[2];

101. mask[j*2 - 1] = cdiff(-3) || cdiff(3) || cdiff(-dstep-3) || cdiff(-dstep) ||

102. cdiff(-dstep+3) || cdiff(dstep-3) || cdiff(dstep) || cdiff(dstep+3);

103. }

104. }

105.  

106. cv::dilate( m, m, cv::Mat() ); //对m膨胀

107. mask = m.data;

108. }

109.  

110. dptr = dst_pyramid[level].data; //dptr指向图像矩阵起始地址

111. dstep = (int)dst_pyramid[level].step; //dstep表示图像矩阵每一行的占内存的字节数

112.  

113. for( i = 0; i < size.height; i++, sptr += sstep - size.width*3,

114. dptr += dstep - size.width*3, //每处理完一行，sptr和dptr都指向下一行的起始地址

115. mask += mstep )

116. {

117. for( j = 0; j < size.width; j++, sptr += 3, dptr += 3 ) //每处理完一列,sptr和dptr都指向同行下一列像素的起始地址，所以sptr和dptr实际就是每个像素点的地址

118. {

119. int x0 = j, y0 = i, x1, y1, iter;

120. int c0, c1, c2;

121.  

122. if( mask && !mask[j] )

123. continue;

124.  

125. c0 = sptr[0], c1 = sptr[1], c2 = sptr[2]; //分别对应像素点三通道的地址

126.  

127. // iterate meanshift procedure

128. for( iter = 0; iter < termcrit.max_iter; iter++ )

129. {

130. uchar* ptr;

131. int x, y, count = 0;

132. int minx, miny, maxx, maxy;

133. int s0 = 0, s1 = 0, s2 = 0, sx = 0, sy = 0; //(x,y)的迭代的坐标值，(s0,s1,s2)是迭代的3通道分量值

134. double icount;

135. int stop_flag;

136.  

137. //mean shift: process pixels in window (p-sigmaSp)x(p+sigmaSp)

138. minx = cvRound(x0 - sp); minx = MAX(minx, 0); //若j-sp>=0,则minx=(j-sp)，否则，minx=0;

139. miny = cvRound(y0 - sp); miny = MAX(miny, 0); //若i-sp>=0,则miny=(i-sp)，否则，miny=0;

140. maxx = cvRound(x0 + sp); maxx = MIN(maxx, size.width-1); //若j+sp<=width+1,则maxx=j+sp，否则,maxx=width-1;

141. maxy = cvRound(y0 + sp); maxy = MIN(maxy, size.height-1); //若i+sp<=height+1,则maxy=i+sp，否则,maxy=height-1;

142. ptr = sptr + (miny - i)*sstep + (minx - j)*3; //sptr指向(i,j)，ptr则指向当前窗口第一个像素点

143.  

144. for( y = miny; y <= maxy; y++, ptr += sstep - (maxx-minx+1)*3 ) //窗口内，每处理完一行，ptr指向下一行首地址

145. {

146. int row_count = 0;

147. x = minx;

148. #if CV_ENABLE_UNROLLED

149. for( ; x + 3 <= maxx; x += 4, ptr += 12 ) //这两次for循环是什么意思？颜色限定和空间限定？

150. {

151. int t0 = ptr[0], t1 = ptr[1], t2 = ptr[2];

152. if( tab[t0-c0+255] + tab[t1-c1+255] + tab[t2-c2+255] <= isr2 )

153. {

154. s0 += t0; s1 += t1; s2 += t2;

155. sx += x; row_count++;

156. }

157. t0 = ptr[3], t1 = ptr[4], t2 = ptr[5];

158. if( tab[t0-c0+255] + tab[t1-c1+255] + tab[t2-c2+255] <= isr2 )

159. {

160. s0 += t0; s1 += t1; s2 += t2;

161. sx += x+1; row_count++;

162. }

163. t0 = ptr[6], t1 = ptr[7], t2 = ptr[8];

164. if( tab[t0-c0+255] + tab[t1-c1+255] + tab[t2-c2+255] <= isr2 )

165. {

166. s0 += t0; s1 += t1; s2 += t2;

167. sx += x+2; row_count++;

168. }

169. t0 = ptr[9], t1 = ptr[10], t2 = ptr[11];

170. if( tab[t0-c0+255] + tab[t1-c1+255] + tab[t2-c2+255] <= isr2 )

171. {

172. s0 += t0; s1 += t1; s2 += t2;

173. sx += x+3; row_count++;

174. }

175. }

176. #endif

177. for( ; x <= maxx; x++, ptr += 3 )

178. {

179. int t0 = ptr[0], t1 = ptr[1], t2 = ptr[2];

180. if( tab[t0-c0+255] + tab[t1-c1+255] + tab[t2-c2+255] <= isr2 )

181. {

182. s0 += t0; s1 += t1; s2 += t2;

183. sx += x; row_count++;

184. }

185. }

186. count += row_count;

187. sy += y*row_count;

188. }

189.  

190. if( count == 0 )

191. break;

192.  

193. icount = 1./count;

194. x1 = cvRound(sx*icount);

195. y1 = cvRound(sy*icount);

196. s0 = cvRound(s0*icount);

197. s1 = cvRound(s1*icount);

198. s2 = cvRound(s2*icount);

199.  

200. stop_flag = (x0 == x1 && y0 == y1) || abs(x1-x0) + abs(y1-y0) +

201. tab[s0 - c0 + 255] + tab[s1 - c1 + 255] +

202. tab[s2 - c2 + 255] <= termcrit.epsilon;

203.  

204. x0 = x1; y0 = y1;

205. c0 = s0; c1 = s1; c2 = s2;

206.  

207. if( stop_flag )

208. break;

209. }

210.  

211. dptr[0] = (uchar)c0;

212. dptr[1] = (uchar)c1;

213. dptr[2] = (uchar)c2;

214. }

215. }

216. }

217. }

218.  

219. void cv::pyrMeanShiftFiltering( InputArray _src, OutputArray _dst,

220. double sp, double sr, int maxLevel,

221. TermCriteria termcrit )

222. {

223. Mat src = _src.getMat();

224.  

225. if( src.empty() )

226. return;

227.  

228. _dst.create( src.size(), src.type() );

229. CvMat c_src = src, c_dst = _dst.getMat();

230. cvPyrMeanShiftFiltering( &c_src, &c_dst, sp, sr, maxLevel, termcrit );

231. }</span><span style="font-size:18px;">

232.