OpenCV分水岭watershed的应用注意

在VS2010，OpenCV进行分水岭的实现时。我遇到了一个问题：
在做好种子图和背景图后，也无法分隔开同一个背景框内的多个种子点。网上给的方法都是将背景点设置为灰度128，种子点设置为255，其他为0。然后前景背景叠加作为mark图。再调用watershed(原图，蒙版)。
但是，按照这个结果做出来的话，同一个128灰度的区域包围的多个255灰度的种子区域无法被区分开。通过阅读OpenCV源码，我发现了原因。

/*-----划分:如果其4邻域只有一个集水区那就归入其中，如果有多个那就标记为分水岭------*/  
        if( q[active_queue].first == 0 )    //当前层级处理完或无像素入队，换下一层  
        {  
            for( i = active_queue+1; i < NQ; i++ )  
                if( q[i].first )    //换这层  
                    break;  
            if( i == NQ )   //退出大循环  
                break;  
            active_queue = i;   //下面给这一层的像素分配标号  
        }  
        //取一个点  
        ws_pop( active_queue, mofs, iofs );   

        m = mask + mofs;    //mask的值  
        ptr = img + iofs;   //img的值  
        //4邻域只有一个区域，确定区域归属，要是处在多个集水区，成为分水岭  
        t = m[-1];  
        if( t > 0 )
            lab = t;  
        t = m[1];  
        if( t > 0 )  
        {  
            if( lab == 0 )
                lab = t;  
            else if( t != lab )
                lab = WSHED;  
        }  
        t = m[-mstep];  
        if( t > 0 )  
        {  
            if( lab == 0 )
                lab = t;  
            else if( t != lab )
                lab = WSHED;  
        }  
        t = m[mstep];  
        if( t > 0 )  
        {  
            if( lab == 0 )
                lab = t;  
            else if( t != lab )
                lab = WSHED;  
        }  
        assert( lab != 0 );  
        m[0] = lab;  
        if( lab == WSHED )  
            continue;   //继续  

在这一部分可以看出，在设置分水岭位置时，是通过比较上下左右四个点，是否有不为0且不同的两个灰度点，如果有，那么设置为WSHED（-1），也就是分水岭位置。
引用一个比较易懂专业的描述：如果mask图像中该像素的四邻域中存在两个不同的非0值，表示该点为两个注水盆地的边缘，即分水岭线，在mark图像中标记该点为-1；否则在mask图像中标记该点为四邻域中大于0的那个值（如果有多个大于0的值，则按照左右上下的顺序取）。
同理，对于同一个128灰度区域包围的几个不同的种子点区域，他们的灰度值要设置的不同，那样在扩张不同水区域接触的时候才能够分开。
其实，什么前景背景对于watershed函数来说根本就没有区别，就是不同灰度值能够表示不同的注水点。以这些注水点为起点开始分水岭算法的实现，而不是之前所想的128的背景图部分是不动的。

同样需要注意的还有，opencv自带的分水岭算法的涨水方式和传统的涨水方式不同，在一片论文上我看到有对OpenCV分水岭算法进行改进，指出OpenCV分水岭算法在判断下一步对何种像素点进行涨水时，依据的是已经被标记为-2，也就是mark边缘和mark边缘的像素差的绝对值的大小。将其放入q队列（0-255个位置分别记录对应灰度值大小的点的位置）中，再从小到大搜寻q队列，对找到的第一个灰度进行涨水。这个就会导致一些过漫水的问题。
那篇论文名称为：《OpenCV分水岭算法的改进及其在细胞分割中的应用_张羽》，我按照这篇论文中提出的方法进行对源码的修改了，效果也不尽人意，如果哪位看到这篇文章并成功修改，效果也不错，请务必联系我。（如果这篇论文的作者看到了，也请帮助一下我这个初学者，具体应该如何修改。）
以下附上不同灰度种子改前改后的结果
：

：

其实最后的结果也不是很好，很有可能是我种子点的给出以及原图像有点太差。不过意思到位了，大家看懂就行。
错误定义mark图展示

这里就是把所有的种子点的像素都设置成了255，就会出现上面第一幅图的结果。也就是没分开。

2018.8.21更新
我突然发现在我2017.5.31的实验室组内汇报的时候，我的总结里有一个根据别人的方法那里学来的新的方法进行改进的想法，就是使用距离变换与分水岭结合的方法，在分水岭出来一部分结果之后，给出局部最小，再加上些新的东西使它能够一定程度上区分开连通域。
在此声明这个方法是我在http://www.it165.net/pro/html/201312/8297.html看到的。我讲的肯定没有人家的好，我在这简单阐述一下我的改进过程和结果。

OpenCV分水岭算法描述如下：
初始化mark矩阵，生成最初的注水区域。
1. 设置mark图像的边框值为-1
2. 标记每个mark区域的边界为-2
3. 对于mark图像一个像素值，如果它本身值为0，但上下左右四邻域有一个像素值大于0，则把img图像中该点按照RGB高度值（实际上是梯度值）放入优先队列q的对应位置。
其中q是这样的一个数据结构：它是一个0到255的数组，每个元素分别是一个队列的头。当往q的同一个位置添加高度值时，就是找到它的位置，再往对应队列的后面添加元素。
之前的OpenCV分水岭结果在上面可以看到。
而Matlab:
Watershed在Matlab中直接输入图片，进行分水岭计算，原理是自动寻找像素值局部最低的地方作为注水口，然后开始注水，基于分水岭算法原理直接输出结果。
读入一张图
filename=’4-HF开运算后.bmp’;
bw = imread(filename);

距离变换反变换
D = -bwdist(~bw);

直接分水岭的话
Ld = watershed(D);
下图为原图和分水岭结果叠加图

这样的结果中包含了不少过分割。原因和上面提到的一样，局部最小太多。对于这种情况，通用的技巧是，在基于watershed的图像分割中使用imextendedmin这个函数过滤掉一些特别小（指区域小）的局部最小。然后，我们修改距离变换的结果，让滤波后的区域不会出现局部最小值，这个操作叫做“强制最小（minima imposition）”，可以用imimposemin这个函数实现。

使用imextendedmin将会只在我们希望分割的区块中间产生小点。然后我们使用imshowpair来将模板叠加到原图上。

mask = imextendedmin(D,2);
imshowpair(bw,mask,’blend’)

imextendedmin函数的介绍：
功能：扩展极小值变换。
用法：BW = imextendedmin(I,h)
D为制定矩阵，2为保留最小的几个值的点，以下分别为2，3，4，6的结果

计算扩展的极小值变换，这是在H -极小区域极小变换。H是一个非负标量。
区域极小是以恒定强度值连通的像素，其外部边界的像素都具有较高的价值。
默认情况下，imextendedmin使用二维图像8邻域，和三维图像26邻域
对于高维，imextendedmin使用conndef(ndims(I),’maximal’).
BW = imextendedmin(I,h,conn)
计算扩展极小值变换，在conn指定连接，conn可以有以下任何标值。
值 含义
二维的连接
4 4邻域
8 8邻域
三维的连接
6 6邻域
18 18邻域
26 26邻域
连接可以被定义为更普遍方式的任何维度通过了conn 3-by-3-by-…-by-3的0’s 和 1’s 矩阵。1值元素定义相对位置附近的conn的中心内容，注意，conn必须和它的中心元素对称。

最后，我们我们修改距离变换的结果，让其只在想要的位置具有局部最小。然后进行watershed。这就是我们想要的结果了。
D2 = imimposemin(D,mask);
Ld2 = watershed(D2);
bw3 = bw;
bw3(Ld2 == 0) = 0;
imshow(bw3)
imimposemin
功能：施加极小值
用法：I2 = imimposemin(I,BW)
使用形态重构修改强度图像I，使得它在BW非零的地方只有区域极小值。
BW是一个二值图像，和I的尺寸相同。(一般选用刚才imextendedmin的输出)
默认情况下，imimposemin对二维图像使用8连通域，对三维图像使用26连通域。
对于更高维数的情况，imimposemin使用conndef( ndims (I), ’minimum’)。
I2 = imimposemin(I,BW,conn)

D2 = imimposemin(D,mask);运算前后D，D2的输出结果

效果已经算是理想，和我能想到的最好结果一样