Opencv大津(otsu)的使用代码

在opencv工程里面使用otsu分割灰度图像类似于matlab里的graythresh, opencv里面提供了otsu threshold的源代码“icvGetThreshVal_Otsu”，但是我们却不能直接调用它，只能把它拷贝出来使用（这就是开源的好处）
//Input:
//pTile_pic: 灰度图像，*IplImage
//output:
//im_th: pTile_pic 的otsu阈值
otsu threshold:
int size = 256;
float range[] = {0, 255};
float *pRange[] = {&range[0]};
CvHistogram *p_imhist = cvCreateHist(1, &size, CV_HIST_ARRAY, pRange);
cvCalcHist(&pTile_pic, p_imhist);
double im_th = icvGetThreshVal_Otsu(p_imhist);

Source code:
static double
icvGetThreshVal_Otsu( const CvHistogram* hist )
{
    double max_val = 0;
   
    CV_FUNCNAME( "icvGetThreshVal_Otsu" );
    __BEGIN__;
    int i, count;
    const float* h;
    double sum = 0, mu = 0;
    bool uniform = false;
    double low = 0, high = 0, delta = 0;
    float* nu_thresh = 0;
    double mu1 = 0, q1 = 0;
    double max_sigma = 0;
    if( !CV_IS_HIST(hist) || CV_IS_SPARSE_HIST(hist) || hist->mat.dims != 1 )
        CV_ERROR( CV_StsBadArg,
        "The histogram in Otsu method must be a valid dense 1D histogram" );
    count = hist->mat.dim[0].size;
    h = (float*)cvPtr1D( hist->bins, 0 );
    if( !CV_HIST_HAS_RANGES(hist) || CV_IS_UNIFORM_HIST(hist) )
    {
        if( CV_HIST_HAS_RANGES(hist) )
        {
            low = hist->thresh[0][0];
            high = hist->thresh[0][1];
        }
        else
        {
            low = 0;
            high = count;
        }
        delta = (high-low)/count;
        low += delta*0.5;
        uniform = true;
    }
    else
        nu_thresh = hist->thresh2[0];
    for( i = 0; i < count; i++ )
    {
        sum += h;
        if( uniform )
            mu += (i*delta + low)*h;
        else
            mu += (nu_thresh[i*2] + nu_thresh[i*2+1])*0.5*h;
    }
   
    sum = fabs(sum) > FLT_EPSILON ? 1./sum : 0;
    mu *= sum;
    mu1 = 0;
    q1 = 0;
    for( i = 0; i < count; i++ )
    {
        double p_i, q2, mu2, val_i, sigma;
        p_i = h*sum;
        mu1 *= q1;
        q1 += p_i;
        q2 = 1. - q1;
        if( MIN(q1,q2) < FLT_EPSILON || MAX(q1,q2) > 1. - FLT_EPSILON )
            continue;
        if( uniform )
            val_i = i*delta + low;
        else
            val_i = (nu_thresh[i*2] + nu_thresh[i*2+1])*0.5;
        mu1 = (mu1 + val_i*p_i)/q1;
        mu2 = (mu - q1*mu1)/q2;
        sigma = q1*q2*(mu1 - mu2)*(mu1 - mu2);
        if( sigma > max_sigma )
        {
            max_sigma = sigma;
            max_val = val_i;
        }
    }
    __END__;
    return max_val;
}