本文轉自：

clear,clc%三種方法進行分水嶺分割

%讀入影象

filename='sar1.bmp';

f=imread(filename);

Info=imfinfo(filename);

if Info.BitDepth>8   

f=rgb2gray(f);

end

figure,

mesh(double(f));%顯示影象，類似集水盆地

%方法1：一般分水嶺分割，從結果可以看出存在過分割問題

b=im2bw(f,graythresh(f));%二值化,注意應保證集水盆地的值較低（為0），否則就要對b取反

d=bwdist(b);             %求零值到最近非零值的距離，即集水盆地到分水嶺的距離  

  l=watershed(-d);         %matlab自帶分水嶺演算法，l中的零值即為風水嶺

w=l==0;                  %取出邊緣

g=b&~w;                  %用w作為mask從二值影象中取值

figure

subplot(2,3,1),

imshow(f);

subplot(2,3,2),

imshow(b);

subplot(2,3,3),

imshow(d);

subplot(2,3,4),

imshow(l);

subplot(2,3,5),

imshow(w);

subplot(2,3,6),

imshow(g);

%方法2：使用梯度的兩次分水嶺分割，從結果可以看出還存在過分割問題（在方法1的基礎上改進）

h=fspecial('sobel');%獲得縱方向的sobel運算元

fd=double(f);

g=sqrt(imfilter(fd,h,'replicate').^2+imfilter(fd,h','replicate').^2);%使用sobel運算元進行梯度運算

l=watershed(g);%分水嶺運算

wr=l==0;      

 g2=imclose(imopen(g,ones(3,3)),ones(3,3));%進行開閉運算對影象進行平滑

l2=watershed(g2);%再次進行分水嶺運算

wr2=l2==0;f2=f;

f2(wr2)=255;

figuresubplot(2,3,1),

imshow(f);

subplot(2,3,2),imshow(g);

subplot(2,3,3),imshow(l);

subplot(2,3,4),imshow(g2);

subplot(2,3,5),imshow(l2);

subplot(2,3,6),imshow(f2);

%方法3：使用梯度加掩模的三次分水嶺演算法（在方法2的基礎上改進）

h=fspecial('sobel');%獲得縱方向的sobel運算元

fd=double(f);

g=sqrt(imfilter(fd,h,'replicate').^2+imfilter(fd,h','replicate').^2);%使用sobel運算元進行梯度運算

l=watershed(g);%分水嶺運算

wr=l==0;

 rm=imregionalmin(g);  %計算影象的區域最小值定位，該函式僅僅是用來觀察為何分水嶺演算法產生這麼多集水盆地

im=imextendedmin(f,2);%上面僅是產生最小值點，而該函式則是得到最小值附近的區域，此處的附近是相差2的區域

fim=f;                 

  fim(im)=175;              %將im在原圖上標識出，用以觀察

 lim=watershed(bwdist(im));%再次分水嶺計算

em=lim==0;

g2=imimposemin(g,im|em);%在梯度圖上標出im和em，im是集水盆地的中心，em是分水嶺

l2=watershed(g2);       %第三次分水嶺計算

f2=f;f2(l2==0)=255;          %從原圖對分水嶺進行觀察

figuresubplot(3,3,1),imshow(f);

subplot(3,3,2),imshow(g);

subplot(3,3,3),imshow(l);

subplot(3,3,4),imshow(im);

subplot(3,3,5),imshow(fim);

subplot(3,3,6),imshow(lim);

subplot(3,3,7),imshow(g2);

subplot(3,3,8),imshow(l2);

subplot(3,3,9),imshow(f2);