matlab中threshold,OpenCV中threshold自动阈值，类似matlab中的graythresh

在opencv工程里面使用otsu分割灰度图像类似于matlab里的graythresh,

opencv里面提供了otsu threshold的源代码“icvGetThreshVal_Otsu”，我们可以根据源代码进行自动阈值处理

icvGetThreshVal_Otsu的源代码如下：

//opencv 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;

}对源代码进行改写，得到我们想要的代码：

double getThreshVal_Otsu_8u( const cv::Mat& _src )

{

cv::Size size = _src.size();

if ( _src.isContinuous() )

{

size.width *= size.height;

size.height = 1;

}

const int N = 256;

int i, j, h[N] = {0};

for ( i = 0; i < size.height; i++ )

{

const uchar* src = _src.data + _src.step*i;

for ( j = 0; j <= size.width - 4; j += 4 )

{

int v0 = src[j], v1 = src[j+1];

h[v0]++; h[v1]++;

v0 = src[j+2]; v1 = src[j+3];

h[v0]++; h[v1]++;

}

for ( ; j < size.width; j++ )

h[src[j]]++;

}

double mu = 0, scale = 1./(size.width*size.height);

for ( i = 0; i < N; i++ )

mu += i*h[i];

mu *= scale;

double mu1 = 0, q1 = 0;

double max_sigma = 0, max_val = 0;

for ( i = 0; i < N; i++ )

{

double p_i, q2, mu2, sigma;

p_i = h[i]*scale;

mu1 *= q1;

q1 += p_i;

q2 = 1. - q1;

if ( std::min(q1,q2) < FLT_EPSILON || std::max(q1,q2) > 1. - FLT_EPSILON )

continue;

mu1 = (mu1 + 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 = i;

}

}

return max_val;

}直接进行调用上述函数，即可：

m_threshold = getThreshVal_Otsu_8u(matMedian);在进行二值化：

threshold(matMedian,matBinary,m_threshold,255,0); // 二值化

详细请参考OpenCV