使用分水岭法分割图像

何为分水岭法？

传统的分水岭分割方法，是一种基于拓扑理论的数学形态学的分割方法，其基本思想是把图像看作是测地学上的拓扑地貌，图像中每一像素的灰度值表示该点的海拔高度，每一个局部极小值及其影响区域称为集水盆地，而集水盆地的边界则形成分水岭。分水岭的概念和形成可以通过模拟浸入过程来说明。在每一个局部极小值表面，刺穿一个小孔，然后把整个模型慢慢浸人水中，随着浸入的加深，每一个局部极小值的影响域慢慢向外扩展，在两个集水盆汇合处构筑大坝如下图所示，即形成分水岭。如图所示：

image.png

然而基于梯度图像的直接分水岭算法容易导致图像的过分割，产生这一现象的原因主要是由于输入的图像存在过多的极小区域而产生许多小的集水盆地，从而导致分割后的图像不能将图像中有意义的区域表示出来。所以必须对分割结果的相似区域进行合并。

因为传统分水岭算法存在过分割的不足，OpenCV提供了一种改进的分水岭算法，使用一系列预定义标记来引导图像分割的定义方式。使用OpenCV的分水岭算法cv::wathershed，需要输入一个标记图像，图像的像素值为32位有符号正数（CV_32S类型），每个非零像素代表一个标签。它的原理是对图像中部分像素做标记，表明它的所属区域是已知的。分水岭算法可以根据这个初始标签确定其他像素所属的区域。传统的基于梯度的分水岭算法和改进后基于标记的分水岭算法示意图。如下图所示：

image.png

从上图可以看出，传统基于梯度的分水岭算法由于局部最小值过多造成分割后的分水岭较多。而基于标记的分水岭算法，水淹过程从预先定义好的标记图像（像素）开始，较好的克服了过度分割的不足。本质上讲，基于标记点的改进算法是利用先验知识来帮助分割的一种方法。因此，改进算法的关键在于如何获得准确的标记图像，即如何将前景物体与背景准确的标记出来。

主要API：

// 分水岭算法
// 参数说明：
// image —— 输入图像矩阵，要求是8位无符号3通道彩色图像。
// markers —— 在执行分水岭函数watershed之前，必须对第二个参数markers进行处理，它应该包含不同区域的轮廓，每个轮廓有一个自己唯一的编号，轮廓的定位可以通过Opencv中findContours方法实现，这个是执行分水岭之前的要求。
void watershed( InputArray image, InputOutputArray markers );

接下来执行分水岭会发生什么呢？算法会根据markers传入的轮廓作为种子（也就是所谓的注水点），对图像上其他的像素点根据分水岭算法规则进行判断，并对每个像素点的区域归属进行划定，直到处理完图像上所有像素点。而区域与区域之间的分界处的值被置为“-1”，以做区分。

简单概括一下就是说第二个入参markers必须包含了种子点信息。Opencv官方例程中使用鼠标划线标记，其实就是在定义种子，只不过需要手动操作，而使用findContours可以自动标记种子点。而分水岭方法完成之后并不会直接生成分割后的图像，还需要进一步的显示处理，如此看来，只有两个参数的watershed其实并不简单。

使用方法：

分水岭法被用来分割图像，分割结果基于markers的标记指导。对markers的标记方法，在日常工作中，根据实际需要，有以下方法：

基于轮廓提取的标记法：

相关步骤包括：

将图像转换为单通道灰度图。
执行模糊滤波，减小细小的边缘。
使用Canny算法或其它梯度算子，求取二值边缘图像。
视情况决定是否执行形态学操作。
对二值边缘图像提取轮廓。
根据轮廓，制作标记掩模。
执行分水岭过程调用。
提取分割结果。

下面举一个例子，边看代码边看执行结果。

// 使用轮廓标记法，执行分水岭算法
{
    cv::Mat matImg = cv::imread("space_shuttle.jpg");

    // 边缘检测
    cv::Mat matImgGray;
    cv::cvtColor(matImg, matImgGray, cv::COLOR_BGR2GRAY);
    cv::GaussianBlur(matImgGray, matImgGray, cv::Size(7, 7), 0);
    cv::Mat matCanny;
    cv::Canny(matImgGray, matCanny, 100, 200);

    // 根据边缘找轮廓
    std::vector<std::vector<cv::Point> > contours;
    std::vector<cv::Vec4i> hierarchy;
    cv::findContours(matCanny, contours, hierarchy, cv::RETR_TREE, cv::CHAIN_APPROX_SIMPLE);

    // 对轮廓进行标记
    cv::Mat matMarks(matImg.size(), CV_32SC1, cv::Scalar(0));
    for (int i = 0; i < contours.size(); ++i)
    {
        cv::drawContours(matMarks, contours, i, cv::Scalar(i + 1), 1, 8, hierarchy);
    }

    // 分水岭区域标记
    cv::watershed(matImg, matMarks);

    // 对分水岭区域标记结果按分块颜色显示
    cv::Mat matCuts(matImg.size(), matImg.type(), cv::Scalar(0, 0, 0));

    for (int y = 0; y < matMarks.rows; ++y)
    {
        for (int x = 0; x < matMarks.cols; ++x)
        {
            int nMarkVal = matMarks.at<int>(y, x);

            // 轮廓区域
            if (nMarkVal == -1)
            {
                matCuts.at<cv::Vec3b>(y, x) = cv::Vec3b(255, 255, 255);
            }
            // 标记区域
            else
            {
                matCuts.at<cv::Vec3b>(y, x) = randomColor(nMarkVal);
            }
        }
    }
}

上面的randomColor(int)函数定义为：

// 根据标记值生成随机颜色
cv::Vec3b randomColor(int n)
{
    static std::map<int, cv::Vec3b> mapNcolor;

    auto iter = mapNcolor.find(n);
    if (iter != mapNcolor.end())
        return iter->second;

    cv::RNG rng(n);
    int B = rng.uniform(1, 256);
    int G = rng.uniform(1, 256);
    int R = rng.uniform(1, 256);

    cv::Vec3b color(B, G, R);

    mapNcolor[n] = color;
    return color;
}

首先加载图像，还是那架航天飞机：

image.png

然后使用Canny算法提取边缘，为了提取较少的边缘，上面代码使用了7x7的高斯算子，同时抬高了Canny边缘检测的高低阀值。提取到的边缘为：

image.png

大家可以尝试改变算子尺寸和Canny的阀值试试其他效果。

接下来使用cv::findContours()函数寻找轮廓，并对轮廓进行标记。注意标记区域必须大于1，所有为0的像素将被视为未知区域，分水岭将进行注水识别。标记掩模图如下：

image.png

看起来与边缘图没有差别，这是由于掩模图格式是32位整数，我们放大其中某个细节看标记情况：

image.png

继续执行分水岭注水操作，完成对所有未知区域的标记。注水标记图如下：

image.png

图中的黑色部分为-1，既原来的掩模标记，而掩模中的0未知区域都被成功标记，放大细节：

image.png

最后是对结果上色，便于查看结果，也可以根据需要，裁出指定部分的分块图。上色结果如下：

image.png

基于前景背景的标记法：

相关步骤包括：

将图像转换为单通道灰度图。
二值化提取前景和背景。
视情况决定是否执行形态学操作。
求取前景的骨架图质心图。
对骨架质心图求取连通域标注。
执行分水岭过程调用。
提取分割结果。

举例代码如下：

// 使用连通域标记法，执行分水岭算法
{
    cv::Mat matImg = cv::imread("image31.png");

    cv::Mat matImgGray;
    cv::cvtColor(matImg, matImgGray, cv::COLOR_BGR2GRAY);

    // 提取暗像素做为前景
    cv::Mat matFg;
    cv::threshold(matImgGray, matFg, 30, 255, cv::THRESH_BINARY_INV);

    // 形态学操作，适当放大前景
    cv::Mat matKernel = cv::getStructuringElement(cv::MORPH_RECT, cv::Size(3, 3));
    cv::morphologyEx(matFg, matFg, cv::MORPH_CLOSE, matKernel, cv::Point(-1, -1), 3);
    cv::dilate(matFg, matFg, matKernel, cv::Point(-1, -1), 3);

    // 求取前景骨架二值图像
    cv::Mat matDistance;
    cv::distanceTransform(matFg, matDistance, cv::DIST_L2, cv::DIST_MASK_5);
    cv::Mat matDistance2;
    cv::normalize(matDistance, matDistance2, 0, 1.0, cv::NORM_MINMAX);
    cv::Mat matFg2;
    double fMax = 0.0;
    cv::minMaxLoc(matDistance, NULL, &fMax);
    cv::threshold(matDistance, matFg2, fMax * 0.4, 255, cv::THRESH_BINARY);
    matFg2.convertTo(matFg2, CV_8UC1);

    // 基于前景骨架二值图像求连通域标注
    cv::Mat matMarks;
    cv::connectedComponents(matFg2, matMarks);

    // 执行分水岭算法
    cv::watershed(matImg, matMarks);

    // 对分水岭区域标记结果按分块颜色显示
    cv::Mat matCuts(matImg.size(), matImg.type(), cv::Scalar(0, 0, 0));

    for (int y = 0; y < matMarks.rows; ++y)
    {
        for (int x = 0; x < matMarks.cols; ++x)
        {
            int nMarkVal = matMarks.at<int>(y, x);

            // 轮廓区域
            if (nMarkVal == -1)
            {
                matCuts.at<cv::Vec3b>(y, x) = cv::Vec3b(255, 255, 255);
            }
            // 标记区域
            else
            {
                matCuts.at<cv::Vec3b>(y, x) = randomColor(nMarkVal);
            }
        }
    }
}

首先还是加载图像，由于航天飞机从灰度值上不太好切分，换个图吧：