控制标记符分水岭分割算法

(2011-01-10 04:44:39)

转载▼

//示例流程

①原图像及默认的分水岭变换：

>> f=imread('rice_tophat.png');
>> h=fspecial('sobel')
h =
     1     2     1
     0     0     0
    -1    -2    -1
>> fd=double(f);
>> g=sqrt(imfilter(fd,h,'replicate').^2+imfilter(fd,h','replicate').^2);%梯度，具体参《梯度》博文
>> L=watershed(g);
>> wr=L==0;%分水岭脊线
>> subplot(121),imshow(f),subplot(122),imshow(wr)

从图中看到，过分割现象特别严重。这样的结果完全没有意义。
②计算图像中大量局部最小区域的位置：

>> rm=imregionalmin(f);whos rm
  Name        Size             Bytes  Class      Attributes

rm        253x251            63503  logical

>> figure,imshow(rm)

其中，白点标记的小区域均为“局部最小区域”——这块区域的亮度值相同且其周围的亮度值都比它大。

③计算图像的“扩展的最小变换”，即通过一个阈值来扩展“局部最小区域”(让更多的与最小区域相似的像素合并到最小区域里边来)：

>> im2=imextendedmin(f,45);%试了多次，阈值为45左右效果较好

>> figure,imshow(im2)
>>fim=f;fim(im2)=175;figure,imshow(fim)%以灰色来显示这些“局部最小区域”

（此时，所有亮的区域均为扩展的“局部最小区域”）
 
可以看到，将扩展的局部最小区域叠加到原图像上后的效果。灰度区域灰度值为175.

这个“扩展的局部最小区域”作为内部标记符集合。

//不对，内部标记应选在要提取的目标上，这幅图像中要提取的是米粒图像，故作如下修改：

>>im2＝~im2;%im2是 二值矩阵，取反即可

>>fim=f;fim(im2)=175;figure,imshow(fim,[])

可见，已将米粒区域修改为对应的“扩展最小区域”——这个作为内部标记符
④寻找外部标记符：

>> temp=bwdist(im2);%距离变换

>>figure,imshow(temp,[])

>>Lim=watershed(temp);

>>em=Lim==0;
>>figure,imshow(em) %em=external marker，外部标记符

结果显示一片黑，进一步查看：

>> whos em
  Name        Size             Bytes  Class      Attributes

em        253x251            63503  logical

>> max(em(:))
ans =
     0

原因：距离变换后得到的一些孤立的亮区域，那些区域包含分水岭脊线，而所有的黑色区域对应汇水盆地。em包含在所有亮区域内，不能连通。求出的分水岭变换全为1：没找到合适的分水岭脊线（不能形成对图像有效分割的分水岭脊线）。
实际上，应该对im2的补求分水岭变换，这样求得的分水岭脊线就在米粒之间，这个作为外部标记符。

下面对“扩展局部最小区域(im2)”的补求分水岭变换：

>> temp=bwdist(~im2);figure,imshow(temp,[])

和预料的一样，目标区域为汇水盆地。最亮的地方对应分水岭脊线。

>> Lim=watershed(temp);em=Lim==0;figure,imshow(em)
 
这些脊线应该作为外部标记才对!

⑤仿二值图像分水岭分割的方法，将此区域叠加到原图像上：

>> f(em)=255;figure,imshow(f,[])

可看到，过分割现象已不存在了，但出现了一点误分割（右下角处）。

⑥对于标记分水岭方法来说，第⑤步应该修改为：给出内部和外部标记后，就可以使用它们来修改梯度图像——使用称为“强制最小”的过程。方法如下：

>> g2=imimposemin(g,im2 | em);%im2为内部标记，em为外部标记
>> L2=watershed(g2);

>> f2=f;f2(L2==0)=255;figure,imshow(f2,[])

分割到的结果比较理想，但为什么会有一部分没有被分割？？
总结一下整个程序：

f=imread('rice_tophat.png');

h=fspecial('sobel');

fd=double(f);

g=sqrt(imfilter(fd,h,'replicate').^2+imfilter(fd,h','replicate').^2);

rm=imregionalmin(f);

im2=imextendedmin(f,45);%阈值要通过观察来定

im2=~im2;

temp=bwdist(im2);

Lim=watershed(temp);em=Lim==0;

g2=imimposemin(g,im2|em);

L2=watershed(g2);

f2=f;f2(L2==0)=255;

figure,imshow(f2,[])

//要根据实际情况来定，如中间不一定要经过im2=~im2;的转换。总体流程是这样的。