opencv编写的二维OTSU算法

/********************************** 二维Otsu 算法，算法见论文 二维Otsu自适应阈值选取算法的快速实现 作者 汪海洋 潘德炉 夏德深 ***********************************/ #include "cv.h" #include "highgui.h" #include<iostream> using namespace std; void TwoDOTSU(IplImage*p1,IplImage*p2); inline float funcp(int m,int n); inline float funcx(int m,int n); inline float funcy(int m,int n); unsigned int graytab[256][256]={0}; float pList[256][256]={0}; float xList[256][256]={0}; float yList[256][256]={0}; float pi[256][256]={0}; float trace[256][256]={0}; int main() { IplImage*p1=0; IplImage*p2=0; p1=cvLoadImage("C:\\Documents and Settings\\i_mandy\\桌面\\image processing\\lena.bmp",0); p2=cvCloneImage(p1); //高斯平滑得到平滑图像 cvSmooth(p2,p2,CV_GAUSSIAN,3,3); TwoDOTSU(p1,p2); cvReleaseImage(&p1); cvReleaseImage(&p2); return 0; } void TwoDOTSU(IplImage*p1,IplImage*p2) { /*得到二维直方图*/ uchar* data; uchar* data1; CvSize ss,ss1; int step,step1; cvGetRawData( p1, &data, &step,&ss); cvGetRawData( p2, &data1, &step1,&ss1); int i,j,k; uchar* lp; uchar*lp1; lp=data; lp1=data1; for(i = 0; i< ss.width * ss.height; i++) { graytab[(unsigned char)*lp][(unsigned char)*lp1]++; lp++; lp1++; } float uy0=0; float uy1=0; for(j=0;j<256;j++) { for(k=0;k<256;k++) { pi[j][k]=(float)graytab[j][k]/p1->imageSize; uy0+=pi[j][k]*j; uy1+=pi[j][k]*k; } } /*建立查找表*/ for(i=0;i<256;i++) { for(j=0;j<256;j++) { pList[i][j]=funcp(i,j); xList[i][j]=funcx(i,j); yList[i][j]=funcy(i,j); } } /*计算迹*/ float w0,w1,x0,x1,y0,y1; int s,t; for(s=0;s<256;s++) { for(t=0;t<256;t++) { w0=pList[s][t]; w1=pList[255][255]-pList[255][t]-pList[s][255]+pList[s][t]; x0=xList[s][t]; y0=yList[s][t]; x1=xList[255][255]-xList[255][t]-xList[s][255]+xList[s][t]; y1=yList[255][255]-yList[255][t]-yList[s][255]+yList[s][t]; float temp0,temp1; temp0=pow((x0-uy0),2)+pow((y0-uy1),2); temp1=pow((x1-uy0),2)+pow((y1-uy1),2); trace[s][t]=temp0*w0+temp1*w1; } } /*排序得到最大值*/ float max=trace[0][0]; int tr1,tr2; for(s=0;s<256;s++) { for(t=0;t<256;t++) { if(max<trace[s][t]) { max=trace[s][t]; tr1=s; tr2=t; } } } cout<<tr1<<endl; cout<<tr2<<endl; } float funcp(int m,int n) { float temp=0; if(m==0||n==0) { if(m==0&&n==0) return pi[0][0]; else if(m==0&&n!=0) { for(int i=0;i<n;i++) temp+=pi[0][i]; return temp; } else { for(int j=0;j<m;j++) temp+=pi[j][0]; return temp; } } else { for(int t=0;t<m;t++) for(int k=0;k<n;k++) temp+=pi[t][k]; return temp; } } float funcx(int m,int n) { float temp=0; if(m==0||n==0) { if(m==0) return 0; if(m!=0) { for(int i=0;i<m;i++) temp+=i*pi[i][0]; return temp; } } else { for(int t=0;t<m;t++) for(int k=0;k<n;k++) temp+=t*pi[t][k]; return temp; } } float funcy(int m,int n) { float temp=0; if(m==0||n==0) { if(n==0) return 0; if(n!=0) { for(int i=0;i<n;i++) temp+=i*pi[0][i]; return temp; } } else { for(int t=0;t<m;t++) for(int k=0;k<n;k++) temp+=k*pi[t][k]; return temp; } }