基于OTSU算法和基本粒子群优化算法的双阈值图像分割

 OTSU自适应阈值求法与粒子群算法的合作，将OTSU算法作为粒子群算法的适应值函数，来计算每个粒子的适应度与最优阈值相比较，经过3000次迭代最后取得优化后的阈值

原图：

经过联合算法优化的双阈值为90 ，140

将背景像素置0：

效果图：

利用所取得的阈值就可以将图像背景和目标区分开来，利用所得阈值二值化后

效果图：

 

通过效果图可知将人这个目标从背景中分割出来了

源代码：

[cpp]  view plain copy

1. #include "stdafx.h"  
2. #include "cv.h"  
3. #include "highgui.h"  
4. #include "cxcore.h"  
5. #include "time.h"  
6. using namespace std;  
7.   
8. #define rnd( low,uper) ((int)(((double)rand()/(double)RAND_MAX)*((double)(uper)-(double)(low))+(double)(low)+0.5))  
9. /*************************************************************8888 
10. 粒子群算法变量的说明 
11. ******************************************************************************/  
12. const int number = 20;  
13. int antThreshold[number][2];//以阈值作为粒子  
14. int vect[number][2];//更新的速度  
15. float pbest[number] = {0.0};;//每个粒子历史最优解  
16. float gbest = 0.0;//全局历史最优解  
17. int pbestThreshold[number][2];//每个粒子的最优历史阈值  
18. int gbestThreshold[2];//全局粒子的最优阈值  
19.   
20. float w = 0.9;//惯性因子  
21. float c1 = 2.0;//加速因子1  
22. float c2 = 2.0;//加速因子2  
23.   
24. //histogram     
25. float histogram[256]={0};    
26. /*********************************************************************8888 
27. 函数名：GetAvgValue 
28. 参数类型：IplImage* src 
29. 实现功能：获得灰度图像的总平均灰度值 
30. *****************************************************************************/  
31. float GetAvgValue(IplImage* src)    
32. {    
33.     int height=src->height;    
34.     int width=src->width;        
35.     
36.   
37.     for(int i=0;i<height;i++)   
38.     {    
39.         unsigned char* p=(unsigned char*)src->imageData+src->widthStep*i;    
40.         for(int j=0;j<width;j++)   
41.         {    
42.             histogram[*p++]++;    
43.         }    
44.     }    
45.     //normalize histogram     
46.     int size=height*width;    
47.     for(int i=0;i<256;i++) {    
48.         histogram[i]=histogram[i]/size;    
49.     }    
50.     
51.     //average pixel value     
52.     float avgValue=0;    
53.     for(int i=0;i<256;i++) {    
54.         avgValue+=i*histogram[i];    
55.     }    
56.     return avgValue;  
57. }  
58. /***************************************************************************** 
59. 函数名：ThresholdOTSU 
60. 参数类型：int threshold1 , int threshold2 , float avgValue 
61. 功能：求得最大类间方差 
62. **********************************************************************************/  
63. float  ThresholdOTSU(int threshold1 , int threshold2 , float avgValue)  
64. {  
65.     
66.     int threshold;      
67.     float maxVariance=0;    
68.     float w=0,u=0;    
69.     for(int i=threshold1;i< threshold2 ;i++)  
70.     {    
71.         w+=histogram[i];    
72.         u+=i*histogram[i];    
73.     }  
74.     
75.         float t=avgValue*w-u;    
76.         float variance=t*t/(w*(1-w));    
77.        /* if(variance>maxVariance) 
78.         {   
79.             maxVariance=variance;   
80.             threshold=i;   
81.         }   
82.          */  
83.     return variance;    
84. }    
85. /***************************************************************** 
86. 函数名：Init 
87. 参数类型：void 
88. 功能：初始化粒子群算法的粒子与速度 
89. ************************************************************************/  
90. void Init()  
91. {  
92.     for(int index=0;index<number;index++)  
93.     {  
94.         antThreshold[index][0] = rnd(10 , 50);  
95.         antThreshold[index][1] = antThreshold[index][0] + 50;  
96.         if(antThreshold[index][1]>255)  
97.             antThreshold[index][1] = 255;  
98.         vect[index][0] = rnd(3 ,5);  
99.         vect[index][1] = rnd(3 ,5);  
100.     }  
101.           
102. }  
103. /****************************************************************** 
104. 函数名：Pso 
105. 参数类型：void 
106. 功能：粒子群算法的实现 
107. ***************************************************************************/  
108.   
109. void Pso(float value)  
110. {  
111.   for(int index=0;index<number;index++)  
112.   {  
113.       float variance;  
114.       variance = ThresholdOTSU(antThreshold[index][0] , antThreshold[index][1] , value);  
115.       if(variance>pbest[index])  
116.       {  
117.           pbest[index] = variance;  
118.           pbestThreshold[index][0] = antThreshold[index][0];  
119.           pbestThreshold[index][1] = antThreshold[index][1];  
120.       }  
121.       if(variance>gbest)  
122.       {  
123.           gbest = variance;  
124.           gbestThreshold[0] = antThreshold[index][0];  
125.           gbestThreshold[1] = antThreshold[index][1];  
126.       }  
127.   }  
128. }  
129. /***************************************************************************************88 
130. 函数名：updateData 
131. 参数类型：void 
132. 功能：更新粒子数据与速度 
133. **************************************************************************************************/  
134. void updateData()  
135. {  
136.     for(int index=0;index<number;index++)  
137.     {  
138.         for(int i=0;i<2;i++)  
139.         {  
140.             vect[index][i] = w*vect[index][i] + c1*((double)(rand())/(double)RAND_MAX)*(pbestThreshold[index][i]-antThreshold[index][i])+  
141.                 c2*c1*((double)(rand())/(double)RAND_MAX)*(gbestThreshold[i]-antThreshold[index][i]);  
142.             if(vect[index][i]>5)  
143.                 vect[index][i] = 5;  
144.             if(vect[index][i]<3)  
145.                 vect[index][i] = 3;  
146.             antThreshold[index][i] = vect[index][i] + antThreshold[index][i];  
147.         }  
148.         if(antThreshold[index][0]>antThreshold[index][1])  
149.             antThreshold[index][1] = antThreshold[index][0] + 50;  
150.         if(antThreshold[index][1]>255)  
151.             antThreshold[index][1] = 255;  
152.         if(antThreshold[index][0]<0)  
153.             antThreshold[index][0] = 0;  
154.     }  
155.   
156. }  
157. /**************************************************************8 
158. 函数名：Threshold 
159. 参数类型：IplImage *src , int lower , int higher 
160. 功能：利用算法得到的双阈值对图像进行阈值分割 
161. ***********************************************************************/  
162. void Threshold(IplImage *src , int lower , int higher)  
163. {  
164.     assert(src->nChannels==1);  
165.     for(int h=0;h<src->height;h++)  
166.         for(int w=0;w<src->width;w++)  
167.         {  
168.             if(*(src->imageData+h*src->widthStep+w)<higher&&*(src->imageData+h*src->widthStep+w)>lower)  
169.                 //*(src->imageData+h*src->widthStep+w) = 255;  
170.                 ;  
171.             else  
172.                 *(src->imageData+h*src->widthStep+w) = 0;  
173.         }  
174. }  
175. int _tmain(int argc, _TCHAR* argv[])  
176. {  
177.     srand((unsigned)time(NULL));  
178.     IplImage *img =0;  
179.     IplImage *ycrcb = 0;  
180.     IplImage *cb = 0;  
181.     cvNamedWindow("cb" , 1);  
182.     img = cvLoadImage("1.jpg" , 1);  
183.     ycrcb = cvCreateImage(cvGetSize(img) , 8 ,3);  
184.     cb = cvCreateImage(cvGetSize(img) , 8 , 1);  
185.   
186.     cvCvtColor(img , ycrcb , CV_BGR2YCrCb);  
187.     cvSplit(ycrcb , 0 ,0,cb , 0);  
188.   
189.     cvSmooth(cb , cb , CV_MEDIAN , 3 , 0,0,0);  
190.     float avgValue = 0.0;  
191.     avgValue = GetAvgValue(cb);  
192.     Init();  
193.     for(int i=0;i<3000;i++)  
194.     {  
195.        Pso(avgValue);  
196.        updateData();  
197.     }  
198.   
199.     //cvThreshold(cb , cb , gbestThreshold[0] , gbestThreshold[1] , CV_THRESH_BINARY);  
200.     Threshold(cb , gbestThreshold[0] , gbestThreshold[1]);  
201.     printf("%d , %d\n" ,  gbestThreshold[0] , gbestThreshold[1]);  
202.     cvShowImage("cb" , cb);  
203.     cvSaveImage("cb1.jpg" ,cb);  
204.     cvWaitKey(0);  
205.   
206.     return 0;  
207. }  

 

网址：http://blog.csdn.net/qinjianganying/article/details/6747980