前面学习了直方图双峰法：图像二值化方法中的阈值法
最大类间方差法（OTSU）是找到自适应阈值的常用方法。原理参考了冈萨雷斯的《数字图像处理》。
以下是自己写的函数：
//----获取灰度图in的OTSU阈值--
int Segment::otsuMat(Mat in)
{
int i,j;
int temp;
//第一类均值，第二类均值，全局均值，mk=p1*m1, 第一类概率，第二类概率
double m1,m2,mG,mk,p1,p2;
int hist[256] = {0};
double pro_hist[256] = {0.0};
double cov;
double maxcov=0.0;
int maxthread=0;
int row = in.rows;
int col = in.cols;
//统计每个灰度的数量
for(i=0; i<row; ++i) {
for(j=0; j<col; ++j) {
temp = in.at<uchar>(i,j);
hist[temp]++;
}
}
//计算每个灰度级占图像的概率
for(i=0; i<256; ++i)
pro_hist[i] = (double)hist[i]/(double)(row*col);
//计算平均灰度值
mG = 0.0;
for(i=0; i<256; ++i)
mG += i * pro_hist[i];
//统计前景和背景的平均灰度值，并计算类间方差
for(i=0; i<256; ++i)
{
m1=0.0; m2=0.0; mk=0.0; p1=0.0; p2=0.0;
for(j=0; j<i; ++j) {
p1 += pro_hist[j];
mk += j*pro_hist[j];
}
m1 = mk/p1; //mk=p1*m1,是一个中间值
p2 = 1 - p1; //p1+p2=1;
m2 = (mG - mk)/p2; //mG=p1*m1+p2*m2;
//计算类间方差
cov = p1*p2*(m1-m2)*(m1-m2);
if(cov>maxcov) {
maxcov = cov;
maxthread = i;
}
}
//cout<<maxthread<<endl;
return maxthread;
}灰度图：
对应的二值图：
图1,2,4的结果都可以接受，但是图3的结果突然就差了很远，即使我在OTSU法算得的阈值基础上加加减减也没办法兼顾到每一张图。
我列出了灰度图的均值和方差，以及OTSU算得的阈值
这四张图的均值和方差很接近，但是OTSU算得的阈值却变化较大，尤其是图3，导致二值化结果不正确。
觉得自己对“最大类间方差”理解有误，至少，它不能解决这一类问题。它的结果 out of my control. 琢磨琢磨，到底是怎么回事？
最后通过均值和标准差设置了阈值，效果可控。
另外，要使用OTSU，也可以直接采用opencv 库函数：
double threshold(InputArray src, OutputArray dst, double thresh, double maxval, int type)当 type = CV_THRESH_OTSU 时, thresh 的值会被忽略，所以可以任取一个值。
在手册上，type 中没有列出OTSU方法，但是在函数源代码中是有的：double cv::threshold( InputArray _src, OutputArray _dst, double thresh, double maxval, int type )
{
Mat src = _src.getMat();
bool use_otsu = (type & THRESH_OTSU) != 0;
type &= THRESH_MASK;
if( use_otsu )
{
CV_Assert( src.type() == CV_8UC1 );
thresh = getThreshVal_Otsu_8u(src);
}
_dst.create( src.size(), src.type() );
Mat dst = _dst.getMat();
if( src.depth() == CV_8U )
{
int ithresh = cvFloor(thresh);
thresh = ithresh;
int imaxval = cvRound(maxval);
if( type == THRESH_TRUNC )
imaxval = ithresh;
imaxval = saturate_cast<uchar>(imaxval);
if( ithresh < 0 || ithresh >= 255 )
{
if( type == THRESH_BINARY || type == THRESH_BINARY_INV ||
((type == THRESH_TRUNC || type == THRESH_TOZERO_INV) && ithresh < 0) ||
(type == THRESH_TOZERO && ithresh >= 255) )
{
int v = type == THRESH_BINARY ? (ithresh >= 255 ? 0 : imaxval) :
type == THRESH_BINARY_INV ? (ithresh >= 255 ? imaxval : 0) :
/*type == THRESH_TRUNC ? imaxval :*/ 0;
dst.setTo(v);
}
else
src.copyTo(dst);
return thresh;
}
thresh = ithresh;
maxval = imaxval;
}
else if( src.depth() == CV_16S )
{
int ithresh = cvFloor(thresh);
thresh = ithresh;
int imaxval = cvRound(maxval);
if( type == THRESH_TRUNC )
imaxval = ithresh;
imaxval = saturate_cast<short>(imaxval);
if( ithresh < SHRT_MIN || ithresh >= SHRT_MAX )
{
if( type == THRESH_BINARY || type == THRESH_BINARY_INV ||
((type == THRESH_TRUNC || type == THRESH_TOZERO_INV) && ithresh < SHRT_MIN) ||
(type == THRESH_TOZERO && ithresh >= SHRT_MAX) )
{
int v = type == THRESH_BINARY ? (ithresh >= SHRT_MAX ? 0 : imaxval) :
type == THRESH_BINARY_INV ? (ithresh >= SHRT_MAX ? imaxval : 0) :
/*type == THRESH_TRUNC ? imaxval :*/ 0;
dst.setTo(v);
}
else
src.copyTo(dst);
return thresh;
}
thresh = ithresh;
maxval = imaxval;
}
else if( src.depth() == CV_32F )
;
else
CV_Error( CV_StsUnsupportedFormat, "" );
parallel_for_(Range(0, dst.rows),
ThresholdRunner(src, dst, thresh, maxval, type),
dst.total()/(double)(1<<16));
return thresh;
}/* Threshold types */
enum
{
CV_THRESH_BINARY =0, /* value = value > threshold ? max_value : 0 */
CV_THRESH_BINARY_INV =1, /* value = value > threshold ? 0 : max_value */
CV_THRESH_TRUNC =2, /* value = value > threshold ? threshold : value */
CV_THRESH_TOZERO =3, /* value = value > threshold ? value : 0 */
CV_THRESH_TOZERO_INV =4, /* value = value > threshold ? 0 : value */
CV_THRESH_MASK =7,
CV_THRESH_OTSU =8 /* use Otsu algorithm to choose the optimal threshold value;
combine the flag with one of the above CV_THRESH_* values */
};static double
getThreshVal_Otsu_8u( const Mat& _src )
{
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;
j = 0;
#if CV_ENABLE_UNROLLED
for( ; 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]++;
}
#endif
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*(double)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;
}