【转】七种常见阈值分割代码（Otsu、最⼤熵、迭代法、⾃适应阀值、⼿动、迭代法、基本全局阈值法） ⼀、⼯具：VC+OpenCV
当飞机飞越彩虹⼆、语⾔：C++
三、原理
otsu法（最⼤类间⽅差法，有时也称之为⼤津算法）使⽤的是聚类的思想，把图像的灰度数按灰度级分成2个部分，使得两个部分之间的灰度值差异最⼤，每个部分之间的灰度差异最⼩，通过⽅差的计算来寻⼀个合适的灰度级别来划分。所以可以在⼆值化的时候采⽤otsu 算法来⾃动选取阈值进⾏⼆值化。otsu算法被认为是图像分割中阈值选取的最佳算法，计算简单，不受图像亮度和对⽐度的影响。因此，使类间⽅差最⼤的分割意味着错分概率最⼩。日本人眼里的中国
设t为设定的阈值。
wo：分开后前景像素点数占图像的⽐例
uo：分开后前景像素点的平均灰度
w1：分开后被景像素点数占图像的⽐例
u1：分开后被景像素点的平均灰度
u=w0*u0 + w1*u1 ：图像总平均灰度
从L个灰度级遍历t，使得t为某个值的时候，前景和背景的⽅差最⼤，则这个 t 值便是我们要求得的阈值。
其中，⽅差的计算公式如下：
g=wo * (uo - u) * (uo - u) + w1 * (u1 - u) * (u1 - u)
[ 此公式计算量较⼤，可以采⽤： g = wo * w1 * (uo - u1) * (uo - u1) ]
由于otsu算法是对图像的灰度级进⾏聚类，so 在执⾏otsu算法之前，需要计算该图像的灰度直⽅图。
迭代法原理：迭代选择法是⾸先猜测⼀个初始阈值，然后再通过对图像的多趟计算对阈值进⾏改进的过程。重复地对图像进⾏阈值操作，将图像分割为对象类和背景类，然后来利⽤每⼀个类中的灰阶级别对阈值进⾏改进。
图像阈值分割---迭代算法
1 .处理流程：
1.为全局阈值选择⼀个初始估计值T(图像的平均灰度)。
2.⽤T分割图像。产⽣两组像素：G1有灰度值⼤于T的像素组成，G2有⼩于等于T像素组成。
3.计算G1和G2像素的平均灰度值m1和m2；
4.计算⼀个新的阈值:T = (m1 + m2) / 2;
5.重复步骤2和4，直到连续迭代中的T值间的差⼩于⼀个预定义参数为⽌。
适合图像直⽅图有明显波⾕
四、程序
主程序(核⼼部分)
[cpp]
1. 阈值分割
2. 1 /*===============================图像分割=====================================*/
3. 2 /*---------------------------------------------------------------------------*/
4. 3 /*⼿动设置阀值*/
5. 4 IplImage* binaryImg = cvCreateImage(cvSize(w, h),IPL_DEPTH_8U, 1);
6. 5 cvThreshold(smoothImgGauss,binaryImg,71,255,CV_THRESH_BINARY);遗传与变异
7. 6 cvNamedWindow("cvThreshold", CV_WINDOW_AUTOSIZE );
8. 7 cvShowImage( "cvThreshold", binaryImg );
9. 8 //cvReleaseImage(&binaryImg);
10. 9 /*---------------------------------------------------------------------------*/
11. 10 /*⾃适应阀值 //计算像域邻域的平均灰度，来决定⼆值化的值*/
12. 11 IplImage* adThresImg = cvCreateImage(cvSize(w, h),IPL_DEPTH_8U, 1);
13. 12 double max_value=255;
14. 13 int adpative_method=CV_ADAPTIVE_THRESH_GAUSSIAN_C;//CV_ADAPTIVE_THRESH_MEAN_C
15. 14 int threshold_type=CV_THRESH_BINARY;
16. 15 int block_size=3;//阈值的象素邻域⼤⼩
17. 16 int offset=5;//窗⼝尺⼨
18. 17 cvAdaptiveThreshold(smoothImgGauss,adThresImg,max_value,adpative_method,threshold_type,block_size,offset);
19. 18 cvNamedWindow("cvAdaptiveThreshold", CV_WINDOW_AUTOSIZE );
20. 19 cvShowImage( "cvAdaptiveThreshold", adThresImg );
21. 20 cvReleaseImage(&adThresImg);
22. 21 /*---------------------------------------------------------------------------*/
23. 22 /*最⼤熵阀值分割法*/
24. 23 IplImage* imgMaxEntropy = cvCreateImage(cvGetSize(imgGrey),IPL_DEPTH_8U,1);
25. 24 MaxEntropy(smoothImgGauss,imgMaxEntropy);
26. 25 cvNamedWindow("MaxEntroyThreshold", CV_WINDOW_AUTOSIZE );
27. 26 cvShowImage( "MaxEntroyThreshold", imgMaxEntropy );//显⽰图像
28. 27 cvReleaseImage(&imgMaxEntropy );
29. 28 /*---------------------------------------------------------------------------*/
30. 29 /*基本全局阀值法*/
31. 30 IplImage* imgBasicGlobalThreshold = cvCreateImage(cvGetSize(imgGrey),IPL_DEPTH_8U,1);
32. 31 cvCopyImage(srcImgGrey,imgBasicGlobalThreshold);
33. 32 int pg[256],i,thre;
34. 33 for (i=0;i<256;i++) pg[i]=0;
电化教育研究35. 34 for (i=0;i<imgBasicGlobalThreshold->imageSize;i++) // 直⽅图统计
36. 35 pg[(BYTE)imgBasicGlobalThreshold->imageData[i]]++;
37. 36 thre = BasicGlobalThreshold(pg,0,256); // 确定阈值
38. 37 cout<<"The Threshold of this Image in BasicGlobalThreshold is:"<<thre<<endl;//输出显⽰阀值
39. 38 cvThreshold(imgBasicGlobalThreshold,imgBasicGlobalThreshold,thre,255,CV_THRESH_BINARY); // ⼆值化
40. 39 cvNamedWindow("BasicGlobalThreshold", CV_WINDOW_AUTOSIZE );
41. 40 cvShowImage( "BasicGlobalThreshold", imgBasicGlobalThreshold);//显⽰图像
42. 41 cvReleaseImage(&imgBasicGlobalThreshold);
43. 42 /*---------------------------------------------------------------------------*/
44. 43 /*OTSU*/
45. 44 IplImage* imgOtsu = cvCreateImage(cvGetSize(imgGrey),IPL_DEPTH_8U,1);
46. 45 cvCopyImage(srcImgGrey,imgOtsu);
47. 46 int thre2;
48. 47 thre2 = otsu2(imgOtsu);
49. 48 cout<<"The Threshold of this Image in Otsu is:"<<thre2<<endl;//输出显⽰阀值
50. 49 cvThreshold(imgOtsu,imgOtsu,thre2,255,CV_THRESH_BINARY); // ⼆值化
51. 50 cvNamedWindow("imgOtsu", CV_WINDOW_AUTOSIZE );
52. 51 cvShowImage( "imgOtsu", imgOtsu);//显⽰图像
53. 52 cvReleaseImage(&imgOtsu);
54. 53 /*---------------------------------------------------------------------------*/
55. 54 /*上下阀值法：利⽤正态分布求可信区间*/
56. 55 IplImage* imgTopDown = cvCreateImage( cvGetSize(imgGrey), IPL_DEPTH_8U, 1 );
57. 56 cvCopyImage(srcImgGrey,imgTopDown);
58. 57 CvScalar mean ,std_dev;//平均值、标准差
59. 58 double u_threshold,d_threshold;
60. 59 cvAvgSdv(imgTopDown,&mean,&std_dev,NULL);
61. 60 u_threshold = mean.val[0] +2.5* std_dev.val[0];//上阀值
62. 61 d_threshold = mean.val[0] -2.5* std_dev.val[0];//下阀值
63. 62 //u_threshold = mean + 2.5 * std_dev; //错误
64. 63 //d_threshold = mean - 2.5 * std_dev;
65. 64 cout<<"The TopThreshold of this Image in TopDown is:"<<d_threshold<<endl;//输出显⽰阀值
66. 65 cout<<"The DownThreshold of this Image in TopDown is:"<<u_threshold<<endl;
67. 66 cvThreshold(imgTopDown,imgTopDown,d_threshold,u_threshold,CV_THRESH_BINARY_INV);//上下阀值
68. 67 cvNamedWindow("imgTopDown", CV_WINDOW_AUTOSIZE );
69. 68 cvShowImage( "imgTopDown", imgTopDown);//显⽰图像
70. 69 cvReleaseImage(&imgTopDown);
71. 70 /*---------------------------------------------------------------------------*/
72. 71 /*迭代法*/
73. 72 IplImage* imgIteration = cvCreateImage( cvGetSize(imgGrey), IPL_DEPTH_8U, 1 );
74. 73 cvCopyImage(srcImgGrey,imgIteration);
75. 74 int thre3,nDiffRec;
76. 75 thre3 =DetectThreshold(imgIteration, 100, nDiffRec);
77. 76 cout<<"The Threshold of this Image in imgIteration is:"<<thre3<<endl;//输出显⽰阀值
78. 77 cvThreshold(imgIteration,imgIteration,thre3,255,CV_THRESH_BINARY_INV);//上下阀值
79. 78 cvNamedWindow("imgIteration", CV_WINDOW_AUTOSIZE );
80. 79 cvShowImage( "imgIteration", imgIteration);
81. 80 cvReleaseImage(&imgIteration);
82. 迭代
83. 1 /*======================================================================*/
84. 2 /* 迭代法*/
85. 3 /*======================================================================*/
86. 4 // nMaxIter：最⼤迭代次数；nDiffRec：使⽤给定阀值确定的亮区与暗区平均灰度差异值
87. 5 int DetectThreshold(IplImage*img, int nMaxIter, int& iDiffRec) //阀值分割：迭代法
88. 6 {
89. 7 //图像信息
90. 8 int height = img->height;
91. 9 int width = img->width;
92. 10 int step = img->widthStep/sizeof(uchar);
93. 11 uchar *data = (uchar*)img->imageData;
94. 12
95. 13 iDiffRec =0;
96. 14 int F[256]={ 0 }; //直⽅图数组
97. 15 int iTotalGray=0;//灰度值和
98. 16 int iTotalPixel =0;//像素数和
99. 17 byte bt;//某点的像素值
100. 18
101. 19 uchar iThrehold,iNewThrehold;//阀值、新阀值
102. 20 uchar iMaxGrayValue=0,iMinGrayValue=255;//原图像中的最⼤灰度值和最⼩灰度值
103. 21 uchar iMeanGrayValue1,iMeanGrayValue2;
104. 22
105. 23 //获取(i,j)的值，存于直⽅图数组F
106. 24 for(int i=0;i<width;i++)
107. 25 {
108. 26 for(int j=0;j<height;j++)
109. 27 {
110. 28 bt = data[i*step+j];
111. 29 if(bt<iMinGrayValue)
112. 30 iMinGrayValue = bt;
113. 31 if(bt>iMaxGrayValue)
114. 32 iMaxGrayValue = bt;
115. 33 F[bt]++;
116. 34 }
117. 35 }
118. 36
119. 37 iThrehold =0;//
120. 38 iNewThrehold = (iMinGrayValue+iMaxGrayValue)/2;//初始阀值
121. 39 iDiffRec = iMaxGrayValue - iMinGrayValue;
122. 40
123. 41 for(int a=0;(abs(iThrehold-iNewThrehold)>0.5)&&a<nMaxIter;a++)//迭代中⽌条件
124. 42 {
125. 43 iThrehold = iNewThrehold;
126. 44 //⼩于当前阀值部分的平均灰度值
127. 45 for(int i=iMinGrayValue;i<iThrehold;i++)
128. 46 {
129. 47 iTotalGray += F[i]*i;//F[]存储图像信息
130. 48 iTotalPixel += F[i];
131. 49 }
132. 50 iMeanGrayValue1 = (uchar)(iTotalGray/iTotalPixel);
133. 51 //⼤于当前阀值部分的平均灰度值
134. 52 iTotalPixel =0;
135. 53 iTotalGray =0;
136. 54 for(int j=iThrehold+1;j<iMaxGrayValue;j++)
138. 56 iTotalGray += F[j]*j;//F[]存储图像信息
139. 57 iTotalPixel += F[j];
140. 58 }
141. 59 iMeanGrayValue2 = (uchar)(iTotalGray/iTotalPixel);
142. 60
143. 61 iNewThrehold = (iMeanGrayValue2+iMeanGrayValue1)/2; //新阀值
144. 62 iDiffRec = abs(iMeanGrayValue2 - iMeanGrayValue1);
145. 63 }
146. 64
147. 65 //cout<<"The Threshold of this Image in imgIteration is:"<<iThrehold<<endl;
148. 66 return iThrehold;
149. 67 }
150. 68
151. Otsu代码⼀
152. 1 /*======================================================================*/ 153. 2 /* OTSU global thresholding routine */
154. 3 /* takes a 2D unsigned char array pointer, number of rows, and */
155. 4 /* number of cols in the array. returns the value of the threshold */
156. 5 /*parameter:
157. 6 *image --- buffer for image
158. 7 rows, cols --- size of image
159. 8 x0, y0, dx, dy --- region of vector used for computing threshold
国税31号文
160. 9 vvv --- debug option, is 0, no debug information outputed
161. 10 */
162. 11 /*
163. 12 OTSU 算法可以说是⾃适应计算单阈值（⽤来转换灰度图像为⼆值图像）的简单⾼效⽅法。164. 13 下⾯的代码最早由 Ryan Dibble提供，此后经过多⼈Joerg.Schulenburg, R.Z.Liu 等修改，补正。165. 14 算法对输⼊的灰度图像的直⽅图进⾏分析，将直⽅图分成两个部分，使得两部分之间的距离最⼤。166. 15 划分点就是求得的阈值。
167. 16 */
168. 17 /*======================================================================*/ 169. 18 int otsu (unsigned char*image, int rows, int cols, int x0, int y0, int dx, int dy, int vvv)
170. 19 {
171. 20
172. 21 unsigned char*np; // 图像指针
173. 22 int thresholdValue=1; // 阈值
174. 23 int ihist[256]; // 图像直⽅图，256个点
175. 24
176. 25 int i, j, k; // various counters
节油剂
177. 26 int n, n1, n2, gmin, gmax;
178. 27 double m1, m2, sum, csum, fmax, sb;
179. 28
180. 29 // 对直⽅图置零
181. 30 memset(ihist, 0, sizeof(ihist));
182. 31
183. 32 gmin=255; gmax=0;
184. 33 // ⽣成直⽅图
185. 34 for (i = y0 +1; i < y0 + dy -1; i++)
186. 35 {
187. 36 np = (unsigned char*)image[i*cols+x0+1];
188. 37 for (j = x0 +1; j < x0 + dx -1; j++)
189. 38 {
190. 39 ihist[*np]++;
191. 40 if(*np > gmax) gmax=*np;
192. 41 if(*np < gmin) gmin=*np;
193. 42 np++; /* next pixel */
194. 43 }
195. 44 }
196. 45
197. 46 // set up everything
198. 47 sum = csum =0.0;
199. 48 n =0;
200. 49
201. 50 for (k =0; k <=255; k++)
202. 51 {
203. 52 sum += (double) k * (double) ihist[k]; /* x*f(x) 质量矩*/
204. 53 n += ihist[k]; /* f(x) 质量 */
206. 55
207. 56 if (!n)
208. 57 {
209. 58 // if n has no value, there
210. 59 fprintf (stderr, "NOT NORMAL thresholdValue = 160\n");
211. 60 return (160);
212. 61 }
213. 62
214. 63 // do the otsu global thresholding method
215. 64 fmax =-1.0;
216. 65 n1 =0;
217. 66 for (k =0; k <255; k++)
218. 67 {
219. 68 n1 += ihist[k];
220. 69 if (!n1)
221. 70 {
222. 71 continue;
223. 72 }
224. 73 n2 = n - n1;
225. 74 if (n2 ==0)
226. 75 {
227. 76 break;
228. 77 }
229. 78 csum += (double) k *ihist[k];
230. 79 m1 = csum / n1;
231. 80 m2 = (sum - csum) / n2;
232. 81 sb = (double) n1 *(double) n2 *(m1 - m2) * (m1 - m2);
233. 82 /* bbg: note: can be optimized. */
234. 83 if (sb > fmax)
235. 84 {
236. 85 fmax = sb;
237. 86 thresholdValue = k;
238. 87 }
239. 88 }
240. 89
241. 90 // at this point we have our thresholding value
242. 91
243. 92 // debug code to display thresholding values
244. 93 if ( vvv &1 )
245. 94 fprintf(stderr,"# OTSU: thresholdValue = %d gmin=%d gmax=%d\n",
246. 95 thresholdValue, gmin, gmax);
247. 96
248. 97 return(thresholdValue);
249. 98 }
250. Otsu代码⼆
251. 1 /*======================================================================*/ 252. 2 /* OTSU global thresholding routine */
253. 3 /*======================================================================*/ 254. 4 int otsu2 (IplImage *image)
255. 5 {
256. 6 int w = image->width;
257. 7 int h = image->height;
258. 8
259. 9 unsigned char*np; // 图像指针
260. 10 unsigned char pixel;
261. 11 int thresholdValue=1; // 阈值
262. 12 int ihist[256]; // 图像直⽅图，256个点
263. 13
264. 14 int i, j, k; // various counters
265. 15 int n, n1, n2, gmin, gmax;
266. 16 double m1, m2, sum, csum, fmax, sb;
267. 17
268. 18 // 对直⽅图置零...
269. 19 memset(ihist, 0, sizeof(ihist));
270. 20
271. 21 gmin=255; gmax=0;
272. 22 // ⽣成直⽅图
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议