阈值分割总结

#include <assert.h>
#include "cv.h"
#include "highgui.h"
#include <math.h> 
#include <iostream>
using namespace std;
// implementation of otsu algorithm
// author: onezeros(@yahoo.cn)
// reference: Rafael C. Gonzalez. Digital Image Processing Using MATLAB
void cvThresholdOtsu(IplImage* src, IplImage* dst)
{
	int height=src->height;
	int width=src->width;	

	//histogram
	float histogram[256]={0};
	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];
	}

	int threshold;	
	float maxVariance=0;
	float w=0,u=0;
	for(int i=0;i<256;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;
		}
	}

	cvThreshold(src,dst,threshold,255,CV_THRESH_BINARY);
}
// implementation of otsu algorithm
// author: onezeros(@yahoo.cn)
// reference: Rafael C. Gonzalez. Digital Image Processing Using MATLAB
int cvThresholdOtsu(IplImage* src)
{
	int height=src->height;
	int width=src->width;	

	//histogram
	float histogram[256]={0};
	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];
	}

	int threshold;	
	float maxVariance=0;
	float w=0,u=0;
	for(int i=0;i<256;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 threshold;
}
#include <cv.h>
#include <cxcore.h>
#include <highgui.h>
#pragma comment(lib,"cv210d.lib")
#pragma comment(lib,"cxcore210d.lib")
#pragma comment(lib,"highgui210d.lib")

#include <iostream>
using namespace std;

int main()
{	


	int threshold=-1;
	IplImage* img =cvLoadImage("E:\\test.jpg");
		cvShowImage("video",img);
		cvCvtColor(img,img,CV_RGB2YCrCb);

		IplImage* imgCb=cvCreateImage(cvGetSize(img),8,1);
		cvSplit(img,NULL,NULL,imgCb,NULL);
		if (threshold<0){
			threshold=cvThresholdOtsu(imgCb);
		}
		//cvThresholdOtsu(imgCb,imgCb);
		cvThreshold(imgCb,imgCb,threshold,255,CV_THRESH_BINARY);
       cvSaveImage("E:\\imgCb.bmp",imgCb);
		cvShowImage("object",imgCb);
		cvReleaseImage(&imgCb);
	return 0;
}

%基于贝叶斯分类算法的图像阈值分割 
clear 
clc; 
Init = imread('E:\\test.jpg'); 
Im=rgb2gray(Init);
figure(1)
imhist(Im),title('直方图') ; 
[x,y]=size(Im);                    % 求出图象大小 
b=double(Im);                   
zd=double(max(Im))                  % 求出图象中最大的灰度 
zx=double(min(Im))                % 最小的灰度  
T=double((zd+zx))/2;                     % T赋初值，为最大值和最小值的平均值 

count=double(0);                         % 记录几次循环 
while 1                   % 迭代最佳阈值分割算法 
    count=count+1; 
    S0=0.0; n0=0.0;                   %为计算灰度大于阈值的元素的灰度总值、个数赋值 
    S1=0.0; n1=0.0;                   %为计算灰度小于阈值的元素的灰度总值、个数赋值 
    for i=1:x
        for j=1:y
            if double(Im(i,j))>=T
                S1=S1+double(Im(i,j));  %大于阈域值图像点灰度值累加
                n1=n1+1;                %大于阈域值图像点个数累加
            else 
                S0=S0+double(Im(i,j));  %小于阈域值图像点灰度值累加
                n0=n0+1;                %小于阀域值图像点个数累加
            end 
        end 
    end  
    T0=S0/n0; %求小于阀域值均值
    T1=S1/n1; %求大于阀域值均值
    if abs(T-((T0+T1)/2))<0.1     %迭代至 前后两次阀域值相差几乎为0时 停止迭代。
        break;
    else
       T=(T0+T1)/2;                 %在阈值T下，迭代阈值的计算过程 
   end 
end 

count                                %显示运行次数
T
i1=im2bw(Im,T/255);               % 图像在最佳阈值下二值化 
figure(2)
imshow(i1,'border','tight','InitialMagnification','fit') 
title('实验结果') ;

*===============================图像分割=====================================*/
*---------------------------------------------------------------------------*/
#include <assert.h>
#include "cv.h"
#include "highgui.h"
#include <math.h> 
#include <iostream>
using namespace std;
int HistogramBins = 256;  
float HistogramRange1[2]={0,255};  
float *HistogramRange[1]={&HistogramRange1[0]};  
typedef  
enum {back,object} entropy_state;  
*======================================================================*/
* 迭代法*/
*======================================================================*/
// nMaxIter：最大迭代次数；nDiffRec：使用给定阀值确定的亮区与暗区平均灰度差异值
int DetectThreshold(IplImage*img, int nMaxIter, int& iDiffRec)  //阀值分割：迭代法
{
	//图像信息
	int height = img->height;
	int width = img->width;
	int step = img->widthStep/sizeof(uchar);
	uchar *data = (uchar*)img->imageData;

	iDiffRec =0;
	int F[256]={ 0 }; //直方图数组
	int iTotalGray=0;//灰度值和
	int iTotalPixel =0;//像素数和
	byte bt;//某点的像素值

	uchar iThrehold,iNewThrehold;//阀值、新阀值
	uchar iMaxGrayValue=0,iMinGrayValue=255;//原图像中的最大灰度值和最小灰度值
	uchar iMeanGrayValue1,iMeanGrayValue2;

	//获取(i,j)的值，存于直方图数组F
	for(int i=0;i<width;i++)
	{
		for(int j=0;j<height;j++)
		{
			bt = data[i*step+j];
			if(bt<iMinGrayValue)
				iMinGrayValue = bt;
			if(bt>iMaxGrayValue)
				iMaxGrayValue = bt;
			F[bt]++;
		}
	}

	iThrehold =0;//
	iNewThrehold = (iMinGrayValue+iMaxGrayValue)/2;//初始阀值
	iDiffRec = iMaxGrayValue - iMinGrayValue;

	for(int a=0;(abs(iThrehold-iNewThrehold)>0.5)&&a<nMaxIter;a++)//迭代中止条件
	{
		iThrehold = iNewThrehold;
		//小于当前阀值部分的平均灰度值
		for(int i=iMinGrayValue;i<iThrehold;i++)
		{
			iTotalGray += F[i]*i;//F[]存储图像信息
			iTotalPixel += F[i];
		}
		iMeanGrayValue1 = (uchar)(iTotalGray/iTotalPixel);
		//大于当前阀值部分的平均灰度值
		iTotalPixel =0;
		iTotalGray =0;
		for(int j=iThrehold+1;j<iMaxGrayValue;j++)
		{
			iTotalGray += F[j]*j;//F[]存储图像信息
			iTotalPixel += F[j];    
		}
		iMeanGrayValue2 = (uchar)(iTotalGray/iTotalPixel);

		iNewThrehold = (iMeanGrayValue2+iMeanGrayValue1)/2;        //新阀值
		iDiffRec = abs(iMeanGrayValue2 - iMeanGrayValue1);
	}

	//cout<<"The Threshold of this Image in imgIteration is:"<<iThrehold<<endl;
	return iThrehold;
}
*======================================================================*/
 /* OTSU global thresholding routine */
 /* takes a 2D unsigned char array pointer, number of rows, and */
 /* number of cols in the array. returns the value of the threshold */
 /*parameter: 
 *image --- buffer for image
 rows, cols --- size of image
 x0, y0, dx, dy --- region of vector used for computing threshold
 vvv --- debug option, is 0, no debug information outputed
 */
 /*
 OTSU 算法可以说是自适应计算单阈值（用来转换灰度图像为二值图像）的简单高效方法。
 下面的代码最早由 Ryan Dibble提供，此后经过多人Joerg.Schulenburg, R.Z.Liu 等修改，补正。
 算法对输入的灰度图像的直方图进行分析，将直方图分成两个部分，使得两部分之间的距离最大。
 划分点就是求得的阈值。
 */
 /*======================================================================*/
 int otsu (unsigned char*image, int rows, int cols, int x0, int y0, int dx, int dy, int vvv)
 {
     
     unsigned char*np; // 图像指针
 int thresholdValue=1; // 阈值
 int ihist[256]; // 图像直方图，256个点
     
 int i, j, k; // various counters
 int n, n1, n2, gmin, gmax;
 double m1, m2, sum, csum, fmax, sb;
     
 // 对直方图置零
     memset(ihist, 0, sizeof(ihist));
     
     gmin=255; gmax=0;
 // 生成直方图
 for (i = y0 +1; i < y0 + dy -1; i++) 
     {
         np = (unsigned char*)image[i*cols+x0+1];
 for (j = x0 +1; j < x0 + dx -1; j++)
         {
             ihist[*np]++;
 if(*np > gmax) gmax=*np;
 if(*np < gmin) gmin=*np;
             np++; /* next pixel */
         }
     }
     
 // set up everything
     sum = csum =0.0;
     n =0;
     
 for (k =0; k <=255; k++) 
     {
         sum += (double) k * (double) ihist[k]; /* x*f(x) 质量矩*/
         n += ihist[k]; /* f(x) 质量 */
     }
     
 if (!n) 
     {
 // if n has no value, there is problems...
         fprintf (stderr, "NOT NORMAL thresholdValue = 160\n");
 return (160);
     }
     
 // do the otsu global thresholding method
     fmax =-1.0;
     n1 =0;
 for (k =0; k <255; k++)
     {
         n1 += ihist[k];
 if (!n1) 
         { 
 continue; 
         }
         n2 = n - n1;
 if (n2 ==0)
         { 
 break; 
         }
         csum += (double) k *ihist[k];
         m1 = csum / n1;
         m2 = (sum - csum) / n2;
         sb = (double) n1 *(double) n2 *(m1 - m2) * (m1 - m2);
 /* bbg: note: can be optimized. */
 if (sb > fmax) 
         {
             fmax = sb;
             thresholdValue = k;
         }
     }
     
 // at this point we have our thresholding value
     
 // debug code to display thresholding values
 if ( vvv &1 )
         fprintf(stderr,"# OTSU: thresholdValue = %d gmin=%d gmax=%d\n",
         thresholdValue, gmin, gmax);
     
 return(thresholdValue);
 }
 /*======================================================================*/
 /* OTSU global thresholding routine */
 /*======================================================================*/
 int otsu2 (IplImage *image)
 {
	 int w = image->width;
	 int h = image->height;

	 unsigned char*np; // 图像指针
	 unsigned char pixel;
	 int thresholdValue=1; // 阈值
	 int ihist[256]; // 图像直方图，256个点

	 int i, j, k; // various counters
	 int n, n1, n2, gmin, gmax;
	 double m1, m2, sum, csum, fmax, sb;

	 // 对直方图置零...
	 memset(ihist, 0, sizeof(ihist));

	 gmin=255; gmax=0;
	 // 生成直方图
	 for (i =0; i < h; i++) 
	 {
		 np = (unsigned char*)(image->imageData + image->widthStep*i);
		 for (j =0; j < w; j++) 
		 {
			 pixel = np[j];
			 ihist[ pixel]++;
			 if(pixel > gmax) gmax= pixel;
			 if(pixel < gmin) gmin= pixel;
		 }
	 }

	 // set up everything
	 sum = csum =0.0;
	 n =0;

	 for (k =0; k <=255; k++) 
	 {
		 sum += k * ihist[k]; /* x*f(x) 质量矩*/
		 n += ihist[k]; /* f(x) 质量 */
	 }

	 if (!n) 
	 {
		 // if n has no value, there is problems...
		 //fprintf (stderr, "NOT NORMAL thresholdValue = 160\n");
		 thresholdValue =160;
		 goto L;
	 }

	 // do the otsu global thresholding method
	 fmax =-1.0;
	 n1 =0;
	 for (k =0; k <255; k++) 
	 {
		 n1 += ihist[k];
		 if (!n1) { continue; }
		 n2 = n - n1;
		 if (n2 ==0) { break; }
		 csum += k *ihist[k];
		 m1 = csum / n1;
		 m2 = (sum - csum) / n2;
		 sb = n1 * n2 *(m1 - m2) * (m1 - m2);
		 /* bbg: note: can be optimized. */
		 if (sb > fmax)
		 {
			 fmax = sb;
			 thresholdValue = k;
		 }
	 }

L:
	 for (i =0; i < h; i++) 
	 {
		 np = (unsigned char*)(image->imageData + image->widthStep*i);
		 for (j =0; j < w; j++) 
		 {
			 if(np[j] >= thresholdValue)
				 np[j] =255;
			 else np[j] =0;
		 }
	 }

	 //cout<<"The Threshold of this Image in Otsu is:"<<thresholdValue<<endl;
	 return(thresholdValue);
 }
 /*============================================================================
 =  代码内容：最大熵阈值分割                                      
 =  修改日期:2009-3-3                                                                                                         
 =  作者:crond123 
 =  博客:http://blog.csdn.net/crond123/
 =  E_Mail:crond123@163.com                                                      
 ===============================================================================*/
 // 计算当前位置的能量熵
 double caculateCurrentEntropy(CvHistogram * Histogram1,int cur_threshold,entropy_state state)
 {
 int start,end;
 int  total =0;
 double cur_entropy =0.0;
 if(state == back)    
     {
         start =0;
         end = cur_threshold;    
     }
 else    
     {
         start = cur_threshold;
         end =256;    
     }    
 for(int i=start;i<end;i++)    
     {
         total += (int)cvQueryHistValue_1D(Histogram1,i);//查询直方块的值 P304
     }
 for(int j=start;j<end;j++)
     {
 if((int)cvQueryHistValue_1D(Histogram1,j)==0)
 continue;
 double percentage = cvQueryHistValue_1D(Histogram1,j)/total;
 /*熵的定义公式*/
         cur_entropy +=-percentage*logf(percentage);
 /*根据泰勒展式去掉高次项得到的熵的近似计算公式
         cur_entropy += percentage*percentage;*/        
     }
 return cur_entropy;
 //    return (1-cur_entropy);
 }
 
 //寻找最大熵阈值并分割
 void  MaxEntropy(IplImage *src,IplImage *dst)
 {
     assert(src != NULL);
     assert(src->depth ==8&& dst->depth ==8);
     assert(src->nChannels ==1);
     CvHistogram * hist  = cvCreateHist(1,&HistogramBins,CV_HIST_ARRAY,HistogramRange);//创建一个指定尺寸的直方图
 //参数含义：直方图包含的维数、直方图维数尺寸的数组、直方图的表示格式、方块范围数组、归一化标志
     cvCalcHist(&src,hist);//计算直方图
 double maxentropy =-1.0;
 int max_index =-1;
 // 循环测试每个分割点，寻找到最大的阈值分割点
 for(int i=0;i<HistogramBins;i++)    
     {
 double cur_entropy = caculateCurrentEntropy(hist,i,object)+caculateCurrentEntropy(hist,i,back);
 if(cur_entropy>maxentropy)
         {
             maxentropy = cur_entropy;
             max_index = i;
         }
     }
     cout<<"The Threshold of this Image in MaxEntropy is:"<<max_index<<endl;
     cvThreshold(src, dst, (double)max_index,255, CV_THRESH_BINARY);
     cvReleaseHist(&hist);
 }
 /*============================================================================
 =  代码内容：基本全局阈值法                              
 ==============================================================================*/
 int BasicGlobalThreshold(int*pg,int start,int end)
 {                                           //  基本全局阈值法
 int  i,t,t1,t2,k1,k2;
 double u,u1,u2;    
     t=0;     
     u=0;
 for (i=start;i<end;i++) 
     {
         t+=pg[i];        
         u+=i*pg[i];
     }
     k2=(int) (u/t);                          //  计算此范围灰度的平均值    
 do 
     {
         k1=k2;
         t1=0;    
         u1=0;
 for (i=start;i<=k1;i++) 
         {             //  计算低灰度组的累加和
             t1+=pg[i];    
             u1+=i*pg[i];
         }
         t2=t-t1;
         u2=u-u1;
 if (t1) 
             u1=u1/t1;                     //  计算低灰度组的平均值
 else 
             u1=0;
 if (t2) 
             u2=u2/t2;                     //  计算高灰度组的平均值
 else 
             u2=0;
         k2=(int) ((u1+u2)/2);                 //  得到新的阈值估计值
     }
 while(k1!=k2);                           //  数据未稳定，继续
 //cout<<"The Threshold of this Image in BasicGlobalThreshold is:"<<k1<<endl;
 return(k1);                              //  返回阈值
 }
int main()
{
	/*手动设置阀值*/
	IplImage* smoothImgGauss = cvLoadImage("E:\\test.jpg", CV_LOAD_IMAGE_GRAYSCALE);
	IplImage* imgGrey = cvCreateImage(cvGetSize(smoothImgGauss),IPL_DEPTH_8U, 1);
	int w = smoothImgGauss->width, h = smoothImgGauss->height;
	IplImage* binaryImg = cvCreateImage(cvGetSize(smoothImgGauss),IPL_DEPTH_8U, 1);
	cvThreshold(smoothImgGauss,binaryImg,71,255,CV_THRESH_BINARY); 
	//cvNamedWindow("cvThreshold", CV_WINDOW_AUTOSIZE );
	//cvShowImage( "cvThreshold", binaryImg );
	cvSaveImage("E:\\手动设置阀值.bmp",binaryImg);
	cvReleaseImage(&binaryImg);    
	/*---------------------------------------------------------------------------*/
	/*自适应阀值  //计算像域邻域的平均灰度，来决定二值化的值*/
	IplImage* adThresImg = cvCreateImage(cvGetSize(smoothImgGauss),IPL_DEPTH_8U, 1);
	double max_value=255;
	int adpative_method=CV_ADAPTIVE_THRESH_GAUSSIAN_C;//CV_ADAPTIVE_THRESH_MEAN_C
	int threshold_type=CV_THRESH_BINARY;
	int block_size=3;//阈值的象素邻域大小
	int offset=5;//窗口尺寸
	cvAdaptiveThreshold(smoothImgGauss,adThresImg,max_value,adpative_method,threshold_type,block_size,offset);
	//cvNamedWindow("cvAdaptiveThreshold", CV_WINDOW_AUTOSIZE );
	cvSaveImage("E:\\自适应阀值.bmp",adThresImg);
	//cvShowImage( "cvAdaptiveThreshold", adThresImg );
	cvReleaseImage(&adThresImg);

	/*---------------------------------------------------------------------------*/
	/*最大熵阀值分割法*/    
	IplImage* imgMaxEntropy = cvCreateImage(cvGetSize(imgGrey),IPL_DEPTH_8U,1);
	MaxEntropy(smoothImgGauss,imgMaxEntropy);
	//cvNamedWindow("MaxEntroyThreshold", CV_WINDOW_AUTOSIZE );
	//cvShowImage( "MaxEntroyThreshold", imgMaxEntropy );//显示图像
	cvSaveImage("E:\\最大熵阀值分割.bmp",imgMaxEntropy);
	cvReleaseImage(&imgMaxEntropy ); 
	
	/*---------------------------------------------------------------------------*/
	/*基本全局阀值法*/
	IplImage* imgBasicGlobalThreshold = cvCreateImage(cvGetSize(imgGrey),IPL_DEPTH_8U,1);
	IplImage* srcImgGrey = cvCreateImage(cvGetSize(imgGrey),IPL_DEPTH_8U,1);
	cvCopyImage(srcImgGrey,imgBasicGlobalThreshold);
	int  pg[256],i,thre;    
	for (i=0;i<256;i++) pg[i]=0;
	for (i=0;i<imgBasicGlobalThreshold->imageSize;i++)      //  直方图统计
		pg[(BYTE)imgBasicGlobalThreshold->imageData[i]]++;    
	thre = BasicGlobalThreshold(pg,0,256);    //  确定阈值
	cout<<"The Threshold of this Image in BasicGlobalThreshold is:"<<thre<<endl;//输出显示阀值
	cvThreshold(imgBasicGlobalThreshold,imgBasicGlobalThreshold,thre,200,CV_THRESH_BINARY);  //  二值化    
	//cvNamedWindow("BasicGlobalThreshold", CV_WINDOW_AUTOSIZE );
	//cvShowImage( "BasicGlobalThreshold", imgBasicGlobalThreshold);//显示图像
    cvSaveImage("E:\\基本全局阀值法.bmp",imgBasicGlobalThreshold);
	cvReleaseImage(&imgBasicGlobalThreshold);
	/*---------------------------------------------------------------------------*/
	/*OTSU*/
	IplImage* imgOtsu = cvCreateImage(cvGetSize(imgGrey),IPL_DEPTH_8U,1);
	cvCopyImage(srcImgGrey,imgOtsu);
	int thre2;
	thre2 = otsu2(imgOtsu);
	cout<<"The Threshold of this Image in Otsu is:"<<thre2<<endl;//输出显示阀值
	cvThreshold(imgOtsu,imgOtsu,thre2,200,CV_THRESH_BINARY);  //  二值化    
	//cvNamedWindow("imgOtsu", CV_WINDOW_AUTOSIZE );
	//cvShowImage( "imgOtsu", imgOtsu);//显示图像    
	cvSaveImage("E:\\OTSU.bmp",imgOtsu);
	cvReleaseImage(&imgOtsu);
	/*---------------------------------------------------------------------------*/
	/*上下阀值法：利用正态分布求可信区间*/
	IplImage* imgTopDown = cvCreateImage( cvGetSize(imgGrey), IPL_DEPTH_8U, 1 );
	cvCopyImage(srcImgGrey,imgTopDown);
	CvScalar mean ,std_dev;//平均值、 标准差
	double u_threshold,d_threshold;
	cvAvgSdv(imgTopDown,&mean,&std_dev,NULL);    
	u_threshold = mean.val[0] +2.5* std_dev.val[0];//上阀值
	d_threshold = mean.val[0] -2.5* std_dev.val[0];//下阀值
	//u_threshold = mean + 2.5 * std_dev; //错误
	//d_threshold = mean - 2.5 * std_dev;
	cout<<"The TopThreshold of this Image in TopDown is:"<<d_threshold<<endl;//输出显示阀值
	cout<<"The DownThreshold of this Image in TopDown is:"<<u_threshold<<endl;
	cvThreshold(imgTopDown,imgTopDown,d_threshold,u_threshold,CV_THRESH_BINARY_INV);//上下阀值
	//cvNamedWindow("imgTopDown", CV_WINDOW_AUTOSIZE );
	//cvShowImage( "imgTopDown", imgTopDown);//显示图像  
	cvSaveImage("E:\\上下阀值法.bmp",imgTopDown);
	cvReleaseImage(&imgTopDown);
}

/************************************************************************/
/* 全局阈值分割  自动求取阈值        */
/************************************************************************/
//自动求取阈值，增加对场景的适应性
//只需求取一次，之后就可以一直使用
#include<cv.h>
#include <highgui.h>
#include <iostream>
#include <math.h>
using namespace std;
int main(){
	IplImage * image,* image2;
	image = cvLoadImage("E:\\111.jpg",0);
	cvNamedWindow("image",1);
	cvShowImage("image",image);
	image2 = cvCreateImage(cvSize(image->width,image->height),image->depth,1);
	double T = 0;
	double dT0 = 1.0;//阈值求取结束标志
	double dT = 255.0;

	//求取平均灰度，作为阈值T的初始值T0
	int i, j;
	double T0 = 0,T1 = 0,T2 = 0;//初始阈值
	int count1,count2;
	unsigned char * ptr,*dst;
	for (i = 0 ; i< image->height ; i++)
	{
		for (j =0 ; j < image->width;j++)
		{
			ptr = (unsigned char *)image->imageData + i*image->widthStep + j;
			T0 += ((double)(*ptr))/image->width/image->height;
		}
	}
	cout<<"T0:     "<<T0<<endl;
	T = (int)(T0 + 0.5); 
	//计算T两侧的灰度平均值，然后把二者的均值赋值给T
	while (dT > dT0)
	{

		T1 = 0;
		T2 = 0;
		count1 = 0;
		count2 = 0;
		for (i = 0 ; i< image->height ; i++)
		{
			for (j =0 ; j < image->width;j++)
			{
				ptr = (unsigned char *)image->imageData + i*image->widthStep + j;
				if (*ptr > T)
				{
					T1 += ((double)(*ptr))/image->width/image->height;
					count1++;
				} 
				else if(*ptr < T)
				{
					T2 +=  ((double)(*ptr))/image->width/image->height;
					count2++;
				}
			}
		}

		T1 = T1*image->width*image->height/count1;
		T2 = T2*image->width*image->height/count2;
		dT = fabs(T - (T1 + T2)/2);

		cout<<"T1"<<T1<<endl;
		cout<<"T2"<<T2<<endl;
		cout<<"dT  " << dT<<endl;
		T = (T1 + T2)/2;
		cout<<"T:     "<<T<<endl;
	}



	//根据求取的阈值进行分割
	for (i = 0 ; i< image2->height ; i++)
	{
		for (j =0 ; j < image2->width;j++)
		{
			ptr = (unsigned char *)image->imageData + i*image->widthStep + j;
			dst = (unsigned char *)image2->imageData+i*image2->widthStep+j;
			if (*ptr > T)
			{
				*dst = 255;
			} 
			else
			{
				*dst =0;
			}
		}
	}

	cvNamedWindow("image2",1);
	cvShowImage("image2",image2);
	cvSaveImage("E:\\image\\dowels2.tif",image2);
	cvWaitKey(0);
	return 0;
}