图像处理中的标记分水岭分割算法

图像处理中的标记分水岭分割算法
如果图像中的目标物体是连接在一起的，则分割起来会更困难，分水岭分割算法经常用于处理这类问题，通常会取得比较好的效果。

分水岭分割算法把图像看成一幅“地形图”，其中亮度比较强的区域像素值较大，而比较暗的区域像素值较小，通过寻找“汇水盆地”和“分水岭界限”，对图像进行分割。

直接应用分水岭分割算法的效果往往并不好，如果在图像中对前景对象和背景对象进行标注区别，再应用分水岭算法会取得较好的分割效果。

有很多图像处理工具箱函数可以用到，如fspecial、imfilter、watershed、lable2rgb、imopen、imclose、imreconstruct、imcomplement、imregionalmax、bwareaopen、graythresh、和imimposemin函数等。

下面进行一个例子，步骤如下。

1、读取图像并求其边界，代码如下。

rgb = imread('pears.png');%读取原图像
I = rgb2gray(rgb);%转化为灰度图像
figure; subplot(121)%显示灰度图像
imshow(I)
text(732,501,'Image courtesy of Corel',...
'FontSize',7,'HorizontalAlignment','right')
hy = fspecial('sobel');%sobel算子
hx = hy';
Iy = imfilter(double(I), hy, 'replicate');%滤波求Y方向边缘
Ix = imfilter(double(I), hx, 'replicate');%滤波求X方向边缘
gradmag = sqrt(Ix.^2 + Iy.^2);%求模
subplot(122); imshow(gradmag,[]), %显示梯度
title('Gradient magnitude (gradmag)')
在这一步骤中，首先读取一套真彩色图像，然后把真色图像转化为灰度图像，结果如图所示：
图1原图和梯度图像
使用sobel 边缘算子对图像进行水平和垂直方向的滤波，然后求取模值，sobel 算子滤波后的图像在边缘处会显示比较大的值，在没有边界处的值会很小，如上图右图所示。

2、 直接使用梯度模值进行分水岭算法，代码如下。

L = watershed(gradmag);%直接应用分水岭算法
Lrgb = label2rgb(L);%转化为彩色图像 figure; imshow(Lrgb), %显示分割后的图像
title(
'Watershed transform of gradient magnitude (Lrgb)')
直接使用梯度模值图像进行分水岭算法得到的结果往往会存在过度分割的现象，如下图所示。

因此通常需要对前景对象和背景对象进行标记，以获得更好的分割效果。

Image courtesy of Corel
Gradient magnitude (gradmag)
图2直接对梯度图像进行分水岭分割
3、 分别对前景和背景进行标记，代码如下。

se = strel('disk', 20);%圆形结构元素
Io = imopen(I, se);%形态学开操作
figure; subplot(121)
imshow(Io), %显示执行开操作后的图像
title('Opening (Io)')
Ie = imerode(I, se);%对图像进行腐蚀
Iobr = imreconstruct(Ie, I);%形态学重建
subplot(122); imshow(Iobr), %显示重建后的图像
title('Opening-by-reconstruction (Iobr)')
Ioc = imclose(Io, se);%形态学关操作
figure; subplot(121)
imshow(Ioc), %显示关操作后的图像
title('Opening-closing (Ioc)')
Iobrd = imdilate(Iobr, se);%对图像进行膨胀
Iobrcbr = imreconstruct(imcomplement(Iobrd), ...
imcomplement(Iobr));%形态学重建
Iobrcbr = imcomplement(Iobrcbr);%图像求反
subplot(122); imshow(Iobrcbr), %显示重建求反后的图像
title('Opening-closing by reconstruction (Iobrcbr)')
fgm = imregionalmax(Iobrcbr);%局部极大值
figure; imshow(fgm), %显示重建后局部极大值图像
title('Regional maxima of opening-closing
by reconstruction
(fgm)')
Watershed transform of gradient magnitude (Lrgb)
I2 = I;
I2(fgm) = 255;%局部极大值处像素设置为255
figure; imshow(I2), %在原图上显示极大值区域
title('Regional maxima superimposed on original image (I2)')
se2 = strel(ones(5,5));%结构元素
fgm2 = imclose(fgm, se2);%关操作
fgm3 = imerode(fgm2, se2);%腐蚀
fgm4 = bwareaopen(fgm3, 20);%开操作
I3 = I;
I3(fgm4) = 255;%前景处设置为255
figure; subplot(121)
imshow(I3)%显示修改后的极大值区域
title('Modified regional maxima')
bw = im2bw(Iobrcbr, graythresh(Iobrcbr));%转化为二值图像
subplot(122); imshow(bw), %显示二值图像
title('Thresholded opening-closing by reconstruction') 可以使用很多算法来对前景对象进行标记，标记的每个对象内部的像素值是连接在一起的。

在这个例子中，
使用形态学重建技术对前景对象进行标记，首先使用imopen
函数对图像进行开操作，如图3左图所示，使用半径为20的圆形结构元素，开操作是膨胀和腐蚀操作的结合。

另外一种方法是先对图像进行腐蚀，然后对图像进行形态学重建，处理后的图像如图3右图所示。

图3开操作和重建操作结果对比
Opening (Io)Opening-by-reconstruction (Iobr)
在开操作之后进行关键操作可以去除一些很小的目标，如苹果上的茎干和一些污点等，
如图4
中左图所示。

另外一种方法是先对图像进行腐蚀，然后对图像进行形态学重建，重建后的图像如图4右图所示，注意在重建之前需要先对图像求反，然后再重建之后再进行一次求反。

图4关操作和重建操作结果对比
对比两幅图像Iobrcbr 和Ioc ，以重建为基础的开关操作（结果为Iobrcbr ）比一般的开关操作（结果为Ioc ）比一般的开关操作（结果为Ioc ）在去除小的污点时会更有效，并且不会影响这个图像的轮廓。

计算Iobrcbr 的局部极大值会得到比较好的前景标记，如图5
所示。

Opening-closing (Ioc)Opening-closing by reconstruction (Iobrcbr)
图5求取局部极大值的图像
为了更好的理解这个结果，可以在原图像的基础上，显示局部极
大值，对前景图像进行标记，如图6所示。

图6在原图上显示局部极大值
Regional maxima of opening-closing by reconstruction (fgm)
Regional maxima superimposed on original image (I2)
注意到图像中还有少部分目标物体，即苹果，未被正确的标记出，如果这些目标物体不能被正确的进行标记，则不能正确的进行分割。

并且，少部分前景目标物体已经拓展到边缘，因此应该收缩一下边缘，可以先对图像进行操作，然后进行腐蚀来达到这样的效果。

这个过程会产生一些孤立的像素点，可以使用
bwareaopen 函数
来达到这样的效果，这个过程将像素点数量较少的孤立像素点去除，如图7左图所示。

图7调整后的局部极大值图像和二值图像
将图像Iobrcbr 使用合适的阈值转化成二值图像，其中淡颜色的值为背景，转化成二值图像如图7右图所示。

4、 进行分水岭变换并显示，代码如下。

D = bwdist(bw);%计算距离
DL = watershed(D);%分水岭变换
Modified regional maxima Thresholded opening-closing by reconstruction
bgm = DL == 0;%求取分割边界
figure; imshow(bgm), %显示分割后的边界
title('Watershed ridge lines (bgm)')
gradmag2 = imimposemin(gradmag, bgm | fgm4);%设置最小值
L = watershed(gradmag2);%分水岭变换
I4 = I;
I4(imdilate(L == 0, ones(3, 3)) | bgm | fgm4) = 255;%前景及
边界处设置为255
figure; subplot(121)
imshow(I4)%突出前景及边界
title('Markers and object boundaries')
Lrgb = label2rgb(L, 'jet', 'w', 'shuffle');%转化为伪色彩图像
subplot(122); imshow(Lrgb)%显示伪色彩图像
title('Colored watershed label matrix')
figure; imshow(I),
hold on
himage = imshow(Lrgb);%在原图上显示伪色彩图像
set(himage, 'AlphaData', 0.3);
title('Lrgb superimposed transparently on original image')
从图7中可以看出，背景像素是黑色的，但在理想情况下，我们不希望背景标记太靠近目标对象的边缘，可以通过“骨骼化”进行细化，对二值图像的距离进行分水岭变换，然后寻找分水岭的界限，分水岭的界限如图8所示。

图8分水岭界限
Watershed ridge lines (bgm)
函数imimpesemin 可以用来对图像进行修改，使图像在一些特定的
区域像素值最小。

在本程序中，
使用imimpesemin 函数修改梯度模值图像，使梯度模值图像在标记的前景对象和背景对象中有最小值，然后进行分水岭变换得到矩阵L 。

还有一种可视化的技术是在原图像中分别标记前景对象，背景对象和边界。

为了使分割的边界更清楚，可以对图像进行膨胀操作，如图9左图所示。

从图9左图中可以看出，对前景对象和背景对象标记后再进行分水岭变换比直接在梯度模值图像上进行分水岭变换的效果要好得多。

另外一种显示分割图像的方法是使用彩色图像显示，使用
label2rgb 函数将分水岭变换后的图像显示为彩色，如图9右图所示。

图9分割图像的显示
Markers and object boundaries Colored watershed label matrix
同样也可以在原图的基础上，使用透明技术将图像显示为伪彩色图像，如图10所示。

图10 分割图像的伪彩色显示
Lrgb superimposed transparently on original image。