--------------------------------------------------------------------------------------------------------------------

附錄A  教程【3】給出的matlab原始碼，附詳細註釋

function [  ] = MarkerControlled_Watershed_tutorial(  )
%標記符控制的分水嶺演算法教程
%本程式按照官方文件的步驟，展示瞭如何使用分水嶺演算法，並測試對於含手勢圖片的分割結果是否理想
%程式修改自matlab官方文件 Marker-Controlled Watershed Segmentation
%網址： http://cn.mathworks.com/help/images/examples/marker-controlled-watershed-segmentation.html?s_tid=gn_loc_drop

%%
%官方文件簡介

%這個例子展示瞭如何使用分水嶺分割演算法來分割一副影象中的物體
%分水嶺變換經常應用到這類問題中。
%分水嶺變換通過將灰度影象看成拓撲表面（亮的畫素看做較高的點，暗畫素看做低點），以此在一幅圖中尋找“匯水盆地”和“分水嶺脊線”

%如果你可以識別或者說標記mark前景物體與背景的位置（提供一些先驗資訊），使用分水嶺變換的分割將會表現的很好。
%當然，人為提供這些標記“mark”，演算法就算有效，很大程度上也失去了意義，於是就有了下面可以自行標記前景背景的分水嶺演算法。

%Marker-Controlled Watershed Segmentation，即標記符控制的分水嶺分割演算法，主要遵循以下幾個基本步驟：
%1.計算一個分割函式（最基本的分水嶺分割）。你嘗試去分割的影象需要滿足：黑暗區域對應著待分割的目標物體
%2.計算前景標記。得到的前景標記是每個物體內部的聯通畫素斑塊。
%3.計算背景標記。背景標記由那些不屬於任何物體的的畫素點組成。
%4.修改分割函式，使得它只能在前景標記和背景標記的位置處達到最小值。
%5.計算修改後分割函式的分水嶺變換


close all;

%%
%步驟一：讀取彩色影象並轉化為灰度圖
HandImage = imread('testImgae/pear.jpg');
I = rgb2gray(HandImage);
imshow(I)

%%
%******************步驟二：使用梯度幅值作為分割函式**************************
%使用Sobel邊緣檢測運算元，imfilter函式，和一些簡單的四則運算來計算梯度幅值。
%梯度值總是在物體的邊緣處高，而總是在物體內部低。
hy = fspecial('sobel');%獲取sobel運算元模板，計算縱向梯度
hx = hy';%橫向模板
Iy = imfilter(double(I), hy, 'replicate');%計算縱向梯度
Ix = imfilter(double(I), hx, 'replicate');
gradmag = sqrt(Ix.^2 + Iy.^2);%計算梯度幅值
figure
imshow(gradmag,[]), title('Gradient magnitude (gradmag)')

%你能直接在梯度幅值矩陣上應用分水嶺變換來分割影象嗎？
L = watershed(gradmag);
Lrgb = label2rgb(L);
figure, imshow(Lrgb), title('Watershed transform of gradient magnitude (Lrgb)')
%從結果可以看到，顯然不能
%如果沒有額外的預處理，比如新增標記符（如下所述），直接使用分水嶺變換通常導致嚴重的過分割

%取消下面3行程式碼的註釋可以觀察label2rgb函式的作用
% testA=[1 1 2 2 2 2 ;1 2 3 3 3 4;0 1 1 3 4 4];
% testArgb = label2rgb(testA);
% figure, imshow(testArgb), title('測試label2rgb函式的作用，使用簡單的標記函式testA')

%%
%******************步驟三：標記前景物體************************************
%我們可以應用多種方法來獲得前景標記，這些標記必須是前景物件內部的連通畫素斑塊。
%這個例子中，將使用形態學技術“基於重建的開操作”和“基於重建的閉操作”來清理影象。
%這些操作將會在每個物件內部建立平坦的極大值小斑塊，這些斑塊可以使用imregionalmax函式來定位。

%我們的目標是

%開運算和閉運算：先腐蝕後膨脹稱為開；先膨脹後腐蝕稱為閉（所謂腐蝕還是膨脹是針對白色而言來理解的）
%開和閉這兩種運算可以除去比結構元素小的特定影象細節，同時保證不產生全域性幾何失真。
%開運算可以把比結構元素小的突刺濾掉，切斷細長搭接而起到分離作用；閉運算可以把比結構元素小的缺口或孔填充上，搭接短的間隔而起到連線作用。

%但是

%開操作的第一步腐蝕能去除白色小物體，隨後的膨脹趨向於恢復保留下來的物體的形狀，這種恢復是不精確的，其精確度取決於形狀和結構體之間的相似性。
%基於重建進行的開操作能夠準確的恢復腐蝕之後的物體形狀。
%所以基於重建的開操作與單純的開操作功能類似，但更好地保留了原物體的形狀。

%基於重建的開操作：I先腐蝕成Ie,後進行形態學重建。

%開操作是腐蝕後膨脹，基於重建的開操作是腐蝕後進行形態學重建。
%下面通過示例比較這兩種方式，直觀展示“基於重建的開操作”如何優於“開操作”
%首先，用imopen做開操作。
se = strel('disk',20);%disk指定構建一個圓形的結構體，第二個引數指定結構體的半徑
Io = imopen(I, se);%開操作
% figure
%imshow(Io), title('Opening (Io)')

%接下來，進行基於重建的開操作。使用imerode和imreconstruct函式實現
Ie = imerode(I, se);%先腐蝕‘erosion’
Iobr = imreconstruct(Ie, I);%再重建
% figure
%imshow(Iobr), title('Opening-by-reconstruction (Iobr)')


%對開操作後的結果進行閉操作，可以移除較暗的斑點和枝幹標記。
%同樣的，我們來對比常規的形態學閉操作和基於重建的閉操作的區別。

%首先，使用imclose對開操作的結果進行常規閉操作：
Ioc = imclose(Io, se);
% figure
%imshow(Ioc), title('Opening-closing 開操作+閉操作 (Ioc)')

%然後我們來對開操作的結果進行基於重建的閉操作
%即首先使用imdilate函式進行膨脹，然後使用imreconstruct進行重建。
%這裡要解釋一下，重建操作是通過膨脹操作和交運算定義的，即 Hk+1=（Hk 膨脹 結構體B）∩G。
%重建定義的具體解釋可以參考《數字影象處理的MATLAB實現（第2版）》P358頁
%重建操作是對標記影象，在模板影象的限制下，進行膨脹。本質上是一種*膨脹*演算法
%而這裡需要的是模擬閉操作對膨脹後的影象進行*腐蝕*，怎麼辦呢？
%可以通過對模板影象和標記影象同時求補，然後進行重建操作，對補的膨脹相當於對原圖的腐蝕。
%當然，重建的結果也要求補，才能得到實際的結果。
%IM2 = imcomplement(IM)計算影象IM的補集。IM可以是二值影象，或者RGB影象。IM2與IM有著相同的資料型別和大小。
Iobrd = imdilate(Iobr, se);%在基於重建的開操作的結果基礎上，進行腐蝕
Iobrcbr = imreconstruct(imcomplement(Iobrd), imcomplement(Iobr));%重建，標記影象為腐蝕後圖像取補，模板為腐蝕前原圖取補。
Iobrcbr = imcomplement(Iobrcbr);%重建結果再取補，得到實際基於重建的閉操作的結果。
figure
imshow(Iobrcbr), title('Opening-closing by reconstruction 基於重建的開+閉操作 (Iobrcbr)')
%如你所見，通過比較Iobrcbr和loc可以看到，在移除小汙點同時不影響物件全域性形狀的應用下，基於開閉操作的重建要比標準的開閉操作更加有效。


%計算Iobrcbr的區域極大值來得到好的前景標記。
%得到的前景標記圖fgm是二值圖，白色對應前景區域
fgm = imregionalmax(Iobrcbr);
figure
imshow(fgm), title('Regional maxima of opening-closing by reconstruction (fgm)')

%為了方便解釋結果，將前景標記圖疊加到原始影象上
I2 = I;
I2(fgm) = 255;%將fgm中的前景區域（畫素值為1）標記到原圖上（置白色）
figure
imshow(I2), title('Regional maxima superimposed on original image (I2)')

%注意到一些大部分重合或被陰影遮擋的物體沒有被標記出來。這意味著這些物體最終可能不會被正確的分割出來。
%並且，有些物體中前景標記正確的到達了物體的邊緣。這意味著你應該清除掉標記斑塊的邊緣，向內收縮一點。
%你可以通過先閉操作，再腐蝕做到這點。
se2 = strel(ones(5,5));
fgm2 = imclose(fgm, se2);
fgm3 = imerode(fgm2, se2);

%這個操作會導致遺留下一些離群的孤立點，這些是需要被移除的。
%你可以通過bwareaopen做到這點，函式將移除那些包含畫素點個數少於指定值的區域。
fgm4 = bwareaopen(fgm3, 20);
I3 = I;
I3(fgm4) = 255;
figure
imshow(I3)
title('Modified regional maxima superimposed on original image (fgm4)')

%%
%*********************第四步：計算背景標記**********************************
%本例中設計出的標記背景演算法的前提假設是：影象中相對亮的是物體，相對暗的區域是背景。
%如果不滿足這條假設，標記結果可能不甚理想

%現在你需要標記背景。在去除噪點後的影象Iobrcbr中，暗畫素屬於背景，所以你可以先進行一下閾值操作。
%bw = imbinarize(Iobrcbr);
bw=im2bw(Iobrcbr,graythresh(Iobrcbr));%我使用上一行程式碼報錯了，故換了一種二值化方法
figure
imshow(bw), title('Thresholded opening-closing by reconstruction (bw)')

%背景畫素現在是黑的了，但是理想情況下我們不希望背景標記太接近我們想要分割的物體的邊界。
%我們要使背景變瘦，通過計算“骨架影響範圍”來“細化”背景，或者SKIZ，bw的前景。這句沒翻通，原文如下
%We'll "thin" the background by computing the "skeleton by influence zones", or SKIZ, of the foreground of bw. 
%這個可以通過對bw的距離變換進行分水嶺變換來實現，然後尋找結果的分水嶺脊線（DL==0）。
D = bwdist(bw);
%D = bwdist(BW)計算二值影象BW的歐幾里得矩陣。對BW的每一個畫素，距離變換指定畫素和最近的BW非零畫素的距離。
%bwdist預設使用歐幾里得距離公式。BW可以由任意維數，D與BW有同樣的大小。

%由於bw中目標物體是白色的1.所以D中對應的目標物體處均是0，隨著進入背景越深，對應畫素值越大。
%這時正好符合我們使用分水嶺演算法的假設（想分出的目標物體數值較低）
%於是得到的背景標記是 物體與背景間的一個圈，能夠包住目標物體
DL = watershed(D);
bgm = DL == 0;%分水嶺變換結果L中，同一區域用同一數字表示，區域間分界線同一由0標識
figure
imshow(bgm), title('Watershed ridge lines (bgm)')


%%
%******************第五步：計算分割函式（修改後）的分水嶺變換*****************
% 函式imimposemin可以被用來修改一副圖片，使得其只在指定的位置處取得區域性最小值
% 這裡你可以使用imimposemin來修改梯度幅值影象，使得區域性最小值只出現在前景標記和背景標記處。
% 從結果來看imimposemin會將指定區域置為-Inf，從而成為極小值
% 且影象變得相當“平整”，一塊塊的相同數值的區域
gradmag2 = imimposemin(gradmag, bgm | fgm4);

%這段測試imimposemin都做了什麼
% gradmag2;
% temp=gradmag2-gradmag;
% 
% maxValue=max(max(temp));
% minValue=min(min(temp));
% temp=temp/maxValue;
% figure
% imshow(temp)

% maxValue=max(max(gradmag));
% temp_gradmag=gradmag/maxValue;
% figure
% imshow(temp_gradmag)
% 
% maxValue=max(max(gradmag2));
% temp_gradmag2=gradmag2/maxValue;
% figure
% imshow(temp_gradmag2)


%Finally we are ready to compute the watershed-based segmentation.
L = watershed(gradmag2);

%%
%******************第六步：視覺化結果************************************
I4 = I;
I4(imdilate(L == 0, ones(3, 3)) | bgm | fgm4) = 255;
figure
imshow(I4)
title('Markers and object boundaries superimposed on original image (I4)')


Lrgb = label2rgb(L, 'jet', 'w', 'shuffle');
figure
imshow(Lrgb)
title('Colored watershed label matrix (Lrgb)')

figure
imshow(I)
hold on
himage = imshow(Lrgb);
himage.AlphaData = 0.3;
title('Lrgb superimposed transparently on original image')
end
 

--------------------------------------------------------------------------------------------------------------------

附錄B  提煉後的原始碼，去掉多餘演示步驟，只可視化關鍵步驟

function [  ] = MarkerControlled_Watershed(  )
%標記符控制的分水嶺演算法
%與教程函式不同，本函式沒有多餘的步驟和註釋，直接實現
%展示的結果圖更少更精煉

close all;

%%
%步驟一：讀取彩色影象並轉化為灰度圖
HandImage = imread('testImgae/hand1.jpg');
I = rgb2gray(HandImage);

%如果原圖中，目標物體是較暗的，即與假設相反，這裡可以取反
% Fan=ones(size(I,1),size(I,2))*255;
% Fan=uint8(Fan);
% Fan=Fan-I;
% I=Fan;

imshow(I), title('原圖I')
%%
%******************步驟二：基於重建的開閉操作************************************

se = strel('disk',20);%disk指定構建一個圓形的結構體，第二個引數指定結構體的半徑

%接下來，進行基於重建的開操作。使用imerode和imreconstruct函式實現
Ie = imerode(I, se);%先腐蝕‘erosion’
Iobr = imreconstruct(Ie, I);%再重建

Iobrd = imdilate(Iobr, se);%在基於重建的開操作的結果基礎上，進行腐蝕
Iobrcbr = imreconstruct(imcomplement(Iobrd), imcomplement(Iobr));%重建，標記影象為腐蝕後圖像取補，模板為腐蝕前原圖取補。
Iobrcbr = imcomplement(Iobrcbr);%重建結果再取補，得到實際基於重建的閉操作的結果。
figure
imshow(Iobrcbr), title('基於重建的開+閉操作 (Iobrcbr)')


%%
%******************步驟三：標記前景物體************************************
%計算Iobrcbr的區域極大值來得到好的前景標記。
%得到的前景標記圖fgm是二值圖，白色對應前景區域
fgm = imregionalmax(Iobrcbr);

%為了方便解釋結果，將前景標記圖疊加到原始影象上
I2 = I;
I2(fgm) = 255;%將fgm中的前景區域（畫素值為1）標記到原圖上（置白色）

%注意到一些大部分重合或被陰影遮擋的物體沒有被標記出來。這意味著這些物體最終可能不會被正確的分割出來。
%並且，有些物體中前景標記正確的到達了物體的邊緣。這意味著你應該清除掉標記斑塊的邊緣，向內收縮一點。
%你可以通過先閉操作，再腐蝕做到這點。
se2 = strel(ones(5,5));
fgm2 = imclose(fgm, se2);
fgm3 = imerode(fgm2, se2);

%這個操作會導致遺留下一些離群的孤立點，這些是需要被移除的。
%你可以通過bwareaopen做到這點，函式將移除那些包含畫素點個數少於指定值的區域。
fgm4 = bwareaopen(fgm3, 20);
I3 = I;
I3(fgm4) = 255;

%%
%*********************第四步：計算背景標記**********************************
%本例中設計出的標記背景演算法的前提假設是：影象中相對亮的是物體，相對暗的區域是背景。
%如果不滿足這條假設，標記結果可能不甚理想

%現在你需要標記背景。在去除噪點後的影象Iobrcbr中，暗畫素屬於背景，所以你可以先進行一下閾值操作。
%bw = imbinarize(Iobrcbr);
bw=im2bw(Iobrcbr,graythresh(Iobrcbr));%我使用上一行程式碼報錯了，故換了一種二值化方法

%背景畫素現在是黑的了，但是理想情況下我們不希望背景標記太接近我們想要分割的物體的邊界。

D = bwdist(bw);
%D = bwdist(BW)計算二值影象BW的歐幾里得矩陣。對BW的每一個畫素，距離變換指定畫素和最近的BW非零畫素的距離。
%bwdist預設使用歐幾里得距離公式。BW可以由任意維數，D與BW有同樣的大小。

%由於bw中目標物體是白色的1.所以D中對應的目標物體處均是0，隨著進入背景越深，對應畫素值越大。
%這時正好符合我們使用分水嶺演算法的假設（想分出的目標物體數值較低）
%於是得到的背景標記是 物體與背景間的一個圈，能夠包住目標物體
DL = watershed(D);
bgm = DL == 0;%分水嶺變換結果L中，同一區域用同一數字表示，區域間分界線同一由0標識

%%
%******************第五步：計算分割函式（修改後）的分水嶺變換*****************
% 函式imimposemin可以被用來修改一副圖片，使得其只在指定的位置處取得區域性最小值
% 這裡你可以使用imimposemin來修改梯度幅值影象，使得區域性最小值只出現在前景標記和背景標記處。
% 從結果來看imimposemin會將指定區域置為-Inf，從而成為極小值
% 且影象變得相當“平整”，一塊塊的相同數值的區域


%使用Sobel邊緣檢測運算元，imfilter函式，和一些簡單的四則運算來計算梯度幅值。
%梯度值總是在物體的邊緣處高，而總是在物體內部低。
hy = fspecial('sobel');%獲取sobel運算元模板，計算縱向梯度
hx = hy';%橫向模板
Iy = imfilter(double(I), hy, 'replicate');%計算縱向梯度
Ix = imfilter(double(I), hx, 'replicate');
gradmag = sqrt(Ix.^2 + Iy.^2);%計算梯度幅值

gradmag2 = imimposemin(gradmag, bgm | fgm4);

%Finally we are ready to compute the watershed-based segmentation.
L = watershed(gradmag2);

%%
%******************第六步：視覺化結果************************************
I4 = I;
I4(imdilate(L == 0, ones(3, 3)) | bgm | fgm4) = 255;
figure
imshow(I4)
title('在原圖上繪製前景、背景標記，以及分割邊界')


Lrgb = label2rgb(L, 'jet', 'w', 'shuffle');
figure
imshow(Lrgb)
title('Colored watershed label matrix (Lrgb)')

figure
imshow(I)
hold on
himage = imshow(Lrgb);
himage.AlphaData = 0.3;
title('Lrgb superimposed transparently on original image')


end