【图像算法】七种常见阈值分割代码(Otsu、最大熵、迭代法、自适应阀值、手动、迭代法、基本全局阈值法)

代码
/*======================================================================*/
/* 迭代法*/
/*======================================================================*/
// 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);                              //  返回阈值
}