阈值分割：基本全局阀值法、上下阀值法、迭代法

/*---------------------------------------------------------------------------*/
/*基本全局阀值法*/
    IplImage* imgBasicGlobalThreshold = 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,255,CV_THRESH_BINARY);  //  二值化    
     cvNamedWindow("BasicGlobalThreshold", CV_WINDOW_AUTOSIZE );
    cvShowImage( "BasicGlobalThreshold", imgBasicGlobalThreshold);//显示图像
     cvReleaseImage(&imgBasicGlobalThreshold);
/*---------------------------------------------------------------------------*/
/*上下阀值法：利用正态分布求可信区间*/
    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);//显示图像    
    cvReleaseImage(&imgTopDown);
/*---------------------------------------------------------------------------*/
/*迭代法*/
    IplImage* imgIteration = cvCreateImage( cvGetSize(imgGrey), IPL_DEPTH_8U, 1 );
    cvCopyImage(srcImgGrey,imgIteration);
int thre3,nDiffRec;
    thre3 =DetectThreshold(imgIteration, 100, nDiffRec);
    cout<<"The Threshold of this Image in imgIteration is:"<<thre3<<endl;//输出显示阀值
    cvThreshold(imgIteration,imgIteration,thre3,255,CV_THRESH_BINARY_INV);//上下阀值
    cvNamedWindow("imgIteration", CV_WINDOW_AUTOSIZE );
    cvShowImage( "imgIteration", imgIteration);
    cvReleaseImage(&imgIteration);

/*======================================================================*/
/* 迭代法*/
/*======================================================================*/
// 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;
}

/*============================================================================
=  代码内容：基本全局阈值法                              
==============================================================================*/
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);                              //  返回阈值
}