OpenCV - 均值迭代分割

【题外话】:之前在博客中写过一篇“区域生长”的博客，区域生长在平时经常用到，也比较容易理解和代码实现，所以在很多情况下大家会选择这种方法。但是区域生长有一个最致命的点就是需要选取一个生长的种子点。
为了交流学习，同时也为了后面查阅方便，准备陆续将基于直方图的几种分割算法加以总结。

1、均值迭代算法的描述

对一幅图像$M$,均值迭代算法就是要迭代计算得到一个灰度值$T$，使得这个灰度值$T$将图像分成的两类$A,B$。满足条件：$A$类的均值和$B$类的均值，再求均值正好等于$T$。
用直方图可以很直观的描述：

图示表示某个图像的灰度直方图，均值迭代就是为找到一个$T$，使得$T$左侧的积分面积中均值$T_{1}$和右侧积分面的均值$T_{2}$相等。

2、算法的步骤

1. 选择一个初始化阈值$T$，通常取整张图灰度值的平均值；
2. 计算$T$分成的两个部分的灰度均值$u_{1}$和$u_{2}$；
3. 更新$T＝（u_{1} + u_{2}）/2$；
4. 重复步骤2～3，直到相邻两次计算的结果相等，或者两次结果的差值小于预先设定的值某个值；
5. 用这个$T$对图像进行分割。

3、OpenCV下的实现

//Mat src： 待分割灰度图像
//int n：   初始阈值
void IsodataSeg(Mat &src, int n)
{
    int threshold = 0;                          //历史阈值
    int MeansO = 0;                             //前景灰度均值
    int nObject = 0;                            //实质像素点个数
    int MeansB = 0;                             //背景灰度均值
    int nBack = 0;                              //背景像素点个数
    int nCol = src.cols * src.channels();       //每行的像素个数
    while (abs(threshold - n) > 10)             //迭代停止条件
    {
        threshold = n;
        for (int i=0; i<src.rows; ++i)
        {
            uchar* pData = src.ptr<uchar>(i);
            for (int j=0; j<nCol; ++j)
            {
                if (pData[j] < threshold)   //背景
                {
                    ++nBack;
                    MeansB += pData[j];
                }
                else                        //物体
                {
                    ++nObject;
                    MeansO += pData[j];
                }
            }
        }
        if (nBack == 0 || nObject == 0)     //防止出现除以0的计算
            continue;
        n = (MeansB/nBack + MeansO/nObject) / 2;
    }
    cv::threshold(src, src, n, 255, 0);     //进行阈值分割
}

4、进一步探究

在实际的应用中常常会发现一幅图像中，我们真正想去分割的并不是整个的矩形。比如下面的图：

假设我们想要将整个奇怪的人的眼睛，嘴巴和他们脸分开，而并不关心整个图的背景部分。（当然有很多种方法可以直接去掉黑色的背景，现在仅仅讨论整个图）在算法的思路中我们可以看到，我们统计的是整个图的像素，试想如果我们不去统计哪些我们已经知道是背景颜色的像素点，那么整个问题不就转化成了我们已经能解决的问题了么？！
所以，我们仅仅在代码的实现中讲哪些我们已知是不要统计的颜色排除在外就可以了。
这里我们假定背景是纯黑的（像素值为0）。这样上述代码仅仅需要作很小的改动就能满足应用的要求了。

void IsodataSeg(Mat &src, int n)
{
    int threshold = 0;                          //历史阈值
    int MeansO = 0;                             //前景灰度均值
    int nObject = 0;                            //实质像素点个数
    int MeansB = 0;                             //背景灰度均值
    int nBack = 0;                              //背景像素点个数
    int nCol = src.cols * src.channels();       //每行的像素个数

    while (abs(threshold - n) > 10)             //迭代停止条件
    {
        threshold = n;
        for (int i=0; i<src.rows; ++i)
        {
            uchar* pData = src.ptr<uchar>(i);
            for (int j=0; j<nCol; ++j)
            {
                //黑色区域为多余的像素，不参与计算
                if (pData[j] == 0)
                    continue;

                if (pData[j] < threshold)   //背景
                {
                    ++nBack;
                    MeansB += pData[j];
                }
                else                        //物体
                {
                    ++nObject;
                    MeansO += pData[j];
                }
            }
        }
        if (nBack == 0 || nObject == 0)
            continue;
        n = (MeansB/nBack + MeansO/nObject) / 2;
    }
    cv::threshold(src, src, n, 255, 0);                             //进行阈值分割
}

5、均值迭代在哪好用

从算法的描述步骤中可以看到，均值迭代算法的收敛速度是很快的。也就意味了在合适的场景下，使用均值迭代算法是具有绝对优势的。但是也需要注意一点均值迭代一般比较使用直方图为典型的“双峰”的图像。这一特征也是基于直方图统计类算法的共同特征。
可以体会到均值迭代算法其实是一种简单的聚类，它每次将样本分成两类，两类的均值作为两类的聚类中心，每一次计算的 $T$ 与两类均值的距离为聚类的半径。