Otsu算法步骤如下：

设图象包含L个灰度级(0,1…,L-1)，灰度值为i的的象素点数为Ni ，图象总的象素点数为N=N0+N1+...+N(L-1)。灰度值为i的点的概为：P(i) = N(i)/N.门限t将整幅图象分为暗区c1和亮区c2两类，则类间方差σ是t的函数：σ=a1*a2(u1-u2)^2 (2)式中，aj 为类cj的面积与图象总面积之比，a1 = sum(P(i)) i->t, a2 = 1-a1; uj为类cj的均值，u1 = sum(i*P(i))/a1 0->t, u2 = sum(i*P(i))/a2, t+1->L-1 该法选择最佳门限t^ 使类间方差最大，即：令Δu=u1-u2，σb = max{a1(t)*a2(t)Δu^2}

代码实现：

parameter: *image --- buffer for image

rows, cols --- size of imagex0, y0, dx, dy --- region of vector used for computing thresholdvvv --- debug option, is 0, no debug information outputed*//*======================================================================*//* 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 *//*======================================================================*/int otsu (IplImage *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 = &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 }}*/for(j=y0;j<dy;j++){     for(i=0;i<dx;i++)    {          unsigned char temp=CV_IMAGE_ELEM(image,uchar,j,i);         ihist[temp]++;     }}

// 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 problemsfprintf (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);

}