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

  

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

原图：

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

将背景像素置0：

效果图：

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

效果图：

 

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

源代码：

#include "stdafx.h"#include "cv.h"#include "highgui.h"#include "cxcore.h"#include "time.h"using namespace std;#define rnd( low,uper) ((int)(((double)rand()/(double)RAND_MAX)*((double)(uper)-(double)(low))+(double)(low)+0.5))/*************************************************************8888粒子群算法变量的说明******************************************************************************/const int number = 20;int antThreshold[number][2];//以阈值作为粒子int vect[number][2];//更新的速度float pbest[number] = {0.0};;//每个粒子历史最优解float gbest = 0.0;//全局历史最优解int pbestThreshold[number][2];//每个粒子的最优历史阈值int gbestThreshold[2];//全局粒子的最优阈值float w = 0.9;//惯性因子float c1 = 2.0;//加速因子1float c2 = 2.0;//加速因子2//histogram   float histogram[256]={0};  /*********************************************************************8888函数名：GetAvgValue参数类型：IplImage* src实现功能：获得灰度图像的总平均灰度值*****************************************************************************/float GetAvgValue(IplImage* src)  {      int height=src->height;      int width=src->width;            for(int i=0;i<height;i++) {          unsigned char* p=(unsigned char*)src->imageData+src->widthStep*i;          for(int j=0;j<width;j++) {              histogram[*p++]++;          }      }      //normalize histogram       int size=height*width;      for(int i=0;i<256;i++) {          histogram[i]=histogram[i]/size;      }        //average pixel value       float avgValue=0;      for(int i=0;i<256;i++) {          avgValue+=i*histogram[i];      }  return avgValue;}/*****************************************************************************函数名：ThresholdOTSU参数类型：int threshold1 , int threshold2 , float avgValue功能：求得最大类间方差**********************************************************************************/float  ThresholdOTSU(int threshold1 , int threshold2 , float avgValue){      int threshold;        float maxVariance=0;      float w=0,u=0;      for(int i=threshold1;i< threshold2 ;i++){          w+=histogram[i];          u+=i*histogram[i];  }          float t=avgValue*w-u;          float variance=t*t/(w*(1-w));         /* if(variance>maxVariance){              maxVariance=variance;              threshold=i;          }           */    return variance;  }  /*****************************************************************函数名：Init参数类型：void功能：初始化粒子群算法的粒子与速度************************************************************************/void Init(){for(int index=0;index<number;index++){antThreshold[index][0] = rnd(10 , 50);antThreshold[index][1] = antThreshold[index][0] + 50;if(antThreshold[index][1]>255)antThreshold[index][1] = 255;vect[index][0] = rnd(3 ,5);vect[index][1] = rnd(3 ,5);}}/******************************************************************函数名：Pso参数类型：void功能：粒子群算法的实现***************************************************************************/void Pso(float value){  for(int index=0;index<number;index++)  {  float variance;  variance = ThresholdOTSU(antThreshold[index][0] , antThreshold[index][1] , value);  if(variance>pbest[index])  {  pbest[index] = variance;  pbestThreshold[index][0] = antThreshold[index][0];          pbestThreshold[index][1] = antThreshold[index][1];  }  if(variance>gbest)  {  gbest = variance;  gbestThreshold[0] = antThreshold[index][0];  gbestThreshold[1] = antThreshold[index][1];  }  }}/***************************************************************************************88函数名：updateData参数类型：void功能：更新粒子数据与速度**************************************************************************************************/void updateData(){for(int index=0;index<number;index++){for(int i=0;i<2;i++){vect[index][i] = w*vect[index][i] + c1*((double)(rand())/(double)RAND_MAX)*(pbestThreshold[index][i]-antThreshold[index][i])+c2*c1*((double)(rand())/(double)RAND_MAX)*(gbestThreshold[i]-antThreshold[index][i]);if(vect[index][i]>5)vect[index][i] = 5;if(vect[index][i]<3)vect[index][i] = 3;antThreshold[index][i] = vect[index][i] + antThreshold[index][i];}if(antThreshold[index][0]>antThreshold[index][1])antThreshold[index][1] = antThreshold[index][0] + 50;if(antThreshold[index][1]>255)antThreshold[index][1] = 255;if(antThreshold[index][0]<0)antThreshold[index][0] = 0;}}/**************************************************************8函数名：Threshold参数类型：IplImage *src , int lower , int higher功能：利用算法得到的双阈值对图像进行阈值分割***********************************************************************/void Threshold(IplImage *src , int lower , int higher){assert(src->nChannels==1);for(int h=0;h<src->height;h++)for(int w=0;w<src->width;w++){if(*(src->imageData+h*src->widthStep+w)<higher&&*(src->imageData+h*src->widthStep+w)>lower)//*(src->imageData+h*src->widthStep+w) = 255;;else*(src->imageData+h*src->widthStep+w) = 0;}}int _tmain(int argc, _TCHAR* argv[]){srand((unsigned)time(NULL));IplImage *img =0;IplImage *ycrcb = 0;IplImage *cb = 0;    cvNamedWindow("cb" , 1);img = cvLoadImage("1.jpg" , 1);ycrcb = cvCreateImage(cvGetSize(img) , 8 ,3);cb = cvCreateImage(cvGetSize(img) , 8 , 1);cvCvtColor(img , ycrcb , CV_BGR2YCrCb);cvSplit(ycrcb , 0 ,0,cb , 0);cvSmooth(cb , cb , CV_MEDIAN , 3 , 0,0,0);float avgValue = 0.0;avgValue = GetAvgValue(cb);Init();for(int i=0;i<3000;i++){       Pso(avgValue);   updateData();}//cvThreshold(cb , cb , gbestThreshold[0] , gbestThreshold[1] , CV_THRESH_BINARY);Threshold(cb , gbestThreshold[0] , gbestThreshold[1]);printf("%d , %d\n" ,  gbestThreshold[0] , gbestThreshold[1]);cvShowImage("cb" , cb);cvSaveImage("cb1.jpg" ,cb);cvWaitKey(0);return 0;}