算法步骤
1、生成灰度直方图，并进行归一化，得到比例直方图。
2、根据比例直方图计算整幅图像的平均灰度$mu_0$。
3、从灰度0迭代到灰度255，每次迭代计算背景（这里将小于当前迭代灰度的部分视为背景）占整幅图像的比例$omega_1$；计算背景的平均灰度$mu_1$；计算前景和背景的类间方差$sigma^2 = frac{omega_1} { (1 - omega_1)} · (mu_1 - mu_0)^2 $。
4、将最大类间方差对应的灰度设置为阈值，并进行二值化。
算法原理
为了表述的方便，这里先定义一些符号：
| 像素数占总图像的比例 | 平均灰度 | 类间方差 | |
| 背景 | $omega_1$ | $mu_1$ | |
| 前景 | $omega_2$ | $mu_2$ | |
| 图像 | $omega_0$ | $mu_0$ | $sigma^2 $ |
类间方差指的是前景和背景之间的差异，显然该差异越大，说明分离度越好。数学上，类间方差的计算方法是：
$sigma^2 = omega_1 · (mu_1 - mu_0)^2 + omega_2 · (mu_2 - mu_0)^2$（式2-1）
但直接采用该公式进行编程，不免有些繁杂，因此需要先进行一定的化简。将（式2-1）展开得：
$sigma^2 = omega_1 · mu_1^2 + omega_2 · mu_2^2 - 2( omega_1 · mu_1 + omega_2 · mu_2)· mu_0 + mu_0^2 $（式2-2）
根据期望的数学定义式$E(X)=sum_{k=1 }^{infty}x_k·p_k$我们可以推知：
$mu_0 = omega_1 · mu_1 + omega_2 · mu_2 $（式2-3）
将（式2-3）代入到（式2-2）中，则$sigma^2 = omega_1 · mu_1^2 + omega_2 · mu_2^2 - mu_0^2 $。然后我们希望可以把前景的成分消除掉，因此再次利用（式2-3）和$omega_2 = 1 - omega_1$的关系进行替换：
$sigma^2 = omega_1 · mu_1^2 + frac{omega_2^2 · mu_2^2}{1 - omega_1} - mu_0^2 = omega_1 · mu_1^2 + frac{(mu_0 - omega_1 · mu_1)^2}{1 - omega_1} - mu_0^2 = frac{omega_1} { (1 - omega_1)} · (mu_1 - mu_0)^2$（式2-4）
使用（式2-4），我们就只需要统计当前迭代灰度以前的像素即可，这大大提升了程序的效率。此外，类间方差的最大点往往处于直方图的低谷处（如下图所示，红色点代表方差值），因此可以比较好地自动进行前景和背景分离。
算法复现
复现需要注意的一点是背景平均灰度的计算，用公式可以表示为：
$mu_1 = frac{ sum_{i=0}^{T} i · p_i}{ sum_{i=0}^{T} p_i}$
其中$i$表示迭代灰度，$p_i$表示$i$灰度级的比例，$T$为当前迭代灰度值。
unsigned char Ostu(int size,unsigned char *image) { int i; unsigned char threshold=0; //阈值 float variance=0; //类间方差 float maxvariance=0; //最大方差 float k=0; float w1=0; //背景像素比例 float u1=0; //背景平均灰度 float u0=0; //图像的平均灰度 float histogram[256]={0}; //直方图 for (i=0;i<size;i++){ histogram[*(image+i)]++; //像素直方图 } for (i=0;i<256;i++){ histogram[i]/=size; //比例直方图 u0+=histogram[i]*i; //计算图像平均灰度 } for (i=0;i<256;i++){ w1+=histogram[i]; //背景比例 k+=i*histogram[i]; u1=k/w1; //背景平均灰度 variance=w1/(1-w1)*(u1-u0)*(u1-u0); //求方差 if (variance>maxvariance){ maxvariance=variance; threshold=i; //将最大g相应的i值作为图像的全局阈值 } } return threshold; }