1. 圖像分割——分水嶺演算法
姓名：謝意遠
學號：19021110366T
嵌牛導讀：圖像中的目標物體是連接在一起的，則分割起來很困難，分水嶺分割演算法經常用於處理這類問題，通常會取得比較好的效果。
嵌牛鼻子：圖像分割、分水嶺演算法
嵌牛提問：分水嶺演算法具體有哪些步驟？
嵌牛正文：
一、綜述
分水嶺分割演算法把圖像看成一幅「地形圖」，其中亮度比較強的區域像素值較大，而比較暗的區域像素值較小，通過尋找「匯水盆地」和「分水嶺界限」，對圖像進行分割。而直接應用分水嶺分割演算法的效果往往並不好，如果在圖像中對前景對象和背景對象進行標注區別，再應用分水嶺演算法會取得較好的分割效果。基於標記控制的分水嶺分割方法有以下基本步驟：
1 綜述
分水嶺分割演算法把圖像看成一幅「地形圖」，其中亮度比較強的區域像素值較大，而比較暗的區域像素值較小，通過尋找「匯水盆地」和「分水嶺界限」，對圖像進行分割。直接應用分水嶺分割演算法的效果往往並不好，如果在圖像中對前景對象和背景對象進行標注區別，再應用分水嶺演算法會取得較好的分割效果。基於標記控制的分水嶺分割方法有以下基本步驟：
1.計算分割函數。圖像中較暗的區域是要分割的對象
2.計算前景標志。這些是每個對象內部連接的斑點像素。
3.計算背景標志。這些是不屬於任何對象的要素。
4.修改分割函數，使其僅在前景和後景標記位置有極小值。
5.對修改後的分割函數做分水嶺變換計算。
使用MATLAB圖像處理工具箱
註：期間用到了很多圖像處理工具箱的函數，例如fspecial、imfilter、watershed、label2rgb、imopen、imclose、imreconstruct、imcomplement、imregionalmax、bwareaopen、graythresh和imimposemin函數等。
2 步驟
第一步：讀入彩色圖像，將其轉化成灰度圖像
clc; clear all; close all;
rgb = imread('pears.png');
if ndims(rgb) == 3
I = rgb2gray(rgb);
else
I = rgb;
end
figure('units', 'normalized', 'position', [0 0 1 1]);
第2步：將梯度幅值作為分割函數
使用Sobel邊緣運算元對圖像進行水平和垂直方向的濾波，然後求取模值，sobel運算元濾波後的圖像在邊界處會顯示比較大的值，在沒有邊界處的值會很小。
hy = fspecial('sobel');
hx = hy';
Iy = imfilter(double(I), hy, 'replicate');
Ix = imfilter(double(I), hx, 'replicate');
gradmag = sqrt(Ix.^2 + Iy.^2);
figure('units', 'normalized', 'position', [0 0 1 1]);
subplot(1, 2, 1); imshow(I,[]), title('灰度圖像')
subplot(1, 2, 2); imshow(gradmag,[]), title('梯度幅值圖像')
可否直接對梯度幅值圖像使用分水嶺演算法？
L = watershed(gradmag);
Lrgb = label2rgb(L);
figure('units', 'normalized', 'position', [0 0 1 1]);
subplot(1, 2, 1); imshow(gradmag,[]), title('梯度幅值圖像')
subplot(1, 2, 2); imshow(Lrgb); title('梯度幅值做分水嶺變換')
直接使用梯度模值圖像進行分水嶺演算法得到的結果往往會存在過度分割的現象。因此通常需要分別對前景對象和背景對象進行標記，以獲得更好的分割效果。
第3步：標記前景對象
有多種方法可以應用在這里來獲得前景標記，這些標記必須是前景對象內部的連接斑點像素。這個例子中，將使用形態學技術「基於開的重建」和「基於閉的重建」來清理圖像。這些操作將會在每個對象內部創建單位極大值，使得可以使用imregionalmax來定位。
開運算和閉運算：先腐蝕後膨脹稱為開；先膨脹後腐蝕稱為閉。開和閉這兩種運算可以除去比結構元素小的特定圖像細節，同時保證不產生全局幾何失真。開運算可以把比結構元素小的突刺濾掉，切斷細長搭接而起到分離作用；閉運算可以把比結構元素小的缺口或孔填充上，搭接短的間隔而起到連接作用。
開操作是腐蝕後膨脹，基於開的重建（基於重建的開操作）是腐蝕後進行形態學重建。下面比較這兩種方式。首先，用imopen做開操作。
se = strel('disk', 20);
Io = imopen(I, se);
figure('units', 'normalized', 'position', [0 0 1 1]);
subplot(1, 2, 1); imshow(I, []); title('灰度圖像');
subplot(1, 2, 2); imshow(Io), title('圖像開操作')
接下來，通過腐蝕後重建來做基於開的重建計算。
Ie = imerode(I,se)
Iobr = imreconstruct(Ie,I);
figure('units', 'normalized', 'position', [0 0 1 1]);
subplot(1, 2, 1); imshow(I, []); title('灰度圖像');
subplot(1, 2, 2); imshow(Iobr, []), title('基於開的重建圖像')
開操作後，接著進行閉操作，可以移除較暗的斑點和枝幹標記。對比常規的形態學閉操作和基於閉的重建操作。首先，使用imclose：
Ioc = imclose(Io, se);
Ic = inclose(I,se);
figure('units', 'normalized', 'position', [0 0 1 1]);
subplot(2, 2, 1); imshow(I, []); title('灰度圖像');
subplot(2, 2, 2); imshow(Io, []); title('開操作圖像');
subplot(2, 2, 3); imshow(Ic, []); title('閉操作圖像');
subplot(2, 2, 4); imshow(Ioc, []), title('開閉操作');
現在使用imdilate，然後使用imreconstruct。注意必須對輸入圖像求補，對imreconstruct輸出圖像求補。IM2 = imcomplement(IM)計算圖像IM的補集。IM可以是二值圖像，或者RGB圖像。IM2與IM有著相同的數據類型和大小。
Iobrd = imdilate(Iobr, se);
Iobrcbr = imreconstruct(imcomplement(Iobrd), imcomplement(Iobr));
Iobrcbr = imcomplement(Iobrcbr);
figure('units', 'normalized', 'position', [0 0 1 1]);
subplot(2, 2, 1); imshow(I, []); title('灰度圖像');
subplot(2, 2, 2); imshow(Ioc, []); title('開閉操作');
subplot(2, 2, 3); imshow(Iobr, []); title('基於開的重建圖像');
subplot(2, 2, 4); imshow(Iobrcbr, []), title('基於閉的重建圖像');
通過比較Iobrcbr和loc可以看到，在移除小污點同時不影響對象全局形狀的應用下，基於重建的開閉操作要比標準的開閉重建更加有效。計算Iobrcbr的局部極大來得到更好的前景標記。
fgm = imregionalmax(Iobrcbr);
figure('units', 'normalized', 'position', [0 0 1 1]);
subplot(1, 3, 1); imshow(I, []); title('灰度圖像');
subplot(1, 3, 2); imshow(Iobrcbr, []); title('基於重建的開閉操作');
subplot(1, 3, 3); imshow(fgm, []); title('局部極大圖像');
為了幫助理解這個結果，疊加前景標記到原圖上。
It1 = rgb(:, :, 1);
It2 = rgb(:, :, 2);
It3 = rgb(:, :, 3);
It1(fgm) = 255; It2(fgm) = 0; It3(fgm) = 0;
I2 = cat(3, It1, It2, It3);
figure('units', 'normalized', 'position', [0 0 1 1]);
subplot(2, 2, 1); imshow(rgb, []); title('原圖像');
subplot(2, 2, 2); imshow(Iobrcbr, []); title('基於重建的開閉操作');
subplot(2, 2, 3); imshow(fgm, []); title('局部極大圖像');
subplot(2, 2, 4); imshow(I2); title('局部極大疊加到原圖像');
注意到大多閉塞處和陰影對象沒有被標記，這就意味著這些對象在結果中將不會得到合理的分割。而且，一些對象的前景標記會一直到對象的邊緣。這就意味著應該清理標記斑點的邊緣，然後收縮它們。可以通過閉操作和腐蝕操作來完成。
se2 = strel(ones(5,5));
fgm2 = imclose(fgm, se2);
fgm3 = imerode(fgm2, se2);
figure('units', 'normalized', 'position', [0 0 1 1]);
subplot(2, 2, 1); imshow(Iobrcbr, []); title('基於重建的開閉操作');
subplot(2, 2, 2); imshow(fgm, []); title('局部極大圖像');
subplot(2, 2, 3); imshow(fgm2, []); title('閉操作');
subplot(2, 2, 4); imshow(fgm3, []); title('腐蝕操作');
這個過程將會留下一些偏離的孤立像素，應該移除它們。可以使用bwareaopen，用來移除少於特定像素個數的斑點。BW2 = bwareaopen(BW,P)從二值圖像中移除所以少於P像素值的連通塊，得到另外的二值圖像BW2。
fgm4 = bwareaopen(fgm3, 20);
It1 = rgb(:, :, 1);
It2 = rgb(:, :, 2);
It3 = rgb(:, :, 3);
It1(fgm4) = 255; It2(fgm4) = 0; It3(fgm4) = 0;
I3 = cat(3, It1, It2, It3);
figure('units', 'normalized', 'position', [0 0 1 1]);
subplot(2, 2, 1); imshow(I2, []); title('局部極大疊加到原圖像');
subplot(2, 2, 2); imshow(fgm3, []); title('閉腐蝕操作');
subplot(2, 2, 3); imshow(fgm4, []); title('去除小斑點操作');
subplot(2, 2, 4); imshow(I3, []); title('修改局部極大疊加到原圖像');
第4步：計算背景標記
現在，需要標記背景。在清理後的圖像Iobrcbr中，暗像素屬於背景，所以可以從閾值操作開始。
bw =im2bw(Iobrcbr, graythresh(Iobrcbr));
figure('units', 'normalized', 'position', [0 0 1 1]);
subplot(1, 2, 1); imshow(Iobrcbr, []); title('基於重建的開閉操作');
subplot(1, 2, 2); imshow(bw, []); title('閾值分割');
背景像素在黑色區域，但是理想情形下，不必要求背景標記太接近於要分割的對象邊緣。通過計算「骨架影響范圍」來「細化」背景，或者SKIZ，bw的前景。這個可以通過計算bw的距離變換的分水嶺變換來實現，然後尋找結果的分水嶺脊線（DL==0）。D = bwdist(BW)計算二值圖像BW的歐幾里得矩陣。對BW的每一個像素，距離變換指定像素和最近的BW非零像素的距離。bwdist默認使用歐幾里得距離公式。BW可以由任意維數，D與BW有同樣的大小。
D = bwdist(bw);
DL = watershed(D);
bgm = DL == 0;
figure('units', 'normalized', 'position', [0 0 1 1]);
subplot(2, 2, 1); imshow(Iobrcbr, []); title('基於重建的開閉操作');
subplot(2, 2, 2); imshow(bw, []); title('閾值分割');
subplot(2, 2, 3); imshow(label2rgb(DL), []); title('分水嶺變換示意圖');
subplot(2, 2, 4); imshow(bgm, []); title('分水嶺變換脊線圖');
第5步：計算分割函數的分水嶺變換
函數imimposemin可以用來修改圖像，使其只是在特定的要求位置有局部極小。這里可以使用imimposemin來修改梯度幅值圖像，使其只在前景和後景標記像素有局部極小。
gradmag2 = imimposemin(gradmag, bgm | fgm4);
figure('units', 'normalized', 'position', [0 0 1 1]);
subplot(2,2,1)imshow(bgm,[]);title('分水嶺變換脊線圖');
subplot(2, 2, 2); imshow(fgm4, []); title('前景標記');
subplot(2, 2, 3); imshow(gradmag, []); title('梯度幅值圖像');
subplot(2, 2, 4); imshow(gradmag2, []); title('修改梯度幅值圖像');
最後，可以做基於分水嶺的圖像分割計算。
第6步：查看結果
一個可視化技術是疊加前景標記、背景標記、分割對象邊界到初始圖像。可以使用膨脹來實現某些要求，比如對象邊界，更加清晰可見。對象邊界定位於L==0的位置。
It1 = rgb(:, :, 1);
It2 = rgb(:, :, 2);
It3 = rgb(:, :, 3);
fgm5 = imdilate(L == 0, ones(3, 3)) | bgm | fgm4;
It1(fgm5) = 255; It2(fgm5) = 0; It3(fgm5) = 0;
I4 = cat(3, It1, It2, It3);
figure('units', 'normalized', 'position', [0 0 1 1]);
subplot(1, 2, 1); imshow(rgb, []); title('原圖像');
subplot(1, 2, 2); imshow(I4, []); title('標記和對象邊緣疊加到原圖像');
可視化說明了前景和後景標記如何影響結果。在幾個位置，部分的較暗對象與它們相鄰的較亮的鄰接對象相融合，這是因為受遮擋的對象沒有前景標記。
另外一個有用的可視化技術是將標記矩陣作為彩色圖像進行顯示。標記矩陣，比如通過watershed和bwlabel得到的，可以使用label2rgb轉換到真彩圖像來顯示。
Lrgb = label2rgb(L,'jet', 'w', 'shuffle');
figure('units', 'normalized', 'position', [0 0 1 1]);
subplot(1, 2, 1); imshow(rgb, []); title('原圖像');
subplot(1, 2, 2); imshow(Lrgb); title('彩色分水嶺標記矩陣');
可以使用透明度來疊加這個偽彩色標記矩陣在原亮度圖像上進行顯示。
figure('units', 'normalized', 'position', [0 0 1 1]);
subplot(1, 2, 1); imshow(rgb, []); title('原圖像');
subplot(1, 2, 2); imshow(rgb, []); hold on;
himage = imshow(Lrgb);
set(himage, 'AlphaData', 0.3);
title('標記矩陣疊加到原圖像');
2. 幫幫我 我不知道分水嶺演算法在圖像分割中的應用…… 代碼,還有別的幫幫我
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;
figure
subplot(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; %從原圖對分水嶺進行觀察
figure
subplot(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);
3. MATLAB的分水嶺演算法
樓上的到處都是這句話，無敵了你
4. MATLAB 分水嶺分割演算法
其實，這涉及到命令和演算法，單一的命令往往不能解決所有的問題，要有前處理或後處理，才能達到目的。另外，也說明，某個命令應該升級或更新了。所以，watershed這個命令，單用達不到所期望的效果，只有加上預處理才行。
5. 求經典MATLAB分水嶺演算法源代碼
%%都不知道你是啥樣的分割 什麼圖片
clear, close all;
clc;
PathName='d:\';%t為自填內容,下面p類似
FileName=[PathName 'test.jpg'];
Image=imread(FileName);
subplot(2,2,1);subimage(Image);title('原圖');;pixval on;
B=[1,1,1;1,1,1;1,1,1];%方形結構元
E8=[-1,0;-1,1;0,1;1,1;1,0;1,-1;0,-1;-1,-1]; % 8－連通結構元坐標
maskLenth=length(E8); % 結構元點的個數
[X,Y]=size(Image);
%原始圖像image 賦值給A1
n=1;
A(:,:,n)=Image;
M=zeros(X,Y);
Mark_Image=zeros(X,Y);
%產生距離圖
while sum(sum(A(:,:,n)))~=0
A(:,:,n+1)= imerode(A(:,:,n),B);
U(:,:,n)= (A(:,:,n)-A(:,:,n+1))*n;
M=M+U(:,:,n);
n=n+1;
end
n=n-1;
subplot(2,2,2);imagesc(M,[0,n]);title('距離圖');
% 搜尋局部最大值，將其放入Deal_Image
Deal_Image=zeros(X,Y);
while n>0
for high=1:X
for width=1:Y
%********************************************************************
Mark_Bool=0;
if M(high,width)==n
%______________________________________________________________
for dot=1:maskLenth
i=E8(dot,1); j=E8(dot,2);
if high+i>=1 & width+j>=1 & high+i<=X & width+j<=Y & M(high+i,width+j)>M(high,width);
Mark_Bool=1;break;
end % if_end
end % for dot_end
%______________________________________________________________
if Mark_Bool==0;
Deal_Image(high,width)=M(high,width);
end %if end
%______________________________________________________________
end %if end
%********************************************************************
end %for-end
end %for-end
n=n-1;
end % while n=0 end
Deal_Image =[Deal_Image>=1]
subplot(2,2,3);subimage(Deal_Image);title('輸出圖像');
Mark_Number=1;
while n>0
for high=1:X
for width=1:Y
Mark_Bool=0;
%********************************************************************
if M(high,width)==n
%______________________________________________________________
for dot=1:maskLenth
i=E8(dot,1); j=E8(dot,2);
if high+i>=1 & width+j>=1 & high+i<=X & width+j<=Y & Mark_Image(high+i,width+j)>0;
Mark_Image(high,width)=Mark_Image(high+i,width+j);
Mark_Bool=1;break;
end % if_end
end % for dot_end
%______________________________________________________________
if Mark_Bool==0;
Mark_Image(high,width)=Mark_Number;
Mark_Number=Mark_Number+1;
end %if end
%______________________________________________________________
pause;
subplot(2,2,2);imagesc(Mark_Image,[0,Mark_Number]);title('輸出圖像');
end %if end
%********************************************************************
end %for-end
end %for-end
n=n-1;
end % while n=0 end
subplot(2,2,3);imagesc(Mark_Image,[0,Mark_Number]);title('分割後的圖像');
uicontrol('Style','edit','string',['分割出區域：',Num2str(Mark_Number-1),'個'],...
'Position', [400 0 150 18],'FontSize',12,'FontWeight','light');