灰度图像的自动阈值分割（Otsu 法）

灰度图像的自动阈值分割（Otsu 法）

机器视觉领域许多算法都要求先对图像进行二值化。这种二值化操作阈值的选取非常重要。阈值选取的不合适，可能得到的结果就毫无用处。今天就来讲讲一种自动计算阈值的方法。这种方法被称之为Otsu法。发明人是个日本人，叫做Nobuyuki Otsu (大津展之）。
简单的说，这种算法假设一副图像由前景色和背景色组成，通过统计学的方法来选取一个阈值，使得这个阈值可以将前景色和背景色尽可能的分开。或者更准确的说是在某种判据下最优。与数理统计领域的 fisher 线性判别算法其实是等价的。

otsu算法中这个判据就是最大类间方差 (intra-class variance or the variance within the class)。下面就来详细说说什么是 intra-class variance。

我们知道一副灰度图像，可以计算它的颜色平均值，或者更进一步。可以计算出灰度直方图。

比如下面的例子图片：

这个图片拍摄的是一个条形码。在这个图中，前景色就是黑色的条形码，背景色是其余部分的灰色。那么我们可以计算出这个图像的灰度直方图。

图中那个大的峰是背景色的部分，小的峰是前景色。

灰度值的均值是 122. 我们称这个均值为 M。

现在任意选取一个灰度值 t，则可以将这个直方图分成前后两部分。我们称这两部分分别为 A 和 B。对应的就是前景色和背景色。这两部分各自的平均值成为 MA 和 MB。
A 部分里的像素数占总像素数的比例记作 PA，B部分里的像素数占总像素数的比例记作 PB。
Nobuyuki Otsu 给出的类间方差定义为：

$$ICV = PA * (MA-M)^2 + PB * (MB - M )^2$$

那么这个最佳的阈值$t$ 就是使得 $ICV$ 最大的那个值。
对于上面的测试图像，我们可以遍历 t 的各种取值，计算 ICV。之后可以画出这样的ICV 曲线（绿色线条）：

可以看出，ICV 取最值的点确实将前景色和背景色分开了。
下面是个例子代码，用到了 Qt 的QImage。

int otsu(const QImage &image)
{
    double hist[256];
    normalizedHistogram(image, hist);

    double omega[256];
    double mu[256];

    omega[0] = hist[0];
    mu[0] = 0;
    for(int i = 1; i < 256; i++)
    {
        omega[i] = omega[i-1] + hist[i]; //累积分布函数
        mu[i] = mu[i-1] + i * hist[i];
    }
    double mean = mu[255];// 灰度平均值
    double max = 0;
    int k_max = 0;
    for(int k = 1; k < 255; k++)
    {
        double PA = omega[k]; // A类所占的比例
        double PB = 1 - omega[k]; //B类所占的比例
        double value = 0;
        if( fabs(PA) > 0.001 && fabs(PB) > 0.001)
        {       
            double MA = mu[k] / PA; //A 类的灰度均值
            double MB = (mean - mu[k]) / PB;//B类灰度均值
            value = PA * (MA - mean) * (MA - mean) + PB * (MB - mean) * (MB - mean);//类间方差

            if (value > max)
            {
                max = value;
                k_max = k;
            }
        }
        //qDebug() <<k << " " << hist[k] << " " << value;
    }
    return k_max;
}
bool normalizedHistogram(const QImage &image, double hist[256])
{
    for(int i = 0; i < 256; i++)
    {
        hist[i] = 0.0;
    }
    int height = image.height();
    int width = image.width();
    int N = height * width;
    if(image.format() != QImage::Format_Indexed8)
    {
        return false;
    }
    for(int i = 0; i < height; i++)
    {
        const quint8 *pData = (const quint8 *)image.constScanLine(i);
        for(int j = 0; j < width; j++)
        {
            ++hist[ pData[j] ];
        }
    }
    for(int i = 0; i < 256; i++)
    {
        hist[i] = hist[i] / N;
    }
    return true;
}

利用这个方法计算出的阈值做了二值化后得到图像如下：

可以看到效果很好。

Otsu 方法也不是万能的。当目标与背景的大小比例悬殊时，类间方差准则函数可能呈现双峰或多峰，此时效果不好。这时就要考虑其他的办法了。

其实，我们可以仔细观察 ICV 的计算公式。

$$ICV = PA * (MA-M)^2 + PB * (MB - M )^2$$

这里面 PA 和 PB 相当于是个前景色和背景色部分做个加权。当前景色或背景色有一个区域很小时。比如 PA 非常的小。那么这时计算出来的 t 就会和 B 区域很接近，这时的分割效果就会比较差。我们可以对ICV的公式进行一点小小的改造。

$$ICV = PA^\alpha * (MA-M)^2 + PB^\alpha * (MB - M )^2$$

这里的 $\alpha$ 可以取一个 0-1之间的值。比如上面的例子图片，如果我们取 $\alpha = 0.8$ 计算出的效果会更好一些。当然这个 $\alpha$ 值就要全凭经验来定了。