ostu算法

//这个是输入jpg格式的，有一行注释掉的是直接输入灰度图的，如果用灰度图的，可以吧上边注释掉

#include <iostream>

#include "cv.h"
#include "highgui.h" 
using namespace std;


int nThreshold;


int cvOtsu2D(CvMat *pGrayMat);






int cvOtsu2D(CvMat *pGrayMat)  
 {  
     double dHistogram[256][256];        //建立二维灰度直方图  
     double dTrMatrix = 0.0;             //离散矩阵的迹
     int height = pGrayMat->rows;
     int width = pGrayMat->cols;
     int N = height*width;               //总像素数
     int i, j;  
     for(i = 0; i < 256; i++)  
     {  
         for(j = 0; j < 256; j++)  
             dHistogram[i][j] = 0.0;      //初始化变量  
     }  
     for(i = 0; i < height; i++)  
     {  
         for(j = 0; j < width; j++)  
         {  
             unsigned char nData1 = (unsigned char)cvGetReal2D(pGrayMat, i, j);//当前的灰度值  
             unsigned char nData2 = 0;                                                                                                                                
             int nData3 = 0;//注意9个值相加可能超过一个字节                                                                                    
             for(int m = i-1; m <= i+1; m++)  
             {  
                 for(int n = j-1; n <= j+1; n++)  
                 {  
                     if((m >= 0) && (m < height) && (n >= 0) && (n < width))  
                         nData3 += (unsigned char)cvGetReal2D(pGrayMat, m, n); //当前的灰度值                                                                        
                 }  
             }  
             nData2 = (unsigned char)(nData3 / 9);    //对于越界的索引值进行补零,邻域均值  
             dHistogram[nData1][nData2]++;  
         }  
     }  
     for(i = 0; i < 256; i++)  
         for(j = 0; j < 256; j++)  
             dHistogram[i][j] /= N;  //得到归一化的概率分布  
   
     double Pai = 0.0;      //目标区均值矢量i分量  
     double Paj = 0.0;      //目标区均值矢量j分量  
     double Pbi = 0.0;      //背景区均值矢量i分量  
     double Pbj = 0.0;      //背景区均值矢量j分量  
     double Pti = 0.0;      //全局均值矢量i分量  
     double Ptj = 0.0;      //全局均值矢量j分量  
     double W0 = 0.0;       //目标区的联合概率密度  
     double W1 = 0.0;       //背景区的联合概率密度  
     double dData1 = 0.0;  
     double dData2 = 0.0;  
     double dData3 = 0.0;  
     double dData4 = 0.0;   //中间变量  
     int nThreshold_s = 0;  
     int nThreshold_t = 0;  
     double temp = 0.0;     //寻求最大值  
     for(i = 0; i < 256; i++)  
     {  
         for(j = 0; j < 256; j++)  
         {  
             Pti += i*dHistogram[i][j];  
             Ptj += j*dHistogram[i][j];  
         }  
     }  
     for(i = 0; i < 256; i++)  
     {  
         for(j = 0; j < 256; j++)  
         {  
             W0 += dHistogram[i][j];  
             dData1 += i*dHistogram[i][j];  
             dData2 += j*dHistogram[i][j];  
   
             W1 = 1-W0;  
             dData3 = Pti-dData1;  
             dData4 = Ptj-dData2;  
           
             Pai = dData1 / W0;  
             Paj = dData2 / W0;  
             Pbi = dData3 / W1;  
             Pbj = dData4 / W1;   // 得到两个均值向量，用4个分量表示  
             dTrMatrix = ((W0 * Pti - dData1) * (W0 * Pti - dData1) + (W0 * Ptj - dData2) * (W0 * Ptj- dData2)) / (W0 * W1);  
             if(dTrMatrix > temp)  
             {  
                 temp = dTrMatrix;  
                 nThreshold_s = i;  
                 nThreshold_t = j;  
             }  
         }  
     }  
     nThreshold = (nThreshold_s + nThreshold_t) / 2;//返回阈值，有多种形式，可以单独某一个，也可                                                         //是两者的均值 
     return nThreshold; 
 }  
 int main( )
 {   
 IplImage *img1 = cvLoadImage("test1.jpg",1);//加载图片
 IplImage *img= cvCreateImage(cvGetSize(img1),img1->depth,1);//灰度图
  cout << img1->nChannels<< endl;
  cvCvtColor(img1,img,CV_RGB2GRAY);


 //IplImage *img = cvLoadImage("copyImage.bmp", 0);//加载图片
cvNamedWindow("Miss", CV_WINDOW_AUTOSIZE);//创建显示窗口
cvShowImage("Miss", img);//显示图片
cvNamedWindow("二维Otsu", CV_WINDOW_AUTOSIZE);
 CvMat *imgMat = cvCreateMat(img->height, img->width, CV_8UC1);//创建矩阵，其type为CV_8UC1
 cvConvert(img, imgMat);//将图像转换成矩阵形式
int t = cvOtsu2D(imgMat);//求二维阈值分割的阈值
cvThreshold(img, img, t, 255, CV_THRESH_BINARY);//用求得的阈值分割图像
cvShowImage("二维Otsu", img);
 cvWaitKey(0);//保持图像的显示窗口
cvReleaseImage(&img);//释放图像资源
cvReleaseMat(&imgMat);//释放矩阵资源
cvDestroyAllWindows();//销毁窗口
return 0;
 }