分水岭算法

%******************************************************** %直接对梯度图像进行分水岭变换

rgb = imread('');%读取原图像

I = rgb2gray(rgb);%转化为灰度图像

imshow(I),title('原始灰度图像')

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);%求模

L = watershed(gradmag);%直接应用分水岭算法

Ig=L==0;

figure; imshow(Ig), %显示分割后的图像边界

title('梯度图像的分水岭变换边界（Ig）')

%******************************************************** %通过对图像进行一系列预处理，解决过分割问题

se = strel('disk', 20);%圆形结构元素

Io = imopen(I, se);%形态学开操作

figure; subplot(121)

imshow(Io), %显示执行开操作后的图像

title('开操作（Io）')

Ie = imerode(I, se);%对图像进行腐蚀

Iobr = imreconstruct(Ie, I);%形态学重建

subplot(122); imshow(Iobr), %显示重建后的图像

title('腐蚀后形态学重建(Iobr)')

Ioc = imclose(Io, se);%形态学关操作

figure; subplot(121)

imshow(Ioc), %显示关操作后的图像

title('关操作(Ioc)')

Iobrd = imdilate(Iobr, se);%对图像进行膨胀

Iobrcbr = imreconstruct(imcomplement(Iobrd), ...

imcomplement(Iobr));%形态学重建

Iobrcbr = imcomplement(Iobrcbr);%图像求反

subplot(122); imshow(Iobrcbr), %显示重建求反后的图像

title('膨胀后形态学重建求反(Iobrcbr)')

fgm = imregionalmax(Iobrcbr);%局部极大值

figure; imshow(fgm), %显示重建后局部极大值图像

title('重建后局部极大值图像(fgm)')

I2 = I;

I2(fgm) = 255;%局部极大值处像素值设置为255

figure; imshow(I2), %在原图上显示极大值区域

title('图上显示极大值区域(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('修改后的极大值区域（I3）')

bw = im2bw(Iobrcbr, graythresh(Iobrcbr));%转化为二值图像

subplot(122); imshow(bw), %显示二值图像

title('二值图像（bw）')

D = bwdist(bw);%计算距离

DL = watershed(D);%分水岭变换

bgm = DL == 0;%求取分割边界

figure; imshow(bgm), %显示分割后的边界

title('分水岭边界(bgm)')

gradmag2 = imimposemin(gradmag, bgm | fgm4);%设置最小值

L = watershed(gradmag2);%分水岭变换

I4 = I;

I4(imdilate(L == 0, ones(3, 3)) | bgm | fgm4) = 255;%前景及边界处设置为255

figure;imshow(I4)%突出前景及边界

title('前景和背景边界')